text:
### utils.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is MIT "Expat".
#
### Commentary
#
# This file contains some utilities for the LombScargle.jl package.
#
### Code:

"""
    power(p::Periodogram)

Return the power vector of Lomb-Scargle periodogram `p`.
"""
power(p::Periodogram) = p.power

"""
    freq(p::Periodogram)

Return the frequency vector of Lomb-Scargle periodogram `p`.
"""
freq(p::Periodogram) = p.freq

"""
    freqpower(p::Periodogram)

Return the 2-tuple `(freq(p), power(p))`, where `freq(p)` and `power(p)` are the
frequency vector and the power vector of Lomb-Scargle periodogram `p`
respectively.
"""
freqpower(p::Periodogram) = (freq(p), power(p))

"""
    findmaxpower(p::Periodogram)

Return the highest power of the periodogram `p`.
"""
findmaxpower(p::Periodogram) = maximum(power(p))

_findmaxfreq{R1<:Real,R2<:Real}(freq::AbstractVector{R1},
                                power::AbstractVector{R2},
                                threshold::Real) =
                                    freq[find(x -> x >= threshold, power)]

"""
    findmaxfreq(p::Periodogram, [interval::AbstractVector{Real}], threshold::Real=findmaxpower(p))

Return the array of frequencies with the highest power in the periodogram `p`.
If a scalar real argument `threshold` is provided, return the frequencies with
power larger than or equal to `threshold`.  If you want to limit the search to a
narrower frequency range, pass as second argument a vector with the extrema of
the interval.
"""
findmaxfreq(p::Periodogram, threshold::Real=findmaxpower(p)) =
    _findmaxfreq(freqpower(p)..., threshold)

function findmaxfreq{R<:Real}(p::Periodogram, interval::AbstractVector{R},
                              threshold::Real=NaN)
    f = freq(p)
    lo, hi = extrema(interval)
    indices = find(x -> lo <= x <= hi, f)
    pow = power(p)[indices]
    if isnan(threshold)
        threshold = maximum(pow)
    end
    return _findmaxfreq(f[indices], pow, threshold)
end

"""
    autofrequency(times::AbstractVector{Real};
                  samples_per_peak::Integer=5,
                  nyquist_factor::Integer=5,
                  minimum_frequency::Real=NaN,
                  maximum_frequency::Real=NaN)

Determine a suitable frequency grid for the given vector of `times`.

Optional keyword arguments are:

* `samples_per_peak`: the approximate number of desired samples across the
  typical peak
* `nyquist_factor`: the multiple of the average Nyquist frequency used to choose
  the maximum frequency if `maximum_frequency` is not provided
* `minimum_frequency`: if specified, then use this minimum frequency rather than
  one chosen based on the size of the baseline
* `maximum_frequency`: if specified, then use this maximum frequency rather than
  one chosen based on the average Nyquist frequency

This is based on prescription given at
https://jakevdp.github.io/blog/2015/06/13/lomb-scargle-in-python/ and uses the
same keywords names adopted in Astropy.
"""
function autofrequency{R<:Real}(times::AbstractVector{R};
                                samples_per_peak::Integer=5,
                                nyquist_factor::Integer=5,
                                minimum_frequency::Real=NaN,
                                maximum_frequency::Real=NaN)
    T = maximum(times) - minimum(times)
    δf = inv(samples_per_peak*T)
    f_min = isfinite(minimum_frequency) ? minimum_frequency : 0.5*δf
    if isfinite(maximum_frequency)
        return f_min:δf:maximum_frequency
    else
        return f_min:δf:0.5*nyquist_factor*length(times)/T
    end
end

"""
    prob(P::Periodogram, p_0::Real)

Return the probability that the periodogram power can exceed the value `p_{0}`.

Its inverse is the `probinv` function.
"""
function prob(P::Periodogram, p_0::Real)
    N = length(P.times)
    normalization = P.norm
    if normalization == "standard"
        return (1.0 - p_0)^((N - 3.0)*0.5)
    elseif normalization == "Scargle"
        return exp(-p_0)
    elseif normalization == "HorneBaliunas"
        return (1.0 - 2p_0/(N - 1))^((N - 3.0)*0.5)
    elseif normalization == "Cumming"
        return (1.0 + 2p_0/(N - 3.0))^((3.0 - N)*0.5)
    else
        error("normalization \"", normalization, "\" not supported")
    end
end

"""
    probinv(P::Periodogram, prob::Real)

Return the `p_0` value of the periodogram power whose probability is `prob`.

This is the inverse of `prob` function.
"""
function probinv(P::Periodogram, prob::Real)
    N = length(P.times)
    normalization = P.norm
    if normalization == "standard"
        return 1.0 - prob^(2.0/(N - 3.0))
    elseif normalization == "Scargle"
        return -log(prob)
    elseif normalization == "HorneBaliunas"
        return 0.5*(N - 1.0)*(1.0 - prob^(2.0/(N - 3.0)))
    elseif normalization == "Cumming"
        return 0.5*(N - 3.0)*(prob^(2.0/(3.0 - N)) - 1.0)
    else
        error("normalization \"", normalization, "\" not supported")
    end
end

"""
    M(P::Periodogram)

Estimates the number of independent frequencies in the periodogram `P`.
"""
function M(P::Periodogram)
    t = P.times
    f = P.freq
    return (maximum(t) - minimum(t))*(maximum(f) - minimum(f))
end

"""
    fap(P::Periodogram, p_0::Real)

Return the false-alarm probability for periodogram `P` and power value `p_0`.

Its inverse is the `fapinv` function.
"""
fap(P::Periodogram, p_0::Real) = 1.0 - (1.0 - prob(P, p_0))^M(P)

"""
    fapinv(P::Periodogram, fap::Real)

Return the `p_0` value of the periodogram power whose false-alarm probability is
`fap`.

This is the inverse of `fap` function.
"""
fapinv(P::Periodogram, fap::Real) = probinv(P, 1.0 - (1.0 - fap)^(inv(M(P))))

"""
    LombScargle.model(times::AbstractVector{Real},
                      signal::AbstractVector{R2},
                      [errors::AbstractVector{R3},]
                      frequency::Real,
                      [times_fit::AbstractVector{R4}];
                      center_data::Bool=true,
                      fit_mean::Bool=true)

Return the best fitting Lomb-Scargle model for the given signal at the given
frequency.

Mandatory arguments are:

* `times`: the observation times
* `signal`: the signal, sampled at `times` (must have the same length as
  `times`)
* `frequency`: the frequency at which to calculate the model

Optional arguments are:

* `errors`: the vector of uncertainties of the signal.  If provided, it must
  have the same length as `signal` and `times`, and be the third argument
* `times_fit`: the vector of times at which the model will be calculated.  It
  defaults to `times`.  If provided, it must come after `frequency`

Optional keyword arguments `center_data` and `fit_mean` have the same meaning as
in `lombscargle` function, which see.
"""
function model{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(t::AbstractVector{R1},
                                                    s::AbstractVector{R2},
                                                    errors::AbstractVector{R3},
                                                    f::Real,
                                                    t_fit::AbstractVector{R4}=t;
                                                    center_data::Bool=true,
                                                    fit_mean::Bool=true)
    @assert length(t) == length(s) == length(errors)
    if center_data
        # Compute weights vector
        w = 1.0./errors.^2
        m = (w⋅s)/sum(w)
        y = (s - m)./errors
    else
        # We don't want to center the data: the mean is 0 and the signal is left
        # unchanged
        m = zero(R2)
        y = s./errors
    end
    # The Lomb-Scargle periodogram looks for the best fitting sinusoidal
    # function
    #     a·cos(ωt) + b·sin(ωt) + c
    # In order to find the coefficients a, b, c for the given frequency we can
    # solve the linear system A·x = y, where A is the matrix with rows:
    # [cos(ωt), sin(ωt), 1]; x is the column vector [a, b, c], and y is the
    # column vector of the signal
    ω = 2pi*f
    if fit_mean
        a, b, c = hcat(cos(ω*t), sin(ω*t), ones(t))./errors \ y
        return a*cos(ω*t_fit) + b*sin(ω*t_fit) + (c + m)
    else
        # If fit_mean==false, the model to be fitted is a·cos(ωt) + b·sin(ωt)
        a, b = hcat(cos(ω*t), sin(ω*t))./errors \ y
        return a*cos(ω*t_fit) + b*sin(ω*t_fit) + m
    end
end

# No uncertainties: errors=ones(s)
model{R1<:Real,R2<:Real,R3<:Real}(t::AbstractVector{R1},
                                  s::AbstractVector{R2},
                                  f::Real,
                                  t_fit::AbstractVector{R3}=t;
                                  kwargs...) =
                                      model(t, s, ones(s), f, t_fit; kwargs...)

# Uncertainties provided via Measurement type
model{R1<:Real,R2<:AbstractFloat,R3<:Real}(t::AbstractVector{R1},
                                           s::AbstractVector{Measurement{R2}},
                                           f::Real,
                                           t_fit::AbstractVector{R3}=t;
                                           kwargs...) =
                                               model(t,
                                                     value(s),
                                                     uncertainty(s),
                                                     f,
                                                     t_fit;
                                                     kwargs...)


text:
using Plots
using DataFrames
using HTTP
using JSON

function get_current_price()
   #REFERENCE: https://developer.oanda.com/rest-live-v20/pricing-ep/#CurrentPrices
   instrument = "EUR_USD"
   account_id = "123-456-7890123-456"
   api_key = "9r32rjijfiewjfoihfoi3u9439243294329432-094328hhrfiehfiheip81724jhfd"
   url = string("https://api-fxpractice.oanda.com/v3/accounts/", account_id, "/pricing?instruments=", instrument)

   r = HTTP.request(
      "GET",
      url,
      ["Accept" => "application/json",
      "Content-Type" => "application/json",
      "Authorization" => string("Bearer ", api_key)]
      )

   res = r.body

   #The body of a request or response can be cast to a string - as below - or you can read and write bytes off of the stream as needed.
   res_in_string = String(res)

   #JSON converts string by curl response to a dictionary object
   data_in_dict = JSON.parse(res_in_string)

   #JSON returns a dictionary of dictionaries. Use [1] to get to the one-level further in the dictionary
   instrument_code = data_in_dict["prices"][1]["instrument"]
   bid = data_in_dict["prices"][1]["closeoutBid"]
   ask = data_in_dict["prices"][1]["closeoutAsk"]

   println("Code: ", instrument_code, "\nBid: ", bid, "\nAsk: ", ask)
end

function get_historical_data()
   #REFERENCE: http://developer.oanda.com/rest-live/rates/#retrieveInstrumentHistory
   # THIS IS V1 AND SHOULD HAVE BEEN ALREADY DEPRICATED???

   instrument = "EUR_USD"
   account_id = "123-456-7890123-456"
   api_key = "9r32rjijfiewjfoihfoi3u9439243294329432-094328hhrfiehfiheip81724jhfd"

   #With count paramter of 5000 and daily ("D") data  (NO start or end date to be included!)
   url = string("https://api-fxpractice.oanda.com/v1/candles?instrument=", instrument, "&granularity=D&count=5000"

   r = HTTP.request(
   "GET",
   url,
   ["Accept" => "application/json",
    "Content-Type" => "application/json",
    "Authorization" => string("Bearer ", api_key) ]
    )
   response = r.body
   data = JSON.parse(String(response))

   #Convert every single array into an instance of Dataframe and then concatenate each one of them with "vcat"
   #Without the splat (...), a warning gets generated, but it executes without a failure
   df = vcat(DataFrame.(data["candles"])...)

   #Load GR backend of Plots (GR is its default)
   gr()

   # Use comprehension to add each datapoint to OHLC type
   y = OHLC[(df[i, :openAsk], df[i, :highAsk], df[i, :lowAsk], df[i, :closeAsk]) for i = 1:nrow(df)]

   title!("EUR_USD")

   graph = ohlc(y)

   #All output is suppressed when running from a file, one needs to "explicitly" display by command
   #I have tried the following methods each, but none of them seems to work on Atom IDE Plot window...
   display(graph)
   gui(graph)
   show = true
end

function get_stream_pricing()
   instrument = "EUR_USD"
   account_id = "123-456-7890123-456"
   api_key = "9r32rjijfiewjfoihfoi3u9439243294329432-094328hhrfiehfiheip81724jhfd"
   url = string("https://stream-fxpractice.oanda.com/v3/accounts/", account_id, "/pricing/stream?instruments=", instrument)

   r = HTTP.request(
       "GET",
       url,
       ["Accept" => "application/json",
        "Content-Type" => "application/json",
        "Authorization" => string("Bearer ", api_key) ]
        )
   response = r.body
end

ext: jl

max_stars_repo_name: tezzytezzy/Oanda-REST-API-in-Julia

max_issues_repo_name: tezzytezzy/Oanda-REST-API-in-Julia

max_forks_repo_name: tezzytezzy/Oanda-REST-API-in-Julia


text:
# Uses a special case of the method in
# Kay, T. L., & Kajiya, J. T. (1986). Ray tracing complex scenes
#
# Assumes the line passes all the way through the AABox if it intersects, which is a
# valid assumption for this ray tracing application.
#
# This version is branchless and is likely faster on the GPU
#function intersect(l::LineSegment, aab::AABox)
#    𝘂⁻¹= 1 ./ l.𝘂
#    𝘁₁ = 𝘂⁻¹*(aab.minima - l.𝘅₁)
#    𝘁₂ = 𝘂⁻¹*(aab.maxima - l.𝘅₁)
#    tmin = maximum(min.(𝘁₁, 𝘁₂))
#    tmax = minimum(max.(𝘁₁, 𝘁₂))
#    return (tmax >= tmin, SVector(l(tmin), l(tmax)))
#end

# This version has branches and is slightly faster on CPU
# Section 5.3.3 in Ericson, C. (2004). Real-time collision detection
function intersect(l::LineSegment{N,T}, aab::AABox{N,T}) where {N,T}
    tmin = typemin(T)
    tmax = typemax(T)
    p_nan = nan(typeof(l.𝘅₁))
    for i = 1:N
        if abs(l.𝘂[i]) < 1e-6
            if l.𝘅₁[i] < aab.minima[i] || aab.maxima[i] < l.𝘅₁[i]
                return (false, SVector(p_nan, p_nan))
            end
        else
            u⁻¹= 1/l.𝘂[i]
            t₁ = (aab.minima[i] - l.𝘅₁[i])*u⁻¹
            t₂ = (aab.maxima[i] - l.𝘅₁[i])*u⁻¹
            if t₁ > t₂
                t₁,t₂ = t₂,t₁
            end
            tmin = max(tmin, t₁)
            tmax = min(tmax, t₂)
            if tmin > tmax
                return (false, SVector(p_nan, p_nan))
            end
        end
    end
    return (true, SVector(l(tmin), l(tmax)))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
@testset ExtendedTestSet "hausdorff" begin
    @testset ExtendedTestSet "hausdorff" begin
        y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = cat(y1, y1; dims=3)
        ŷ = copy(y)

        ŷ_dtm1 = [
            3.0  2.0  1.0  0.0
            3.0  2.0  1.0  0.0
            3.0  2.0  1.0  0.0
            3.0  2.0  1.0  0.0
        ]
        ŷ_dtm = cat(ŷ_dtm1, ŷ_dtm1, dims=3)
        y_dtm = copy(ŷ_dtm)

        @test hausdorff(ŷ, y, ŷ_dtm, y_dtm) == 0.0
    end

    @testset ExtendedTestSet "hausdorff" begin
        y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = cat(y1, y1; dims=3)
        ŷ = copy(y)

        ŷ_dtm1 = [
            0.0  0.0  0.0  3.0
            0.0  0.0  0.0  3.0
            0.0  0.0  0.0  3.0
            0.0  0.0  0.0  3.0
        ]
        ŷ_dtm = cat(ŷ_dtm1, ŷ_dtm1, dims=3)
        y_dtm = copy(ŷ_dtm)

        @test hausdorff(ŷ, y, ŷ_dtm, y_dtm) == 0.0
    end
end


text:
function solcCompiler()

    # getVersion = `./solcv/solc-macosx-amd64-v0.8.9+commit.e5eed63a solc --version`
    runCommand = `./solcv/solc-macosx-amd64-v0.8.9+commit.e5eed63a solc --help`
    run(runCommand)
end


text:
immutable Exponential <: ContinuousUnivariateDistribution
    scale::Float64 # note: scale not rate
    function Exponential(sc::Real)
        if sc > 0.0
            new(float64(sc))
        else
            error("scale must be positive")
        end
    end
end

Exponential() = Exponential(1.0)

cdf(d::Exponential, q::Real) = q <= 0.0 ? 0.0 : -expm1(-q / d.scale)

function logcdf(d::Exponential, q::Real)
    if q <= 0.0
        return -Inf
    else
        qs = -q / d.scale
        if qs > log(0.5)
            return log(-expm1(qs))
        else
            return log1p(-exp(qs))
        end
    end
end

function ccdf(d::Exponential, q::Real)
    q <= 0.0 ? 1.0 : exp(-q / d.scale)
end

function logccdf(d::Exponential, q::Real)
    q <= 0.0 ? 0.0 : -q / d.scale
end

function invlogcdf(d::Exponential, lp::Real)
    if lp <= 0.0
        -d.scale * (lp > log(0.5) ? log(-expm1(lp)) : log1p(-exp(lp)))
    else
        return NaN
    end
end

function invlogccdf(d::Exponential, lp::Real)
    lp <= 0.0 ? -d.scale * lp : NaN
end

entropy(d::Exponential) = 1.0 - log(d.scale)

insupport(d::Exponential, x::Number) = isreal(x) && isfinite(x) && 0.0 <= x

kurtosis(d::Exponential) = 6.0

mean(d::Exponential) = d.scale

median(d::Exponential) = d.scale * log(2.0)

function mgf(d::Exponential, t::Real)
    s = d.scale
    return (1.0 - t * s)^(-1)
end

function cf(d::Exponential, t::Real)
    s = d.scale
    return (1.0 - t * im * s)^(-1)
end

modes(d::Exponential) = [0.0]

function pdf(d::Exponential, x::Real)
    x <= 0.0 ? 0.0 : exp(-x / d.scale) / d.scale
end

function logpdf(d::Exponential, x::Real)
    x <= 0.0 ? -Inf : -x / d.scale - log(d.scale)
end

function quantile(d::Exponential, p::Real)
    0.0 <= p <= 1.0 ? -d.scale * log1p(-p) : NaN
end

function cquantile(d::Exponential, p::Real)
    0.0 <= p <= 1.0 ? -d.scale * log(p) : NaN
end

rand(d::Exponential) = d.scale * Random.randmtzig_exprnd()

rand!(d::Exponential, A::Array{Float64}) = d.scale * Random.randmtzig_fill_exprnd!(A)

skewness(d::Exponential) = 2.0

var(d::Exponential) = d.scale * d.scale


text:
println("Julia script executing!")

@assert snakemake.config["test"] == true
@assert snakemake.params["integer"] == 123
@assert snakemake.output[1] == "julia.out"
@assert snakemake.input[1] == "test.in"
@assert snakemake.input["named_input"] == "test.in"

f = open(snakemake.output[1], "w")
println(f, "Julia test succeded!")
close(f)

max_stars_repo_name: graingert/snakemake

max_issues_repo_name: graingert/snakemake

max_forks_repo_name: graingert/snakemake


text:
## Dimension Test
## ==============
using Sigma
import Sigma: SimplexBenchmark, weighted_mid_split, rand_partial_split,
              weighted_partial_split
import Sigma: dreal, z3, dreal3
import Sigma: SigmaSMT, SigmaAI
import Sigma: cond_sample_tlmh
import Sigma: runbenchmarks

SplitBenchmarks = [Simplex(i,[:sample_distribution,
                                        :accumulative_KL,
                                        :total_time,],400,.1)
                   for i = [3,10]]

mh_captures = [:start_loop, :refinement_depth]
all_splits = [weighted_partial_split, rand_partial_split]

# SMT algorithms
SMTAlgorithms = [SigmaSMT(mh_captures, solver, sampler, nprocs, split)
  for nprocs = [1],
      solver = [dreal3],
      split = all_splits,
      sampler = [cond_sample_tlmh]][:]

AIAlgorithms = [SigmaAI(mh_captures, sampler, nprocs, split)
  for nprocs = [1],
      split = all_splits,
      sampler = [cond_sample_tlmh]][:]

dimbenchmarks = SplitBenchmarks
dimalgorithms = vcat(AIAlgorithms,SMTAlgorithms)

runbenchmarks(AIAlgorithms,dimbenchmarks;runname = "depth")

# # function extractdata(data)
#     xs = Int64[]
#     ys = Vector[]
#     for (key,value) in data
#       if (key[1].split == Sigma.rand_partial_split) & !(isa(value,Exception))
#         push!(xs,key[2].ndims)
#         push!(ys,value[:refinement_depth])
#       end
#     end


text:
# # [Using GXBeam with DifferentialEquations.jl](@id diffeq)
#
# While the capabilities provided by GXBeam are probably sufficient for most users,
# advanced users may wish to make use of some of the features of the
# [`DifferentialEquations`](https://github.com/SciML/DifferentialEquations.jl) package.
# For this reason, we have created an interface in GXBeam to allow users to model the
# differential algebraic equations encountered in GXBeam in DifferentialEquations.
#
#md # ```@contents
#md # Pages = ["diffeq.md"]
#md # Depth = 3
#md # ```
#
#-
#md # !!! tip
#md #     This guide is also available as a Jupyter notebook:
#md #     [`diffeq.ipynb`](@__NBVIEWER_ROOT_URL__/examples/diffeq.ipynb).
#-
#
# ## Interface Functions
#
# The following constructors are available for modeling the differential algebraic
# equations from GXBeam in DifferentialEquations.
#
# ```@docs
# GXBeam.ODEFunction(system::System, assembly; kwargs...)
# GXBeam.ODEProblem(system::System, assembly, tspan; kwargs...)
# GXBeam.DAEFunction(system::System, assembly; kwargs...)
# GXBeam.DAEProblem(system::System, assembly, tspan; kwargs...)
# ```
#
# ## Example Usage
#
# For this example we demonstrate how to solve the [Nonlinear Dynamic Analysis of a Wind
# Turbine Blade](@ref) problem using DifferentialEquations.
#
# We start by setting up the problem as if we were solving the problem using GXBeam's
# internal solver.

using GXBeam, LinearAlgebra

L = 60 # m

## create points
nelem = 10
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]

## index of endpoints of each beam element
start = 1:nelem
stop = 2:nelem+1

## stiffness matrix for each beam element
stiffness = fill(
    [2.389e9  1.524e6  6.734e6 -3.382e7 -2.627e7 -4.736e8
     1.524e6  4.334e8 -3.741e6 -2.935e5  1.527e7  3.835e5
     6.734e6 -3.741e6  2.743e7 -4.592e5 -6.869e5 -4.742e6
    -3.382e7 -2.935e5 -4.592e5  2.167e7 -6.279e5  1.430e6
    -2.627e7  1.527e7 -6.869e5 -6.279e5  1.970e7  1.209e7
    -4.736e8  3.835e5 -4.742e6  1.430e6  1.209e7  4.406e8],
    nelem)

## mass matrix for each beam element
mass = fill(
    [258.053      0.0        0.0      0.0      7.07839  -71.6871
       0.0      258.053      0.0     -7.07839  0.0        0.0
       0.0        0.0      258.053   71.6871   0.0        0.0
       0.0       -7.07839   71.6871  48.59     0.0        0.0
       7.07839    0.0        0.0      0.0      2.172      0.0
     -71.6871     0.0        0.0      0.0      0.0       46.418],
     nelem)

## create assembly of interconnected nonlinear beams
assembly = Assembly(points, start, stop; stiffness=stiffness, mass=mass)

## prescribed conditions
prescribed_conditions = (t) -> begin
    Dict(
        ## fixed left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        ## force on right side
        nelem+1 => PrescribedConditions(Fz = 1e5*sin(20*t))
    )
end

#!jl nothing #hide

# At this point if we wanted to use GXBeam's internal solver, we would choose a time
# discretization and call the `time_domain_analysis` function.

## simulation time
t = 0:0.001:2.0

system, gxbeam_history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions = prescribed_conditions)

#!jl nothing #hide

# To instead use the capabilities of the DifferentialEquations package we can do the
# following.

using DifferentialEquations

## define simulation time
tspan = (0.0, 2.0)

## run initial condition analysis to get consistent set of initial conditions
system, converged = initial_condition_analysis(assembly, tspan[1]; prescribed_conditions)

## construct a DAEProblem
prob = DAEProblem(system, assembly, tspan; prescribed_conditions)

## solve DAEProblem
sol = solve(prob, DABDF2())

#!jl nothing #hide

# We can extract the outputs from the solution in a easy to understand format using the
# [`AssemblyState`](@ref) constructor.

diffeq_history = [AssemblyState(system, assembly, sol[it]; prescribed_conditions)
    for it in eachindex(sol)]

#!jl nothing #hide

# Let's now compare the solutions from GXBeam's internal solver and the `DABDF2` solver
# from DifferentialEquations.

using Plots
#md using Suppressor #hide
pyplot()

#md @suppress_err begin #hide

point = vcat(fill(nelem+1, 6), fill(1, 6))
field = [:u, :u, :u, :theta, :theta, :theta, :F, :F, :F, :M, :M, :M]
direction = [1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]
ylabel = ["\$u_x\$ (\$m\$)", "\$u_y\$ (\$m\$)", "\$u_z\$ (\$m\$)",
    "Rodriguez Parameter \$\\theta_x\$ (degree)",
    "Rodriguez Parameter \$\\theta_y\$ (degree)",
    "Rodriguez Parameter \$\\theta_z\$ (degree)",
    "\$F_x\$ (\$N\$)", "\$F_y\$ (\$N\$)", "\$F_z\$ (\$N\$)",
    "\$M_x\$ (\$Nm\$)", "\$M_y\$ (\$Nm\$)", "\$M_z\$ (\$N\$)"]
#nb ph = Vector{Any}(undef, 12)

for i = 1:12
    #md local y #hide
#nb    ph[i] = plot(
#!nb    plot(
        xlim = (0, 2.0),
        xticks = 0:0.5:2.0,
        xlabel = "Time (s)",
        ylabel = ylabel[i],
        grid = false,
        overwrite_figure=false
        )
    y_gxbeam = [getproperty(state.points[point[i]], field[i])[direction[i]]
        for state in gxbeam_history]

    y_diffeq = [getproperty(state.points[point[i]], field[i])[direction[i]]
        for state in diffeq_history]

    if field[i] == :theta
        ## convert to Rodriguez parameter
        @. y_gxbeam = 4*atan(y_gxbeam/4)
        @. y_diffeq = 4*atan(y_diffeq/4)
        ## convert to degrees
        @. y_gxbeam = rad2deg(y_gxbeam)
        @. y_diffeq = rad2deg(y_diffeq)
    end

    if field[i] == :F || field[i] == :M
        y_gxbeam = -y_gxbeam
        y_diffeq = -y_diffeq
    end

    plot!(t, y_gxbeam, label="GXBeam")
    plot!(sol.t, y_diffeq, label="DifferentialEquations")
#!nb     plot!(show=true)
#md     savefig("../assets/diffeq-"*string(field[i])*string(direction[i])*".svg"); #hide
#md     closeall() #hide
end

#md end #hide
#md nothing #hide

#nb ph[1]
#nb #-
#nb ph[2]
#nb #-
#nb ph[3]
#nb #-
#nb ph[4]
#nb #-
#nb ph[5]
#nb #-
#nb ph[6]
#nb #-
#nb ph[7]
#nb #-
#nb ph[8]
#nb #-
#nb ph[9]
#nb #-
#nb ph[10]
#nb #-
#nb ph[11]
#nb #-
#nb ph[12]
#nb #-

#md # ![](../assets/diffeq-u1.svg)
#md # ![](../assets/diffeq-u2.svg)
#md # ![](../assets/diffeq-u3.svg)
#md # ![](../assets/diffeq-theta1.svg)
#md # ![](../assets/diffeq-theta2.svg)
#md # ![](../assets/diffeq-theta3.svg)
#md # ![](../assets/diffeq-F1.svg)
#md # ![](../assets/diffeq-F2.svg)
#md # ![](../assets/diffeq-F3.svg)
#md # ![](../assets/diffeq-M1.svg)
#md # ![](../assets/diffeq-M2.svg)
#md # ![](../assets/diffeq-M3.svg)

# As can be seen, the solutions provided by GXBeam and DifferentialEquations track closely
# with each other.

root_chord = 1.9000
tip_chord =  0.4540
airfoil = [ # MH-104
    1.00000000  0.00000000;
    0.99619582  0.00017047;
    0.98515158  0.00100213;
    0.96764209  0.00285474;
    0.94421447  0.00556001;
    0.91510964  0.00906779;
    0.88074158  0.01357364;
    0.84177999  0.01916802;
    0.79894110  0.02580144;
    0.75297076  0.03334313;
    0.70461763  0.04158593;
    0.65461515  0.05026338;
    0.60366461  0.05906756;
    0.55242353  0.06766426;
    0.50149950  0.07571157;
    0.45144530  0.08287416;
    0.40276150  0.08882939;
    0.35589801  0.09329359;
    0.31131449  0.09592864;
    0.26917194  0.09626763;
    0.22927064  0.09424396;
    0.19167283  0.09023579;
    0.15672257  0.08451656;
    0.12469599  0.07727756;
    0.09585870  0.06875796;
    0.07046974  0.05918984;
    0.04874337  0.04880096;
    0.03081405  0.03786904;
    0.01681379  0.02676332;
    0.00687971  0.01592385;
    0.00143518  0.00647946;
    0.00053606  0.00370956;
    0.00006572  0.00112514;
    0.00001249 -0.00046881;
    0.00023032 -0.00191488;
    0.00079945 -0.00329201;
    0.00170287 -0.00470585;
    0.00354717 -0.00688469;
    0.00592084 -0.00912202;
    0.01810144 -0.01720842;
    0.03471169 -0.02488211;
    0.05589286 -0.03226730;
    0.08132751 -0.03908459;
    0.11073805 -0.04503763;
    0.14391397 -0.04986836;
    0.18067874 -0.05338180;
    0.22089879 -0.05551392;
    0.26433734 -0.05636585;
    0.31062190 -0.05605816;
    0.35933893 -0.05472399;
    0.40999990 -0.05254383;
    0.46204424 -0.04969990;
    0.51483073 -0.04637175;
    0.56767889 -0.04264894;
    0.61998250 -0.03859653;
    0.67114514 -0.03433153;
    0.72054815 -0.02996944;
    0.76758733 -0.02560890;
    0.81168064 -0.02134397;
    0.85227225 -0.01726049;
    0.88883823 -0.01343567;
    0.92088961 -0.00993849;
    0.94797259 -0.00679919;
    0.96977487 -0.00402321;
    0.98607009 -0.00180118;
    0.99640466 -0.00044469;
    1.00000000  0.00000000;
]

sections = zeros(3, size(airfoil, 1), length(points))
for ip = 1:length(points)
    chord = root_chord * (1 - x[ip]/L) + tip_chord * x[ip]/L
    sections[1, :, ip] .= 0
    sections[2, :, ip] .= chord .* (airfoil[:,1] .- 0.5)
    sections[3, :, ip] .= chord .* airfoil[:,2]
end

write_vtk("dynamic-wind-turbine", assembly, gxbeam_history, sol.t; sections = sections)

# ![](../assets/wind-turbine-blade-simulation.gif)


text:
# License for this file: MIT
# Copyright 2018-2019, DLR Institute of System Dynamics and Control
#

"""
   module Fluid1_TestConstantPropertyLiquidWater

Date: Jan. 22, 2019

Test module to evaluate how fluid components can be modeled in Julia. The test models
in this module are for ModiaMedia ConstantPropertyLiquidWater. Since Modia does not yet support
automatic state selection, the test models of different media are slightly different,
because state=false has to be manually set at some components.

The components are modelled based on the new approach described in

- Zimmer, Bender, Pollok (2018): Robust Modeling of Directed Thermofluid Flows in Complex Networks.
  2nd Japanese Modelica Conference, Tokyo, Proceedings, pp 39-48. Download:
  https://www.modelica.org/events/modelica2018japan/conference-proceedings/modelica-final-proceedings-2018-Japan.pdf

# Main Author of this module
```

      /|      Dirk Zimmer, Martin Otter
  ___/_|___   German Aerospace Center (DLR)
 /  /  /  /   Institute of System Dynamics and Control
/__/__/__/    https://www.dlr.de/sr/en/
   | /
   |/

```

# Issues with this test
The differences between Fluid1_TestConstantPropertyLiquidWater with respect to TestFluid1_SimpleAir are
(a) the used medium model,
(b) the explicit setting of a non-state in examples TestLine2 and TestLoop.
(c) that TestLoop fails in Fluid1_TestConstantPropertyLiquidWater.

Model TestLoop fails:
- if tearing=true: Translation fails with an error (pipe1.inPlug.m_flow not defined)
- if tearing=false: Translation is successful and initialization fails with a KINSol error
  (It might be that the Jacobian is singular (= there are redundant equations))
- Various settings of removeSingularities and aliasElimination do not change the situation
"""
module Fluid1_TestConstantPropertyLiquidWater

include("Components.jl")
using .Fluid1_Components

using  Modia
using  ModiaMedia
import ModiaMedia.ModiaMath   # included via ModiaMedia, to avoid requirement to add it in the standard environment

const MediumUsedInExamples = getMedium("ConstantPropertyLiquidWater")

@model TestLine begin
    fixedSource = FixedSource_pT(p0=2e5, T0=400.0, Medium=MediumUsedInExamples)
    fixedSink   = FixedSink(p0=1.5e5)
    pipe        = ShortPipe(l=1.0, A=0.01)

    @equations begin
        connect(fixedSource.outPlug, pipe.inPlug)
        connect(pipe.outPlug, fixedSink.inPlug)
    end
end

result = simulate(TestLine, 0.02; logSimulation=true)
add!(result, "pipe.outPlug", ["T", "d"])

ModiaMath.plot(result, ["pipe.m_flow", "pipe.outPlug.p", "pipe.outPlug.r", "pipe.outPlug.T", "pipe.outPlug.d", "pipe.eta"], heading="Fluid1_TestConstantPropertyLiquidWater.jl: TestLine", figure=1)

@model TestJunction begin
    fixedSource = FixedSource_pT(p0=1.0e5, T0=400.0, Medium=MediumUsedInExamples)
    fixedSink   = FixedSink(  p0=0.8e5)
    pipe1       = ShortPipe(k=3e5, l=1.0, A=0.01)
    pipe2       = ShortPipe(       l=1.0, A=0.01)
    junction    = Junction()
    splitter    = Splitter()

    @equations begin
        connect(pipe2.outPlug      , junction.inPlugA)
        connect(pipe1.outPlug      , junction.inPlugB)
        connect(junction.outPlugC  , fixedSink.inPlug)
        connect(splitter.outPlugB  , pipe1.inPlug)
        connect(splitter.outPlugC  , pipe2.inPlug)
        connect(fixedSource.outPlug, splitter.inPlugA)
    end
end

result = simulate(TestJunction, 0.02; logSimulation=true)
ModiaMath.plot(result, [("junction.inPlugA.m_flow", "junction.inPlugB.m_flow" , "junction.outPlugC.m_flow"),
                        ("splitter.inPlugA.m_flow", "splitter.outPlugB.m_flow", "splitter.outPlugC.m_flow"),
                        ("pipe1.outPlug.p", "pipe2.outPlug.p")], heading="Fluid1_TestConstantPropertyLiquidWater.jl: TestJunction", figure=2)

@model TestLine2 begin
    fixedSource = FixedSource_pT(p0=2.0e5, T0=400.0, Medium=MediumUsedInExamples)
    fixedSink   = FixedSink(p0=1.0e5)
    pipe1       = ShortPipe(l=1.0, A=0.01, m_flow=MassFlowRate(state=false))
    pipe2       = ShortPipe(l=1.0, A=0.01)
    volume      = ClosedVolume(p0=1.5e5, T0=500.0, Medium=MediumUsedInExamples)

    @equations begin
        connect(fixedSource.outPlug, pipe1.inPlug)
        connect(pipe1.outPlug      , volume.inPlug)
        connect(volume.outPlug     , pipe2.inPlug)
        connect(pipe2.outPlug      , fixedSink.inPlug)
    end
end

result = simulate(TestLine2, 1.0, logTranslation=true, logSimulation=true, tearing=true)
add!(result, "pipe2.outPlug", ["T", "d"])
add!(result, "fixedSource.outPlug", ["T", "d"])
ModiaMath.plot(result, [("volume.inPlug.m_flow", "volume.outPlug.m_flow"),
                         "volume.medium.p",
                        ("fixedSource.outPlug.T", "volume.medium.T"),
                        ("volume.h", "volume.medium.h")],
                         heading="Fluid1_TestConstantPropertyLiquidWater.jl: TestLine2", figure=3)

@model TestCompressor begin
    fixedSource = FixedSource_pT(p0=1.0e5, T0=400, Medium=MediumUsedInExamples)
    fixedSink   = FixedSink(  p0=1.5e5)
    pipe1       = ShortPipe(l=1.0, A=0.01, m_flow=MassFlowRate(state=false))
    pipe2       = ShortPipe(l=1.0, A=0.01, m_flow=MassFlowRate(state=false))
    compressor  = Compressor(l=1.0, A=0.01)

    @equations begin
        connect(fixedSource.outPlug, pipe1.inPlug)
        connect(pipe1.outPlug      , compressor.inPlug)
        connect(compressor.outPlug , pipe2.inPlug)
        connect(pipe2.outPlug      , fixedSink.inPlug)
    end
end

result = simulate(TestCompressor, 0.02; logTranslation=true, logSimulation=true, tearing=false)
ModiaMath.plot(result, [ "compressor.inPlug.m_flow",
                        ("compressor.inPlug.p", "compressor.outPlug.p", "pipe2.outPlug.p")],
                         heading="Fluid1_TestConstantPropertyLiquidWater.jl: TestCompressor", figure=4)

@model TestLoop begin
    pipe1       = ShortPipe(       l=1.0, A=0.01, m_flow=MassFlowRate(state=false))
    pipe2       = ShortPipe(k=3e5, l=1.0, A=0.01, m_flow=MassFlowRate(state=false))
    compressor  = Compressor(      l=1.0, A=0.01, m_flow=MassFlowRate(state=false))
    volume      = ClosedVolume(p0=1e5, T0=500, Medium=MediumUsedInExamples)

    @equations begin
        connect(volume.outPlug    , pipe1.inPlug)
        connect(pipe1.outPlug     , compressor.inPlug)
        connect(compressor.outPlug, pipe2.inPlug)
        connect(pipe2.outPlug     , volume.inPlug)
    end
end

result = simulate(TestLoop, 0.02; logTranslation=true, logSimulation=true, tearing=false, removeSingularities=true, aliasElimination=true)
ModiaMath.plot(result, [ "compressor.inPlug.m_flow",
                        ("compressor.inPlug.p", "compressor.outPlug.p", "pipe2.outPlug.p", "volume.medium.p")],
                         heading="Fluid1_TestConstantPropertyLiquidWater.jl: TestLoop", figure=5)

end


text:
module MaxHelpingHandCustomSandwichLava

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/CustomSandwichLava" CustomSandwichLava(x::Integer, y::Integer,
    direction::String="CoreModeBased", speed::Number=20.0, sandwichGap::Number=160.0,
    hotSurfaceColor::String="ff8933", hotEdgeColor::String="f25e29", hotCenterColor::String="d01c01",
    coldSurfaceColor::String="33ffe7", coldEdgeColor::String="4ca2eb", coldCenterColor::String="0151d0")

const directions = String["AlwaysUp", "AlwaysDown", "CoreModeBased"]

const placements = Ahorn.PlacementDict(
    "Sandwich Lava (Customizable) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomSandwichLava
    )
)

Ahorn.editingOptions(entity::CustomSandwichLava) = Dict{String, Any}(
    "direction" => directions
)

function Ahorn.selection(entity::CustomSandwichLava)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSandwichLava, room::Maple.Room)
    direction = get(entity.data, "direction", "CoreModeBased")

    if direction == "AlwaysUp"
        Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/lava_sandwich_up", 0, 0)
    elseif direction == "AlwaysDown"
        Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/lava_sandwich_down", 0, 0)
    else
        Ahorn.drawImage(ctx, Ahorn.Assets.lavaSanwitch, -12, -12)
    end
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using Images, ImageView, ImageCore, BenchmarkTools
include("src/SLIC.jl")

input_image = load("lenna.bmp")
#input_image = load("dog.png")
#imshow(input_image)
#mosaicview(input_image)

#@btime out_image = slic($input_image, 100, 10, 10)
#save("foo.png", out_image)
@btime out_image_1 = slic($input_image, 100, 10, 10, true)
#save("bar.png", out_image_1)


text:
using SparseArrays
using LinearAlgebra
using MatrixDepot

matrices = [
    "LPnetlib/lpi_itest6",
    "HB/west0132",
    "LPnetlib/lp_etamacro",
    "LPnetlib/lp_blend",
    "Pajek/GD99_c",
    "HB/can_292",
]

open("matrices.jl", "w") do f
    println(f, "using SparseArrays")
    println(f, "using LinearAlgebra")
    println(f, "matrices = Dict{String, Any}(")
    for mtx in matrices
        A = mdopen(mtx).A
        if issymmetric(A)
            n, n = size(A)
            I, J, V = findnz(SparseMatrixCSC(A.data))
            println(f, "\"$mtx\" => Symmetric(sparse($I, $J, $V, $n, $n), Symbol(\"$(A.uplo)\")),")
        else
            m, n = size(A)
            I, J, V = findnz(SparseMatrixCSC(A))
            println(f, "\"$mtx\" => sparse($I, $J, $V, $m, $n),")
        end
    end
    println(f, ")")
end


text:
module Aggregator

# Base Mangrove types and conversions.
include("MangroveBase.jl")
using .MangroveBase: ReviewInfo, normalize, generate_data, convert, mean, hash
export ReviewInfo, normalize, generate_data

# Statistical model and inference.
include("Model.jl")
using .Model: mangrove_model, get_chains, subs, qualities
export mangrove_model, get_chains, subs, qualities

# Mangrove API access.
include("data/Api.jl")
# Mangrove DB reading and writing.
include("data/Db.jl")

# Script to be ran periodically to compute statistics and store them.
include("Job.jl")
using .Job: store_inferred
export store_inferred

export convert, hash, mean

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
# If adding models from MLJModels for testing purposes, then do the
# following in the interface file (eg, DecisionTree.jl):

# - change `import ..DecisionTree` to `import DecisionTree`
# - remove wrapping as module

module Models

using MLJBase
import MLJModelInterface: @mlj_model, metadata_model, metadata_pkg
import MLJModelInterface

include("_models/Constant.jl")
include("_models/NearestNeighbors.jl")

end


text:
using GeoPhyInv
using Test
using Random

# # Intro

#md # ```@docs
#md # Data
#md # ```

# # Examples
# Define an acquisition geometry.
mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

# Need a time grid.
tgrid=range(0, stop=1.0, step=0.1);

# Lets initialize records for `:P` and `:Vz` fields.
data=Data(tgrid, ageom, [:P,:Vz]);

# Fill the `:P` field of 3rd supersource with random numbers.
Random.randn!(data[3][:P]);

# # Methods
# Methods listed for `SrcWav` can used on instances of `Data` too.


text:
Base.show(io::IO, ::MIME"text/org", org::Org) = org(io, org)

function Base.string(component::Union{Org, OrgComponent})
    b = IOBuffer()
    org(b, component)
    String(take!(b))
end

function org(io::IO, o::Org)
    for component in o.contents
        (component isa Heading && component !== first(o.contents)) && print(io, '\n')
        org(io, component)
        print(io, '\n')
    end
end

function org(io::IO, component::OrgComponent, indent::Integer)
    print(io, ' '^indent)
    org(io, component)
end

org(o::Union{Org, OrgComponent}) = org(stdout, o)

function org(_::IO, component::C) where { C <: OrgComponent}
    @warn "No method for converting $C to a string representation currently exists"
end

# ---------------------
# Sections
# ---------------------

function org(io::IO, heading::Heading, indent::Integer=0)
    print(io, ' '^indent)
    print(io, '*'^heading.level, ' ')
    if !isnothing(heading.keyword)
        print(io, heading.keyword, ' ')
    end
    if !isnothing(heading.priority)
        print(io, "[#", heading.priority, "] ")
    end
    org.(Ref(io), heading.title)
    if length(heading.tags) > 0
        print(io, " :", join(heading.tags, ':'), ":")
    end
    if !isnothing(heading.planning)
        print(io, '\n')
        org(io, heading.planning, indent)
    end
    if !isnothing(heading.properties)
        print(io, '\n')
        org(io, heading.properties, indent)
    end
    if !isnothing(heading.section)
        print(io, "\n\n")
        org(io, heading.section, indent)
    end
end

function org(io::IO, section::Section, indent::Integer=0)
    for component in section.contents
        org(io, component, indent)
        component === last(section.contents) || print(io, "\n\n")
    end
end

function org(io::IO, afkw::AffiliatedKeywordsWrapper, indent::Integer=0)
    for afk in afkw.keywords
        org(io, afk, indent)
        print(io, '\n')
    end
    org(io, afkw.element, indent)
end

# ---------------------
# Greater Elements
# ---------------------

function org(io::IO, specialblock::SpecialBlock, indent::Integer=0)
    print(io, ' '^indent, "#+begin_", specialblock.name)
    if !isnothing(specialblock.parameters)
        print(io, ' ', specialblock.parameters)
    end
    print(io, '\n')
    for el in specialblock.contents
        org(io, el, indent)
    end
    print(io, ' '^indent, "#+end_", specialblock.name)
end

function org(io::IO, centerb::CenterBlock, indent::Integer=0)
    org(io, SpecialBlock("center", centerb.parameters, centerb.contents), indent)
end

function org(io::IO, quoteb::QuoteBlock, indent::Integer=0)
    org(io, SpecialBlock("quote", quoteb.parameters, quoteb.contents), indent)
end

function org(io::IO, drawer::Drawer, indent::Integer=0)
    print(io, ' '^indent, ':', drawer.name, ":\n")
    for node in drawer.contents
        org(io, node, indent)
        print(io, '\n')
    end
    print(io, ' '^indent, ":END:")
end

function org(io::IO, o::Org, dynblock::DynamicBlock, indent::Integer=0)
    print(io, ' '^indent, "#+begin: ", dynblock.name)
    if !isnothing(dynblock.parameters)
        print(io, ' ', dynblock.parameters)
    end
    for el in dynblock.contents
        org(io, o, el, indent)
    end
    print(io, ' '^indent, "#+end: ")
end

function org(io::IO, fn::FootnoteDefinition, indent::Integer=0)
    print(io, ' '^indent, "[fn:", fn.label, "] ")
    contentbuf = IOContext(IOBuffer(), :color => get(io, :color, false),
                           :displaysize => (displaysize(io)[1],
                                            displaysize(io)[2] - indent - 2))
    parlines = if fn.definition[1] isa Paragraph
        for obj in fn.definition[1]; org(contentbuf, obj) end
        contents = String(take!(contentbuf.io))
        wraplines(contents, displaysize(io)[2] - indent, 6 + ncodeunits(fn.label))
    else
        [""]
    end
    components = @view fn.definition[if fn.definition[1] isa Paragraph 2 else 1 end:end]
    for component in components
        org(contentbuf, component, indent)
        component === last(components) || print(contentbuf, '\n')
    end
    otherlines = split(String(take!(contentbuf.io)), '\n')
    lines = vcat(parlines, if otherlines == [""]; [] else otherlines end)
    for line in lines
        print(io, line)
        line === last(lines) || print(io, '\n', ' '^indent)
    end
end

# InlineTask

function org(io::IO, item::Item, indent::Integer=0, offset::Integer=0)
    print(io, ' '^indent, item.bullet)
    offset += length(item.bullet)
    if !isnothing(item.counterset)
        print(io, " [@", item.counterset, "]")
        offset += length(item.counterset) + 4
    end
    if !isnothing(item.checkbox)
        print(io, " [", item.checkbox, ']')
        offset += 4
    end
    if !isnothing(item.tag)
        print(io, item.tag, " ::")
        offset += length(item.tag) + 3
    end
    if length(item.contents) > 0
        print(io, ' ')
        offset += 1
        contentbuf = IOContext(IOBuffer(), :color => get(io, :color, false),
                               :displaysize => (displaysize(io)[1],
                                                displaysize(io)[2] - indent - 2))
        parlines = if item.contents[1] isa Paragraph
            for obj in item.contents[1]; org(contentbuf, obj) end
            contents = String(take!(contentbuf.io))
            wraplines(contents, displaysize(io)[2] - indent - 2, offset)
        else
            [""]
        end
        components = @view item.contents[if item.contents[1] isa Paragraph 2 else 1 end:end]
        for component in components
            org(contentbuf, component, indent)
            component === last(components) || print(contentbuf, '\n')
        end
        otherlines = split(String(take!(contentbuf.io)), '\n')
        lines = vcat(parlines, if otherlines == [""]; [] else otherlines end)
        for line in lines
            print(io, line)
            line === last(lines) || print(io, '\n', ' '^(indent+2))
        end
    end
end

function org(io::IO, list::List, indent::Integer=0)
    for item in list.items
        org(io, item, indent)
        item === last(list.items) || print(io, '\n')
    end
end

function org(io::IO, propdrawer::PropertyDrawer, indent::Integer=0)
    print(io, ' '^indent, ":PROPERTIES:\n")
    for node in propdrawer.contents
        org(io, node, indent)
        print(io, '\n')
    end
    print(io, ' '^indent, ":END:")
end

const table_charset_org =
    Dict('|' => '|',
         '>' => '|',
         '<' => '|',
         '[' => ' ',
         ']' => ' ',
         '-' => '-',
         '+' => '+')

org(io::IO, table::Table, indent::Integer=0) =
    layouttable(io, org, table, table_charset_org, indent)

# ---------------------
# Elements
# ---------------------

org(io::IO, bcall::BabelCall) = print(io, "#+call: ", bcall.name)

function org(io::IO, block::Block, indent::Integer=0)
    name, data = if block isa CommentBlock
        ("comment", nothing)
    elseif block isa VerseBlock
        ("verse", nothing)
    elseif block isa ExampleBlock
        ("example", nothing)
    elseif block isa ExportBlock
        ("example", block.backend)
    elseif block isa SourceBlock
        ("src", if isnothing(block.arguments)
             block.lang
         else
             string(block.lang, " ", block.arguments)
         end)
    end
    print(io, ' '^indent, "#+begin_", name)
    if !isnothing(data)
        print(io, ' ', data)
    end
    print(io, '\n')
    if block isa VerseBlock
        print(io, "Oh noes! A verse block...\n")
    else
        for line in block.contents
            print(io, ' '^indent)
            if startswith(line, '*')
                print(io, ',')
            end
            print(io, line, '\n')
        end
    end
    print(io, ' '^indent, "#+end_", name)
end

function org(io::IO, clock::Clock{TimestampInactive})
    print(io, "clock: ")
    org(io, clock.timestamp)
end

function org(io::IO, clock::Clock{TimestampInactiveRange})
    print(io, "clock: ")
    org(io, clock.timestamp)
    print(io, " => ", clock.duration[1], ':',
          lpad(string(clock.duration[2]), 2, '0'))
end

org(io::IO, diarysexp::DiarySexp) = print(io, "%%", diarysexp.sexp)

function org(io::IO, planning::Planning)
    values = [(type, getproperty(planning, type))
              for type in (:deadline, :scheduled, :closed)] |>
                  vals -> filter(v -> !isnothing(v[2]), vals)
    for val in values
        print(io, uppercase(string(val[1])), ' ')
        org(io, val[2])
        val === last(values) || print(io, ' ')
    end
end

function org(io::IO, comment::Comment, indent::Integer=0)
    for line in comment.contents
        print(io, ' '^indent, "# ", line)
        line === last(comment.contents) || print(io, '\n')
    end
end

function org(io::IO, fw::FixedWidth, indent::Integer=0)
    for line in fw.contents
        print(io, ' '^indent, ": ", line)
        line === last(fw.contents) || print(io, '\n')
    end
end

org(io::IO, ::HorizontalRule) = print(io, "-----")

org(io::IO, keyword::Keyword{<:AbstractString}) =
    print(io, "#+", keyword.key, ": ", keyword.value)

function org(io::IO, keyword::Keyword{Vector{Object}})
    print(io, "#+", keyword.key, ": ")
    for obj in keyword.value
        org(io, obj)
    end
end

function org(io::IO, afk::AffiliatedKeyword{<:AbstractString})
    print(io, "#+", afk.key)
    if !isnothing(afk.optval)
        print(io, '[', afk.optval, ']')
    end
    print(io, ':', afk.value)
    if !isnothing(afk.value)
        print(io, ' ', afk.value)
    end
end

function org(io::IO, afk::AffiliatedKeyword{Vector{Object}})
    print(io, "#+", afk.key)
    if !isnothing(afk.optval)
        print(io, '[')
        for obj in afk.optval
            org(io, obj)
        end
        print(io, ']')
    end
    print(io, ':')
    if !isnothing(afk.optval)
        print(io, ' ')
        for obj in afk.value
            org(io, obj)
        end
    end
end

function org(io::IO, env::LaTeXEnvironment, indent::Integer=0)
    print(io, ' '^indent, "\\begin{", env.name, '}', env.contents[1], '\n')
    for line in env.contents[2:end]
        print(io, ' '^(2+indent), line, '\n')
    end
    print(io, ' '^indent, "\\end{", env.name, '}')
end

org(io::IO, node::NodeProperty) =
    print(io, ':', node.name, if node.additive "+:" else ":" end, if !isnothing(node.value) node.value else "" end)

function org(io::IO, objs::Vector{Object})
    for obj in objs
        org(io, obj)
    end
end

function Base.string(par::Paragraph)
    b = IOBuffer()
    org(b, par.contents)
    String(take!(b))
end

function org(io::IO, par::Paragraph, indent::Integer=0)
    lines = wraplines(string(par), displaysize(io)[2] - indent)
    for line in lines
        print(io, ' '^indent, line)
        line === last(lines) || print(io, '\n')
    end
end

# Table Row

# ---------------------
# Objects
# ---------------------

org(io::IO, entity::Entity) = print(io, '\\', entity.name)

function org(io::IO, latex::LaTeXFragment)
    if isnothing(latex.delimiters)
        print(io, latex.contents)
    else
        print(io, latex.delimiters[1], latex.contents, latex.delimiters[2])
    end
end

org(io::IO, snippet::ExportSnippet) =
    print(io, "@@", snippet.backend, ':', snippet.snippet, "@@")

function org(io::IO, fn::FootnoteReference)
    print(io, "[fn:")
    if !isnothing(fn.label)
        print(io, fn.label)
    end
    if !isnothing(fn.definition)
        print(io, ':')
        contentbuf = IOContext(IOBuffer(), :color => get(io, :color, false))
        for obj in fn.definition
            org(contentbuf, obj)
        end
        print(io, String(take!(contentbuf.io)))
    end
    print(io, ']')
end

function org(io::IO, keycite::CitationReference)
    org(io, keycite.prefix)
    print(io, '@', keycite.key)
    org(io, keycite.suffix)
end

function org(io::IO, cite::Citation)
    print(io, "[cite")
    if !isnothing(cite.style[1])
        print(io, '/', cite.style[1])
    end
    if !isnothing(cite.style[2])
        print(io, '/', cite.style[2])
    end
    print(io, ':')
    if !isnothing(cite.globalprefix)
        org(io, cite.globalprefix)
        print(io, ';')
    end
    for keycite in cite.citerefs
        org(io, keycite)
        keycite === last(cite.citerefs) || print(io, ';')
    end
    if !isnothing(cite.globalsuffix)
        print(io, ';')
        org(io, cite.globalsuffix)
    end
    print(io, ']')
end

function org(io::IO, bcall::InlineBabelCall)
    print(io, "call_", bcall.name)
    if !isnothing(bcall.header)
        print(io, '[', bcall.header, ']')
    end
    print(io, '(', bcall.arguments, ')')
end

function org(io::IO, src::InlineSourceBlock)
    print(io, "src_", src.lang)
    if !isnothing(src.options)
        print(io, '[', src.options, ']')
    end
    print(io, '{', src.body, '}')
end

org(io::IO, ::LineBreak) = print(io, "\\\\\n")

function org(io::IO, link::RadioLink)
    for obj in link.radio.contents
        org(io, obj)
    end
end

const org_link_protocol_prefixes =
    Dict(:coderef => p -> "($p)",
         :custom_id => p -> "#$p",
         :heading => p -> "*$p",
         :fuzzy => identity)

const org_link_protocol_prefer_slashprefix =
    ("https", "http", "ftp")

function org(io::IO, path::LinkPath)
    if path.protocol in keys(org_link_protocol_prefixes)
        print(io, org_link_protocol_prefixes[path.protocol](path.path))
    else
        print(io, path.protocol, ':',
              if path.protocol in org_link_protocol_prefer_slashprefix
                  "//" else "" end,
              path.path)
    end
end

org(io::IO, link::PlainLink) = org(io, link.path)

function org(io::IO, link::AngleLink)
    print(io, '<')
    org(io, link.path)
    print(io, '>')
end

function org(io::IO, link::RegularLink)
    print(io, "[[")
    org(io, link.path)
    if isnothing(link.description)
        print(io, "]]")
    else
        print(io, "][")
        for obj in link.description
            org(io, obj)
        end
        print(io, "]]")
    end
end

org(io::IO, mac::Macro) = print(io, "{{{", mac.name, '(', join(mac.arguments, ","), ")}}}")

function org(io::IO, radio::RadioTarget)
    print(io, "<<<")
    org.(Ref(io), radio.contents)
    print(io, ">>>")
end

org(io::IO, target::Target) = print(io, "<<", target.target, ">>")

org(io::IO, statscookie::StatisticsCookiePercent) =
    print(io, '[', statscookie.percentage, "%]")
org(io::IO, statscookie::StatisticsCookieFraction) =
    print(io, '[', if !isnothing(statscookie.complete) string(statscookie.complete) else "" end,
          '/', if !isnothing(statscookie.total) string(statscookie.total) else "" end, ']')

org(io::IO, script::Superscript) = print(io, script.char, '^', script.script)
org(io::IO, script::Subscript) = print(io, script.char, '_', script.script)

org(io::IO, cell::TableCell) = print(io, "| ", cell.contents, " |")

const org_timestamp_repeaterordelay_marks =
    Dict(:cumulative => "+",
         :catchup => "++",
         :restart => ".+",
         :warningall => "-",
         :warningfirst => "--")

org(io::IO, tsrod::TimestampRepeaterOrDelay) =
    print(io, org_timestamp_repeaterordelay_marks[tsrod.type], tsrod.value, tsrod.unit)

org(io::IO, tsd::TimestampDiary) = print(io, "<%%", tsd.sexp, '>')

function org(io::IO, ts::TimestampInstant)
    bra, ket = if ts isa TimestampActive; ('<', '>') else ('[', ']') end
    print(io, bra)
    print(io, ts.date, ' ', dayabbr(ts.date))
    if !isnothing(ts.time)
        print(io, ' ', hour(ts.time), ':', minute(ts.time))
    end
    if !isnothing(ts.repeater)
        print(io, ' ')
        org(io, ts.repeater)
    end
    if !isnothing(ts.warning)
        print(io, ' ')
        org(io, ts.warning)
    end
    print(io, ket)
end

function org(io::IO, tsr::TimestampRange)
    if tsr.start.date == tsr.stop.date &&
        tsr.start.repeater == tsr.stop.repeater
        bra, ket = if tsr isa TimestampActiveRange; ('<', '>') else ('[', ']') end
        print(io, bra)
        print(io, tsr.start.date, ' ', dayabbr(tsr.start.date))
        print(io, ' ', hour(tsr.start.time), ':', minute(tsr.start.time),
              '-', hour(tsr.stop.time), ':', minute(tsr.stop.time))
        if !isnothing(tsr.start.repeater)
            print(io, ' ')
            org(io, tsr.start.repeater)
        end
        if !isnothing(tsr.start.warning)
            print(io, ' ')
            org(io, tsr.start.warning)
        end
        print(io, ket)
    else
        org(io, tsr.start)
        print(io, "--")
        org(io, tsr.stop)
    end
end

org(io::IO, text::TextPlain) = print(io, text.text)

const TextMarkupMarkers =
    Dict(:bold => '*',
         :italic => '/',
         :underline => '_',
         :verbatim => '=',
         :code => '~',
         :strikethrough => '+')

function org(io::IO, markup::TextMarkup)
    print(io, TextMarkupMarkers[markup.formatting])
    if markup.contents isa AbstractString
        print(io, markup.contents)
    else
        org.(Ref(io), markup.contents)
    end
    print(io, TextMarkupMarkers[markup.formatting])
end


text:
UV = UncertainValueDataset(example_uvals)

n = 3
@test resample(UV) isa Vector{T} where T <: Real
@test resample(UV, n) isa Vector{Vector{T}} where T <: Real


text:
using RobustAndOptimalControl, ControlSystems
using RobustAndOptimalControl: @check_unique, @check_all_unique

@test :x^3 == expand_symbol(:x, 3) == [:x1, :x2, :x3]

G1 = ss(1,1,1,1)
G2 = ss(1,1,1,1)
s1 = named_ss(G1, x = [:x], u = :u, y=:y) # test providing only symbol
s2 = named_ss(G2, x = [:z], u = [:u], y=[:y])

@test s1[:y, :u] == s1
@test s1[[:y], [:u]] == s1

G3 = ControlSystems.ssrand(1,2,3)
s3 = named_ss(G3, x = [:x1, :x2, :x3], u = [:u1, :u2], y=[:y])

@test s3[:y, :u1] == named_ss(G3[1,1], x = [:x1, :x2, :x3], u = [:u1], y=[:y])
@test s3[:y, :u2] == named_ss(G3[1,2], x = [:x1, :x2, :x3], u = [:u2], y=[:y])
@test s3[1, 2] == named_ss(G3[1,2], x = [:x1, :x2, :x3], u = [:u2], y=[:y]) # indexing with integer preserves names

@test s3[:y, [:u1, :u2]] == s3

G3 = partition(ssrand(2,3,4), 1, 1)
s3 = named_ss(G3, x = :x, u = :u, y=:y, z=:z, w=:w)
@test s3[:z, :w].sys == G3[1,1]
@test length(s3.u) == 3
@test length(s3.y) == 2

for op in (+, -)
    @testset "$op" begin
        @info "Testing $op"

        G12 = op(G1, G2)
        s12 = op(s1, s2)

        @test s12.x == [:x, :z]
        @test s12.u == [:u]
        @test s12.y == [:y]

        @test s12.A == G12.A
        @test s12.B == G12.B
        @test s12.C == G12.C
        @test s12.D == G12.D
    end
end

s1 = named_ss(G1, x = [:x], u = [:u1], y=[:y1])
s2 = named_ss(G2, x = [:z], u = [:u2], y=[:y2])

@testset "*" begin
    G12 = *(G1, G2)
    s12 = *(s1, s2)

    @test s12.x == [:x, :z]
    @test s12.u == [:u2]
    @test s12.y == [:y1]

    @test s12.A == G12.A
    @test s12.B == G12.B
    @test s12.C == G12.C
    @test s12.D == G12.D

    @test_throws ArgumentError s1*s1
end

@testset "Concatenation" begin
    @info "Testing Concatenation"

    # same output name
    s1 = named_ss(G1, x = :x, u = :u1, y=:y1)
    s2 = named_ss(G2, x = :z, u = :u2, y=:y1)
    s3 = named_ss(G2, x = :w, u = :u3, y=:y1)
    s1s2 = [s1 s2]
    @test s1s2.sys == [G1 G2]

    @test s1s2.x == [:x, :z]
    @test s1s2.u == [:u1, :u2]
    @test s1s2.y == [:y1]

    @test_throws ArgumentError [s1 s2 s1 s2]
    @test_nowarn [s1 s2 s3]

    # same input name
    s1 = named_ss(G1, x = :x, u = :u1, y=:y1)
    s2 = named_ss(G2, x = :z, u = :u1, y=:y2)
    s3 = named_ss(G2, x = :w, u = :u1, y=:y3)
    s1s2 = [s1; s2]
    @test s1s2.sys == [G1; G2]

    @test s1s2.x == [:x, :z]
    @test s1s2.u == [:u1]
    @test s1s2.y == [:y1, :y2]

    @test_throws ArgumentError [s1; s2; s1; s2]
    @test_nowarn [s1; s2; s3]
end

G1 = ss(1,1,1,0)
G2 = ss(1,1,1,0)
s1 = named_ss(G1, x = :x, u = :u1, y=:y1)
s2 = named_ss(G2, x = :z, u = :u2, y=:y2)
@test_nowarn @check_all_unique s1 s2

s1d = c2d(s1, 1.0)
@test s1d.sys == c2d(s1.sys, 1.0)
@test s1d.x == s1.x
@test s1d.u == s1.u
@test s1d.y == s1.y
@test_nowarn plot(step(s1d))

fb = feedback(s1, s2)
# @test fb.x == [:x, :z]
@test fb.u == s1.u
@test fb.y == [:y1]

s1 = named_ss(G1, x = [:x], u = [:u1], y=[:y1])
s2 = named_ss(G2, x = [:z], u = [:u1], y=[:y2])
@test_throws ArgumentError @check_all_unique s1 s2

## Promotion and conversion
@testset "Promotion and conversion" begin
    @info "Testing Promotion and conversion"
    # Named and not named
    G1 = named_ss(ssrand(1,1,1))
    G2 = HeteroStateSpace(ssrand(1,1,2))
    @test promote_type(typeof.((G1, G2))...) == NamedStateSpace{Continuous, HeteroStateSpace{Continuous, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}}

    G1p, G2p = promote(G1, G2)
    @test G1p isa NamedStateSpace{Continuous, HeteroStateSpace{Continuous, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}}
    @test G2p isa NamedStateSpace{Continuous, HeteroStateSpace{Continuous, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}}

    # Named and named
    @test promote_type(typeof.((G1, G2p))...) == NamedStateSpace{Continuous, HeteroStateSpace{Continuous, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}}
    G1p, G2p2 = promote(G1, G2p)
    @test G1p isa NamedStateSpace{Continuous, HeteroStateSpace{Continuous, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}}
    @test G2p == G2p
end

# P = ssrand(1,2,3)
# addP1 = named_ss(ss([1  1], P.timeevol), u=[:vf_a, :yL], y=:uP)
# addL1 = named_ss(ss([I(P.nx) -I(P.nx)], P.timeevol), u=[:xr^P.nx; :xh^P.nx], y=:x_diff)
# addP = sumblock(:(uP = vf_a + yL))
# addL = sumblock(:(x_diff = xr - xh); n=P.nx)
# @test addP1 == addP
# @test addL1 == addL

# P = ssrand(1,2,3; Ts=1)
# addP1 = named_ss(ss([1  1], P.timeevol), u=[:vf_a, :yL], y=:uP)
# addL1 = named_ss(ss([I(P.nx) -I(P.nx)], P.timeevol), u=[:xr^P.nx; :xh^P.nx], y=:x_diff)
# addP = sumblock(:(uP = vf_a + yL); P.Ts)
# addL = sumblock(:(x_diff = xr - xh); n=P.nx, P.Ts)
# @test addP1 == addP
# @test addL1 == addL

P = ssrand(1,2,3)
addP1 = named_ss(ss([1  1], P.timeevol), u=[:vf_a, :yL], y=:uP)
addL1 = named_ss(ss([I(P.nx) -I(P.nx)], P.timeevol), u=[:xr^P.nx; :xh^P.nx], y=:x_diff)
addP = sumblock("uP = vf_a + yL")
addL = sumblock("x_diff = xr - xh"; n=P.nx)
@test addP1 == addP
@test addL1 == addL

P = ssrand(1,2,3; Ts=1)
addP1 = named_ss(ss([1  1], P.timeevol), u=[:vf_a, :yL], y=:uP)
addL1 = named_ss(ss([I(P.nx) -I(P.nx)], P.timeevol), u=[:xr^P.nx; :xh^P.nx], y=:x_diff)
addP = sumblock("uP = vf_a + yL"; P.Ts)
addL = sumblock("x_diff = xr - xh"; n=P.nx, P.Ts)
@test addP1 == addP
@test addL1 == addL

# Without spaces
addP = sumblock("uP=vf_a+yL"; P.Ts)
addL = sumblock("x_diff=xr-xh"; n=P.nx, P.Ts)
@test addP1 == addP
@test addL1 == addL

# More spaces
addP = sumblock("uP  =  vf_a  +  yL"; P.Ts)
addL = sumblock("x_diff  =  xr  -  xh"; n=P.nx, P.Ts)
@test addP1 == addP
@test addL1 == addL

# many
addmany = sumblock("uP  =  vf_a  +  yL + v"; P.Ts)
@test addmany.y == [:uP]
@test addmany.u == [:vf_a, :yL, :v]
@test addmany.D == [1 1 1]

# many
addmany = sumblock("uP  =  vf_a  +  yL - v"; P.Ts, n=2)
@test addmany.y == :uP^2
@test addmany.u == [:vf_a^2; :yL^2; :v^2]
@test addmany.D == [I(2) I(2) -I(2)]

##

G1 = ss(1,[1 1],1,0)
G2 = ss(1,1,1,0)
s1 = named_ss(G1, x = [:x], u = :u, y=:y1) # test providing only symbol
s2 = named_ss(G2, x = [:z], u = [:u], y=[:y2])

connections = [:y1 => :u, :y2 => :u1]
w1 = [:u2]
G = connect([s1,s2], connections; w1)
@test G.u == w1
@test G.y == :y^2
@test G.sys == ss(ones(2,2), [1,0], I(2), 0)

sp = splitter(:u2, 2)
@test all(sp.D .== [1, 1])
@test sp.u == [:u2]
@test sp.y == :u2^2

@testset "complicated_feedback" begin
    @info "Testing complicated_feedback"
    include("../examples/complicated_feedback.jl")
end


text:
##### Entrainment-Detrainment models

abstract type EntrDetrModel end

struct BOverW2{FT} <: EntrDetrModel
    ε_factor::FT
    δ_factor::FT
end

struct MoistureDeficit{FT} <: EntrDetrModel
    entr_factor::FT
    detr_factor::FT
    detr_RH_power::FT
    sigmoid_slope_param::FT
    upd_mixing_frac::FT
    turb_entr_fac::FT
    entr_tke_fac::FT
end

"""
    compute_entrainment_detrainment!

Define entrainment and detrainment fields

 - `aux[:ε_model, k, i]`
 - `aux[:δ_model, k, i]`

 for all `k` and all `i`
"""
function compute_entrainment_detrainment! end

function compute_entrainment_detrainment!(
    grid::Grid{FT},
    UpdVar,
    aux,
    q,
    params,
    model::BOverW2,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    Δzi = grid.Δzi
    k_1 = first_interior(grid, Zmin())
    @inbounds for i in ud
        zi = UpdVar[i].cloud.base
        @inbounds for k in over_elems_real(grid)
            buoy = aux[:buoy, k, i]
            w = q[:w, k, i]
            if grid.zc[k] >= zi
                detr_sc = 4.0e-3 + 0.12 * abs(min(buoy, 0.0)) / max(w * w, 1e-2)
            else
                detr_sc = FT(0)
            end
            entr_sc = 0.12 * max(buoy, FT(0)) / max(w * w, 1e-2)
            aux[:ε_model, k, i] = entr_sc * model.ε_factor
            aux[:δ_model, k, i] = detr_sc * model.δ_factor
        end
        aux[:ε_model, k_1, i] = 2 * Δzi
        aux[:δ_model, k_1, i] = FT(0)
    end
end

function compute_entrainment_detrainment!(
    grid::Grid{FT},
    UpdVar,
    aux,
    q,
    params,
    model::MoistureDeficit,
) where {FT}
    @unpack param_set = params
    gm, en, ud, sd, al = allcombinations(q)
    Δzi = grid.Δzi
    k_1 = first_interior(grid, Zmin())
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            # get parameters
            c_ε = model.entr_factor
            c_turb = model.turb_entr_fac
            if aux[:q_liq, k, en] + aux[:q_liq, k, i] > 0.0
                c_δ = model.detr_factor
            else
                c_δ = 0.0
            end
            β = model.detr_RH_power
            μ_0 = model.sigmoid_slope_param
            χ = model.upd_mixing_frac
            c_λ = model.entr_tke_fac

            # get subdomain properties
            b_up = aux[:buoy, k, i]
            b_en = aux[:buoy, k, en]
            w_up = max(q[:w, k, i], 1e-4)
            w_en = q[:w, k, en]
            dw = max(w_up - w_en, 1e-4)
            db = b_up - b_en
            sqrt_tke = sqrt(max(q[:tke, k, en], 0.0))
            ts = ActiveThermoState(param_set, q, aux, k, en)
            RH_en = relative_humidity(ts)
            ts = ActiveThermoState(param_set, q, aux, k, i)
            RH_up = relative_humidity(ts)

            # compute aux functions
            D_ϵ =
                1 / (
                    1 + exp(
                        -db / dw / μ_0 *
                        (χ - q[:a, k, i] / (q[:a, k, i] + q[:a, k, en])),
                    )
                )
            D_δ =
                1 / (
                    1 + exp(
                        db / dw / μ_0 *
                        (χ - q[:a, k, i] / (q[:a, k, i] + q[:a, k, en])),
                    )
                )
            M_δ = (max((RH_up^β - RH_en^β), 0.0))^(1 / β)
            M_ϵ = (max((RH_en^β - RH_up^β), 0.0))^(1 / β)
            λ = min(abs(db / dw), c_λ * abs(db / (sqrt_tke + 1e-8)))

            # compute entrainment/detrainmnet components
            aux[:εt_model, k, i] =
                2 * q[:a, k, i] * c_turb * sqrt_tke /
                (w_up * q[:a, k, i] * UpdVar[i].cloud.updraft_top)
            aux[:ε_model, k, i] = λ / w_up * (c_ε * D_ϵ + c_δ * M_ϵ)
            aux[:δ_model, k, i] = λ / w_up * (c_ε * D_δ + c_δ * M_δ)
        end
        aux[:εt_model, k_1, i] = FT(0)
        aux[:ε_model, k_1, i] = 2 * Δzi
        aux[:δ_model, k_1, i] = FT(0)
    end
end

function compute_cv_entr!(
    grid::Grid{FT},
    q,
    aux,
    aux_O2,
    ϕ,
    ψ,
    cv,
    tke_factor,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems_real(grid)
        aux_O2[cv][:entr_gain, k] = FT(0)
        @inbounds for i in ud
            Δϕ = q[ϕ, k, i] - q[ϕ, k, en]
            Δψ = q[ψ, k, i] - q[ψ, k, en]
            aux_O2[cv][:entr_gain, k] +=
                tke_factor *
                q[:a, k, i] *
                abs(q[:w, k, i]) *
                aux[:δ_model, k, i] *
                Δϕ *
                Δψ
        end
        aux_O2[cv][:entr_gain, k] *= aux[:ρ_0, k]
    end
end


text:
using Pkg

# Set Julia environment path
juliaenv = ".julia"

# Activate environment
Pkg.activate(juliaenv)

# Set Julia development dir
ENV["JULIA_PKG_DEVDIR"] = abspath(juliaenv)

# Add Spine package registry
pkg"registry add https://github.com/Spine-project/SpineJuliaRegistry.git"

# Add packages from Project.toml
Pkg.instantiate()

### Other examples #############################################################
## Download the source for development (will go in `JULIA_PKG_DEVDIR`)
# Pkg.develop(PackageSpec(url="https://github.com/Spine-project/SpineOpt.jl"))
#
## Specify commit hashes
# Pkg.add(PackageSpec(url="https://github.com/Spine-project/SpineInterface.jl.git",
#                     rev="ebe87b8c3ccbf046d44cf30b66a9ad8830b0e477"))
# Pkg.add(PackageSpec(url="https://github.com/Spine-project/SpineOpt.jl",
#                     rev="71ff747f64ae9d318a18f223c67c6b409928cf34"))
################################################################################