text:
#----------------------------------------------------------------------------------------------#
#                                    Amount Type Interface                                     #
#----------------------------------------------------------------------------------------------#

"""
`function deco end`\n
Interface to return a unique decorative `Symbol` from a method's argument type.
"""
function deco end

"""
`function ppu end`\n
Interface to pretty-print units.
"""
function ppu end

# A 191113-212130 benchmark showed amt(x) is faster than x.amt:
#
# ```julia-repl
# julia> u1 = u(1.0f0 ± 0.1f0, MO)
# ū₃₂: (1 ± 0.1) kJ/kmol
#
# julia> typeof(u1)
# uAmt{Float32,MM,MO}
#
# julia> @benchmark u1.amt
# ✂ ✂ ✂   median time:      26.918 ns (0.00% GC)   ✂ ✂ ✂
#
# julia> @benchmark amt(u1)
# ✂ ✂ ✂   median time:      16.710 ns (0.00% GC)   ✂ ✂ ✂
#
# ```

"""
`function amt end`\n
Interface to get an `AMOUNTS`' `:amt` field in a type-stable manner.
"""
function amt end

export deco, ppu, amt

#----------------------------------------------------------------------------------------------#
#                                     Generic Amount Type                                      #
#----------------------------------------------------------------------------------------------#

import Base: cp, convert
import Unicode: normalize
import Base: +, -, *, /

"""
Generic Amount type factory.
"""
function mkGenAmt(TYPE::Symbol,         # Type name:            :_Amt
                  SUPT::Symbol,         # Supertype:            :GenerAmt
                  SYMB::AbstractString, # Printing symbol:      "?"
                  WHAT::AbstractString, # Description:          "generic amounts"
                  DELT::Bool=false,     # Whether a Δ quantity
                 )
    # Constants
    i, f = DELT ? (3, 4) : (1, 2)
    𝑠SY = normalize((DELT ? "Δ" : "") * string(SYMB))
    # Documentation
    hiStr = tyArchy(eval(SUPT))
    dcStr = """
`struct $TYPE{𝗽<:PREC,𝘅<:EXAC} <: $SUPT{𝗽,𝘅}`\n
Precision-, and Exactness- parametric $WHAT amounts based in arbitrary units.\n
`$TYPE{𝗽,𝘅}` parameters are:\n
- Precision `𝗽<:Union{Float16,Float32,Float64,BigFloat}`;\n
- Exactness `𝘅<:Union{EX,MM}`, i.e., either a single, precise value or an uncertainty-bearing
  measurement, respectively;\n
A `$TYPE` can be natively constructed from the following argument types:\n
- A plain, unitless float;\n
- A plain, unitless `Measurement`; hence, any `AbstractFloat`;\n
- A `Quantity{AbstractFloat}` with any units.\n
## Hierarchy\n
`$(TYPE) <: $(hiStr)`
    """
    # @eval block
    @eval begin
        # Concrete type definition
        struct $TYPE{𝗽,𝘅} <: $SUPT{𝗽,𝘅}
            amt::UATY{𝗽} where 𝗽<:PREC
            # Inner, non-converting, parameter-determining constructors
            $TYPE(x::$TYPE{𝗽,𝘅}) where {𝗽<:PREC,𝘅<:EXAC} = new{𝗽,𝘅}(amt(x))
            $TYPE(x::Union{𝗽,UETY{𝗽}}) where 𝗽<:PREC = new{𝗽,EX}(_qty(x))
            $TYPE(x::Union{PMTY{𝗽},UMTY{𝗽}}) where 𝗽<:PREC = new{𝗽,MM}(_qty(x))
            # Inner, non-converting, fully-specified constructors
            (::Type{$TYPE{𝗽,EX}})(x::𝗽) where 𝗽<:PREC = new{𝗽,EX}(_qty(x))
            (::Type{$TYPE{𝗽,EX}})(x::PMTY{𝗽}) where 𝗽<:PREC = new{𝗽,EX}(_qty(x.val))
            (::Type{$TYPE{𝗽,MM}})(x::𝗽) where 𝗽<:PREC = new{𝗽,MM}(_qty(measurement(x)))
            (::Type{$TYPE{𝗽,MM}})(x::PMTY{𝗽}) where 𝗽<:PREC = new{𝗽,MM}(_qty(x))
        end
        # Type documentation
        @doc $dcStr $TYPE
        # Precision-changing external constructors
        (::Type{$TYPE{𝘀}})(x::$TYPE{𝗽,EX}) where {𝘀<:PREC,𝗽<:PREC} = begin
            $TYPE(𝘀(amt(x).val) * unit(amt(x)))
        end
        (::Type{$TYPE{𝘀}})(x::$TYPE{𝗽,MM}) where {𝘀<:PREC,𝗽<:PREC} = begin
            $TYPE(Measurement{𝘀}(amt(x).val) * unit(amt(x)))
        end
        # Precision+Exactness-changing external constructors
        (::Type{$TYPE{𝘀,EX}})(x::$TYPE{𝗽,EX}) where {𝘀<:PREC,𝗽<:PREC} = begin
            $TYPE(𝘀(amt(x).val) * unit(amt(x)))
        end
        (::Type{$TYPE{𝘀,EX}})(x::$TYPE{𝗽,MM}) where {𝘀<:PREC,𝗽<:PREC} = begin
            $TYPE(𝘀(amt(x).val.val) * unit(amt(x)))
        end
        (::Type{$TYPE{𝘀,MM}})(x::$TYPE{𝗽,EX},
                              e::𝘀=𝘀(max(eps(𝘀), eps(amt(x).val)))) where {𝘀<:PREC,
                                                                           𝗽<:PREC} = begin
            $TYPE(measurement(𝘀(amt(x).val), e) * unit(amt(x)))
        end
        (::Type{$TYPE{𝘀,MM}})(x::$TYPE{𝗽,MM}) where {𝘀<:PREC,𝗽<:PREC} = begin
            $TYPE(Measurement{𝘀}(amt(x).val) * unit(amt(x)))
        end
        # Type export
        export $TYPE
        # Type-stabler wrapped amount obtaining function
        amt(x::$TYPE{𝗽,EX}) where 𝗽<:PREC = x.amt::Quantity{𝗽}
        amt(x::$TYPE{𝗽,MM}) where 𝗽<:PREC = x.amt::Quantity{Measurement{𝗽}}
        # Type-specific functions
        deco(x::$TYPE{𝗽,𝘅} where {𝗽,𝘅}) = Symbol($𝑠SY)
        ppu(x::$TYPE{𝗽,𝘅} where {𝗽,𝘅}) = string(unit(amt(x)))
        # Conversions
        convert(::Type{$TYPE{𝘀,𝘅}},
                y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘅<:EXAC} = begin
            $TYPE{promote_type(𝘀,𝗽),𝘅}(y)
        end
        convert(::Type{$TYPE{𝘀,𝘆}},
                y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE{promote_type(𝘀,𝗽),promote_type(𝘆,𝘅)}(y)
        end
        # Promotion rules
        promote_rule(::Type{$TYPE{𝘀,𝘆}},
                     ::Type{$TYPE{𝗽,𝘅}}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE{promote_type(𝘀,𝗽),promote_type(𝘆,𝘅)}
        end
        # same-type sum,sub with Unitful promotion
        +(x::$TYPE{𝘀,𝘆}, y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE(+(amt(x), amt(y)))
        end
        -(x::$TYPE{𝘀,𝘆}, y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE(-(amt(x), amt(y)))
        end
        # scalar mul,div with Unitful promotion
        *(y::plnF{𝘀}, x::$TYPE{𝗽}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(*(amt(x), y))
        *(x::$TYPE{𝗽}, y::plnF{𝘀}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(*(amt(x), y))
        /(x::$TYPE{𝗽}, y::plnF{𝘀}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(/(amt(x), y))
        # Type-preserving scalar mul,div
        *(y::REAL, x::$TYPE{𝗽}) where 𝗽<:PREC = $TYPE(*(amt(x), 𝗽(y)))
        *(x::$TYPE{𝗽}, y::REAL) where 𝗽<:PREC = $TYPE(*(amt(x), 𝗽(y)))
        /(x::$TYPE{𝗽}, y::REAL) where 𝗽<:PREC = $TYPE(/(amt(x), 𝗽(y)))
    end
end

#----------------------------------------------------------------------------------------------#
#                                 Generic Amount Declarations                                  #
#----------------------------------------------------------------------------------------------#

# The fallback generic amount
mkGenAmt(:_Amt  , :GenerAmt , "?"   , "generic amounts"     , false )

#----------------------------------------------------------------------------------------------#
#                                  Whole Amount Type Factory                                   #
#----------------------------------------------------------------------------------------------#

"""
Whole Amount type factory.
"""
function mkWhlAmt(TYPE::Symbol,         # Type name:            :sysT
                  SUPT::Symbol,         # Supertype:            :WProperty
                  FNAM::Symbol,         # Function Name:        :T
                  SYMB::AbstractString, # Printing symbol:      "T"
                  UNIT::Unitful.Units,  # SY quantity units:    u"K"
                  USTR::AbstractString, # PrettyPrinting units: "K"
                  WHAT::AbstractString, # Description:          "temperature"
                  DELT::Bool=false,     # Whether a Δ quantity
                 )
    # Constants
    uSY = UNIT
    𝑢SY = typeof(uSY)
    𝑑SY = dimension(uSY)
    i, f = DELT ? (3, 4) : (1, 2)
    𝑠SY = normalize((DELT ? "Δ" : "") * string(SYMB))
    # Documentation
    hiStr = tyArchy(eval(SUPT))
    dcStr = """
`struct $TYPE{𝗽<:PREC,𝘅<:EXAC} <: $SUPT{𝗽,𝘅}`\n
Precision-, and Exactness- parametric $WHAT amounts based in $USTR.\n
`$TYPE{𝗽,𝘅}` parameters are:\n
- Precision `𝗽<:Union{Float16,Float32,Float64,BigFloat}`;\n
- Exactness `𝘅<:Union{EX,MM}`, i.e., either a single, precise value or an uncertainty-bearing
  measurement, respectively;\n
A `$TYPE` can be natively constructed from the following argument types:\n
- A plain, unitless float;\n
- A plain, unitless `Measurement`; hence, any `AbstractFloat`;\n
- A `Quantity{AbstractFloat}` with compatible units.\n
Constructors determine all parameters from their arguments.\n
## Hierarchy\n
`$(TYPE) <: $(hiStr)`
    """
    # @eval block
    @eval begin
        # Concrete type definition
        struct $TYPE{𝗽,𝘅} <: $SUPT{𝗽,𝘅}
            amt::UATY{𝗽,$𝑑SY,$𝑢SY} where 𝗽<:PREC
            # Inner, non-converting, parameter-determining constructors
            # ---------------------------------------------------------
            # Copy constructor
            $TYPE(x::$TYPE{𝗽,𝘅}) where {𝗽<:PREC,𝘅<:EXAC} = new{𝗽,𝘅}(amt(x))
            # Plain constructors enforce default units & avoid unit conversion
            $TYPE(x::𝗽) where 𝗽<:PREC = new{𝗽,EX}(_qty(x * $uSY))
            $TYPE(x::PMTY{𝗽}) where 𝗽<:PREC = new{𝗽,MM}(_qty(x * $uSY))
            # Quantity constructors have to perform unit conversion despite matching dimensions
            $TYPE(x::UETY{𝗽,$𝑑SY}) where 𝗽<:PREC = new{𝗽,EX}(_qty(uconvert($uSY, x)))
            $TYPE(x::UMTY{𝗽,$𝑑SY}) where 𝗽<:PREC = new{𝗽,MM}(_qty(uconvert($uSY, x)))
            # Inner, non-converting, fully-specified constructors
            # ---------------------------------------------------
            (::Type{$TYPE{𝗽,EX}})(x::𝗽) where 𝗽<:PREC = new{𝗽,EX}(_qty(x * $uSY))
            (::Type{$TYPE{𝗽,EX}})(x::PMTY{𝗽}) where 𝗽<:PREC = new{𝗽,EX}(_qty(x.val * $uSY))
            (::Type{$TYPE{𝗽,MM}})(x::𝗽) where 𝗽<:PREC = new{𝗽,MM}(_qty(measurement(x) * $uSY))
            (::Type{$TYPE{𝗽,MM}})(x::PMTY{𝗽}) where 𝗽<:PREC = new{𝗽,MM}(_qty(x * $uSY))
        end
        # Type documentation
        @doc $dcStr $TYPE
        # Precision-changing external constructors
        (::Type{$TYPE{𝘀}})(x::$TYPE{𝗽,EX}
                          ) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(𝘀(amt(x).val))
        (::Type{$TYPE{𝘀}})(x::$TYPE{𝗽,MM}
                          ) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(Measurement{𝘀}(amt(x).val))
        # Precision+Exactness-changing external constructors
        (::Type{$TYPE{𝘀,EX}})(x::$TYPE{𝗽,EX}
                             ) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(𝘀(amt(x).val))
        (::Type{$TYPE{𝘀,EX}})(x::$TYPE{𝗽,MM}
                             ) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(𝘀(amt(x).val.val))
        (::Type{$TYPE{𝘀,MM}})(x::$TYPE{𝗽,EX},
                              e::𝘀=𝘀(max(eps(𝘀),eps(amt(x).val)))
                             ) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(measurement(𝘀(amt(x).val), e))
        (::Type{$TYPE{𝘀,MM}})(x::$TYPE{𝗽,MM}
                             ) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(Measurement{𝘀}(amt(x).val))
        # Type export
        export $TYPE
        # Type-stable wrapped amount obtaining function
        amt(x::$TYPE{𝗽,EX}) where 𝗽<:PREC = x.amt::Quantity{𝗽,$𝑑SY,$𝑢SY}
        amt(x::$TYPE{𝗽,MM}) where 𝗽<:PREC = x.amt::Quantity{Measurement{𝗽},$𝑑SY,$𝑢SY}
        # Type-specific functions
        deco(x::$TYPE{𝗽,𝘅} where {𝗽,𝘅}) = Symbol($𝑠SY)
        ppu(x::$TYPE{𝗽,𝘅} where {𝗽,𝘅}) = $USTR
        # Indirect construction from plain
        $FNAM(x::plnF) = $TYPE(x)
        $FNAM(x::REAL) = $TYPE(float(x))
        # Indirect construction from quantity
        $FNAM(x::UATY{𝗽,$𝑑SY}) where 𝗽<:PREC = $TYPE(x)
        $FNAM(x::uniR{𝗽,$𝑑SY}) where 𝗽<:REAL = $TYPE(float(x.val) * unit(x))
        export $FNAM
        # Conversions
        convert(::Type{$TYPE{𝘀,𝘅}},
                y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘅<:EXAC} = begin
            $TYPE{promote_type(𝘀,𝗽),𝘅}(y)
        end
        convert(::Type{$TYPE{𝘀,𝘆}},
                y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE{promote_type(𝘀,𝗽),promote_type(𝘆,𝘅)}(y)
        end
        # Promotion rules
        promote_rule(::Type{$TYPE{𝘀,𝘆}},
                     ::Type{$TYPE{𝗽,𝘅}}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE{promote_type(𝘀,𝗽),promote_type(𝘆,𝘅)}
        end
        # same-type sum,sub with Unitful promotion
        +(x::$TYPE{𝘀,𝘆}, y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE(+(amt(x), amt(y)))
        end
        -(x::$TYPE{𝘀,𝘆}, y::$TYPE{𝗽,𝘅}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC} = begin
            $TYPE(-(amt(x), amt(y)))
        end
        # scalar mul,div with Unitful promotion
        *(y::plnF{𝘀}, x::$TYPE{𝗽}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(*(amt(x), y))
        *(x::$TYPE{𝗽}, y::plnF{𝘀}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(*(amt(x), y))
        /(x::$TYPE{𝗽}, y::plnF{𝘀}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(/(amt(x), y))
        # Type-preserving scalar mul,div
        *(y::REAL, x::$TYPE{𝗽}) where 𝗽<:PREC = $TYPE(*(amt(x), 𝗽(y)))
        *(x::$TYPE{𝗽}, y::REAL) where 𝗽<:PREC = $TYPE(*(amt(x), 𝗽(y)))
        /(x::$TYPE{𝗽}, y::REAL) where 𝗽<:PREC = $TYPE(/(amt(x), 𝗽(y)))
    end
end

#----------------------------------------------------------------------------------------------#
#                           Thermodynamic Whole Amount Declarations                            #
#----------------------------------------------------------------------------------------------#

# Regular properties -- \bb#<TAB> velocity/speed function names
mkWhlAmt(:sysT  , :WProperty, :T    , "T"   , u"K"          , "K"       , "temperature"         , false )
mkWhlAmt(:sysP  , :WProperty, :P    , "P"   , u"kPa"        , "kPa"     , "pressure"            , false )
mkWhlAmt(:VELO  , :WProperty, :velo , "𝕍"   , u"√(kJ/kg)"   , "√kJ/kg"  , "velocity"            , false )
mkWhlAmt(:SPEE  , :WProperty, :spee , "𝕧"   , u"m/s"        , "m/s"     , "speed"               , false )

# Regular unranked -- \sans#<TAB> function names
mkWhlAmt(:TIME  , :WUnranked, :TIME , "𝗍"   , u"s"          , "s"       , "time"                , false )
mkWhlAmt(:grav  , :WUnranked, :grav , "𝗀"   , u"m/s^2"      , "m/s²"    , "gravity"             , false )
mkWhlAmt(:alti  , :WUnranked, :alti , "𝗓"   , u"m"          , "m"       , "altitude"            , false )

#----------------------------------------------------------------------------------------------#
#                                  Based Amount Type Factory                                   #
#----------------------------------------------------------------------------------------------#

"""
Based Amount type factory.
"""
function mkBasAmt(TYPE::Symbol,         # Type Name:            :uAmt
                  SUPT::Symbol,         # Supertype:            :BProperty
                  FNAM::Symbol,         # Function Name:        :u
                  SYMB::AbstractString, # Printing symbol:      "U"
                  UNIT::Unitful.Units,  # SY quantity units:    u"kJ"
                  USTR::AbstractString, # PrettyPrinting units: "K"
                  WHAT::AbstractString, # Description:          "internal energy"
                  DELT::Bool=false;     # Whether a Δ quantity
                  bsym::NTuple{4,Symbol}=(:none,:none,:none,:none)
                 )
    # Constants
    uSY = UNIT
    uDT = UNIT / u"s"
    uMA = UNIT / u"kg"
    uMO = UNIT / u"kmol"
    𝑢SY = typeof(uSY)
    𝑢DT = typeof(uDT)
    𝑢MA = typeof(uMA)
    𝑢MO = typeof(uMO)
    𝑑SY = dimension(uSY)
    𝑑DT = dimension(uDT)
    𝑑MA = dimension(uMA)
    𝑑MO = dimension(uMO)
    i, f = DELT ? (3, 4) : (1, 2)
    𝑠SY = bsym[1] == :none ?
        normalize((DELT ? "Δ" : "") * uppercase(string(SYMB))) :
        string(bsym[1])
    𝑠DT = bsym[2] == :none ?
        normalize(string(𝑠SY[1:i], "\u0307", 𝑠SY[f:end])) :
        string(bsym[2])
    𝑠MA = bsym[3] == :none ?
        normalize((DELT ? "Δ" : "") * lowercase(string(SYMB))) :
        string(bsym[3])
    𝑠MO = bsym[4] == :none ?
        normalize(string(𝑠MA[1:i], "\u0304", 𝑠MA[f:end])) :
        string(bsym[4])
    # Documentation
    hiStr = tyArchy(eval(SUPT))
    dcStr = """
`struct $TYPE{𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE} <: $SUPT{𝗽,𝘅,𝗯}`\n
Precision-, Exactness-, and Base- parametric $WHAT amounts based in $USTR.\n
`$TYPE{𝗽,𝘅,𝗯}` parameters are:\n
- Precision `𝗽<:Union{Float16,Float32,Float64,BigFloat}`;\n
- Exactness `𝘅<:Union{EX,MM}`, i.e., either a single, precise value or an uncertainty-bearing
  measurement, respectively;\n
- Thermodynamic base `𝗯<:Union{SY,DT,MA,MO}` respectively for system, rate, mass, or molar
  quantities, respectively in units of $(uSY), $(uDT), $(uMA), or $(uMO).\n
A `$TYPE` can be natively constructed from the following argument types:\n
- A plain, unitless float;\n
- A plain, unitless `Measurement`; hence, any `AbstractFloat`;\n
- A `Quantity{AbstractFloat}` with compatible units.\n
Constructors determine parameters from their arguments. `Quantity` constructors do not need a
base argument. Plain, `AbstractFloat` ones require the base argument.\n
## Hierarchy\n
`$(TYPE) <: $(hiStr)`
    """
    # @eval block
    @eval begin
        # Concrete type definition
        struct $TYPE{𝗽,𝘅,𝗯} <: $SUPT{𝗽,𝘅,𝗯}
            amt::Union{UATY{𝗽,$𝑑SY,$𝑢SY},UATY{𝗽,$𝑑DT,$𝑢DT},
                       UATY{𝗽,$𝑑MA,$𝑢MA},UATY{𝗽,$𝑑MO,$𝑢MO}} where 𝗽<:PREC
            # Inner, non-converting, parameter-determining constructors
            # ---------------------------------------------------------
            # Copy constructor
            $TYPE(x::$TYPE{𝗽,𝘅,𝗯}) where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE} = new{𝗽,𝘅,𝗯}(amt(x))
            # Plain constructors enforce default units & avoid unit conversion
            # Plain Exact (𝗽<:PREC) float constructors
            $TYPE(x::𝗽, ::Type{SY}) where 𝗽<:PREC = new{𝗽,EX,SY}(_qty(x * $uSY))
            $TYPE(x::𝗽, ::Type{DT}) where 𝗽<:PREC = new{𝗽,EX,DT}(_qty(x * $uDT))
            $TYPE(x::𝗽, ::Type{MA}) where 𝗽<:PREC = new{𝗽,EX,MA}(_qty(x * $uMA))
            $TYPE(x::𝗽, ::Type{MO}) where 𝗽<:PREC = new{𝗽,EX,MO}(_qty(x * $uMO))
            # Plain Measurement (PMTY) constructors
            $TYPE(x::PMTY{𝗽}, ::Type{SY}) where 𝗽<:PREC = new{𝗽,MM,SY}(_qty(x * $uSY))
            $TYPE(x::PMTY{𝗽}, ::Type{DT}) where 𝗽<:PREC = new{𝗽,MM,DT}(_qty(x * $uDT))
            $TYPE(x::PMTY{𝗽}, ::Type{MA}) where 𝗽<:PREC = new{𝗽,MM,MA}(_qty(x * $uMA))
            $TYPE(x::PMTY{𝗽}, ::Type{MO}) where 𝗽<:PREC = new{𝗽,MM,MO}(_qty(x * $uMO))
            # Quantity constructors have to perform unit conversion despite matching dimensions
            # United Exact (UETY) constructors
            $TYPE(x::UETY{𝗽,$𝑑SY}) where 𝗽<:PREC = new{𝗽,EX,SY}(_qty(uconvert($uSY, x)))
            $TYPE(x::UETY{𝗽,$𝑑DT}) where 𝗽<:PREC = new{𝗽,EX,DT}(_qty(uconvert($uDT, x)))
            $TYPE(x::UETY{𝗽,$𝑑MA}) where 𝗽<:PREC = new{𝗽,EX,MA}(_qty(uconvert($uMA, x)))
            $TYPE(x::UETY{𝗽,$𝑑MO}) where 𝗽<:PREC = new{𝗽,EX,MO}(_qty(uconvert($uMO, x)))
            # United Measurement (UMTY) constructors
            $TYPE(x::UMTY{𝗽,$𝑑SY}) where 𝗽<:PREC = new{𝗽,MM,SY}(_qty(uconvert($uSY, x)))
            $TYPE(x::UMTY{𝗽,$𝑑DT}) where 𝗽<:PREC = new{𝗽,MM,DT}(_qty(uconvert($uDT, x)))
            $TYPE(x::UMTY{𝗽,$𝑑MA}) where 𝗽<:PREC = new{𝗽,MM,MA}(_qty(uconvert($uMA, x)))
            $TYPE(x::UMTY{𝗽,$𝑑MO}) where 𝗽<:PREC = new{𝗽,MM,MO}(_qty(uconvert($uMO, x)))
            # Inner, non-converting, fully-specified constructors
            # ---------------------------------------------------
            # SY-based constructors
            (::Type{$TYPE{𝗽,EX,SY}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,EX,SY}(_qty(             x * $uSY))
            end
            (::Type{$TYPE{𝗽,EX,SY}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,EX,SY}(_qty(         x.val * $uSY))
            end
            (::Type{$TYPE{𝗽,MM,SY}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,MM,SY}(_qty(measurement(x) * $uSY))
            end
            (::Type{$TYPE{𝗽,MM,SY}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,MM,SY}(_qty(             x * $uSY))
            end
            # DT-based constructors
            (::Type{$TYPE{𝗽,EX,DT}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,EX,DT}(_qty(             x * $uDT))
            end
            (::Type{$TYPE{𝗽,EX,DT}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,EX,DT}(_qty(         x.val * $uDT))
            end
            (::Type{$TYPE{𝗽,MM,DT}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,MM,DT}(_qty(measurement(x) * $uDT))
            end
            (::Type{$TYPE{𝗽,MM,DT}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,MM,DT}(_qty(             x * $uDT))
            end
            # MA-based constructors
            (::Type{$TYPE{𝗽,EX,MA}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,EX,MA}(_qty(             x * $uMA))
            end
            (::Type{$TYPE{𝗽,EX,MA}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,EX,MA}(_qty(         x.val * $uMA))
            end
            (::Type{$TYPE{𝗽,MM,MA}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,MM,MA}(_qty(measurement(x) * $uMA))
            end
            (::Type{$TYPE{𝗽,MM,MA}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,MM,MA}(_qty(             x * $uMA))
            end
            # MO-based constructors
            (::Type{$TYPE{𝗽,EX,MO}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,EX,MO}(_qty(             x * $uMO))
            end
            (::Type{$TYPE{𝗽,EX,MO}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,EX,MO}(_qty(         x.val * $uMO))
            end
            (::Type{$TYPE{𝗽,MM,MO}})(x::𝗽) where 𝗽<:PREC = begin
                new{𝗽,MM,MO}(_qty(measurement(x) * $uMO))
            end
            (::Type{$TYPE{𝗽,MM,MO}})(x::PMTY{𝗽}) where 𝗽<:PREC = begin
                new{𝗽,MM,MO}(_qty(             x * $uMO))
            end
        end
        # Type documentation
        @doc $dcStr $TYPE
        # Precision-changing external constructors
        (::Type{$TYPE{𝘀}})(x::$TYPE{𝗽,EX,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝗯<:BASE} = begin
            $TYPE(𝘀(amt(x).val), 𝗯)
        end
        (::Type{$TYPE{𝘀}})(x::$TYPE{𝗽,MM,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝗯<:BASE} = begin
            $TYPE(Measurement{𝘀}(amt(x).val), 𝗯)
        end
        # Precision+Exactness-changing external constructors
        (::Type{$TYPE{𝘀,EX}})(x::$TYPE{𝗽,EX,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝗯<:BASE} = begin
            $TYPE(𝘀(amt(x).val), 𝗯)
        end
        (::Type{$TYPE{𝘀,EX}})(x::$TYPE{𝗽,MM,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝗯<:BASE} = begin
            $TYPE(𝘀(amt(x).val.val), 𝗯)
        end
        (::Type{$TYPE{𝘀,MM}})(x::$TYPE{𝗽,EX,𝗯},
                            e::𝘀=𝘀(max(eps(𝘀),eps(amt(x).val)))
                            ) where {𝘀<:PREC,𝗽<:PREC,𝗯<:BASE} = begin
            $TYPE(measurement(𝘀(amt(x).val), e), 𝗯)
        end
        (::Type{$TYPE{𝘀,MM}})(x::$TYPE{𝗽,MM,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝗯<:BASE} = begin
            $TYPE(Measurement{𝘀}(amt(x).val), 𝗯)
        end
        # Type export
        export $TYPE
        # Type-stable wrapped amount obtaining function
        amt(x::$TYPE{𝗽,EX,SY}) where 𝗽<:PREC = x.amt::Quantity{𝗽,$𝑑SY,$𝑢SY}
        amt(x::$TYPE{𝗽,EX,DT}) where 𝗽<:PREC = x.amt::Quantity{𝗽,$𝑑DT,$𝑢DT}
        amt(x::$TYPE{𝗽,EX,MA}) where 𝗽<:PREC = x.amt::Quantity{𝗽,$𝑑MA,$𝑢MA}
        amt(x::$TYPE{𝗽,EX,MO}) where 𝗽<:PREC = x.amt::Quantity{𝗽,$𝑑MO,$𝑢MO}
        amt(x::$TYPE{𝗽,MM,SY}) where 𝗽<:PREC = x.amt::Quantity{Measurement{𝗽},$𝑑SY,$𝑢SY}
        amt(x::$TYPE{𝗽,MM,DT}) where 𝗽<:PREC = x.amt::Quantity{Measurement{𝗽},$𝑑DT,$𝑢DT}
        amt(x::$TYPE{𝗽,MM,MA}) where 𝗽<:PREC = x.amt::Quantity{Measurement{𝗽},$𝑑MA,$𝑢MA}
        amt(x::$TYPE{𝗽,MM,MO}) where 𝗽<:PREC = x.amt::Quantity{Measurement{𝗽},$𝑑MO,$𝑢MO}
        # Type-specific functions
        deco(x::$TYPE{𝗽,𝘅,SY} where {𝗽,𝘅}) = Symbol($𝑠SY)
        deco(x::$TYPE{𝗽,𝘅,DT} where {𝗽,𝘅}) = Symbol($𝑠DT)
        deco(x::$TYPE{𝗽,𝘅,MA} where {𝗽,𝘅}) = Symbol($𝑠MA)
        deco(x::$TYPE{𝗽,𝘅,MO} where {𝗽,𝘅}) = Symbol($𝑠MO)
        ppu(x::$TYPE{𝗽,𝘅,SY} where {𝗽,𝘅}) = $USTR
        ppu(x::$TYPE{𝗽,𝘅,DT} where {𝗽,𝘅}) = $USTR * "/s"
        ppu(x::$TYPE{𝗽,𝘅,MA} where {𝗽,𝘅}) = $USTR * "/kg"
        ppu(x::$TYPE{𝗽,𝘅,MO} where {𝗽,𝘅}) = $USTR * "/kmol"
        # Indirect construction from plain
        $FNAM(x::plnF, b::Type{𝗯}=DEF[:IB]) where 𝗯<:BASE = $TYPE(x, b)
        $FNAM(x::REAL, b::Type{𝗯}=DEF[:IB]) where 𝗯<:BASE = $TYPE(float(x), b)
        # Indirect construction from quantity
        $FNAM(x::Union{UATY{𝗽,$𝑑SY},UATY{𝗽,$𝑑DT},
                       UATY{𝗽,$𝑑MA},UATY{𝗽,$𝑑MO}}) where 𝗽<:PREC = begin
            $TYPE(x)
        end
        $FNAM(x::Union{uniR{𝗽,$𝑑SY},uniR{𝗽,$𝑑DT},
                       uniR{𝗽,$𝑑MA},uniR{𝗽,$𝑑MO}}) where 𝗽<:REAL = begin
            $TYPE(float(x.val) * unit(x))
        end
        export $FNAM
        # Conversions - Change of base is _not_ a conversion
        # Same {EXAC,BASE}, {PREC}- conversion
        convert(::Type{$TYPE{𝘀,𝘅,𝗯}},
                y::$TYPE{𝗽,𝘅,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE} = begin
            $TYPE{promote_type(𝘀,𝗽),𝘅}(y)
        end
        # Same {BASE}, {PREC,EXAC}- conversion
        convert(::Type{$TYPE{𝘀,𝘆,𝗯}},
                y::$TYPE{𝗽,𝘅,𝗯}) where {𝘀<:PREC,𝗽<:PREC,𝘆<:EXAC,𝘅<:EXAC,𝗯<:BASE} = begin
            $TYPE{promote_type(𝘀,𝗽),promote_type(𝘆,𝘅)}(y)
        end
        # Promotion rules
        promote_rule(::Type{$TYPE{𝘀,𝘆,𝗯}},
                     ::Type{$TYPE{𝗽,𝘅,𝗯}}) where {𝘀<:PREC,𝗽<:PREC,
                                                  𝘆<:EXAC,𝘅<:EXAC,𝗯<:BASE} = begin
            $TYPE{promote_type(𝘀,𝗽),promote_type(𝘆,𝘅),𝗯}
        end
        # same-type sum,sub with Unitful promotion
        +(x::$TYPE{𝘀,𝘆,𝗯}, y::$TYPE{𝗽,𝘅,𝗯}) where {𝘀<:PREC,𝗽<:PREC,
                                                   𝘆<:EXAC,𝘅<:EXAC,
                                                   𝗯<:BASE} = begin
            $TYPE(+(amt(x), amt(y)))
        end
        -(x::$TYPE{𝘀,𝘆,𝗯}, y::$TYPE{𝗽,𝘅,𝗯}) where {𝘀<:PREC,𝗽<:PREC,
                                                   𝘆<:EXAC,𝘅<:EXAC,
                                                   𝗯<:BASE} = begin
            $TYPE(-(amt(x), amt(y)))
        end
        # scalar mul,div with Unitful promotion
        *(y::plnF{𝘀}, x::$TYPE{𝗽}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(*(amt(x), y))
        *(x::$TYPE{𝗽}, y::plnF{𝘀}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(*(amt(x), y))
        /(x::$TYPE{𝗽}, y::plnF{𝘀}) where {𝘀<:PREC,𝗽<:PREC} = $TYPE(/(amt(x), y))
        # Type-preserving scalar mul,div
        *(y::REAL, x::$TYPE{𝗽}) where 𝗽<:PREC = $TYPE(*(amt(x), 𝗽(y)))
        *(x::$TYPE{𝗽}, y::REAL) where 𝗽<:PREC = $TYPE(*(amt(x), 𝗽(y)))
        /(x::$TYPE{𝗽}, y::REAL) where 𝗽<:PREC = $TYPE(/(amt(x), 𝗽(y)))
    end
end

#----------------------------------------------------------------------------------------------#
#                           Thermodynamic Amount Group Declarations                            #
#----------------------------------------------------------------------------------------------#

# Mass / Mass fraction anomalous
mkBasAmt(:mAmt  , :BProperty, :m    , "m"   , u"kg"         , "kg"      , "mass"                , false , bsym=(:m , :ṁ , :mf, :M))
# Chemical amount / Molar fraction anomalous
mkBasAmt(:nAmt  , :BProperty, :N    , "N"   , u"kmol"       , "kmol"    , "chemical amount"     , false , bsym=(:N , :Ṅ , :n , :y))
# Gas constant / System constant anomalous
mkBasAmt(:RAmt  , :BProperty, :R    , "mR"  , u"kJ/K"       , "kJ/K"    , "gas constant"        , false , bsym=(:mR, :ṁR, :R , :R̄))
# Plank function anomalous
mkBasAmt(:rAmt  , :BProperty, :r    , "mr"  , u"kJ/K"       , "kJ/K"    , "Planck function"     , false , bsym=(:mr, :ṁr, :r , :r̄))

# Regular properties
mkBasAmt(:vAmt  , :BProperty, :v    , "V"   , u"m^3"        , "m³"      , "volume"              , false )
mkBasAmt(:uAmt  , :BProperty, :u    , "U"   , u"kJ"         , "kJ"      , "internal energy"     , false )
mkBasAmt(:hAmt  , :BProperty, :h    , "H"   , u"kJ"         , "kJ"      , "enthalpy"            , false )
mkBasAmt(:gAmt  , :BProperty, :g    , "G"   , u"kJ"         , "kJ"      , "Gibbs energy"        , false )
mkBasAmt(:aAmt  , :BProperty, :a    , "A"   , u"kJ"         , "kJ"      , "Helmholtz energy"    , false )
mkBasAmt(:eAmt  , :BProperty, :e    , "E"   , u"kJ"         , "kJ"      , "total energy"        , false )
mkBasAmt(:ekAmt , :BProperty, :ek   , "Ek"  , u"kJ"         , "kJ"      , "kinetic energy"      , false )
mkBasAmt(:epAmt , :BProperty, :ep   , "Ep"  , u"kJ"         , "kJ"      , "potential energy"    , false )
mkBasAmt(:sAmt  , :BProperty, :s    , "S"   , u"kJ/K"       , "kJ/K"    , "entropy"             , false )
mkBasAmt(:cpAmt , :BProperty, :cp   , "Cp"  , u"kJ/K"       , "kJ/K"    , "iso-P specific heat" , false )
mkBasAmt(:cvAmt , :BProperty, :cv   , "Cv"  , u"kJ/K"       , "kJ/K"    , "iso-v specific heat" , false )
mkBasAmt(:jAmt  , :BProperty, :j    , "J"   , u"kJ/K"       , "kJ/K"    , "Massieu function"    , false )

# Regular interactions
mkBasAmt(:qAmt  , :BInteract, :q    , "Q"   , u"kJ"         , "kJ"      , "heat"                , false )
mkBasAmt(:wAmt  , :BInteract, :w    , "W"   , u"kJ"         , "kJ"      , "work"                , false )
mkBasAmt(:ΔeAmt , :BInteract, :Δe   , "E"   , u"kJ"         , "kJ"      , "energy variation"    , true  )
mkBasAmt(:ΔsAmt , :BInteract, :Δs   , "S"   , u"kJ/K"       , "kJ/K"    , "entropy variation"   , true  )

#----------------------------------------------------------------------------------------------#
#                                      AMOUNT Type Unions                                      #
#----------------------------------------------------------------------------------------------#

# Unions of amounts of like units and thermodynamic classification, for same-unit operations

# --- energy
"""
`ENERGYP{𝗽,𝘅,𝗯} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}`\n
Energy property type union.
"""
ENERGYP{𝗽,𝘅,𝗯} = Union{uAmt{𝗽,𝘅,𝗯},hAmt{𝗽,𝘅,𝗯},
                       gAmt{𝗽,𝘅,𝗯},aAmt{𝗽,𝘅,𝗯},
                       eAmt{𝗽,𝘅,𝗯},ekAmt{𝗽,𝘅,𝗯},
                       epAmt{𝗽,𝘅,𝗯}} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}

"""
`ENERGYI{𝗽,𝘅,𝗯} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}`\n
Energy interaction type union.
"""
ENERGYI{𝗽,𝘅,𝗯} = Union{qAmt{𝗽,𝘅,𝗯},wAmt{𝗽,𝘅,𝗯},
                       ΔeAmt{𝗽,𝘅,𝗯}} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}

"""
`ENERGYA{𝗽,𝘅,𝗯} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}`\n
Energy amount type union.
"""
ENERGYA{𝗽,𝘅,𝗯} = Union{ENERGYP{𝗽,𝘅,𝗯},ENERGYI{𝗽,𝘅,𝗯}} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}

# --- entropy
"""
`NTROPYP{𝗽,𝘅,𝗯} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}`\n
Entropy property type union.
"""
NTROPYP{𝗽,𝘅,𝗯} = Union{RAmt{𝗽,𝘅,𝗯},rAmt{𝗽,𝘅,𝗯},
                       sAmt{𝗽,𝘅,𝗯},jAmt{𝗽,𝘅,𝗯},
                       cpAmt{𝗽,𝘅,𝗯},cvAmt{𝗽,𝘅,𝗯}} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}

"""
`NTROPYI{𝗽,𝘅,𝗯} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}`\n
Entropy interaction type union.
"""
NTROPYI{𝗽,𝘅,𝗯} = Union{ΔsAmt{𝗽,𝘅,𝗯}} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}

"""
`NTROPYA{𝗽,𝘅,𝗯} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}`\n
Entropy amount type union.
"""
NTROPYA{𝗽,𝘅,𝗯} = Union{NTROPYP{𝗽,𝘅,𝗯},NTROPYI{𝗽,𝘅,𝗯}} where {𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE}

# --- velocity
"""
`VELOCYP{𝗽,𝘅} where {𝗽<:PREC,𝘅<:EXAC}`\n
Velocity property type union.
"""
VELOCYP{𝗽,𝘅,𝗯} = Union{VELO{𝗽,𝘅},SPEE{𝗽,𝘅}} where {𝗽<:PREC,𝘅<:EXAC}

#----------------------------------------------------------------------------------------------#
#                                       Pretty Printing                                        #
#----------------------------------------------------------------------------------------------#

import Base: show
import Formatting: sprintf1

# Auxiliar method
function subscript(x::Int)
    asSub(c::Char) = Char(Int(c) - Int('0') + Int('₀'))
    map(asSub, "$(x)")
end

# Precision decoration
pDeco(::Type{Float16})  = DEF[:showPrec] ? subscript(16) : ""
pDeco(::Type{Float32})  = DEF[:showPrec] ? subscript(32) : ""
pDeco(::Type{Float64})  = DEF[:showPrec] ? subscript(64) : ""
pDeco(::Type{BigFloat}) = DEF[:showPrec] ? subscript(precision(BigFloat)) : ""

# Custom printing
Base.show(io::IO, x::AMOUNTS{𝗽,EX}) where 𝗽<:PREC = begin
    if DEF[:pprint]
        print(io,
            "$(string(deco(x)))$(pDeco(𝗽)): ",
            sprintf1("%.$(DEF[:showSigD])g", amt(x).val),
            " ", ppu(x))
    else
        Base.show_default(io, x)
    end
end

Base.show(io::IO, x::AMOUNTS{𝗽,MM}) where 𝗽<:PREC = begin
    if DEF[:pprint]
        print(io,
            "$(string(deco(x)))$(pDeco(𝗽)): (",
            sprintf1("%.$(DEF[:showSigD])g", amt(x).val.val),
            " ± ",
            sprintf1("%.2g", amt(x).val.err),
            ") ", ppu(x))
    else
        Base.show_default(io, x)
    end
end


text:
# # Matrix Formalism
#
# In this example we will consider the matrix formalism approach for a geometry of cylinders.

# LSP indexing solution                                                          #src
# https://github.com/julia-vscode/julia-vscode/issues/800#issuecomment-650085983 #src
if isdefined(@__MODULE__, :LanguageServer)                                       #src
    include("../src/SpinDoctor.jl")                                              #src
    using .SpinDoctor                                                            #src
end                                                                              #src

using SpinDoctor
using LinearAlgebra

if haskey(ENV, "GITHUB_ACTIONS")
    using CairoMakie
else
    using GLMakie
end

setup = CylinderSetup(;
    name = "Slice",
    ncell = 3,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    height = 1.0,
    bend = 0.0,
    twist = 0.0,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.5,
)

# We also define coefficients for the different cell compartments `:in` (axon), `:out`
# (myelin), and `:ecs` (ECS).

coeffs = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)

# We then proceed to build the geometry and finite element mesh.

mesh, = create_geometry(setup; recreate = true)
plot_mesh(mesh)

# The mesh looks good, so we may then proceed to assemble the biological model and the
# associated finite element matrices.

model = Model(; mesh, coeffs...);
matrices = assemble_matrices(model);

# We may also compute some useful quantities, including a scalar diffusion coefficient from
# the diffusion tensors.

volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)

# The eigenfunctions of the diffusive part of the Bloch-Torrey operator forms a good basis
# for the finite element function space. The basis may be truncated at a certain level, thus
# reducing the number of degrees of freedom. We here opt for 400 eigenfunctions.

laplace = Laplace(; model, matrices, neig_max = 400)
lap_eig = solve(laplace)

# The resulting eigenvalues may be represented as length scales, describing the wavelength
# of the eigenfunctions.

length_scales = eig2length.(lap_eig.values, D_avg)

# We may also further truncate the eigenfunction basis, if we are satisfied skipping
# features below a threshold length scale of 3 micrometers.

length_scale = 3
λ_max = length2eig(length_scale, D_avg)
lap_eig = limit_lengthscale(lap_eig, λ_max)

# Each of the resulting eigenfunctions is represented in the same way as the initial
# magnetization field `ρ`.

ncompartment, nboundary = size(mesh.facets)
fig = Figure()
for i = 1:3, j = 1:4
    ieig = 6(i - 1) + j
    ϕ_cmpts = split_field(mesh, lap_eig.funcs[:, ieig])
    ax = Axis3(fig[i, j]; title = "n = $ieig, ℓ = $(length_scales[ieig])", aspect = :data)
    nboundary = size(mesh.facets, 2)
    scene = nothing
    first = true
    for icmpt = 1:ncompartment, iboundary = 1:nboundary
        facets = mesh.facets[icmpt, iboundary]
        points = mesh.points[icmpt]
        mesh!(ax, points', facets', color = ϕ_cmpts[icmpt], shading = false)
    end
end
fig

# We observe that the first functions have large features, while the higher-index functions
# have more rapidly varying features. We may now choose a gradient and compute the
# projection of magnetization field onto the truncated basis.

dir = [1.0, 0.0, 0.0]
profile = CosOGSE(5000.0, 5000.0, 2)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)

# The matrix formalism problem is solved in the same way as the [`BTPDE`](@ref). The time
# profile is approximated on 500 points, since it is non-constant.

mf = MatrixFormalism(; model, matrices, lap_eig)
ξ = solve(mf, gradient; ninterval = 500)

# The resulting magnetization field may be plotted.

plot_field(model.mesh, ξ)

ext: jl


text:
import LinearAlgebra: adjoint, transpose
adjoint(Inv::Inverse) = inverse(adjoint(Inv.parent))
adjoint(Inv::PseudoInverse) = pseudoinverse(adjoint(Inv.parent))

transpose(Inv::Inverse) = inverse(transpose(Inv.parent))
transpose(Inv::PseudoInverse) = pseudoinverse(transpose(Inv.parent))

import LinearAlgebra: ishermitian, issymmetric
ishermitian(Inv::AbstractInverse) = ishermitian(Inv.parent)
issymmetric(Inv::AbstractInverse) = issymmetric(Inv.parent)

import LinearAlgebra: diag
# IDEA: allows for stochastic approximation:
# A Probing Method for Cοmputing the Diagonal of the Matrix Inverse
diag(Inv::AbstractInverse) = diag(AbstractMatrix(Inv))

#################### Basic multiplication and division #########################
import LinearAlgebra: *, /, \

*(L::AbstractInverse, B::Number) = pseudoinverse(L.parent / B) # returns AbstractInverse
*(L::AbstractInverse{<:Any, <:Number}, B::Number) = L.parent \ B
*(L::AbstractInverse, B::AbstractVector) = L.parent \ B
*(L::AbstractInverse, B::AbstractMatrix) = L.parent \ B
*(L::AbstractInverse, B::Adjoint{Any, <:AbstractMatrix}) = L.parent \ B

*(B::Number, L::AbstractInverse) = pseudoinverse(B \ L.parent) # returns AbstractInverse
*(B::Number, L::AbstractInverse{<:Any, <:Number}) = B / L.parent
*(B::Adjoint{<:Any, <:AbstractVector}, L::AbstractInverse) = B / L.parent
*(B::Adjoint{<:Any, <:AbstractMatrix}, L::AbstractInverse) = B / L.parent
*(B::AbstractMatrix, L::AbstractInverse) = B / L.parent

\(L::AbstractInverse, B::AbstractVector) = L.parent * B
\(L::AbstractInverse, B::AbstractMatrix) = L.parent * B
\(L::AbstractInverse, B::Adjoint{<:Any, <:AbstractMatrix}) = L.parent * B

/(B::Adjoint{<:Any, <:AbstractVector}, L::AbstractInverse) = B * L.parent
/(B::Adjoint{<:Any, <:AbstractMatrix}, L::AbstractInverse) = B * L.parent
/(B::AbstractMatrix, L::AbstractInverse) = B * L.parent

# scaling of AbstractInverse by Number
/(L::AbstractInverse, B::Number) = pseudoinverse(L.parent * B)
\(B::Number, L::AbstractInverse) = pseudoinverse(B * L.parent)

# *(L1::Inverse, L2::Inverse) =  Inverse(L1.parent * L2.parent) IDEA: LazyMatrixProduct to avoid O(n^3) multiply
# IDEA: could have check for L.parent ≡ B in multiply, to return identity with O(1) operations
##################### in-place multiplication and solving ######################
# IDEA: add rdiv!, ldiv! with Number types
# TODO: add further tests for mul!, and div! methods (e.g. involving scalar)
import LinearAlgebra: ldiv!, rdiv!, mul!, rmul!, lmul!
function ldiv!(Y::AbstractVecOrMat, A::AbstractInverse, B::AbstractVecOrMat)
	mul!(Y, A.parent, B)
end
# left multiplying with inverse
mul!(Y, A::AbstractInverse, B) = ldiv!(Y, A.parent, B) # 5 arg?
function mul!(Y, A::AbstractInverse, B, α::Real)
	mul!(Y, A, B)
	@. Y *= α
end
function mul!(Y, A::AbstractInverse, B, α::Real, β::Real)
	A.parent isa Matrix && throw("in place mul! only works if this only works if A.parent is Factorization OR a special matrix type like Diagonal, Bidiagonal, etc.")
	Z = copy(Y) # IDEA: pre-allocate somewhere?
	mul!(Y, A, B)
	@. Y = α*Y + β*Z
	return Y
end

# right multiplying with inverse
function mul!(Y, A, B::AbstractInverse)
	copy!(Y, A)
	rdiv!(Y, B.parent)
end
function mul!(Y, A, B::AbstractInverse, α::Real)
	mul!(Y, A, B)
	@. Y *= α
end
function mul!(Y, A, B::AbstractInverse, α::Real, β::Real)
	Z = copy(Y) # IDEA: pre-allocate somewhere?
	mul!(Y, A, B)
	if !(α == one(α) && β == zero(β))
		@. Y = α*Y + β*Z
	end
	return Y
end

# A \ b in place, overwriting B
lmul!(A::AbstractInverse, B) = ldiv!(A.parent, B) # these are usuall only defined for numbers
rmul!(A, B::AbstractInverse) = rdiv!(A, B.parent)

function ldiv!(A::AbstractInverse, B)
	Y = zero(B)
	mul!(Y, A.parent, B)
end
function rdiv!(A, B::AbstractInverse)
	Y = zero(A)
	mul!(Y, A, B.parent)
end

############################# ternary dot product ##############################
LinearAlgebra.dot(x, A::AbstractInverse, y) = dot(x, A*y)


text:
using ArgParse
parse_setting = ArgParseSettings()
@add_arg_table parse_setting begin
    "--io"
        help = "IO mode, `socket`, `file` default is `file`. "
        arg_type = String
        default = "file"
    "--host"
        help = "another option with an argument"
        arg_type = String
        default = "127.0.0.1"
    "--port"
        help = "another option with an argument"
        arg_type = Int
        default = 6543
    "--clienthost"
        help = "another option with an argument"
        arg_type = String
        default = "127.0.0.1"
    "--clientport"
        help = "another option with an argument"
        arg_type = Int
        default = 6544
    "--keep_open"
        help = "whether shutdown calculation after finish"
        action = :store_true
    "input_file"
        help = "a positional argument"
        required = false
end

args = parse_args(ARGS, parse_setting, as_symbols=true)

# Load MPI
import MPI
MPI.Init()
comm = MPI.COMM_WORLD
mpi_kwargs = Dict(:MPI => MPI, :comm => comm, :root => 0,
                  :nprc => MPI.Comm_size(comm), :rank => MPI.Comm_rank(comm),
                  :mFloat64 => MPI.Datatype(Float64),
                  :mUInt8 => MPI.Datatype(UInt8))

# Merge args and mpi
merge!(args, mpi_kwargs)

# Load TAPS IO module
if args[:io] == "file"
    IO = include("./file.jl")
    io = getfield(IO, :FileIO)(;args...)
elseif args[:io] == "socket"
    include("./socket.jl")
    # io = getfield(IO, :SocketIO)(;args...)
    io = SocketIO(;args...)
end

run_taps(io)
MPI.Finalize()

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
molecular_weight(model::EoSModel,z = @SVector [1.]) = 0.001*mapreduce(+,*,mw(model),z)
mw(model::EoSModel) = paramvals(model.params.Mw)

include("initial_guess.jl")
include("differentials.jl")
include("vt.jl")
include("property_solvers/property_solvers.jl")
include("pt.jl")
include("unitful_base.jl")
include("unitful_methods.jl")


text:
@testset "Poisson 2D" begin

    eta1_min = 0.0
    eta1_max = 2π
    eta2_min = 0.0
    eta2_max = 2π

    nc_eta1 = 128
    nc_eta2 = 128

    ex = zeros(Float64, nc_eta1, nc_eta2)
    ey = zeros(Float64, nc_eta1, nc_eta2)

    rho = zeros(Float64, nc_eta1, nc_eta2)
    phi = zeros(Float64, nc_eta1, nc_eta2)

    ex_exact = zeros(nc_eta1, nc_eta2)
    ey_exact = zeros(nc_eta1, nc_eta2)
    phi_exact = zeros(nc_eta1, nc_eta2)

    grid = TwoDGrid(eta1_min, eta1_max, nc_eta1, eta2_min, eta2_max, nc_eta2)

    poisson = TwoDPoissonPeriodic(grid)

    mode = 2
    for i = 1:nc_eta1, j = 1:nc_eta2
        x1 = eta1_min + (i - 1) * (eta1_max - eta1_min) / nc_eta1
        x2 = eta2_min + (j - 1) * (eta2_max - eta2_min) / nc_eta2
        phi_exact[i, j] = -mode * sin(mode * x1) * cos(mode * x2)
        ex_exact[i, j] = mode^2 * cos(mode * x1) * cos(mode * x2)
        ey_exact[i, j] = -mode^2 * sin(mode * x1) * sin(mode * x2)
        rho[i, j] = -2 * mode^3 * sin(mode * x1) * cos(mode * x2)
    end

    rhs = copy(rho)
    solve!(phi, poisson, rhs)
    @test phi_exact ≈ phi

    rhs = copy(rho)
    solve!(ex, ey, poisson, rhs)
    @test ex_exact ≈ ex
    @test ey_exact ≈ ey

    x1_min = 0
    x1_max = 1

    x2_min = 0
    x2_max = 1

    nc_x1 = 32
    nc_x2 = 64

    phi = zeros(nc_x1, nc_x2)
    e1 = zeros(nc_x1, nc_x2)
    e2 = zeros(nc_x1, nc_x2)
    rho = ones(nc_x1, nc_x2)

    grid = TwoDGrid(x1_min, x1_max, nc_x1, x2_min, x2_max, nc_x2)

    poisson = TwoDPoissonPeriodic(grid)

    solve!(phi, poisson, rho)

    solve!(e1, e2, poisson, rho)

    nx, ny = 32, 32
    xmin, xmax = 0.0, 2π
    ymin, ymax = 0.0, 2π
    degree_smoother = 3
    n_particles = 1000

    particles = ParticleGroup{2,2}(n_particles, charge = 1.0, mass = 1.0, n_weights = 1)

    grid = TwoDGrid(xmin, xmax, nx, ymin, ymax, ny)

    poisson = TwoDPoissonPeriodic(grid)

    # Initialize the kernel smoother
    kernel_smoother = ParticleMeshCoupling2D(particles, grid, degree_smoother, :collocation)

    @test true

    # Initialize the PIC field solver
    pic_poisson = (poisson, kernel_smoother)

    @test true

end


text:
module HydrostaticBoussinesq

export HydrostaticBoussinesqModel, AbstractHydrostaticBoussinesqProblem, OceanDGModel,
       LinearHBModel, calculate_dt

using StaticArrays
using LinearAlgebra: I, dot, Diagonal
using ..VariableTemplates
using ..MPIStateArrays
using ..DGmethods: init_ode_state
using ..PlanetParameters: grav
using ..Mesh.Filters: CutoffFilter, apply!, ExponentialFilter
using ..Mesh.Grids: polynomialorder, VerticalDirection, HorizontalDirection, min_node_distance

using ..DGmethods.NumericalFluxes: Rusanov, CentralNumericalFluxGradient,
                                   CentralNumericalFluxDiffusive,
                                   CentralNumericalFluxNonDiffusive

import ..DGmethods.NumericalFluxes: update_penalty!, numerical_flux_diffusive!,
                                    NumericalFluxNonDiffusive

import ..DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                    vars_diffusive, flux_nondiffusive!,
                    flux_diffusive!, source!, wavespeed,
                    boundary_state!, update_aux!, update_aux_diffusive!,
                    gradvariables!, init_aux!, init_state!,
                    LocalGeometry, DGModel, nodal_update_aux!, diffusive!,
                    copy_stack_field_down!, create_state, calculate_dt,
                    vars_integrals, vars_reverse_integrals,
                    indefinite_stack_integral!,
                    reverse_indefinite_stack_integral!,
                    integral_load_aux!, integral_set_aux!,
                    reverse_integral_load_aux!,
                    reverse_integral_set_aux!

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
∘(a::SVector, b::SVector) = StaticArrays.dot(a, b)

abstract type AbstractHydrostaticBoussinesqProblem end

"""
    HydrostaticBoussinesqModel <: BalanceLaw

A `BalanceLaw` for ocean modeling.

write out the equations here

# Usage

    HydrostaticBoussinesqModel(problem)

"""
struct HydrostaticBoussinesqModel{P,T} <: BalanceLaw
  problem::P
  cʰ::T
  cᶻ::T
  αᵀ::T
  νʰ::T
  νᶻ::T
  κʰ::T
  κᶻ::T
  function HydrostaticBoussinesqModel{FT}(problem;
                                      cʰ = FT(0),     # m/s
                                      cᶻ = FT(0),     # m/s
                                      αᵀ = FT(2e-4),  # (m/s)^2 / K
                                      νʰ = FT(5e3),   # m^2 / s
                                      νᶻ = FT(5e-3),  # m^2 / s
                                      κʰ = FT(1e3),   # m^2 / s
                                      κᶻ = FT(1e-4),  # m^2 / s
                                      ) where {FT <: AbstractFloat}
    return new{typeof(problem),FT}(problem, cʰ, cᶻ, αᵀ, νʰ, νᶻ, κʰ, κᶻ)
  end
end
HBModel = HydrostaticBoussinesqModel

"""
    calculate_dt(grid, model::HBModel)

calculates the time step based on grid spacing and model parameters
takes minimum of gravity wave, diffusive, and viscous CFL
factor of 1000 in the diffusive CFL is to handle the convective adjustment

"""
function calculate_dt(grid, model::HBModel, Courant_number)
    minΔx = min_node_distance(grid, HorizontalDirection())
    minΔz = min_node_distance(grid, VerticalDirection())

    CFL_gravity = minΔx / model.cʰ
    CFL_diffusive = minΔz^2 / (1000 * model.κᶻ)
    CFL_viscous = minΔz^2 / model.νᶻ

    dt = 1//2 * minimum([CFL_gravity, CFL_diffusive, CFL_viscous])

    return dt
end

"""
    OceanDGModel()

helper function to add required filtering
not used in the Driver+Config setup
"""
function OceanDGModel(bl::HBModel, grid, numfluxnondiff, numfluxdiff,
                      gradnumflux; kwargs...)
  vert_filter = CutoffFilter(grid, polynomialorder(grid)-1)
  exp_filter  = ExponentialFilter(grid, 1, 8)

  modeldata = (vert_filter = vert_filter, exp_filter=exp_filter)

  return DGModel(bl, grid, numfluxnondiff, numfluxdiff, gradnumflux;
                 kwargs..., modeldata=modeldata)
end

"""
    vars_state(::HBModel)

prognostic variables evolved forward in time

u = (u,v) = (zonal velocity, meridional velocity)
η = sea surface height
θ = temperature
"""
# If this order is changed check the filter usage!
function vars_state(m::HBModel, T)
  @vars begin
    u::SVector{2, T}
    η::T # real a 2-D variable TODO: should be 2D
    θ::T
  end
end

"""
    init_state!(::HBModel)

sets the initial value for state variables
dispatches to ocean_init_state! which is defined in a problem file such as SimpleBoxProblem.jl
"""
function ocean_init_state! end
function init_state!(m::HBModel, Q::Vars, A::Vars, coords, t)
  return ocean_init_state!(m.problem, Q, A, coords, t)
end

"""
    vars_aux(::HBModel)
helper variables for computation

first half is because there is no dedicated integral kernels
these variables are used to compute vertical integrals
w = vertical velocity
w_reverse = w but integrated in the opposite direction
wz0 = w at z = 0
pkin = bulk hydrostatic pressure contribution
pkin_reverse = pkin but integrated in the opposite direction

second half of these are fields that are used for computation
θʳ = relaxation value of the sea surface temperature
f = coriolis force
τ = wind stress
ν = vector of viscosities (zonal, meridional, vertical)
κ = vector of diffusivities (zonal, meridional, vertical)
"""
# If this order is changed check update_aux!
function vars_aux(m::HBModel, T)
  @vars begin
    w::T
    pkin::T         # ∫(-αᵀ θ)
    wz0::T          # w at z=0
    θʳ::T           # SST given    # TODO: Should be 2D
    f::T            # coriolis
    τ::T            # wind stress  # TODO: Should be 2D
    ν::SVector{3, T}
    κ::SVector{3, T}
  end
end

"""
    init_aux!(::HBModel)

sets the initial value for auxiliary variables (those that aren't related to vertical integrals)
dispatches to ocean_init_aux! which is defined in a problem file such as SimpleBoxProblem.jl
"""
function ocean_init_aux! end
function init_aux!(m::HBModel, A::Vars, geom::LocalGeometry)
  return ocean_init_aux!(m, m.problem, A, geom)
end

"""
    vars_gradient(::HBModel)

variables that you want to take a gradient of
these are just copies in our model
"""
function vars_gradient(m::HBModel, T)
  @vars begin
    u::SVector{2, T}
    θ::T
  end
end

"""
    gradvariables!(::HBModel)

copy u and θ to var_gradient
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `G`: array of gradient variables
- `Q`: array of state variables
- `A`: array of aux variables
- `t`: time, not used
"""
@inline function gradvariables!(m::HBModel, G::Vars, Q::Vars, A, t)
  G.u = Q.u
  G.θ = Q.θ

  return nothing
end

"""
    vars_diffusive(::HBModel)

the output of the gradient computations
once again just copies, we don't do any transforms or reductions
"""
function vars_diffusive(m::HBModel, T)
  @vars begin
    ∇u::SMatrix{3, 2, T, 6}
    ∇θ::SVector{3, T}
  end
end

"""
    diffusive!(::HBModel)

copy ∇u and ∇θ to var_diffusive
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `D`: array of diffusive variables
- `G`: array of gradient variables
- `Q`: array of state variables
- `A`: array of aux variables
- `t`: time, not used
"""
@inline function diffusive!(m::HBModel, D::Vars, G::Grad, Q::Vars,
                            A::Vars, t)
  D.∇u = G.u
  D.∇θ = G.θ

  return nothing
end

"""
    vars_integral(::HBModel)

location to store integrands for bottom up integrals
∇hu = the horizontal divegence of u, e.g. dw/dz
"""
function vars_integrals(m::HBModel, T)
  @vars begin
    ∇hu::T
  end
end

"""
    integral_load_aux!(::HBModel)

copy w to var_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
I -> array of integrand variables
Q -> array of state variables
A -> array of aux variables
"""
@inline function integral_load_aux!(m::HBModel, I::Vars, Q::Vars, A::Vars)
  I.∇hu = A.w # borrow the w value from A...

  return nothing
end

"""
    integral_set_aux!(::HBModel)

copy integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function integral_set_aux!(m::HBModel, A::Vars, I::Vars)
  A.w = I.∇hu

  return nothing
end

"""
    vars_reverse_integral(::HBModel)

location to store integrands for top down integrals
αᵀθ = density perturbation
"""
function vars_reverse_integrals(m::HBModel, T)
  @vars begin
    αᵀθ::T
  end
end

"""
    reverse_integral_load_aux!(::HBModel)

copy αᵀθ to var_reverse_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
I -> array of integrand variables
A -> array of aux variables
"""
@inline function reverse_integral_load_aux!(m::HBModel, I::Vars, Q::Vars, A::Vars)
  I.αᵀθ = A.pkin

  return nothing
end

"""
    reverse_integral_set_aux!(::HBModel)

copy reverse integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function reverse_integral_set_aux!(m::HBModel, A::Vars, I::Vars)
  A.pkin = I.αᵀθ

  return nothing
end

"""
    flux_nondiffusive!(::HBModel)

calculates the hyperbolic flux contribution to state variables
this computation is done pointwise at each nodal point

# arguments:
m -> model in this case HBModel
F -> array of fluxes for each state variable
Q -> array of state variables
A -> array of aux variables
t -> time, not used

# computations
∂ᵗu = ∇∘(g*η + g∫αᵀθdz + v∘u)
∂ᵗθ = ∇∘(vθ) where v = (u,v,w)
"""
@inline function flux_nondiffusive!(m::HBModel, F::Grad, Q::Vars,
                                    A::Vars, t::Real)
  @inbounds begin
    u = Q.u # Horizontal components of velocity
    η = Q.η
    θ = Q.θ
    w = A.w   # vertical velocity
    pkin = A.pkin

    v = @SVector [u[1], u[2], w]
    Ih = @SMatrix [ 1 -0;
                   -0  1;
                   -0 -0]

    # ∇h • (g η)
    F.u += grav * η * Ih

    # ∇h • (- ∫(αᵀ θ))
    F.u += grav * pkin * Ih

    # ∇h • (v ⊗ u)
    # F.u += v * u'

    # ∇ • (u θ)
    F.θ += v * θ
  end

  return nothing
end

"""
    flux_diffusive!(::HBModel)

calculates the parabolic flux contribution to state variables
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `F`: array of fluxes for each state variable
- `Q`: array of state variables
- `D`: array of diff variables
- `A`: array of aux variables
- `t`: time, not used

# computations
∂ᵗu = -∇∘(ν∇u)
∂ᵗθ = -∇∘(κ∇θ)
"""
@inline function flux_diffusive!(m::HBModel, F::Grad, Q::Vars, D::Vars,
                                 A::Vars, t::Real)
  F.u -= Diagonal(A.ν) * D.∇u
  F.θ -= Diagonal(A.κ) * D.∇θ

  return nothing
end

"""
    source!(::HBModel)
    calculates the source term contribution to state variables
    this computation is done pointwise at each nodal point

    arguments:
    m -> model in this case HBModel
    F -> array of fluxes for each state variable
    Q -> array of state variables
    A -> array of aux variables
    t -> time, not used

    computations
    ∂ᵗu = -f×u
    ∂ᵗη = w|(z=0)
"""
@inline function source!(m::HBModel{P}, source::Vars, Q::Vars,
                         diffusive::Vars, A::Vars, t::Real) where P
  @inbounds begin
    u,v = Q.u # Horizontal components of velocity
    f = A.f
    wz0 = A.wz0

    # f × u
    source.u -= @SVector [-f * v, f * u]

    source.η += wz0
  end

  return nothing
end

"""
    wavespeed(::HBModel)

calculates the wavespeed for rusanov flux
"""
@inline wavespeed(m::HBModel, n⁻, _...) = abs(SVector(m.cʰ, m.cʰ, m.cᶻ)' * n⁻)

"""
    update_penalty(::HBModel)
    set Δη = 0 when computing numerical fluxes
"""
# We want not have jump penalties on η (since not a flux variable)
function update_penalty!(::Rusanov, ::HBModel, n⁻, λ, ΔQ::Vars,
                         Q⁻, A⁻, Q⁺, A⁺, t)
  ΔQ.η = -0

  return nothing
end

"""
    update_aux!(::HBModel)

    applies the vertical filter to the zonal and meridional velocities to preserve numerical incompressibility
    applies an exponential filter to θ to anti-alias the non-linear advective term

    doesn't actually touch the aux variables any more, but we need a better filter interface than this anyways
"""
function update_aux!(dg::DGModel, m::HBModel, Q::MPIStateArray, t::Real)
  MD = dg.modeldata

  # required to ensure that after integration velocity field is divergence free
  vert_filter = MD.vert_filter
  # Q[1] = u[1] = u, Q[2] = u[2] = v
  apply!(Q, (1, 2), dg.grid, vert_filter, VerticalDirection())

  exp_filter = MD.exp_filter
  # Q[4] = θ
  apply!(Q, (4,), dg.grid, exp_filter, VerticalDirection())

  return true
end

"""
    update_aux_diffusive!(::HBModel)

    ∇hu to w for integration
    performs integration for w and pkin (should be moved to its own integral kernels)
    copies down w and wz0 because we don't have 2D structures

    now for actual update aux stuff
    implements convective adjustment by bumping the vertical diffusivity up by a factor of 1000 if dθdz < 0
"""
function update_aux_diffusive!(dg::DGModel, m::HBModel, Q::MPIStateArray, t::Real)
  A  = dg.auxstate

  # store ∇ʰu as integrand for w
  # update vertical diffusivity for convective adjustment
  function f!(m::HBModel, Q, A, D, t)
    @inbounds begin
      A.w = -(D.∇u[1,1] + D.∇u[2,2])
      A.pkin = -m.αᵀ * Q.θ

      D.∇θ[3] < 0 ? A.κ = (m.κʰ, m.κʰ, 1000 * m.κᶻ) : A.κ = (m.κʰ, m.κʰ, m.κᶻ)
    end

    return nothing
  end
  nodal_update_aux!(f!, dg, m, Q, t; diffusive=true)

  # compute integrals for w and pkin
  indefinite_stack_integral!(dg, m, Q, A, t) # bottom -> top
  reverse_indefinite_stack_integral!(dg, m, Q, A, t) # top -> bottom

  # project w(z=0) down the stack
  # Need to be consistent with vars_aux
  # A[1] = w, A[3] = wz0
  copy_stack_field_down!(dg, m, A, 1, 3)

  return true
end

"""
    boundary_state!(nf, ::HBModel, Q⁺, A⁺, Q⁻, A⁻, bctype)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, m::HBModel, Q⁺::Vars, A⁺::Vars, n⁻,
                                 Q⁻::Vars, A⁻::Vars, bctype, t, _...)
  return ocean_boundary_state!(m, m.problem, bctype, nf, Q⁺, A⁺, n⁻, Q⁻, A⁻, t)
end

"""
    boundary_state!(nf, ::HBModel, Q⁺, D⁺, A⁺, Q⁻, D⁻, A⁻, bctype)

applies boundary conditions for the parabolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, m::HBModel,
                                 Q⁺::Vars, D⁺::Vars, A⁺::Vars,
                                 n⁻,
                                 Q⁻::Vars, D⁻::Vars, A⁻::Vars,
                                 bctype, t, _...)
  return ocean_boundary_state!(m, m.problem, bctype, nf, Q⁺, D⁺, A⁺, n⁻, Q⁻, D⁻, A⁻, t)
end

include("OceanBoundaryConditions.jl")
include("LinearHBModel.jl")

end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text: n = 922337203685477580900000000000

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
const colnames_style = Dict(:color=>:yellow,)
const rownames_style = Dict(:color=>:cyan,)

# Base.show
function Base.show(io::IO, colored::ColoredString)
    if Base.have_color && colored.style isa Dict
         printstyled(io, colored.data, color=colored.style[:color])
    else
         print(io, colored.data)
    end
end

function Base.show(io::IO, v::Vertical)
    print(io, v.data)
end

function Base.show(io::IO, dash::Dash)
    print(io, repeat(dash.data, dash.repeat))
end

function Base.show(io::IO, connector::Connector)
    print(io, connector.data)
end

function print_plate(io::IO, horizon::Horizontal, tablemode::TableMode, islast::Bool)
    for j=1:length(horizon)
        hr = horizon[j]
        print(io, hr)
    end
    !islast && println(io)
end

function print_plate(io::IO, linear::Linear, tablemode::TableMode, islast::Bool)
    isempty(linear) && return
    firstcell = first(linear)
    height = firstcell.height
    @inbounds for i=1:height
        for j=1:length(linear)
            cv = linear[j]
            if isa(cv, Vertical)
                print(io, cv)
            else
                cel = cv
                rows,cols = arraysize(cel.data)
                if rows >= i
                    for n=1:cols
                        print(io, cel.data[i,n])
                    end
                else
                    print(io, repeat(" ", cel.width))
                end
            end
        end
        if tablemode.style==:markdown && i==height
            !islast && println(io)
        else
            println(io)
        end
    end
end

function print_plates(io::IO, plates::PlateVector, tablemode::TableMode)
    len = length(plates)
    @inbounds for i=1:len
        plate = plates[i]
        print_plate(io, plate, tablemode, i==len)
    end
end

function Base.show(io::IO, mill::Mill)
    prepend_newline = display_style[:prepend_newline]
    prepend_newline && println(io)
    tablemode = display_style[:table_mode]
    plates = decking(mill, tablemode)
    print_plates(io, plates, tablemode)
end

# Base.isless
Base.isless(::String, ::Symbol) = false

# arraysize
arraysize(a::AbstractArray) = size(a[:,:])

# precell
precell(::Nothing, style=nothing) = DataCell([], 0, 0, style)

function precell(el::String, style=nothing)
    if occursin("\n", el)
        a = map(split(el, "\n")) do x
            width = textwidth(x)
            string(x, repeat(" ", width-textwidth(x)))
        end
        m,n = arraysize(a)
        DataCell(a, maximum(map(length, a)), m, style)
    else
        width = textwidth(el)
        DataCell([el], width, 1, style)
    end
end

function precell(el::AbstractArray, style=nothing)
    rows,cols = arraysize(el)
    if 0==rows
        DataCell([""], 0, 0, style)
    else
        prelines = similar(el, String)
        widths = zeros(Int, rows, cols)
        for i=1:rows
            for j=1:cols
                s = replace(string(el[i,j]), "\n" => " ")
                widths[i,j] = length(s)
                prelines[i,j] = s
            end
        end
        maxwidths = zeros(Int, cols)
        for j=1:cols
            maxwidths[j] = maximum(widths[:,j])
        end

        lines = AbstractArray{String}([])
        for i=1:rows
            line = AbstractArray{String}([])
            for (j,x) in enumerate(prelines[i,:])
                push!(line, lpad(x, maxwidths[j]))
            end
            push!(lines, join(line, " "))
        end
        DataCell(lines, sum(maxwidths) + (cols-1), rows, style)
    end
end

function precell(el::Any, style=nothing)
    DataCell([el], length(repr(el)), 1, style)
end

# postcell
function postcell(cell::DataCell, width::Int, height::Int, margin::Margin)
    data = cell.data
    rows,cols = arraysize(data)
    A = Array{ColoredString}(undef, rows, cols)
    @inbounds for i=1:rows
        for j=1:cols
            el = data[i,j]
            x = isa(el, String) ? el : repr(el)
            prep = repeat(" ", margin.leftside + width - textwidth(x))
            A[i,j] = ColoredString(string(prep, x, repeat(" ", margin.rightside)), cell.style)
        end
    end
    Cell(A, width+margin.leftside+margin.rightside, height)
end

function horizontal(maxwidths::Vector{Int}, cols::Int, margin::Margin, tablemode::TableMode; dash="-")
    corner_corner = string(tablemode.corner_corner)
    horizon = Horizontal{Union{Dash,Connector}}([])
    push!(horizon, Connector(corner_corner))
    for j=1:cols
        width = maxwidths[j]
        push!(horizon, Dash(dash, width + margin.leftside + margin.rightside))
        push!(horizon, Connector(corner_corner))
    end
    horizon
end

function marginal(option::Dict{Symbol,Any})
    Margin(1, 1) # fixed
end

function decking(mill::Mill, tablemode::TableMode)
    board = mill.board
    option = mill.option
    header_fillchar = string(tablemode.header_fillchar)

    rows,cols = arraysize(board)
    plates = PlateVector{Union{Linear,Horizontal}}([])
    if 0==rows
        return plates
    end

    widths = zeros(Int, rows+1, cols+1)
    heights = zeros(Int, rows+1, cols+1)
    # prepare ells
    preplates = Vector{Vector{DataCell}}([])
    headlinear = Vector{Union{DataCell,Vertical}}([])
    zerocolname = nothing
    if haskey(option, :zerocolname)
        zerocolname = option[:zerocolname]
    end
    cel = precell(zerocolname, colnames_style)
    widths[1,1] = cel.width
    heights[1,1] = cel.height
    push!(headlinear, cel)
    ols = Vector{DataCell}(undef, cols)
    for j=1:cols
        ols[j] = precell(j, colnames_style)
    end
    if haskey(option, :colnames)
        @inbounds for (j,name) in enumerate(option[:colnames])
            if cols >= j
                ols[j] = precell(name, colnames_style)
            end
        end
    end
    for j=1:cols
        cel = ols[j]
        widths[1,j+1] = cel.width
        heights[1,j+1] = cel.height
        push!(headlinear, cel)
    end
    push!(preplates, headlinear)

    prerows = Vector{DataCell}(undef, rows)
    for i=1:rows
        prerows[i] = precell(i, rownames_style)
    end
    if haskey(option, :rownames)
        @inbounds for (i,name) in enumerate(option[:rownames])
            if rows >= i
                prerows[i] = precell(name, rownames_style)
            end
        end
    end
    for i=1:rows
        linear = Vector{Union{DataCell,Vertical}}([])
        rownamecell = prerows[i]
        widths[1+i,1] = rownamecell.width
        heights[1+i,1] = rownamecell.height
        push!(linear, rownamecell)
        for j=1:cols
            cel = precell(board[i,j])
            widths[1+i,1+j] = cel.width
            heights[1+i,1+j] = cel.height
            push!(linear, cel)
        end
        push!(preplates, linear)
    end

    maxwidths = zeros(Int, cols+1)
    for j=1:cols+1
        maxwidths[j] = maximum(widths[:,j])
    end
    maxheights = zeros(Int, rows+1)
    for i=1:rows+1
        maxheights[i] = maximum(heights[i,:])
    end
    margin = marginal(option)
    if 0==rows
        cols = 1
    else
        rows += 1
        cols += 1
    end

    isgridmode = tablemode.style == :grid
    isgridmode && push!(plates, horizontal(maxwidths, cols, margin, tablemode, dash="-"))
    @inbounds for i=1:rows
        linear = Linear{Union{Cell,Vertical}}([])
        preplate = preplates[i]
        push!(linear, Vertical(maxheights[i]))
        maxheight = maxheights[i]
        for j=1:cols
            el = preplate[j]
            cel = postcell(el, maxwidths[j], maxheight, margin)
            if 1==j && cel.height > 1
                prep = [repeat(" ", cel.width) for x in 1:(cel.height-length(cel.data))]
                cel.data = vcat(prep, cel.data)
            end
            push!(linear, cel)
            push!(linear, Vertical(maxheight))
        end
        push!(plates, linear)
        (1==i || isgridmode) && push!(plates, horizontal(maxwidths, cols, margin, tablemode, dash = 1==i ? header_fillchar : "-"))
    end
    plates
end

# table
table(board::Any; options...) = Mill([board]; options...)
table(board::AbstractArray; options...) = Mill(board; options...)
table(board::Tuple; options...) = Mill(collect(board); options...)

function table(board::Dict; options...)
    dict = Dict{Symbol, Any}(options)
    if !haskey(dict, :colnames)
        dict[:zerocolname] = "KEY"
        dict[:colnames] = ["VALUE"]
    end
    rownames = sort(collect(keys(board)))
    dict[:rownames] = rownames
    Mill(getindex.(Ref(board), rownames); tuple(dict...)...)
end

table(board::T; options...) where {T<:Base.Enums.Enum} = table(typeof(board), options...)
function table(board::Type{T}; options...) where {T<:Base.Enums.Enum}
    dict = Dict{Symbol, Any}(options)
    enums = instances(board)
    dict[:rownames] = string.(enums)
    if !haskey(dict, :colnames)
        dict[:zerocolname] = string(board)
        dict[:colnames] = ["VALUE"]
    end
    Mill(collect(Int.(enums)); tuple(dict...)...)
end


text:
function flux(u, q, w, v, x, t, mu, eta)
    f = mu[1]*q;
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = (3*pi*pi)*sin(pi*x[1])*sin(pi*x[2])*sin(pi*x[3]);
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = mu[2];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[1]*n[1] + f[2]*n[2] + f[3]*n[3] + tau[1]*(u[1]-0.0);
    return fb;
end
function initu(x, mu, eta)
    u0 = 0.0;
    return u0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
function FODEdG!(λ::T, f::Function, u0::T, t::OffsetVector{T},
                 r::Integer, M::Integer,
                 store::Store{T}) where T <: AbstractFloat
    @argcheck r ≤ store.rmax
    @argcheck M ≤ store.Mmax
    N = length(t) - 1
    U = Vector{Vector{T}}(undef, N)
    G = coef_G(T, r)
    K = coef_K(T, r, r)
    Hn = OffsetArray{Matrix{T}}(undef, 0:N-1)
    for n = 0:N-1
        Hn[n] = Array{T}(undef, r, r)
    end
    τ, wτ = rule(store.legendre[M])
    Ψ = view(store.Ψ, 1:r, 1:M)
    legendre_polys!(Ψ, τ)
    Fn = Array{T}(undef, r)
    load_vector!(Fn, (t[0],t[1]), f, τ, wτ, Ψ)
    Ψ_at_minus1 = Vector{T}(undef, r)
    legendre_polys!(Ψ_at_minus1, -one(T))
    coef_Hn!(Hn, 1, 0, t, M, store)
    A = G + λ * Hn[0]
    b = copy(Fn) + Ψ_at_minus1 * u0
    Fact = LinearAlgebra.lu(A)
    U[1] = Fact \ b
    for n = 2:N
        coef_Hn!(Hn, n, n-1, t, M, store)
        load_vector!(Fn, (t[n-1],t[n]), f, τ, wτ, Ψ)
        fill!(b, zero(T))
        for ℓ = 1:n-1
            b .= b .+ ( Hn[n-ℓ] * U[ℓ] )
        end
        b .= Fn - λ * b + K * U[n-1]
        A = G + λ * Hn[0]
        Fact = LinearAlgebra.lu(A)
        U[n] = Fact \ b
    end
    return U
end

function FODEdG!(λ::T, tmax::T, f::Function, u0::T,
                 N::Integer, r::Integer, M::Integer,
                 store::Store{T}) where T <: AbstractFloat
    @argcheck r ≤ store.rmax
    @argcheck M ≤ store.Mmax
    G = coef_G(T, r)
    K = coef_K(T, r, r)
    H = coef_H_uniform!(r, N, M, store)
    U = Vector{Vector{T}}(undef, N)
    k = tmax / N
    α = store.α
    t_ = range(zero(T), tmax, length=N+1)
    t = OffsetArray(t_, 0:N)
    τ, wτ = rule(store.legendre[M])
    Ψ = view(store.Ψ, 1:r, 1:M)
    legendre_polys!(Ψ, τ)
    Fn = Array{T}(undef, r)
    load_vector!(Fn, (t[0],t[1]), f, τ, wτ, Ψ)
    Ψ_at_minus1 = Vector{T}(undef, r)
    legendre_polys!(Ψ_at_minus1, -one(T))
    A = G + ( λ * k^α ) * H[0]
    b = copy(Fn) + Ψ_at_minus1 * u0
    Fact = LinearAlgebra.lu(A)
    U[1] = Fact \ b
    for n = 2:N
        load_vector!(Fn, (t[n-1],t[n]), f, τ, wτ, Ψ)
        fill!(b, zero(T))
        for ℓ = 1:n-1
            b .= b .+ ( H[n-ℓ] * U[ℓ] )
        end
        b .= Fn - ( λ * k^α ) .* b + K * U[n-1]
        U[n] = Fact \ b
    end
    return t, U
end

function load_vector!(Fn::Vector{T}, In::Tuple{T,T}, f::Function,
                      τ::Vector{T}, wτ::Vector{T}, Ψ::AbstractMatrix{T}
                     ) where T <: AbstractFloat
    r, M = size(Ψ)
    tnm1, tn = In
    fill!(Fn, zero(T))
    for m = 1:M
        tm = ( (1-τ[m]) * tnm1 + (1+τ[m]) * tn ) / 2
        for i = 1:r
            Fn[i] += wτ[m] * f(tm) * Ψ[i,m]
        end
    end
    kn = tn - tnm1
    for i = 1:r
        Fn[i] *= kn / 2
    end
end

function graded_mesh(max_t::T, cutoff_t::T,
                     N_outer::Integer, mg::T) where { T <: AbstractFloat }
    k = ( max_t - cutoff_t ) / N_outer
    N_inner = ceil(Integer, -1/expm1(log1p(-k/cutoff_t)/mg))
    N = N_inner + N_outer
    t = OffsetVector{T}(undef, 0:N)
    t[0] = zero(T)
    for n = 1:N_inner
        t[n] = (n/N_inner)^mg * cutoff_t
    end
    for n = N_inner+1:N
        t[n] = cutoff_t + (n-N_inner)*k
    end
    return t
end

function ODEdG!(λ::T, tmax::T, f::Function, u0::T,
                N::Integer, r::Integer, M::Integer,
                ) where T <: AbstractFloat
    G = coef_G(T, r)
    K = coef_K(T, r, r)
    H = coef_H(T, r)
    U = Vector{Vector{T}}(undef, N)
    k = tmax / N
    t_ = range(zero(T), tmax, length=N+1)
    t = OffsetArray(t_, 0:N)
    τ, wτ = GaussQuadrature.legendre(M)
    Ψ = Matrix{T}(undef, r, M)
    legendre_polys!(Ψ, τ)
    Fn = Array{T}(undef, r)
    load_vector!(Fn, (t[0],t[1]), f, τ, wτ, Ψ)
    Ψ_at_minus1 = Vector{T}(undef, r)
    legendre_polys!(Ψ_at_minus1, -one(T))
    A = G + ( λ * k ) * H
    b = copy(Fn) + Ψ_at_minus1 * u0
    Fact = LinearAlgebra.lu(A)
    U[1] = Fact \ b
    for n = 2:N
        load_vector!(Fn, (t[n-1],t[n]), f, τ, wτ, Ψ)
        b .= Fn + K * U[n-1]
        U[n] = Fact \ b
    end
    return t, U
end


text:
int_rules_8_4 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.4*Trig*integral*functions =#
    @apply_utils Antiderivative(SinIntegral(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * SinIntegral(~a + ~b * ~x)) / ~b + cos(~a + ~b * ~x) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(CosIntegral(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * CosIntegral(~a + ~b * ~x)) / ~b - sin(~a + ~b * ~x) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(SinIntegral(~(b') * ~x) / ~x, ~x) => (1 / 2) * ~b * ~x * HypergeometricPFQ([1, 1, 1], [2, 2, 2], -im * ~b * ~x) + (1 / 2) * ~b * ~x * HypergeometricPFQ([1, 1, 1], [2, 2, 2], im * ~b * ~x) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative(CosIntegral(~(b') * ~x) / ~x, ~x) => (-1 / 2) * im * ~b * ~x * HypergeometricPFQ([1, 1, 1], [2, 2, 2], -im * ~b * ~x) + (1 / 2) * im * ~b * ~x * HypergeometricPFQ([1, 1, 1], [2, 2, 2], im * ~b * ~x) + EulerGamma * log(~x) + (1 / 2) * log(~b * ~x) ^ 2 <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * SinIntegral(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * SinIntegral(~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(((~c + ~d * ~x) ^ (~m + 1) * sin(~a + ~b * ~x)) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * CosIntegral(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * CosIntegral(~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(((~c + ~d * ~x) ^ (~m + 1) * cos(~a + ~b * ~x)) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(SinIntegral(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * SinIntegral(~a + ~b * ~x) ^ 2) / ~b - 2 * Antiderivative(sin(~a + ~b * ~x) * SinIntegral(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(CosIntegral(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * CosIntegral(~a + ~b * ~x) ^ 2) / ~b - 2 * Antiderivative(cos(~a + ~b * ~x) * CosIntegral(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * SinIntegral(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * SinIntegral(~b * ~x) ^ 2) / (~m + 1) - (2 / (~m + 1)) * Antiderivative((~x) ^ ~m * sin(~b * ~x) * SinIntegral(~b * ~x), ~x) <-- FreeQ(~b, ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * CosIntegral(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * CosIntegral(~b * ~x) ^ 2) / (~m + 1) - (2 / (~m + 1)) * Antiderivative((~x) ^ ~m * cos(~b * ~x) * CosIntegral(~b * ~x), ~x) <-- FreeQ(~b, ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * SinIntegral(~a + ~(b') * ~x) ^ 2, ~x) => (((~a + ~b * ~x) * (~c + ~d * ~x) ^ ~m * SinIntegral(~a + ~b * ~x) ^ 2) / (~b * (~m + 1)) - (2 / (~m + 1)) * Antiderivative((~c + ~d * ~x) ^ ~m * sin(~a + ~b * ~x) * SinIntegral(~a + ~b * ~x), ~x)) + (((~b * ~c - ~a * ~d) * ~m) / (~b * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * SinIntegral(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * CosIntegral(~a + ~(b') * ~x) ^ 2, ~x) => (((~a + ~b * ~x) * (~c + ~d * ~x) ^ ~m * CosIntegral(~a + ~b * ~x) ^ 2) / (~b * (~m + 1)) - (2 / (~m + 1)) * Antiderivative((~c + ~d * ~x) ^ ~m * cos(~a + ~b * ~x) * CosIntegral(~a + ~b * ~x), ~x)) + (((~b * ~c - ~a * ~d) * ~m) / (~b * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * CosIntegral(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    #= Antiderivative((~x)^(~m')*SinIntegral((~a)+(~b')*(~x))^2,~x) :=  b*x^(m+2)*SinIntegral(a+b*x)^2/(a*(m+1)) + x^(m+1)*SinIntegral(a+b*x)^2/(m+1) - 2*b/(a*(m+1))*Antiderivative(x^(m+1)*sin(a+b*x)*SinIntegral(a+b*x),x) - b*(m+2)/(a*(m+1))*Antiderivative(x^(m+1)*SinIntegral(a+b*x)^2,x) <--  FreeQ([a,b],x) && ILtQ(m,-2) =#
    #= Antiderivative((~x)^(~m')*CosIntegral((~a)+(~b')*(~x))^2,~x) :=  b*x^(m+2)*CosIntegral(a+b*x)^2/(a*(m+1)) + x^(m+1)*CosIntegral(a+b*x)^2/(m+1) - 2*b/(a*(m+1))*Antiderivative(x^(m+1)*cos(a+b*x)*CosIntegral(a+b*x),x) - b*(m+2)/(a*(m+1))*Antiderivative(x^(m+1)*CosIntegral(a+b*x)^2,x) <--  FreeQ([a,b],x) && ILtQ(m,-2) =#
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * SinIntegral(~(c') + ~(d') * ~x), ~x) => (-(cos(~a + ~b * ~x)) * SinIntegral(~c + ~d * ~x)) / ~b + (~d / ~b) * Antiderivative((cos(~a + ~b * ~x) * sin(~c + ~d * ~x)) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * CosIntegral(~(c') + ~(d') * ~x), ~x) => (sin(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x)) / ~b - (~d / ~b) * Antiderivative((sin(~a + ~b * ~x) * cos(~c + ~d * ~x)) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) * SinIntegral(~(c') + ~(d') * ~x), ~x) => (-((~e + ~f * ~x) ^ ~m) * cos(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x)) / ~b + (~d / ~b) * Antiderivative(((~e + ~f * ~x) ^ ~m * cos(~a + ~b * ~x) * sin(~c + ~d * ~x)) / (~c + ~d * ~x), ~x) + ((~f * ~m) / ~b) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * cos(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x) * CosIntegral(~(c') + ~(d') * ~x), ~x) => (((~e + ~f * ~x) ^ ~m * sin(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x)) / ~b - (~d / ~b) * Antiderivative(((~e + ~f * ~x) ^ ~m * sin(~a + ~b * ~x) * cos(~c + ~d * ~x)) / (~c + ~d * ~x), ~x)) - ((~f * ~m) / ~b) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * sin(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~m * sin(~(a') + ~(b') * ~x) * SinIntegral(~(c') + ~(d') * ~x), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * sin(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x)) / (~f * (~m + 1)) - (~d / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * sin(~a + ~b * ~x) * sin(~c + ~d * ~x)) / (~c + ~d * ~x), ~x)) - (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * cos(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x) * CosIntegral(~(c') + ~(d') * ~x), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * cos(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x)) / (~f * (~m + 1)) - (~d / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * cos(~a + ~b * ~x) * cos(~c + ~d * ~x)) / (~c + ~d * ~x), ~x)) + (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * sin(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative(cos(~(a') + ~(b') * ~x) * SinIntegral(~(c') + ~(d') * ~x), ~x) => (sin(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x)) / ~b - (~d / ~b) * Antiderivative((sin(~a + ~b * ~x) * sin(~c + ~d * ~x)) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative(sin(~(a') + ~(b') * ~x) * CosIntegral(~(c') + ~(d') * ~x), ~x) => (-(cos(~a + ~b * ~x)) * CosIntegral(~c + ~d * ~x)) / ~b + (~d / ~b) * Antiderivative((cos(~a + ~b * ~x) * cos(~c + ~d * ~x)) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x) * SinIntegral(~(c') + ~(d') * ~x), ~x) => (((~e + ~f * ~x) ^ ~m * sin(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x)) / ~b - (~d / ~b) * Antiderivative(((~e + ~f * ~x) ^ ~m * sin(~a + ~b * ~x) * sin(~c + ~d * ~x)) / (~c + ~d * ~x), ~x)) - ((~f * ~m) / ~b) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * sin(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * sin(~(a') + ~(b') * ~x) * CosIntegral(~(c') + ~(d') * ~x), ~x) => (-((~e + ~f * ~x) ^ ~m) * cos(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x)) / ~b + (~d / ~b) * Antiderivative(((~e + ~f * ~x) ^ ~m * cos(~a + ~b * ~x) * cos(~c + ~d * ~x)) / (~c + ~d * ~x), ~x) + ((~f * ~m) / ~b) * Antiderivative((~e + ~f * ~x) ^ (~m - 1) * cos(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~(m') * cos(~(a') + ~(b') * ~x) * SinIntegral(~(c') + ~(d') * ~x), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * cos(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x)) / (~f * (~m + 1)) - (~d / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * cos(~a + ~b * ~x) * sin(~c + ~d * ~x)) / (~c + ~d * ~x), ~x)) + (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * sin(~a + ~b * ~x) * SinIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative((~(e') + ~(f') * ~x) ^ ~m * sin(~(a') + ~(b') * ~x) * CosIntegral(~(c') + ~(d') * ~x), ~x) => (((~e + ~f * ~x) ^ (~m + 1) * sin(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x)) / (~f * (~m + 1)) - (~d / (~f * (~m + 1))) * Antiderivative(((~e + ~f * ~x) ^ (~m + 1) * sin(~a + ~b * ~x) * cos(~c + ~d * ~x)) / (~c + ~d * ~x), ~x)) - (~b / (~f * (~m + 1))) * Antiderivative((~e + ~f * ~x) ^ (~m + 1) * cos(~a + ~b * ~x) * CosIntegral(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative(SinIntegral(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * SinIntegral(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ~b * ~d * ~n * Antiderivative(sin(~d * (~a + ~b * log(~c * (~x) ^ ~n))) / (~d * (~a + ~b * log(~c * (~x) ^ ~n))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(CosIntegral(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * CosIntegral(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ~b * ~d * ~n * Antiderivative(cos(~d * (~a + ~b * log(~c * (~x) ^ ~n))) / (~d * (~a + ~b * log(~c * (~x) ^ ~n))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~F)(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))) / ~x, ~x) => (1 / ~n) * Subst((~F)(~d * (~a + ~b * ~x)), ~x, log(~c * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && MemberQ([SinIntegral, CosIntegral], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * SinIntegral(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * SinIntegral(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((~b * ~d * ~n) / (~m + 1)) * Antiderivative(((~e * ~x) ^ ~m * sin(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~d * (~a + ~b * log(~c * (~x) ^ ~n))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * CosIntegral(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * CosIntegral(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((~b * ~d * ~n) / (~m + 1)) * Antiderivative(((~e * ~x) ^ ~m * cos(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~d * (~a + ~b * log(~c * (~x) ^ ~n))), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
end

max_stars_repo_path: src/rules/8 Special functions/8.4 Trig integral functions.jl


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

@testset "Initializers" begin
  @testset "Bilinear initializer" begin
    # Setup a filter with scale = 2
    expectedFilter = Float32[
                       0.0625 0.1875 0.1875 0.0625;
                       0.1875 0.5625 0.5625 0.1875;
                       0.1875 0.5625 0.5625 0.1875;
                       0.0625 0.1875 0.1875 0.0625]
    filter = mx.zeros(Float32, 4, 4, 1, 4)
    mx.init(mx.XavierInitializer(), :upsampling0_weight, filter)

    mx.@nd_as_jl ro=filter begin
      for s in 1:size(filter, 4)
        @test all(filter[:, :, 1, s] .== expectedFilter)
      end
    end
  end
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
# Defining our own `_unsafe_getindex!` means we don't have to rely on Base internals.
# If there is a performance cost, we may want to switch to a copy of the Base version,
# but we would then need to call `to_indices` earlier to avoid passing in `Colon`.
function _unsafe_getindex!(
    dest::AbstractArray,
    src::AbstractArray,
    I::Vararg{Union{Real, AbstractArray, Colon}, N}
) where N  # `where N` used to force specialisation
    @inbounds for (i, j) in zip(eachindex(dest), Iterators.product(to_indices(src, I)...))
        dest[i] = src[j...]
    end
    return dest
end

#=
Optimized version (from Base) of the `_unsafe_getindex!` function above.
Can use this if it improves performance, but ideally we'd avoid an `@generated` function.

# This code is copied from Base/multidimensional.jl
# https://github.com/JuliaLang/julia/blob/f12cde0e8c7f4268838d3201dd9e13b90ea44633/base/multidimensional.jl#L761-L777
# Does not support `Colon`, as it assumes `I == to_indices(src, I)`.
@generated function _unsafe_getindex!(dest::AbstractArray, src::AbstractArray, I::Vararg{Union{Real, AbstractArray}, N}) where N
    quote
        Base.@_inline_meta
        D = eachindex(dest)
        Dy = iterate(D)
        @inbounds Base.Cartesian.@nloops $N j d->I[d] begin
            # This condition is never hit, but at the moment
            # the optimizer is not clever enough to split the union without it
            Dy === nothing && return dest
            (idx, state) = Dy
            dest[idx] = Base.Cartesian.@ncall $N getindex src j
            Dy = iterate(D, state)
        end
        return dest
    end
end
=#


text:
# _sequential_tasks_param
concept_sequential_tasks(x) = x[1] ≤ x[2] || x[3] ≤ x[4]

"""
    sum_equal_param
Constraint ensuring that the start and completion times of two tasks are not intersecting: `x[1] ≤ x[2] || x[3] ≤ x[4]`.
"""
const sequential_tasks = Constraint(
    args_length = 4,
    concept = concept_sequential_tasks,
    error = make_error(:sequential_tasks),
)


text:
module NeuronBuilder

using ModelingToolkit

const t = Num(Sym{Real}(:t))
const D = Differential(t)
export t, D

include("channel_types.jl")
export Synapse, IonChannel, Conductance, Soma, Compartment

function channel_dynamics() end

include("assets/Liu/Liu.jl")
using .Liu
export Liu

include("assets/Prinz/Prinz.jl")
using .Prinz
export Prinz

export channel_dynamics

include("synapses.jl")
include("build_neuron.jl")
include("compartments.jl")

export build_channel, build_synapse, build_neuron, build_group
export add_connection
#export Liu_NaV, Liu_CaS, Liu_CaT, Liu_H, Liu_Ka, Liu_KCa, Liu_Leak, Liu_Kdr
#export Prinz_NaV, Prinz_CaS, Prinz_CaT, Prinz_H, Prinz_Ka, Prinz_KCa, Prinz_Leak, Prinz_Kdr

#
# include("assets/Prinz/Prinz.jl")

end