text:
"""
    struct Torus{T,CO,TR,TP,TT,TT1,TT2} <: AbstractVolumePrimitive{T,CO}

Volume primitive describing a [Torus](@ref). It is defined as all points that are
within a given radius to a circle, parallel to the `xy` plane, with constant
radius around a given origin.

## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
* `TR`: Type of `r_tube`.
    * `TR == T`: Full tube without cutout (constant radius `r_tube`).
* `TP`: Type of the azimuthial angle `φ`.
    * `TP == Nothing`: Full 2π in `φ`.
    * `TP == T`: Partial Torus ranging from `0` to `φ`.
* `TT`: Type of the polar angle `θ`.
    * `TT == Nothing`: Full 2π in `θ`.
    * `TP == Tuple{T,T}`: Partial Torus ranging from `θ[1]` to `θ[2]`.
* `TT1`: Type of the surface at `θ[1]` (`:inwards`, `:outwards` or `:flat`).
* `TT2`: Type of the surface at `θ[2]` (`:inwards`, `:outwards` or `:flat`).

## Fields
* `r_torus::T`: Distance of the center of the `Torus` to the center of the tube (in m).
* `r_tube::TR`: Radius of the tube of the `Torus` (in m).
* `φ::TP`: Range in azimuthial angle `φ` of the `Torus`.
* `θ::TT`: Range in polar angle `θ` of the `Torus`.
* `origin::CartesianPoint{T}`: The position of the center of the `Torus`.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Torus` around its `origin`.

## Definition in Configuration File

A `Torus` is defined in the configuration file as part of the `geometry` field
of an object through the field `torus`.

Example definitions of a `Torus` looks like this:
```yaml
torus:
  r_torus: 10.0   # => r_torus = 10.0
  r_tube: 2       # => r_tube = 2.0
  phi:
      from: 0.0°
      to: 360.0°  # => φ = nothing
  theta:
      from: 0.0°
      to: 360.0°  # => θ = nothing
```
The fields `phi` and `theta` do not need to defined if they are full 2π.

Example definition of a `Torus` with an inner cut-out:
```yaml
torus:
  r_torus: 10.0   # => r_torus = 10.0
  r_tube:
      from: 1.0
      to: 2.0     # => r_tube = (1.0, 2.0)
  phi:
      from: 0.0°
      to: 360.0°  # => φ = nothing
  theta:
      from: 0.0°
      to: 360.0°  # => θ = nothing
```
This is a `Torus` with `r_tube` having an inner radius of 1 and an outer radius of 2.

See also [Constructive Solid Geometry (CSG)](@ref).
"""
@with_kw struct Torus{T,CO,TR,TP<:Union{Nothing,T},TT,TT1,TT2} <: AbstractVolumePrimitive{T,CO}
    r_torus::T = 1
    r_tube::TR = 1  # (r_tube_in, r_tube_out)
    φ::TP = nothing
    θ::TT = nothing

    origin::CartesianPoint{T} = zero(CartesianPoint{T})
    rotation::SMatrix{3,3,T,9} = one(SMatrix{3, 3, T, 9})
end

Torus{T,CO,TR,TP,TT,TT1,TT2}( t::Torus{T,CO,TR,TP,TT}; COT = CO,
            origin::CartesianPoint{T} = t.origin,
            rotation::SMatrix{3,3,T,9} = t.rotation) where {T,CO<:Union{ClosedPrimitive, OpenPrimitive},TR,TP<:Union{Nothing,T},TT,TT1,TT2} =
    Torus{T,COT,TR,TP,TT,TT1,TT2}(t.r_torus, t.r_tube, t.φ, t.θ, origin, rotation)

const FullTorus{T,CO} = Torus{T,CO,T,Nothing,Nothing,Nothing,Nothing}
const FullPhiTorus{T,CO,TT1,TT2} = Torus{T,CO,T,Nothing,Tuple{T,T},TT1,TT2}
const FullThetaTorus{T,CO} = Torus{T,CO,T,T,Nothing,Nothing,Nothing}
const HollowTorus{T,CO} = Torus{T,CO,Tuple{T,T},Nothing,Nothing,Nothing,Nothing}
const HollowPhiTorus{T,CO,TT1,TT2} = Torus{T,CO,Tuple{T,T},Nothing,Tuple{T,T},TT1,TT2}
const HollowThetaTorus{T,CO} = Torus{T,CO,Tuple{T,T},T,Nothing,Nothing,Nothing}

function _get_conemantle_type(θ::Tuple{T,T})::Tuple{Symbol, Symbol} where {T}
    θ1::T, θ2::T = θ
    return ( mod(θ1, π) == 0 ? :flat : (mod(θ1, 2π) in 0..π ? :inwards : :outwards),
             mod(θ2, π) == 0 ? :flat : (mod(θ2, 2π) in 0..π ? :outwards : :inwards))
end

_get_conemantle_type(::Nothing) = (Nothing, Nothing)

function Geometry(::Type{T}, ::Type{Torus}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
    length_unit = input_units.length
    angle_unit = input_units.angle
    origin = get_origin(T, dict, length_unit)
    rotation = get_rotation(T, dict, angle_unit)

    r_torus = _parse_value(T, dict["r_torus"], length_unit)
    r_tube = _parse_radial_interval(T, dict["r_tube"], length_unit)
    r_tube isa Tuple{T,T} && iszero(r_tube[1]) ? r_tube = r_tube[2] : nothing
    φ_interval = parse_φ_of_primitive(T, dict, angle_unit)
    φ = if φ_interval isa Tuple{T,T}
        rotation = rotation * RotZ{T}(φ_interval[1])
        φ_interval[2] - φ_interval[1]
    elseif isnothing(φ_interval)
        nothing
    else
        throw(ConfigFileError("Error when trying to parse φ from configuration file."))
    end
    θ = parse_θ_of_primitive(T, dict, angle_unit)
    TT1, TT2 = _get_conemantle_type(θ)
    if haskey(dict, "z")
        @warn "Deprecation warning: Field `z` for `Torus` is deprecated.
                Use instead (only) `origin` to specify the origin of the primitive.
                There might be a conflict with the possible field `origin`:
                The `z` component of the origin of the primitive is overwritten by the `z`."
        z = _parse_value(T, dict["z"], length_unit)
        origin = CartesianPoint{T}(origin[1], origin[2], z)
    end

    t = Torus{T,ClosedPrimitive,typeof(r_tube),typeof(φ),typeof(θ),TT1,TT2}(
            r_torus = r_torus, r_tube = r_tube, φ = φ, θ = θ,
            origin = origin, rotation = rotation)
    transform(t, transformations)
end

function Dictionary(t::Torus{T,<:Any,TR})::OrderedDict{String, Any} where {T,TR}
    dict = OrderedDict{String, Any}()
    dict["r_torus"] = t.r_torus
    dict["r_tube"] = TR <: Real ? t.r_tube : OrderedDict{String, Any}("from" => t.r_tube[1], "to" => t.r_tube[2])
    if !isnothing(t.φ) dict["phi"]   = OrderedDict("from" => "0°", "to" => string(rad2deg(t.φ))*"°") end
    if !isnothing(t.θ) dict["theta"] = OrderedDict("from" => string(rad2deg(t.θ[1]))*"°", "to" => string(rad2deg(t.θ[2]))*"°") end
    if t.origin != zero(CartesianVector{T}) dict["origin"] = t.origin end
    if t.rotation != one(SMatrix{3,3,T,9})
        d = Dictionary(t.rotation)
        if unique(keys(d)) == ["Z"]
            φ0 = mod2pi(_parse_value(T, d["Z"], internal_angle_unit))
            dict["phi"] = OrderedDict("from" => string(rad2deg(φ0))*"°", "to" => string(rad2deg(φ0 + t.φ))*"°")
        else
            dict["rotation"] = d
        end
    end
    OrderedDict{String, Any}("torus" => dict)
end

function surfaces(t::FullTorus{T,ClosedPrimitive}) where {T}
    tm = FullTorusMantle{T,:inwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    (tm, )
end
function surfaces(t::FullTorus{T,OpenPrimitive}) where {T}
    tm = FullTorusMantle{T,:outwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    (tm, )
end
function surfaces(t::FullPhiTorus{T,ClosedPrimitive,TT1,TT2}) where {T,TT1,TT2}
    tm = FullPhiTorusMantle{T,:inwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ1)), t.φ, t.r_tube * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ1)/2), t.rotation, Val{TT1}())
    cm2 = TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ2)), t.φ, t.r_tube * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ2)/2), t.rotation, Val{TT2}())
    (tm, cm1, cm2)
end
function surfaces(t::FullPhiTorus{T,OpenPrimitive,TT1,TT2}) where {T,TT1,TT2}
    tm = FullPhiTorusMantle{T,:outwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = flip(TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ1)), t.φ, t.r_tube * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ1)/2), t.rotation, Val{TT1}()))
    cm2 = flip(TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ2)), t.φ, t.r_tube * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ2)/2), t.rotation, Val{TT2}()))
    (tm, cm1, cm2)
end
function surfaces(t::FullThetaTorus{T,ClosedPrimitive}) where {T}
    tm = FullThetaTorusMantle{T,:inwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    es1 = CircularArea{T}(t.r_tube, t.θ, t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(-π/2) )
    es2 = CircularArea{T}(t.r_tube, t.θ, t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(π/2) )
    (tm, es1, es2)
end
function surfaces(t::FullThetaTorus{T,OpenPrimitive}) where {T}
    tm = FullThetaTorusMantle{T,:outwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    φ1, φ2 = T(0), t.φ
    es1 = CircularArea{T}(t.r_tube, t.θ, t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(π/2) )
    es2 = CircularArea{T}(t.r_tube, t.θ, t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(-π/2) )
    (tm, es1, es2)
end
function surfaces(t::Torus{T,ClosedPrimitive,T,T,Tuple{T,T},TT1,TT2}) where {T,TT1,TT2}
    tm = TorusMantle{T,T,Tuple{T,T},:inwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ1)), t.φ, t.r_tube * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ1)/2), t.rotation, Val{TT1}())
    cm2 = TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ2)), t.φ, t.r_tube * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ2)/2), t.rotation, Val{TT2}())
    es1 = flip(PartialCircularArea{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(π/2) * RotZ{T}(t.θ[1])))
    es2 = PartialCircularArea{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(π/2) * RotZ{T}(t.θ[1]))
    (tm, cm1, cm2, es1, es2)
end
function surfaces(t::Torus{T,OpenPrimitive,T,T,Tuple{T,T},TT1,TT2}) where {T,TT1,TT2}
    tm = TorusMantle{T,T,Tuple{T,T},:outwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = flip(TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ1)), t.φ, t.r_tube * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ1)/2), t.rotation, Val{TT1}()))
    cm2 = flip(TorusThetaSurface((t.r_torus, t.r_torus + t.r_tube * cos(θ2)), t.φ, t.r_tube * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,t.r_tube * sin(θ2)/2), t.rotation, Val{TT2}()))
    es1 = PartialCircularArea{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(π/2) * RotZ{T}(t.θ[1]))
    es2 = flip(PartialCircularArea{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(π/2) * RotZ{T}(t.θ[1])))
    (tm, cm1, cm2, es1, es2)
end

function surfaces(t::HollowTorus{T,ClosedPrimitive}) where {T}
    tm_in = FullTorusMantle{T,:outwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = FullTorusMantle{T,:inwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    (tm_in, tm_out)
end
function surfaces(t::HollowTorus{T,OpenPrimitive}) where {T}
    tm_in = FullTorusMantle{T,:inwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = FullTorusMantle{T,:outwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    (tm_in, tm_out)
end
function surfaces(t::HollowPhiTorus{T,ClosedPrimitive,TT1,TT2}) where {T,TT1,TT2}
    tm_in = FullPhiTorusMantle{T,:outwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = FullPhiTorusMantle{T,:inwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ1), t.φ, abs(-(t.r_tube...)) * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ1)), t.rotation, Val{TT1}())
    cm2 = TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ2), t.φ, abs(-(t.r_tube...)) * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ2)), t.rotation, Val{TT2}())
    (tm_in, tm_out, cm1, cm2)
end
function surfaces(t::HollowPhiTorus{T,OpenPrimitive,TT1,TT2}) where {T,TT1,TT2}
    tm_in = FullPhiTorusMantle{T,:inwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = FullPhiTorusMantle{T,:outwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = flip(TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ1), t.φ, abs(-(t.r_tube...)) * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ1)), t.rotation, Val{TT1}()))
    cm2 = flip(TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ2), t.φ, abs(-(t.r_tube...)) * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ2)), t.rotation, Val{TT2}()))
    (tm_in, tm_out, cm1, cm2)
end
function surfaces(t::HollowThetaTorus{T,ClosedPrimitive}) where {T}
    tm_in = FullThetaTorusMantle{T,:outwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = FullThetaTorusMantle{T,:inwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    es1 = Annulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(-π/2)  * RotZ{T}(t.θ[1]))
    es2 = Annulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(π/2)  * RotZ{T}(t.θ[1]))
    (tm_in, tm_out, es1, es2)
end
function surfaces(t::HollowThetaTorus{T,OpenPrimitive}) where {T}
    tm_in = FullThetaTorusMantle{T,:inwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = FullThetaTorusMantle{T,:outwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    es1 = Annulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(π/2)  * RotZ{T}(t.θ[1]))
    es2 = Annulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(-π/2)  * RotZ{T}(t.θ[1]))
    (tm_in, tm_out, es1, es2)
end
function surfaces(t::Torus{T,ClosedPrimitive,Tuple{T,T},T,Tuple{T,T},TT1,TT2}) where {T,TT1,TT2}
    tm_in = TorusMantle{T,T,Tuple{T,T},:outwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = TorusMantle{T,T,Tuple{T,T},:inwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ1), t.φ, abs(-(t.r_tube...)) * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ1)), t.rotation, Val{TT1}())
    cm2 = TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ2), t.φ, abs(-(t.r_tube...)) * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ2)), t.rotation, Val{TT2}())
    es1 = flip(PartialAnnulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(π/2) * RotZ{T}(t.θ[1])))
    es2 = PartialAnnulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(π/2) * RotZ{T}(t.θ[1]))
    (tm_in, tm_out, cm1, cm2, es1, es2)
end
function surfaces(t::Torus{T,OpenPrimitive,Tuple{T,T},T,Tuple{T,T},TT1,TT2}) where {T,TT1,TT2}
    tm_in = TorusMantle{T,T,Tuple{T,T},:inwards}(t.r_torus, t.r_tube[1], t.φ, t.θ, t.origin, t.rotation)
    tm_out = TorusMantle{T,T,Tuple{T,T},:outwards}(t.r_torus, t.r_tube[2], t.φ, t.θ, t.origin, t.rotation)
    θ1, θ2 = t.θ
    cm1 = flip(TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ1), t.φ, abs(-(t.r_tube...)) * sin(θ1)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ1)), t.rotation, Val{TT1}()))
    cm2 = flip(TorusThetaSurface(t.r_torus .+ t.r_tube .* cos(θ2), t.φ, abs(-(t.r_tube...)) * sin(θ2)/2, t.origin + t.rotation * CartesianVector{T}(0,0,mean(t.r_tube) * sin(θ2)), t.rotation, Val{TT2}()))
    es1 = PartialAnnulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus, 0, 0), t.rotation * RotX(π/2) * RotZ{T}(t.θ[1]))
    es2 = flip(PartialAnnulus{T}(t.r_tube, abs(-(t.θ...)), t.origin + t.rotation * CartesianVector{T}(t.r_torus * cos(t.φ), t.r_torus * sin(t.φ), 0), t.rotation * RotZ(t.φ) * RotX(π/2) * RotZ{T}(t.θ[1])))
    (tm_in, tm_out, cm1, cm2, es1, es2)
end

function _in(pt::CartesianPoint{T}, t::Torus{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    _r = hypot(hypot(pt.x, pt.y) - t.r_torus, pt.z)
    rmin::T, rmax::T = _radial_endpoints(t.r_tube)
    return rmin - csgtol <= _r <= rmax + csgtol &&
        (isnothing(t.φ) || _in_angular_interval_closed(atan(pt.y, pt.x), t.φ, csgtol = csgtol)) &&
        (isnothing(t.θ) || _in_angular_interval_closed(atan(pt.z, hypot(pt.x, pt.y) - t.r_torus), t.θ, csgtol = csgtol))
end
function _in(pt::CartesianPoint{T}, t::Torus{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    _r = hypot(hypot(pt.x, pt.y) - t.r_torus, pt.z)
    rmin::T, rmax::T = _radial_endpoints(t.r_tube)
    return rmin + csgtol < _r < rmax - csgtol &&
        (isnothing(t.φ) || _in_angular_interval_open(atan(pt.y, pt.x), t.φ, csgtol = csgtol)) &&
        (isnothing(t.θ) || _in_angular_interval_open(atan(pt.z, hypot(pt.x, pt.y) - t.r_torus), t.θ, csgtol = csgtol))
end

# #Constructors
# function Torus(;r_torus = 1, r_tubeMin = 0, r_tubeMax = 1, φMin = 0, φMax = 2π, θMin = 0, θMax = 2π, z = 0)
#     T = float(promote_type(typeof.((r_torus, r_tubeMin, r_tubeMax, φMin, φMax, θMin, θMax, z))...))
#     r_tube = r_tubeMin == 0 ? T(r_tubeMax) : T(r_tubeMin)..T(r_tubeMax)
#     φ = mod(T(φMax) - T(φMin), T(2π)) == 0 ? nothing : T(φMin)..T(φMax)
#     θ = mod(T(θMax) - T(θMin), T(2π)) == 0 ? nothing : T(θMin)..T(θMax)
#     Torus( T, T(r_torus), r_tube, φ, θ, T(z))
# end

# Torus(r_torus, r_tubeMin, r_tubeMax, φMin, φMax, θMin, θMax, z) = Torus(;r_torus = r_torus, r_tubeMin = r_tubeMin, r_tubeMax = r_tubeMax, φMin = φMin, φMax = φMax, θMin = θMin, θMax = θMax, z = z)

# function Torus(r_torus::R1, r_tube::R2, z::TZ) where {R1<:Real, R2<:Real, TZ<:Real}
#     T = float(promote_type(R1, R2, TZ))
#     Torus( T, T(r_torus), T(r_tube), nothing, nothing, T(z))
# end

# function RoundChamfer(r_torus::R1, r_tube::R2, z::TZ) where {R1<:Real, R2<:Real, TZ<:Real}
#     T = float(promote_type(R1, R2, TZ))
#     Torus( T, T(r_torus), T(r_tube), nothing, T(0)..T(π/2), T(z))
# end

extremum(t::Torus{T}) where {T} = t.r_torus + max(t.r_tube...)

# function Dictionary(t::Torus{T}) where {T}
#     dict = OrderedDict{String,Any}()
#     dict["r_torus"] = t.r_torus
#     dict["r_tube"] = typeof(t.r_tube) == T ? t.r_tube : OrderedDict{String,Any}("from" => t.r_tube.left, "to" => t.r_tube.right)
#     if !isnothing(t.φ) dict["phi"] = OrderedDict{String,Any}("from" => t.φ.left, "to" => t.φ.right) end
#     if !isnothing(t.θ) dict["theta"] = OrderedDict{String,Any}("from" => t.θ.left, "to" => t.θ.right) end
#     if t.z != 0 dict["z"] = t.z end
#     OrderedDict{String,Any}("torus" => dict)
# end

# get_r_tube_limits(t::Torus{T}) where {T} = _radial_endpoints(t.r_tube)
# get_φ_limits(t::Torus{T, <:Any, Nothing, <:Any}) where {T} = (T(0), T(2π), true)
# get_φ_limits(t::Torus{T, <:Any, <:AbstractInterval, <:Any}) where {T} = (t.φ.left, t.φ.right, false)

# get_θ_limits(t::Torus{T, <:Any, <:Any, Nothing}) where {T} = (T(0), T(2π), true)
# get_θ_limits(t::Torus{T, <:Any, <:Any, <:AbstractInterval}) where {T} = (t.θ.left, t.θ.right, false)

# function _is_torus_collapsed(t::Torus{T}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     isapprox(r_tubeMin, r_tubeMax, atol = geom_atol_zero(T))
# end

# function _get_decomposed_surfaces_torus(t::Torus{T, T}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     AbstractSurfacePrimitive[TorusMantle(t, r_tube = r_tubeMax)]
# end

# function _get_decomposed_surfaces_torus(t::Torus{T, <:AbstractInterval{T}}) where {T}
#     surfaces = AbstractSurfacePrimitive[]
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     if isapprox(r_tubeMin, r_tubeMax, atol = geom_atol_zero(T))
#         push!(surfaces, TorusMantle(t, r_tube = r_tubeMax))
#     else
#         push!(surfaces, TorusMantle(t, r_tube = r_tubeMin), TorusMantle(t, r_tube = r_tubeMax))
#     end
#     surfaces
# end

# get_decomposed_surfaces(t::Torus{T, <:Any, Nothing, Nothing}) where {T} = _get_decomposed_surfaces_torus(t)

# function get_decomposed_surfaces(t::Torus{T, <:Any, <:AbstractInterval, Nothing}) where {T}
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     surfaces = _get_decomposed_surfaces_torus(t)
#     if !_is_torus_collapsed(t)
#         push!(surfaces, ToroidalAnnulus(t, φ = φMin), ToroidalAnnulus(t, φ = φMax))
#     end
#     surfaces
# end

# function get_decomposed_surfaces(t::Torus{T, <:Any, Nothing, <:AbstractInterval}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     surfaces = _get_decomposed_surfaces_torus(t)
#     if !_is_torus_collapsed(t)
#         if r_tubeMin == T(0) && minmax(mod.((θMin, θMax), T(2π))...) == (T(0),T(π))
#             rMin = t.r_torus - r_tubeMax
#             rMax = t.r_torus + r_tubeMax
#             r = rMin == T(0) ? T(rMax) : T(rMin)..T(rMax)
#             push!(surfaces, CylindricalAnnulus(T, r, t.φ, t.z))
#         else
#             for θ in [θMin, θMax]
#                 mod(θ, T(2π)) in [T(0),T(π)] ? push!(surfaces, CylindricalAnnulus(t, θ = θ)) : push!(surfaces, ConeMantle(t, θ = θ))
#             end
#         end
#     end
#     surfaces
# end

# function get_decomposed_surfaces(t::Torus{T, <:Any, <:AbstractInterval, <:AbstractInterval}) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     surfaces = _get_decomposed_surfaces_torus(t)
#     if !_is_torus_collapsed(t)
#         if r_tubeMin == T(0) && minmax(mod.((θMin, θMax), T(2π))...) == (T(0),T(π))
#             rMin = t.r_torus - r_tubeMax
#             rMax = t.r_torus + r_tubeMax
#             r = rMin == T(0) ? T(rMax) : T(rMin)..T(rMax)
#             push!(surfaces, CylindricalAnnulus(T, r, t.φ, t.z))
#         else
#             for θ in [θMin, θMax]
#                 mod(θ, T(2π)) in [T(0),T(π)] ? push!(surfaces, CylindricalAnnulus(t, θ = θ)) : push!(surfaces, ConeMantle(t, θ = θ))
#             end
#         end
#         push!(surfaces, ToroidalAnnulus(t, φ = φMin), ToroidalAnnulus(t, φ = φMax))
#     end
#     surfaces
# end

# function sample(t::Torus{T}, step::Real) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     samples = [
#         CylindricalPoint{T}(t.r_torus+r_tube*cos(θ),φ,r_tube*sin(θ))
#         for r_tube in r_tubeMin:step:r_tubeMax
#         for θ in (r_tube == 0 ? θMin : θMin:step/r_tube:θMax)
#         for φ in (r_tube == 0 ? φMin : φMin:step/r_tube:φMax)
#     ]
# end

# function sample(t::Torus{T}, Nsamps::NTuple{3,Int} = (2,5,3)) where {T}
#     r_tubeMin::T, r_tubeMax::T = get_r_tube_limits(t)
#     θMin::T, θMax::T, _ = get_θ_limits(t)
#     φMin::T, φMax::T, _ = get_φ_limits(t)
#     samples = [
#         CylindricalPoint{T}(t.r_torus+r_tube*cos(θ),φ,r_tube*sin(θ))
#         for r_tube in (Nsamps[1] ≤ 1 ? r_tubeMin : range(r_tubeMin, r_tubeMax, length = Nsamps[1]))
#         for θ in (Nsamps[3] ≤ 1 ? θMin : range(θMin, θMax, length = Nsamps[3]))
#         for φ in (Nsamps[2] ≤ 1 ? φMin : range(φMin, φMax, length = Nsamps[2]))
#     ]
# end

ext: jl


text:
default_metric_dispatch(::Stiefel, ::EuclideanMetric) = Val(true)

function decorator_transparent_dispatch(::typeof(distance), ::Stiefel, args...)
    return Val(:intransparent)
end

@doc raw"""
    exp(M::Stiefel, p, X)

Compute the exponential map on the [`Stiefel`](@ref)`{n,k,𝔽}`() manifold `M`
emanating from `p` in tangent direction `X`.

````math
\exp_p X = \begin{pmatrix}
   p\\X
 \end{pmatrix}
 \operatorname{Exp}
 \left(
 \begin{pmatrix} p^{\mathrm{H}}X & - X^{\mathrm{H}}X\\
 I_n & p^{\mathrm{H}}X\end{pmatrix}
 \right)
\begin{pmatrix}  \exp( -p^{\mathrm{H}}X) \\ 0_n\end{pmatrix},
````

where $\operatorname{Exp}$ denotes matrix exponential,
$\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian, and $I_k$ and
$0_k$ are the identity matrix and the zero matrix of dimension $k × k$, respectively.
"""
exp(::Stiefel, ::Any...)

function exp!(::Stiefel{n,k}, q, p, X) where {n,k}
    A = p' * X
    B = exp([A -X'*X; I A])
    @views begin
        r = p * B[1:k, 1:k]
        mul!(r, X, B[(k + 1):(2 * k), 1:k], true, true)
    end
    mul!(q, r, exp(-A))
    return q
end

@doc raw"""
    get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis) where {n,k}

Create the default basis using the parametrization for any $X ∈ T_p\mathcal M$.
Set $p_\bot \in ℝ^{n\times(n-k)}$ the matrix such that the $n\times n$ matrix of the common
columns $[p\ p_\bot]$ is an ONB.
For any skew symmetric matrix $a ∈ ℝ^{k\times k}$ and any $b ∈ ℝ^{(n-k)\times k}$ the matrix

````math
X = pa + p_\bot b ∈ T_p\mathcal M
````

and we can use the $\frac{1}{2}k(k-1) + (n-k)k = nk-\frac{1}{2}k(k+1)$ entries
of $a$ and $b$ to specify a basis for the tangent space.
using unit vectors for constructing both
the upper matrix of $a$ to build a skew symmetric matrix and the matrix b, the default
basis is constructed.

Since $[p\ p_\bot]$ is an automorphism on $ℝ^{n\times p}$ the elements of $a$ and $b$ are
orthonormal coordinates for the tangent space. To be precise exactly one element in the upper
trangular entries of $a$ is set to $1$ its symmetric entry to $-1$ and we normalize with
the factor $\frac{1}{\sqrt{2}}$ and for $b$ one can just use unit vectors reshaped to a matrix
to obtain orthonormal set of parameters.
"""
function get_basis(M::Stiefel{n,k,ℝ}, p, B::DefaultOrthonormalBasis{ℝ}) where {n,k}
    V = get_vectors(M, p, B)
    return CachedBasis(B, V)
end

function get_coordinates!(
    M::Stiefel{n,k,ℝ},
    c,
    p,
    X,
    B::DefaultOrthonormalBasis{ℝ},
) where {n,k}
    V = get_vectors(M, p, B)
    c .= inner.(Ref(M), Ref(p), V, Ref(X))
    return c
end

function get_vector!(M::Stiefel{n,k,ℝ}, X, p, c, B::DefaultOrthonormalBasis{ℝ}) where {n,k}
    V = get_vectors(M, p, B)
    zero_tangent_vector!(M, X, p)
    length(c) < length(V) && error(
        "Coordinate vector too short. Excpected $(length(V)), but only got $(length(c)) entries.",
    )
    @inbounds for i in 1:length(V)
        X .+= c[i] .* V[i]
    end
    return X
end

function get_vectors(::Stiefel{n,k,ℝ}, p, ::DefaultOrthonormalBasis{ℝ}) where {n,k}
    p⊥ = nullspace([p zeros(n, n - k)])
    an = div(k * (k - 1), 2)
    bn = (n - k) * k
    V = vcat(
        [p * vec2skew(1 / sqrt(2) .* _euclidean_unit_vector(an, i), k) for i in 1:an],
        [p⊥ * reshape(_euclidean_unit_vector(bn, j), (n - k, k)) for j in 1:bn],
    )
    return V
end

_euclidean_unit_vector(n, i) = [k == i ? 1.0 : 0.0 for k in 1:n]

@doc raw"""
    project(M::Stiefel,p)

Projects `p` from the embedding onto the [`Stiefel`](@ref) `M`, i.e. compute `q`
as the polar decomposition of $p$ such that $q^{\mathrm{H}q$ is the identity,
where $\cdot^{\mathrm{H}}$ denotes the hermitian, i.e. complex conjugate transposed.
"""
project(::Stiefel, ::Any, ::Any)

function project!(::Stiefel, q, p)
    s = svd(p)
    mul!(q, s.U, s.Vt)
    return q
end

@doc raw"""
    project(M::Stiefel, p, X)

Project `X` onto the tangent space of `p` to the [`Stiefel`](@ref) manifold `M`.
The formula reads

````math
\operatorname{proj}_{\mathcal M}(p, X) = X - p \operatorname{Sym}(p^{\mathrm{H}}X),
````

where $\operatorname{Sym}(q)$ is the symmetrization of $q$, e.g. by
$\operatorname{Sym}(q) = \frac{q^{\mathrm{H}}+q}{2}$.
"""
project(::Stiefel, ::Any...)

function project!(::Stiefel, Y, p, X)
    A = p' * X
    T = eltype(Y)
    copyto!(Y, X)
    mul!(Y, p, A + A', T(-0.5), true)
    return Y
end

function vector_transport_to!(M::Stiefel, Y, ::Any, X, q, ::ProjectionTransport)
    return project!(M, Y, q, X)
end


text:
export get_weather_darksky

function fetch_weather_darksky(api_key, time, longitude, latitude, altitude)
    timestr = Dates.format(astimezone(time, tz"UTC"), "yyyy-mm-ddTHH:MM:SS")

    url = "https://api.darksky.net/forecast/$( api_key )/$( latitude ),$( longitude ),$( timestr )Z?exclude=currently,alerts&units=si"
    req = HTTP.get(url)

    j = JSON.parse(String(req.body))

    weather = vcat(DataFrame.(merge.(Ref(Dict("precipType" => missing, "precipAccumulation" => missing, "windBearing" => missing, "pressure" => missing)), j["hourly"]["data"]))...)
    weather[:time] = ZonedDateTime.(unix2datetime.(weather[:time]), tz"UTC")
    weather[:pressure_absolute] = pressure_absolute.(weather[:pressure], altitude)
    weather[:density] = air_density.(weather[:pressure_absolute] .* 100, weather[:temperature] .+ 273.15, weather[:humidity])

    weather
end

fetch_weather_darksky(time, longitude, latitude, altitude) = fetch_weather_darksky(ENV["DARKSKY_API_KEY"], time, longitude, latitude, altitude)

function get_weather_darksky(time_date, ampm, longitude, latitude, altitude)
    gametime = create_time(time_date, ampm)

    w = fetch_weather_darksky(time_date, longitude, latitude, altitude)
    w2 = filter(r -> (round(time_date, Hour, RoundDown) + Hour(4)) >= r[Symbol("time")] >= round(time_date, Hour, RoundDown), w)

    wind_bearing = mean(skipmissing(w2[:windBearing]))
    pressure_absolute = mean(skipmissing(w2[:pressure_absolute]))

    (
        density           = round(mean(w2[:density]), digits = 4),
        temperature       = round(mean(w2[:temperature]), digits = 1),
        dew_point         = round(mean(w2[:dewPoint]), digits = 1),
        humidity          = round(mean(w2[:humidity]), digits = 3),
        pressure_absolute = ismissing(pressure_absolute) ? pressure_absolute : round(pressure_absolute, digits = 2),
        wind_bearing      = ismissing(wind_bearing) ? wind_bearing : round(wind_bearing, digits = 1),
        wind_speed        = round(mean(w2[:windSpeed]), digits = 2),
        wind_gust         = round(mean(w2[:windGust]), digits = 2),
    )
end


text:
################
# addition (+) #
################

# in-place version (necessary for nested differentiation to work)

function plus!(out, x, y)
    for i in eachindex(out)
        out[i] = x[i] + y[i]
    end
    return out
end

@inline plus!(out::TrackedArray, x::TrackedArray, y::TrackedArray) = record_plus!(out, x, y)

for A in ARRAY_TYPES
    @eval @inline plus!(out::TrackedArray, x::TrackedArray, y::$(A)) = record_plus!(out, x, y)
    @eval @inline plus!(out::TrackedArray, x::$(A), y::TrackedArray) = record_plus!(out, x, y)
end

function record_plus!(out::TrackedArray, x, y)
    copyto!(value(out), value(x) + value(y))
    record!(tape(x, y), SpecialInstruction, +, (x, y), out)
    return out
end

# Base allocating version

@inline Base.:+(x::TrackedArray{X,D}, y::TrackedArray{Y,D}) where {X,Y,D} = record_plus(x, y, D)

for A in ARRAY_TYPES
    @eval @inline Base.:+(x::TrackedArray{V,D}, y::$(A)) where {V,D} = record_plus(x, y, D)
    @eval @inline Base.:+(x::$(A), y::TrackedArray{V,D}) where {V,D} = record_plus(x, y, D)
end

function record_plus(x, y, ::Type{D}) where D
    tp = tape(x, y)
    out = track(value(x) + value(y), D, tp)
    record!(tp, SpecialInstruction, +, (x, y), out)
    return out
end

# reverse/forward passes

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(+)})
    a, b = instruction.input
    output = instruction.output
    output_deriv = deriv(output)
    istracked(a) && increment_deriv!(a, output_deriv)
    istracked(b) && increment_deriv!(b, output_deriv)
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(+)})
    a, b = instruction.input
    pull_value!(a)
    pull_value!(b)
    plus!(value(instruction.output), value(a), value(b))
    return nothing
end

###################
# subtraction (-) #
###################

# in-place version (necessary for nested differentiation to work)

function minus!(out, x)
    for i in eachindex(out)
        out[i] = -(x[i])
    end
end

function minus!(out, x, y)
    for i in eachindex(out)
        out[i] = x[i] - y[i]
    end
end

@inline minus!(out::TrackedArray, x::TrackedArray, y::TrackedArray) = record_minus!(out, x, y)
@inline minus!(out::TrackedArray, x::TrackedArray) = record_minus!(out, x, y)

for A in ARRAY_TYPES
    @eval @inline minus!(out::TrackedArray, x::TrackedArray, y::$(A)) = record_minus!(out, x, y)
    @eval @inline minus!(out::TrackedArray, x::$(A), y::TrackedArray) = record_minus!(out, x, y)
    @eval @inline minus!(out::TrackedArray, x::$(A)) = record_minus!(out, x)
end

function record_minus!(out::TrackedArray, x)
    copyto!(value(out), -(value(x)))
    record!(tape(x), SpecialInstruction, -, x, out)
    return out
end

function record_minus!(out::TrackedArray, x, y)
    copyto!(value(out), value(x) - value(y))
    record!(tape(x, y), SpecialInstruction, -, (x, y), out)
    return out
end

# Base allocating version

Base.:-(x::TrackedArray{X,D}, y::TrackedArray{Y,D}) where {X,Y,D} = record_minus(x, y, D)

for A in ARRAY_TYPES
    @eval Base.:-(x::TrackedArray{V,D}, y::$(A)) where {V,D} = record_minus(x, y, D)
    @eval Base.:-(x::$(A), y::TrackedArray{V,D}) where {V,D} = record_minus(x, y, D)
end

function Base.:-(x::TrackedArray{V,D}) where {V,D}
    tp = tape(x)
    out = track(-(value(x)), D, tp)
    record!(tp, SpecialInstruction, -, x, out)
    return out
end

function record_minus(x, y, ::Type{D}) where D
    tp = tape(x, y)
    out = track(value(x) - value(y), D, tp)
    record!(tp, SpecialInstruction, -, (x, y), out)
    return out
end

# reverse/forward passes

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(-)})
    input = instruction.input
    output = instruction.output
    output_deriv = deriv(output)
    if istracked(input)
        decrement_deriv!(input, output_deriv)
    else
        a, b = input
        istracked(a) && increment_deriv!(a, output_deriv)
        istracked(b) && decrement_deriv!(b, output_deriv)
    end
    unseed!(output)
    return nothing
end

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(-)})
    input = instruction.input
    output = instruction.output
    if istracked(input)
        minus!(value(output), value(input))
    else
        a, b = input
        pull_value!(a)
        pull_value!(b)
        minus!(value(output), value(a), value(b))
    end
    return nothing
end

######################
# multiplication (*) #
######################

mulargvalue(x) = value(x)
mulargvalue(x::Adjoint) = adjoint(value(adjoint(x)))
mulargvalue(x::Transpose) = transpose(value(transpose(x)))

mulargpullvalue!(x) = pull_value!(x)
mulargpullvalue!(x::Adjoint) = pull_value!(adjoint(x))
mulargpullvalue!(x::Transpose) = pull_value!(transpose(x))

# recording pass #
#----------------#

@inline function record_mul(x, y, ::Type{D}) where D
    tp = tape(x, y)
    out = track(*(mulargvalue(x), mulargvalue(y)), D, tp)
    cache = (similar(x, D), similar(y, D))
    record!(tp, SpecialInstruction, *, (x, y), out, cache)
    return out
end

@inline function record_mul!(out::TrackedArray{V,D}, x, y) where {V,D}
    copyto!(mulargvalue(out), *(mulargvalue(x), mulargvalue(y)))
    cache = (similar(x, D), similar(y, D))
    record!(tape(x, y), SpecialInstruction, *, (x, y), out, cache)
    return out
end

for S1 in (:TrackedArray, :TrackedVector, :TrackedMatrix)
    for S2 in (:TrackedArray, :TrackedVector, :TrackedMatrix)
        @eval begin
            LinearAlgebra.:*(x::$(S1){X,D}, y::$(S2){Y,D}) where {X,Y,D} = record_mul(x, y, D)

            LinearAlgebra.:*(x::Transpose{<:TrackedReal,<:$(S1){X,D}}, y::Transpose{<:TrackedReal,<:$(S2){Y,D}}) where {X,Y,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::Adjoint{<:TrackedReal,<:$(S1){X,D}}, y::Adjoint{<:TrackedReal,<:$(S2){Y,D}}) where {X,Y,D} = record_mul(x, y, D)

            LinearAlgebra.:*(x::Transpose{<:TrackedReal,<:$(S1){X,D}}, y::$(S2){Y,D}) where {X,Y,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(S1){X,D}, y::Transpose{<:TrackedReal,<:$(S2){Y,D}}) where {X,Y,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::Adjoint{<:TrackedReal,<:$(S1){X,D}}, y::$(S2){Y,D}) where {X,Y,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(S1){X,D}, y::Adjoint{<:TrackedReal,<:$(S2){Y,D}}) where {X,Y,D} = record_mul(x, y, D)

            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::$(S1){X,D}, y::$(S2){Y,D}) where {V,X,Y,D} = record_mul!(out, x, y)

            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::Transpose{<:TrackedReal,<:$(S1){X,D}}, y::Transpose{<:TrackedReal,<:$(S2){Y,D}}) where {V,X,Y,D} = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::Adjoint{<:TrackedReal,<:$(S1){X,D}}, y::Adjoint{<:Number,<:$(S2){Y,D}}) where {V,X,Y,D} = record_mul!(out, x, y)

            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::Transpose{<:TrackedReal,<:$(S1){X,D}}, y::$(S2){Y,D}) where {V,X,Y,D} = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::$(S1){X,D}, y::Transpose{<:TrackedReal,<:$(S2){Y,D}}) where {V,X,Y,D} = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::Adjoint{<:TrackedReal,<:$(S1){X,D}}, y::$(S2){Y,D}) where {V,X,Y,D} = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray{V,D}, x::$(S1){X,D}, y::Adjoint{<:TrackedReal,<:$(S2){Y,D}}) where {V,X,Y,D} = record_mul!(out, x, y)
        end
    end

    for T in ARRAY_TYPES
        @eval begin
            LinearAlgebra.:*(x::$(S1){V,D}, y::$(T)) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(T), y::$(S1){V,D}) where {V,D} = record_mul(x, y, D)

            LinearAlgebra.:*(x::Transpose{<:Number,<:$(T)}, y::$(S1){V,D}) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(S1){V,D}, y::Transpose{<:Number,<:$(T)}) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::Adjoint{<:Number,<:$(T)}, y::$(S1){V,D}) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(S1){V,D}, y::Adjoint{<:Number,<:$(T)}) where {V,D} = record_mul(x, y, D)

            LinearAlgebra.:*(x::Transpose{<:Number,<:$(S1){V,D}}, y::$(T)) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(T), y::Transpose{<:Number,<:$(S1){V,D}}) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::Adjoint{<:Number,<:$(S1){V,D}}, y::$(T)) where {V,D} = record_mul(x, y, D)
            LinearAlgebra.:*(x::$(T), y::Adjoint{<:Number,<:$(S1){V,D}}) where {V,D} = record_mul(x, y, D)

            LinearAlgebra.mul!(out::TrackedArray, x::$(S1), y::$(T)) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::$(T), y::$(S1)) = record_mul!(out, x, y)

            LinearAlgebra.mul!(out::TrackedArray, x::$(S1), y::Transpose{<:Number,<:$(T)}) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::Transpose{<:Number,<:$(T)}, y::$(S1)) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::$(S1), y::Adjoint{<:Number,<:$(T)}) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::Adjoint{<:Number,<:$(T)}, y::$(S1)) = record_mul!(out, x, y)

            LinearAlgebra.mul!(out::TrackedArray, x::Transpose{<:Number,<:$(S1)}, y::$(T)) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::$(T), y::Transpose{<:Number,<:$(S1)}) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::Adjoint{<:Number,<:$(S1)}, y::$(T)) = record_mul!(out, x, y)
            LinearAlgebra.mul!(out::TrackedArray, x::$(T), y::Adjoint{<:Number,<:$(S1)}) = record_mul!(out, x, y)
        end
    end
end

# forward pass #
#--------------#

@noinline function special_forward_exec!(instruction::SpecialInstruction{typeof(*)})
    a, b = instruction.input
    mulargpullvalue!(a)
    mulargpullvalue!(b)
    if instruction.output isa Number
        value!(instruction.output, mulargvalue(a) * mulargvalue(b))
    else
        mul!(mulargvalue(instruction.output), mulargvalue(a), mulargvalue(b))
    end
    return nothing
end

# reverse pass #
#--------------#

@noinline function special_reverse_exec!(instruction::SpecialInstruction{typeof(*)})
    a, b = instruction.input
    a_tmp, b_tmp = instruction.cache
    output = instruction.output
    output_deriv = deriv(output)
    reverse_mul!(output, output_deriv, a, b, a_tmp, b_tmp)
    unseed!(output)
    return nothing
end

# a * b

function reverse_mul!(output, output_deriv, a, b, a_tmp, b_tmp)
    istracked(a) && increment_deriv!(a, mul!(a_tmp, output_deriv, transpose(value(b))))
    istracked(b) && increment_deriv!(b, mul!(b_tmp, transpose(value(a)), output_deriv))
end

for (f, F) in ((:transpose, :Transpose), (:adjoint, :Adjoint))
    @eval begin
        # a * f(b)
        function reverse_mul!(output, output_deriv, a, b::$F, a_tmp, b_tmp)
            _b = ($f)(b)
            istracked(a) && increment_deriv!(a, mul!(a_tmp, output_deriv, mulargvalue(b)))
            istracked(_b) && increment_deriv!(_b, ($f)(mul!(b_tmp, ($f)(output_deriv), value(a))))
        end
           # f(a) * b
        function reverse_mul!(output, output_deriv, a::$F, b, a_tmp, b_tmp)
            _a = ($f)(a)
            istracked(_a) && increment_deriv!(_a, ($f)(mul!(a_tmp, output_deriv, value(b))))
            istracked(b) && increment_deriv!(b, mul!(b_tmp, ($f)(mulargvalue(a)), ($f)(output_deriv)))
        end
        # f(a) * f(b)
        function reverse_mul!(output, output_deriv, a::$F, b::$F, a_tmp, b_tmp)
            _a = ($f)(a)
            _b = ($f)(b)
            istracked(_a) && increment_deriv!(_a, ($f)(mul!(a_tmp, ($f)(mulargvalue(b)), ($f)(output_deriv))))
            istracked(_b) && increment_deriv!(_b, ($f)(mul!(b_tmp, ($f)(output_deriv), ($f)(mulargvalue(a)))))
        end
    end
end

# adjoint(a) * transpose(b)

function reverse_mul!(output, output_deriv, a::Adjoint, b::Transpose, a_tmp, b_tmp)
    _a = adjoint(a)
    _b = transpose(b)
    if istracked(_a)
        reverse_mul!(output, output_deriv, transpose(_a), b, a_tmp, b_tmp)
    elseif istracked(_b)
        increment_deriv!(_b, transpose(mul!(b_tmp, adjoint(output_deriv), adjoint(mulargvalue(a)))))
    end
end

# transpose(a) * adjoint(b)

function reverse_mul!(output, output_deriv, a::Transpose, b::Adjoint, a_tmp, b_tmp)
    _a = transpose(a)
    _b = adjoint(b)
    if istracked(_b)
        reverse_mul!(output, output_deriv, a, transpose(_b), a_tmp, b_tmp)
    elseif istracked(_a)
        increment_deriv!(_a, transpose(mul!(a_tmp, adjoint(mulargvalue(b)), adjoint(output_deriv))))
    end
end

## zero

Base.zero(x::ReverseDiff.TrackedArray) = track(zero(x.value))


text:
using IncoherentScatter
using ElectronDisplay,PlotlyJS
using Tullio,LoopVectorization

import Base: Threads.@threads

function testspectra(;threads=true)
    sv = ScatterVolume(tfeff=230e6, system=1e-6, b=50000e-9, alpha=deg2rad(13));
    freqs = 0.0:20:6400

    tr = Vector{GenericTrace}(undef, 3)

    for k ∈ 1:3
        p     = Plasma(ne=1.5e11, te=(k+1)*500.0, ti=1000.0, sv)
        testr = "$(Int(p.te))"
        tr[k] = scatter(x=freqs, y=pwrspec(p, sv, freqs; threads), showlegend=true,
                        name="Te=$testr K")
    end
    lyo = Layout(;title="IS spectrum, varying Te/Ti
Ne=1.5*1011, Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    for k ∈ 1:3
        p     = Plasma(ne=1.5e11, te=1500.0, ti=1000.0, sv, nion=2, fractn=((k-1)*0.5, 0.0))
        ofstr = "$(Int(p.fractn1*100))"
        mfstr = "$(Int(p.fractn[1]*100))"
        tr[k] = scatter(x=freqs, y=pwrspec(p, sv, freqs), showlegend=true,
                        name="O+ $ofstr %, M+ $mfstr %")
    end
    lyo = Layout(;title="IS spectrum, varying composition
"*
                        "Ne=1.5*1011, Te=1500 K, Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    for k ∈ 1:3
        p     = Plasma(ne=1.0*10^(9+k), te=2000.0, ti=1000.0, sv)
        tr[k] = scatter(x=freqs, y=10^(3-k)*pwrspec(p, sv, freqs), showlegend=true,
                        name="scale 10$(3-k), Ne=10$(9+k) m-3")
    end
    lyo = Layout(;title="IS spectrum, varying Debye length
Te=2000 K, Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    for k ∈ 1:3
        svv = ScatterVolume(tfeff=230e6, system=1e-6, b=50000e-9,
                                             alpha=deg2rad((k-1)*30.0))
        p   = Plasma(ne=1.5e11, te=2000.0, ti=1000.0, sv)
        tr[k] = scatter(x=freqs, y=pwrspec(p, svv, freqs), showlegend=true,
                        name="∠(k, B) = $((k-1)*30) deg")
    end
    lyo = Layout(;title="IS spectrum, varying magnetic field
"*
                        "Ne=1.5*1011, Te=2000 K Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    p = Plasma(ne=1.5e11, te=300.0, ti=300.0, sv)
    for k ∈ 1:3
        cp    = CollisionalPlasma(p; ν=(k-1)*300.0)
        cfstr = "$(Int(cp.ν))"
        tr[k] = scatter(x=freqs, y=pwrspec(cp, sv, freqs), showlegend=true, name="ν=$cfstr Hz")
    end
    lyo = Layout(;title="IS spectrum, ν
"*
                        "Ne=1.5*1011, Te=300 K Ti=300 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))
end

function testft()
    freqs = vcat(0.0:15:4800, 4850:50:7600)
    τ     = 0.0:2e-5:144e-5
    cft   = IncoherentScatter.cosine(freqs, τ)

    sv    = ScatterVolume(tfeff=230e6, system=1e-6, b=50000e-9, alpha=deg2rad(13));
    tr = Vector{GenericTrace}(undef, 3)

    for k ∈ 1:3
        p     = Plasma(ne=1.5e11, te=(k+1)*500.0, ti=1000.0, sv)
        acf   = autocor(cft, pwrspec(p, sv, freqs))
        nm    = "Te=$(Int(p.te)) K"
        tr[k] = scatter(x=τ, y=acf, mode="lines+markers", showlegend=true, name=nm)
    end
    lyo = Layout(;title="ACF, varying Te/Ti
Ne=1.5*1011, Ti=1000 K",
        xaxis_title="Lag in s", yaxis_title="Power", legend_x=0.75)
    display(plot(tr, lyo))

    for k ∈ 1:3
        p     = Plasma(ne=1.5e11, te=1500.0, ti=1000.0, sv, nion=2, fractn=((k-1)*0.5, 0.0))
        nm    = "O+ $(Int(p.fractn1*100)) %, M+ %$(Int(p.fractn[1]*100)) %"
        acf   = autocor(cft, pwrspec(p, sv, freqs))
        tr[k] = scatter(x=freqs, y=acf, mode="lines+markers", name=nm)
    end
    lyo = Layout(;title="IS spectrum, varying composition
"*
                        "Ne=1.5*1011, Te=1500 K, Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    for k ∈ 1:3
        p     = Plasma(ne=1.0*10^(9+k), te=2000.0, ti=1000.0, sv)
        nm    = "scale 10$(3-k), Ne=10$(9+k) m-3"
        acf   = 10^(3-k)*autocor(cft, pwrspec(p, sv, freqs))
        tr[k] = scatter(x=freqs, y=nm, mode="lines+markers", showlegend=true, name=nm)
    end
    lyo = Layout(;title="IS spectrum, varying Debye length
Te=2000 K, Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    for k ∈ 1:3
        svv   = ScatterVolume(tfeff=230e6, system=1e-6, b=50000e-9, alpha=deg2rad((k-1)*30.0))
        p     = Plasma(ne=1.5e11, te=2000.0, ti=1000.0, sv)
        acf   = autocor(cft, pwrspec(p, svv, freqs))
        nm    = "∠(k, B) = $((k-1)*30) deg"
        tr[k] = scatter(x=freqs, y=autocor(cft, pwrspec(p, svv, freqs)), mode="lines+markers",
                        showlegend=true, name=nm)
    end
    lyo = Layout(;title="IS spectrum, varying magnetic field
"*
                        "Ne=1.5*1011, Te=2000 K Ti=1000 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))

    p = Plasma(ne=1.5e11, te=300.0, ti=300.0, sv)
    for k ∈ 1:3
        cp    = CollisionalPlasma(p; ν=(k-1)*300.0)
        acf   = autocor(cft, pwrspec(cp, sv, freqs))
        nm    = "ν=$(Int(cp.ν)) Hz"
        tr[k] = scatter(x=freqs, y=acf, mode="lines+markers", showlegend=true, name=nm)
    end
    lyo = Layout(;title="IS spectrum, ν
"*
                        "Ne=1.5*1011, Te=300 K Ti=300 K",
        xaxis_title="Frequency in Hz", yaxis_title="Power",
        legend_x=0.75)
    display(plot(tr, lyo))
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export feffery_antdtabs

"""
    feffery_antdtabs(;kwargs...)
    feffery_antdtabs(children::Any;kwargs...)
    feffery_antdtabs(children_maker::Function;kwargs...)

An AntdTabs component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): tab children, i.e. contents
- `id` (String; optional)
- `activeKey` (String; optional)
- `animated` (Bool; optional): whether to animate change of tab
- `className` (String; optional)
- `defaultActiveKey` (String; optional)
- `disabled` (Bool; optional)
- `forceRender` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `nClicksAdd` (Real; optional): number of times add button has been clicked
- `nClicksRemove` (Real; optional): number of times add button has been clicked
- `nTabClicks` (Real; optional): number of times add button has been clicked
- `size` (String; optional)
- `style` (Dict; optional): new tab to add on click
- `tabPosition` (String; optional)
- `targetKey` (String; optional)
"""
function feffery_antdtabs(; kwargs...)
        available_props = Symbol[:children, :id, :activeKey, :animated, :className, :defaultActiveKey, :disabled, :forceRender, :loading_state, :nClicksAdd, :nClicksRemove, :nTabClicks, :size, :style, :tabPosition, :targetKey]
        wild_props = Symbol[]
        return Component("feffery_antdtabs", "AntdTabs", "feffery_antd_components", available_props, wild_props; kwargs...)
end

feffery_antdtabs(children::Any; kwargs...) = feffery_antdtabs(;kwargs..., children = children)
feffery_antdtabs(children_maker::Function; kwargs...) = feffery_antdtabs(children_maker(); kwargs...)

max_stars_repo_name: nodelab-aps/feffery-antd-components

max_issues_repo_name: nodelab-aps/feffery-antd-components

max_forks_repo_name: nodelab-aps/feffery-antd-components


text:
# While Loop

i = 0
while i <= 5
  println(i)
  global i += 1
end

# For Loop

for i in 0:5
  println(i)
end

for i in [1,3,"five"]
  println(i)
end

# Comprehensions

x = [i for i=1:50]
println(x)

# Generator Expressions
sum(x^2 for x=1:5)

# Multiple Iterables
for i =0:2, j = 1:2
  println("($i,$j)")
end

# break and continue

i = 0
while true
  println(i)
  if (i >= 5)
    break
  end
  global i += 1
end

for i = -2:5
  if (i < 0)
    println("skipped $i")
    continue
  end
  println(i)
end

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
module ClusterTest

using Test
using ForneyLab

using ForneyLab: Cluster

@testset "Cluster" begin
    g = FactorGraph()

    m = Variable(id=:m)
    v = Variable(id=:v)
    y = Variable(id=:y)
    nd = GaussianMeanVariance(y, m, v)
    em = nd.i[:m].edge
    ev = nd.i[:v].edge

    cluster = Cluster(nd, [em, ev])

    @test cluster.id == :m_v
    @test cluster.node == nd
    @test cluster.edges[1] == em
    @test cluster.edges[2] == ev
end

end # module


text:
using ZigZagBoomerang
using LinearAlgebra
using Random
using SparseArrays
using Distributions
using Statistics
using ZigZagBoomerang: idot, idot_moving!
using ReverseDiff
using Test
using FileIO
using GLM

println("Sparse logistic regression")

# Design matrix
Random.seed!(2)

include("exampledesign.jl")
# create mock design matrix with 2 categorical explanatory variables
# and their interaction effects and 2 continuous explanatory variable
A = example_design_matrix(;num_rows = 50_000)

n, p = size(A)
@show n, p

const m = 1 # m Bernoulli samples for each set of covariates
println("Av. number of regressors per column ", mean(sum(A .!= 0, dims=1)), ", row ", mean(sum(A .!= 0, dims=2)))

# Data from the model
xtrue = randn(p) .* (rand(p) .< 0.2)
sigmoid(x) = inv(one(x) + exp(-x))
lsigmoid(x) = -log(one(x) + exp(-x))
y_ = Matrix(rand(n, m) .< sigmoid.(A*xtrue))
y = sum(y_, dims=2)[:]

# gaussian mixture prior and target potential
pri(x, w1 = 0.5, a=(0.1)^2, w2 = 0.5, b=4.0^2) = -log(w1*exp(-x^2/(2*a))/sqrt(2pi*a) + w2*exp(-x^2/(2*b))/sqrt(2pi*b))
∇pri(x, w1 = 0.5, a=(0.1)^2, w2 = 0.5, b=4.0^2) = -(-(w1*x*exp(-x^2/(2*a)))/a^(3/2) - (w2*x*exp(-x^2/(2*b)))/b^(3/2))/(w1*exp(-x^2/(2*a))/sqrt(a) + w2*exp(-x^2/(2*b))/sqrt(b))
ϕ(x, A, y::Vector) = pri.(x) - sum(y .* lsigmoid.(A*x) + (m .- y) .* lsigmoid.(-A*x))

sigmoidn(x) = sigmoid(-x)
nsigmoid(x) = -sigmoid(x)

# Gradient (found using http://www.matrixcalculus.org/ S. Laue, M. Mitterreiter, and J. Giesen. A Simple and Efficient Tensor Calculus, AAAI 2020.)
∇ϕ(x, A, y) = ∇pri.(x) - A'*(y .* sigmoidn.(A*x)) - A'*((m .- y).*nsigmoid.(A*x))
function ∇ϕ!(out, x, A, y)
    out .= ∇pri.(x) - A'*(y .* sigmoidn.(A*x)) - A'*((m .- y).*nsigmoid.(A*x))
    out
end

function ∇ϕ(x, i, A, y) # partial derivative
    s = ∇pri(x[i])
    n, p = size(A)
    for j in 1:n
        z = sum(A[j,r]*x[r] for r in 1:p)
        s += -A[j,i]*(y[j]*sigmoidn(z) + (m - y[j])*nsigmoid(z))
    end
    return s
end

function ∇ϕr(x, i, A, y, control, μ, bias, k) # with subsampling and optional control variate
    s = ∇pri(x[i]) + control*(bias[i] - pri(μ[i]))
    n, p = size(A)
    sampler = Random.SamplerRangeNDL(1:n)
    for _ in 1:k
        j = rand(sampler)
        z = sum(A[j,r]*x[r] for r in 1:p)
        s += -n/k*A[j,i]*(y[j]*sigmoidn(z) + (m - y[j])*nsigmoid(z))
        if control
            z2 = sum(A[j,r]*μ[r] for r in 1:p)
            s -= -n/k*A[j,i]*(y[j]*sigmoidn(z2) + (m - y[j])*nsigmoid(z2))
        end
    end
    return s
end

# Some Newton steps towards the mode as starting point
# The mode can be used as control variate
println("Newton steps")
x0 = zeros(p)
@time for i in 1:8
    println(i)
    global x0
    H = ReverseDiff.jacobian(x -> ∇ϕ(x, A, y), x0)
    x0 = x0 - cholesky(Symmetric(I + H))\∇ϕ(x0, A, y)
end
@show norm(xtrue - x0)

μ = copy(x0)
Γfull = ReverseDiff.jacobian(x -> ∇ϕ(x, A, y), μ)
Γ = sparse(1.0I, p, p)
bias = ∇ϕ(μ, A, y)
@show norm(bias)

xtest = randn(p)
i = 8
@test ∇ϕ(xtest, A, y)[i] ≈ ∇ϕ(xtest, i, A, y)
#@test ∇ϕ(xtest, i, A, y) ≈ mean(∇ϕr(xtest, i, A, y, false, μ, bias, 100) for k in 1:20000)
#@test ∇ϕ(xtest, i, A, y) ≈ mean(∇ϕr(xtest, i, A, y, true, μ, bias, 100) for k in 1:20000)

# Initial values
t0 = 0.0

# Define ZigZag

#σ = sqrt.(diag(inv(Matrix(Γfull))))
#extrema(σ)
σ = 0.1*(Vector(diag(Γ))).^(-0.5) # cheaper

Z = ZigZag(Γ, μ, σ)

θ0 = rand([-1.0,1.0], p) .* σ

# Run sparse ZigZag for T time units and collect trajectory
println("Distance starting point")
@show norm(x0 - xtrue)

println("Run pdmp")
T = 100.0
c = ones(Float64, p) # rejection bounds
c0 = copy(c)
traj, u, (acc,num), c = @time pdmp(∇ϕr, t0, x0, θ0, T, c, Z, A, y, true, μ, bias, 40; adapt=true)
@show acc, num, acc/num
@show extrema(c ./ c0)

dt = T/4000
x̂ = mean(x for (t,x) in discretize(traj, dt))
X = Float64[]
ts = Float64[]
for (t,x) in discretize(traj, dt)
    append!(ts, t)
    append!(X, x)
end
X = reshape(X, p, length(X)÷p)

[xtrue x̂ abs.(xtrue - x̂)]

@show norm(xtrue - μ)
@show norm(xtrue - x̂)

println("Plot")
using Makie
using Colors
using GoldenSequences

ps = [3, 22, 50, p-3, p-2, p-1]
p0 = length(ps)
cs = map(x->RGB(x...), (Iterators.take(GoldenSequence(3), p0)))
pis = []
for i in 1:p0
    p_i = lines(ts, X[ps[i], :], color=cs[i])
    lines!(p_i, [0, T], [xtrue[ps[i]], xtrue[ps[i]]], color=cs[i], linewidth=2.0)
    ylabel!(p_i, "var$(ps[i])")
    xlabel!(p_i, "t")
    push!(pis, p_i)
end
p1 = title(hbox(pis...), "Spike and slab Bayesian logistic regression n=$(m*n), p=$p", textsize=20)
save(joinpath("figures","spikeandslab$(typeof(Z).name).png"), p1)

p2 = title(vbox(hbox([text("$p", show_axis=false, textsize=10) for p in ps]...),
    [hbox([i==j ? lines(ts, X[i,:]) : lines(X[i, :], X[j, :]) for i in ps]...) for j in ps]...), "$ps")

r = -3:0.01:3
p3 = lines(r, ∇pri.(r))
lines!(p3, r, exp.(-pri.(r)))

p1


text:
# ## Introduction to vector graphics with Luxor.jl

# This course is a short introduction to making vector graphics. By the end of this course you will have learnt quite a bit about how to use Luxor to make vector graphics, and started making some cool images.

# Luxor.jl is a simple package designed to make *vector graphics*, which are pictures made out of colored lines and curves, plus any text. Luxor produces PNG, SVG, and PDF files; PDF and SVG are resolution-independent formats, so you can scale them up or down as much as you like without losing detail. PNG files are, like photos and JPEG images, made at a specific resolution.

# Luxor sits in the space between the heavy-duty Julia plotting packages such as Plots.jl, Gadfly, and Pyplot (among others), and applications such as Inkscape and Adobe Illustrator. Like Plots.jl et al, you specify what you want to draw with Julia functions and parameters using code, rather than by clicking and pointing at the screen with a mouse or trackpad. Like Inkscape/Illustrator, you have the task of specifying every shape and color you want in your final drawing, rather than having the software do everything for you.

# If you've learnt a bit about Julia already, and if you can remember some of the geometry you learned at school (angles, coordinates, and so on), you know enough to get started. To install Luxor, use the Julia package manager:

## using Pkg
## Pkg.add("Luxor")

# # Chapter 1: The first task: design a logo

# The new (and currently fictitious) organization JuliaFission has just asked you to design a new logo for them. They're something to do with atoms and particles, perhaps? So we'll design a new logo for them using some basic shapes; perhaps colored circles in some kind of spiral formation would look suitably "atomic". Let's try out some ideas.

# ### A first drawing

using Luxor

Drawing(500, 500, "my-drawing.png")
origin()
setcolor("red")
circle(Point(0, 0), 100, :fill)
finish()
preview()

# This short piece of code does the following things:

# - Makes a new drawing 500 units square, and saves it in "my-drawing.png" in PNG format.

# - Moves the zero point from the top left to the center. Graphics applications often start measuring from the top left (or bottom left), but it's easier to work out the positions of things if you start at the center.)

# - Selects one of the 200 or so named colors (defined [here](http://juliagraphics.github.io/Colors.jl/stable/))

# - Draws a circle at x = 0, y = 0, with radius 100 units, and fills it with the current color.)

# - Finishes the drawing and displays it on the screen.

# In case you're wondering, the units are *points* (as in font sizes), and there are 72 points in an inch, just over 3 per millimeter. The y-axis points down the page, by the way. If you want to be reminded of where the x and y axes are, uses the `rulers()` function.

# The `:fill` at the end of `circle()` is one of a set of symbols that lets you display the shape in different ways. There's also `:stroke`, which draws around the edges but doesn't fill the shape with color. You might also meet `:path`, `fillpreserve`, `strokepreserve`, `clip`, and `:none`.

# ### Circles in a spiral

# We want more than just one circle. We'll define a triangular shape, and place circles at each corner. The `ngon()` function creates regular polygon (eg triangles, squares, etc.), and the `vertices=true` keyword returns the corner points - just what we want.

Drawing(500, 500, "my-drawing.png")
origin()
setcolor("red")
corners = ngon(Point(0, 0), 80, 3, vertices=true)
circle.(corners, 10, :fill)
finish()
preview()

# Notice the "." after `circle`. This broadcasts the `circle()` function over the `corners`, drawing a red filled circle at every point.

# The arguments to `ngon` are centerpoint, radius, and number of sides. Try changing the third argument from 3 to 4 or 33.

# To create a spiral of circles, we want to repeat this `ngon`...`circle` part more than once. A simple loop will do: we'll rotate everything by `i * ` 5° (`deg2rad(5)` radians) each time (so 5°, 10°, 15°, 20°, 25°, and 30°), and increase the size of the shape by multiples of 10:

Drawing(500, 500, "my-drawing.png")
origin()
setcolor("red")
for i in 1:6
    rotate(i * deg2rad(5))
    corners = ngon(Point(0, 0), 80 + 10i, 3, vertices=true)
    circle.(corners, 10, :fill)
end
finish()
preview()

# ### Just add color

# The Julia colors are available as constants in Luxor, so we can make two changes that cycle through them. The first line creates the set of colors; the `setcolor()` function then works through them. `mod1()` is the 1-based `mod` function, essential for working with Julia and its 1-based indexing.

const colors = (Luxor.julia_green, Luxor.julia_red,Luxor.julia_purple, Luxor.julia_blue) ### <-
Drawing(500, 500, "my-drawing.png")
origin()
for i in 1:6
    setcolor(colors[mod1(i, end)]) ### <-
    rotate(i * deg2rad(5))
    corners = ngon(Point(0, 0), 80 + 10i, 3, vertices=true)
    circle.(corners, 10, :fill)
end
finish()
preview()

# ## Taking particles seriously

# The flat circles are a bit dull, so let's write a function that takes circles seriously. The `drawcircle()` function draws lots of circles, but each one is drawn with a slightly smaller radius and a slightly lighter shade of the incoming color. The `rescale()` function in Luxor provides an easy way to map or adjust values from one range to another; here, numbers between 5 and 1 are mapped to numbers between 0.5 and 3. And the radius is scaled to run between `radius` and `radius/6`. Also, let's make them get larger as they spiral outwards, by adding `4i` to the radius when called by `drawcircle()`.

function drawcircle(pos, radius, n) ## <-
    c = colors[mod1(n, end)]
    for i in 5:-0.1:1
        setcolor(rescale(i, 5, 1, 0.5, 3) .* c)
        circle(pos + (i/2, i/2), rescale(i, 5, 1, radius, radius/6), :fill)
    end
end

Drawing(500, 500, "my-drawing.png")
    origin()
    for i in 1:6
        rotate(i * deg2rad(5))
        corners = ngon(Point(0, 0), 80 + 10i, 3, vertices=true)
        drawcircle.(corners, 10 + 4i, i) ## <-
    end
    finish()
    preview()

# This is looking quite promising. But here's the thing: in a parallel universe, you might already have made this in no time at all using Adobe Illustrator or Inkscape. But with this Luxor code, you can try all kinds of different variations with almost immediate results - just walk through the parameter space, either manually or via code, and see what effects you get. You don't have to redraw everything with different angles and radii...

# So here's what a pentagonal theme with more circles looks like:

Drawing(500, 500, "my-drawing.png")
origin()
for i in 1:12
    rotate(i * deg2rad(1.5))
    corners = ngon(Point(0, 0), 10 + 12i, 5, vertices=true)
    drawcircle.(corners, 5 + 2i, i)
end
finish()
preview()

# To tidy up, it's a good idea to move the code into functions, and do a bit of housekeeping. The final script looks like this:

using Luxor

const colors = (Luxor.julia_green, Luxor.julia_red,Luxor.julia_purple, Luxor.julia_blue)

function drawcircle(pos, radius, n)
    c = colors[mod1(n, end)]
    for i in 5:-.1:1
        setcolor(rescale(i, 5, 1, 0.5, 3) .* c)
        circle(pos + (i/2, i/2), rescale(i, 5, 1, radius, radius/6), :fill)
    end
end

function main(filename)
    Drawing(500, 500, filename)
    origin()
    for i in 1:12
        rotate(i * deg2rad(1.5))
        corners = ngon(Point(0, 0), 10 + 12i, 5, vertices=true)
        drawcircle.(corners, 5 + 2i, i)
    end
    finish()
    preview()
end

main("my-drawing.svg")

# To create a high-quality resolution-independent SVG, just change the suffix from ".png" to ".svg"

# So, did the JuliaFission organization like their logo? Who knows!? But if not, we can always recycle some of these ideas for future projects.

# ### Challenges

# Before you move on to the next chapter, perhaps you'd like to take an extra challenge or two?

# #### 1. Remember the random values

# Using random numbers is a great way to find new patterns and shapes; but unless you know what the values are, you can't reproduce them. So modify the code so that the random numbers are remembered, and drawn on the screen (you can use the `text(str, position)` function),

# #### 2. Backgrounds

# Because there's no background, the SVG or PNG image created will have a transparent background. This is usually what you want for an icon or logo. But this design might look good against a darker colored background. Use `background()` or `paint()` and experiment.

# #### 3. Randomize

# Try refactoring your code so that you can automatically run through various parameter ranges.

# # Chapter 2: Interactive geometry

# Circles, lines, arcs, and polygons are the fundamental particles of Luxor. One way to explore some of the features is to use the interactive abilities of a Jupyter notebook.

# You'll need to add the Interact.jl and WebIO.jl packages if you haven't already done so:

# using Pkg

# Pkg.add("Interact")
# Pkg.add("WebIO")

using WebIO
import Interact: @manipulate

# Interact.jl has a `@manipulate()` macro that provides a slider. This feeds a changing numeric value into the Luxor code that's generating the image. Because we don't want to generate thousands of images for the slider positions, we'll use the in-memory (file-less) version of the `Drawing()` function with the `:svg` option, assign the drawing to `d`, and then return `d` as the result. (You must run this in Jupyter, not Atom/Juno or VS Code.)

@manipulate for r1 in 80:200, r2 in 80:200
    d = Drawing(500, 500, :svg)
    origin()
    sethue("black")
    circle(Point(-100, 0), r1, :stroke)
    circle(Point(100, 0), r2, :stroke)
    flag, ip1, ip2 = intersectioncirclecircle(Point(-100, 0), r1, Point(100, 0), r2)

    if flag
        sethue("red")
        circle(ip1, 12, :fill)
        sethue("green")
        circle(ip2, 10, :fill)
    end
    finish()
    return d
end

# The mathematics going on behind the scenes here avoids the situation when there is a single intersection point, probably because it would demand too much accuracy from the floating-point routines involved. Simple graphics aren't perfectly accurate!

# ## Sample size

# Now let's investigate the `polysample()` function, a useful tool that creates new versions of polygons by "sampling" them a given number of times. First we create a regular polygon, using `ngonside()` rather than `ngon()`, with `nsides` sides and a side length of 100. Then we sample the polygon `s` times.

@manipulate for s in 3:100, nsides=3:12
    d = Drawing(500, 500, :svg)
    origin()
    sethue("black")
    pgon = ngonside(Point(0, 0), 100, nsides, vertices=true)

    ## resample it
    psampled = polysample(pgon, s)

    ## draw the new polygon
    prettypoly(psampled, :stroke, close=true)

    ## show perimeter values
    text(string(polyperimeter(psampled)), halign=:center)
    text(string(polyperimeter(pgon)), Point(0, 40), halign=:center)

    finish()
    return d
end

# We used the `prettypoly()` function to apply a stroke to the sampled polygon. It's called `prettypoly()` because it automatically highlights the vertices for you with small circles - very useful when you're working with unfamiliar polygons.

# We've also introduced the `text(str, pt)` function, that draws a string of text at a point.

# Notice that at low sampling rates (low `s` values) the shape is rendered approximately rather than accurately.

# ## Pretty poly

# With a few changes we can see some more of `prettypoly()` in action. Let's get rid of the sampling, and just vary the size of the polygons drawn at the vertices of the original `pgon`. We pass an *anonymous function* that is applied at each vertex, allowing us to apply any sequence of graphics instructions at each vertex in turn. In this case we select a color at random, draw another polygon centered at the vertex, then outline it in black.

@manipulate for nsides in 3:12, sidelength in 40:80
    d = Drawing(500, 500, :svg)
    origin()
    setopacity(0.6)
    sethue("black")
    pgon = ngonside(Point(0, 0), sidelength, nsides, 0, vertices=true)
    prettypoly(pgon, :none, close=true, () ->
        begin
            randomhue()
            ngon(O, sidelength, nsides, 0, :fillpreserve)
            sethue("black")
            strokepath()
        end)
    finish()
    d
end

# The letter `O` is a short-cut for `Point(0, 0)`. The anonymous function runs after the x and y axes have been moved to the position of the vertex, and therefore all graphics are applied relative to the vertex position.

# ## Specifying locations

# It's useful to be able to refer to specific points on a drawing. To help you do this, there are functions such as `midpoint()` and `between()`, among many others. There are also functions for working with *bounding boxes*, which rectangular areas defined by two opposite corners.

# `BoundingBox()` is a constructor which returns the rectangular extent of various things, such as, in this case, the entire drawing. Multiplying by `k` reduces its size.

# The `between()` function takes two points and a number between 0.0 and 1.0, and finds the point at that interval between them. Here, it's broadcasting over a range, finding 11 points between 0.0 and 1.0, ready for `circle()` to use.

@manipulate for k in 0:0.1:1.0
    d = Drawing(500, 500, :png)
    origin()

    bbox = BoundingBox() * k
    box(bbox, :stroke) ## quick way to draw around a bounding box

    ## draw circles between top left and bottom right corner
    circle.(between.(bbox..., 0:0.1:1.0), 5, :stroke)

    sethue("red")
    ## refer the top left corner of the bounding box
    circle(boxtopleft(bbox), 10, :fill)
    finish()
    return d
end

# A bounding box provides nine useful 'named' positioning functions: `boxtopleft()`, `boxmiddleright()`, and so on.

# ### Challenges

# #### 1. More exploring

# Other interesting Luxor functions to explore might be `rule()`, which can be told to keep the lines inside a bounding box, `pointinverse()`, `intersectionlinecircle()`, and `intersectboundingboxes()`, to name a few.

# # Chapter 3: Conference badges; placing text

# Remember when people used to attend conferences in person, and they'd wear badges to help identify themselves to others? Some day soon, perhaps, we'll be doing it again. Until then, we can work on a useful application for Luxor that will generate these badges automatically. Given a DataFrame containing the details of attendees, Luxor can process the rows and create a unique badge for each person.

# If you want to follow along, install DataFrames.jl, and download the free [Barlow Condensed](https://fonts.google.com/specimen/Barlow+Condensed) font from Google Fonts, because we want a strong but narrow font for the confined space of a badge.

using DataFrames, Colors

# Let's quickly define a few arbitrary attendees (any resemblance to real individuals or companies is unintended):

attendees = DataFrame(
    FirstName      =   ["Alice", "Bob", "Chuck", "Damian", "Elizabeth", "Jacquelinetta", "Gus"],
    LastName       =   ["Blackwood", "Greybeard", "Whitehat", "Gryffindor", "van Purple", "Magickata", "McTurquoisepherson" ],
    Organization   =   ["Avocado Corp",
                        "Sirius Cybernetics Corporation",
                        "Institute of Conference Badge Designers",
                        "Umbrella Corporation",
                        "Abstergo Industries",
                        "Cyberdyne Systems",
                        "Acme Industries"],
    URL            =   ["www.avocadocorp.com"
                         "www.siriuscyberneticscorporation.com"
                         "www.instituteofconferencebadgedesigners.com"
                         "www.umbrellacorporation.com"
                         "www.abstergoindustries.com"
                         "www.cyberdynesystemscorporation.com"
                         "www.acme.com"])

# Some useful information:

const conference_name = "Dark Matters"
const conference_year = "2020"
const card_dimensions = (298, 420)

struct Attendee
    firstname:: String
    lastname:: String
    organization:: String
    url:: String
end

# Luxor provides a few ways to divide up the drawing area, and two useful ones are Tilers and Tables. A Tiler constructor divides up a specified area into identical rectangular tiles; a Table constructor makes rows and columns, but each row can be a different height, and each column can be a different width. The Table constructor can accepts rows/columns, or vectors defining row and column heights. For this task, we'll define a table with 1 column, and 8 rows, and a bit of sketching and experimentation suggests the following row heights.

table = Table([100, 5, 90, 50, 40, 40, 5, 90], [298])

# This quick test shows the basic structure:

Drawing(card_dimensions..., "testcard.svg")
origin()
sethue("black")
for (pos, n) in table
    box(pos, table.colwidths[1], table.rowheights[n], :stroke)
    circle(pos, 4, :fill)
end
finish()
preview()

# The four pieces of information for each attendee are to be placed in table rows 3, 4, 5, and 6. Let's try it out. The throwaway `test_draw_attendee()` function takes a single row of the DataFrame and draws the four pieces of text in these four table cells.

function test_draw_attendee(a, table)
    person = Attendee(a...)
    Drawing(card_dimensions..., "/tmp/testcard-$(person.firstname)-$(person.lastname).png")
    origin()
    sethue("green")
    for (pos, k) in table
        box(pos, table.colwidths[1], table.rowheights[k], :stroke)
    end
    sethue("red")
    text(person.firstname,    table[3])
    text(person.lastname,     table[4])
    text(person.organization, table[5])
    text(person.url,          table[6])
    finish()
    preview()
end

test_draw_attendee(attendees[1, :], table)

# This looks structurally OK. But it's not very appealing visually, so the main job now is to place the text more carefully and use some color.

# We'll make a custom text-drawing function. This draws the text in `str`, and tries to use the font and font size supplied to fit in a space `w` by `h`. But, if the text would be too wide, the font size is reduced until it fits. We can use the `textextents()` function for this, because it returns a useful list of all the dimensions of the text with the current settings.

function drawtext(str, pos, w, h;
        margin=20,
        fsize=20,
        fface="BarlowCondensed-Bold",
        fground = colorant"black",
        bground = colorant"white",
        halignment = :left)
    fontface(fface)
    fontsize(fsize)
    ## check and recalculate if necessary
    textwidth = textextents(str)
    if textwidth[3] > (w - 2margin)
        fs = fsize
        while(textwidth[3] > (w - 2margin))
            fs -= 2
            fontsize(fs)
            textwidth = textextents(str)
        end
        fontsize(fs)
    end
    ## paint the box
    sethue(bground)
    box(pos, w, h, :fillstroke)
    sethue(fground)
    ## allow for centered text
    if halignment == :left
        text(str, pos - (w/2 - margin, 0), halign=:left, valign=:middle)
    else
        text(str, pos - (w/2 - margin, 0), valign=:middle)
    end
end

# Finally, this function draws a single badge. It's mostly calls to `drawtext()`, adding some graphics before and after:

function draw_attendee(a, table)
    person = Attendee(a...)
    Drawing(card_dimensions..., "/tmp/testbadge-$(person.firstname)-$(person.lastname).svg")
    origin()
    background("white")
    sethue("black")
    w = table.colwidths[1]

    @layer begin
        sethue("black")
        box(table, 1, :fill)

        ## dividers
        sethue("red")
        box(table, 2, :fill)
        box(table, 7, :fill)
    end

    drawtext(person.firstname, table[3], w, table.rowheights[3],
        fsize=70, fface="BarlowCondensed-Black",
        fground = "white",
        bground = HSV(230, .8, .2))
    drawtext(person.lastname, table[4], w, table.rowheights[4],
        fsize=40, fface="BarlowCondensed-Bold",
        fground = "white",
        bground = HSV(230, .8, .2))
    drawtext(person.organization, table[5], w, table.rowheights[5],
        fsize=30, fface="BarlowCondensed-Bold",
        fground = "white",
        bground = HSV(230, .6, .5))
    drawtext(person.url, table[6], w, table.rowheights[6],
        fsize=20, fface="BarlowCondensed-Bold",
        fground = "white",
        bground = HSV(230, .5, .7))

    ## what conference is this anyway?
    drawtext(uppercase(conference_name) * "∞" * conference_year, table[8], table.colwidths[1], table.rowheights[8],
        fsize=40, fface="BarlowCondensed-Black",
        fground = "black",
        bground = "darkorange",
        halignment=:center)

    ## final decoration
    @layer begin
        translate(table[1])
        setopacity(0.3)
        setline(0.5)
        sethue("white")
        epitrochoid(200, 31, 250, :stroke)
        end

    finish()
    preview()
end

# Each one can be checked individually:

#-

function testone()
    draw_attendee(attendees[6, :], table)
end

testone()

#-

# And then we can generate them all:

function testall()
    for r in eachrow(attendees)
        println("printing $(r.FirstName * " " * r.LastName)")
        draw_attendee(r, table)
    end
end

testall()

# ### Challenges

# #### 1. Add a QR code

# Adding a QR code to each badge sounds like fun. You can find a nice Julia QR code package [here](https://github.com/jiegillet/QRCode.jl). It's easy enough to generate a grid of black and white squares using Luxor's Table or Tiler objects.

# #### 2. Logos

# People love adding logos to conference-related material. Try adding some to the design. If the logo is provided as a PNG image, you can use Luxor's `readpng()` and `placeimage()` functions. Otherwise you might have to incorporate Luxor code to reproduce the desired logo, because there's no PDF/EPS/SVG import tools. Luckily the Julia logo has already been done.

# #### 3. Provide printer-friendly masters

# We've made a separate image for each badge. But often you would send larger sheets containing four to six badges (possibly complete with crop marks) to the print shop, so that they could run them through their machines and then trim them. So you could refactor the code so that it fits four or six badges on a suitably larger sheet. Perhaps the PDF format is going to be more useful here.

# # Chapter 4: Generative art

# You can make art with anything - these days you might not even need a pencil. Luxor has plenty of tools to help you make some art. The phrase "generative art" suggests that the computer code does most of the work generating the image; perhaps the presence of your guiding hand justifies you claiming to be co-creator.

# You need only have a basic idea to get started. If it doesn't go anywhere, then you just start over. The idea for this chapter is based on a random walk: you place a simple line and circle ("blob") graphic next to other ones, and use randomness to see what emerges.

using Colors

# A simple Blob structure holds start and end points, and keeps track of a color and a radius.

struct Blob
   st::Point
   fi::Point
   col::Colorant
   radius::Float64
end

# This is a place to store Blobs:

todolist = Blob[]

# This function takes an existing Blob object, and returns a new one that's slightly different. Its color, direction, and radius are randomly shifted. A call to `rand(1:10)` introduces some different behaviour: the new Blob is either placed alongside an existing one, or it's attached to the end:

function mutate(blob::Blob)
   currentcol = convert(HSV, blob.col)
   newcol = HSV(mod(currentcol.h + rand(-4:4), 360),
        mod(currentcol.s + 0.1, 1),
        mod(currentcol.v + 0.1, 1))
   newradius = mod(rand(Bool) ? blob.radius * 1.1 : blob.radius * 0.9, 20)
   newdirection = slope(blob.st, blob.fi) + rand(-1:0.1:1) * (π/5 * rand())
   if rand(1:10) < 5
      return Blob(blob.fi, blob.fi + polar(rand(5:10), newdirection), newcol, newradius)
   else
      return Blob(blob.st, blob.st + polar(rand(5:10), newdirection), newcol, newradius)
   end
end

# The `draw()` function draws the blob, and adds new mutated blobs, but only if they're still on the drawing.

function draw(blob::Blob, todolist)
    move(blob.st)
    line(blob.fi)
    setline(4)
    sethue(blob.col)
    strokepath()
    circle(blob.fi, blob.radius, :fill)
    ## add new blobs if they're on the drawing
    if isinside(blob.fi, BoundingBox())
        push!(todolist, mutate(blob))
    end
end

# Because each walk only goes in one direction, a better balanced image is achieved by stacking a dozen images with varying opacities.

# The results are unpredictable, and - this being art - it's difficult to say whether the images are pleasing or horrible. The next time you run it, you might get a masterpiece! Or perhaps you won't...

function gen()
    Drawing(500, 500, "blobart.svg")
    origin()
    background("antiquewhite")
    for i in 1:12
      setopacity(i/12)
      seed = Blob(Point(0, 0), Point(rand(-5:5, 2)...), HSB(180, 0.5, 0.5), 10)
      push!(todolist, seed)
      while !isempty(todolist)
         draw(pop!(todolist), todolist)
      end
    end
    finish()
    preview()
end

gen()

# Fortunately "art is in the eye of the beholder", so if you like these, that's great. And if you don't, you can start thinking of different ways to coax beautiful images from Julia code!

# ## Challenges

# #### 1. Create more masterpieces

# And post your creations online!

max_stars_repo_name: radovankavicky/graphics

max_issues_repo_name: radovankavicky/graphics

max_forks_repo_name: cormullion/graphics


text:
using Oceananigans
using JLD2
using GLMakie
using Printf

xz_file = jldopen("free_convection_xz_slice.jld2")
yz_file = jldopen("free_convection_yz_slice.jld2")
xy_file = jldopen("free_convection_xy_slice.jld2")

Nx = xz_file["grid/Nx"]
Ny = xz_file["grid/Ny"]
Nz = xz_file["grid/Nz"]
Lx = xz_file["grid/Lx"]
Ly = xz_file["grid/Ly"]
Lz = xz_file["grid/Lz"]

grid = RectilinearGrid(size=(Nx, Ny, Nz), x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

# Coordinates

x, y, z = nodes((Center, Center, Face), grid)

slice_height = -8
k_xy_slice = searchsortedfirst(grid.zF[:], slice_height)

z = z[1:k_xy_slice]
Nz′ = length(z)

x_xz = repeat(x, 1, Nz′)
y_xz = 511 * ones(grid.Nx, Nz′)
z_xz = repeat(reshape(z, 1, Nz′), grid.Nx, 1)

x_yz = 511 * ones(grid.Ny, Nz′)
y_yz = repeat(y, 1, Nz′)
z_yz = repeat(reshape(z, 1, Nz′), grid.Ny, 1)

# Slight displacements to "stitch" the cube together
x_xy = x
y_xy = y
z_xy = (slice_height - 0.5) * ones(grid.Nx, grid.Ny)

# Analyze

iterations = parse.(Int, keys(xz_file["timeseries/t"]))

w_max = maximum([maximum(abs, xz_file["timeseries/w/$iter"]) for iter in iterations])
@show w_max

w_lim = 0.5 * w_max

# Animate

iter = Node(0)

wxz = @lift reverse(xz_file["timeseries/w/" * string($iter)][:, 1, 1:k_xy_slice], dims=1)
wyz = @lift reverse(yz_file["timeseries/w/" * string($iter)][1, :, 1:k_xy_slice], dims=1)
wxy = @lift rot180(xy_file["timeseries/w/" * string($iter)][:, :, 1])

fig = Figure(resolution=(1000, 800))

ax = fig[1:10, 1:3] = LScene(fig, aspect=:data, xlabel="x (meters)", ylabel = "y (meters)", zlabel = "height (meters)")

# Three surfaces
pl = surface!(ax, x_xz, y_xz, z_xz, color=wxz, colormap=:bwr, colorrange=(-w_lim, w_lim))
surface!(ax, x_yz, y_yz, z_yz, color=wyz, colormap=:bwr, colorrange=(-w_lim, w_lim))
surface!(ax, x_xy, y_xy, z_xy, color=wxy, colormap=:bwr, colorrange=(-w_lim, w_lim))

# Colorbar
cb = fig[3:8, 4] = Colorbar(fig, pl) #, label="Vertical velocity", height=Relative(0.5))

# Title
title_str = @lift begin
    primitive_title = @sprintf("t = %-16s", prettytime(xz_file["timeseries/t/" * string($iter)]))
    @sprintf("% 52s", primitive_title)
end

title_label = fig[1, 1:3] = Label(fig, title_str, textsize=30)

display(fig)

record(fig, "convection.mp4", iterations; framerate = 16) do i
    iter[] = i
end


text:
using Documenter
using EmojiSymbols

makedocs(
    build = joinpath(@__DIR__, "local" in ARGS ? "build_local" : "build"),
    modules = [EmojiSymbols],
    clean = false,
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        assets = ["assets/custom.css"],
    ),
    sitename = "EmojiSymbols.jl 🤔",
    authors = "WooKyoung Noh",
    pages = Any[
        "Home" => "index.md",
        "additional Emoji symbols" => [
            "additional_symbols_170.md",
            "additional_symbols_below_170.md",
        ],
        "Emoji symbols in REPL" => [
            "emoji_symbols_in_repl_170.md",
            "emoji_symbols_in_repl_below_170.md",
        ],
        "LaTeX symbols in REPL" => "latex_symbols_in_repl.md",
    ],
)


text:
using Test
using CMPFit

## testlinfit
x = [-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
     -2.8394297E-01,1.3416969E+00,1.3757038E+00,
     -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
     8.2065094E-01]

y = [1.9000429E-01,6.5807428E+00,1.4582725E+00,
     2.7270851E+00,5.5969253E+00,5.6249280E+00,
     0.787615,3.2599759E+00,2.9771762E+00,
     4.5936475E+00]

e = fill(0., size(y)) .+ 0.07
param = [1., 1.]

function linfunc(x::Vector{Float64}, p::Vector{Float64})
    return @. p[1] + p[2]*x
end

res = cmpfit(x, y, e, linfunc, param)

@test res.bestnorm ≈ 2.756285  rtol=1.e-5
@test res.status   == 1
@test res.npar     == 2
@test res.nfree    == 2
@test res.npegged  == 0
@test res.dof      == 8

zero = res.perror .- [0.0222102, 0.018938];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

zero = res.param .- [3.209966, 1.770954];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

## testquadfit
x = [-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
	 -2.8394297E-01,1.3416969E+00,1.3757038E+00,
	 -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
	 8.2065094E-01]
y = [2.3095947E+01,2.6449392E+01,1.0204468E+01,
	 5.40507,1.5787588E+01,1.6520903E+01,
	 1.5971818E+01,4.7668524E+00,4.9337711E+00,
	 8.7348375E+00]

e = fill(0., size(y)) .+ 0.2
param = [1., 1., 1.]

function quadfunc(x::Vector{Float64}, p::Vector{Float64})
    return @. p[1] + p[2]*x + p[3]*(x*x)
end

res = cmpfit(x, y, e, quadfunc, param)

@test res.bestnorm ≈ 5.679323  rtol=1.e-5
@test res.status   == 1
@test res.npar     == 3
@test res.nfree    == 3
@test res.npegged  == 0
@test res.dof      == 7

zero = res.perror .- [0.097512, 0.054802, 0.054433];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

zero = res.param .- [4.703829, 0.062586, 6.163087];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

## testquadfix
param = [1., 0., 1.]
pinfo = CMPFit.Parinfo(length(param))
pinfo[2].fixed = 1

res = cmpfit(x, y, e, quadfunc, param, parinfo=pinfo)

@test res.bestnorm ≈ 6.983588  rtol=1.e-5
@test res.status   == 1
@test res.npar     == 3
@test res.nfree    == 2
@test res.npegged  == 0
@test res.dof      == 8

zero = res.perror .- [0.097286, 0., 0.053743];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

zero = res.param .- [4.696254, 0., 6.172954];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

## testgaussfit
x = [-1.7237128E+00,1.8712276E+00,-9.6608055E-01,
	 -2.8394297E-01,1.3416969E+00,1.3757038E+00,
	 -1.3703436E+00,4.2581975E-02,-1.4970151E-01,
	 8.2065094E-01]
y = [-4.4494256E-02,8.7324673E-01,7.4443483E-01,
	 4.7631559E+00,1.7187297E-01,1.1639182E-01,
	 1.5646480E+00,5.2322268E+00,4.2543168E+00,
	 6.2792623E-01]

e = fill(0., size(y)) .+ 0.5
param = [0.0, 1.0, 1.0, 1.0]

function gaussfunc(x::Vector{Float64}, p::Vector{Float64})
    sig2 = p[4] * p[4];
    xc = @. x - p[3];
    return @. p[2] * exp(-0.5 * xc *xc / sig2) + p[1]
end

res = cmpfit(x, y, e, gaussfunc, param)

@test res.bestnorm ≈ 10.350032  rtol=1.e-5
@test res.status   == 1
@test res.npar     == 4
@test res.nfree    == 4
@test res.npegged  == 0
@test res.dof      == 6

zero = res.perror .- [0.232234, 0.395434, 0.074715, 0.089997];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

zero = res.param .- [0.480441, 4.550754, -0.062562, 0.397473];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

## testgaussfit (with explicit derivatives)
gaussfunc_deriv(x::Vector{Float64}, p::Vector{Float64}) = gaussfunc(x, p)
function gaussfunc_deriv(x::Vector{Float64}, p::Vector{Float64}, pderiv::Vector{Int}, vderiv::Vector{Vector{Float64}})
    sig2 = p[4] * p[4];
    xc = x .- p[3];
    ee = exp.(-0.5 .* xc.^2 ./ sig2)

    for i in 1:length(pderiv)
        ipar = pderiv[i]
        if ipar == 1
            vderiv[i] .= 1.
        elseif ipar == 2
            vderiv[i] .= ee
        elseif ipar == 3
            vderiv[i] .= p[2] .* ee .* (x .- p[3]) ./ sig2
        elseif ipar == 4
            vderiv[i] .= p[2] .* ee .* (x .- p[3]).^2 ./ sig2 ./ p[4]
        end
    end

    return p[2] .* ee .+ p[1]
end

pinfo = CMPFit.Parinfo(length(param))
setfield!.(pinfo, :side, Int32(3))
res = cmpfit(x, y, e, gaussfunc_deriv, param, parinfo=pinfo)

## testgaussfix
param = [0.0, 1.0, 0.0, 0.1]
pinfo = CMPFit.Parinfo(length(param))
pinfo[1].fixed = 1
pinfo[3].fixed = 1

res = cmpfit(x, y, e, gaussfunc, param, parinfo=pinfo)

@test res.bestnorm ≈ 15.516134  rtol=1.e-5
@test res.status   == 1
@test res.npar     == 4
@test res.nfree    == 2
@test res.npegged  == 0
@test res.dof      == 8

zero = res.perror .- [0., 0.329307, 0., 0.053804];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5

zero = res.param .- [0., 5.059244, 0., 0.479746];
@test minimum(zero) ≈ 0 atol=1.e-5
@test maximum(zero) ≈ 0 atol=1.e-5


text:
pygui(true)

export useFilterDesignGUI
function useFilterDesignGUI()
win = GtkWindow("Filter Designer")
g = GtkGrid()

LDriveDesc = GtkEntry()  # a widget for entering text
set_gtk_property!(LDriveDesc, :text, "Load Inductance [H]")
makeDescBox(LDriveDesc)
LDriveVal =GtkEntry()
makeValBox(LDriveVal,100e-6)
g[1,1] = LDriveDesc
g[2,1] = LDriveVal

RDriveDesc = GtkEntry()
set_gtk_property!(RDriveDesc, :text, "Load Resistance [Ω]")
makeDescBox(RDriveDesc)
RDriveVal =GtkEntry()
makeValBox(RDriveVal,0.1)
g[1,2] = RDriveDesc
g[2,2] = RDriveVal

CDriveDesc = GtkEntry()
set_gtk_property!(CDriveDesc, :text, "Load Capacitance [F]")
makeDescBox(CDriveDesc)
CDriveVal =GtkEntry()
makeValBox(CDriveVal,1)
g[1,3] = CDriveDesc
g[2,3] = CDriveVal

TargetZDesc = GtkEntry()
set_gtk_property!(TargetZDesc, :text, "Target Impedance [Ω]")
makeDescBox(TargetZDesc)
TargetZVal =GtkEntry()
makeValBox(TargetZVal,20)
g[1,4] = TargetZDesc
g[2,4] = TargetZVal

DriveFreqDesc = GtkEntry()
set_gtk_property!(DriveFreqDesc, :text, "Drive Freq [Hz]")
makeDescBox(DriveFreqDesc)
DriveFreqVal =GtkEntry()
makeValBox(DriveFreqVal,25000)
g[1,5] = DriveFreqDesc
g[2,5] = DriveFreqVal

FilterZDesc = GtkEntry()
set_gtk_property!(FilterZDesc, :text, "Filter Characteristic Impedance [Ω]")
makeDescBox(FilterZDesc)
FilterZVal =GtkEntry()
makeValBox(FilterZVal,25)
g[1,6] = FilterZDesc
g[2,6] = FilterZVal

DampRValDesc = GtkEntry()
set_gtk_property!(DampRValDesc, :text, "Damping Resistor [Ω]")
makeDescBox(DampRValDesc)
DampRVal =GtkEntry()
makeValBox(DampRVal,50)
g[1,7] = DampRValDesc
g[2,7] = DampRVal

OptimButton = GtkCheckButton("Run Optimization?")
Optim_cb = GtkComboBoxText()
choices = ["Dist to Trough", "Magnitude Drift"]
for choice in choices
  push!(Optim_cb,choice)
end
g[1,8] = OptimButton
g[2,8] = Optim_cb

SaveButton = GtkCheckButton("Save LTSPICE File?")

g[1,9] = SaveButton
SavePathBox = GtkEntry()
set_gtk_property!(SavePathBox, :text, "Not saving")
makeDescBox(SavePathBox)
g[2,9] = SavePathBox

id = signal_connect(SaveButton, "button-press-event") do widget, event
    if ~(get_gtk_property(widget, :active, Bool))

        set_gtk_property!(SavePathBox, :has_frame,true)
        set_gtk_property!(SavePathBox, :sensitive,true)
        set_gtk_property!(SavePathBox, :editable,true)
        set_gtk_property!(SavePathBox, :text, "FilterDesignOutput.asc")
    else
        set_gtk_property!(SavePathBox, :has_frame,false)
        set_gtk_property!(SavePathBox, :sensitive,false)

        set_gtk_property!(SavePathBox, :text, "Not saving")
    end
end

Notch1Button = GtkCheckButton("Add notch #1")

g[1,10] = Notch1Button
Notch1Val = GtkEntry()
set_gtk_property!(Notch1Val, :text, "No notch")
makeDescBox(Notch1Val)
g[2,10] = Notch1Val

signal_connect(Notch1Button, "button-press-event") do widget, event
    if ~(get_gtk_property(widget, :active, Bool))

        set_gtk_property!(Notch1Val, :has_frame,true)
        set_gtk_property!(Notch1Val, :sensitive,true)
        set_gtk_property!(Notch1Val, :editable,true)
        set_gtk_property!(Notch1Val, :text, "$(2*getVal(DriveFreqVal))")
    else
        set_gtk_property!(Notch1Val, :has_frame,false)
        set_gtk_property!(Notch1Val, :sensitive,false)

        set_gtk_property!(Notch1Val, :text, "Not notching 1")
    end
end

Notch2Button = GtkCheckButton("Add notch #2")

g[1,11] = Notch2Button
Notch2Val = GtkEntry()
set_gtk_property!(Notch2Val, :text, "No notch")
makeDescBox(Notch2Val)
g[2,11] = Notch2Val

signal_connect(Notch2Button, "button-press-event") do widget, event
    if ~(get_gtk_property(widget, :active, Bool))

        set_gtk_property!(Notch2Val, :has_frame,true)
        set_gtk_property!(Notch2Val, :sensitive,true)
        set_gtk_property!(Notch2Val, :editable,true)
        set_gtk_property!(Notch2Val, :text, "$(3*getVal(DriveFreqVal))")
    else
        set_gtk_property!(Notch2Val, :has_frame,false)
        set_gtk_property!(Notch2Val, :sensitive,false)

        set_gtk_property!(Notch2Val, :text, "Not notching 2")
    end
end

FindFreqButton = GtkCheckButton("Determine ideal freq?")
g[1,12] = FindFreqButton

RunScriptButton = GtkButton()
set_gtk_property!(RunScriptButton, :label, "RUN DESIGNER")
g[1:2,13] = RunScriptButton

set_gtk_property!(Optim_cb,:active,1)

id = signal_connect(RunScriptButton, "button-press-event") do widget, event
    println("Starting filter designer")

    if  (get_gtk_property(SaveButton, :active, Bool))
            Savepath = get_gtk_property(SavePathBox, :text,String)
            println("The SPICE file will be saved at location:"*get_gtk_property(SavePathBox, :text,String))
            if ~(Savepath[end-3:end]==".asc")
                Savepath = Savepath*".asc"
            end
            set_gtk_property!(SavePathBox, :text, Savepath)
    else
        Savepath = nothing
        println("Will not save schematic")
    end

    NotchArray = [1e9,1e9]
    if  (get_gtk_property(Notch1Button, :active, Bool))
        NotchArray[1] = getVal(Notch1Val)
    end

    if  (get_gtk_property(Notch2Button, :active, Bool))
        NotchArray[2] = getVal(Notch2Val)
    end

    if  (get_gtk_property(OptimButton, :active, Bool))

        MinimizeDistToTrough = parse(Int,get_gtk_property(Optim_cb,:active,String))==0

        MinZ,DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray = DesignDriveFilter_OptimDrift(getVal(LDriveVal),getVal(RDriveVal),getVal(TargetZVal),getVal(DriveFreqVal); MinimizeDistToTrough=MinimizeDistToTrough,
        CDrive = getVal(CDriveVal),
        NumDriveElements = 1,
        WireDiam = 2e-3,
        WireFillFac = 0.75,
        PlotSrcR = 0.0001,
        PlotFTs = true,
        VSrc = 2,
        DetermineFreq = (get_gtk_property(FindFreqButton, :active, Bool)),
        AddNotchFreq = NotchArray,
        FilterZ = getVal(FilterZVal),
        RDampVal = getVal(DampRVal),
        PerturbTxReactance = nothing,
        WriteFileName=Savepath
        )
    else
    DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray = DesignDriveFilter(getVal(LDriveVal),getVal(RDriveVal),getVal(TargetZVal),getVal(DriveFreqVal);
        CDrive = getVal(CDriveVal),
        NumDriveElements = 1,
        WireDiam = 2e-3,
        WireFillFac = 0.75,
        PlotSrcR = 0.0001,
        PlotFTs = true,
        VSrc = 2,
        DetermineFreq = (get_gtk_property(FindFreqButton, :active, Bool)),
        AddNotchFreq = NotchArray,
        FilterZ = getVal(FilterZVal),
        RDampVal = getVal(DampRVal),
        PerturbTxReactance = nothing,
        AxesArray = nothing,
        WriteFileName=Savepath
        )
    end
        println(SPICE_DF)

        if  (get_gtk_property(SaveButton, :active, Bool))

            Savepath_CSV = Savepath[1:end-4]*".csv"
            println(Savepath_CSV)
            CSV.write(Savepath_CSV, SPICE_DF)

        end
end

set_gtk_property!(g, :column_homogeneous, true)
set_gtk_property!(g, :column_spacing, 15)  # introduce a 15-pixel gap between columns
push!(win, g)
showall(win)

end

function makeDescBox(Box)
    set_gtk_property!(Box, :has_frame,false)
    set_gtk_property!(Box, :sensitive,false)
    set_gtk_property!(Box, :editable,false)
end

function makeValBox(Box,DefaultVal)
    set_gtk_property!(Box, :text, "$DefaultVal")
end

function getVal(Box)
    parse(Float64,get_gtk_property(Box, :text, String))
end


text:
using Test
using Py2Jl
include("test_case1.jl")