text:
include("../utils.jl")
include("group_utils.jl")

@testset "Product group" begin
    SOn = SpecialOrthogonal(3)
    Tn = TranslationGroup(2)
    Rn = Rotations(3)
    M = ProductManifold(SOn, Tn)
    G = ProductGroup(M)
    @test_throws ErrorException ProductGroup(ProductManifold(Rotations(3), Stiefel(3, 2)))
    @test G isa ProductGroup
    @test submanifold(G, 1) === SOn
    @test submanifold(G, 2) === Tn
    @test base_manifold(G) === M
    @test sprint(show, G) == "ProductGroup($(SOn), $(Tn))"
    @test sprint(show, "text/plain", G) == "ProductGroup with 2 subgroups:\n $(SOn)\n $(Tn)"
    x = Matrix{Float64}(I, 3, 3)
    t = Vector{Float64}.([1:2, 2:3, 3:4])
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
    tuple_pts = [(exp(Rn, x, hat(Rn, x, ωi)), ti) for (ωi, ti) in zip(ω, t)]
    tuple_v = (hat(Rn, x, [1.0, 0.5, -0.5]), [-1.0, 2.0])
    eA = [x, zeros(2)]

    @testset "Product Identity" begin
        i = Identity(G)
        @test submanifold_components(G, i) == (Identity(SOn), Identity(Tn))
        @test submanifold_component(G, i, 1) == Identity(SOn)
        @test submanifold_component(G, i, 2) == Identity(Tn)
    end

    @testset "product repr" begin
        pts = [ProductRepr(tp...) for tp in tuple_pts]
        X_pts = [ProductRepr(tuple_v...)]

        @test compose(G, pts[1], Identity(G)) == pts[1]
        @test compose(G, Identity(G), pts[1]) == pts[1]
        test_group(G, pts, X_pts, X_pts; test_diff=true)
        @test isapprox(
            G,
            Identity(G),
            exp_lie(G, X_pts[1]),
            ProductRepr(exp_lie(SOn, X_pts[1].parts[1]), exp_lie(Tn, X_pts[1].parts[2])),
        )
        @test isapprox(
            G,
            Identity(G),
            log_lie(G, pts[1]),
            ProductRepr(log_lie(SOn, pts[1].parts[1]), log_lie(Tn, pts[1].parts[2])),
        )
        X = log_lie(G, pts[1])
        Z = zero_vector(G, pts[1])
        log_lie!(G, Z, pts[1])
        @test isapprox(G, pts[1], X, Z)
        p = exp_lie(G, X)
        q = identity_element(G)
        @test is_identity(G, q)
        @test isapprox(G, q, Identity(G))
        @test isapprox(G, Identity(G), q)
        exp_lie!(G, q, X)
        @test isapprox(G, p, q)
        log_lie!(G, Z, Identity(G))
        @test isapprox(G, Identity(G), Z, zero_vector(G, identity_element(G)))
        @test isapprox(
            G,
            Identity(G),
            log_lie(G, Identity(G)),
            zero_vector(G, identity_element(G)),
        )

        @test compose(G, pts[1], Identity(G)) == pts[1]
        @test compose(G, Identity(G), pts[1]) == pts[1]
        test_group(G, pts, X_pts, X_pts; test_diff=true, test_mutating=false)
        @test isapprox(
            G,
            exp_lie(G, X_pts[1]),
            ProductRepr(exp_lie(SOn, X_pts[1].parts[1]), exp_lie(Tn, X_pts[1].parts[2])),
        )
        @test isapprox(
            G,
            log_lie(G, pts[1]),
            ProductRepr(log_lie(SOn, pts[1].parts[1]), log_lie(Tn, pts[1].parts[2])),
        )
    end
    @test sprint(show, "text/plain", G) === """
    ProductGroup with 2 subgroups:
     SpecialOrthogonal(3)
     TranslationGroup(2; field = ℝ)"""
end

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
z = to_code(SUnion((Trans(Sphere(5*sqrt(2),4),space(0.0,-10.0,0.0)),
    Trans(Ball(5*sqrt(2),5,1),space(0.0,5.0,8.66)),
    Trans(Ball(5*sqrt(2),6,0.5),space(0.0,5.0,-8.66)),
    Plane(space(1.0,0.0,0.0),2),
    RepQ(Trans(Ball(0.5,1,3.5),space(2.5,2.5,2.5)),5.0),
    Trans(Ball(5,3,2.5),space(0.0,0.0,0.0)))))

scene = func(z,:space)

ext: jl


text:
module GitUtils

"""
    updating(f, url, branch; [sshkey, commit_message])

Checkout `branch` from `url` run `f(ctx)` while `cd`ing into the local
clone and then push any changes made.
"""
function updating(
    f,
    url::AbstractString,
    branch::AbstractString;
    sshkey::Union{AbstractString,Nothing} = nothing,
    commit_message::AbstractString = "",
)
    mktempdir() do tmpd
        repodir = joinpath(tmpd, "repo")
        if sshkey !== nothing
            sshdir = joinpath(tmpd, "ssh")
            GIT_SSH_COMMAND = prepare_ssh_command(sshdir, sshkey)
            git = function git_with_ssh(c)
                env = copy(ENV)
                env["GIT_SSH_COMMAND"] = GIT_SSH_COMMAND
                setenv(`git $c`, env)
            end
        else
            git = c -> `git $c`
        end

        try
            run(git(`clone --branch=$branch --depth=1 $url $repodir`))
            cd(repodir) do
                setup_git_user()
            end
        catch
            run(git(`init $repodir`))
            cd(repodir) do
                setup_git_user()
                run(git(`remote add origin $url`))
                run(git(`checkout -b $branch`))
                run(git(`commit --allow-empty --allow-empty-message --message=""`))
                _push_with_retry(git, branch)
            end
        end
        cd(repodir) do
            f((git = git,))
        end
        cd(repodir) do
            if !isempty(read(`git status --porcelain=v1`))
                @info(
                    "Committing uncommitted files.",
                    var"git status --short" = Text(read(`git status --short`, String)),
                )
                run(git(`add --all`))
                run(git(`commit --allow-empty-message --message $commit_message`))
            end
            _push_with_retry(git, branch)
        end
    end
end

function prepare_ssh_command(sshdir::AbstractString, sshkey::AbstractString)
    mkpath(sshdir)
    chmod(sshdir, 0o700)

    keyfile = joinpath(sshdir, "key")
    write(keyfile, "")
    chmod(keyfile, 0o600)
    write(keyfile, sshkey)

    return "ssh -i $keyfile"
end

function setup_git_user()
    GITHUB_ACTOR = get(ENV, "GITHUB_ACTOR", nothing)
    GITHUB_ACTOR === nothing && return
    # WARNING: This function is run via tests. Do NOT use `--global`.
    run(`git config user.email "$GITHUB_ACTOR@users.noreply.github.com"`)
    run(`git config user.name $GITHUB_ACTOR`)
end

printable(cmd::Cmd) = Cmd(cmd.exec::Vector{String})

function saferun(cmd::Cmd; log = true)
    @info "Executing $(printable(cmd))..."
    try
        return run(cmd)
    catch err
        log && @info "Failed to run: $(printable(cmd))"
        @debug("Failed to run: $(printable(cmd))", exception = (err, catch_backtrace()),)
        return nothing
    end
end

function _push_with_retry(git, branch; timeout = 5 * 60, tries = 10)
    t0 = time_ns()

    push_cmd = git(`push origin "HEAD:refs/heads/$branch"`)
    for i in 1:tries
        if  saferun(push_cmd; log = false) !== nothing
            @info "Successfully pushed branch `$branch` to remote."
            return
        end

        if (time_ns() - t0) / 1e9 > timeout
            @info "Timeout ($timeout seconds) reached."
            break
        end
        @info "$i-th $(printable(push_cmd)) failed; retrying...."

        saferun(git(`fetch origin $branch`)) == nothing && break
        saferun(git(`merge --no-edit "origin/$branch"`)) == nothing && break
    end
    error("Failed to: $(printable(push_cmd))")
end

end  # module


text:
module FactorAnalysis

using Distributions

export CFADistribution, rand, loglikelihood, dloglikelihood, fit_mle, fit_map

include("sigma.jl")
include("theta.jl")
include("em.jl")
include("utils.jl")
#include("pgradient.jl")
#include("em.jl")

end # module


text:
ToeplitzOperator(f::Fun{Laurent{DD,RR}}) where {DD,RR}=ToeplitzOperator(f.coefficients[2:2:end],
                                                    f.coefficients[1:2:end])

##Taylor

Multiplication(f::Fun{Taylor{DD,RR}},sp::Taylor{DD,RR}) where {DD,RR} =
    MultiplicationWrapper(f,SpaceOperator(ToeplitzOperator(f.coefficients[2:end],
                                                           [f.coefficients[1]]),
                                          sp,sp))

## Evaluation

getindex(T::ConcreteEvaluation{Taylor{DD,RR},Complex{Float64},Int,Complex{Float64}},cols::AbstractRange) where {DD,RR} =
    mappoint(domain(T),Circle(),T.x).^(cols-1)

## Multiplication

Multiplication(f::Fun{Laurent{DD,RR}},sp::Laurent{DD,RR}) where {DD,RR} = ConcreteMultiplication(cfstype(f),f,sp)

function laurent_getindex(negative::AbstractVector{T},nonnegative::AbstractVector{T},k::Integer,j::Integer) where T
    # switch to double-infinite indices
    k=iseven(k) ? -k÷2 : (k-1)÷2
    j=iseven(j) ? -j÷2 : (j-1)÷2

    if 0<k-j≤length(negative)
        negative[k-j]
    elseif 0≤j-k≤length(nonnegative)-1
        nonnegative[j-k+1]
    else
        zero(T)
    end
end

rangespace(T::ConcreteMultiplication{Laurent{DD,RR},Laurent{DD,RR}}) where {DD,RR} = domainspace(T)
function getindex(T::ConcreteMultiplication{Laurent{DD,RR},Laurent{DD,RR}},k::Integer,j::Integer) where {DD,RR}
    isempty(T.f.coefficients) && return zero(eltype(T))
    laurent_getindex(T.f.coefficients[3:2:end],T.f.coefficients[[1;2:2:end]],k,j)
end

function bandwidths(T::ConcreteMultiplication{Laurent{DD,RR},Laurent{DD,RR}}) where {DD,RR}
    bbi = blockbandwidths(T)
    (2bbi[1]+1,2bbi[2]+1)
end

function blockbandwidths(T::ConcreteMultiplication{Laurent{DD,RR},Laurent{DD,RR}}) where {DD,RR}
    m = ncoefficients(T.f)÷2
    (m,m)
end

Multiplication(f::Fun{Fourier{DD,R1}},sp::Laurent{DD,R2}) where {DD,R1,R2} = Multiplication(Fun(f,sp),sp)
Multiplication(f::Fun{Laurent{DD,R1}},sp::Fourier{DD,R2}) where {DD,R1,R2} = Multiplication(Fun(f,sp),sp)

# override SumSpace default
coefficienttimes(f::Fun{Laurent{DD,RR}},g::Fun{Laurent{DD,RR}}) where {DD,RR} = Multiplication(f,space(g))*g

## Derivative

# override map definition
Derivative(S::Hardy{s,DD},k::Integer) where {s,DD<:Circle} = ConcreteDerivative(S,k)
Derivative(S::Hardy{s,DD},k::Integer) where {s,DD<:PeriodicSegment} = ConcreteDerivative(S,k)
Derivative(S::Laurent{DD,RR},k::Integer) where {DD<:Circle,RR} =
    DerivativeWrapper(InterlaceOperator(Diagonal([map(s->Derivative(s,k),S.spaces)...]),SumSpace),k)

bandwidths(D::ConcreteDerivative{Hardy{s,DD,RR}}) where {s,DD<:PeriodicSegment,RR}=(0,0)
bandwidths(D::ConcreteDerivative{Hardy{s,DD,RR}}) where {s,DD<:Circle,RR}=s ? (0,D.order) : (D.order,0)

rangespace(D::ConcreteDerivative{S}) where {S<:Hardy}=D.space

function taylor_derivative_getindex(d::PeriodicSegment,m,k::Integer,j::Integer)
    C=2/complexlength(d)*π*im
    k==j ? (C*(k-1))^m : zero(C)
end

function hardyfalse_derivative_getindex(d::PeriodicSegment,m,k::Integer,j::Integer)
    C=2/complexlength(d)*π*im
    k==j ? (-C*k)^m : zero(C)
end

function taylor_derivative_getindex(d::Circle,m,k::Integer,j::Integer)
    C=d.radius^(-m)

    if j==k+m
        D=k
        for j=k+1:k+m-1
          D*=j
        end
        C*D
    else
        zero(C)
    end
end

function hardyfalse_derivative_getindex(d::Circle,m,k::Integer,j::Integer)
    C=(-d.radius)^(-m)

    if j == k-m
        D=k-m
        for j=k-m+1:k-1
          D*=j
        end
        C*D
    else
        zero(C)
    end
end

getindex(D::ConcreteDerivative{Taylor{DD,RR},OT,T},k::Integer,j::Integer) where {DD,RR,OT,T}  =
    convert(T,taylor_derivative_getindex(domain(D),D.order,k,j))
getindex(D::ConcreteDerivative{Hardy{false,DD,RR},OT,T},k::Integer,j::Integer) where {DD,RR,OT,T} =
    convert(T,hardyfalse_derivative_getindex(domain(D),D.order,k,j))

## Integral

# Integral{D<:Circle}(S::Taylor{D},m)=Integral{Taylor,typeof(m),Complex{Float64}}(S,m)
#
# function Integral{D<:PeriodicSegment}(S::Hardy{false,D},m)
#     Integral{Hardy{false,D},typeof(m),Complex{Float64}}(S,m)
# end

Integral(S::Hardy{s,DD,RR},k::Integer) where {s,DD<:Circle,RR} = ConcreteIntegral(S,k)

bandwidths(D::ConcreteIntegral{Taylor{DD,RR}}) where {DD<:Circle,RR} = (D.order,0)
rangespace(Q::ConcreteIntegral{Hardy{s,DD,RR}}) where {s,DD<:Circle,RR} = Q.space

function getindex(D::ConcreteIntegral{Taylor{DD,RR}},k::Integer,j::Integer) where {DD<:Circle,RR}
    d=domain(D)
    m=D.order
    T=eltype(D)
    C=convert(T,d.radius^m)

    if j==k-m
        D=k-m
        for j=k-m+1:k-1
          D*=j
        end
        C/D
    else
        zero(T)
    end
end

function Integral(S::SubSpace{<:Hardy{false,<:Circle}, <:AbstractInfUnitRange{Int}},k::Integer)
    if first(S.indexes) == k+1
        ConcreteIntegral(S,k)
    else
        @assert first(S.index) > k+1
        S2=S.space|(k+1:∞)
        IntegralWrapper(ConcreteIntegral(S2,k)*Converion(S,S2),k)
    end
end

bandwidths(D::ConcreteIntegral{<:SubSpace{<:Hardy{false,<:Circle}, <:AbstractInfUnitRange{Int}}}) =
    (0,0)

rangespace(D::ConcreteIntegral{<:SubSpace{<:Hardy{false,<:Circle}, <:AbstractInfUnitRange{Int}}}) =
    D.space.space

function getindex(D::ConcreteIntegral{<:SubSpace{<:Hardy{false,<:Circle}, <:AbstractInfUnitRange{Int}},OT,TT},
                 k::Integer,j::Integer) where {OT,TT}
    d=domain(D)
    m=D.order

    C=convert(TT,(-d.radius)^m)

    if k==j
        D=k
        for j=k+1:k+m-1
          D*=j
        end
        C/D
    else
        zero(TT)
    end
end

bandwidths(D::ConcreteIntegral{Hardy{false,DD,RR}}) where {DD<:PeriodicSegment,RR} = (0,0)
rangespace(D::ConcreteIntegral{Taylor{DD,RR}}) where {DD<:PeriodicSegment,RR} = D.space

function getindex(D::ConcreteIntegral{Hardy{false,DD,RR}},k::Integer,j::Integer) where {DD<:PeriodicSegment,RR}
    d=domain(D)
    m=D.order
    T=eltype(D)
    C=2/complexlength(d)*π*im
    if k==j
        convert(T,(-C*k)^(-m))
    else
        zero(T)
    end
end

bandwidths(D::ConcreteIntegral{<:SubSpace{Taylor{DD,RR},<:AbstractInfUnitRange{Int},DD,RR}}) where {RR,DD<:PeriodicSegment} =
    (0,0)
rangespace(D::ConcreteIntegral{<:SubSpace{Taylor{DD,RR},<:AbstractInfUnitRange{Int},DD,RR}}) where {RR,DD<:PeriodicSegment} =
    D.space

function getindex(D::ConcreteIntegral{<:SubSpace{Taylor{DD,RR},<:AbstractInfUnitRange{Int},DD,RR}},
                 k::Integer,j::Integer) where {RR,DD<:PeriodicSegment}
    d=domain(D)
    m=D.order
    TT=eltype(D)
    C=2/complexlength(d)*π*im
    if k==j
        convert(TT,(C*(k+n-1))^(-m))
    else
        zero(TT)
    end
end

## Definite integral

for SP in (:Taylor,:(Hardy{false}),:Laurent)
    @eval begin
        DefiniteIntegral(S::$SP{D,R}) where {D,R} =
            ConcreteDefiniteIntegral{typeof(S),prectype(S)}(S)
        DefiniteLineIntegral(S::$SP{D,R}) where {D,R} =
            ConcreteDefiniteLineIntegral{typeof(S),real(prectype(S))}(S)
    end
end

getindex(Σ::ConcreteDefiniteIntegral{Taylor{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    k == 1 ? convert(T,complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Hardy{false,D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Laurent{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    k == 1 ? convert(T,complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Taylor{D,R},T},k::Integer) where {T,D<:Circle,R} =
    zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Hardy{false,D,R},T},k::Integer) where {T,D<:Circle,R} =
    k == 1 ? convert(T,complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Laurent{D,R},T},k::Integer) where {T,D<:Circle,R} =
    k == 2 ? convert(T,complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{Taylor{D,R},T},k::Integer) where {T,D,R} =
    k == 1 ? convert(T,arclength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{Hardy{false,D,R},T},k::Integer) where {T,D,R} =
    zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{Laurent{D,R},T},k::Integer) where {T,D,R} =
    k == 1 ? convert(T,arclength(domain(Σ))) : zero(T)

bandwidths(Σ::ConcreteDefiniteIntegral{Taylor{D,R}}) where {D<:PeriodicSegment,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Hardy{false,D,R}}) where {D<:PeriodicSegment,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Laurent{D,R}}) where {D<:PeriodicSegment,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Taylor{D,R}}) where {D<:Circle,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Hardy{false,D,R}}) where {D<:Circle,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Laurent{D,R}}) where {D<:Circle,R} = 0,1
bandwidths(Σ::ConcreteDefiniteLineIntegral{Taylor{D,R}}) where {D,R} = 0,0
bandwidths(Σ::ConcreteDefiniteLineIntegral{Hardy{false,D,R}}) where {D,R} = 0,0
bandwidths(Σ::ConcreteDefiniteLineIntegral{Laurent{D,R}}) where {D,R} = 0,0

## reverse orientation

conversion_type(A::Laurent{DD,RR},B::Laurent{DD,RR}) where {DD<:Circle,RR}=domain(A).orientation ? A : B
function Conversion(A::Laurent{DD,RR},B::Laurent{DD,RR}) where {DD,RR}
    @assert domain(A) == reverseorientation(domain(B))
    ConversionWrapper(SpaceOperator(
        InterlaceOperator(Diagonal([Matrix(I,1,1),PermutationOperator([2,1])]))
    ,A,B))
end


text:
import Test

import Thermodynamics
import CloudMicrophysics
import CLIMAParameters

const TT = Test
const TD = Thermodynamics
const CM0 = CloudMicrophysics.Microphysics_0M
const CM1 = CloudMicrophysics.Microphysics_1M
const CP = CLIMAParameters
const CP_planet = CLIMAParameters.Planet
const CP_micro0 = CLIMAParameters.Atmos.Microphysics_0M
const CP_micro1 = CLIMAParameters.Atmos.Microphysics

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const prs = EarthParameterSet()

const liquid = CM1.LiquidType()
const ice = CM1.IceType()
const rain = CM1.RainType()
const snow = CM1.SnowType()

TT.@testset "τ_relax" begin

    TT.@test CM1.τ_relax(prs, liquid) ≈ 10
    TT.@test CM1.τ_relax(prs, ice) ≈ 10

end

TT.@testset "0M_microphysics" begin

    _τ_precip = CP_micro0.τ_precip(prs)
    _qc_0 = CP_micro0.qc_0(prs)
    _S_0 = CP_micro0.S_0(prs)

    q_vap_sat = 10e-3
    qc = 3e-3
    q_tot = 13e-3
    frac = [0.0, 0.5, 1.0]

    # no rain if no cloud
    q = TD.PhasePartition(q_tot, 0.0, 0.0)
    TT.@test CM0.remove_precipitation(prs, q) ≈ 0
    TT.@test CM0.remove_precipitation(prs, q, q_vap_sat) ≈ 0

    # rain based on qc threshold
    for lf in frac
        q_liq = qc * lf
        q_ice = (1 - lf) * qc

        q = TD.PhasePartition(q_tot, q_liq, q_ice)

        TT.@test CM0.remove_precipitation(prs, q) ≈
                 -max(0, q_liq + q_ice - _qc_0) / _τ_precip
    end

    # rain based on supersaturation threshold
    for lf in frac
        q_liq = qc * lf
        q_ice = (1 - lf) * qc

        q = TD.PhasePartition(q_tot, q_liq, q_ice)

        TT.@test CM0.remove_precipitation(prs, q, q_vap_sat) ≈
                 -max(0, q_liq + q_ice - _S_0 * q_vap_sat) / _τ_precip
    end
end

TT.@testset "RainFallSpeed" begin
    # eq. 5d in [Grabowski1996](@cite)
    function terminal_velocity_empir(
        q_rai::FT,
        q_tot::FT,
        ρ::FT,
        ρ_air_ground::FT,
    ) where {FT <: Real}
        rr = q_rai / (1 - q_tot)
        vel = FT(14.34) * ρ_air_ground^FT(0.5) * ρ^-FT(0.3654) * rr^FT(0.1346)
        return vel
    end

    # some example values
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    ρ_air, q_tot, ρ_air_ground = 1.2, 20 * 1e-3, 1.22

    for q_rai in q_rain_range
        TT.@test CM1.terminal_velocity(prs, rain, ρ_air, q_rai) ≈
                 terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground) atol =
            0.2 * terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground)

    end
end

TT.@testset "CloudLiquidCondEvap" begin

    q_liq_sat = 5e-3
    frac = [0.0, 0.5, 1.0, 1.5]

    _τ_cond_evap = CM1.τ_relax(prs, liquid)

    for fr in frac
        q_liq = q_liq_sat * fr

        TT.@test CM1.conv_q_vap_to_q_liq_ice(
            prs,
            liquid,
            TD.PhasePartition(0.0, q_liq_sat, 0.0),
            TD.PhasePartition(0.0, q_liq, 0.0),
        ) ≈ (1 - fr) * q_liq_sat / _τ_cond_evap
    end
end

TT.@testset "CloudIceCondEvap" begin

    q_ice_sat = 2e-3
    frac = [0.0, 0.5, 1.0, 1.5]

    _τ_cond_evap = CM1.τ_relax(prs, ice)

    for fr in frac
        q_ice = q_ice_sat * fr

        TT.@test CM1.conv_q_vap_to_q_liq_ice(
            prs,
            ice,
            TD.PhasePartition(0.0, 0.0, q_ice_sat),
            TD.PhasePartition(0.0, 0.0, q_ice),
        ) ≈ (1 - fr) * q_ice_sat / _τ_cond_evap
    end
end

TT.@testset "RainAutoconversion" begin

    _q_liq_threshold = CP_micro1.q_liq_threshold(prs)
    _τ_acnv_rai = CP_micro1.τ_acnv_rai(prs)

    q_liq_small = 0.5 * _q_liq_threshold
    TT.@test CM1.conv_q_liq_to_q_rai(prs, q_liq_small) == 0.0

    q_liq_big = 1.5 * _q_liq_threshold
    TT.@test CM1.conv_q_liq_to_q_rai(prs, q_liq_big) ==
             0.5 * _q_liq_threshold / _τ_acnv_rai
end

TT.@testset "SnowAutoconversionNoSupersat" begin

    _q_ice_threshold = CP_micro1.q_ice_threshold(prs)
    _τ_acnv_sno = CP_micro1.τ_acnv_sno(prs)

    q_ice_small = 0.5 * _q_ice_threshold
    TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(prs, q_ice_small) == 0.0

    q_ice_big = 1.5 * _q_ice_threshold
    TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(prs, q_ice_big) ≈
             0.5 * _q_ice_threshold / _τ_acnv_sno
end

TT.@testset "SnowAutoconversion" begin

    ρ = 1.0

    # above freezing temperatures -> no snow
    q = TD.PhasePartition(15e-3, 2e-3, 1e-3)
    T = 273.15 + 30
    TT.@test CM1.conv_q_ice_to_q_sno(prs, q, ρ, T) == 0.0

    # no ice -> no snow
    q = TD.PhasePartition(15e-3, 2e-3, 0.0)
    T = 273.15 - 30
    TT.@test CM1.conv_q_ice_to_q_sno(prs, q, ρ, T) == 0.0

    # no supersaturation -> no snow
    T = 273.15 - 5
    q_sat_ice = TD.q_vap_saturation_generic(prs, T, ρ, TD.Ice())
    q = TD.PhasePartition(q_sat_ice, 2e-3, 3e-3)
    TT.@test CM1.conv_q_ice_to_q_sno(prs, q, ρ, T) == 0.0

    # TODO - coudnt find a plot of what it should be from the original paper
    # just chacking if the number stays the same
    T = 273.15 - 10
    q_vap = 1.02 * TD.q_vap_saturation_generic(prs, T, ρ, TD.Ice())
    q_liq = 0.0
    q_ice = 0.03 * q_vap
    q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
    TT.@test CM1.conv_q_ice_to_q_sno(prs, q, ρ, T) ≈ 1.8512022335645584e-9

end

TT.@testset "RainLiquidAccretion" begin

    # eq. 5b in [Grabowski1996](@cite)
    function accretion_empir(q_rai::FT, q_liq::FT, q_tot::FT) where {FT <: Real}
        rr = q_rai / (FT(1) - q_tot)
        rl = q_liq / (FT(1) - q_tot)
        return FT(2.2) * rl * rr^FT(7 / 8)
    end

    # some example values
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    ρ_air, q_liq, q_tot = 1.2, 5e-4, 20e-3

    for q_rai in q_rain_range
        TT.@test CM1.accretion(prs, liquid, rain, q_liq, q_rai, ρ_air) ≈
                 accretion_empir(q_rai, q_liq, q_tot) atol =
            (0.1 * accretion_empir(q_rai, q_liq, q_tot))
    end
end

TT.@testset "Accretion" begin
    # TODO - coudnt find a plot of what it should be from the original paper
    # just chacking if the number stays the same

    # some example values
    ρ = 1.2
    q_tot = 20e-3
    q_ice = 5e-4
    q_sno = 5e-4
    q_liq = 5e-4
    q_rai = 5e-4

    TT.@test CM1.accretion(prs, liquid, rain, q_liq, q_rai, ρ) ≈
             1.4150106417043544e-6
    TT.@test CM1.accretion(prs, ice, snow, q_ice, q_sno, ρ) ≈
             2.453070979562392e-7
    TT.@test CM1.accretion(prs, liquid, snow, q_liq, q_sno, ρ) ≈
             2.453070979562392e-7
    TT.@test CM1.accretion(prs, ice, rain, q_ice, q_rai, ρ) ≈
             1.768763302130443e-6

    TT.@test CM1.accretion_rain_sink(prs, q_ice, q_rai, ρ) ≈
             3.085229094251214e-5

    TT.@test CM1.accretion_snow_rain(prs, snow, rain, q_sno, q_rai, ρ) ≈
             2.1705865794293408e-4
    TT.@test CM1.accretion_snow_rain(prs, rain, snow, q_rai, q_sno, ρ) ≈
             6.0118801860768854e-5
end

TT.@testset "RainEvaporation" begin

    # eq. 5c in [Grabowski1996](@cite)
    function rain_evap_empir(
        prs::CP.AbstractParameterSet,
        q_rai::FT,
        q::TD.PhasePartition,
        T::FT,
        p::FT,
        ρ::FT,
    ) where {FT <: Real}

        q_sat = TD.q_vap_saturation_generic(prs, T, ρ, TD.Liquid())
        q_vap = q.tot - q.liq
        rr = q_rai / (1 - q.tot)
        rv_sat = q_sat / (1 - q.tot)
        S = q_vap / q_sat - 1

        ag, bg = FT(5.4 * 1e2), FT(2.55 * 1e5)
        G = FT(1) / (ag + bg / p / rv_sat) / ρ

        av, bv = FT(1.6), FT(124.9)
        F =
            av * (ρ / FT(1e3))^FT(0.525) * rr^FT(0.525) +
            bv * (ρ / FT(1e3))^FT(0.7296) * rr^FT(0.7296)

        return 1 / (1 - q.tot) * S * F * G
    end

    # example values
    T, p = 273.15 + 15, 90000.0
    ϵ = 1.0 / CP_planet.molmass_ratio(prs)
    p_sat = TD.saturation_vapor_pressure(prs, T, TD.Liquid())
    q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    q_tot = 15e-3
    q_vap = 0.15 * q_sat
    q_ice = 0.0
    q_liq = q_tot - q_vap - q_ice
    q = TD.PhasePartition(q_tot, q_liq, q_ice)
    R = TD.gas_constant_air(prs, q)
    ρ = p / R / T

    for q_rai in q_rain_range
        TT.@test CM1.evaporation_sublimation(prs, rain, q, q_rai, ρ, T) ≈
                 rain_evap_empir(prs, q_rai, q, T, p, ρ) atol =
            -0.5 * rain_evap_empir(prs, q_rai, q, T, p, ρ)
    end

    # no condensational growth for rain
    T, p = 273.15 + 15, 90000.0
    ϵ = 1.0 / CP_planet.molmass_ratio(prs)
    p_sat = TD.saturation_vapor_pressure(prs, T, TD.Liquid())
    q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))
    q_rai = 1e-4
    q_tot = 15e-3
    q_vap = 1.15 * q_sat
    q_ice = 0.0
    q_liq = q_tot - q_vap - q_ice
    q = TD.PhasePartition(q_tot, q_liq, q_ice)
    R = TD.gas_constant_air(prs, q)
    ρ = p / R / T

    TT.@test CM1.evaporation_sublimation(prs, rain, q, q_rai, ρ, T) ≈ 0.0

end

TT.@testset "SnowSublimation" begin
    # TODO - coudnt find a plot of what it should be from the original paper
    # just chacking if the number stays the same

    cnt = 0
    ref_val = [
        -1.9756907119482267e-7,
        1.9751292385808357e-7,
        -1.6641552112891826e-7,
        1.663814937710236e-7,
    ]
    # some example values
    for T in [273.15 + 2, 273.15 - 2]
        p = 90000.0
        ϵ = 1.0 / CP_planet.molmass_ratio(prs)
        p_sat = TD.saturation_vapor_pressure(prs, T, TD.Ice())
        q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))

        for eps in [0.95, 1.05]
            cnt += 1

            q_tot = eps * q_sat
            q_ice = 0.0
            q_liq = 0.0
            q = TD.PhasePartition(q_tot, q_liq, q_ice)

            q_sno = 1e-4

            R = TD.gas_constant_air(prs, q)
            ρ = p / R / T

            TT.@test CM1.evaporation_sublimation(prs, snow, q, q_sno, ρ, T) ≈
                     ref_val[cnt]

        end
    end
end

TT.@testset "SnowMelt" begin

    # TODO - find a good reference to compare with
    T = 273.15 + 2
    ρ = 1.2
    q_sno = 1e-4
    TT.@test CM1.snow_melt(prs, q_sno, ρ, T) ≈ 9.518235437405256e-6

    # no snow -> no snow melt
    T = 273.15 + 2
    ρ = 1.2
    q_sno = 0.0
    TT.@test CM1.snow_melt(prs, q_sno, ρ, T) ≈ 0

    # T < T_freeze -> no snow melt
    T = 273.15 - 2
    ρ = 1.2
    q_sno = 1e-4
    TT.@test CM1.snow_melt(prs, q_sno, ρ, T) ≈ 0

end


text:
module Fatou
using SyntaxTree,Reduce,LaTeXStrings,Requires,Base.Threads

#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

export fatou, juliafill, mandelbrot, newton, basin, plot

"""
    Fatou.Define(E::Any;                  # primary map, (z, c) -> F
      Q::Expr     = :(abs2(z)),           # escape criterion, (z, c) -> Q
      C::Expr     = :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      newt::Bool  = false,                # toggle Newton mode
      m::Number   = 0,                    # Newton multiplicity factor
      O::String   = F,                    # original Newton map
      mandel::Bool= false,                # toggle Mandelbrot mode
      seed::Number= 0.0+0.0im,            # Mandelbrot seed value
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` the metadata for a `Fatou.FilledSet`.
"""
struct Define{FT<:Function,QT<:Function,CT<:Function,M,N,P,D}
    E::Any # input expression
    F::FT # primary map
    Q::QT # escape criterion
    C::CT # complex fixed point coloring
    ∂::Array{Float64,1} # bounds
    n::UInt16 # number of grid points
    N::UInt8 # number of iterations
    ϵ::Float64 # epsilon Limit criterion
    iter::Bool # toggle iteration mode
    p::Float64 # iteration color exponent
    newt::Bool # toggle Newton mode
    m::Number # newton multiplicity factor
    mandel::Bool # toggle Mandelbrot mode
    seed::Number # Mandelbrot seed value
    x0 # orbit starting point
    orbit::Int # orbit cobweb depth
    depth::Int # depth of function composition
    cmap::String # imshow color map
    plane::Bool # convert input disk to half-plane
    disk::Bool # convert output half-plane to disk
    function Define(E;
            Q=:(abs2(z)),
            C=:((angle(z)/(2π))*n^p),
            ∂=π/2,
            n::Integer=176,
            N::Integer=35,
            ϵ::Number=4,
            iter::Bool=false,
            p::Number=0,
            newt::Bool=false,
            m::Number=0,
            mandel::Bool=false,
            seed::Number=0.0+0.0im,
            x0=nothing,
            orbit::Int=0,
            depth::Int=1,
            cmap::String="",
            plane::Bool=false,
            disk::Bool=false)
        !(typeof(∂) <: Array) && (∂ = [-float(∂),∂,-∂,∂])
        length(∂) == 2 && (∂ = [∂[1],∂[2],∂[1],∂[2]])
        !newt ? (f = genfun(E,[:z,:c]); q = genfun(Q,[:z,:c])) :
        (f = genfun(newton_raphson(E,m),[:z,:c]); q = genfun(Expr(:call,:abs,E),[:z,:c]))
        c = genfun(C,[:z,:n,:p])
        e = typeof(E) == String ? parse(E) : E
        return new{typeof(f),typeof(q),typeof(c),mandel,newt,plane,disk}(e,f,q,c,convert(Array{Float64,1},∂),UInt16(n),UInt8(N),float(ϵ),iter,float(p),newt,m,mandel,seed,x0,orbit,depth,cmap,plane,disk)
    end
end

"""
    Fatou.FilledSet(::Fatou.Define)

Compute the `Fatou.FilledSet` set using `Fatou.Define`.
"""
struct FilledSet{FT,QT,CT,M,N,P,D}
    meta::Define{FT,QT,CT,M,N,P,D}
    set::Matrix{Complex{Float64}}
    iter::Matrix{UInt8}
    mix::Matrix{Float64}
    function FilledSet{FT,QT,CT,M,N,P,D}(K::Define{FT,QT,CT,M,N,P,D}) where {FT,QT,CT,M,N,P,D}
        (i,s) = Compute(K)
        return new{FT,QT,CT,M,N,P,D}(K,s,i,broadcast(K.C,s,broadcast(float,i./K.N),K.p))
    end
end

  """
      fatou(::Fatou.Define)

Compute the `Fatou.FilledSet` set using `Fatou.Define`.

# Examples
```Julia
julia> fatou(K)
```
"""
fatou(K::Define{FT,QT,CT,M,N,P,D}) where {FT,QT,CT,M,N,P,D} = FilledSet{FT,QT,CT,M,N,P,D}(K)

"""
    juliafill(::Expr;                     # primary map, (z, c) -> F
      Q::Expr     = :(abs2(z)),           # escape criterion, (z, c) -> Q
      C::Expr     = :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` filled Julia basin in `Fatou`

# Exmaples
```Julia
julia> juliafill("z^2-0.06+0.67im",∂=[-1.5,1.5,-1,1],N=80,n=1501,cmap="RdGy")
```
"""
function juliafill(E;
        Q=:(abs2(z)),
        C=:((angle(z)/(2π))*n^p),
        ∂=π/2,
        n::Integer=176,
        N::Integer=35,
        ϵ::Number=4,
        iter::Bool=false,
        p::Number=0,
        newt::Bool=false,
        m::Number=0,
        x0=nothing,
        orbit::Int=0,
        depth::Int=1,
        cmap::String="",
        plane::Bool=false,
        disk::Bool=false)
    return Define(E,Q=Q,C=C,∂=∂,n=n,N=N,ϵ=ϵ,iter=iter,p=p,newt=newt,m=m,x0=x0,orbit=orbit,depth=depth,cmap=cmap,plane=plane,disk=disk)
end

"""
    mandelbrot(::Expr;                    # primary map, (z, c) -> F
      Q::Expr     = :(abs2(z)),           # escape criterion, (z, c) -> Q
      C::Expr     = :(exp(-abs(z))*n^p),  # coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      m::Number   = 0,                    # Newton multiplicity factor
      seed::Number= 0.0+0.0im,            # Mandelbrot seed value
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` Mandelbrot basin in `Fatou`

# Examples
```Julia
mandelbrot(:(z^2+c),n=800,N=20,∂=[-1.91,0.51,-1.21,1.21],cmap="nipy_spectral")
```
"""
function mandelbrot(E;
        Q=:(abs2(z)),
        ∂=π/2,
        C=:(exp(-abs(z))*n^p),
        n::Integer=176,
        N::Integer=35,
        ϵ::Number=4,
        iter::Bool=false,
        p::Number=0,
        newt::Bool=false,
        m::Number=0,
        seed::Number=0.0+0.0im,
        x0=nothing,
        orbit::Int=0,
        depth::Int=1,
        cmap::String="",
        plane::Bool=false,
        disk::Bool=false)
    m ≠ 0 && (newt = true)
    return Define(E,Q=Q,C=C,∂=∂,n=n,N=N,ϵ=ϵ,iter=iter,p=p,newt=newt,m=m,mandel=true,seed=seed,x0=x0,orbit=orbit,depth=depth,cmap=cmap,plane=plane,disk=disk)
end

"""
    newton(::Expr;                        # primary map, (z, c) -> F
      C::Expr     = :((angle(z)/(2π))*n^p)# coloring, (z, n=iter., p=exp.) -> C
      ∂    = π/2, # Array{Float64,1}      # Bounds, [x(a),x(b),y(a),y(b)]
      n::Integer  = 176,                  # horizontal grid points
      N::Integer  = 35,                   # max. iterations
      ϵ::Number   = 4,                    # basin ϵ-Limit criterion
      iter::Bool  = false,                # toggle iteration mode
      p::Number   = 0,                    # iteration color exponent
      m::Number   = 0,                    # Newton multiplicity factor
      mandel::Bool= false,                # toggle Mandelbrot mode
      seed::Number= 0.0+0.0im,            # Mandelbrot seed value
      x0          = nothing,              # orbit starting point
      orbit::Int  = 0,                    # orbit cobweb depth
      depth::Int  = 1,                    # depth of function composition
      cmap::String= "",                   # imshow color map
      plane::Bool = false,                # convert input disk to half-plane
      disk::Bool  = false)                # convert output half-plane to disk

`Define` Newton basin in `Fatou`

# Examples
```Julia
julia> newton("z^3-1",n=800,cmap="brg")
```
"""
function newton(E;
        C=:((angle(z)/(2π))*n^p),
        ∂=π/2,
        n::Integer=176,
        N::Integer=35,
        ϵ::Number=0.01,
        iter::Bool=false,
        p::Number=0,
        m::Number=1,
        mandel::Bool=false,
        seed::Number=0.0+0.0im,
        x0=nothing,
        orbit::Int=0,
        depth::Int=1,
        cmap::String="",
        plane::Bool=false,
        disk::Bool=false)
    return Define(E,C=C,∂=∂,n=n,N=N,ϵ=ϵ,iter=iter,p=p,newt=true,m=m,mandel=mandel,seed=seed,x0=x0,orbit=orbit,depth=depth,cmap=cmap,plane=plane,disk=disk)
end

# load additional functionality
include("internals.jl"); include("orbitplot.jl")

"""
    basin(::Fatou.Define, ::Integer)

Output the `j`-th basin of `Fatou.Define` as LaTeX.
Each subsequent iteration of the Newton-Raphson method will yield a more complicated set.

# Examples
```Julia
julia> basin(newton("z^3-1"),2)
L"\$\\displaystyle D_2(\\epsilon) = \\left\\{z\\in\\mathbb{C}:\\left|\\,z - \\frac{\\left(z - \\frac{z^{3} - 1}{3 z^{2}}\\right)^{3} - 1}{3 \\left(z - \\frac{z^{3} - 1}{3 z^{2}}\\right)^{2}} - \\frac{z^{3} - 1}{3 z^{2}} - r_i\\,\\right|<\\epsilon,\\,\\forall r_i(\\,f(r_i)=0 )\\right\\}\$"
```
"""
basin(K::Define,j) = K.newt ? nrset(K.E,K.m,j) : jset(K.E,j)

plane(z::Complex) = (2z.re/(z.re^2+(1-z.im)^2))+im*(1-z.re^2-z.im^2)/(z.re^2+(1-z.im)^2)
disk(z::Complex) = (2z.re/(z.re^2+(1+z.im)^2))+im*(z.re^2+z.im^2-1)/(z.re^2+(1+z.im)^2)

# define function for computing orbit of a z0 input
function orbit(K::Define{FT,QT,CT,M,N,P,D},z0::Complex{Float64}) where {FT,QT,CT,M,N,P,D}
    M ? (z = K.seed) : (z = P ? plane(z0) : z0)
    zn = 0x00
    while (N ? (K.Q(z,z0)::Float64>K.ϵ)::Bool : (K.Q(z,z0)::Float64<K.ϵ))::Bool && K.N>zn
        z = K.F(z,z0)::Complex{Float64}
        zn+=0x01
    end; #end
    # return the normalized argument of z or iteration count
    return (zn::UInt8,(D ? disk(z) : z)::Complex{Float64})
end

"""
    Compute(::Fatou.Define)::Union{Matrix{UInt8},Matrix{Float64}}

`Compute` the `Array` for `Fatou.FilledSet` as specefied by `Fatou.Define`.
"""
function Compute(K::Define{FT,QT,CT,M,N,D})::Tuple{Matrix{UInt8},Matrix{Complex{Float64}}} where {FT,QT,CT,M,N,D}
    # generate coordinate grid
    Kyn = round(UInt16,(K.∂[4]-K.∂[3])/(K.∂[2]-K.∂[1])*K.n)
    x = range(K.∂[1]+0.0001,stop=K.∂[2],length=K.n)
    y = range(K.∂[4],stop=K.∂[3],length=Kyn)
    Z = x' .+ im*y # apply Newton-Orbit function element-wise to coordinate grid
    (matU,matF) = (Array{UInt8,2}(undef,Kyn,K.n),Array{Complex{Float64},2}(undef,Kyn,K.n))
    @time @threads for j = 1:length(y); for k = 1:length(x);
        (matU[j,k],matF[j,k]) = orbit(K,Z[j,k])::Tuple{UInt8,Complex{Float64}}
    end; end
    return (matU,matF)
end

function __init__()
    println("Fatou detected $(Threads.nthreads()) julia threads.")
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("pyplot.jl")
    @require ImageInTerminal="d8c32880-2388-543b-8c61-d9f865259254" include("term.jl")
    @require UnicodePlots="b8865327-cd53-5732-bb35-84acbb429228" include("uniplots.jl")
end

end # module


text:
"""
    struct MultiFieldFEFunction <: CellField
      # private fields
    end
"""
struct MultiFieldFEFunction{T<:MultiFieldCellField} <: FEFunction
  single_fe_functions::Vector{<:SingleFieldFEFunction}
  free_values::AbstractArray
  cell_dof_values::AbstractArray
  fe_space::MultiFieldFESpace
  multi_cell_field::T

  function MultiFieldFEFunction(
    free_values::AbstractVector,
    space::MultiFieldFESpace,
    single_fe_functions::Vector{<:SingleFieldFEFunction})

    multi_cell_field = MultiFieldCellField(map(i->i.cell_field,single_fe_functions))
    T = typeof(multi_cell_field)

    blocks = map(get_cell_dof_values,single_fe_functions)
    cell_axes_blocks = map(i->lazy_map(axes,i),blocks)
    cell_axes = lazy_map(_multifield_axes_dofs,cell_axes_blocks...)
    blockids = [(i,) for i in 1:length(single_fe_functions)]
    bsize = (length(single_fe_functions),)
    cell_dof_values = lazy_map(BlockArrayCooMap(bsize,blockids),cell_axes,blocks...)

    new{T}(
      single_fe_functions,
      free_values,
      cell_dof_values,
      space,
      multi_cell_field)
  end
end

CellData.get_data(f::MultiFieldFEFunction) = get_data(f.multi_cell_field)
CellData.get_triangulation(f::MultiFieldFEFunction) = get_triangulation(f.multi_cell_field)
CellData.DomainStyle(::Type{MultiFieldFEFunction{T}}) where T = DomainStyle(T)
FESpaces.get_free_values(f::MultiFieldFEFunction) = f.free_values
FESpaces.get_fe_space(f::MultiFieldFEFunction) = f.fe_space
FESpaces.get_cell_dof_values(f::MultiFieldFEFunction) = f.cell_dof_values

"""
    num_fields(m::MultiFieldFEFunction)
"""
num_fields(m::MultiFieldFEFunction) = length(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction) = iterate(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction,state) = iterate(m.single_fe_functions,state)
Base.getindex(m::MultiFieldFEFunction,field_id::Integer) = m.single_fe_functions[field_id]


text:
using MultivariateFunctions
using Test

# Run tests

println("Test of Univariate Functions.")
@time @test include("1_test_univariate.jl")
println("Test of MultivariateFunctions taking a date.")
@time @test include("2_test_dates.jl")
println("Test of Multivariate Functions.")
@time @test include("3_test_multivariate.jl")
println("Test of Univariate Piecewise Functions.")
@time @test include("4_piecewise_tests.jl")
println("Test of Multivariate Piecewise Functions.")
@time @test include("5_test_piecewise_multivariate.jl")
println("Test of Rootfinding and optimisation spline.")
@time @test include("6_optimisation_and_rootfinder.jl")
println("Test of interpolation with a univariate function.")
@time @test include("7_schumaker_test.jl")
println("Test of interpolation with a univariate function.")
@time @test include("8_interpolation_test.jl")
println("Test of chebyshev fitting with univariate and multivariate functions.")
@time @test include("9_chebyshev_tests.jl")
println("Test of ols regressions with univariate and multivariate functions.")
@time @test include("10_test_regressions.jl")
println("Test of high dimensional algorithms.")
@time @test include("11_test_Piecewise_ML_Algos.jl")


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

using ModelBlocks
using Test

println("Combo Test starting");

variables1 = Variables([
    Variable("X", 1, (0, 100), "kg", ""),
]);

parameters1 = Variables([
    Variable("p1", 0.0, (0, 100), "kg/s", ""),
    Variable("p2", 0.0, (0, 100), "kg/s", ""),
    Variable("p3", 0.0, (0, 100), "kg/s", ""),
]);

reactions1 = [
    GeneralRateReaction("R", [], ["X"], (t, v, p) -> p.p1 + p.p2 * v.X + p.p3 * t),
];

block1 = Block(variables1, parameters1, reactions1);

variables2 = Variables([
    Variable("Y", 1, (0, 100), "kg", ""),
]);

parameters2 = Variables([
    Variable("p3", 0.0, (0, 100), "kg/s", ""),
    Variable("p4", 2.0, (0, 100), "kg/s", ""),
]);

reactions2 = [
    GeneralRateReaction("R", [], ["Y"], (t, v, p) -> mod(t, 10) < 0.05 ? p.p4 : 0),
];

block2 = Block(variables2, parameters2, reactions2);

extraParameters = Variables([
    Variable("p5", 1.0, (0, 100), "kg/s", ""),
]);

combo = BlockCombo(
    [("first", block1), ("second", block2)],
    [
        ("first", "p1", (t, v, p) -> v.Y - v.X * p.p5)
    ],
    extraParameters
);

setTimeRange!(combo, 0:100);
setDiscontinuities!(combo, (p, timeRange) -> vcat([[t * 10, t * 10 + 0.05] for t in 0:10]...));

@time solution = runBlock(combo);

println("Combo Test finished");

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
#=
Author: Shunsuke Hori
=#

"""
This holds the results for Harrison Kreps.  In particular, it
accepts two matrices Qa and Qb and compares the single belief,
optimistic belief, and pessimistic belief prices
"""
type PriceHolder
    qa
    qb
    qpess
    qopt
    dividend_payoff
    qaprice
    qbprice
    qpessprice
    qoptprice
    optimisticprice
    pessimisticprice
    phat_a
    phat_b
end

function PriceHolder(qa, qb, dividend_payoff;
        beta=.75, max_iters = 10000, tolerance = 1e-16)

    # Create the Pessimistic and Optimistic Beliefs
    qpess = Array{Float64}(2, 2)
    qpess[1, :] = ifelse(qa[1, 2] < qb[1, 2],qa[1, :],qb[1, :])
    qpess[2, :] = ifelse(qa[2, 2] < qb[2, 2],qa[2, :],qb[2, :])
    qopt = Array{Float64}(2, 2)
    qopt[1, :] = ifelse(qa[1, 2] > qb[1, 2], qa[1, :], qb[1, :])
    qopt[2, :] = ifelse(qa[2, 2] > qb[2, 2], qa[2, :], qb[2, :])

    # Price everything
    p_singlebelief, p_optimistic, phat_a, phat_b, p_pessimistic=
        create_prices(qa, qb, qpess, qopt, dividend_payoff,
                                    beta,max_iters,tolerance)

    qaprice = p_singlebelief[1]
    qbprice = p_singlebelief[2]
    qpessprice = p_singlebelief[3]
    qoptprice = p_singlebelief[4]

    return PriceHolder(qa,
                        qb,
                        qpess,
                        qopt,
                        dividend_payoff,
                        qaprice,
                        qbprice,
                        qpessprice,
                        qoptprice,
                        p_optimistic,
                        p_pessimistic,
                        phat_a,
                        phat_b)
end

"""
Computes prices under all belief systems
"""
function create_prices(qa, qb, qpess, qopt, dividend_payoff,
                                bet,max_iters,tolerance)
    transitionmatrix = [qa, qb, qpess, qopt]
    # Single Belief Prices
    p_singlebelief = [price_singlebeliefs(q, dividend_payoff,beta=bet)
                         for q in transitionmatrix]

    # Compute Optimistic and Pessimistic beliefs
    p_optimistic, phat_a, phat_b =
        price_optimisticbeliefs([qa, qb],dividend_payoff,
                    beta=bet, max_iter=max_iters, tol=tolerance)
    p_pessimistic =
        price_pessimisticbeliefs([qa, qb], dividend_payoff,
                    beta=bet, max_iter=max_iters, tol=tolerance)

    return p_singlebelief, p_optimistic, phat_a, phat_b, p_pessimistic
end

function print_prices(ph::PriceHolder)
    qaprice, qbprice, qpessprice, qoptprice =
        ph.qaprice, ph.qbprice, ph.qpessprice, ph.qoptprice
    optimisticprice, pessimisticprice =
        ph.optimisticprice, ph.pessimisticprice
    phata, phatb = ph.phat_a, ph.phat_b
    println("The Single Belief Price Vectors are:")
    println(" P(Qa) = $(qaprice)\n P(Qb) = $(qbprice)\n P(Qopt) = $(qoptprice)\n P(Qpess) = $(qpessprice)\n")
    println("The Optimistic Belief Price Vector is:")
    println(" P(Optimistic) = $(optimisticprice)\n Phat(a) = $(phata)\n Phat(b) = $(phatb)\n")
    println("The Pessimistic Belief Price Vector is:")
    println(" P(Pessimistic) = $(pessimisticprice)")
end

ph=PriceHolder(qa, qb, [0, 1], beta=0.75)

print_prices(ph)


text:
# This script generates catalogues of objects in disk galaxies
# with the model of uniform distribution in the symmetry plane
# and Laplace distribution on the vertical axis, and it fits
# the same model to the generated data, retrieving parameters
# and evaluating confidence intervals

println('\n', " "^4, "> Loading the packages...")

using Distributions # Probability distributions
using LaTeXStrings # LaTeX strings
using Optim # Optimization
using Plots # Plotting
using Random # Random numbers
using Roots # Finding roots
using Zygote # Derivatives

# Use the PGFPlotsX backend for plots
# (this `@eval` trick is necessary to fool `DaemonMode`)
@eval Main begin
    using Plots
    pgfplotsx()
end

# Change some of the default parameters for plots
default(fontfamily = "Computer Modern", dpi = 300, size = (400, 400), legend = :topright)

# Define the paths to output directories
CURRENT_DIR = @__DIR__
ROOT_DIR = basename(CURRENT_DIR) == "scripts" ? dirname(CURRENT_DIR) : CURRENT_DIR
PLOTS_DIR = joinpath(ROOT_DIR, "plots")
TRACES_DIR = joinpath(ROOT_DIR, "traces")

# Make sure the needed directories exist
mkpath(PLOTS_DIR)
mkpath(TRACES_DIR)

"Probability density function of the Laplace distribution"
laplace_pdf(x, μ, b) = 1 / (2 * b) * exp(-abs(x - μ) / b)

"""
Negative log-likelihood function of the model
(the model is a Laplace distribution)
"""
function nll(z::Vector{F}, θ::Vector{F})::F where {F <: Real}
    μ, b = θ
    return -sum(@. log(laplace_pdf(z, μ, b)))
end

"""
Calculate the minimum of the negative log-likelihood
function while one of the parameters is frozen
"""
function nll_frozen(
    idx::Int,
    value::F,
    z::Vector{F},
    θ₀::Vector{F},
    θₗ::Vector{F},
    θᵤ::Vector{F},
) where {F <: Real}
    # Exclude the frozen parameter from the active ones
    θ₀ = [θ₀[1:idx - 1]; θ₀[idx + 1:end]]
    θₗ = [θₗ[1:idx - 1]; θₗ[idx + 1:end]]
    θᵤ = [θᵤ[1:idx - 1]; θᵤ[idx + 1:end]]

    # Recreate the function, handling a frozen parameter
    function nll_frozen_inner(θ::Vector{F}) where {F <: Real}
        return nll(z, [θ[1:idx - 1]; value; θ[idx:end]])
    end

    # Optimize the new negative log likelihood function
    res = Optim.optimize(
        nll_frozen_inner,
        θ -> Zygote.gradient(nll_frozen_inner, θ)[1],
        θₗ,
        θᵤ,
        θ₀,
        Fminbox(LBFGS()),
        Optim.Options(
            show_trace=false,
            extended_trace=true,
            store_trace=true,
        );
        inplace=false,
    )

    return res.minimum
end

"Plot a histogram of the height sample"
function plot_histogram_height(z, μₜ, bₜ, μ, b, dir)
    # Plot the histogram of the height sample
    histogram(
        z;
        label = "",
        title = "",
        xlabel = L"z",
        ylabel = L"f(z)",
        normalize = true,
        yminorticks = 5,
    )

    # Plot the probability functions over the histogram
    plot!((x) -> laplace_pdf(x, μₜ, bₜ), label = "true", lw = 1.5)
    plot!((x) -> laplace_pdf(x, μ, b), label = "fit", lw = 1.5)

    # Save the figure as PDF and PNG
    savefig(joinpath(dir, "Histogram Z.pdf"))
end

"Create the objective function profiles for all of the parameters"
function profiles(
    z::Vector{F},
    L₀::F,
    θ::Vector{F},
    θ₀::Vector{F},
    θₗ::Vector{F},
    θᵤ::Vector{F},
    labels::Vector{Vector{String}},
    plots_dir::String,
) where {F <: Real}
    for (i, x₀, (label, tex_label)) in zip(1:length(θ), θ, labels)
        println(" "^6, "> Evaluating the confidence intervals of the $label parameter...")

        # Create an alias function for the frozen parameter
        nll_frozen_with(x) = nll_frozen(i, x, z, θ₀, θₗ, θᵤ)

        # Define a function whose roots define confidence intervals
        confidence_intervals(x) = nll_frozen_with(x) - L₀ - 0.5

        # Evaluate the confidence intervals
        r1, r2 = find_zeros(confidence_intervals, x₀ - 0.2, x₀ + 0.2)
        σr1, σr2 = x₀ - r1, r2 - x₀

        # Compute the mean of the sigmas
        σr = (σr1 + σr2) / 2

        # Print the results
        println(
            '\n',
            " "^8, "$label = $(x₀)\n",
            " "^8, "$(label)₋ₕ = $(r1)\n",
            " "^8, "$(label)₊h = $(r2)\n",
            " "^8, "σ$(label)₋ₕ = $(σr1)\n",
            " "^8, "σ$(label)₊ₕ = $(σr2)\n",
        )

        println(" "^6, "> Plotting the objective function profile for the $label parameter...")

        # Plot the profile
        plot(nll_frozen_with, r1 - σr1, r2 + σr2; label="", xlabel=latexstring(tex_label), ylabel=latexstring("L_p($tex_label)"))

        # Add vertical lines to the plot
        plot!([r1, r1], [L₀ - 0.2, L₀ + 0.5]; label="", linestyle=:dash)
        plot!([x₀, x₀], [L₀ - 0.2, L₀]; label="", linestyle=:dash)
        plot!([r2, r2], [L₀ - 0.2, L₀ + 0.5]; label="", linestyle=:dash)

        # Add points to the plot
        scatter!([r1], [L₀ + 0.5]; label="")
        scatter!([x₀], [L₀]; label="")
        scatter!([r2,], [L₀ + 0.5]; label="")

        # Add horizontal lines to the plot
        hline!([L₀ + 0.5]; label="", linestyle=:dash)
        hline!([L₀]; label="", linestyle=:dash)

        # Add annotations to the plot
        rounded = round.([r1, x₀, r2]; digits=2)
        annotate!(
            [
                (r1 + σr * 0.2, L₀ - 0.1, text("$(rounded[1])", 8, "Computer Modern")),
                (x₀, L₀ + 0.1, text("$(rounded[2])", 8, "Computer Modern")),
                (r2 + σr * 0.2, L₀ - 0.1, text("$(rounded[3])", 8, "Computer Modern")),
                (r1 - σr * 0.22, L₀ - 0.1, text("$(round(rounded[1] - rounded[2]; digits=2))", 8, "Computer Modern")),
                (r2 - σr * 0.25, L₀ - 0.103, text("+$(round(rounded[3] - rounded[2]; digits=2))", 8, "Computer Modern")),
                (x₀ - σr1 * 1.80, L₀ + 0.06, text(L"L_0", 8, "Computer Modern")),
                (x₀ - σr1 * (1.80 - 0.21), L₀ + 0.5 + 0.06, text(L"L_0 + 1/2", 8, "Computer Modern")),
            ]
        )

        # Save the figure
        savefig(joinpath(plots_dir, "$label profile.pdf"))
    end
end

# Define a set of initial parameters
θ₀ = [0.0, 0.5]

# Define the lower and upper boundaries, respectively
θₗ = [-Inf, 0]
θᵤ = [Inf, Inf]

"""
Generate a catalogue for a disk with uniform distribution in
the symmetry plane and Laplace distribution on the vertical axis
"""
function fit_uniform_disk(seed; N::Int = 1000, μₜ::F = 2.0, bₜ::F = 0.9) where {F <: Real}
    println(" "^4, "> Fitting a sample (seed = $seed, N = $N) of the uniform disk distribution...")

    # Set a seed for the random generator
    Random.seed!(seed)

    # Define the output directories
    name = "$seed, $N"
    plots_dir = joinpath(PLOTS_DIR, "fit_uniform_disk", name)
    traces_dir = joinpath(TRACES_DIR, name)

    # Make sure the needed directories exists
    mkpath(plots_dir)
    mkpath(traces_dir)

    # Generate a sample of heights from the Laplace distribution
    z = rand(Laplace(μₜ, bₜ), N)

    # Link the negative log-likelihood function with data
    nll_inner(θ) = nll(z, θ)

    # Optimize the negative log likelihood function
    res = Optim.optimize(
        nll_inner,
        θ -> Zygote.gradient(nll_inner, θ)[1],
        θₗ,
        θᵤ,
        θ₀,
        Fminbox(LBFGS()),
        Optim.Options(
            show_trace=false,
            extended_trace=true,
            store_trace=true,
        );
        inplace=false,
    )

    # Unpack the results
    θ = res.minimizer
    μ, b = res.minimizer
    L₀ = res.minimum

    # Save the trace and the results
    open(joinpath(traces_dir, "fit.trace"), "w") do io
        println(io, res.trace)
        println(
            io,
            " * Parameters:\n",
            "    μ = $(μ)\n",
            "    b = $(b)\n",
        )
        show(io, res)
    end

    # Print results
    println(
        '\n',
        " "^6, "μₜ = $(μₜ)\n",
        " "^6, "bₜ = $(bₜ)\n",
        " "^6, "μ = $(μ)\n",
        " "^6, "b = $(b)\n",
    )

    println(" "^6, "> Plotting the histogram...")

    # Plot the histogram of the height sample
    plot_histogram_height(z, μₜ, bₜ, μ, b, plots_dir)

    # Create the objective function profiles for all of the parameters
    profiles(z, L₀, θ, θ₀, θₗ, θᵤ, [["μ", "\\mu"], ["b", "b"]], plots_dir)

    print("\n")
    sleep(0.5)
end

# Plot samples from the uniform disk distribution
fit_uniform_disk(1, N=100)
fit_uniform_disk(1)
fit_uniform_disk(1, N=10000)
fit_uniform_disk(2)
fit_uniform_disk(3)
fit_uniform_disk(4)
fit_uniform_disk(5)


text:
"""
Atomic updateBTG function

Incorporate new data into btg object, namely triple (x_new, Fx_new, y_new)
"""
function updateBTG!(btg::btg, x_new::Array{T, 2} where T<:Real, Fx_new::Array{T, 2} where T<:Real, y_new:: Union{T, G} where T<:Real where G<:Array{F, 1} where F<:Real)
    @assert typeof(btg.trainingData) <: extensible_trainingData
    update!(btg.trainingData, x_new, Fx_new, y_new) #update training data once, and importantly, first

    for key in keys(btg.train_buffer_dict)
        update!(btg.train_buffer_dict[key], btg.trainingData) #update train_buffer_dict, so we don't recompute cholesky decompositions from scratch
    end

    for key in keys(btg.λbuffer_dict)
        update!(btg.λbuffer_dict[key], y_new, btg.g)
    end

    for key in keys(btg.θλbuffer_dict)
        update!(btg.θλbuffer_dict[key], btg.λbuffer_dict[key[2]], btg.train_buffer_dict[key[1]], btg.trainingData) #update θλbuffer
    end
    btg.n = btg.n + 1
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# ------------------------------------------------------------------------------------------
# # 多重派发
#
# **多重派发**是这篇notebook要探索的Julia的一个关键特性。
#
# 它有助于快速开发。它也让软件可拓展、可设计和更好把玩。
#
# 它可能是并行计算取得重大突破的预兆。
#
# ## 大纲
# 1. 罗马数字
# 2. 函数
# 3. 并行计算
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 1. 罗马数字（玩玩）
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 我们来定义一个**新结构**表示罗马数字。为了方便写代码，我们将只考虑 0 到 9。
#
# **练习**：拓展到更大的数字。（别忘了罗马数字是十进制计数系统！）
# ------------------------------------------------------------------------------------------

struct Roman
    n::Int
end

Base.show(io::IO, r::Roman) = print(io, 'ⅰ' + (r.n - 1) % 10 )  # 漂亮的显示；'ⅰ'是一个Unicode罗马数字

# ------------------------------------------------------------------------------------------
# 我们可以像这样创建一个这个类型的对象：
# ------------------------------------------------------------------------------------------

Roman(4)

typeof.([5 5.0 Roman(5) "Five" '5'  5//1])

# ------------------------------------------------------------------------------------------
# 用罗马数字漂亮地打印数组：
# ------------------------------------------------------------------------------------------

x = [7 1 2 5 8 9]
Roman.(x)   # equivalent to map(Roman, x)  or  [Roman(w) for w in x]

# ------------------------------------------------------------------------------------------
# 要是能像普通数字那样给罗马数字做加法就好了：
# ------------------------------------------------------------------------------------------

Roman(4) + Roman(5)

# ------------------------------------------------------------------------------------------
# 但 Julia 不知道该怎么做。
# 我们可以通过 `import` 函数 `+` 并拓展它的定义来教 Julia 给罗马数字做加法：
# ------------------------------------------------------------------------------------------

import Base: +, *

+(a::Roman, b::Roman) = Roman(a.n + b.n)

Roman(4) + Roman(5)

# ------------------------------------------------------------------------------------------
# 这样就给函数 `+` **添加了一个新方法**：
# ------------------------------------------------------------------------------------------

methods(+)

@which Roman(4) + Roman(5)

import Base.*
*(i::Roman, j::Roman) = Roman(i.n * j.n)                     # Multiply like a Roman

Roman(3) * Roman(2)

Roman.(1:3) .* [Roman(1) Roman(2) Roman(3)]

# ------------------------------------------------------------------------------------------
# 但是我们的乘法仍有一些问题
# ------------------------------------------------------------------------------------------

Roman(3) * 2

# mytimes 函数目前基于类型来判断该做什么，很复杂
# 不用担心，更好的方法马上就来
function mytimes(i,j)
  if isa(i,Roman) & isa(j,Number)
        return  fill(1, i.n, j)   # i by j matrix with ones
    elseif    isa(i,Number) & isa(j,Roman)
        return "😄"^ (i*j.n)   #  i * j happy faces
    else
        return("I Don't know")
    end
end

mytimes(4,Roman(3)) # 12个笑脸

mytimes(Roman(4),3) # 4x3的全1矩阵

# ------------------------------------------------------------------------------------------
# 最简单的实现方式是明确定义一个 `Roman` 和一个数字的乘法。我们可以按我们的想法来定义：
# ------------------------------------------------------------------------------------------

*(i::Number, j::Roman) = "😄"^ (i*j.n)    # i * j个笑脸

*(i::Roman, j::Number) =   fill(1, i.n, j) # i * j的矩阵

3 * Roman(3) # 9个笑脸

Roman(3) * 5  # 3*5的全1矩阵

t(x::Roman,y::Roman) = x.n * y.n

t(Roman(5),Roman(4))

# 注意它的汇编代码是多紧凑！
@code_native t(Roman(2),Roman(4))

# ------------------------------------------------------------------------------------------
# ## 函数
# ------------------------------------------------------------------------------------------

import Base: *, +, ^

*(α::Number,   g::Function) = x -> α * g(x)   # 标量乘以函数

*(f::Function, λ::Number)   = x -> f(λ * x)   # Scale the argument

*(f::Function, g::Function) = x -> f(g(x))    # 复合函数  -- 滥用符号！  在Julia0.6中使用 \circ

^(f::Function, n::Integer) = n == 1 ? f : f*f^(n-1) # 一个通过递归乘法实现的天真的求幂算法

+(f::Function, g::Function) = x -> f(x) + g(x)

# ------------------------------------------------------------------------------------------
# 举个例子，定义成这样的指数函数
#
# $$\exp(x) = \sum_{n=0}^\infty \frac{1}{n!} x^n.$$
#
# 我们可以把它看成是这样的函数：
#
# $$\exp = \sum_{n=0}^\infty \frac{1}{n!} \mathrm{pow}_n,$$
#
# 其中 $\mathrm{pow}_n(x) = x^n$.
#
# （开始用数字模糊符号！）
# ------------------------------------------------------------------------------------------

pow(n) = x -> x^n

myexp = sum(1/factorial(big(n)) * pow(n) for n in 0:100)   # 效率低的泰勒级数！

[myexp(1); exp(1); exp(big(1))]

f = x -> x^2
f(10)

g = 3f
g(10)

(f^2)(10)  # 因为我们已经定义了函数乘法为符合函数

using Plots;
gr()

x = pi*(0:0.001:4)

plot( x,  sin.(x),        c="black",            label="Fun")
plot!(x, (12*sin).(x),    c="green",            label="Num * Fun")
plot!(x, (sin*12).(x),    c="red",   alpha=0.9, label="Fun * Num")
plot!(x, (5*sin*exp).(x), c="blue",  alpha=0.2, label="Num * Fun * Fun")

plot([12*sin, sin*12, 5*sin*exp], 0:.01:4π, α=[1 .9 .2], c=[:green :red :blue])

# ------------------------------------------------------------------------------------------
# > “我很讨厌 $sin^2 \phi$，虽然 Laplace 用过它；
# > 应该要担心 $sin^2 \phi$ 可能会引起歧义，如果说 $sin(\phi^2)$ 就不会或者说几乎不会引起歧义了，那么我们应该写 $(sin \phi)^2$，而不是
# $sin^2 \phi$，以此类推 $sin^2 \phi$应该指的是 $sin(sin \phi)$。”
# >
# > —— Gauss
# ------------------------------------------------------------------------------------------

x=(0:.01:2) * pi;

plot(x, (sin^2).(x), c="blue")     # 乘方能工作，y=sin(sin(x))，Gauss会开心的！
plot!(x, sin.(x).^2, c="red")

# ------------------------------------------------------------------------------------------
# ## 练习
# ------------------------------------------------------------------------------------------

h(a, b::Any) = "fallback"
h(a::Number, b::Number) = "a and b are both numbers"
h(a::Number, b) = "a is a number"
h(a, b::Number) = "b is a number"
h(a::Integer, b::Integer) = "a and b are both integers"

# 试着把玩 h，尝试使用 h 的 5 种方法
# 在此作答

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
##############################################################################
#
# REFERENCES: Wasserman, "All of Statistics"
#
# NOTES: Fails CDF/Quantile matching test
#
##############################################################################

immutable Bernoulli <: DiscreteUnivariateDistribution
    p0::Float64
    p1::Float64

    function Bernoulli(p::Real)
        zero(p) <= p <= one(p) || error("prob must be in [0,1]")
        new(1.0 - p, float64(p))
    end
end

Bernoulli() = Bernoulli(0.5)

min(d::Bernoulli) = 0
max(d::Bernoulli) = 1

cdf(d::Bernoulli, q::Real) = q >= zero(q) ? (q >= one(q) ? 1.0 : d.p0) : 0.

function entropy(d::Bernoulli)
    p0 = d.p0
    p1 = d.p1
    p0 == 0. || p0 == 1. ? 0. : -(p0 * log(p0) + p1 * log(p1))
end

insupport(::Bernoulli, x::Real) = (x == zero(x)) || (x == one(x))
insupport(::Type{Bernoulli}, x::Real) = (x == zero(x)) || (x == one(x))

mean(d::Bernoulli) = d.p1

var(d::Bernoulli) = d.p0 * d.p1

skewness(d::Bernoulli) = (d.p0 - d.p1) / sqrt(d.p0 * d.p1)

kurtosis(d::Bernoulli) = 1.0 / (d.p0 * d.p1) - 6.0

median(d::Bernoulli) = d.p1 < 0.5 ? 0.0 : 1.0

mgf(d::Bernoulli, t::Real) = d.p0 + d.p1 * exp(t)

cf(d::Bernoulli, t::Real) = d.p0 + d.p1 * exp(im * t)

mode(d::Bernoulli) = d.p1 > 0.5 ? 1 : 0

function modes(d::Bernoulli)
    d.p1 < 0.5 ? [0] :
    d.p1 > 0.5 ? [1] : [0, 1]
end

pdf(d::Bernoulli, x::Real) = x == zero(x) ? d.p0 : x == one(x) ? d.p1 : 0.0

quantile(d::Bernoulli, p::Real) = zero(p) <= p <= one(p) ? (p <= d.p0 ? 0 : 1) : NaN

rand(d::Bernoulli) = rand() > d.p1 ? 0 : 1

## MLE fitting

immutable BernoulliStats <: SufficientStats
    cnt0::Float64
    cnt1::Float64

    BernoulliStats(c0::Real, c1::Real) = new(float64(c0), float64(c1))
end

fit_mle(::Type{Bernoulli}, ss::BernoulliStats) = Bernoulli(ss.cnt1 / (ss.cnt0 + ss.cnt1))

function suffstats{T<:Integer}(::Type{Bernoulli}, x::Array{T})
    n0 = 0
    n1 = 0
    for xi in x
        if xi == 0
            n0 += 1
        elseif xi == 1
            n1 += 1
        else
            throw(DomainError())
        end
    end
    BernoulliStats(n0, n1)
end

function suffstats{T<:Integer}(::Type{Bernoulli}, x::Array{T}, w::Array{Float64})
    n = length(x)
    if length(w) != n
        throw(ArgumentError("Inconsistent argument dimensions."))
    end

    n0 = 0.0
    n1 = 0.0
    for i = 1:n
        xi = x[i]
        if xi == 0
            n0 += w[i]
        elseif xi == 1
            n1 += w[i]
        else
            throw(DomainError())
        end
    end
    BernoulliStats(n0, n1)
end


text:
using SQLite
using Test, Dates, Random, WeakRefStrings, Tables, DBInterface

import Base: +, ==

mutable struct Point{T}
    x::T
    y::T
end

==(a::Point, b::Point) = a.x == b.x && a.y == b.y

mutable struct Point3D{T<:Number}
    x::T
    y::T
    z::T
end

==(a::Point3D, b::Point3D) = a.x == b.x && a.y == b.y && a.z == b.z
+(a::Point3D, b::Point3D) = Point3D(a.x + b.x, a.y + b.y, a.z + b.z)

triple(x) = 3x

function add4(q)
    q + 4
end

mult(args...) = *(args...)
str2arr(s) = Vector{UInt8}(s)
doublesum_step(persist, current) = persist + current
doublesum_final(persist) = 2 * persist
mycount(p, c) = p + 1
bigsum(p, c) = p + big(c)
sumpoint(p::Point3D, x, y, z) = p + Point3D(x, y, z)

"""
    setup_clean_test_db() do db
        @test db isa SQLite.DB
    end

Copy test sqlite file to temp directory, path `test.sqlite`,
overwrite file if necessary and set permissions.
"""
function setup_clean_test_db(f::Function, args...)
    dbfile = joinpath(dirname(pathof(SQLite)), "../test/Chinook_Sqlite.sqlite")
    tmp_dir = mktempdir()
    test_dbfile = joinpath(tmp_dir, "test.sqlite")

    cp(dbfile, test_dbfile; force = true)
    chmod(test_dbfile, 0o777)

    db = SQLite.DB(test_dbfile)
    try
        f(db)
    finally
        close(db)
        rm(tmp_dir, recursive = true)
    end
end

@testset "basics" begin

    @testset "Julia to SQLite3 type conversion" begin
        @test SQLite.sqlitetype(Int) == "INT NOT NULL"
        @test SQLite.sqlitetype(Union{Float64,Missing}) == "REAL"
        @test SQLite.sqlitetype(String) == "TEXT NOT NULL"
        @test SQLite.sqlitetype(Symbol) == "BLOB NOT NULL"
        @test SQLite.sqlitetype(Missing) == "NULL"
        @test SQLite.sqlitetype(Nothing) == "NULL"
        @test SQLite.sqlitetype(Any) == "BLOB"
    end

    @testset "DB Connection" begin
        dbfile = joinpath(dirname(pathof(SQLite)), "../test/Chinook_Sqlite.sqlite")
        @test SQLite.DB(dbfile) isa SQLite.DB
        con = DBInterface.connect(SQLite.DB, dbfile)
        @test con isa SQLite.DB
        DBInterface.close!(con)
    end

    @testset "Issue #207: 32 bit integers" begin
        setup_clean_test_db() do db
            ds =
                DBInterface.execute(db, "SELECT RANDOM() as a FROM Track LIMIT 1") |>
                columntable
            @test ds.a[1] isa Int64

        end
    end

    @testset "Regular SQLite Tests" begin
        setup_clean_test_db() do db
            @test_throws SQLiteException DBInterface.execute(db, "just some syntax error")
            # syntax correct, table missing
            @test_throws SQLiteException DBInterface.execute(
                db,
                "SELECT name FROM sqlite_nomaster WHERE type='table';",
            )
        end
    end

    @testset "close!(query)" begin
        setup_clean_test_db() do db
            qry = DBInterface.execute(db, "SELECT name FROM sqlite_master WHERE type='table';")
            DBInterface.close!(qry)
            DBInterface.close!(qry) # test it doesn't throw on double-close
        end
    end

    @testset "Query tables" begin
        setup_clean_test_db() do db
            ds =
                DBInterface.execute(db, "SELECT name FROM sqlite_master WHERE type='table';") |>
                columntable
            @test length(ds) == 1
            @test keys(ds) == (:name,)
            @test length(ds.name) == 11

            results1 = SQLite.tables(db)
            @test isequal(ds, results1)
        end
    end

    @testset "DBInterface.execute([f])" begin
        setup_clean_test_db() do db

            # pipe approach
            results = DBInterface.execute(db, "SELECT * FROM Employee;") |> columntable
            @test length(results) == 15
            @test length(results[1]) == 8
            # callable approach
            @test isequal(
                DBInterface.execute(columntable, db, "SELECT * FROM Employee"),
                results,
            )
            employees_stmt = DBInterface.prepare(db, "SELECT * FROM Employee")
            @test isequal(columntable(DBInterface.execute(employees_stmt)), results)
            @test isequal(DBInterface.execute(columntable, employees_stmt), results)
            @testset "throwing from f()" begin
                f(::SQLite.Query) = error("I'm throwing!")
                @test_throws ErrorException DBInterface.execute(f, employees_stmt)
                @test_throws ErrorException DBInterface.execute(f, db, "SELECT * FROM Employee")
            end
            DBInterface.close!(employees_stmt)
        end
    end

    @testset "isempty(::Query)" begin
        setup_clean_test_db() do db

            @test !DBInterface.execute(isempty, db, "SELECT * FROM Employee")
            @test DBInterface.execute(
                isempty,
                db,
                "SELECT * FROM Employee WHERE FirstName='Joanne'",
            )
        end
    end

    @testset "empty query has correct schema and return type" begin
        setup_clean_test_db() do db
            empty_scheme = DBInterface.execute(
                Tables.schema,
                db,
                "SELECT * FROM Employee WHERE FirstName='Joanne'",
            )
            all_scheme = DBInterface.execute(
                Tables.schema,
                db,
                "SELECT * FROM Employee WHERE FirstName='Joanne'",
            )
            @test empty_scheme.names == all_scheme.names
            @test all(ea -> ea[1] <: ea[2], zip(empty_scheme.types, all_scheme.types))

            empty_tbl = DBInterface.execute(
                columntable,
                db,
                "SELECT * FROM Employee WHERE FirstName='Joanne'",
            )
            all_tbl = DBInterface.execute(columntable, db, "SELECT * FROM Employee")
            @test propertynames(empty_tbl) == propertynames(all_tbl)
            @test all(
                col ->
                    eltype(empty_tbl[col]) == Missing ||
                        eltype(empty_tbl[col]) >: eltype(all_tbl[col]),
                propertynames(all_tbl),
            )
        end
    end

    @testset "Create table, run commit/rollback tests" begin
        setup_clean_test_db() do db
            DBInterface.execute(db, "create table temp as select * from album")
            DBInterface.execute(db, "alter table temp add column colyear int")
            DBInterface.execute(db, "update temp set colyear = 2014")
            r = DBInterface.execute(db, "select * from temp limit 10") |> columntable
            @test length(r) == 4 && length(r[1]) == 10
            @test all(==(2014), r[4])

            @test_throws SQLiteException SQLite.rollback(db)
            @test_throws SQLiteException SQLite.commit(db)

            SQLite.transaction(db)
            DBInterface.execute(db, "update temp set colyear = 2015")
            SQLite.rollback(db)
            r = DBInterface.execute(db, "select * from temp limit 10") |> columntable
            @test all(==(2014), r[4])

            SQLite.transaction(db)
            DBInterface.execute(db, "update temp set colyear = 2015")
            SQLite.commit(db)
            r = DBInterface.execute(db, "select * from temp limit 10") |> columntable
            @test all(==(2015), r[4])
        end
    end

    @testset "Dates" begin
        setup_clean_test_db() do db
            DBInterface.execute(db, "create table temp as select * from album")
            DBInterface.execute(db, "alter table temp add column dates blob")
            stmt = DBInterface.prepare(db, "update temp set dates = ?")
            DBInterface.execute(stmt, (Date(2014, 1, 1),))

            r = DBInterface.execute(db, "select * from temp limit 10") |> columntable
            @test length(r) == 4 && length(r[1]) == 10
            @test typeof(r[4][1]) == Date
            @test all(Bool[x == Date(2014, 1, 1) for x in r[4]])
            DBInterface.execute(db, "drop table temp")

            rng = Dates.Date(2013):Dates.Day(1):Dates.Date(2013, 1, 5)
            dt = (i = collect(rng), j = collect(rng))
            tablename = dt |> SQLite.load!(db, "temp")
            r = DBInterface.execute(db, "select * from $tablename") |> columntable
            @test length(r) == 2 && length(r[1]) == 5
            @test all([i for i in r[1]] .== collect(rng))
            @test all([typeof(i) for i in r[1]] .== Dates.Date)
            SQLite.drop!(db, "$tablename")
        end
    end

    @testset "Prepared Statements" begin
        setup_clean_test_db() do db

            DBInterface.execute(db, "CREATE TABLE temp AS SELECT * FROM Album")
            r = DBInterface.execute(db, "SELECT * FROM temp LIMIT ?", [3]) |> columntable
            @test length(r) == 3 && length(r[1]) == 3
            r =
                DBInterface.execute(db, "SELECT * FROM temp WHERE Title LIKE ?", ["%time%"]) |>
                columntable
            @test r[1] == [76, 111, 187]
            DBInterface.execute(
                db,
                "INSERT INTO temp VALUES (?1, ?3, ?2)",
                [0, 0, "Test Album"],
            )
            r = DBInterface.execute(db, "SELECT * FROM temp WHERE AlbumId = 0") |> columntable
            @test r[1][1] == 0
            @test r[2][1] == "Test Album"
            @test r[3][1] == 0
            SQLite.drop!(db, "temp")

            DBInterface.execute(db, "CREATE TABLE temp AS SELECT * FROM Album")
            r = DBInterface.execute(db, "SELECT * FROM temp LIMIT :a", (a = 3,)) |> columntable
            @test length(r) == 3 && length(r[1]) == 3
            r = DBInterface.execute(db, "SELECT * FROM temp LIMIT :a", a = 3) |> columntable
            @test length(r) == 3 && length(r[1]) == 3
            r =
                DBInterface.execute(
                    db,
                    "SELECT * FROM temp WHERE Title LIKE @word",
                    (word = "%time%",),
                ) |> columntable
            @test r[1] == [76, 111, 187]
            DBInterface.execute(
                db,
                "INSERT INTO temp VALUES (@lid, :title, \$rid)",
                (rid = 1, lid = 0, title = "Test Album"),
            )
            DBInterface.execute(
                db,
                "INSERT INTO temp VALUES (@lid, :title, \$rid)",
                rid = 3,
                lid = 400,
                title = "Test2 Album",
            )
            r =
                DBInterface.execute(db, "SELECT * FROM temp WHERE AlbumId IN (0, 400)") |>
                columntable
            @test r[1] == [0, 400]
            @test r[2] == ["Test Album", "Test2 Album"]
            @test r[3] == [1, 3]
            SQLite.drop!(db, "temp")

            SQLite.register(db, SQLite.regexp, nargs = 2, name = "regexp")
            r =
                DBInterface.execute(
                    db,
                    SQLite.@sr_str(
                        "SELECT LastName FROM Employee WHERE BirthDate REGEXP '^\\d{4}-08'"
                    )
                ) |> columntable
            @test r[1][1] == "Peacock"

            SQLite.register(db, identity, nargs = 1, name = "identity")
            r =
                DBInterface.execute(
                    db,
                    """SELECT identity("abc") as x, "abc" == identity("abc") as cmp""",
                ) |> columntable
            @test first(r.x) == "abc"
            @test first(r.cmp) == 1

            @test_throws AssertionError SQLite.register(db, triple, nargs = 186)
            SQLite.register(db, triple, nargs = 1)
            r =
                DBInterface.execute(
                    db,
                    "SELECT triple(Total) FROM Invoice ORDER BY InvoiceId LIMIT 5",
                ) |> columntable
            s =
                DBInterface.execute(
                    db,
                    "SELECT Total FROM Invoice ORDER BY InvoiceId LIMIT 5",
                ) |> columntable
            for (i, j) in zip(r[1], s[1])
                @test abs(i - 3 * j) < 0.02
            end
        end
    end

    @testset "Register functions" begin
        setup_clean_test_db() do db
            SQLite.@register db add4
            r = DBInterface.execute(db, "SELECT add4(AlbumId) FROM Album") |> columntable
            s = DBInterface.execute(db, "SELECT AlbumId FROM Album") |> columntable
            @test r[1][1] == s[1][1] + 4

            SQLite.@register db mult
            r = DBInterface.execute(db, "SELECT GenreId, UnitPrice FROM Track") |> columntable
            s =
                DBInterface.execute(db, "SELECT mult(GenreId, UnitPrice) FROM Track") |>
                columntable
            @test (r[1][1] * r[2][1]) == s[1][1]
            t =
                DBInterface.execute(db, "SELECT mult(GenreId, UnitPrice, 3, 4) FROM Track") |>
                columntable
            @test (r[1][1] * r[2][1] * 3 * 4) == t[1][1]

            SQLite.@register db sin
            u =
                DBInterface.execute(db, "select sin(milliseconds) from track limit 5") |>
                columntable
            @test all(-1 .< convert(Vector{Float64}, u[1]) .< 1)

            SQLite.register(db, hypot; nargs = 2, name = "hypotenuse")
            v =
                DBInterface.execute(
                    db,
                    "select hypotenuse(Milliseconds,bytes) from track limit 5",
                ) |> columntable
            @test [round(Int, i) for i in v[1]] == [11175621, 5521062, 3997652, 4339106, 6301714]

            SQLite.@register db str2arr
            r =
                DBInterface.execute(db, "SELECT str2arr(LastName) FROM Employee LIMIT 2") |>
                columntable
            @test r[1][2] == UInt8[0x45, 0x64, 0x77, 0x61, 0x72, 0x64, 0x73]

            SQLite.@register db big
            r = DBInterface.execute(db, "SELECT big(5)") |> columntable
            @test r[1][1] == big(5)
            @test typeof(r[1][1]) == BigInt

            SQLite.register(db, 0, doublesum_step, doublesum_final, name = "doublesum")
            r = DBInterface.execute(db, "SELECT doublesum(UnitPrice) FROM Track") |> columntable
            s = DBInterface.execute(db, "SELECT UnitPrice FROM Track") |> columntable
            @test abs(r[1][1] - 2 * sum(convert(Vector{Float64}, s[1]))) < 0.02

            SQLite.register(db, 0, mycount)
            r =
                DBInterface.execute(db, "SELECT mycount(TrackId) FROM PlaylistTrack") |>
                columntable
            s =
                DBInterface.execute(db, "SELECT count(TrackId) FROM PlaylistTrack") |>
                columntable
            @test r[1][1] == s[1][1]

            SQLite.register(db, big(0), bigsum)
            r =
                DBInterface.execute(db, "SELECT bigsum(TrackId) FROM PlaylistTrack") |>
                columntable
            s = DBInterface.execute(db, "SELECT TrackId FROM PlaylistTrack") |> columntable
            @test r[1][1] == big(sum(convert(Vector{Int}, s[1])))

            DBInterface.execute(db, "CREATE TABLE points (x INT, y INT, z INT)")
            DBInterface.execute(db, "INSERT INTO points VALUES (?, ?, ?)", (1, 2, 3))
            DBInterface.execute(db, "INSERT INTO points VALUES (?, ?, ?)", (4, 5, 6))
            DBInterface.execute(db, "INSERT INTO points VALUES (?, ?, ?)", (7, 8, 9))

            SQLite.register(db, Point3D(0, 0, 0), sumpoint)
            r = DBInterface.execute(db, "SELECT sumpoint(x, y, z) FROM points") |> columntable
            @test r[1][1] == Point3D(12, 15, 18)
            SQLite.drop!(db, "points")

            db2 = DBInterface.connect(SQLite.DB)
            DBInterface.execute(db2, "CREATE TABLE tab1 (r REAL, s INT)")

            @test_throws SQLiteException SQLite.drop!(db2, "nonexistant")
            # should not throw anything
            SQLite.drop!(db2, "nonexistant", ifexists = true)
            # should drop "tab2"
            SQLite.drop!(db2, "tab2", ifexists = true)
            @test !in("tab2", SQLite.tables(db2)[1])

            SQLite.drop!(db, "sqlite_stat1", ifexists = true)
            tables = SQLite.tables(db)
            @test length(tables[1]) == 11
        end
    end

    @testset "Remove Duplicates" begin
        setup_clean_test_db() do db
            db = SQLite.DB() #In case the order of tests is changed
            dt = (ints = Int64[1, 1, 2, 2, 3], strs = ["A", "A", "B", "C", "C"])
            tablename = dt |> SQLite.load!(db, "temp")
            SQLite.removeduplicates!(db, "temp", ["ints", "strs"]) #New format
            dt3 = DBInterface.execute(db, "Select * from temp") |> columntable
            @test dt3[1][1] == 1
            @test dt3[2][1] == "A"
            @test dt3[1][2] == 2
            @test dt3[2][2] == "B"
            @test dt3[1][3] == 2
            @test dt3[2][3] == "C"
        end
    end

    @testset "Issue #104: bind!() fails with named parameters" begin
        db = SQLite.DB() #In case the order of tests is changed
        DBInterface.execute(db, "CREATE TABLE IF NOT EXISTS tbl(a  INTEGER);")
        stmt = DBInterface.prepare(db, "INSERT INTO tbl (a) VALUES (@a);")
        SQLite.bind!(stmt, "@a", 1)
        SQLite.clear!(stmt)
    end

    @testset "SQLite to Julia type conversion" begin
        binddb = SQLite.DB()
        DBInterface.execute(
            binddb,
            "CREATE TABLE temp (n NULL, i1 INT, i2 integer,
                                f1 REAL, f2 FLOAT, f3 NUMERIC,
                                s1 TEXT, s2 CHAR(10), s3 VARCHAR(15), s4 NVARCHAR(5),
                                d1 DATETIME, ts1 TIMESTAMP,
                                b BLOB,
                                x1 UNKNOWN1, x2 UNKNOWN2)",
        )
        DBInterface.execute(
            binddb,
            "INSERT INTO temp VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)",
            [
                missing,
                Int64(6),
                Int64(4),
                6.4,
                6.3,
                Int64(7),
                "some long text",
                "short text",
                "another text",
                "short",
                "2021-02-21",
                "2021-02-12::1532",
                b"bytearray",
                "actually known",
                435,
            ],
        )
        rr = (;) # just to have the var declared
        @test_logs(
            (
                :warn,
                "Unsupported SQLite declared type UNKNOWN1, falling back to String type",
            ),
            (
                :warn,
                "Unsupported SQLite declared type UNKNOWN2, falling back to $(Int64) type",
            ),
            rr = DBInterface.execute(rowtable, binddb, "SELECT * FROM temp")
        )
        @test length(rr) == 1
        r = first(rr)
        @test typeof.(Tuple(r)) == (
            Missing,
            Int64,
            Int64,
            Float64,
            Float64,
            Int64,
            String,
            String,
            String,
            String,
            String,
            String,
            Base.CodeUnits{UInt8,String},
            String,
            Int64,
        )
    end

    @testset "Issue #158: Missing DB File" begin
        @test_throws SQLiteException SQLite.DB("nonexistentdir/not_there.db")
    end

    @testset "Issue #180, Query" begin
        param = "Hello!"
        query = DBInterface.execute(SQLite.DB(), "SELECT ?1 UNION ALL SELECT ?1", [param])
        param = "x"
        for row in query
            @test row[1] == "Hello!"
            GC.gc() # this must NOT garbage collect the "Hello!" bound value
        end

        db = SQLite.DB()
        DBInterface.execute(db, "CREATE TABLE T (a TEXT, PRIMARY KEY (a))")

        q = DBInterface.prepare(db, "INSERT INTO T VALUES(?)")
        DBInterface.execute(q, ["a"])

        SQLite.bind!(q, 1, "a")
        @test_throws SQLiteException DBInterface.execute(q)
    end

    @testset "Enable extension" begin
        db = SQLite.DB()
        @test SQLite.@OK SQLite.enable_load_extension(db)
    end

    @testset "show(DB)" begin
        io = IOBuffer()
        db = SQLite.DB()

        show(io, db)
        @test String(take!(io)) == "SQLite.DB(\":memory:\")"

        DBInterface.close!(db)
    end

    @testset "SQLite.execute()" begin
        db = SQLite.DB()
        DBInterface.execute(db, "CREATE TABLE T (x INT UNIQUE)")

        q = DBInterface.prepare(db, "INSERT INTO T VALUES(?)")
        SQLite.execute(q, (1,))
        r = DBInterface.execute(db, "SELECT * FROM T") |> columntable
        @test r[1] == [1]

        SQLite.execute(q, [2])
        r = DBInterface.execute(db, "SELECT * FROM T") |> columntable
        @test r[1] == [1, 2]

        q = DBInterface.prepare(db, "INSERT INTO T VALUES(:x)")
        SQLite.execute(q, Dict(:x => 3))
        r = DBInterface.execute(columntable, db, "SELECT * FROM T")
        @test r[1] == [1, 2, 3]

        SQLite.execute(q, x = 4)
        r = DBInterface.execute(columntable, db, "SELECT * FROM T")
        @test r[1] == [1, 2, 3, 4]

        SQLite.execute(db, "INSERT INTO T VALUES(:x)", x = 5)
        r = DBInterface.execute(columntable, db, "SELECT * FROM T")
        @test r[1] == [1, 2, 3, 4, 5]

        r = DBInterface.execute(db, strip("   SELECT * FROM T  ")) |> columntable
        @test r[1] == [1, 2, 3, 4, 5]

        SQLite.createindex!(db, "T", "x", "x_index"; unique = false)
        inds = SQLite.indices(db)
        @test last(inds.name) == "x"
        SQLite.dropindex!(db, "x")
        @test length(SQLite.indices(db).name) == 1

        cols = SQLite.columns(db, "T")
        @test cols.name == ["x"]

        @test SQLite.last_insert_rowid(db) == 5

        r = DBInterface.execute(db, "SELECT * FROM T")
        @test Tables.istable(r)
        @test Tables.rowaccess(r)
        @test Tables.rows(r) === r
        @test Base.IteratorSize(typeof(r)) == Base.SizeUnknown()
        @test eltype(r) == SQLite.Row
        row = first(r)
        SQLite.reset!(r)
        row2 = first(r)
        @test row[:x] == row2[:x]
        @test propertynames(row) == [:x]
        @test DBInterface.lastrowid(r) == 5

        r = DBInterface.execute(db, "SELECT * FROM T") |> columntable
        SQLite.load!(nothing, Tables.rows(r), db, "T2", SQLite.tableinfo(db, "T2"))
        r2 = DBInterface.execute(db, "SELECT * FROM T2") |> columntable
        @test r == r2
    end

    @testset "Escaping" begin
        @test SQLite.esc_id(["1", "2", "3"]) == "\"1\",\"2\",\"3\""
    end

    @testset "Issue #193: Throw informative error on duplicate column names" begin
        db = SQLite.DB()
        @test_throws SQLiteException SQLite.load!((a = [1, 2, 3], A = [1, 2, 3]), db)
    end

    @testset "Issue #216: Table should map by name" begin
        db = SQLite.DB()

        tbl1 = (a = [1, 2, 3], b = [4, 5, 6])
        tbl2 = (b = [7, 8, 9], a = [4, 5, 6])
        SQLite.load!(tbl1, db, "data")
        SQLite.load!(tbl2, db, "data")

        res = DBInterface.execute(db, "SELECT * FROM data") |> columntable
        expected = (a = [1, 2, 3, 4, 5, 6], b = [4, 5, 6, 7, 8, 9])
        @test res == expected
    end

    @testset "Issue #216: Table should error if names don't match" begin
        db = SQLite.DB()

        tbl1 = (a = [1, 2, 3], b = [4, 5, 6])
        SQLite.load!(tbl1, db, "data")
        tbl3 = (c = [7, 8, 9], a = [4, 5, 6])
        @test_throws SQLiteException SQLite.load!(tbl3, db, "data")
    end

    @testset "Test busy_timeout" begin
        db = SQLite.DB()
        @test SQLite.busy_timeout(db, 300) == 0
    end

    @testset "Issue #253: Ensure query column names are unique by default" begin
        db = SQLite.DB()
        res =
            DBInterface.execute(db, "select 1 as x2, 2 as x2, 3 as x2, 4 as x2_2") |>
            columntable
        @test res == (x2 = [1], x2_1 = [2], x2_2 = [3], x2_2_1 = [4])
    end

    @testset "load!() / drop!() table name escaping" begin
        db = SQLite.DB()
        tbl = (a = [1, 2, 3], b = ["a", "b", "c"])
        SQLite.load!(tbl, db, "escape 10.0%")
        r =
            DBInterface.execute(db, "SELECT * FROM $(SQLite.esc_id("escape 10.0%"))") |>
            columntable
        @test r == tbl
        SQLite.drop!(db, "escape 10.0%")
    end

    @testset "load!() column names escaping" begin
        db = SQLite.DB()
        tbl = NamedTuple{(:a, Symbol("50.0%"))}(([1, 2, 3], ["a", "b", "c"]))
        SQLite.load!(tbl, db, "escape_colnames")
        r = DBInterface.execute(db, "SELECT * FROM escape_colnames") |> columntable
        @test r == tbl
        SQLite.drop!(db, "escape_colnames")
    end

    @testset "Bool column data" begin
        db = SQLite.DB()
        tbl = (a = [true, false, false], b = [false, missing, true])
        SQLite.load!(tbl, db, "bool_data")
        r = DBInterface.execute(db, "SELECT * FROM bool_data") |> columntable
        @test isequal(r, (a = [1, 0, 0], b = [0, missing, 1]))
        SQLite.drop!(db, "bool_data")
    end

    @testset "serialization edgecases" begin
        db = SQLite.DB()
        r = DBInterface.execute(db, "SELECT zeroblob(2) as b") |> columntable
        @test first(r.b) == [0, 0]
        r = DBInterface.execute(db, "SELECT zeroblob(0) as b") |> columntable
        @test first(r.b) == []
    end

    @testset "Stmt scope" begin
        dbfile = joinpath(tempdir(), "test_stmt_scope.sqlite")
        db = SQLite.DB(dbfile)
        tbl = (a = [1, 2, 3], b = ["a", "b", "c"])

        @testset "explicit finalization by finalize_statements!(db)" begin
            SQLite.load!(tbl, db, "test_table")
            stmt = SQLite.Stmt(db, "SELECT a, b FROM test_table")
            @test SQLite.isready(stmt)
            @test SQLite.execute(stmt) == 100
            # test cannot drop the table locked by the statement
            @test_throws SQLiteException SQLite.drop!(db, "test_table")
            SQLite.finalize_statements!(db)
            @test !SQLite.isready(stmt)
            SQLite.drop!(db, "test_table")
            DBInterface.close!(stmt) # test can call close!() 2nd time
        end

        @testset "explicit finalization by close!(stmt)" begin
            SQLite.load!(tbl, db, "test_table2")
            stmt = SQLite.Stmt(db, "SELECT a, b FROM test_table2")
            @test SQLite.isready(stmt)
            @test SQLite.execute(stmt) == 100
            # test cannot drop the table locked by the statement
            @test_throws SQLiteException SQLite.drop!(db, "test_table2")
            DBInterface.close!(stmt)
            @test !SQLite.isready(stmt)
            SQLite.drop!(db, "test_table2")
            DBInterface.close!(stmt) # test can call close!() 2nd time
        end

        @testset "automatic close of implicit prepared statement" begin
            @testset "SQLite.execute() call" begin
                SQLite.load!(tbl, db, "test_table3")
                @test SQLite.execute(db, "SELECT a, b FROM test_table3") == 100
                # test can immediately drop the table, since no locks by the statement
                SQLite.drop!(db, "test_table3")
            end

            @testset "DBInterface.execute() call" begin
                SQLite.load!(tbl, db, "test_table4")
                @test SQLite.execute(db, "SELECT a, b FROM test_table4") == 100
                GC.gc() # close implicitly created statement
                # test can immediately drop the table, since no locks by the GC-ed statement
                SQLite.drop!(db, "test_table4")
            end
        end

        close(db)
        rm(dbfile)
    end

end # @testset

struct UnknownSchemaTable end

Tables.isrowtable(::Type{UnknownSchemaTable}) = true
Tables.rows(x::UnknownSchemaTable) = x
Base.length(x::UnknownSchemaTable) = 3
Base.iterate(::UnknownSchemaTable, st = 1) =
    st == 4 ? nothing : ((a = 1, b = 2 + st, c = 3 + st), st + 1)

@testset "misc" begin

    # https://github.com/JuliaDatabases/SQLite.jl/issues/259
    db = SQLite.DB()
    SQLite.load!(UnknownSchemaTable(), db, "tbl")
    tbl = DBInterface.execute(db, "select * from tbl") |> columntable
    @test tbl == (a = [1, 1, 1], b = [3, 4, 5], c = [4, 5, 6])

    # https://github.com/JuliaDatabases/SQLite.jl/issues/251
    q = DBInterface.execute(db, "select * from tbl")
    row, st = iterate(q)
    @test row.a == 1 && row.b == 3 && row.c == 4
    row2, st = iterate(q, st)
    @test_throws ArgumentError row.a

    # https://github.com/JuliaDatabases/SQLite.jl/issues/243
    db = SQLite.DB()
    DBInterface.execute(
        db,
        "create table tmp ( a INTEGER NOT NULL PRIMARY KEY, b INTEGER, c INTEGER )",
    )
    @test_throws SQLite.SQLiteException SQLite.load!(UnknownSchemaTable(), db, "tmp")
    SQLite.load!(UnknownSchemaTable(), db, "tmp"; replace = true)
    tbl = DBInterface.execute(db, "select * from tmp") |> columntable
    @test tbl == (a = [1], b = [5], c = [6])

    db = SQLite.DB()
    DBInterface.execute(db, "create table tmp ( x TEXT )")
    DBInterface.execute(db, "insert into tmp values (?)", (nothing,))
    DBInterface.execute(db, "insert into tmp values (?)", (:a,))
    tbl = DBInterface.execute(db, "select x from tmp") |> columntable
    @test isequal(tbl.x, [missing, :a])

    db = SQLite.DB()
    DBInterface.execute(db, "create table tmp (a integer, b integer, c integer)")
    stmt = DBInterface.prepare(db, "INSERT INTO tmp VALUES(?, ?, ?)")
    tbl = (a = [1, 1, 1], b = [3, 4, 5], c = [4, 5, 6])
    DBInterface.executemany(stmt, tbl)
    tbl2 = DBInterface.execute(db, "select * from tmp") |> columntable
    @test tbl == tbl2

end


text:
steps = AbstractStep[]

f() = nothing
push!(steps, PerPulseStep(f, [:pulse, :pre_samples, :frame_time],
	[:pretrig_mean, :pretrig_rms, :pulse_average, :pulse_rms, :rise_time, :postpeak_deriv, :peak_index, :peak_value, :min_value]))
push!(steps, PerPulseStep(f, [:pretrig_rms, :pretrig_rms_criteria, :peak_index, :peak_index_criteria, :postpeak_deriv, :postpeak_deriv_criteria], [:selection_good]))
push!(steps, PerPulseStep(f, [:pulse, :whitenoise_filter], [:filt_value]))
push!(steps, ThresholdStep(f, [:pretrig_mean, :filt_value, :selection_good], [:ptm_dc],:filt_value_hist, counted, 3000, true))

g = Graphs.graph(steps)


text:
module SwitchOnSafety

using LinearAlgebra, SparseArrays, Statistics

using DynamicPolynomials
using MultivariatePolynomials
using SemialgebraicSets
using MultivariateMoments

import Reexport
Reexport.@reexport using SetProg
Reexport.@reexport using MathematicalSystems
Reexport.@reexport using HybridSystems

export ρ, quicklb, quickub, quickb

function ρ(A::Matrix)
    if isempty(A)
        throw(ArgumentError("Cannot compute the spectral radius of 0x0 matrix"))
    end
    maximum(abs.(eigvals(A)))
end
# eigvals is not defined for SparseMatrixCSC
# eigvals is not defined for SMatrix in StaticArrays for non-Hermitian
ρ(A::AbstractMatrix) = ρ(Matrix(A))

# TODO move to StaticArrays
using StaticArrays
function ρ(A::SMatrix{2, 2, <:Real})
    trace = tr(A)
    deter = det(A)
    # Need to solve λ^2 - trace * λ + deter
    Δ = trace^2 - 4deter
    if Δ >= 0
        sΔ = √Δ
        return max(abs(trace + sΔ), abs(trace - sΔ)) / 2
    else
        return √(trace^2 + abs(Δ)) / 2
    end
end

#abstract type AbstractSwitchedSystem end

include("scaled.jl")

const AbstractDiscreteSwitchedSystem = Union{DiscreteSwitchedLinearSystem, ConstrainedDiscreteSwitchedLinearSystem}
const AbstractSwitchedSystem = AbstractDiscreteSwitchedSystem
integratorfor(s::AbstractDiscreteSwitchedSystem, t) = dynamicfort(s, t)
#integratorfor(s::AbstractContinuousSwitchedSystem, mode::Tuple{Int,Float64}) = expm(dynamicfor(s, mode[1]) * mode[2])

ρA(s::ConstrainedDiscreteSwitchedLinearSystem) = ρ(adjacency_matrix(s.automaton.G))

state(s::AbstractSwitchedSystem, t, forward=true) = forward ? target(s, t) : source(s, t)

nlabels(s::AbstractSwitchedSystem, t) = 1

#modes(s::AbstractSwitchedSystem, v::Int, forward=true) = 1:length(s.A)

function measurefor(μ, s::HybridSystems.DiscreteSwitchingSequence)
    μ[first(s.seq)]
end

# Only makes sense for discrete
function dynamicfort(s::AbstractDiscreteSwitchedSystem, sw::HybridSystems.DiscreteSwitchingSequence)
    sw.A
end

dynamicforσ(s::AbstractDiscreteSwitchedSystem, σ) = s.resetmaps[σ].A
dynamicfort(s::AbstractDiscreteSwitchedSystem, t) = dynamicforσ(s, symbol(s, t))
dynamicfort(s::AbstractDiscreteSwitchedSystem, t, γ) = dynamicfort(s, t) / γ

io_transitions(s, st, forward::Bool) = forward ? out_transitions(s, st) : in_transitions(s, st)

function _eyet(s, t)
    # transpose needed for rectangle system
    n, m = size(integratorfor(s, t)')
    Matrix(1.0I, n, m)
end

include("periodic.jl")
include("quick.jl")
include("gripenberg.jl")

include("polytopes.jl")

include("sosdata.jl")
include("veronese.jl")
include("kronecker.jl")
include("pradius.jl")
include("sos.jl")
include("sosseq.jl")
include("sosext.jl")

include("debruijn.jl")

include("plot.jl")

include("invariant.jl")

end # module