text:
#---Begin-Functions-of-the-physical-model---
"""
    pressure(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")

Calculate the pressure at position `x` at time `t`.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `Fpin_itp`: Interpolated (linear) flow `F(t)` resp. inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m. Can be a function of position `x`.
* `gas`: Name of the mobile phase gas.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`,
    eq. 2)
    or with a gradient (`ng = false`, eq. 1).
* `vis`: used model for viscosity "Blumberg" or "HP"
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``p(x,t) =
\\sqrt(p_{in}(t)^2-\\frac{κ(x,t)}{κ_L(t)}\\left(p_{in}^2-p_{out}^2\\right))``
Eq. 1

``p(x,t) =
\\sqrt(p_{in}(t)^2-\\frac{x}{L}\\left(p_{in}^2-p_{out}^2\\right))`` Eq. 2

with ``κ(x,t)`` the flow restriction up to position `x` at time `t` and
``κ_L(t) = κ(x=L,t)`` the flow restriction of the whole column at
time `t`.

See also: [`flow_restriction`](@ref)
"""
function pressure(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")
    if ng==true
        if control == "Pressure"
            pin_itp = Fpin_itp
            pp = sqrt(pin_itp(t)^2 - x/L*(pin_itp(t)^2-pout_itp(t)^2))
        elseif control == "Flow"
            F_itp = Fpin_itp
            pp = sqrt(pout_itp(t)^2 + 256/π * pn/Tn * viscosity(x, t, T_itp, gas, vis=vis)* T_itp(x,t)/d(x)^4 * F_itp(t) * (L - x))
        end
    else
        if control == "Pressure"
            pin_itp = Fpin_itp
            pp = sqrt(pin_itp(t)^2 - flow_restriction(x, t, T_itp, d, gas, vis=vis)/flow_restriction(L, t, T_itp, d, gas, vis=vis)*(pin_itp(t)^2-pout_itp(t)^2))
        elseif control == "Flow"
            F_itp = Fpin_itp
            κL = flow_restriction(L, t, T_itp, d, gas, vis=vis)
            κx = flow_restriction(x, t, T_itp, d, gas, vis=vis)
            pp = sqrt(pout_itp(t)^2 + 256/π * pn/Tn * F_itp(t) * (κL - κx))
        end
    end
    return pp
end

"""
    flow_restriction(x, t, T_itp, d, gas; ng=false, vis="Blumberg")

Calculate the flow restriction ``κ`` up to position `x` at time `t`.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `d`: Diameter of the GC column, in m. Can be a function of position `x`.
* `gas`: Name of the mobile phase gas.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`,
    eq. 2)
    or with a gradient (`ng = false`, eq. 1).
* `vis`: used model for viscosity "Blumberg" or "HP"

``κ(x,t) = \\int_0^x \\frac{η(y,t) T(y,t)}{d(y)^4}dy``
Eq. 1

``κ(x,t) = \\frac{η(t) T(t) x}{d^4}`` Eq. 2

with ``η(x,t)`` the viscosity of the mobile phase gas.

See also: [`viscosity`](@ref)
"""
function flow_restriction(x, t, T_itp, d, gas; ng=false, vis="Blumberg")
    if ng==true
        κ = x*viscosity(x, t, T_itp, gas, vis=vis)*T_itp(x, t)*d(x)^-4
    else
        κ = quadgk(y -> viscosity(y, t, T_itp, gas, vis=vis)*T_itp(y, t)*d(y)^-4, 0, x)[1]
    end
    return κ
end

"""
    viscosity(x, t, T_itp, gas; vis="Blumberg")

Calculate the (dynamic) viscosity of the mobile phase gas at position `x`
at time `t` in Pa s.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `gas`: Name of the mobile phase gas.
* `vis`: used model

`vis = "HP"`

Simple model used in the HP Flow calculator

``η(x,t) = C_1 * \\left(T(x,t) + T_{st}\\right) + C_2`

`vis = "Blumberg"`

``η(x,t) = η_{st}\\left(\\frac{T(x,t)}{T_{st}}\\right)^{(ξ_0 + ξ_1 \\frac{T(x,t)-T_{st}}{T_{st}})}``

with ``η_{st}``, ``ξ_0`` and ``ξ_1`` parameters dependent on the
mobile phase gas [1].

[1] Blumberg, Leonid M., Temperature-Programmed Gas Chromatography,
Wiley-VCH, 2010.
"""
function viscosity(x, t, T_itp, gas; vis="Blumberg")
    if vis == "Blumberg"
        if gas=="He"
            ηst = 18.63e-6
            ξ₀ = 0.6958
            ξ₁ = -0.0071
        elseif gas=="H2"
            ηst = 8.382e-6
            ξ₀ = 0.6892
            ξ₁ = 0.005
        elseif gas=="N2"
            ηst = 16.62e-6
            ξ₀ = 0.7665
            ξ₁ = -0.0378
        elseif gas=="Ar"
            ηst = 21.04e-6
            ξ₀ = 0.8131
            ξ₁ = -0.0426
        else
            error("Unknown selection of gas. Choose one of these: He, H2, N2 or Ar.")
        end
        T = T_itp(x, t)
        η = ηst*(T/Tst)^(ξ₀ + ξ₁*(T-Tst)/Tst)
    elseif vis == "HP"
        if gas=="He"
            C₁ = 4.28e-8
            C₂ = 6.968e-6
        elseif gas=="H2"
            C₁ = 3.5e-8
            C₂ = 7.994e-6
        elseif gas=="N2"
            C₁ = 1.83e-8
            C₂ = 4.416e-6
        else
            error("Unknown selection of gas. Choose one of these: He, H2 or N2.")
        end
        T = T_itp(x, t)
        η = C₁*T + C₂
    else
        error("Unknown selection for the viscosity model. Choose one of these options: 'Blumberg' or 'HP'.")
    end
    return η
end

"""
    viscosity(T, gas; vis="Blumberg")

Calculate the (dynamic) viscosity of the mobile phase gas at temperature `T` in Pa s.

# Arguments
* `T`: Temperature in K.
* `gas`: Name of the mobile phase gas.
* `vis`: used model for viscosity "Blumberg" or "HP".

`vis = "HP"`

Simple model used in the HP Flow calculator

``η(x,t) = C_1 * \\left(T(x,t) + T_{st}\\right) + C_2`

`vis = "Blumberg"``

``η(x,t) = η_{st}\\left(\\frac{T)}{T_{st}}\right)^{(ξ_0 + ξ_1 \\frac{T-T_{st}}{T_{st}})}``

with ``η_{st}``, ``ξ_0`` and ``ξ_1`` parameters dependent on the
mobile phase gas [1].

[1] Blumberg, Leonid M., Temperature-Programmed Gas Chromatography,
Wiley-VCH, 2010.
"""
function viscosity(T::Float64, gas::String; vis="Blumberg")
    if vis == "Blumberg"
        if gas=="He"
            ηst = 18.63e-6
            ξ₀ = 0.6958
            ξ₁ = -0.0071
        elseif gas=="H2"
            ηst = 8.382e-6
            ξ₀ = 0.6892
            ξ₁ = 0.005
        elseif gas=="N2"
            ηst = 16.62e-6
            ξ₀ = 0.7665
            ξ₁ = -0.0378
        elseif gas=="Ar"
            ηst = 21.04e-6
            ξ₀ = 0.8131
            ξ₁ = -0.0426
        else
            error("Unknown selection of gas. Choose one of these: He, H2, N2 or Ar.")
        end
        η = ηst*(T/Tst)^(ξ₀ + ξ₁*(T-Tst)/Tst)
    elseif vis == "HP"
        if gas=="He"
            C₁ = 4.28e-8
            C₂ = 6.968e-6
        elseif gas=="H2"
            C₁ = 3.5e-8
            C₂ = 7.994e-6
        elseif gas=="N2"
            C₁ = 1.83e-8
            C₂ = 4.416e-6
        else
            error("Unknown selection of gas. Choose one of these: He, H2 or N2.")
        end
        η = C₁*T + C₂
    else
        error("Unknown selection for the viscosity model. Choose one of these options: 'Blumberg' or 'HP'.")
    end
    return η
end

"""
    inlet_pressure(t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")

Calculate the inlet pressure for a given column, temperature, flow and outlet pressure
"""
function inlet_pressure(t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")
    if control == "Pressure"
        pin = Fpin_itp(t)
    elseif control == "Flow"
        F_itp = Fpin_itp
        if ng == true
            pin = sqrt(pout_itp(t)^2 + 256/π * pn/Tn * viscosity(0.0, t, T_itp, gas; vis="Blumberg") * T_itp(0.0, t) * L / d(0.0)^4 * F_itp(t))
        else
            κL = flow_restriction(L, t, T_itp, d, gas, vis=vis)
            pin = sqrt(pout_itp(t)^2 + 256/π * pn/Tn * κL * F_itp(t))
        end
    end
    return pin
end

"""
    holdup_time(T, Fpin, pout, L, d, gas; vis="Blumberg", control="Pressure")

Calculate the hold-up time in s without a gradient.

# Arguments
* `T`: Temperature in K.
* `Fpin`: Flow in m³/s resp. inlet pressure in Pa(a).
* `pout`: Outlet pressure in Pa(g).
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``t_M = \\frac{128}{3}\\frac{L^2}{d^2}η\\frac{p_{in}^3-p_{out}^3}{(p_{in}^2-p_{out}^2)^2}``
"""
function holdup_time(T::Float64, Fpin::Float64, pout::Float64, L::Float64, d::Float64, gas::String; vis="Blumberg", control="Pressure")
    # hold-up time at temperature T (non-gradient)
	η = viscosity(T, gas; vis=vis)
    if control == "Pressure"
        pin = Fpin
    elseif control == "Flow"
        F = Fpin
        pin = sqrt(pout^2 + 256/π * pn/Tn * viscosity(T, gas; vis="Blumberg") * T * L / d^4 * F)
    end
	tM = 128/3*L^2/d^2*η*(pin^3-pout^3)/(pin^2-pout^2)^2
	return tM
end

"""
    holdup_time(t, T_itp, pin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")

Calculate the hold-up time in s at time `t` with a gradient.

# Arguments
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `pin_itp`: Interpolated (linear) inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`,
    eq. 2)
    or with a gradient (`ng = false`, eq. 1).
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``t_M(t) = 64\\frac{κ_L(t)}{p_{in}(t)^2-p_{out}(t)^2} \\int_0^L
d(y)^2\\frac{p(y,t)}{T(y,t)}dy`` Eq. 1

``t_M(t) =
\\frac{128}{3}\\frac{L^2}{d^2}η(t)\\frac{p_{in}(t)^3-p_{out}(t)^3}{(p_{in}(t)^2-p_{out}(t)^2)^2}``
Eq. 2
"""
function holdup_time(t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")
    # hold-up time at time t in a temperature program with potential thermal gradient
    if control == "Pressure"
        #pin(t) = Fpin_itp(t)
        if ng==true
            η = GasChromatographySimulator.viscosity(L, t, T_itp, gas; vis=vis)
            tM = 128/3*L^2/d(L)^2*η*(Fpin_itp(t)^3-pout_itp(t)^3)/(Fpin_itp(t)^2-pout_itp(t)^2)^2
        else
            κL = flow_restriction(L, t, T_itp, d, gas; ng=false, vis=vis)
            integral = quadgk(y -> d(y)^2*pressure(y, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis=vis, control="Pressure")/T_itp(y, t), 0, L)[1]
            tM = 64*κL/(Fpin_itp(t)^2-pout_itp(t)^2)*integral
        end
    elseif control == "Flow"
        if ng==true
            pin = GasChromatographySimulator.inlet_pressure(t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=true, vis=vis, control="Flow")
            η = GasChromatographySimulator.viscosity(L, t, T_itp, gas; vis=vis)
            tM = 128/3*L^2/d(L)^2*η*(pin^3-pout_itp(t)^3)/(pin^2-pout_itp(t)^2)^2
        else
            integral = quadgk(y -> d(y)^2*pressure(y, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis=vis, control="Flow")/T_itp(y, t), 0, L)[1]
            tM = π/4 * Tn/pn * integral/Fpin_itp(t)
        end
    end

    return tM
end

"""
    flow(T, Fpin, pout, L, d, gas; vis="Blumberg", control="Pressure")

Calculate the normalized flow through the GC column in m³/s without a gradient.

# Arguments
* `T`: Temperature in K.
* `Fpin`: Flow in m³/s resp. inlet pressure in Pa(a).
* `pout`: Outlet pressure in Pa(g).
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``F =
\\frac{π}{256}\\frac{T_n}{p_n}\\frac{d^4}{L}\\frac{p_{in}^2-p_{out}^2}{η T}``

with ``T_n`` the normalized temperature (``T_n=(25 + 273.15)``K), ``p_n``
the normalized pressure (``p_n = 101300`` Pa(a)) and ``η`` the viscosity
the mobile phase gas at temperature ``T``.
"""
function flow(T::Float64, Fpin::Float64, pout::Float64, L::Float64, d::Float64, gas::String; vis="Blumberg", control="Pressure")
	# normalized Flow at temperature T (non-gradient)
    if control == "Pressure"
        pin = Fpin
	    η = viscosity(T, gas; vis=vis)
	    F = π/256 * Tn/pn * d^4/L * (pin^2-pout^2)/(η*T)
    elseif control == "Flow"
        F = Fpin
    end
	return F
end

"""
    flow(t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")

Calculate the normalized flow through the GC column in m³/s at time `t`.

# Arguments
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `Fpin_itp`: Interpolated (linear) Flow F(t) resp. inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`,
    eq. 2)
    or with a gradient (`ng = false`, eq. 1).
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``F(t) =
\\frac{π}{256}\\frac{T_n}{p_n}\\frac{p_{in}(t)^2-p_{out}(t)^2}{κ_L(t)}``
Eq. 1

``F(t) =
\\frac{π}{256}\\frac{T_n}{p_n}\\frac{d^4}{L}\\frac{p_{in}(t)^2-p_{out}(t)^2}{η(t)
T(t)}``
Eq. 2

with ``T_n`` the normalized temperature (``T_n=(25 + 273.15)``K), ``p_n``
the normalized pressure (``p_n = 101300`` Pa(a)), ``κ_L`` the flow
restriction of the column and ``η`` the viscosity
the mobile phase gas at temperature ``T``.
"""
function flow(t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")
	# normalized Flow at time t in a temperature program with potential thermal
	# gradient
    # TODO: test for gradient in d(x)
    if control == "Pressure"
        pin_itp = Fpin_itp
        if ng==true
            η = GasChromatographySimulator.viscosity(L, t, T_itp, gas, vis=vis)
            F = π/256 * Tn/pn * d(L)^4/L * (pin_itp(t)^2-pout_itp(t)^2)/(η*T_itp(L,t))
        else
            κL = flow_restriction(L, t, T_itp, d, gas; ng=false, vis=vis)
            F = π/256 * Tn/pn * (pin_itp(t)^2-pout_itp(t)^2)/κL
        end
    elseif control == "Flow"
        F = Fpin_itp(t)
    end
	return F
end

"""
    mobile_phase_residency(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")

Calculate the residency (the inverse velocity) of the mobile phase at
position `x` at time `t`.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `Fpin_itp`: Interpolated (linear) Flow F(t) resp. inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`)
or with a gradient (`ng = false`).
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``r_M(x,t) = 64 \\frac{d^2 κ_L}{T(x,t)}\\frac{p(x,t)}{p_{in}^2-p_{out}^2}``

with ``T_n`` the normalized temperature (``T_n=(25 + 273.15)``K), ``p_n``
the normalized pressure (``p_n = 101300`` Pa(a)), ``κ_L`` the flow
restriction of the column and ``p(x,t)`` the local pressure.

See also: [`pressure`](@ref), [`flow_restriction`](@ref)
"""
function mobile_phase_residency(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=false, vis="Blumberg", control="Pressure")
    if control == "Pressure"
        pin_itp = Fpin_itp
        pp = pressure(x, t, T_itp, pin_itp, pout_itp, L, d, gas; ng=ng, vis=vis, control="Pressure")
        κL = flow_restriction(L, t, T_itp, d, gas; ng=ng, vis=vis)
        rM = 64*(pp*(d(x))^2)/T_itp(x, t)*κL/(pin_itp(t)^2-pout_itp(t)^2)
    elseif control == "Flow"
        F_itp = Fpin_itp
        pp = pressure(x, t, T_itp, F_itp, pout_itp, L, d, gas; ng=ng, vis=vis, control="Flow")
        rM = π/4 * Tn/pn * d(x)^2/F_itp(t) * pp/T_itp(x,t)
    end
    return rM
end

"""
    residency(x, t, T_itp, Fpin_itp, pout_itp, L, d, df, gas, Tchar, θchar, ΔCp,  φ₀; ng=false, vis="Blumberg", control="Pressure")

Calculate the residency (the inverse velocity) of the solute at
position `x` at time `t`.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `Fpin_itp`: Interpolated (linear) Flow F(t) resp. inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `df`: Film thickness of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `Tchar`: Characteristic temperature of the solute, in K.
* `θchar`: Characteristic parameters of the solute, in °C.
* `ΔCp`: Change of the isobaric heat capacity of the solute moving from the mobile to the
stationary phase, in J mol⁻¹ K⁻¹.
* `φ₀`: Dimensionless film thickness (φ ≈ df/d) of the column for which the
thermodynamic parameters (Tchar, θchar, ΔCp) were estimated.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`)
or with a gradient (`ng = false`).
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``r(x,t) = r_M(x,t) \\left(1+k(x,t)\\right)``

with ``r_M`` the residency of the mobile phase and ``k(x,t)`` the retention
factor of the solute on the stationary phase.

See also: [`mobile_phase_residency`](@ref), [`retention_factor`](@ref)
"""
function residency(x, t, T_itp, Fpin_itp, pout_itp, L, d, df, gas, Tchar, θchar, ΔCp, φ₀; ng=false, vis="Blumberg", control="Pressure")
    r = mobile_phase_residency(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=ng, vis=vis, control=control)*(1 + retention_factor(x, t, T_itp, d, df, Tchar, θchar, ΔCp, φ₀))
    return r
end

"""
    retention_factor(x, t, T_itp, d, df, Tchar, θchar, ΔCp, φ₀)

Calculate the retention factor of the solute in the stationary phase at
position `x` at time `t`.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `d`: Diameter of the GC column, in m.
* `df`: Film thickness of the GC column, in m.
* `Tchar`: Characteristic temperature of the solute, in K.
* `θchar`: Characteristic parameters of the solute, in °C.
* `ΔCp`: Change of the isobaric heat capacity of the solute moving from the mobile to the
stationary phase, in J mol⁻¹ K⁻¹.
* `φ₀`: Dimensionless film thickness (φ ≈ df/d) of the column for which the
thermodynamic parameters (Tchar, θchar, ΔCp) were estimated.

``k(x,t) = \\frac{φ}{φ₀}
\\exp{\\left((\\frac{ΔC_p}{R}+\\frac{T_{char}}{θ_{char}})(\\frac{T_{char}}{T}+-1)
    \\frac{ΔC_p}{R}\\ln{(\\frac{T}{T_{char}})}\\right)}``

with ``R`` the molar gas constant and ``φ`` the dimensionless film thickness
of the simulated GC Column (``φ = d_f/d``).

**TODO**: add option for the retention model ('ABC', 'K-centric')
"""
function retention_factor(x, t, T_itp, d, df, Tchar, θchar, ΔCp, φ₀)
    # this version of the function, where every parameter is
    # given to the function separatly seems to be the fastest
    # version
    # for now only the ideal thermodynamic model
    T = T_itp(x, t)
    φ = df(x)/d(x)
    C = ΔCp/R
    lnk₀ = (C + Tchar/θchar) * (Tchar/T - 1) + C*log(T/Tchar)
    k = φ/φ₀*exp(lnk₀)
    return k
end

"""
    plate_height(x, t, T_itp, Fpin_itp, pout_itp, L, d, df, gas, Tchar, θchar, ΔCp, φ₀, Dag; ng=false, vis="Blumberg", control="Pressure")

Calculate the plate height of the solute at position `x` at time `t`
according to the Golay equation.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `Fpin_itp`: Interpolated (linear) Flow F(t) resp. inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `df`: Film thickness of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `Tchar`: Characteristic temperature of the solute, in K.
* `θchar`: Characteristic parameters of the solute, in °C.
* `ΔCp`: Change of the isobaric heat capacity of the solute moving from the mobile to the
stationary phase, in J mol⁻¹ K⁻¹.
* `φ₀`: Dimensionless film thickness (φ ≈ df/d) of the column for which the
thermodynamic parameters (Tchar, θchar, ΔCp) were estimated.
* `Dag`: diffusivity of solute `a` in gas `g`.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`)
or with a gradient (`ng = false`).
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``H(x,t) = 2 \\frac{D_M}{u_M} + \\frac{d^2}{96}\\left(6 μ^2-16 μ +11
\\right) \\frac{u_M}{D_M} + \\frac{2}{3} d_f^2 μ(1-μ) \\frac{u_M}{D_S}``

with ``D_M`` the diffusion coefficient of the solute in the mobile phase,
``D_S`` the diffusion coefficient of the solute in the stationary phase,
``u_M`` the velocity of the mobile phase and μ the mobility of the solute.

``D_S`` is correlated to ``D_M`` by:

``D_S = \\frac{D_M}{10000}``

**TODO**: alternative correlations?

``u_M`` is realated to the residency of the mobile phase ``r_M``:

``u_M = \\frac{1}{r_M}``

μ is correlated to the retention factor ``k``:

``μ = \\frac{1}{1 + k}``

See also: [`diffusion_mobile`](@ref), [`mobile_phase_residency`](@ref), [`retention_factor`](@ref)
"""
function plate_height(x, t, T_itp, Fpin_itp, pout_itp, L, d, df, gas, Tchar, θchar, ΔCp, φ₀, Dag; ng=false, vis="Blumberg", control="Pressure")
    id = d(x)# - 2.0*df(x)
    uM = 1/mobile_phase_residency(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=ng, vis=vis, control=control)
    μ = 1/(1 + retention_factor(x, t, T_itp, d, df, Tchar, θchar, ΔCp, φ₀))
    DM = diffusion_mobile(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas, Dag; ng=ng, vis=vis, control=control)
    DS = DM/10000
    H1 = 2*DM/uM
    H2 = id^2/96*(6*μ^2-16*μ+11)*uM/DM
    H3 = 2/3*df(x)^2*μ*(1-μ)*uM/DS
    H = H1 + H2 + H3
    return H
end

"""
    diffusion_mobile(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas, Dag; ng=false, vis="Blumberg", control="Pressure")

Calculate the diffusion coefficient of the solute in the mobile phase at
position `x` at time `t`.

# Arguments
* `x`: Position along the GC column, in m.
* `t`: Time in s.
* `T_itp`: Interpolated (linear) temperature `T(x,t)`.
* `Fpin_itp`: Interpolated (linear) Flow F(t) resp. inlet pressure `pin(t)`.
* `pout_itp`: Interpolated (linear) outlet pressure `pout(t)`.
* `L`: Length of the capillary measured in m (meter).
* `d`: Diameter of the GC column, in m.
* `gas`: Name of the mobile phase gas.
* `Dag`: diffusivity of solute `a` in gas `g`.
* `ng`: Option to calculate the simulation without a gradient (`ng = true`)
or with a gradient (`ng = false`).
* `vis`: used model for viscosity "Blumberg" or "HP".
* `control`: Control of the "Flow" or of the "Pressure" (at column inlet) during the program

``D_M(x,t) = D_{ag} \\frac{T(x,t)^{1.75}}{p(x,t)}``
"""
function diffusion_mobile(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas, Dag; ng=false, vis="Blumberg", control="Pressure")
    DM = T_itp(x, t)^1.75/pressure(x, t, T_itp, Fpin_itp, pout_itp, L, d, gas; ng=ng, vis=vis, control=control)*Dag
    return DM
end

"""
    diffusivity(M, Cn, Hn, On, Nn, Rn, gas)

Calculate the diffusivity `Dag` of solute `a` in gas `g` using the
emperical Fuller-Schettler-Giddings model [1].

[1] Fuller, Edward N.; Ensley, Keith; Giddings, J. Calvin, Diffusion of
Halogenated Hydrocarbons in Helium. The Effect of Structure on Collision
Cross Sections, The Journal of Physical Chemistry, Volume 73, Issue 11,
1969, 3679–3685

# Arguments
* `M`: Molar mass of the solute.
* `Cn`: Number of carbon atoms of the solute.
* `Hn`: Number of hydrogen atoms of the solute.
* `On`: Number of oxygen atoms of the solute.
* `Nn`: Number of nitrogen atoms of the solute.
* `Rn`: Number of closed rings of the structure of the solute.
* `gas`: The name of the mobile phase. Allowed values: He, H2 or N2.
"""
function diffusivity(M, Cn, Hn, On, Nn, Rn, gas)
    # calculates diffusitivity Dag of an analyte in a gas
    # from the (simplified) molecular formula of the solute
    # using the Fuller-Schettler-Giddings equation
    #
    # M ... molar mass
    # Cn ... number of carbons
    # Hn ... number of hydrogens
    # On ... number of oxygens
    # Nn ... number of nitrogens
    # Rn ... number of ring structures
    # gas ... the used mobile phase gas
    # **TODO**: add other atoms (P, Cl, Si, ...)
    if gas=="H2"
        Vg = 6.12
        Mg = 2.02
    elseif gas=="He"
        Vg = 2.67
        Mg = 4
    elseif gas=="N2"
        Vg = 18.5
        Mg = 28.01
    elseif gas=="Ar"
        Vg = 16.2
        Mg = 39.95
    else
        error("Unknown selection of gas. Choose one of these: He, H2, N2 or Ar.")
    end
    Va = 15.9*Cn + 2.31*Hn + 6.11*On + 4.54*Nn - 18.3*Rn
    Dag = pn*sqrt(1/M+1/Mg)/(Vg^(1/3)+Va^(1/3))^2*1e-7 # m²/s
    return Dag
end
#---End-Functions-of-the-physical-model---


text:
module Detect_test

using SpikeSorting

if VERSION > v"0.7-"
    using Test
else
    using Base.Test
end

cluster=ClusterOSort();
align=AlignMax();
feature=FeatureTime();
reduce=ReductionNone();
thres=ThresholdMeanN();
num_channels=1;

(buf,nums)=output_buffer(num_channels);
v=rand(1:1000, 1000, num_channels);

#=
Power

Rutishauser et al 2006
=#
#=
d=DetectPower()
s1=create_multi(d,cluster,align,feature,reduce,thres,num_channels);

cal!(s1,v,buf,nums,0)
cal!(s1,v,buf,nums,1)
cal!(s1,v,buf,nums,2)
onlinesort!(s1,v,buf,nums);

    @fact s1[1].d.a --> less_than(s1[1].d.b)
    @fact s1[1].d.a --> greater_than(s1[1].d.c)

=#

#=
Raw Signal

=#

d=DetectNeg()
s1=create_multi(d,cluster,align,feature,reduce,thres,num_channels);

cal!(s1,v,buf,nums,0)
cal!(s1,v,buf,nums,1)
cal!(s1,v,buf,nums,2)
onlinesort!(s1,v,buf,nums);

end

ext: jl


text:
using RegistryTools: gitcmd, Compress
using Pkg: Pkg, TOML

# Read `project_file` and create or update a corresponding bare bones
# package directory under `package_dir`. Commit the changes to git.
#
# Strictly speaking it is only the "Project.toml" file and the src
# file that matter. Additional files are added to make the set up
# somewhat more normal looking.
function prepare_package(packages_dir, project_file, subdir = "";
                         use_julia_project = false)
    project_file = joinpath(@__DIR__, "project_files", project_file)
    project_data = TOML.parsefile(project_file)
    name = project_data["name"]
    version = project_data["version"]
    # Fake repository URL.
    repo = "git@example.com:Julia/$(name).jl.git"
    top_dir = joinpath(packages_dir, name)
    package_dir = joinpath(top_dir, subdir)
    mkpath(package_dir)
    mkpath(joinpath(package_dir, "src"))
    mkpath(joinpath(package_dir, "test"))
    if !isempty(subdir)
        write(joinpath(top_dir, "README.md"), "# Top Level README")
    end
    git = gitcmd(top_dir, TEST_GITCONFIG)
    if !isdir(joinpath(top_dir, ".git"))
        run(`$(git) init -q`)
        run(`$git remote add origin $repo`)
    end
    if use_julia_project
        project = read(project_file, String)
        write(joinpath(package_dir, "JuliaProject.toml"), project)
        # Write a corrupted Project.toml as well to test that the
        # right file is chosen. (This makes an assumption of the
        # package name used for this case.)
        write(joinpath(package_dir, "Project.toml"),
              replace(project, "JuliaProjectTest" => "ProjectTest"))
    else
        write(joinpath(package_dir, "Project.toml"), read(project_file, String))
    end
    write(joinpath(package_dir, "README.md"), "# $(name)\n")
    write(joinpath(package_dir, "LICENSE"), "$(name) is in the public domain\n")
    write(joinpath(package_dir, "src", "$(name).jl"), "module $(name)\nend\n")
    write(joinpath(package_dir, "test", "runtests.jl"),
          "using Test\n@test true\n")
    run(`$git add --all`)
    vers = "Version $(version)"
    @show vers
    if !occursin("nothing to commit, working tree clean",
                read(`$git status`, String))
        @info("nothing to commit")
        run(`$git commit -qm $(vers)`)
    end
    return
end

function readdir_excluding_git(dir)
    return filter(!isequal(".git"), readdir(dir))
end

function check_result(actual_result_dir, expected_result_dir)
    check1 = compare_file_trees(actual_result_dir,
                                joinpath(@__DIR__, "expected_results",
                                         expected_result_dir))
    check2 = sanity_check_registry(actual_result_dir)
    return check1 && check2
end

function compare_file_trees(path1, path2)
    if isdir(path1) && isdir(path2)
        dir1 = readdir_excluding_git(path1)
        if sort(dir1) != sort(readdir_excluding_git(path2))
            println("Directories $(path1) and $(path2) differ.")
            return false
        end
        for path in dir1
            if !compare_file_trees(joinpath(path1, path), joinpath(path2, path))
                return false
            end
        end
    elseif isfile(path1) && isfile(path2)
        if !compare_files(path1, path2)
            println("Files $(path1) and $(path2) differ.")
            return false
        end
    else
        println("Mismatch between paths $(path1) and $(path2)")
        return false
    end

    return true
end

# Compress.save writes files differently on e.g. Julia 1.1 and
# 1.3. Compare the Deps.toml and Compat.toml files after reading them
# with Compress.load.
function compare_files(path1, path2)
    if endswith(path1, "Deps.toml") || endswith(path1, "Compat.toml")
        return Compress.load(path1) == Compress.load(path2)
    end
    return read_normalize_line_end(path1) == read_normalize_line_end(path2)
end

function read_normalize_line_end(path)
    return replace(read(path, String), "\r\n" => "\n")
end

# Check that the deps and compat files can be read at all by a
# function used in the internals of Pkg. Return true whenever there is
# no error.
function sanity_check_registry(path)
    registry = TOML.parsefile(joinpath(path, "Registry.toml"))
    if VERSION >= v"1.7-"
        # TODO: Check if the `parse_packages` keyword is available
        # when the final Julia 1.7 has been released and if so set it
        # explicitly to true.
        registry = Pkg.Registry.RegistryInstance(path)  # , parse_packages = true)
        return true
    end

    for (uuid, package) in registry["packages"]
        package_path = joinpath(path, package["path"])
        deps_file = joinpath(package_path, "Deps.toml")
        compat_file = joinpath(package_path, "Compat.toml")
        if isdefined(Pkg.Operations, :load_package_data_raw)
            deps_data = Pkg.Operations.load_package_data_raw(Pkg.Types.UUID,
                                                             deps_file)
            compat_data = Pkg.Operations.load_package_data_raw(Pkg.Types.VersionSpec,
                                                               compat_file)
        else
            ctx = Pkg.Types.Context()
            version_info = Pkg.Operations.load_versions(ctx, package_path;
                                                        include_yanked = false)
            versions = sort!(collect(keys(version_info)))
            if applicable(Pkg.Operations.load_package_data,
                          Pkg.Types.UUID, deps_file, versions)
                deps_data = Pkg.Operations.load_package_data(Pkg.Types.UUID,
                                                             deps_file, versions)
                compat_data = Pkg.Operations.load_package_data(Pkg.Types.VersionSpec,
                                                               compat_file, versions)
            else
                deps_data = Pkg.Operations.load_package_data(ctx, Pkg.Types.UUID,
                                                             deps_file, versions)
                compat_data = Pkg.Operations.load_package_data(ctx, Pkg.Types.VersionSpec,
                                                               compat_file, versions)
            end
        end
    end

    return true
end


text:
using DiscreteEvents, Printf, Random

function watchdog(name)
    delay!(until, 6 + rand())
    now!(SF(println, @sprintf("%5.2f %s: yawn!, bark!, yawn!", tau(), name)))
    wait!(((@val :hunger :≥ 7),(@tau :≥ 6.5)))
    while 5 ≤ hunger ≤ 10
        now!(SF(println, @sprintf("%5.2f %s: %s", tau(), name, repeat("wow ", Int(trunc(hunger))))))
        delay!(rand()/2)
        if scuff
            now!(SF(println, @sprintf("%5.2f %s: smack smack smack", tau(), name)))
            global hunger = 2
            global scuff = false
        end
    end
    delay!(rand())
    now!(SF(println, @sprintf("%5.2f %s: snore ... snore ... snore", tau(), name)))
end

hunger = 0
scuff = false
reset!(𝐶)
Random.seed!(1122)

sample!(SF(()-> global hunger += rand()), 0.5)
event!(SF(()-> global scuff = true ), 7+rand())
process!(SP(1, watchdog, "Snoopy"), 1)

run!(𝐶, 10)


text:
using JSON

include("./bundlepaths.jl")
include("./jupyter.jl")

download_js_bundles()

# See https://github.com/JuliaGizmos/WebIO.jl/issues/314 for the rational behind
# why we're not installing Jupyter packages by default anymore.
@warn(
    "WebIO no longer installs Jupyter extensions automatically; please run "
    * "`WebIO.install_jupyter_notebook()` or `WebIO.install_jupyter_lab()` if "
    * "needed."
)


text:
"""
    tensorsvd(A::AbstractTensor{T,2}; svdtrunc)
returns the `svd` of an AbstractTensor. `svdtrunc` is a function that has the
singular values as input and returns a number specifying how many of them to keep.
The default `svdtrunc` is `svdtrunc_default` and keeps all of them.

Other options are:
svdtrunc_discardzero
svdtrunc_maxχ
svdtrunc_maxcumerror
svdtrunc_maxerror
"""
function tensorsvd end

function _tensorsvd(A::Array{<:Any,2}; svdtrunc = svdtrunc_default,
        helper::Bool = false)
    F = LA.svd(A)
    svals = F.S
    cutoff = svdtrunc(svals)
    U = F.U[:, 1:cutoff]
    S = LA.diagm(0=>svals[1:cutoff])
    Vt = F.Vt[1:cutoff, :]
    helper && return (U, S, Vt, cutoff)
    return (U, S, Vt)
end

"""
    tensorsvd(A::AbstractTensor, indexes; svdtrunc)
works like `tensorsvd` except that `A` can have arbitrary rank.
`indexes` specifies which indices the fuse for `A` to be a rank-2 tensor
as in `fuselegs`. The tensor is then fused, `tensorsvd` applied and split again.
"""
function tensorsvd(A::AbstractTensor, indexes; svdtrunc = svdtrunc_default)
    fA, rs = fuselegs(A, indexes)
    U, S, Vt = tensorsvd(fA, svdtrunc=svdtrunc)
    li1 = length(indexes[1])
    li2 = length(indexes[2])
    if indexes[1] isa Tuple
        indxs = (ntuple(x -> (1,1,x), li1)..., 2)
        U = splitlegs(U, indxs, rs...)
    end
    if indexes[2] isa Tuple
        indxs = (1, ntuple(x -> (2,2,x), li2)...)
        Vt = splitlegs(Vt, indxs, rs...)
    end
    return (U,S,Vt)
end

svdtrunc_default(x) = length(x)
"""
    svdtrunc_discardzero(s)
return the number of non-zero values in s
"""
svdtrunc_discardzero(x) = length(filter(!iszero, x))

"""
    svdtrunc_maxχ(χ)
return a function that given singular values `s` returns the min of the
length of `s` and χ.
"""
svdtrunc_maxχ(χ) = x -> min(length(x), χ)

"""
    svdtrunc_maxcumerror(ϵ::Real; χ::Int = typemax(Int))
return a function that returns then min of χ and l where l is the number
of singular values that need to be kept to have the truncated sum up to not
more than ϵ.
"""
svdtrunc_maxcumerror(ϵ::Real; χ::Int = typemax(Int)) = x -> _maxcumerror(x, ϵ, χ)

"""
    svdtrunc_maxerror(ϵ::Real; χ::Int = typemax(Int))
return a function that returns then min of χ and l where l is the number
of singular values that need to be kept such that the largest discarded singular
value is below ϵ.
"""
svdtrunc_maxerror(ϵ::Real; χ::Int = typemax(Int)) = x -> _maxerror(x, ϵ, χ)

function _maxcumerror(xs, ϵ, χ)
    acc, lxs = 0, length(xs)
    for i in lxs:-1:1
        acc += xs[i]/xs[1]
        acc >= ϵ && return min(i+1,χ)
    end
    return min(lxs, χ)
end

function _maxerror(xs, ϵ, χ)
    index = findfirst(x -> x/xs[1] < ϵ, xs)
    index == nothing && return min(length(xs), χ)
    return min(index-1, χ)
end

#DTensor
function tensorsvd(A::DTensor{T,2}; svdtrunc = svdtrunc_default) where T
    DTensor.(_tensorsvd(A.array, svdtrunc = svdtrunc))
end

#DASTensor
function tensorsvd(A::DASTensor{T,N,SYM,CHARGES,SIZES,CHARGE};
                   svdtrunc = svdtrunc_default) where
                   {T,N,SYM,CHARGES,SIZES,CHARGE}
    N == 2 || throw(ArgumentError("SVD only works on rank 2 tensors"))
    lch = connectingcharge(A)
    ld =  sizes(A,2)[[chargeindex(c,charges(A,2)) for c in lch]]

    U  = DASTensor{T,N}(SYM, (charges(A,1), lch),
            deepcopy.((sizes(A,1), ld)),
            in_out(A), charge(A))
    S  = DASTensor{real(T),N}(SYM, (lch, lch),
            deepcopy.((ld, ld)),
            vcat(inv(in_out(A,2)), in_out(A,2)))
    Vd = DASTensor{T,N}(SYM, (lch, charges(A,2)),
            deepcopy((ld,sizes(A,2))),
            vcat(inv(in_out(A,2)), in_out(A,2)))

    for (k, degen) in tensor(A)
        in, out = k[1], k[2]
        U[DASSector(in, out)], S[DASSector(out, out)], Vd[DASSector(out, out)], cutoff =
                _tensorsvd(degen, svdtrunc = svdtrunc, helper = true)
        U.dims[2][chargeindex(out,charges(U,2))] = cutoff
        S.dims[1][chargeindex(out,charges(S,1))] = cutoff
        S.dims[2][chargeindex(out,charges(S,2))] = cutoff
        Vd.dims[1][chargeindex(out,charges(Vd,1))] = cutoff
    end
    return (U, S, Vd)
end

"""
    connectingcharge(A)
where A is a rank-2 DASTensor,
returns the charges on the second index that can be realized
given the charges on the first index and the charge of A.
"""
function connectingcharge(A::DASTensor{<:Any,2})
    chs1, chs2 = charges(A)
    io = in_out(A)
    io[1] == InOut(-1) && (chs1 = inv(chs1))
    io[2] == InOut( 1) && (chs2 = inv(chs2))
    #                         (1)ch
    #                          ↓
    #canonical form: chs1 →(1)→A→(-1)→ lch = chs1 + ch ∩ chs2
    chs1 += charge(A)
    lch = chs1 ∩ chs2
    io[2] == InOut(1) && return inv(lch)
    return lch
end

"""
    tensorqr(A::AbstractTensor)

returns tensor Q,R such that A = QR and Q is an orthogonal/unitary matrix
and R is an upper triangular matrix.
"""
function tensorqr end

"""
    tensorqr(A::AbstractTensor, inds)

returns the `tensorqr` of `A` fused according to `inds`.
"""
function tensorqr(A::AbstractTensor, inds)
    fA, rs = fuselegs(A, inds)
    Q,R = tensorqr(fA)
    li1, li2 = length.(inds)
    if inds[1] isa Tuple
        indxs = (ntuple(x -> (1,1,x), li1)..., 2)
        Q = splitlegs(Q, indxs, rs...)
    end
    if inds[2] isa Tuple
        indxs = (1, ntuple(x -> (2,2,x), li2)...)
        R = splitlegs(R, indxs, rs...)
    end
    return (Q,R)
end

function tensorqr(a::DTensor{<:Any,2})
    f = LA.qr(a.array)
    return DTensor(Matrix(f.Q)), DTensor(Matrix(f.R))
end

function tensorqr(a::DASTensor{T,2,SYM}) where {T,SYM}
    lch = connectingcharge(a)
    ld =  sizes(a,2)[[chargeindex(c,charges(a,2)) for c in lch]]

    Q = DASTensor{T,2}(SYM, (charges(a,1), lch),
            deepcopy.((sizes(a,1), ld)),
            in_out(a), charge(a))
    R = DASTensor{T,2}(SYM, (lch, charges(a,2)),
            deepcopy((ld,sizes(a,2))),
            vcat(inv(in_out(a,2)), in_out(a,2)))

    for (k, degen) in tensor(a)
        in, out = k[1], k[2]
        Q[DASSector(in, out)], R[DASSector(out, out)] = LA.qr(degen)
    end
    return (Q,R)
end

"""
    tensorrq(A::AbstractTensor)

returns tensor R,Q such that A = RQ and Q is
obeys Q*Q' = 1 and R is triangular.
This is simply a wrapper for `tensorqr`, useful for e.g. MPS canonicalization.
"""
function tensorrq end

"""
    tensorrq(A::AbstractTensor, inds)

returns the `tensorrq` of `A` fused according to `inds`.
"""
function tensorrq(c, inds = ((1,),(2,)))
    q, r = tensorqr(c, reverse(inds))
    nq, nr = ndims.((q,r))
    q = TO.tensorcopy(q, circshift(1:nq,-1),  1:nq)
    r = TO.tensorcopy(r, circshift(1:nr,1), 1:nr)
    return r, q
end

function tensoreig end

function tensoreig(A::AbstractTensor, indexes; truncfun = svdtrunc_default)
    fA, rs = fuselegs(A, indexes)
    E, D, Ep = tensoreig(fA, truncfun=svdtrunc)
    li1, li2 = length.(indexes)

    if indexes[1] isa Tuple
        indxs = (ntuple(x -> (1,1,x), li1)..., 2)
        E = splitlegs(E, indxs, rs...)
    end
    if indexes[2] isa Tuple
        indxs = (1, ntuple(x -> (2,2,x), li2)...)
        Ep = splitlegs(Ep, indxs, rs...)
    end
    return (E, D, Ep)
end

function _tensoreig(a::Matrix; truncfun = length, helper=false)
    if LA.ishermitian(a)
        evals, evecs = LA.eigen(LA.Hermitian(a))
    else
        evals, evecs = LA.eigen(a)
    end

    p = sortperm(evals, by=abs, rev=true)
    evals = evals[p]
    evecs = evecs[:,p]
    cutoff = truncfun(evals)
    E = evecs[:,1:cutoff]
    d = LA.diagm(0 => evals[1:cutoff])
    helper && return E, d, E', cutoff
    return E, d, E'
end

function tensoreig(a::DTensor{T,2}; truncfun = length) where T
    E, d, Ep = _tensoreig(a.array, truncfun = truncfun)
    DTensor(E), DTensor(d), DTensor(collect(Ep))
end

function tensoreig(A::DASTensor{T,N,SYM,CHARGES,SIZES,CHARGE};
                   truncfun = svdtrunc_default) where
                   {T,N,SYM,CHARGES,SIZES,CHARGE}
    N == 2 || throw(ArgumentError("EIG only works on rank 2 tensors"))
    lch = connectingcharge(A)
    ld =  sizes(A,2)[[chargeindex(c,charges(A,2)) for c in lch]]

    E  = DASTensor{T,N}(SYM, (charges(A,1), lch),
            deepcopy.((sizes(A,1), ld)),
            in_out(A), charge(A))
    D  = DASTensor{T,N}(SYM, (lch, lch),
            deepcopy.((ld, ld)),
            vcat(inv(in_out(A,2)), in_out(A,2)))
    Ep = DASTensor{T,N}(SYM, (lch, charges(A,2)),
            deepcopy((ld,sizes(A,2))),
            vcat(inv(in_out(A,2)), in_out(A,2)))

    for (k, degen) in tensor(A)
        in, out = k[1], k[2]
        E[DASSector(in, out)], D[DASSector(out, out)], Ep[DASSector(out, out)], cutoff =
                _tensoreig(degen, truncfun = truncfun, helper = true)
        E.dims[2][chargeindex(out,charges(E,2))] = cutoff
        D.dims[1][chargeindex(out,charges(D,1))] = cutoff
        D.dims[2][chargeindex(out,charges(D,2))] = cutoff
        Ep.dims[1][chargeindex(out,charges(Ep,1))] = cutoff
    end
    return (E, D, Ep)
end


text:
using EffectiveWaves, Plots

pyplot()

Maxtime=100.
T=Float64

"Derivative of Hankel function of the first kind"
function diffhankelh1(n,z)
  if n!=0
    0.5*(hankelh1(-1 + n, z) - hankelh1(1 + n, z))
  else
    - hankelh1(1, z)
  end
end

"Derivative of Bessel function of first kind"
function diffbesselj(n,z)
  if n!=0
    0.5*(besselj(-1 + n, z) - besselj(1 + n, z))
  else
    - besselj(1, z)
  end
end

# background medium
medium = Medium(1.0,1.0+0.0im)
radius_multiplier = 1.005
max_hankel_order=15

incident_medium = medium
θin = 0.2
ω = 15.
k = ω/incident_medium.c

## Strong scatterers
species = [
    Specie(ρ=0.8,r=0.1, c=0.2, volfrac=0.1),
    Specie(ρ=0.2, r=0.2, c=0.1, volfrac=0.1)
]

S = length(species)
# @memoize Z_l_n(l,n) = Zn(ω,species[l],medium,n)

as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
function M(keff,j,l,m,n)
    (n==m ? 1.0+im*0.0:0.0+im*0.0)*(j==l ? 1.0+im*0.0:0.0+im*0.0) + 2.0pi*species[l].num_density*Z_l_n[l,n]*
        kernelN(n-m,k*as[j,l],keff*as[j,l])/(k^2.0-keff^2.0)
end

ho = -1 + sum([ tol .< norm([M(0.9*k + 0.1im,j,l,1,n) for j = 1:S, l = 1:S]) for n=0:max_hankel_order ])

# this matrix is needed to calculate the eigenvectors
MM(keff::Complex{T}) = reshape(
    [M(keff,j,l,m,n) for m in -ho:ho, j = 1:S, n in -ho:ho, l = 1:S]
, ((2ho+1)*S, (2ho+1)*S))
detMM2(keff_vec::Array{T}) = map(x -> real(x*conj(x)), det(MM(keff_vec[1]+im*keff_vec[2])))

function detMM!(F,x)
    F[1] = abs(det(MM(x[1]+im*x[2])))
end

k0 = real(k)
x = k0 .* LinRange(0.0,1.8,100)/10.
y = k0 .* LinRange(0.,0.22,100)

X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))

Z0 = map( (x,y) -> detMM2([x,y]),X,Y)
Z = map( z -> (abs(z)> 2.) ? NaN : z, Z0)
contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M - high ω - low φ")

x = k0 .* LinRange(0.18,0.23,50)
y = k0 .* LinRange(0.01,0.17,50)

x = kx; y = ky;

X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map( (x,y) -> (z = detMM2([x,y]); (abs(z)> 0.5) ? NaN: z),X,Y)

contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M")

x = LinRange(3.,6.,40)
y = LinRange(-1.5,14.,40)

X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map( (x,y) -> (z = detMM2([x,y]); (abs(z)> 1.) ? NaN: z),X,Y)
# Z = map( (x,y) -> detMM2([x,y]),X,Y)

# contour(x,y,Z,fill=true, xlab = "Re k*", ylab = "Im k*", title="Roots of secular det M")
heatmap(x,y,Z, xlab = "Re k*", ylab = "Im k*", title="det M - high ω and almost dirichlet")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Wine is a bit of a pain to build; it requires a two-and-a-half step process;
# we first build the 64-bit version of Wine part-way, use its BUILD DIRECTORY
# (not the result of `make install`, the result of `make`!) to build 32-bit
# Wine, install and bundle that version up, then resume the 64-bit build.

# We also don't support passing explicit targets in
SAFE_ARGS = [a for a in ARGS if startswith(a, "-") && !startswith(a, "--deploy")]
platform_spec = [a for a in ARGS if !startswith(a, "-")]
is_meta = any(a->startswith(a, "--meta-json"), ARGS)
requested_platforms = map(p->parse(Platform, p; validate_strict=true), platform_spec)

version = v"7.0-rc3"
sources = Any[
    GitSource("https://github.com/JuliaComputing/wine-staging.git", "10ff6d27133070cff6369c48d7ee4c56ee8da7aa"),
    DirectorySource("./bundled")
]

# The products that we will ensure are always built
products = Product[
    ExecutableProduct("wine64", :wine64, "wine64/loader"),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("libpng_jll"),
    Dependency("JpegTurbo_jll"),
    Dependency("GnuTLS_jll"),
]

platform32_musl = Platform("i686", "linux"; libc="musl")
platform64_musl = Platform("x86_64", "linux"; libc="musl")

platform32_glibc = Platform("i686", "linux")
platform64_glibc = Platform("x86_64", "linux")

platform64_darwin = Platform("x86_64", "macos")

all_platforms = [
    platform64_musl,
    platform64_glibc,
    platform64_darwin
]

platform_map = Dict(
    platform64_musl => platform32_musl,
    platform64_glibc => platform32_glibc
)

# In meta mode, note all supported platforms.
platforms = all_platforms

if !is_meta
    # TODO: Make it possible to build multiple platforms in one invocation.
    if isempty(requested_platforms)
        @info "No platform specified; defaulting to musl."
        platform64 = platform64_musl
    elseif length(requested_platforms) != 1
        error("This build script is special and only supports one platform at a time.")
    else
        platform64 = requested_platforms[1]
    end

    if platform64 != platform64_darwin
        # No wine32 build for darwin at the moment.

        platform32 = platform_map[platform64]

        wine64_build_script = raw"""
        cd $WORKSPACE/srcdir/wine-staging
        atomic_patch -p1 $WORKSPACE/srcdir/patches/hwcap2.patch

        # First, build wine64, making the actual build directory itself the thing we will install.
        mkdir $WORKSPACE/destdir/wine64
        cd $WORKSPACE/destdir/wine64
        # N.B.: The --build=${target} is technically wrong here, but
        # because we support both glibc and musl exectuables in the
        # build environment, everything goes through ok and doesn't
        # trigger wine's cross compile detection (which would require
        # a more complicated bootstrap process).
        $WORKSPACE/srcdir/wine-staging/configure --prefix=${prefix} --build=${target} --host=${target} --without-x --without-freetype --enable-win64
        make -j${nproc}
        """

        # package up the wine64 build directory itself:
        product_hashes = build_tarballs(copy(SAFE_ARGS), "BuildTmpWine64Build", version, sources, wine64_build_script, [platform64], products, dependencies; skip_audit=true, preferred_gcc_version=v"11.1")

        # Include that tarball as one of the sources we need to include:
        wine64_build_path, wine64_build_hash = product_hashes[platform64]
        push!(sources, ArchiveSource(joinpath("products", wine64_build_path), wine64_build_hash))

        # Next, build wine32, without wine64, then use those tools to build wine32 WITH wine64
        wine32_script = raw"""
        cd $WORKSPACE/srcdir/wine-staging
        atomic_patch -p1 $WORKSPACE/srcdir/patches/hwcap2.patch

        # Next, build wine32, linking against the previously included wine64 stuff:
        mkdir $WORKSPACE/srcdir/wine32_only
        cd $WORKSPACE/srcdir/wine32_only
        $WORKSPACE/srcdir/wine-staging/configure --host=${target} --without-x --without-freetype
        make -j${nproc}

        mkdir $prefix/wine32
        cd $prefix/wine32
        $WORKSPACE/srcdir/wine-staging/configure --prefix=${prefix}/wine32 --host=${target} --without-x --without-freetype --with-wine64=$WORKSPACE/srcdir/wine64 --with-wine-tools=$WORKSPACE/srcdir/wine32_only
        make -j${nproc}
        make -j${nproc} install
        """

        wine32_products = Product[
            ExecutableProduct("wine", :wine, "wine32/loader"),
        ]

        product_hashes = build_tarballs(copy(SAFE_ARGS), "BuildTmpWine32", version, sources, wine32_script, [platform32], wine32_products, dependencies; skip_audit=true, preferred_gcc_version=v"11.1")
        wine32_path, wine32_hash = product_hashes[platform32]
        push!(sources, ArchiveSource(joinpath("products", wine32_path), wine32_hash))
    end

    platforms = [platform64]
end

# Finally, install both:
script = raw"""
cd $WORKSPACE/srcdir/wine-staging
atomic_patch -p1 $WORKSPACE/srcdir/patches/hwcap2.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/darwin.patch

if [[ "${target}" == *darwin* ]]; then
    # On macOS, we currently do not support 32-bit wine. It is unclear whether
    # this is supported in upstream wine (though it is supported in CrossOver,
    # the commercial wine version). Because of that, we can skip the complicated
    # boostrap procedure above and just build everything right here. However,
    # we do have to do a bit of extra work to sucessfully cross-compile for
    # darwin, since we can't rely on natively running the tools.

    # Go to wine_tools directory
    mkdir $WORKSPACE/wine_tools
    cd $WORKSPACE/wine_tools

    # winebuild will use {nm, ld, as} indiscriminately and expect those to be
    # the host tools, so put those on the PATH.
    mkdir aliases
	ln -s $(which x86_64-linux-musl-nm) aliases/nm
	ln -s $(which x86_64-linux-musl-ld) aliases/ld
	ln -s $(which x86_64-linux-musl-as) aliases/as
    export OLD_PATH=$PATH
    export PATH=$PWD/aliases:$PATH
    ${WORKSPACE}/srcdir/wine-staging/configure --build=x86_64-linux-musl --host=x86_64-linux-musl CC=x86_64-linux-musl-gcc LD=x86_64-linux-musl-ld AS=x86_64-linux-musl-as --enable-win64 --without-x --without-freetype
    make -j${nproc} CC=x86_64-linux-musl-gcc LD=x86_64-linux-musl-ld AS=x86_64-linux-musl-as
    export PATH=$OLD_PATH

    # Ok, now we can actually do the proper build of wine.
    cd $WORKSPACE/srcdir/wine-staging
    ./configure --build=${MACHTYPE} --prefix=${prefix} --host=${target} --without-x --without-freetype --with-wine-tools=$WORKSPACE/wine_tools --enable-win64
    make -j${nproc}
else
    cp -r $WORKSPACE/srcdir/wine32/* ${prefix}/
    cd $WORKSPACE/srcdir/wine64
fi

make -j${nproc} install
install_license $WORKSPACE/srcdir/wine-staging/LICENSE
install_license $WORKSPACE/srcdir/wine-staging/COPYING.LIB
"""

final_products = Product[
    ExecutableProduct("wine64", :wine64),
# Currently the macos build is wine64 only, so let's leave it at that.
#    ExecutableProduct("wine", :wine),
]

build_tarballs(ARGS, "Wine", version, sources, script, platforms, final_products, dependencies, preferred_gcc_version=v"11.1")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ---
# title: 1648. Sell Diminishing-Valued Colored Balls
# id: problem1648
# author: AquaIndigo
# date: 2020-11-08
# difficulty: Medium
# categories: Math, Greedy, Sort
# link: <https://leetcode.com/problems/sell-diminishing-valued-colored-balls/description/>
# hidden: true
# ---
#
# You have an `inventory` of different colored balls, and there is a customer
# that wants `orders` balls of **any** color.
#
# The customer weirdly values the colored balls. Each colored ball's value is
# the number of balls **of that color  **you currently have in your `inventory`.
# For example, if you own `6` yellow balls, the customer would pay `6` for the
# first yellow ball. After the transaction, there are only `5` yellow balls
# left, so the next yellow ball is then valued at `5` (i.e., the value of the
# balls decreases as you sell more to the customer).
#
# You are given an integer array, `inventory`, where `inventory[i]` represents
# the number of balls of the `ith` color that you initially own. You are also
# given an integer `orders`, which represents the total number of balls that the
# customer wants. You can sell the balls **in any order**.
#
# Return _the **maximum** total value that you can attain after selling
# _`orders` _colored balls_. As the answer may be too large, return it
# **modulo**`109 + 7`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/jj.gif)
#
#
#
#     Input: inventory = [2,5], orders = 4
#     Output: 14
#     Explanation: Sell the 1st color 1 time (2) and the 2nd color 3 times (5 + 4 + 3).
#     The maximum total value is 2 + 5 + 4 + 3 = 14.
#
#
# **Example 2:**
#
#
#
#     Input: inventory = [3,5], orders = 6
#     Output: 19
#     Explanation: Sell the 1st color 2 times (3 + 2) and the 2nd color 4 times (5 + 4 + 3 + 2).
#     The maximum total value is 3 + 2 + 5 + 4 + 3 + 2 = 19.
#
#
# **Example 3:**
#
#
#
#     Input: inventory = [2,8,4,10,6], orders = 20
#     Output: 110
#
#
# **Example 4:**
#
#
#
#     Input: inventory = [1000000000], orders = 1000000000
#     Output: 21
#     Explanation: Sell the 1st color 1000000000 times for a total value of 500000000500000000. 500000000500000000 modulo 109 + 7 = 21.
#
#
#
#
# **Constraints:**
#
#   * `1 <= inventory.length <= 105`
#   * `1 <= inventory[i] <= 109`
#   * `1 <= orders <= min(sum(inventory[i]), 109)`
#
#
## @lc code=start
using LeetCode

function max_profit(inventory::Vector{Int}, orders::Int)
    rest = sum(inventory) - orders
    sort!(inventory)
    len = length(inventory)
    for n in inventory
        if rest ÷ len >= n
            len -= 1
            rest -= n
        else
            break
        end
    end
    q = rest ÷ len
    r = rest - q * len
    ret = 0
    for i in (length(inventory) - len + 1):length(inventory)
        ret += sum((q + 2):inventory[i])
    end
    ret += (q + 1) * (len - r)
    return ret % (10^9 + 7)
end
## @lc code=end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# qcp.jl
#
# A simple quadratically constrained probgram
# Based on http://www.gurobi.com/documentation/5.5/example-tour/node25
#############################################################################

using JuMP

# Will require either Gurobi.jl, CPLEX.jl, or Mosek.jl to run
m = Model()

# Need nonnegativity for (rotated) second-order cone
@variable(m, x)
@variable(m, y >= 0)
@variable(m, z >= 0)

# Maximize x
@objective(m, Max, x)

# Subject to 1 linear and 2 nonlinear constraints
@constraint(m, x + y + z == 1)
@constraint(m, x*x + y*y - z*z <= 0)
@constraint(m, x*x - y*z <= 0)

# Print the model to check correctness
print(m)

# Solve with Gurobi
status = solve(m)

# Solution
println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))


text:
function mgc2mgc!(dst_ceps::AbstractVector, dst_α, dst_γ,
                  src_ceps::AbstractVector, src_α, src_γ)
    dst_order = length(dst_ceps) - 1
    α = (dst_α-src_α) / (1.0-dst_α*src_α)

    if α == zero(α)
        dst_ceps = gnorm(src_ceps, src_γ)
        dst_ceps = gc2gc(dst_ceps, src_γ, dst_order, dst_γ)
        ignorm!(dst_ceps, dst_γ)
    else
        freqt!(dst_ceps, src_ceps, α)
        gnorm!(dst_ceps, src_γ)
        dst_ceps = gc2gc(dst_ceps, src_γ, dst_order, dst_γ)
        ignorm!(dst_ceps, dst_γ)
    end

    dst_ceps
end

function mgc2mgc(src_ceps::AbstractVector,
                 src_α=0.0,
                 src_γ=0.0,
                 dst_order=length(src_ceps) - 1,
                 dst_α=0.0,
                 dst_γ=0.0)
    mgc2mgc!(Array(eltype(src_ceps), dst_order+1), dst_α, dst_γ,src_ceps, src_α, src_γ)
end

function mgc2mgc{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T},
                                            dst_order=param_order(paramdef(state)),
                                            dst_α=allpass_alpha(paramdef(state)),
                                            dst_γ=glog_gamma(paramdef(state)))
    assert_not_ready_to_filt(state)

    def = paramdef(state)
    src_α = allpass_alpha(def)
    src_γ = glog_gamma(def)
    normalized = (dst_γ == zero(dst_γ))
    data = rawdata(state)
    newdef = MelGeneralizedCepstrum(dst_order, dst_α, dst_γ)
    SpectralParamState(newdef, mgc2mgc(data, src_α, src_γ, dst_order, dst_α, dst_γ),
                       has_loggain(state), normalized)
end


text:
using PyCall, Conda

using Conda
Conda.add_channel("conda-forge")
Conda.add("quadpy")

cmd = `pip3 install quadpy --user`
run(cmd)

quadpy = pyimport("quadpy")


text:
function show(io::IO, object::testOutput)
    if object.method != nothing
        @printf("%s\n", object.method)
    end
    if object.h0 != nothing
        @printf(" Alternative hypothesis:   % s\n", object.h0)
    end
    if object.n != nothing
        @printf(" Sample size:         ")
        for i = 1:length(object.n)
            @printf("% d  ", object.n[i])
        end
        @printf("\n")
    end
    if object.df != nothing
            @printf(" Degrees of freedom:  ")
        for i = 1:length(object.df)
            if isa(object.df[i], Integer)
                @printf(" %d  ", object.df[i])
            elseif isa(object.df, Number)
                @printf("% .2f  ", object.df[i])
            end
        end
       @printf("\n")
    end
    if object.estimate != nothing
        @printf(" Estimate:            % .6f\n", object.estimate)
    end
    if object.se != nothing
        @printf(" Standard error:      % .6f\n", object.se)
    end
    if object.statistic != nothing
        if ndims(object.statistic)==0
            @printf(" Statistic:           % .6f\n", object.statistic)
        else ndims(object.statistic)==2
            @printf(" Statistic:            % 7s  % 7s  % 9s\n", "Group 1", "Group 2", "Statistic")
            for i = 1:size(object.statistic, 1)
                @printf("                       % 7d  % 7d  %  7.6f\n",
                        object.statistic[i,1], object.statistic[i,2], object.statistic[i,3])
            end
        end
    end
    if object.crit != nothing
        @printf(" Critical value:      % .6f\n", object.crit)
    end
    if object.ci != nothing
        if isa(object.ci, Integer)
            @printf(" Confidence interval: % d      % d\n", object.ci[1], object.ci[2])
        else
            @printf(" Confidence interval: % .6f      % .6f\n", object.ci[1], object.ci[2])
        end
    end
    if object.p != nothing
        if object.p < 10e-16
            @printf(" p value:             < 10e-16\n")
        elseif object.p < 1
            @printf(" p value:             % .6f\n", object.p)
        else
            @printf(" p value:              1.0")
        end
    end
end

function show(io::IO, object::outOutput)
    if object.method!=nothing
        println(io, "$(object.method)")
    end
    @printf("Outlier ID:         ")
    if length(object.outid) == 0
        @printf("\n")
    else
        for i = 1:length(object.outid)
            if i < length(object.outid)
                @printf("%d, ", object.outid[i])
            else
                @printf("%d\n", object.outid[i])
            end
        end
    end
    @printf("Outlier value:      ")
    if length(object.outid) == 0
        @printf("\n")
    else
        if isa(object.outval, Integer)
            for i = 1:length(object.outval)
                if i < length(object.outval)
                    @printf("%d, ", object.outval[i])
                else
                    @printf("%d\n", object.outval[i])
                end
            end
        else
            for i=1:length(object.outval)
                if i < length(object.outval)
                    @printf("%.5f, ", object.outval[i])
                else
                    @printf("%.5f\n", object.outval[i])
                end
            end
        end
    end
    @printf("Number of outliers: %d\n", object.nout)
    @printf("Non-outlier ID:     ")
    for i=1:length(object.keepid)
        if i < length(object.keepid)
            @printf("%d, ", object.keepid[i])
        else
            @printf("%d\n", object.keepid[i])
        end
    end
end

function show(io::IO, object::binomciOutput)
   @printf(" p_hat:               %.4f\n", object.p_hat)
   @printf(" confidence interval: %.4f   %.4f\n", object.confint[1], object.confint[2])
   @printf(" Sample size          %d\n", object.n)
end


text:
using Test, IonSim

@testset "lasers -- Laser" begin
    # should not be able to set unnormalized k-vectors or polarizations
    @test_throws AssertionError Laser(ϵ=(x=2, y=0, z=0))
    @test_throws AssertionError Laser(k=(x=2, y=0, z=0))
    # should not be able to initialize with non-orthogonal ϵ̂, k̂
    @test_throws AssertionError Laser(ϵ=(x=1, y=0, z=0), k=(x=1, y=0, z=0))
    # pointing cannot be overspecified
    @test_throws AssertionError Laser(pointing=[(1, 0.5), (1, 0.5)])
    # pointing magnitude must be within [0,1]
    @test_throws AssertionError Laser(pointing=[(1, 2.0)])

    # should be able to set L.E, L.ϕ to a constant Real value
    t = 0:1:100
    L = Laser(E=1, ϕ=1)
    ones = [1 for _ in t]
    # should be able to set L.E, L.ϕ to a function of time
    @test L.E.(t) == ones
    @test L.ϕ.(t) == ones
    L = Laser(E=sin, ϕ=sin)
    @test L.E.(t) == sin.(t)
    @test L.ϕ.(t) == sin.(t)

    # test that normalization is enforced for altered L.ϵ/L.k
    @test_throws AssertionError L.ϵ = (x=2, y=0, z=0)
    @test_throws AssertionError L.k = (x=2, y=0, z=0)
    # test that pointing constraints are enforced for altered values
    @test_throws AssertionError L.pointing = [(1, 0.5), (1, 0.5)]
    @test_throws AssertionError L.pointing = [(1, 2.0)]
    L.λ = 1
    @test L.λ == 1

    # just make sure print(L) isn't broken
    @suppress print(L)

    # test non-orthogonal warnings
    @test_logs (:warn, "!(ϵ ⟂ k)") L.ϵ = (x̂+ŷ+ẑ)/√3
    @test_logs (:warn, "!(ϵ ⟂ k)") L.k = (x̂+ŷ+ẑ)/√3

    # test comparison
    L2 = copy(L)
    @test L == L2
    L.E = 7
    @test L ≠ L2
end


text:
using PlotlyJS

function datetimestrings()
    x = ["2013-10-04 22:23:00", "2013-11-04 22:23:00", "2013-12-04 22:23:00"]
    plot(scatter(x=x, y=[1 ,3, 6]))
end


text:
using StatsBase: mean
using .NonadiabaticDistributions: NonadiabaticDistribution
using NQCDynamics: masses

export FSSH

"""
    FSSH{T} <: SurfaceHopping

Type for fewest-switches surface hopping

```jldoctest
Simulation{FSSH}(Atoms(:H), Free())

# output

Simulation{FSSH{Float64}}:
  Atoms{1, Float64}([:H], UInt8[0x01], [1837.4715941070515])
  Free(1)
```
"""
mutable struct FSSH{T} <: SurfaceHopping
    density_propagator::Matrix{Complex{T}}
    hopping_probability::Vector{T}
    state::Int
    new_state::Int
    rescaling::Symbol
    function FSSH{T}(states::Integer, rescaling::Symbol) where {T}
        density_propagator = zeros(states, states)
        hopping_probability = zeros(states)
        new{T}(density_propagator, hopping_probability, 0, 0, rescaling)
    end
end

function Simulation{FSSH}(atoms::Atoms{S,T}, model::Model; rescaling=:standard, kwargs...) where {S,T}
    Simulation(atoms, model, FSSH{T}(NQCModels.nstates(model), rescaling); kwargs...)
end

function DynamicsMethods.DynamicsVariables(
    sim::AbstractSimulation{<:SurfaceHopping}, v, r, electronic::NonadiabaticDistribution
)
    σ = NonadiabaticDistributions.initialise_adiabatic_density_matrix(electronic, sim.calculator, r)
    state = sample(Weights(diag(real.(σ))))
    return SurfaceHoppingVariables(ComponentVector(v=v, r=r, σreal=σ, σimag=zero(σ)), state)
end

function acceleration!(dv, v, r, sim::AbstractSimulation{<:FSSH}, t, state)
    for I in eachindex(dv)
        dv[I] = -sim.calculator.adiabatic_derivative[I][state, state]
    end
    DynamicsUtils.divide_by_mass!(dv, masses(sim))
    return nothing
end

"""
    evaluate_hopping_probability!(sim::Simulation{<:FSSH}, u, dt)

Evaluates the probability of hopping from the current state to all other states

# Implementation

- `σ` is Hermitan so the choice `σ[m,s]` or `σ[s,m]` is irrelevant; we take the real part.
- 'd' is skew-symmetric so here the indices are important.
"""
function evaluate_hopping_probability!(sim::AbstractSimulation{<:FSSH}, u, dt)
    v = get_hopping_velocity(sim, u)
    σ = DynamicsUtils.get_quantum_subsystem(u)
    s = sim.method.state
    d = get_hopping_nonadiabatic_coupling(sim)

    fewest_switches_probability!(sim.method.hopping_probability, v, σ, s, d, dt)
end

function get_hopping_nonadiabatic_coupling(sim::Simulation{<:SurfaceHopping})
    sim.calculator.nonadiabatic_coupling
end

function get_hopping_velocity(::Simulation{<:SurfaceHopping}, u)
    DynamicsUtils.get_velocities(u)
end

function get_hopping_eigenvalues(sim::Simulation{<:SurfaceHopping})
    sim.calculator.eigen.values
end

function fewest_switches_probability!(probability, v, σ, s, d, dt)
    probability .= 0 # Set all entries to 0
    for m in axes(σ, 1)
        if m != s
            for I in eachindex(v)
                probability[m] += 2v[I]*real(σ[m,s]/σ[s,s])*d[I][s,m] * dt
            end
        end
    end

    clamp!(probability, 0, 1) # Restrict probabilities between 0 and 1
    cumsum!(probability, probability)
end

function select_new_state(sim::AbstractSimulation{<:FSSH}, u)

    random_number = rand()
    for (i, prob) in enumerate(sim.method.hopping_probability)
        if i != sim.method.state # Avoid self-hops
            if prob > random_number
                return i
            end
        end
    end
    return sim.method.state
end

"""
    rescale_velocity!(sim::AbstractSimulation{<:FSSH}, u)::Bool

Rescale the velocity in the direction of the nonadiabatic coupling.

# References

[HammesSchiffer1994](@cite)
"""
function rescale_velocity!(sim::AbstractSimulation{<:FSSH}, u)::Bool
    sim.method.rescaling === :off && return true

    new_state, old_state = unpack_states(sim)
    velocity = get_hopping_velocity(sim, u)
    eigs = get_hopping_eigenvalues(sim)

    d = extract_nonadiabatic_coupling(get_hopping_nonadiabatic_coupling(sim), new_state, old_state)
    a = calculate_a(sim, d)
    b = calculate_b(d, velocity)
    c = calculate_potential_energy_change(eigs, new_state, old_state)

    discriminant = b^2 - 4a * c
    discriminant < 0 && return false # Frustrated hop with insufficient kinetic energy

    root = sqrt(discriminant)
    if b < 0
        γ = (b + root) / 2a
    else
        γ = (b - root) / 2a
    end
    perform_rescaling!(sim, DynamicsUtils.get_velocities(u), γ, d)

    return true
end

"""
    unpack_states(sim)

Get the two states that we are hopping between.
"""
function unpack_states(sim::AbstractSimulation{<:FSSH})
    return (sim.method.new_state, sim.method.state)
end

"""
    extract_nonadiabatic_coupling(coupling, new_state, old_state)

Extract the nonadiabatic coupling vector between states `new_state` and `old_state`
"""
function extract_nonadiabatic_coupling(coupling, new_state, old_state)
    [coupling[I][new_state, old_state] for I ∈ CartesianIndices(coupling)]
end

"""
    calculate_a(sim::AbstractSimulation{<:SurfaceHopping}, coupling::AbstractMatrix)

Equation 40 from [HammesSchiffer1994](@cite).
"""
function calculate_a(sim::AbstractSimulation{<:SurfaceHopping}, coupling::AbstractMatrix)
    a = zero(eltype(coupling))
    for I in CartesianIndices(coupling)
        a += coupling[I]^2 / masses(sim, I)
    end
    return a / 2
end

"""
    calculate_b(coupling::AbstractMatrix, velocity::AbstractMatrix)

Equation 41 from [HammesSchiffer1994](@cite).
"""
function calculate_b(coupling::AbstractMatrix, velocity::AbstractMatrix)
    return dot(coupling, velocity)
end

"""
    perform_rescaling!(
        sim::AbstractSimulation{<:SurfaceHopping}, velocity, velocity_rescale, d
    )

Equation 33 from [HammesSchiffer1994](@cite).
"""
function perform_rescaling!(
    sim::AbstractSimulation{<:SurfaceHopping}, velocity, γ, d
)
    for I in CartesianIndices(d)
        velocity[I] -= γ * d[I] / masses(sim, I)
    end
    return nothing
end

function calculate_potential_energy_change(eigs::AbstractVector, new_state::Integer, current_state::Integer)
    return eigs[new_state] - eigs[current_state]
end

function Estimators.diabatic_population(sim::AbstractSimulation{<:FSSH}, u)
    r = DynamicsUtils.get_positions(u)
    U = NonadiabaticDistributions.evaluate_transformation(sim.calculator, r)

    σ = copy(DynamicsUtils.get_quantum_subsystem(u).re)
    σ[diagind(σ)] .= 0
    σ[u.state, u.state] = 1

    return diag(U * σ * U')
end

function Estimators.adiabatic_population(sim::AbstractSimulation{<:FSSH}, u)
    population = zeros(NQCModels.nstates(sim.calculator.model))
    population[u.state] = 1
    return population
end

function DynamicsUtils.classical_hamiltonian(sim::Simulation{<:FSSH}, u)
    kinetic = DynamicsUtils.classical_kinetic_energy(sim, DynamicsUtils.get_velocities(u))
    Calculators.evaluate_potential!(sim.calculator, DynamicsUtils.get_positions(u))
    Calculators.eigen!(sim.calculator)
    potential = sim.calculator.eigen.values[sim.method.state]
    return kinetic + potential
end


text:
using .DataFrames

"""

# DataFrame

DataFrame constructor from a Chains object.
Returns either a DataFrame or an Array{DataFrame}

### Method
```julia
  DataFrame(
    chn::Chains,
    sections::Vector{Symbo);
    append_chains::Bool,
    remove_missing_union::Bool
  )
```

### Required arguments
```julia
* `chn` : Chains object to convert to an DataFrame
```

### Optional arguments
```julia
* `sections = Symbol[]` : Sections from the Chains object to be included
* `append_chains = true`  : Append chains into a single column
* `remove_missing_union = true`  : Remove Union{Missing, Real}
```

### Examples
```julia
* `DataFrame(chns)` : DataFrame with chain values are appended
* `DataFrame(chns[:par])` : DataFrame with single parameter (chain values are appended)
* `DataFrame(chns, [:parameters])`  : DataFrame with only :parameter section
* `DataFrame(chns, [:parameters, :internals])`  : DataFrame includes both sections
* `DataFrame(chns, append_chains=false)`  : Array of DataFrames, each chain in its own array
* `DataFrame(chns, remove_missing_union=false)` : No conversion to remove missing values
```

"""
function DataFrames.DataFrame(chain::Chains,
    sections::Union{Symbol, Vector{Symbol}}=Symbol[:parameters];
    append_chains=true,
    remove_missing_union=true,
    sorted=false,
    showall=false
)
    sections = _clean_sections(chain, sections)
    sections = sections isa AbstractArray ? sections : [sections]
    sections = showall ? [] : sections
    section_list = length(sections) == 0 ? sort_sections(chain) : sections

    if remove_missing_union
        chn = concretize(Chains(chain, section_list))
    else
        chn = Chains(chain, section_list)
    end
    d, p, c = size(chn.value.data)

    local b
    if append_chains
        b = DataFrame()
        for section in section_list
            names = sorted ?
                sort(chn.name_map[section],
                    by=x->string(x), lt = natural) :
                chn.name_map[section]
            for par in names
                x = get(chn, Symbol(par))
                d, c = size(x[Symbol(par)])
                b = hcat(b, DataFrame(Symbol(par) => reshape(x[Symbol(par)], d*c)[:, 1]))
            end
        end
    else
        b = Vector{DataFrame}(undef, c)
        for i in 1:c
            b[i] = DataFrame()
            for section in section_list
                names = sorted ?
                    sort(chn.name_map[section],
                        by=x->string(x), lt = natural) :
                    chn.name_map[section]
                for par in names
                    x = get(chn[:,:,i], Symbol(par))
                    d, c = size(x[Symbol(par)])
                    b[i] = hcat(b[i], DataFrame(Symbol(par) => x[Symbol(par)][:,1]))
                end
            end
        end
    end
    return b
end

DataFrames.DataFrame(cdfs::Vector{<:ChainDataFrame}) = map(DataFrame, cdfs)
function DataFrames.DataFrame(cdf::ChainDataFrame)
    colnames = collect(keys(cdf.nt))
    cols = collect(values(cdf.nt))

    return DataFrame(cols, colnames)
end

function concretize(df::AbstractDataFrame)
    if all(isconcretetype_recursive ∘ typeof, eachcol(df))
        df
    else
        mapcols(concretize, df)
    end
end


text:
# (1) パッケージをインストールする（初回のみ）
## ] + addで

# (2) ライブラリの読み出し
using DataFrames
using DataFramesMeta
using GLM
using RCall
@rlibrary readr
@rimport tidyr
@rlibrary broom

# (3) データの読み込み
email_data =
  rcopy(R"readr::read_csv('http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv')")

# (4) 女性向けメールが配信されたデータを削除したデータを作成
@linq male_df =
  email_data |> where(:segment .!= "Womens E-Mail") |> # 女性向けメールが配信されたデータを削除
  transform(treatment = ifelse.(:segment .== "Mens E-Mail", 1, 0)) #介入を表すtreatment変数を追加

# (5) セレクションバイアスのあるデータを作成
## seedを固定
R"set.seed(1)"

## 条件に反応するサンプルの量を半分にする
obs_rate_c = 0.5
obs_rate_t = 0.5

## バイアスのあるデータを作成
@linq biased_data =
  male_df |> transform(
    obs_rate_c = ifelse.(
      (:history .> 300) .| (:recency .< 6) .| (:channel .== "Multichannel"),
      obs_rate_c,
      1,
    ),
    obs_rate_t = ifelse.(
      (:history .> 300) .| (:recency .< 6) .| (:channel .== "Multichannel"),
      1,
      obs_rate_t,
    ),
    random_number = rcopy(R"runif($(nrow(male_df)))"), # 乱数を合わせるためRに作らせる
  ) |> where(
    ((:treatment .== 0) .& (:random_number .< :obs_rate_c)) .|
    ((:treatment .== 1) .& (:random_number .< :obs_rate_t)),
  )

# (6) バイアスのあるデータでの回帰分析
## 回帰分析の実行
## JuliaのGLM.jlの結果をRobjectに変換することは現状不可能。逆はOrderedDictになる
biased_reg = lm(@formula(spend ~ treatment + history), biased_data)
# robject(biased_reg) # エラーになる
biased_reg_R = R"summary(lm(formula = spend ~ treatment + history, $biased_data))"
rcopy(biased_reg_R) # OrderedDictになる

## 分析結果のレポート
## lmのsummaryに比べ情報量が少ない。
println(biased_reg)
println(R"summary($biased_reg_R)")

## 推定されたパラメーターの取り出し
## coeftableでできるがbroomのようにデータフレームにはしてくれない。
biased_reg_coef = coeftable(biased_reg)

biased_reg_coef_R = rcopy(R"broom::tidy($biased_reg_R)")

## broom::tidyと同じ働きをする関数を自作
function broom_tidy_julia(model)
  table = coeftable(model)
  return DataFrame(
    term = coefnames(model),
    estimate = table.cols[1],
    stderror = table.cols[2],
    statistic = table.cols[3],
    pvalue = table.cols[4],
  )
end

broom_tidy_julia(biased_reg)

# (7) RCTデータでの回帰分析とバイアスのあるデータでの回帰分析の比較
## RCTデータでの単回帰
rct_reg = lm(@formula(spend ~ treatment), male_df)
rct_reg_coef = broom_tidy_julia(rct_reg)

## バイアスのあるデータでの単回帰
nonrct_reg = lm(@formula(spend ~ treatment), biased_data)
nonrct_reg_coef = broom_tidy_julia(nonrct_reg)

## バイアスのあるデータでの重回帰
nonrct_mreg =
  lm(@formula(spend ~ treatment + recency + channel + history), biased_data)
nonrct_mreg_coef = broom_tidy_julia(nonrct_mreg)

# (8) OVBの確認
## (a) history抜きの回帰分析とパラメーターの取り出し
short_coef =
  biased_data |>
  data ->
    lm(@formula(spend ~ treatment + recency + channel), data) |>
    broom_tidy_julia

## aの結果から介入効果に関するパラメーターのみを取り出す
## pullがないので少し面倒
@linq alpha_1_p =
  short_coef |> where(:term .== "treatment") |> select(:estimate)
alpha_1 = alpha_1_p[1, 1]

## (b) historyを追加した回帰分析とパラメーターの取り出し
long_coef =
  biased_data |>
  data ->
    lm(@formula(spend ~ treatment + recency + channel + history), data) |>
    broom_tidy_julia

## bの結果から介入とhistoryに関するパラメーターを取り出す
@linq beta_1_p = long_coef |> where(:term .== "treatment") |> select(:estimate)
beta_1 = beta_1_p[1, 1]
@linq beta_2_p = long_coef |> where(:term .== "history") |> select(:estimate)
beta_2 = beta_2_p[1, 1]

## (c) 脱落した変数と介入変数での回帰分析
omitted_coef =
  biased_data |>
  data ->
    lm(@formula(history ~ treatment + channel + recency), data) |>
    broom_tidy_julia
## cの結果から介入変数に関するパラメーターを取り出す
@linq gamma_1_p =
  omitted_coef |> where(:term .== "treatment") |> select(:estimate)
gamma_1 = gamma_1_p[1, 1]

## OVBの確認
beta_2 * gamma_1
alpha_1 - beta_1

# (9) OVBの確認(broomを利用した場合)
## formulaをデータフレーム化しまとめて処理してみる
## モデル式のベクトルを用意
formula_vec = [
  @formula(spend ~ treatment + recency + channel), # モデルA
  @formula(spend ~ treatment + recency + channel + history), # モデルB
  @formula(history ~ treatment + channel + recency)
] # モデルC

## formulaに名前を付ける
names_formula_vec = ["reg" * "_" * letter for letter = 'A':'C']

## モデル式のデータフレーム化
models = DataFrame(model_index = names_formula_vec, formula = formula_vec)

## まとめて回帰分析を実行
@linq df_models =
  models |>
  transform(model = map(formula -> lm(formula, biased_data), :formula)) |>
  transform(lm_result = map(model -> broom_tidy_julia(model), :model))

## モデルの結果を整形
## formulaは元から文字列のように表示されるのでそのまま
## unnestはJuliaにはないのでRCallで
@linq df_results =
  df_models |> select(:formula, :model_index, :lm_result) |>
  tidyr.unnest(cols = :lm_result) |> rcopy()

## モデルA,B,Cでのtreatmentのパラメータを抜き出す
@linq treatment_coef =
  df_results |> where(:term .== "treatment") |> select(:estimate)

## モデルBからhistoryのパラメータを抜き出す
@linq history_coef =
  df_results |> where(:model_index .== "reg_B", :term .== "history") |>
  select(:estimate)

## OVBの確認
OVB = history_coef[1, 1] * treatment_coef[3, 1]
coef_gap = treatment_coef[1, 1] - treatment_coef[2, 1]
OVB # beta_2*gamma_1
coef_gap # alpha_1 - beta_1

# (10) 入れてはいけない変数を入れてみる
#visitとtreatmentとの相関
cor_visit_treatment =
  lm(@formula(treatment ~ visit + channel + recency + history), biased_data) |>
  broom_tidy_julia

# visitを入れた回帰分析を実行
bad_control_reg =
  lm(
    @formula(spend ~ treatment + channel + recency + history + visit),
    biased_data,
  ) |> broom_tidy_julia


text:
module OptimalTrees
    using Clustering
    using DataFrames
    using DocStringExtensions
    using JuMP
    using MathOptInterface
    using Missings
    using Parameters
    using ProgressMeter
    using Random

    include("helpers.jl")
    include("binarynode.jl")
    include("tree.jl")
    include("tree_utils.jl")
    include("training.jl")
    include("ensemble.jl")

    const MOI = MathOptInterface
    const OPTIMALTREES_ROOT = dirname(dirname(@__FILE__))
    const DATA_DIR = OPTIMALTREES_ROOT * "\\data\\"

    # Structs
    export MIOTree, BinaryNode, TreeEnsemble,

    # Tree building
        MIOTree_defaults,
        leftchild, rightchild, children,
        get_lower_child, get_upper_child, get_parent,
        printnode, generate_binary_tree,
        delete_children!,
        populate_nodes!, prune!,
        chop_down!,

    # More advanced tree building
        deepen_to_max_depth!,
        deepen_one_level!, clone,

    # Tree querying
        is_leaf, depth,
        lineage, alloffspring,
        get_classification_label,
        set_classification_label!,
        get_split_values, set_split_values!,
        get_split_threshold, get_split_weights,
        get_regression_constant, get_regression_weights,
        apply, predict,

    # Tree training
        fit!,
        SVM, hyperplane_cart,
        warmstart, generate_MIO_model,
        sequential_train!,

    # Tree utilities
        find_leaves, check_if_trained,
        trust_region_data, pwl_constraint_data,

    # Scoring functions
        score,

    # Helper functions
        set_param!, set_params!, get_param,
        normalize, pairwise_distances, mode,
        allnodes, allleaves,
        clean_model!,
        split_data,

    # Tree ensembles
        TreeEnsemble_defaults,
        plant_trees, weigh_trees,

    # Debugging
        debug_if_trained
end


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Performance testing for the parsing modules.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-12

include("../hdr.jl")

require("structs.jl")
require("conllx.jl")
require("dep.jl")
require("hybrid.jl")

using Structs
using CoNLLX
using DepGraph
using Hybrid

# Featurise along the gold path and measure the time consumption.
open(talpath, "r") do data_f
    sents = collect(Sentence, conllxparse(data_f, useproj=true))
    coder = Coder()

    fidbyfstruct = Identifier(FeatStruct, Uint)
    featids = Array(Uint, NUMFEATS)
    numcalls = 0
    duration = 0
    for sent in sents
        conf = config(sent, coder)
        while !isterminal(conf)
            before = time()
            featurise!(featids, fidbyfstruct, conf)
            duration += time() - before
            numcalls += 1

            trans = oracle(conf)
            apply!(conf, trans)
        end
    end

    throughput = numcalls / duration
    @printf("Featurisation throughput: %.1f configuration(s) per second\n",
        throughput)
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
using Test
using LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators
import ClimaCore.Utilities: half

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

# 3D deformation flow (DCMIP 2012 Test 1-1)
# Reference: http://www-personal.umich.edu/~cjablono/DCMIP-2012_TestCaseDocument_v1.7.pdf, Section 1.1

const R = 6.37122e6 # radius
const grav = 9.8 # gravitational constant
const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const z_top = 1.2e4 # height position of the model top
const p_top = 25494.4 # pressure at the model top
const T_0 = 300 # isothermal atmospheric temperature
const H = R_d * T_0 / grav # scale height
const p_0 = 1.0e5 # reference pressure
const τ = 1036800.0 # period of motion
const ω_0 = 23000 * pi / τ # maxium of the vertical pressure velocity
const b = 0.2 # normalized pressure depth of divergent layer
const λ_c1 = 150.0 # initial longitude of first tracer
const λ_c2 = 210.0 # initial longitude of second tracer
const ϕ_c = 0.0 # initial latitude of tracers
const z_c = 5.0e3 # initial altitude of tracers
const R_t = R / 2 # horizontal half-width of tracers
const Z_t = 1000.0 # vertical half-width of tracers

# set up function space
function sphere_3D(
    R = 6.37122e6,
    zlim = (0, 12.0e3),
    helem = 6,
    zelem = 12,
    npoly = 3,
)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_names = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = zelem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.SphereDomain(R)
    horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
    horztopology = Topologies.Topology2D(horzmesh)
    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

# set up 3D domain
hv_center_space, hv_face_space = sphere_3D()

# initial conditions
coords = Fields.coordinate_field(hv_center_space)
face_coords = Fields.coordinate_field(hv_face_space)

r1(λ, ϕ) = R * acos(sind(ϕ_c) * sind(ϕ) + cosd(ϕ_c) * cosd(ϕ) * cosd(λ - λ_c1))
r2(λ, ϕ) = R * acos(sind(ϕ_c) * sind(ϕ) + cosd(ϕ_c) * cosd(ϕ) * cosd(λ - λ_c2))

p(z) = p_0 * exp(-z / H)
ρ(z) = p(z) / R_d / T_0

y0 = map(coords) do coord
    z = coord.z
    zd = z - z_c
    λ = coord.long
    ϕ = coord.lat
    rd1 = r1(λ, ϕ)
    rd2 = r2(λ, ϕ)

    d1 = min(1, (rd1 / R_t)^2 + (zd / Z_t)^2)
    d2 = min(1, (rd2 / R_t)^2 + (zd / Z_t)^2)

    q1 = 1 / 2.0 * (1 + cos(pi * d1)) + 1 / 2.0 * (1 + cos(pi * d2))
    q2 = 0.9 - 0.8 * q1^2
    q3 = 0.0
    if d1 < 0.5 || d2 < 0.5
        q3 = 1.0
    else
        q3 = 0.1
    end

    if (z > z_c) && (ϕ > ϕ_c - rad2deg(1 / 8)) && (ϕ < ϕ_c + rad2deg(1 / 8))
        q3 = 0.1
    end
    q4 = 1 - 3 / 10 * (q1 + q2 + q3)

    ρq1 = ρ(z) * q1
    ρq2 = ρ(z) * q2
    ρq3 = ρ(z) * q3
    ρq4 = ρ(z) * q4

    return (ρq1 = ρq1, ρq2 = ρq2, ρq3 = ρq3, ρq4 = ρq4)
end

function rhs!(dydt, y, (coords, face_coords), t)

    ϕ = coords.lat
    λ = coords.long
    zc = coords.z
    zf = face_coords.z
    λp = λ .- 360 * t / τ
    k = 10 * R / τ

    ϕf = face_coords.lat
    λf = face_coords.long
    λpf = λf .- 360 * t / τ

    sp =
        @. 1 + exp((p_top - p_0) / b / p_top) - exp((p(zf) - p_0) / b / p_top) -
           exp((p_top - p(zf)) / b / p_top)
    ua = @. k * sind(λp)^2 * sind(2 * ϕ) * cos(pi * t / τ) +
       2 * pi * R / τ * cosd(ϕ)
    ud = @. ω_0 * R / b / p_top *
       cosd(λp) *
       cosd(ϕ)^2 *
       cos(2 * pi * t / τ) *
       (-exp((p(zc) - p_0) / b / p_top) + exp((p_top - p(zc)) / b / p_top))
    uu = @. ua + ud
    uv = @. k * sind(2 * λp) * cosd(ϕ) * cos(pi * t / τ)
    ω = @. ω_0 * sind(λpf) * cosd(ϕf) * cos(2 * pi * t / τ) * sp
    uw = @. -ω / ρ(zf) / grav

    uₕ = Geometry.Covariant12Vector.(Geometry.UVVector.(uu, uv))
    w = Geometry.Covariant3Vector.(Geometry.WVector.(uw))

    ρq1 = y.ρq1
    ρq2 = y.ρq2
    ρq3 = y.ρq3
    ρq4 = y.ρq4

    dρq1 = dydt.ρq1
    dρq2 = dydt.ρq2
    dρq3 = dydt.ρq3
    dρq4 = dydt.ρq4

    If2c = Operators.InterpolateF2C()
    Ic2f = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    vdivf2c = Operators.DivergenceF2C(
        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
    )
    hdiv = Operators.Divergence()

    cw = If2c.(w)
    cuvw = Geometry.Covariant123Vector.(uₕ) .+ Geometry.Covariant123Vector.(cw)

    @. dρq1 = -hdiv(cuvw * ρq1)
    @. dρq1 -= vdivf2c(w * Ic2f(ρq1))
    @. dρq1 -= vdivf2c(Ic2f(uₕ * ρq1))

    @. dρq2 = -hdiv(cuvw * ρq2)
    @. dρq2 -= vdivf2c(w * Ic2f(ρq2))
    @. dρq2 -= vdivf2c(Ic2f(uₕ * ρq2))

    @. dρq3 = -hdiv(cuvw * ρq3)
    @. dρq3 -= vdivf2c(w * Ic2f(ρq3))
    @. dρq3 -= vdivf2c(Ic2f(uₕ * ρq3))

    @. dρq4 = -hdiv(cuvw * ρq4)
    @. dρq4 -= vdivf2c(w * Ic2f(ρq4))
    @. dρq4 -= vdivf2c(Ic2f(uₕ * ρq4))

    Spaces.weighted_dss!(dydt)

    return dydt
end

dydt = similar(y0)
rhs!(dydt, y0, (coords, face_coords), 0.0)

# run!
T = 86400 * 12
dt = 60 * 60
prob = ODEProblem(rhs!, y0, (0.0, T), (coords, face_coords))
sol = solve(
    prob,
    SSPRK33(),
    dt = dt,
    saveat = dt,
    progress = true,
    adaptive = false,
    progress_message = (dt, u, p, t) -> t,
)

q1_error =
    norm(sol.u[end].ρq1 ./ ρ.(coords.z) .- y0.ρq1 ./ ρ.(coords.z)) /
    norm(y0.ρq1 ./ ρ.(coords.z))
@test q1_error ≈ 0.0 atol = 0.9

q2_error =
    norm(sol.u[end].ρq2 ./ ρ.(coords.z) .- y0.ρq2 ./ ρ.(coords.z)) /
    norm(y0.ρq2 ./ ρ.(coords.z))
@test q2_error ≈ 0.0 atol = 0.04

q3_error =
    norm(sol.u[end].ρq3 ./ ρ.(coords.z) .- y0.ρq3 ./ ρ.(coords.z)) /
    norm(y0.ρq3 ./ ρ.(coords.z))
@test q3_error ≈ 0.0 atol = 0.6

q4_error =
    norm(sol.u[end].ρq4 ./ ρ.(coords.z) .- y0.ρq4 ./ ρ.(coords.z)) /
    norm(y0.ρq4 ./ ρ.(coords.z))
@test q4_error ≈ 0.0 atol = 0.04

# visualization artifacts
ENV["GKSwstype"] = "nul"
using ClimaCorePlots, Plots
Plots.GRBackend()
dirname = "deformation_flow"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

Plots.png(
    Plots.plot(
        sol.u[trunc(Int, end / 2)].ρq3 ./ ρ.(coords.z),
        level = 5,
        clim = (-1, 1),
    ),
    joinpath(path, "q3_6day.png"),
)

Plots.png(
    Plots.plot(sol.u[end].ρq3 ./ ρ.(coords.z), level = 5, clim = (-1, 1)),
    joinpath(path, "q3_12day.png"),
)


text:
include("MLMML.jl")

include("../lib/Newton.jl")
using ForwardDiff
using .MLMML
using PyPlot

PERIOD = 86400.0 * 360.0
J0 = 100.0 * MLMML.α * MLMML.g / MLMML.ρ / MLMML.c_p
J_func(t) = J0 / 90.0 / 86400.0 * t
E_func(t) = - J_func(t) * t/2.0
#J_func(t) = J0
#E_func(t) = - J_func(t) * t

ω  = 2π / PERIOD
h0 = 50.0
s  = 30.0 / 5000.0 * MLMML.α * MLMML.g
k  = MLMML.getTKE(fric_u=MLMML.getFricU(ua=0.0))
n  = 0.2
Δb0 = 10.0 * MLMML.α * MLMML.g *0
#sol = t -> (m * h0)^(-1.0) * (√(k^2.0 + (2.0*m*h0^2.0*(-J0)*n)/ω * (1.0 - cos(ω*t))) - k)
ts = collect(Float64, range(0.0, stop=90.0*86400.0, length=1000))
#ana_h = h0 .+ sol.(ts)

using DifferentialEquations
f(h, p, t) = h * n * J_func(t) / ( s * (h^2.0 - h0^2.0) / 2.0 + k + Δb0 * h0 + E_func(t))
prob = ODEProblem(f, h0, (ts[1], ts[end]))
num_h = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
num_dhdt = [f(num_h.u[i], 0 , num_h.t[i]) for i = 1:length(num_h.u)]

# dh/dt
#ana_dhdt = - ana_h * J0 .* sin.(ω*ts) * n ./ (m/2.0 * (ana_h.^2.0 .- h0^2.0) .+ k)
#num_dhdt = (ana_h[2:end] - ana_h[1:end-1]) / (ts[2] - ts[1])
#ana_dhdt = [f(h, 0, ts[i]) for i = 1:length(ts)]

fig, ax = plt[:subplots](2, 1, sharex=true)
ax[1][:set_title]("solution for MLD")
ax[1][:plot](num_h.t/86400.0, - (num_h.u ), "--", label="numerical")
ax[1][:legend]()

#ax[2][:plot](ts/86400.0, ana_dhdt, label="analytical")
#ax[2][:plot]((ts[2:end] + ts[1:end-1])/2.0/86400.0, num_dhdt, "--", label="numerical")
ax[2][:set_title]("dh/dt")
ax[2][:plot](num_h.t/86400.0, num_dhdt, "--", label="numerical")
ax[2][:legend]()
plt[:show]()

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
#workspace()

using Test

#Arithmetic library imports
using IntervalArithmetic
using IntervalContractors

#Preamble
setprecision(53)
setprecision(Interval, Float64)

# Using approximate checks for validaty update later?
import Base.isapprox
isapprox(x::Interval,y::Interval) = isapprox(x.lo,y.lo,atol=1E-4) && isapprox(x.hi,y.hi,atol=1E-4)

c = ∅
d = entireinterval(Float64)
a,b = exp2_rev(c, d)


text:
function FieldApsideDetector(arg0::FieldOrbit)
    return FieldApsideDetector((FieldOrbit,), arg0)
end

function FieldApsideDetector(arg0::RealFieldElement, arg1::FieldOrbit)
    return FieldApsideDetector((RealFieldElement, FieldOrbit), arg0, arg1)
end

function g(obj::FieldApsideDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end