text:
# Julia wrapper for header: /usr/include/scip/cons_nonlinear.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrNonlinear(scip)
    ccall((:SCIPincludeConshdlrNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPincludeNonlinconsUpgrade(scip, nonlinconsupgd, nodereform, priority, active, conshdlrname)
    ccall((:SCIPincludeNonlinconsUpgrade, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Cvoid}, Ptr{Cvoid}, Cint, UInt32, Cstring), scip, nonlinconsupgd, nodereform, priority, active, conshdlrname)
end

function SCIPcreateConsNonlinear(scip, cons, name, nlinvars, linvars, lincoefs, nexprtrees, exprtrees, nonlincoefs, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, Ptr{Ptr{SCIP_EXPRTREE}}, Ptr{Cdouble}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nlinvars, linvars, lincoefs, nexprtrees, exprtrees, nonlincoefs, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicNonlinear(scip, cons, name, nlinvars, linvars, lincoefs, nexprtrees, exprtrees, nonlincoefs, lhs, rhs)
    ccall((:SCIPcreateConsBasicNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cint, Ptr{Ptr{SCIP_EXPRTREE}}, Ptr{Cdouble}, Cdouble, Cdouble), scip, cons, name, nlinvars, linvars, lincoefs, nexprtrees, exprtrees, nonlincoefs, lhs, rhs)
end

function SCIPcreateConsNonlinear2(scip, cons, name, nlinvars, linvars, lincoefs, exprgraphnode, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsNonlinear2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{SCIP_EXPRGRAPHNODE}, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nlinvars, linvars, lincoefs, exprgraphnode, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicNonlinear2(scip, cons, name, nlinvars, linvars, lincoefs, exprgraphnode, lhs, rhs)
    ccall((:SCIPcreateConsBasicNonlinear2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Ptr{SCIP_EXPRGRAPHNODE}, Cdouble, Cdouble), scip, cons, name, nlinvars, linvars, lincoefs, exprgraphnode, lhs, rhs)
end

function SCIPaddLinearVarNonlinear(scip, cons, var, coef)
    ccall((:SCIPaddLinearVarNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, coef)
end

function SCIPsetExprtreesNonlinear(scip, cons, nexprtrees, exprtrees, coefs)
    ccall((:SCIPsetExprtreesNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cint, Ptr{Ptr{SCIP_EXPRTREE}}, Ptr{Cdouble}), scip, cons, nexprtrees, exprtrees, coefs)
end

function SCIPaddExprtreesNonlinear(scip, cons, nexprtrees, exprtrees, coefs)
    ccall((:SCIPaddExprtreesNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cint, Ptr{Ptr{SCIP_EXPRTREE}}, Ptr{Cdouble}), scip, cons, nexprtrees, exprtrees, coefs)
end

function SCIPgetNlRowNonlinear(scip, cons, nlrow)
    ccall((:SCIPgetNlRowNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_NLROW}}), scip, cons, nlrow)
end

function SCIPgetNLinearVarsNonlinear(scip, cons)
    ccall((:SCIPgetNLinearVarsNonlinear, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinearVarsNonlinear(scip, cons)
    ccall((:SCIPgetLinearVarsNonlinear, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinearCoefsNonlinear(scip, cons)
    ccall((:SCIPgetLinearCoefsNonlinear, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetNExprtreesNonlinear(scip, cons)
    ccall((:SCIPgetNExprtreesNonlinear, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetExprtreesNonlinear(scip, cons)
    ccall((:SCIPgetExprtreesNonlinear, libscip), Ptr{Ptr{SCIP_EXPRTREE}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetExprtreeCoefsNonlinear(scip, cons)
    ccall((:SCIPgetExprtreeCoefsNonlinear, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetExprgraphNodeNonlinear(scip, cons)
    ccall((:SCIPgetExprgraphNodeNonlinear, libscip), Ptr{SCIP_EXPRGRAPHNODE}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLhsNonlinear(scip, cons)
    ccall((:SCIPgetLhsNonlinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRhsNonlinear(scip, cons)
    ccall((:SCIPgetRhsNonlinear, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPcheckCurvatureNonlinear(scip, cons)
    ccall((:SCIPcheckCurvatureNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetCurvatureNonlinear(scip, cons, checkcurv, curvature)
    ccall((:SCIPgetCurvatureNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32, Ptr{SCIP_EXPRCURV}), scip, cons, checkcurv, curvature)
end

function SCIPgetExprtreeCurvaturesNonlinear(scip, cons, checkcurv, curvatures)
    ccall((:SCIPgetExprtreeCurvaturesNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32, Ptr{Ptr{SCIP_EXPRCURV}}), scip, cons, checkcurv, curvatures)
end

function SCIPgetViolationNonlinear(scip, cons, sol, violation)
    ccall((:SCIPgetViolationNonlinear, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{Cdouble}), scip, cons, sol, violation)
end

function SCIPgetLinvarMayDecreaseNonlinear(scip, cons)
    ccall((:SCIPgetLinvarMayDecreaseNonlinear, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinvarMayIncreaseNonlinear(scip, cons)
    ccall((:SCIPgetLinvarMayIncreaseNonlinear, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetExprgraphNonlinear(scip, conshdlr)
    ccall((:SCIPgetExprgraphNonlinear, libscip), Ptr{SCIP_EXPRGRAPH}, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}), scip, conshdlr)
end

function SCIPcomputeHyperplaneThreePoints(scip, a1, a2, a3, b1, b2, b3, c1, c2, c3, alpha, beta, gamma_, delta)
    ccall((:SCIPcomputeHyperplaneThreePoints, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), scip, a1, a2, a3, b1, b2, b3, c1, c2, c3, alpha, beta, gamma_, delta)
end


text:
# From https://github.com/JuliaIntervals/TaylorModels.jl/issues/127
@testset "Taylor model" begin
    using ReachableControls.TaylorModels

    n(z, a) = a * z
    d(z, a) = 1 + a * z
    G(z, a) = n(z, a) / d(z, a)

    dom_a = 0.5 ± 0.125 # domain for a
    dom_ω = 1 ± (1 / 16)  # domain for ω

    set_variables("δa δω", order = 6)
    aN = TaylorModelN(1, 3, Interval(mid(dom_a)) × Interval(mid(dom_ω)), dom_a × dom_ω)
    ωN = TaylorModelN(2, 3, Interval(mid(dom_a)) × Interval(mid(dom_ω)), dom_a × dom_ω)

    gNr, gNi = G(as_real(ωN)*im, as_real(aN))

    gNi(IntervalBox(dom_a - mid(dom_a), dom_ω - mid(dom_ω)))
end

ext: jl


text:
# adapted from Winston.jl

# TODO: implement M,k,... prefixes; see http://jkorpela.fi/c/eng.html

import Base.range

function _magform(x)
    # Given x, returns (a,b), where x = a*10^b [a >= 1., b integral].
    if x == 0
        return 0., 0
    end
    a, b = modf(log10(abs(x)))
    a, b = 10^a, Int(b)
    if a < 1.
        a, b = a * 10, b - 1
    end
    if x < 0.
        a = -a
    end
    return a, b
end

grisu(a,b,c) = ((w,x,y) = Base.Grisu.grisu(a,b,c); (y,Base.Grisu.DIGITS[1:w],x))

function _format_ticklabel(x, range=0.; min_pow10=4)
    if x == 0
        return "0"
    end
    neg, digits, b = grisu(x, Base.Grisu.SHORTEST, Int32(0))
    if length(digits) > 5
        neg, digits, b = grisu(x, Base.Grisu.PRECISION, Int32(6))
        n = length(digits)
        while digits[n] == UInt32('0')
            n -= 1
        end
        digits = digits[1:n]
    end
    b -= 1
    if abs(b) >= min_pow10
        s = IOBuffer()
        if neg write(s, '-') end
        if digits != [0x31]
            write(s, Char(digits[1]))
            if length(digits) > 1
                write(s, '.')
                for i = 2:length(digits)
                    write(s, Char(digits[i]))
                end
            end
            #write(s, "\\times ")
        end
        #write(s, "10^{")
        #write(s, dec(b))
        #write(s, '}')
        write(s, "e")
        write(s, string(b))
        return String(take!(s))
    end
    # XXX: @sprint doesn't implement %.*f
    #if range < 1e-6
    #    a, b = _magform(range)
    #    return @sprintf "%.*f" (abs(b),x)
    #end
    s = sprint(show, x; context = :compact => true)
    endswith(s, ".0") ? s[1:end-2] : s
end

range(a::Real, b::Real) = (a <= b) ? (ceil(Int, a):floor(Int, b)) :
                                     (floor(Int, a):-1:ceil(Int, b))

function _ticklist_linear(lo, hi, sep, origin=0.)
    a = (lo - origin)/sep
    b = (hi - origin)/sep
    [ origin + i*sep for i in range(a,b) ]
end

function _ticks_default_linear(lim)
    a, b = _magform(abs(lim[2] - lim[1])/5.)
    if a < (1 + 2)/2.
        x = 1
    elseif a < (2 + 5)/2.
        x = 2
    elseif a < (5 + 10)/2.
        x = 5
    else
        x = 10
    end

    major_div = x * 10.0^b
    return _ticklist_linear(lim[1], lim[2], major_div)
end

function _ticks_default_log(lim)
    a = log10(lim[1])
    b = log10(lim[2])
    r = range(a, b)
    nn = length(r)

    if nn >= 10
        return 10.0 .^ _ticks_default_linear((a,b))
    elseif nn >= 2
        return 10.0 .^ r
    else
        return _ticks_default_linear(lim)
    end
end

_ticks_num_linear(lim, num) = linspace(lim[1], lim[2], num)
_ticks_num_log(lim, num) = logspace(log10(lim[1]), log10(lim[2]), num)

function _subticks_linear(lim, ticks, num=nothing)
    major_div = abs(ticks[end] - ticks[1])/float(length(ticks) - 1)
    if num === nothing
        _num = 4
        a, b = _magform(major_div)
        if 1. < a < (2 + 5)/2.
            _num = 3
        end
    else
        _num = num
    end
    minor_div = major_div/float(_num+1)
    return _ticklist_linear(lim[1], lim[2], minor_div, ticks[1])
end

function _subticks_log(lim, ticks, num=nothing)
    a = log10(lim[1])
    b = log10(lim[2])
    r = range(a, b)
    nn = length(r)

    if nn >= 10
        return 10.0 .^ _subticks_linear((a,b), map(log10,ticks), num)
    elseif nn >= 2
        minor_ticks = Float64[]
        for i in (minimum(r)-1):maximum(r)
            for j in 1:9
                z = j * 10.0^i
                if (lim[1] <= z <= lim[2]) || (lim[1] >= z >= lim[2])
                    push!(minor_ticks, z)
                end
            end
        end
        return minor_ticks
    else
        return _subticks_linear(lim, ticks, num)
    end
end


text:
# ---
# title: 367. Valid Perfect Square
# id: problem367
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/valid-perfect-square/description/>
# hidden: true
# ---
#
# Given a **positive** integer _num_ , write a function which returns True if
# _num_ is a perfect square else False.
#
# **Follow up:** **Do not** use any built-in library function such as `sqrt`.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 16
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: num = 14
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= num <= 2^31 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
############################################################################################
"Abstract super type for sorting data in different dimensions."
abstract type DataSorting end
"Data is sorted by rows."
struct ByRows <: DataSorting end
"Data is sorted by columns."
struct ByCols <: DataSorting end

"""
$(TYPEDEF)
Array configuration struct.

Checks if data is propagated row or column wise.

# Fields
$(TYPEDFIELDS)
"""
struct ArrayConfig{K,T<:DataSorting}
    "Dimension of data."
    size::NTuple{K,Int64}
    "Boolean if data is sorted by row"
    sorted::T
    function ArrayConfig(arraysize::NTuple{K,Int64}) where {K}
        ArgCheck.@argcheck length(arraysize) < 3 "Array dimension only allowed to be in 1 or 2 dimensions."
        ## Check if data is per row or column
        if length(arraysize) > 1 && arraysize[1] < arraysize[2]
            sorted = ByCols()
        else
            sorted = ByRows()
        end
        return new{K,typeof(sorted)}(arraysize, sorted)
    end
end
ArrayConfig(array::AbstractArray{T}) where {T} = ArrayConfig(size(array))

############################################################################################
"Abstract super type to define valid data structures."
abstract type DataStructure end

"""
$(TYPEDEF)
Determines whether data can be subsampled for inference.

# Fields
$(TYPEDFIELDS)
"""
struct SubSampled <: DataStructure
    "Buffer index for subsampled indices."
    buffer::Vector{Int64}
end
SubSampled(n::Integer) = SubSampled(repeat(1:n, 1))

"""
$(TYPEDEF)
Determines if data size for sampler is increasing over time.

# Fields
$(TYPEDFIELDS)
"""
struct Expanding <: DataStructure
    "maximal visible index"
    index::Iterator
end
Expanding(t::Integer) = Expanding(Iterator(t))

"""
$(TYPEDEF)
Determines if data size for sampler is rolled forward over time.

# Fields
$(TYPEDFIELDS)
"""
struct Rolling <: DataStructure
    "Rolling window."
    length::Int64
    "maximal visible index"
    index::Iterator
end
Rolling(t::Integer) = Rolling(t, Iterator(t))

"""
$(TYPEDEF)
Determines whether only full data can be used.

# Fields
$(TYPEDFIELDS)
"""
struct Batch <: DataStructure end

############################################################################################
"""
$(TYPEDEF)
Data tuning struct.

Contains information about data structure.

# Fields
$(TYPEDFIELDS)
"""
struct DataTune{D<:DataStructure,C<:Union{Nothing, ArrayConfig},M}
    "Subtype of DataStructure"
    structure::D
    "Data dimension configuration."
    config::C
    "Treatment of missing data."
    miss::M
    function DataTune(
        structure::D, config::C, miss::M
    ) where {D<:DataStructure,C<:Union{Nothing, ArrayConfig}, M}
        return new{D,C,M}(structure, config, miss)
    end
end

function DataTune(data, structure::Batch, miss=nothing)
    return DataTune(structure, nothing, miss)
end
function DataTune(data::AbstractArray{T}, structure::D, miss=nothing
) where {T, D<:Union{SubSampled, Expanding, Rolling}}
    return DataTune(structure, ArrayConfig(data), miss)
end

############################################################################################
"""
$(SIGNATURES)
Update datetune struct based on datastructure.

# Examples
```julia
```

"""
function update!(_rng::Random.AbstractRNG, tune::DataTune, structure::Batch)
    return nothing
end
function update!(_rng::Random.AbstractRNG, tune::DataTune, structure::Expanding)
    update!(tune.structure.index)
    return nothing
end
function update!(_rng::Random.AbstractRNG, tune::DataTune, structure::Rolling)
    update!(tune.structure.index)
    return nothing
end

function update!(_rng::Random.AbstractRNG, tune::DataTune, structure::SubSampled)
    idx = 1:maximum(tune.config.size)
    @inbounds for iter in eachindex(tune.structure.buffer)
        tune.structure.buffer[iter] = rand(_rng, idx)
    end
    return nothing
end
update!(_rng::Random.AbstractRNG, tune::DataTune) = update!(_rng, tune, tune.structure)

############################################################################################
"""
$(SIGNATURES)
Adjust data with current datatune configuration.

# Examples
```julia
```

"""
function adjust(tune::DataTune{<:Batch}, data)
    return data
end
function adjust(tune::DataTune{<:Expanding}, data)
    return grab(data, 1:(tune.structure.index.current), tune.config)
end
function adjust(tune::DataTune{<:Rolling}, data)
    return grab(
        data,
        (tune.structure.index.current - tune.structure.length + 1):(tune.structure.index.current),
        tune.config,
    )
end
function adjust(tune::DataTune{<:SubSampled}, data)
    return grab(data, tune.structure.buffer, tune.config)
end

############################################################################################
#export
export ArrayConfig, DataTune, update!, DataStructure, SubSampled, Expanding, Rolling, Batch


text:
code_0 = 0
code_1 = 1
code_a = "a"
code_b = "b"
l_b = Set([code_a])
l_c = Dict(code_b => code_0)
function main()
    @assert("a" in l_b)
    println(join(["OK"], " "))
end

main()


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Frames (ModiaMath/Frames/_module.jl)
#
const seps = sqrt(eps())

"""
    path = ModiaMath.Path(r::Vector{Vector3D},
                          q::Vector{Quaternion} = Quaternion[];
                          v = ones(length(r)))

Return an instance of a new Path object. The Path object consists of n frames
defined by the position vectors of their origins (`r[i]` for frame `i`)
and optionally of their absolute rotation quaternions (`q[i]` for frame `i`)
describing the rotation from the world frame to the respective frame.

A path parameter `t` is defined in the following way on these frames:

- `t[1] = 0.0`.
- `t[i] = t[i-1] + pathLength_i/(v[i]+v[i-1])/2` if the origins of frames `i-1` and `i` do not coincide.
- `t[i] = t[i-1] + pathAngle_i/(v[i]+v[i-1])/2` if the origins of frames `i-1` and `i` do coincide.

Hereby `pathLength_i` is the distance between the origins of frames `i-1` and `i` in [m] and
`pathAngle_i` is the planar rotation angle between frames `i-1` and `i` in [rad].

If `v[i]` is the desired velocity or angular velocity at frame `i`, then path parameter
`t` is approximately the time to move along the path. The time instant `t_end` of the last frame can
be inquired with `ModiaMath.t_pathEnd(path)`. For example, if a simulation shall be performed in such
a way that the simulation should start with the first frame and end at `stopTime` at the last frame,
then the path parameter should be selected as `t = time*t_end/stopTime`.

Given the actual path parameter, typically `0 <= t <= t_end` (if `t` is outside of this interval,
the frame at `t` is determined by extrapolation through the first two or the last two frames),
the corresponding frame is determined by linear interpolation in the following way:

```julia
(rt, qt) = interpolate(  path,t)
 rt      = interpolate_r(path,t)
```

where `rt` is the position vector to the origin of the frame at path parameter `t`
and `qt` is the absolute quaternion of the frame at path parameter `t`.

# Example

```julia
import ModiaMath
using Unitful

r = [ ModiaMath.Vector3D(1,0,0),
      ModiaMath.Vector3D(0,1,0),
      ModiaMath.Vector3D(0,0,1) ]
q = [ ModiaMath.NullQuaternion,
      ModiaMath.qrot1(45u"°"),
      ModiaMath.qrot2(60u"°")]

path     = ModiaMath.Path(r,q)
t_end    = ModiaMath.t_pathEnd(path)
dt       = 0.1
stopTime = 2.0
time     = 0.0

while time <= stopTime
   (rt, qt) = ModiaMath.interpolate(path, time*t_end/stopTime)
   time += dt
end
```
"""
struct Path
    t::Vector{Float64}        # path parameter; t=0 is (r[1],q[1])
    r::Vector{Vector3D}       # Position vectors from world frame to origin of Frames
    q::Vector{Quaternion}     # Quaternions describing rotation from world frame to Frames

    function Path(r::AbstractVector, q::AbstractVector=Quaternion[];
                  v::AbstractVector=ones(length(r)),
                  seps=sqrt(eps()) )
        nframes = size(r, 1)
        @assert(seps > 0.0)
        @assert(nframes > 1)
        @assert(length(v) == nframes)
        @assert(size(q, 1) == 0 || size(q, 1) == nframes)
        @assert(v[1]   >= 0.0)
        @assert(v[end] >= 0.0)
        for i in 2:length(v) - 1
            @assert(v[i] > 0.0)
        end
        vv = convert(Vector{Float64}, v)

        # Determine path length for every segment
        t = zeros(nframes)
        for i in 2:nframes
            slen = norm(r[i] - r[i - 1])
            if slen > seps
                # Use distance between r[i] and r[i-1] as path parameter (in [m]), scale with vv
                t[i] = t[i - 1] + slen / ((vv[i] + vv[i - 1]) / 2)

            elseif length(q) > 0
                # Use planar rotation angle between q[i] and q[i-1] as path parameter (in [rad])
                q_rel = relativeRotation(q[i - 1], q[i])
                q4    =  q_rel[4]
                q4    = q4 >  1 + seps ?  1.0 :
                    q4 < -1 - seps ? -1.0 : q4
                absAngle = 2 * acos(q4)

                if absAngle < seps
                    error("ModiaMath.Path(..): r[i] == r[i-1] and q[i] == q[i-1] (i = ", i, ").")
                end
                t[i] = t[i - 1] + absAngle / ((vv[i] + vv[i - 1]) / 2)

            else
                error("ModiaMath.Path(..): r[i] == r[i-1] (i = ", i, ").")
            end
        end
        new(t, r, q)
    end
end

"""
    t_end = ModiaMath.t_pathEnd(path::[`ModiaMath.Path`](@ref))

Return the final path parameter `t`of the last frame in path
(path parameter of first frame = 0.0).
"""
t_pathEnd(path::Path)::Float64 = path.t[end]

function get_interval(path::Path, t::Number)
    # returns index i, such that path.t[i] <= t < path.t[i+1]

    tvec = path.t
    if t <= 0.0
        return 1
    elseif t >= tvec[end]
        return length(tvec) - 1
    end

    low  = 1
    high = length(tvec)
    mid  = round(Int, (low + high) / 2)

    while t < tvec[mid] || t >= tvec[mid + 1]
        if t < tvec[mid]
            high = mid
        else
            low = mid
        end
        mid = round(Int, (low + high) / 2, RoundDown)
    end

    return mid
end

"""
    (rt, qt) = ModiaMath.interpolate(path, t)

Return position `rt`and Quaternion `qt` of `path::`[`ModiaMath.Path`](@ref) at path parameter `t::Number`.
"""
function interpolate(path::Path, t::Number)
    i    = get_interval(path, t)
    tvec = path.t
    tt::Float64  = convert(Float64, t)
    fac::Float64 = (tt - tvec[i]) / (tvec[i + 1] - tvec[i])
    rt = path.r[i] + fac * (path.r[i + 1] - path.r[i])
    qt = length(path.q) > 0 ? normalize(path.q[i] + fac * (path.q[i + 1] - path.q[i])) : NulLQuaternion

    return (rt, qt)
end

"""
    rt = ModiaMath.interpolate_r(path, t)

Return position `r` of `path::`[`ModiaMath.Path`](@ref) at path parameter `t::Number`.
"""
function interpolate_r(path::Path, t::Number)::Vector3D
    i    = get_interval(path, t)
    tvec = path.t
    tt::Float64  = convert(Float64, t)
    fac::Float64 = (tt - tvec[i]) / (tvec[i + 1] - tvec[i])
    return path.r[i] + fac * (path.r[i + 1] - path.r[i])
end


text:
# update_G0.jl

function rpost_β0star_ηy(βstar_ηy::Array{T,2}, δ_y::Array{T,1}, Λ0star_ηy::PDMat{T},
                            β0star_ηy_mean::Array{T, 1}, β0star_ηy_Prec::PDMat{T}) where T <: Real

    ## Collect sufficient statistics
    nstar, Kstar = size(βstar_ηy)
    sumδinv = 0.0
    sumyδinv = zeros(T, Kstar)
    for i = 1:nstar
        sumδinv += 1.0 / δ_y[i]
        sumyδinv += βstar_ηy[i,:] ./ δ_y[i]
    end

    ## Updated parameters
    Prec_out = β0star_ηy_Prec + sumδinv * Λ0star_ηy # result is of type PDMat
    a = β0star_ηy_Prec * β0star_ηy_mean + Λ0star_ηy * sumyδinv
    mean_out = Prec_out \ a

    ## Sample
    return mean_out + Prec_out.chol.U \ randn(Kstar)
end

function rpost_Λ0star_ηy(βstar_ηy::Array{T,2}, δ_y::Array{T,1}, β0star_ηy::Array{T,1},
                                 df::T, invSc::PDMat{T}) where T <: Real
    ## Collect sufficient statistic
    nstar, Kstar = size(βstar_ηy)
    SS = zeros(T, Kstar, Kstar)
    for i = 1:nstar
        dev = (βstar_ηy[i,:] - β0star_ηy) ./ δ_y[i]
        SS +=  dev * dev'
    end

    ## Updated parameters
    df1 = df + nstar
    invSc1 = PDMat_adj(invSc + SS)
    Sc1 = inv(invSc1) # is there some way around this?

    ## Sample
    return PDMat_adj(rand(Distributions.Wishart(df1, Sc1)))
end

function rpost_IGs0_gammaPri(δ::Array{T,1}, ν::T, n0::T, s00::T) where T <: Real
    # δ is a vector of IG variates
    # ν is the degrees of freedom of the IG
    # the gamma prior has shape ν*n0/2 and rate ν*n0/(2*s00)
    cc = 0.5 * ν
    sh_out = cc * (n0 + length(δ))
    rate_out = cc * (n0 / s00 + sum( 1.0 ./ δ) )
    rand( Distributions.Gamma(sh_out, 1.0 / rate_out) )
end

function update_G0!(model::Model_BNP_WMReg_Joint)
    # ii = sort(unique(model.state.S)) # this may be the source of some numerical instabilities
    ii = 1:model.H

    #### We can't even consider updating G0 for η_y with these functions unless we have a current draw for η_y.
    # βstar_ηy = hcat(model.state.μ_y[ii], model.state.β_y[ii,:])

    # model.state.β0star_ηy = rpost_β0star_ηy(βstar_ηy,
    #         model.state.δ_y[ii], model.state.Λ0star_ηy,
    #         model.prior.β0star_ηy_mean, model.prior.β0star_ηy_Prec)

    # model.state.Λ0star_ηy = rpost_Λ0star_ηy(βstar_ηy, model.state.δ_y[ii],
    #                                  model.state.β0star_ηy,
    #                                  model.prior.Λ0star_ηy_df,
    #                                  model.prior.Λ0star_ηy_df * model.prior.Λ0star_ηy_S0)

    # model.state.s0_δy = rpost_IGs0_gammaPri(model.state.δ_y[ii], model.state.ν_δy,
    #                             model.prior.s0_δy_df, model.prior.s0_δy_s0)

    ###### the updates below appear to need all clusters because of the likelihood denominator (the above might too since we are conditioning on config variables, but it appears that they are okay)

    model.state.μ0_μx = BayesInference.rpost_MvN_knownPrec(model.state.μ_x[:,:],
                                model.state.Λ0_μx,
                                model.prior.μ0_μx_mean, model.prior.μ0_μx_Prec)

    model.state.Λ0_μx = PDMat_adj(BayesInference.rpost_MvNprec_knownMean(model.state.μ_x[:,:],
                                model.state.μ0_μx,
                                model.prior.Λ0_μx_df, model.prior.Λ0_μx_df * model.prior.Λ0_μx_S0))

    if model.Σx_type == :full && model.K > 1
        for k = 1:(model.K - 1)
            model.state.β0_βx[k] = BayesInference.rpost_MvN_knownPrec(
                model.state.β_x[k][:,:],
                model.state.Λ0_βx[k],
                model.prior.β0_βx_mean[k], model.prior.β0_βx_Prec[k])

            model.state.Λ0_βx[k] = PDMat_adj(BayesInference.rpost_MvNprec_knownMean(
                model.state.β_x[k][:,:],
                model.state.β0_βx[k],
                model.prior.Λ0_βx_df[k], model.prior.Λ0_βx_df[k] * model.prior.Λ0_βx_S0[k]))
        end
    end

    for k = 1:model.K
        model.state.s0_δx[k] = rpost_IGs0_gammaPri(model.state.δ_x[:,k], model.state.ν_δx[k],
            model.prior.s0_δx_df[k], model.prior.s0_δx_s0[k])
    end

    return nothing
end


text:
__precompile__(true)
module Hargreaves

using Cairo, Colors, LightGraphs, LinearAlgebra

function _get_weight(weight_matrix::Array{T,2} where T <: Number, e::AbstractEdge{T} where T <: Integer)
#This function only supports upper triangular matrices.
#(only unidirected graphs)

    weigth_matrix = UpperTriangular(weight_matrix)
    weight = 0
    if src(e) < dst(e) && src(e) != dst(e)
        weight = weight_matrix[src(e), dst(e)]
    elseif dst(e) < src(e) && src(e) != dst(e)
        weight = weight_matrix[dst(e), src(e)]
    end
    return weight
end

function _get_width(value::Number,
                    actual_min::Number,
                    actual_max::Number,
                    min_width::Number,
                    max_width::Number,
                    static_width_value::Number,
                    width_scale::Symbol,
                    static_widths::Bool)

    output_width = static_width_value
    normrange = max_width - min_width
    weight_range = actual_max - actual_min

    if value == 0
        output_width = 0
    else
        if !static_widths && !(normrange == 0) && !(weight_range == 0) && width_scale == :linear
            output_width = (value - actual_min) / (weight_range) * normrange
        elseif !static_widths && !(normrange == 0) && !(weight_range == 0) && width_scale == :quadratic
            # A proper quadratic function would be better.
            # The value max_width is exceeded by min_width * (actual_min / value) ^ 8.
            output_width = ((value - actual_min) / weight_range) ^ 2 * normrange + min_width * (actual_min / value) ^ 8
        end
        return output_width
    end
end

function _compute_node_positions(res_x, res_y, node_count, ratio=0.8, switch=1)
    center = [res_x / 2, res_y / 2]
    degrees = 2 * pi
    degree_sections = degrees / node_count
    main_radius = res_x / 2 * ratio
    positions = zeros(node_count, 2)
    for i in 1:node_count
        positions[i,1] = cos(degree_sections * (i - 1)) * main_radius + center[1]
        positions[i,2] = sin(degree_sections * (i - 1)) * main_radius + center[2]
    end
    return positions
end

function _get_color(value::Number,
                    actual_min::Number,
                    actual_max::Number,
                    static_colors::Bool,
                    color_scale,
                    scheme::Symbol = :generic)

    # This function only supports decreasing color values
    # while weights increase: IF-clause needed!

    (r_min, r_max) = (200, 112)
    (g_min, g_max) = (255, 30)
    (b_min, b_max) = (225, 112)
    if scheme == :sky
        (r_min, r_max) = (254, 8)
        (g_min, g_max) = (255, 29)
        (b_min, b_max) = (217, 88)
    elseif scheme == :blue
        (r_min, r_max) = (230, 0)
        (g_min, g_max) = (238, 55)
        (b_min, b_max) = (255, 160)
    elseif scheme == :purple
        (r_min, r_max) = (255, 75)
        (g_min, g_max) = (246, 0)
        (b_min, b_max) = (242, 106)
    end

    weight_range = actual_max - actual_min

    @inline function output_color(color_min, color_max)

        output = color_min / 255
        color_range = abs(color_max - color_min)

        if color_scale == :linear
            output = color_min
            if color_max < color_min
                output = (color_min - ((color_range) / weight_range) * value) / 255
            elseif color_max > color_min
                output = (color_max - ((color_range) / weight_range) * value) / 255
            end
        elseif color_scale == :quadratic
            output = color_min
            if color_max < color_min
                output = ((color_range / weight_range) * value) / 255
                output = color_min / 255 - output ^ 2
            elseif color_max > color_min
                output = (((color_range / weight_range) * value) ^ 2) / 255
            end
        end

        return output
    end # inline function

    r_out = output_color(r_min, r_max)
    g_out = output_color(g_min, g_max)
    b_out = output_color(b_min, b_max)

    if static_colors
    # Override min color used in legend's gradient
        r_min = r_max
        g_min = g_max
        b_min = b_max
    end
    return [
        [r_out,g_out,b_out],
        [0.0,0.0,0.0],
        [r_min / 255,g_min / 255,b_min / 255],
        [r_max / 255,g_max / 255,b_max / 255]
    ]
end

function wireplot(g::AbstractGraph, basefn = "wireplot";
    distmx = weights(g),
    res_x = 4096, res_y = 4096,
    plot_to_surface_ratio = 0.9,
    wireplot_bg_color = [1.0,1.0,1.0],
    wireplot_jitter = 0,
    wireplot_node_diameter = 8.0,
    wireplot_node_label_offset = 16.0,
    wireplot_node_label_font_color = [0.0,0.0,0.0],
    wireplot_node_label_font_size = 48.0,
    wireplot_edge_color_scheme =:generic,
    wireplot_legend_font_color = [0.0,0.0,0.0],
    wireplot_legend_heading = "Edges' weights",
    wireplot_legend_show_values = true,
    static_widths = false,
    static_colors = false,
    color_scale = :linear,
    width_scale = :linear,
    min_width = 1,  # minimum normalized arc width
    max_width = 8, # maximum normalized arc width
    static_width_value = 1, # static arc width (if static_widths)
    export_type = :svg
    )

    (actual_minwidth, actual_maxwidth) = extrema(distmx)
    actual_widthrange = actual_maxwidth - actual_minwidth
    @inline function _norm_width(w)
        if actual_widthrange == 0
            return 1
        else
            w_pct = (w - actual_minwidth) / actual_widthrange
            normrange = max_width - min_width
            if w_pct <= 0 || normrange <= 0
                return 1
            end
        end
        return output_min
    end

    actual_min = minimum(distmx)#_drop_zeros(distmx, actual_max)
    #println(actual_min)

    ## SETUP SURFACE
    #-----------------------

    center = [res_x / 2,res_y / 2]

    if export_type == :pdf
        c = CairoPDFSurface("$(basefn).pdf", res_x, res_y)
    else
        c = CairoSVGSurface("$(basefn).svg", res_x, res_y)
    end
    cr = CairoContext(c)
    select_font_face(cr, "Times", 1, 1)
    set_font_size(cr, wireplot_node_label_font_size)

    ## DRAW BACKGROUND
    #-----------------------

    set_source_rgb(cr, wireplot_bg_color...)
    rectangle(cr, 0.0, 0.0, res_x, res_y)
    fill(cr)

    ## PLOT EDGES
    #-----------------------

    plot_border_top = res_y
    plot_border_bottom = 0
    plot_border_left = res_x
    plot_border_right = 0

    node_pos = _compute_node_positions(res_x, res_y, nv(g), plot_to_surface_ratio, 1)

    degs = degree(g)

    for (i, e) in enumerate(edges(g))
        (s, d) = (src(e), dst(e))
        # indeg_d = indegs[d]
        # deg_d = indeg_d + outdeg_d

        if distmx == weights(g) || static_colors
            set_source_rgb(cr, _get_color(actual_maxwidth, actual_minwidth, actual_maxwidth, static_colors, color_scale, wireplot_edge_color_scheme)[1]...)
        else
            set_source_rgb(cr, _get_color(_get_weight(distmx, e), actual_minwidth, actual_maxwidth, static_colors, color_scale, wireplot_edge_color_scheme)[1]...)
        end
        line_width = _get_width(_get_weight(distmx, e),
                                      actual_minwidth,
                                      actual_maxwidth,
                                      min_width,
                                      max_width,
                                      static_width_value,
                                      width_scale,
                                      static_widths)
        set_line_width(cr, line_width)
        if line_width > 0
            curve_to(cr, node_pos[s, 1], node_pos[s, 2],
                center[1] + rand() * wireplot_jitter - rand() * wireplot_jitter,
                center[2] + rand() * wireplot_jitter - rand() * wireplot_jitter,
                node_pos[d, 1], node_pos[d, 2])
            stroke(cr)
        end
    end
    set_line_width(cr, 1)

    ## PLOT NODES BASED ON COMPUTED POSITIONS
    #----------------------------------------

    check = false

    for i in vertices(g)

        (pos_x, pos_y) = node_pos[i,:]
        set_source_rgb(cr, [0.0,0.0,0.0]...)
        move_to(cr, pos_x, pos_y) # Prevents artifacts in the exported pdf file
        circle(cr, pos_x, pos_y, wireplot_node_diameter)
        fill(cr)
        text = string(i)
        label_extents = text_extents(cr, text)

        set_source_rgb(cr, wireplot_node_label_font_color...)

        if round(Int, pos_x) == center[1] && pos_y > center[2]
            label_origin_x = pos_x - label_extents[3] / 2 - 1
            label_origin_y = pos_y + wireplot_node_diameter / 2 + label_extents[4] + wireplot_node_label_offset
            label_border_left = pos_x - label_extents[3] / 2 - 1
            label_border_right = label_origin_x + label_extents[3]
            label_border_top = label_origin_y + label_extents[4]
            label_border_bottom = label_origin_y
        elseif round(Int, pos_x) == center[1] && pos_y < center[2]
            label_origin_x = pos_x - label_extents[3] / 2 - 1
            label_origin_y = pos_y - wireplot_node_diameter / 2 - wireplot_node_label_offset
            label_border_left = pos_x - label_extents[3] / 2 - 1
            label_border_right = label_origin_x + label_extents[3]
            label_border_top = label_origin_y + label_extents[4]
            label_border_bottom = label_origin_y
        elseif pos_x < center[1]
            label_origin_x = pos_x - wireplot_node_diameter / 2 - label_extents[3] - wireplot_node_label_offset
            label_origin_y = pos_y + label_extents[4] / 2 + sin(((i - 1) * 2 * pi) / nv(g)) * (wireplot_node_label_offset + wireplot_node_diameter / 2 + label_extents[4] / 2)
            label_border_left = label_origin_x
            label_border_right = label_origin_x + label_extents[3]
            label_border_top = label_origin_y + label_extents[4]
            label_border_bottom = label_origin_y
        elseif pos_x > center[1]
            label_origin_x = pos_x + wireplot_node_diameter / 2 + wireplot_node_diameter / 2 + wireplot_node_label_offset
            label_origin_y = pos_y + label_extents[4] / 2 + sin(((i - 1) * 2 * pi) / nv(g)) * (wireplot_node_label_offset + wireplot_node_diameter / 2 + label_extents[4] / 2)
            label_border_left = label_origin_x
            label_border_right = label_origin_x + label_extents[3]
            label_border_top = label_origin_y + label_extents[4]
            label_border_bottom = label_origin_y
        end

        move_to(cr, label_origin_x, label_origin_y)
        show_text(cr, text)
        plot_border_bottom = max(plot_border_bottom, label_border_bottom)
        plot_border_right = max(plot_border_right, label_border_right)
        plot_border_top = min(plot_border_top, label_border_top)
        plot_border_left = min(plot_border_left, label_border_left)
    end
    ## PLOT LEGEND
    #-----------------------

    ## HEADING

    heading_extents = text_extents(cr, wireplot_legend_heading)

    set_source_rgb(cr, wireplot_legend_font_color...)
    move_to(cr, plot_border_right - res_x * 0.25, plot_border_bottom)
    show_text(cr, wireplot_legend_heading)
    line_to(cr, plot_border_right, plot_border_bottom)

    stroke(cr)

    ## COLOR DESCRIPTOR

    rectangle(cr, plot_border_right - res_x * 0.23 + heading_extents[3],  plot_border_bottom - res_y * 0.015,  res_x * 0.23 - heading_extents[3], res_y * 0.01)

    gradient = pattern_create_linear(plot_border_right - res_x * 0.23 + heading_extents[3], plot_border_bottom,  plot_border_right,  plot_border_bottom)

    pattern_add_color_stop_rgb(gradient, 0, _get_color(0, actual_minwidth, actual_maxwidth, static_colors, color_scale, wireplot_edge_color_scheme)[3]...)
    pattern_add_color_stop_rgb(gradient, 1, _get_color(0, actual_minwidth, actual_maxwidth, static_colors, color_scale, wireplot_edge_color_scheme)[4]...)
    set_source(cr, gradient)

    fill_preserve(cr)
    destroy(gradient)

    ## COLOR DESCRIPTOR FONT

    if wireplot_legend_show_values
        set_source_rgb(cr, wireplot_legend_font_color...)
        text_min = string("[", 1, ", ...")
        text_max = string(round(Int, actual_maxwidth), "]")

        move_to(cr, plot_border_right - res_x * 0.23 + heading_extents[3], plot_border_bottom - res_y * 0.02)
        show_text(cr, text_min)

        move_to(cr, plot_border_right - text_extents(cr, text_max)[3] - text_extents(cr, text_max)[1], plot_border_bottom - res_y * 0.02)
        show_text(cr, text_max)
    end

    ## EXPORT SVG AND PNG
    #--------------------

    if export_type == :png
        write_to_png(c, "$(basefn).png")
    else
        finish(c)
    end
end

export wireplot

end # module


text:
#    #function g(d, fm, to)
#    function g(d)
#      h = FramedPlot(title="Equity Curve")
#  #    x = lterows(d, fm)   #fm
#  #    x = gterows(d, to)   #to
#      x = [0:size(d, 1)]
#      y = equity(d["Close"])
#      y[1] = 1.0
#      pos = y .> 1
#      neg = y .<= 1
#      add(h, FillBetween(x[neg],y[neg],x,y,color=0x006bab5b))
#      add(h, FillBetween(x[pos],y[pos],x,y,color=0x00d66661))
#  #     add(h, Curve(y,x))
#      return h
#    end

# function plotequity(df::DataFrame, col::String, fm::Int,
#                                                 fd::Int,
#                                                 fy::Int,
#                                                 tm::Int,
#                                                 td::Int,
#                                                 ty::Int)
#   # shape the time series
#   d = gterows(df, fm,fd,fy)   #from
#   d = lterows(d, tm,td,ty)   #to
#   x = [0:size(d, 1)-1]
#   y = equity(d[col]) # take the column passed into function
#   y[1] = 1.0    # replace NA with starting equity
#
#   # useful plot variables
#   neg  = y .< 1.
#   y1   = ones(length(y))
#   onlybad = Float64[]
#   for i in 1:length(neg)
#     if neg[i] == false
#       push!(onlybad, 1.0)
#     else
#       j=i
#       push!(onlybad, y[j])
#     end
#   end
#
#   # build plot with canonical Julia colors
#   h = FramedPlot(title="Equity Curve")
#   add(h, FillBetween(x,y1, x,y, color=0x6bab5b))
#   add(h, FillBetween(x,onlybad, x,y1, color=0xd66661))
#   add(h, Curve(x,y, color=0x3b972e))
#   add(h, Curve(x,onlybad, color=0xc93d39))
#   add(h, Curve(x,y1, color=0xffffff))
#   return h
# end
#
# # default to using the entire date range
# plotequity(df::DataFrame) =
# plotequity(df::DataFrame, "Close", int(month(df[1,1])),
#                                    int(day(df[1,1])),
#                                    int(year(df[1,1])),
#                                    int(month(df[size(df, 1),1])),
#                                    int(day(df[size(df, 1),1])),
#                                    int(year(df[size(df, 1),1])))


text:
function pk_plot(pk::Vector{Float64}; title="PSIS diagnostic plot.",
    leg=:topleft, kwargs...)
    scatter(pk, xlab="Datapoint", ylab="Pareto shape k",
        marker=2.5, lab="Pk points", leg=leg)
  hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
  hline!([1], lab="pk = 1.0")
  title!(title)
end

export
    pk_plot


text:
using ASCIIPlots, Test

x = randn(20)
y = 21.0 * x + randn(20)
lineplot(x, y)


text:
struct Continuity3dModel{M} <: AbstractOceanModel
    ocean::M
    function Continuity3dModel(ocean::M) where {M}
        return new{M}(ocean)
    end
end
vars_state(cm::Continuity3dModel, ::Prognostic, FT) =
    vars_state(cm.ocean, Prognostic(), FT)

# Continuity3dModel is used to compute the horizontal divergence of u

vars_state(cm::Continuity3dModel, ::Auxiliary, T) = @vars()
vars_state(cm::Continuity3dModel, ::Gradient, T) = @vars()
vars_state(cm::Continuity3dModel, ::GradientFlux, T) = @vars()
vars_state(cm::Continuity3dModel, ::UpwardIntegrals, T) = @vars()
init_state_auxiliary!(cm::Continuity3dModel, _...) = nothing
init_state_prognostic!(cm::Continuity3dModel, _...) = nothing
@inline flux_second_order!(cm::Continuity3dModel, _...) = nothing
@inline source!(cm::Continuity3dModel, _...) = nothing
@inline update_penalty!(::RusanovNumericalFlux, ::Continuity3dModel, _...) =
    nothing

# This allows the balance law framework to compute the horizontal gradient of u
# (which will be stored back in the field θ)
@inline function flux_first_order!(
    m::Continuity3dModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    @inbounds begin
        u = state.u # Horizontal components of velocity
        v = @SVector [u[1], u[2], -0]

        # ∇ • (v)
        # Just using θ to store w = ∇h • u
        flux.θ += v
    end

    return nothing
end

# This is zero because when taking the horizontal gradient we're piggy-backing
# on θ and want to ensure we do not use it's jump
@inline wavespeed(cm::Continuity3dModel, n⁻, _...) = -zero(eltype(n⁻))

boundary_state!(
    ::CentralNumericalFluxSecondOrder,
    cm::Continuity3dModel,
    _...,
) = nothing

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

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""

@inline function boundary_state!(nf, cm::Continuity3dModel, args...)
    # hack for handling multiple boundaries for now
    # will fix with a future update
    boundary_conditions = (
        cm.ocean.problem.boundary_conditions[1],
        cm.ocean.problem.boundary_conditions[1],
        cm.ocean.problem.boundary_conditions[1],
    )
    return ocean_boundary_state!(nf, boundary_conditions, cm, args...)
end
#=
@inline function boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Continuity3dModel,
    Q⁺,
    A⁺,
    n⁻,
    Q⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    Q⁺.u = Q⁻.u

    return nothing
end
=#

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
function δ(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E<:AllowedSpeciesTypes}
    for s in species(N)
        @assert haskey(L, s)
    end

    this_l = filter(p -> p.first in species(N), L)
    tl = [this_l[s] for s in species(N; dims=1)]
    bl = [this_l[s] for s in species(N; dims=2)]
    D = tl .== bl'
    return D
end

"""
    Q(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E<:AllowedSpeciesTypes}

Modularity of a network and its partition. The second argument is a dictionary
where every species of `N` is associated to an `Int64` value representing the
identity of the module. This function returns the same value of bipartite
networks and their unipartite projection.
"""
function Q(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E<:AllowedSpeciesTypes}
    for s in species(N)
        @assert haskey(L, s)
    end

    # Degrees
    dkin, dkout = degree_in(N), degree_out(N)
    kin = map(x -> dkin[x], species(N; dims=2))
    kout = map(x -> dkout[x], species(N; dims=1))

    # Value of m -- sum of weights, total number of int, ...
    m = links(N)

    # Null model
    kikj = (kout .* kin')
    Pij = kikj ./ m

    # Difference
    diff = N.A .- Pij

    # Diff × delta
    dd = diff .* δ(N,L)

    return sum(dd)/m
end

"""
    Qr(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E<:AllowedSpeciesTypes}

Realized modularity -- this function returns a value giving the proportion of
all links that are within the same module. Higher values reflect a more strongly
modular partition (whereas `Q` represents the deviation of modularity from the
random expectation).
"""
function Qr(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E<:AllowedSpeciesTypes}
    for s in species(N)
        @assert haskey(L, s)
    end
    W = sum(N.A .* δ(N, L))
    B = links(N)
    return 2.0 * (W/B) - 1.0
end

function tidy_modules!(L::Dict{E,Int64}) where {E<:AllowedSpeciesTypes}
  l_values = sort(unique(collect(values(L))))
  for i in 1:length(l_values)
    for (k,v) in L
      if v == l_values[i]
        L[k] = i
      end
    end
  end
end


text:
#
#	Local Operator is a sparse operator structure only storing in and outgoing indices
#	and the coefficient associated with eacht in-out index pair:
#
#
#	σʸ = (1,2),(2,1) index pairs with coeffcients -1.0im and 1.0im
#
"""
    localOp

Local Operator is a sparse operator structure only storing in and outgoing indices
and the coefficient associated with each in-out index pair:

σʸ = (1,2),(2,1) index pairs with coeffcients -1.0im and 1.0im

To Do: !!! Essentially could be replaced by sparse array with coo indexing.
"""
struct localOp{DT<:Number, D}
	c::Vector{DT}
	in::Vector{Int}
 	out::Vector{Int}

 	function localOp(A::AbstractArray{<:Number,2})
		idx = findall(x->x!=0, A)
		DT = eltype(A)

		in = [x[1] for x in idx]
		out = [x[2] for x in idx]
		c = A[idx]

		return new{DT, size(A)[1]}(c, in, out)
	end

	function localOp(c::Vector{<:Number}, in::Vector{Int}, out::Vector{Int}, D::Int)
		# check dimensions
		D == maximum(in) || throw(ArgumentError("ingoing index exceeds the local operator dimension"))
		D == maximum(out) || throw(ArgumentError("outgoing index exceeds the local operator dimension"))
		length(c) <= D*D || throw(ArgumentError("specified entires of local operator exceeds dimension"))
		# check if each element is only referenced ones
		allunique(in) || throw(ArgumentError("ingoing indicies contains duplicate"))
		allunique(out) || throw(ArgumentError("outgoing indicies contains duplicate"))

		return new{eltype(c), D}(c, in, out)
	end
end

physdim(op::localOp{DT, D}) where{DT<:Number, D} = D

import Base.convert
function convert(::Type{Array{DT, 2}}, op::localOp{DT, D}) where {DT<:Number, D}
	OP = zeros(DT, D, D)
	@inbounds for x in eachoperation(op)
		OP[x[2],x[3]] = x[1]
	end

	return OP
end

import Base.Array
Array(op::localOp{DT, D}) where {DT<:Number, D} = convert(Array{DT,2}, op)

import Base.show
function Base.show(io::IO, op::localOp)
	print(io, "local operator: ", Array(op))
end

import Base.*
*(a::Number, b::localOp{DT,D}) where {DT<:Number, D}= localOp(a.*b.c, b.in, b.out, D)
*(b::localOp{DT,D}, a::Number) where {DT<:Number, D}= localOp(a.*b.c, b.in, b.out, D)

import Base.+
function +(a::localOp{DT1, D}, b::localOp{DT2, D}) where {DT1<:Number, DT2<:Number, D}

end

import Base.eltype
eltype(op::localOp{DT,D}) where {DT<:Number, D} = DT

# get list with tuples (coef, in_index, out_index)
function eachoperation(a::localOp)
	return zip(a.c, a.in, a.out)
end

# check if two localOP are parallel
function (∥)(a::localOp, b::localOp)
	parallel = false
	scaling = 0
	if a.in == b.in && a.out == b.out
		if iszero(diff(a.c./b.c))
			parallel = true
			scaling = (a.c./b.c)[1]
		end
	end

	return parallel, scaling
end

max_stars_repo_name: Gawatz/MPOmodule

max_issues_repo_name: Gawatz/MPOmodule

max_forks_repo_name: Gawatz/MPOmodule


text:
@inbounds @inline dist2(p, q)::Float64 = sqrt((p[1]-q[1])^2+(p[2]-q[2])^2)
# FMTree class
mutable struct FMTree
    s_init::Vec4f
    s_goal::Vec4f
    N #number of samples
    Pset::Vector{Vec4f} # Point set
    cost::Vector{Float64} #cost
    time::Vector{Float64} #optimal time to connect one node to its parent node
    parent::Vector{Int64} #parent node
    bool_unvisit::BitVector #logical value for Vunvisit
    bool_open::BitVector #logical value for Open
    bool_closed::BitVector #logical value for Closed
    world::World # simulation world config
    itr::Int64 # iteration num

    function FMTree(s_init::Vec4f, s_goal::Vec4f, N, world)
        # constructer: sampling valid point from the configurationspace
        println("initializing fmt ...")
        Pset = Vec4f[]
        push!(Pset, s_init) #inply idx_init = 1
        myrn(min, max) = min + (max-min)*rand()
        for n = 1:N-2
            while(true)
                p = Vec4f(myrn(world.x_min[1], world.x_max[1]),
                          myrn(world.x_min[2], world.x_max[2]),
                          myrn(world.v_min[1], world.v_max[1]),
                          myrn(world.v_min[1], world.v_max[1]))
                if isValid(world, p)
                    push!(Pset, p)
                    break
                end
            end
        end
        push!(Pset, s_goal) #inply idx_goal = N [last]
        cost = zeros(N)
        time = zeros(N)
        parent = ones(Int, N)
        bool_unvisit = trues(N)
        bool_unvisit[1] = false
        bool_closed = falses(N)
        bool_open = falses(N)
        bool_open[1] = true
        println("finish initializing")
        new(s_init, s_goal,
            N, Pset, cost, time, parent, bool_unvisit, bool_open, bool_closed, world, 0)
    end
end

function show(this::FMTree)
    println("drawing...")
    show(this.world)
    N = length(this.Pset)
    mat = zeros(2, N)
    for idx = 1:N
        mat[:, idx] = this.Pset[idx][1:2]
    end
    idxset_open = findall(this.bool_open)
    idxset_closed = findall(this.bool_closed)
    idxset_unvisit = findall(this.bool_unvisit)
    idxset_tree = setdiff(union(idxset_open, idxset_closed), [1])
    scatter(mat[1, 1], mat[2, 1], c=:blue, s=10, zorder = 100)
    scatter(mat[1, end], mat[2, end], c=:blue, s=10, zorder = 101)
    scatter(mat[1, idxset_open], mat[2, idxset_open], c=:gray, s=2)
    scatter(mat[1, idxset_closed], mat[2, idxset_closed], c=:gray, s=2)
    #scatter(mat[1, idxset_unvisit], mat[2, idxset_unvisit], c=:orange, s=5)
    for idx in idxset_tree
        s0 = this.Pset[this.parent[idx]]
        s1 = this.Pset[idx]
        tau = this.time[idx]
        show_trajectory(s0, s1, tau)
    end

    scatter(mat[1, 1], mat[2, 1], c=:blue, s=20, zorder = 100)
    scatter(mat[1, end], mat[2, end], c=:blue, s=20, zorder = 101)

    xlim(this.world.x_min[1]-0.05, this.world.x_max[1]+0.05)
    ylim(this.world.x_min[2]-0.05, this.world.x_max[2]+0.05)
    println("finish drawing")
end

function solve(this::FMTree, with_savefig = false)
    # keep extending the node until the tree reaches the goal
    println("please set with_savefig=false if you want to measure the computation time" )
    println("start solving")
    while(true)
        extend(this)
        if with_savefig
            if ((this.itr<100) & (this.itr % 20 == 1)) || (this.itr % 200==1)
                close()
                show(this)
                savefig("../fig/"*string(this.itr)*".png")
            end
        end
        !this.bool_unvisit[end] && break
    end
    idx = this.N
    idx_solution = Int64[idx]
    while(true)
        idx = this.parent[idx]
        push!(idx_solution, idx)
        idx == 1 && break
    end
    println("finish solving")
    return idx_solution
end

function extend(this::FMTree)
    # extend node
    this.itr += 1
    r = 1.0

    idxset_open = findall(this.bool_open)
    idxset_unvisit = findall(this.bool_unvisit)

    idx_lowest = idxset_open[findmin(this.cost[idxset_open])[2]]
    idxset_near, _, _ = filter_reachable(this.Pset, idxset_unvisit,
                                      this.Pset[idx_lowest], r, :F)

    for idx_near in idxset_near
        idxset_cand, distset_cand, timeset_cand = filter_reachable(this.Pset, idxset_open,
                                                     this.Pset[idx_near], r, :B)
        isempty(idxset_cand) && return
        cost_new, idx_costmin = findmin(this.cost[idxset_cand] + distset_cand)
        time_new = timeset_cand[idx_costmin] # optimal time for new connection
        idx_parent = idxset_cand[idx_costmin]
        waypoints = gen_trajectory(this.Pset[idx_parent], this.Pset[idx_near], time_new, 10)
        if isValid(this.world, waypoints)
            this.bool_unvisit[idx_near] = false
            this.bool_open[idx_near] = true
            this.cost[idx_near] = cost_new
            this.time[idx_near] = time_new
            this.parent[idx_near] = idx_parent
        end
    end
    this.bool_open[idx_lowest] = false
    this.bool_closed[idx_lowest] = true
end


text:
function compute_routing(rd::RoutingData, ::Load, ::MinimumMaximum, ::PathFormulation, ::Val{true}, ::DualReformulation, ::ObliviousUncertainty, ::UncertainDemand)
    # For column generation, more dual values are required for the pricing.
    old_tm = rd.model_robust_reformulation_traffic_matrices
    rd.model_robust_reformulation_traffic_matrices = true

    # Create the master model (same as without column generation).
    start = time_ns()
    rm = master_formulation(rd, PathFormulation())
    time_create_master_model_ms = (time_ns() - start) / 1_000_000.

    # Data structures to hold intermediate results.
    n_new_paths = 0
    routings = Routing[]
    objectives = Float64[]
    times_solve = Float64[]

    # Main loop.
    while true
        start = time_ns()

        # Solve the current master problem.
        m = rm.model
        optimize!(m)
        status = termination_status(m)
        push!(routings, Routing(rd, value.(rm.routing)))
        push!(objectives, objective_value(m))

        # Export if needed. TODO: use iteration number.
        _export_lp_if_allowed(rd, rm.model, "lp_master")
        _export_lp_if_failed(rd, status, rm.model, "error_master", "Subproblem could not be solved!")

        # Check if there are still columns to add.
        dual_value_convexity = dual.(rm.constraints_convexity) # Demand -> dual value
        dual_values_uncer = Dict{Edge{Int}, Dict{Edge{Int}, Float64}}(
            d => Dict(e => NaN for e in edges(rd))
            for d in demands(rd)
        ) # Demand -> edge -> dual value. All NaNs should be rewritten by the next loop.
        for (e, d_cstr) in rm.constraints_uncertainty_set
            for (d, constraint) in d_cstr
                dual_values_uncer[d][e] = dual(constraint)
            end
        end

        new_paths = _oblivious_reformulation_solve_pricing_problem(rd, dual_values_uncer, dual_value_convexity)
        n_new_paths += length(new_paths)

        # Add the new columns to the formulation.
        for (d, path) in new_paths
            # Add the new columns to the data object and retrieve a path ID.
            path_id = add_path!(rd, d, path)

            # Get your hands dirty in the formulation.
            # - Create a new variable for this path, add it into the
            #   convexity constraint.
            add_routing_var!(rm, demand, path_id)

            for e in edges(rd)
                # - Add it into the uncertainty sets.
                set_normalized_coefficient(rm.constraints_uncertainty_set[e][d], rm.routing[d, path_id], 1)

                # - Create the new constraint.
                @constraint(m, rm.dual_alpha[e, d] + sum(rm.dual_beta[e, e2] for e2 in path) >= 0)
            end
        end

        push!(times_solve, (time_ns() - start) / 1_000_000.)

        if length(new_paths) == 0
            break
        end
    end

    # Restore the original value.
    rd.model_robust_reformulation_traffic_matrices = old_tm

    # Export if needed.
    status = termination_status(rm.model)
    _export_lp_if_allowed(rd, rm.model, "lp_master")
    _export_lp_if_failed(rd, status, rm.model, "error_master", "Subproblem could not be solved!")

    # Retrieve the traffic matrices, if asked.
    matrices = if rd.model_robust_reformulation_traffic_matrices
        get_traffic_matrices(rm)
    else
        Dict{Edge{Int}, Float64}[]
    end

    stop = time_ns()

    # TODO: Export things like the normal case.

    return RoutingSolution(rd,
                           time_precompute_ms=rd.time_precompute_ms,
                           time_create_master_model_ms=time_create_master_model_ms,
                           time_solve_ms=times_solve,
                           objectives=objectives,
                           matrices=matrices,
                           routings=routings,
                           master_model=rm)
end

function _oblivious_reformulation_solve_pricing_problem(rd::RoutingData, dual_values_uncer, dual_value_convexity)
    paths = Dict{Edge{Int}, Vector{Edge}}()

    for demand in demands(rd)
        # Determine the weight matrix to use for the computations:
        #     for d: \sum_{e\in p} dual_{e, d}
        rd.logmessage("Pricing for $demand")
        weight_matrix = spzeros(Float64, n_nodes(rd), n_nodes(rd))
        for e in edges(rd)
            if dual_values_uncer[demand][e] <= - CPLEX_REDUCED_COST_TOLERANCE
                weight_matrix[src(e), dst(e)] = dual_values_uncer[demand][e]
            end
        end
        weight_matrix .+= maximum(abs.(weight_matrix))

        # Check that at least some weights are nonzero.
        if nnz(weight_matrix) == 0
            # No need for the potentially expensive _check_weight_matrix,
            # as no value in weight_matrix can be below the tolerance by
            # construction (no addition, unlike iterative version of oblivious).
            continue
        end

        # Compute the shortest path for this demand.
        state = desopo_pape_shortest_paths(rd.g, src(demand), weight_matrix)

        # Is there a solution? Has this path a negative reduced cost?
        if all(state.parents .== 0)
            continue
        end
        reduced_cost = state.dists[dst(demand)] - dual_value_convexity[demand]
        if reduced_cost >= - CPLEX_REDUCED_COST_TOLERANCE
            continue
        end
        rd.logmessage(reduced_cost)

        # Build the path.
        path = _build_path(state, demand)
        if path in rd.paths_edges
            rd.logmessage("New path $(path) is already in the formulation!")
        else
            paths[demand] = path
        end
    end

    return paths
end


text:
import InteractiveUtils, Functors
using Functors: @functor
function implementations(T::Type)
    return mapreduce(t -> isabstracttype(t) ? implementations(t) : t, vcat, InteractiveUtils.subtypes(T), init=[])
end

issame(v1::AbstractVertex, v2::AbstractVertex, visited=AbstractVertex[]) = false
function issame(v1::T, v2::T, visited=AbstractVertex[]) where T<:AbstractVertex
    v1 in visited && return true
    v2 in visited && return true
    push!(visited, v1,v2)
    for n in fieldnames(typeof(v1))
        issame(getfield(v1, n), getfield(v2, n), visited) || return false
    end
    return true
end
function issame(a1::AbstractArray{T,1}, a2::AbstractArray{T,1}, visited=AbstractVertex[]) where T<:AbstractVertex
    length(a1) != length(a2) && return false
    return all(map(vs -> issame(vs..., visited), zip(a1,a2)))
end
issame(g1::CompGraph, g2::CompGraph, visited=AbstractVertex[]) = issame(g1.outputs, g2.outputs, visited)
issame(d1, d2, visited=AbstractVertex[]) = d1 == d2

function showstr(f, v)
    buffer = IOBuffer()
    f(buffer, v)
    return String(take!(buffer))
end

# Testing mocks
mutable struct IndMem{F}
    wrapped::F
    lastins::Vector{Union{Missing, Vector{Int}}}
    lastouts::Vector{Int}
end
IndMem(w) = IndMem(w, convert(Vector{Union{Missing, Vector{Int}}}, map(i -> collect(1:i), nin(w))), collect(1:nout(w)))
IndMem(w, is::AbstractVector{<:Integer}, os::Integer) = IndMem(w, convert(Vector{Union{Missing, Vector{Int}}}, map(i -> collect(1:i), is)), collect(1:os))
IndMem(w, is::Tuple{Vararg{Integer}}, os::Integer) = IndMem(w, collect(is), os)
IndMem(w, is::Tuple{Vararg{AbstractVertex}}, os) = IndMem(w, nout.(is), os)
IndMem(w, is::Integer, os) = IndMem(w, [is], os)
IndMem(w, is::AbstractVertex, os) = IndMem(w, nout(is), os)

@functor IndMem

(im::IndMem)(x...) = im.wrapped(x...)
NaiveNASlib.nout(im::IndMem) = nout(im.wrapped)
NaiveNASlib.nin(im::IndMem) = nin(im.wrapped)

lastins(v::AbstractVertex) = lastins(computation(v))
lastins(f) = missing
lastins(im::IndMem) = im.lastins

lastouts(v::AbstractVertex) = lastouts(computation(v))
lastouts(f) = missing
lastouts(im::IndMem) = im.lastouts

function NaiveNASlib.Δsize!(im::IndMem, ins::AbstractVector, outs::AbstractVector)
    im.lastins = ins
    im.lastouts = outs
    Δsize!(im.wrapped, ins, outs)
end

mutable struct SizeDummy
    nin::Vector{Int}
    nout::Int
end
SizeDummy(insize::Integer, outsize::Integer) = SizeDummy([insize], outsize)

@functor SizeDummy

function NaiveNASlib.Δsize!(sd::SizeDummy, ins::AbstractVector{<:Integer}, out::Integer)
    sd.nin .= ins
    sd.nout = out
    nothing
end
NaiveNASlib.nout(sd::SizeDummy) = sd.nout
NaiveNASlib.nin(sd::SizeDummy) = sd.nin

NaiveNASlib.Δsizetype(::SizeDummy) = NaiveNASlib.ScalarSize()

issame(sd1::SizeDummy, sd2::SizeDummy, visited=AbstractVertex[]) = nout(sd1) == nout(sd2) && nin(sd1) == nin(sd2)

mutable struct MatMul{M<:AbstractMatrix}
    W::M
end
MatMul(nin, nout) = MatMul(reshape(collect(1:nin*nout), nout,nin))

@functor MatMul

(mm::MatMul)(x) = mm.W * x

function NaiveNASlib.Δsize!(mm::MatMul, ins::AbstractVector, outs::AbstractVector)
    mm.W = NaiveNASlib.parselect(mm.W, 1=>outs, 2=>ins[1])
    nothing
end
NaiveNASlib.nout(mm::MatMul) = size(mm.W, 1)
NaiveNASlib.nin(mm::MatMul) = [size(mm.W, 2)]

lastins(im::IndMem{<:MatMul}) = im.lastins[1]

computation(v::AbstractVertex) = computation(base(v))
computation(v::CompVertex) = v.computation
computation(::InputVertex) = identity


text:
module BufferedArrays

import Base: size, IndexStyle, getindex, setindex

using ..TiledArrays

export BufferedArray

struct BufferedArray{T, N, A} <: AbstractTiledArray{T, N}
    parent::A
    buffer::Array{T, N}
    tile::Ref{Tile{NTuple{N, Int}}}
    function BufferedArray{T, N, A}(parent, buffer, tile) where {T, N, A}
        new{T, N, A}(parent, buffer, Ref(tile))
end

BufferedArray{T, N, A}(parent, dims::NTuple{N, Integer}) where {T, N, A} =
    BufferedArray{T, N, A}(parent, Array{T, N}(undef, dims))

BufferedArray(parent::AbstractArray{T, N}, args...) where {T, N} =
    BufferedArray{T, N, typeof(A)}(parent, args...)

size(B::BufferedArray) = size(B.parent)

# Indexing

tile(B::BufferedArray) = B.tile[]

function getindex(B::BufferedArray, t::Tile)
    if t != tile(B)
        writeblock!(B, t)
    end
    B.buffer
end

end # module


text:
# Directory listing function test
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.4.0
# Created: 01.11.2013
# Last Edit: 01.04.2020

"""# dirlist(directory::ASCIIString;recur_depth::Int,regular_expression::Regex)

`(files,folders) = dirlist(\"K:\\\\Code\\\\\")` List all the files in the Code directory and all its subdirectories

`(files,folders) = dirlist([\"K:\\\\Code\\\\Matlab\\\\\",\"K:\\\\Code\\\\Julia\\\\\"])` List all the files in the K:\\Code\\Matlab\\ and K:\\Code\\Julia\\ directories and all of their subdirectories

`(files,folders) = dirlist(\"K:\\\\Code\\\\\",recur=2)` List all the files in the Code directory and one levels of directories below that.

`(files,folders) = dirlist([\"K:\\\\Code\\\\\"],recur=2)` List all the files and folders in the Code directory and one levels below that.

`(files,folders) = dirlist([\"K:\\\\Code\\\\\"],recur=2,regex=r\"\\.jl\$\")` Limit the files to those that match the regular expression `r\"\\.jl\$\"`

---

#### Keywords:\n
* recur::Int = Subdirectory recursion depth. 1 is the root directory.
* regex::Regex = Regular expression\n\n

---

\n
#### Regular Expression Examples\n
* Regular expression reference: http://www.regular-expressions.info
* `r\"(AirTemp).*\\d{8}-\\d{6}\"` Search for names with AirTemp followed by any number of characters then followed by 8 numbers, a dash, and 6 more numbers
* `r\"\\d{8}-\\d{6}\"` Search for names with eight digits followed by a dash and another six digits (12345678-654321)
* `r\"(example)\"` Search for all files with \"example\" in the name.
* `r\"\\.jl\$\"` Search for .jl at the end of the file. The . must be escaped with a \\
* `r\"(slt|cfg)\$\"` Search for files that end in slt or cfg\n\n

---

\n
#### Related Functions\n
* joinpath
* isdir
* isfile
* splitdir
* splitdrive
* basename
* dirname"""
function dirlist(directory::String;args...) #function dirlist{T <: String}(directory::Vector{T};args...)
    return dirlist([directory];args...)
end

#function dirlist{T <: String}(directories::Array{T,1};recur=typemax(Int),regex=r"")
function dirlist(directories::Array{T,1};recur::Int=typemax(Int),regex::Regex=r"") where T <: String
	###################################
	##  Regular Expression Examples  ##
	###################################
	# Regular expression reference: http://www.regular-expressions.info
	# reg = r"(AirTemp).*\d{8}-\d{6}" # Search for names with AirTemp followed by any number of characters then followed by 8 numbers, a dash, and 6 more numbers
	# r"\d{8}-\d{6}" # Search for names with eight digits followed by a dash and another six digits (12345678-654321)
	# r"(example)" # Search for all files with "example" in the name
	# r"\.jl$" # Search for .jl at the end of the file. The . must be escaped with a \

    # Parse Inputs
    rcr_max = Int[]
    try
		rcr_max = Int(recur)
	catch e
		error("recur must be convertable to Int")
	end

	# Initialize variables to be used inside the loop
	flist = Vector{String}() # File list
	dlist = Vector{String}() # Directory list
	rcr = 1
	lst = Vector{String}()
	nextdirs = Vector{String}()

	# Iterate through the given directory and its subdirectories
	while rcr <= rcr_max && ~isempty(directories)
		# While the maximum recursion level hasn't been reached and "directories" is not empty
		nextdirs = Vector{String}() # List of directories for the next iteration
		for i=1:1:length(directories)
			lst = Vector{String}() # Reset the file list each cycle
			try # May fail on certain folders like K:\__Papierkorb__\
				lst = joinpath.(directories[i],readdir(directories[i])) # List contents of specified directory and reformat
			catch
				println("FAILED READING " * directories[i])
				println("Continuing...")
			end
			fDirs = isdir.(lst)

			append!(flist,lst[fDirs .== false]) # Add the new files
			append!(dlist,lst[fDirs]) # Add the new directories
			append!(nextdirs,lst[fDirs]) # Add the directories to a list of future directories to search
		end
		directories = nextdirs
		rcr += 1 # Next level of recursion
	end

    if regex != r""
		# Find specific files if regular expression exists
		filter!(x -> occursin(regex,x),flist)
	end

    return flist,dlist # Return the results of the file and directory listing
end


text:
module PlutoProfile

using AbstractTrees
using FlameGraphs
using Profile
using ProfileSVG

export @plutoprofview

function pluto_subtree(g)
    for n in PreOrderDFS(g)
        if (
            endswith(string(n.data.sf.file), "PlutoRunner.jl") &&
            endswith(string(n.data.sf.func), "run_inside_trycatch")
        )
            return n
        end
    end
    return g
end

"""
    @plutoprofview f(args...)
Clear the Profile buffer, profile `f(args...)`, and view the result graphically excluding artefacts due to Pluto's use of multiprocessing.
"""
macro plutoprofview(ex)
    return quote
        Profile.clear()
        @profile $(esc(ex))
        ProfileSVG.view(pluto_subtree(flamegraph()))
    end
end

end


text:
@require NLopt="76087f3c-5699-56af-9a33-bf431cd00edd" begin
    function __map_optimizer_args(prob::OptimizationProblem, opt::MultistartOptimization.TikTak;
        cb=nothing,
        maxiters::Union{Number, Nothing}=nothing,
        maxtime::Union{Number, Nothing}=nothing,
        abstol::Union{Number, Nothing}=nothing,
        reltol::Union{Number, Nothing}=nothing,
        local_method::Union{NLopt.Algorithm, Nothing} = nothing,
        local_maxiters::Union{Number, Nothing} = nothing,
        local_maxtime::Union{Number, Nothing} = nothing,
        local_options::Union{NamedTuple,Nothing} = nothing,
        kwargs...)

        if !isnothing(abstol)
            @warn "common abstol is currently not used by the global method $(opt). Set the absolute tolerance of the optimized states for the local method via local_options using the keyword xtol_abs"
        end

        if !isnothing(reltol)
            @warn "common reltol is currently not used by the global method $(opt). Set the relative tolerance of the optimized states for the local method via local_options using the keyword xtol_rel"
        end

        if !isnothing(maxtime)
            @warn "common maxtime is currently not used by the global method $(opt). Set maxtime of local optimiser via local_maxtime."
        end

        if !isnothing(maxiters)
            @warn "common maxiters is currently not used by the global method $(opt). Set maxiters of local optimiser via local_maxiters."
        end

        if !isnothing(cb)
            @warn "callbacks are currently not used by the global method $(opt)"
        end

        global_meth = opt

        if isa(local_method, NLopt.Algorithm)
            local_kwargs = (; )
            if !(isnothing(local_maxiters))
                local_kwargs = (; local_kwargs..., maxeval=local_maxiters)
            end

            if !(isnothing(local_maxtime))
                local_kwargs = (; local_kwargs..., maxtime=local_maxtime)
            end

            if !(isnothing(local_options))
                local_kwargs = (; local_kwargs..., local_options...)
            end

            local_meth = MultistartOptimization.NLoptLocalMethod(local_method; local_kwargs...)
        else
            error("A local method has to be defined using algorithm of type NLopt.Algorithm e.g. 'NLopt.LN_NELDERMEAD()'.")
        end

        return global_meth, local_meth
    end

    function __solve(prob::OptimizationProblem, opt::MultistartOptimization.TikTak;
                    maxiters::Union{Number, Nothing}=nothing,
                    maxtime::Union{Number, Nothing}=nothing,
                    abstol::Union{Number, Nothing}=nothing,
                    reltol::Union{Number, Nothing}=nothing,
                    local_method::Union{NLopt.Algorithm, Nothing} = nothing,
                    local_maxiters::Union{Number, Nothing} = nothing,
                    local_maxtime::Union{Number, Nothing} = nothing,
                    local_options::Union{NamedTuple,Nothing} = nothing,
                    progress = false,
                    use_threads=true,
                    kwargs...)
        local x, _loss

        @warn "This method might soon be deprecated; if you are on MultistartOptimization v0.1.2#master or higher consider using the new method solve(prob::OptimizationProblem, multiopt::MultistartOptimization.TikTak, opt; use_threads=true,kwargs...) instead."

        maxiters = _check_and_convert_maxiters(maxiters)
        maxtime = _check_and_convert_maxtime(maxtime)

        local_maxiters = _check_and_convert_maxiters(local_maxiters)
        local_maxtime = _check_and_convert_maxtime(local_maxtime)

        _loss = function(θ)
            x = prob.f(θ, prob.p)
            return first(x)
        end

        opt_setup = MultistartOptimization.MinimizationProblem(_loss, prob.lb, prob.ub)
        multistart_method, local_method = _map_optimizer_args(prob, opt, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol,  local_method=local_method, local_maxiters =local_maxiters, local_maxtime=local_maxtime, local_options=local_options; kwargs...)

        t0 = time()
        opt_res = MultistartOptimization.multistart_minimization(multistart_method, local_method, opt_setup; use_threads=use_threads)
        t1 = time()
        opt_ret = hasproperty(opt_res, :ret) ? opt_res.ret : nothing

        SciMLBase.build_solution(prob, opt, opt_res.location, opt_res.value; (isnothing(opt_ret) ? (; original=opt_res) : (; original=opt_res, retcode=opt_ret))... )
    end
end

function __solve(prob::OptimizationProblem, multiopt::MultistartOptimization.TikTak, opt;
                use_threads=true,
                kwargs...)
    local x, _loss

    _loss = function(θ)
        x = prob.f(θ, prob.p)
        return first(x)
    end

    opt_setup = MultistartOptimization.MinimizationProblem(_loss, prob.lb, prob.ub)

    _local_optimiser = function(pb, θ0, prob)
        prob_tmp = remake(prob, u0=θ0)
        res = solve(prob_tmp, opt;
            kwargs...)
        return (value = res.minimum, location = res.minimizer, ret = res.retcode)
    end

    local_optimiser(pb, θ0) = _local_optimiser(pb, θ0, prob)

    t0 = time()
    opt_res = MultistartOptimization.multistart_minimization(multiopt, local_optimiser, opt_setup; use_threads=use_threads)
    t1 = time()
    opt_ret = hasproperty(opt_res, :ret) ? opt_res.ret : nothing

    SciMLBase.build_solution(prob, (multiopt, opt), opt_res.location, opt_res.value; (isnothing(opt_ret) ? (; original=opt_res) : (; original=opt_res, retcode=opt_ret))... )
end


text:
function day_test_functional()
    test_courier_bike_monday()
    test_prior_failure()
end

function test_courier_bike_monday()
    ok, schedule = schedule_day(bike_monday_employees, bike_monday_env)
    @test ok == true
    @test day_sum_coverage(bike_monday_env) * TEST_THRESHOLD > day_hours_scheduled(schedule)
end

function test_prior_failure()
    ok, schedule = schedule_day(test_prior_failure_employees, test_prior_failure_env)
    @test ok == true
    @test day_sum_coverage(test_prior_failure_env) * TEST_THRESHOLD > day_hours_scheduled(schedule)
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
#=
1. Create a right cylinder C:
	(a) the base of C is a regular polygon with N sides.
	(b) the base is parallel to the xy-plane.
2. Create a copy C′ of C.
3. Rotate C′ around its vertical centerline by α degrees.
4. Measure time for C ∪ C′.
=#

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation;


text:
"""
This file istarts the ClimateERA module by creating the root directory and by
specifying whether the data is to be downloaded or analyzed.  Functionalities
include:
    - Creation of root directory
    - Specifying whether purpose is to download or analyse data

"""

## ClimateERA Folders

function eraregfolder(ereg::Dict,eroot::Dict)

    fol = joinpath(eroot["era"],ereg["fol"]);

    if !isdir(fol)
        @info "$(Dates.now()) - The folder for the $(ereg["name"]) region at does not exist.  Creating now ..."
        mkpath(fol);
    end

    return fol

end

function eravarfolder(epar::Dict,ereg::Dict,eroot::Dict)

    fol = joinpath(eroot["era"],ereg["fol"],epar["ID"]);

    if !isdir(fol)
        @info "$(Dates.now()) - The folder for the $(epar["name"]) parameter in the $(ereg["name"]) region does not exist.  Creating now ..."
        mkpath(fol);
    end

    return fol

end

function erarawfolder(epar::Dict,ereg::Dict,eroot::Dict)

    if !haskey(epar,"level")

        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"raw");
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for raw data of the $(epar["name"]) parameter in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    else

        if epar["level"] == "sfc";

            fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"raw");
            if !isdir(fol)
                @info "$(Dates.now()) - The folder for raw data of the $(epar["name"]) parameter in the $(ereg["name"]) region does not exist.  Creating now ..."
                mkpath(fol);
            end

        else

            phPa = "$(epar["ID"])-$(epar["level"])hPa"
            fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],phPa,"raw");
            if !isdir(fol)
                @info "$(Dates.now()) - The folder for raw data of the $(epar["name"]) parameter at pressure level $(epar["level"])hPa in the $(ereg["name"]) region does not exist.  Creating now ..."
                mkpath(fol);
            end
        end

    end

    return fol

end

function erarawfolder(epar::Dict,ereg::Dict,eroot::Dict,date::TimeType)

    yr = "$(yr2str(date))";

    if epar["level"] == "sfc";

        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"raw",yr);
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for raw data of the $(epar["name"]) parameter in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    else

        phPa = "$(epar["ID"])-$(epar["level"])hPa"
        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],phPa,"raw",yr);
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for raw data of the $(epar["name"]) parameter at pressure level $(epar["level"])hPa in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    end

    return fol

end

function eraanafolder(epar::Dict,ereg::Dict,eroot::Dict)

    if epar["level"] == "sfc";

        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"ana");
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for analyzed data of the $(epar["name"]) parameter in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    else

        phPa = "$(epar["ID"])-$(epar["level"])hPa"
        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],phPa,"ana");
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for analyzed data of the $(epar["name"]) parameter at pressure level $(epar["level"])hPa in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    end

    return fol

end

function eraimgfolder(epar::Dict,ereg::Dict,eroot::Dict)

    if epar["level"] == "sfc";

        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],"img");
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for images/maps of the $(epar["name"]) parameter in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    else

        phPa = "$(epar["ID"])-$(epar["level"])hPa"
        fol = joinpath(eroot["era"],ereg["fol"],epar["ID"],phPa,"img");
        if !isdir(fol)
            @info "$(Dates.now()) - The folder for images/maps of the $(epar["name"]) parameter at pressure level $(epar["level"])hPa in the $(ereg["name"]) region does not exist.  Creating now ..."
            mkpath(fol);
        end

    end

    return fol

end

## ClimateERA Folder Setup

function erafolder(emod::Dict,epar::Dict,ereg::Dict,etime::Dict,eroot::Dict)

    yrbeg = etime["Begin"]; yrend = etime["End"];

    folreg = joinpath(eroot["era"],ereg["fol"]);
    if !isdir(folreg)
        @debug "$(Dates.now()) - Creating folder for the $(ereg["name"]) region at $(folreg) ..."
        mkpath(folreg);
    else; @debug "$(Dates.now()) - The folder for the $(ereg["name"]) region $(folreg) exists."
    end

    folvar = joinpath(eroot["era"],ereg["fol"],epar["ID"]);
    if !isdir(folvar)
        @debug "$(Dates.now()) - Creating variable folder for the $(epar["name"]) parameter at $(folvar) ..."
        mkpath(folvar);
    else; @debug "$(Dates.now()) - The folder for the $(epar["name"]) parameter $(folvar) exists."
    end

    if epar["level"] != "sfc"
        folvar = joinpath(folvar,"$(epar["ID"])-$(epar["level"])hPa")
    end
    folraw = joinpath(folvar,"raw"); foltmp = joinpath(folvar,"tmp");
    folana = joinpath(folvar,"ana"); folimg = joinpath(folvar,"img");

    @debug "$(Dates.now()) - Creating relevant subdirectories for data downloading, temporary storage, analysis and image creation."
    if !isdir(folraw);
        @debug "$(Dates.now()) - Creating folder for downloaded raw data: $(folraw)";
        for yr = yrbeg : yrend; mkpath(joinpath(folraw,"$(yr)")); end
    end
    if !isdir(folana);
        @debug "$(Dates.now()) - Creating folder for data analysis ouput: $(folana)"; mkpath(folana)
    end
    if !isdir(folimg);
        @debug "$(Dates.now()) - Creating folder for data images: $(folimg)"; mkpath(folimg)
    end
    if emod["actionID"] == 1 && !isdir(foltmp);
        @debug "$(Dates.now()) - Creating folder for temporary data storage: $(foltmp)";
        mkpath(foltmp)
    end

    return Dict(
        "reg"=>folreg,"var"=>folvar,"raw"=>folraw,
        "tmp"=>foltmp,"ana"=>folana,"img"=>folimg
    );

end

function erafolder(emod::Dict,epar::Dict,ereg::Dict,etime::Dict,eroot::Dict,pre)

    yrbeg = etime["Begin"]; yrend = etime["End"];

    folreg = joinpath(eroot["era"],ereg["fol"]);
    if !isdir(folreg)
        @debug "$(Dates.now()) - Creating folder for the $(ereg["name"]) region at $(folreg) ..."
        mkpath(folreg);
    else; @debug "$(Dates.now()) - The folder for the $(ereg["name"]) region $(folreg) exists."
    end

    folvar = joinpath(eroot["era"],ereg["fol"],epar["ID"]);
    if !isdir(folvar)
        @debug "$(Dates.now()) - Creating variable folder for the $(epar["name"]) parameter at $(folvar) ..."
        mkpath(folvar);
    else; @debug "$(Dates.now()) - The folder for the $(epar["name"]) parameter $(folvar) exists."
    end

    if pre != "sfc"; folvar = joinpath(folvar,"$(epar["ID"])-$(pre)hPa"); end
    folraw = joinpath(folvar,"raw"); foltmp = joinpath(folvar,"tmp");
    folana = joinpath(folvar,"ana"); folimg = joinpath(folvar,"img");

    @info "$(Dates.now()) - Creating relevant subdirectories for data downloading, temporary storage, analysis and image creation."
    if !isdir(folraw);
        @debug "$(Dates.now()) - Creating folder for downloaded raw data: $(folraw)";
        for yr = yrbeg : yrend; mkpath(joinpath(folraw,"$(yr)")); end
    end
    if !isdir(folana);
        @debug "$(Dates.now()) - Creating folder for data analysis ouput: $(folana)"; mkpath(folana)
    end
    if !isdir(folimg);
        @debug "$(Dates.now()) - Creating folder for data images: $(folimg)"; mkpath(folimg)
    end
    if emod["actionID"] == 1 && !isdir(foltmp);
        @debug "$(Dates.now()) - Creating folder for temporary data storage: $(foltmp)";
        mkpath(foltmp)
    end

    return Dict(
        "reg"=>folreg,"var"=>folvar,"raw"=>folraw,
        "tmp"=>foltmp,"ana"=>folana,"img"=>folimg
    );

end

## ClimateERA NetCDF Names

function erarawname(emod::Dict,epar::Dict,ereg::Dict,date::TimeType)

    if !(emod["levels"][1] == "sfc")
          fname = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-$(epar["level"])hPa";
    else; fname = "$(emod["prefix"])-$(ereg["fol"])-$(epar["ID"])-sfc";
    end

    return "$(fname)-$(yrmo2str(date)).nc"

end

function eraananame(emod::Dict,epar::Dict,ereg::Dict,date::TimeType)

    prefix = replace(emod["prefix"],"era"=>"eraa")

    if !(emod["levels"][1] == "sfc")
          fname = "$(prefix)-$(ereg["fol"])-$(epar["ID"])-$(epar["level"])hPa";
    else; fname = "$(prefix)-$(ereg["fol"])-$(epar["ID"])-sfc";
    end

    return "$(fname)-$(yr2str(date)).nc"

end

function eracmpname(emod::Dict,epar::Dict,ereg::Dict)

    prefix = replace(emod["prefix"],"era"=>"erac")

    if !(emod["levels"][1] == "sfc")
          return "$(prefix)-$(ereg["fol"])-$(epar["ID"])-$(epar["level"])hPa.nc";
    else; return "$(prefix)-$(ereg["fol"])-$(epar["ID"])-sfc.nc";
    end

end

erancread(ncname::AbstractString,fol::AbstractString="") = Dataset(joinpath(fol,ncname))

function erarawread(
    emod::Dict,epar::Dict,ereg::Dict,eroot::Dict,
    date::TimeType
)

    ebase = erarawfolder(epar,ereg,eroot,date);
    enc = erarawname(emod,epar,ereg,date);
    eds = erancread(enc,ebase);
    if haskey(eds,epar["ID"]); ID = epar["ID"]; else; ID = epar["IDnc"]; end
    return eds,eds[ID]

end

function eraanaread(
    ID::AbstractString,
    emod::Dict, epar::Dict, ereg::Dict, eroot::Dict,
    date::TimeType
)

    ebase = eraanafolder(epar,ereg,eroot);
    enc = eraananame(emod,epar,ereg,date);
    eds = erancread(enc,ebase);
    return eds,eds[ID]

end

function eracmpread(ID::AbstractString,emod::Dict,epar::Dict,ereg::Dict,eroot::Dict)

    ebase = eraanafolder(epar,ereg,eroot);
    enc = eracmpname(emod,epar,ereg);
    eds = erancread(enc,ebase);
    return eds,eds[ID]

end

function erancoffsetscale(data::Array{<:Real})

    dmax = maximum(data); dmin = minimum(data);
    scale = (dmax-dmin) / 65533;
    offset = (dmax+dmin-scale) / 2;

    return scale,offset

end

function erarawsave(
    data::Array{<:Real,3},
    emod::Dict, epar::Dict, ereg::Dict, date::TimeType, eroot::Dict
)

    @info "$(Dates.now()) - Saving raw $(uppercase(emod["dataset"])) $(epar["name"]) data in $(ereg["name"]) (Horizontal Resolution: $(ereg["step"])) for $(year(date)) $(Dates.monthname(date)) ..."

    fnc = joinpath(erarawfolder(epar,ereg,eroot,date),erarawname(emod,epar,ereg,date));
    if isfile(fnc)
        @info "$(Dates.now()) - Stale NetCDF file $(fnc) detected.  Overwriting ..."
        rm(fnc);
    end
    ds = NCDataset(fnc,"c",attrib = Dict(
        "Conventions" => "CF-1.6",
        "history"     => "Created on $(Dates.now())"
    ));

    ehr = hrindy(emod); nhr = ehr * daysinmonth(date);
    scale,offset = erancoffsetscale(data);

    ds.dim["longitude"] = ereg["size"][1];
    ds.dim["latitude"] = ereg["size"][2];
    ds.dim["time"] = nhr

    nclongitude = defVar(ds,"longitude",Float32,("longitude",),attrib = Dict(
        "units"     => "degrees_east",
        "long_name" => "longitude",
    ))

    nclatitude = defVar(ds,"latitude",Float32,("latitude",),attrib = Dict(
        "units"     => "degrees_north",
        "long_name" => "latitude",
    ))

    nctime = defVar(ds,"time",Int32,("time",),attrib = Dict(
        "units"     => "hours since $(date) 00:00:00.0",
        "long_name" => "time",
        "calendar"  => "gregorian",
    ))

    ncvar = defVar(ds,epar["ID"],Int16,("longitude","latitude","time"),attrib = Dict(
        "long_name"     => epar["era5"],
        "full_name"     => epar["name"],
        "units"         => epar["unit"],
        "scale_factor"  => scale,
        "add_offset"    => offset,
        "_FillValue"    => Int16(-32767),
        "missing_value" => Int16(-32767),
    ))

    nclongitude[:] = ereg["lon"]; nclatitude[:] = ereg["lat"]
    nctime[:] = (collect(1:nhr).-1) * ehr; ncvar[:] = data;

    close(ds)

    @info "$(Dates.now()) - Raw $(uppercase(emod["dataset"])) $(epar["name"]) in $(ereg["name"]) (Horizontal Resolution: $(ereg["step"])) for $(year(date)) $(Dates.monthname(date)) has been saved into $(fnc)."

end

function erasubregion!(
    emod::Dict, epar::Dict, ereg::Dict, etime::Dict, eroot::Dict,
    preg::Dict
)

    elon,elat,reginfo = gregiongridvec(ereg["region"],preg["lon"],preg["lat"])
    ereg["lon"]  = elon
    ereg["lat"]  = elat
    ereg["step"] = preg["step"]
    ereg["size"] = [length(elon),length(elat)]
    ereg["fol"]  = "$(ereg["region"])x$(@sprintf("%.2f",ereg["step"]))"

    for yr = etime["Begin"] : etime["End"], mo = 1:12

        @info "$(Dates.now()) - Extracting $(uppercase(emod["dataset"])) $(epar["name"]) data in $(gregionfullname(ereg["region"])) (Horizontal Resolution: $(ereg["step"])) from $(gregionfullname(preg["region"])) for $yr $(Dates.monthname(mo)) ..."

        date  = Date(yr,mo); nt = daysinmonth(date) * hrindy(emod)
        pdata = Array{Float32,2}(undef,preg["size"][1],preg["size"][2])
        tdata = Array{Float32,2}(undef,preg["size"][1],preg["size"][2])
        edata = Array{Float32,3}(undef,ereg["size"][1],ereg["size"][2],nt)
        pds,pvar = erarawread(emod,epar,preg,eroot,date);

        for it = 1 : nt
            pdata = pvar[:,:,it]*1;
            edata[:,:,it] .= regionextractgrid!(pdata,reginfo,tdata);
        end

        close(pds);
        erarawsave(edata,emod,epar,ereg,date,eroot);

    end

    return

end

function erasubregion(
    init::AbstractDict, eroot::AbstractDict;
    modID::AbstractString, parID::AbstractString,
    regID::AbstractString, timeID::Union{Integer,Vector}=0,
    gres::Real=0,
    pregID::AbstractString=""
)

    if pregID == ""; pregID = gregionparent(regID); end

    _,_,preg,_ = erainitialize(
        init;
        modID=modID,parID=parID,regID=pregID,timeID=timeID,
        gres=gres
    );

    emod,epar,ereg,etime = erainitialize(
        init;
        modID=modID,parID=parID,regID=regID,timeID=timeID,
        gres=gres
    );

    erasubregion!(emod,epar,ereg,etime,eroot,preg)
    putinfo(emod,epar,ereg,etime,eroot);

end

function putinfo(emod::Dict,epar::Dict,ereg::Dict,etime::Dict,eroot::Dict)

    rfol = pwd(); efol = erafolder(emod,epar,ereg,etime,eroot);

    @info "$(Dates.now()) - Saving emodule, eparameter, and GeoRegion information ..."

    cd(efol["var"]); @save "info_par.jld2" emod epar;
    cd(efol["reg"]); @save "info_reg.jld2" ereg;
    cd(rfol);

end


text:
export SymReal

struct SymReal <: Real
    ex::Expression
end

SymReal(name::Symbol) = SymReal(Variable(name))
SymReal(x::SymReal) = x
SymReal(x::Real) = SymReal(Numeric(x))
SymReal(::Irrational{sym}) where sym = SymReal(Constant(sym))
DOMAIN_TYPES[:Real] = SymReal

term(x::SymReal) = x.ex
Base.show(io::IO, t::SymReal) = print(io, t.ex)

for (M, f, arity) in DiffRules.diffrules()
    op = :($M.$f)

    if arity == 1
        @eval $op(x::SymReal) = track($op, x)
    elseif arity == 2
        @eval $op(x::SymReal, y::SymReal) = track($op, x, y)
    end
end

Base.promote_rule(::Type{<:SymReal}, ::Type{<:Real}) = SymReal
Base.convert(::Type{<:SymReal}, x::SymReal) = x
Base.convert(::Type{<:SymReal}, x::Real) = SymReal(x)
Base.convert(::Type{<:SymReal}, x::Irrational) = SymReal(x)


text:
module SCAConst

export PenConflict, PenMinSep, PenCloseness, PenAction, StdDist, MinSepSq, InvVar
export DT, DTI, G, Xmin, Xmax, Ymin, Ymax, Bearingmin, Bearingmax, Speedmin, Speedmax
export Xdim, Ydim, Bearingdim, Speeddim, COCdim, NStates, NActions
export Xs, Ys, Bearings, Speeds, Actions
export SigmaSpeed, SigmaBank, SigmaBankCOC
export SigmaDim, SigmaWeightNominal, SigmaWeightOffNominal

const PenConflict = 1.0
const PenMinSep = 1000.0
const PenCloseness = 10.0
const PenAction = 0.02

const StdDist = 500.0  # [m]
const MinSepSq = StdDist^2  # [m^2]
const InvVar = 1 / MinSepSq  # [1/m^2]

const DT = 5.0  # [s]
const DTI = 1.0  # [s]
const G = 9.8  # [m/s^2]

const Xmin = -2000.0  # [m]
const Xmax = 2000.0  # [m]
const Ymin = -2000.0  # [m]
const Ymax = 2000.0  # [m]
const Bearingmin = 0.0  # [rad]
const Bearingmax = 2 * pi  # [rad]
const Speedmin = 10  # [m/s]
const Speedmax = 20  # [m/s]

const Xdim = 11  # 51
const Ydim = 11  # 51
const Bearingdim = 5  # 37
const Speeddim = 3  # 3

const NStates = Xdim * Ydim * Bearingdim * Speeddim^2 + 1
const NActions = 36

const Xs = linspace(Xmin, Xmax, Xdim)
const Ys = linspace(Ymin, Ymax, Ydim)
const Bearings = linspace(Bearingmin, Bearingmax, Bearingdim)
const Speeds = linspace(Speedmin, Speedmax, Speeddim)

const Actions = [:right20, :right10, :straight, :left10, :left20, :clearOfConflict]

const SigmaSpeed = 2.0  # [m/s]
const SigmaBank = deg2rad(4.0)  # [rad]
const SigmaBankCOC = deg2rad(10.0)  # [rad]

const SigmaDim = 4  # number of dimensions for sigma-point sampling
const SigmaWeightNominal = 1 / 3
const SigmaWeightOffNominal = (1 - SigmaWeightNominal) / (2 * SigmaDim)

end # module SCAConst


text:
using ExtensibleUnions
using Test

# Most of the code in this file is taken from:
# 1. https://github.com/NHDaly/DeepcopyModules.jl (license: MIT)
# 2. https://github.com/perrutquist/CodeTransformation.jl (license: MIT)

@test ExtensibleUnions.getmodule(typeof(sin)) === Base
@test ExtensibleUnions.getmodule(sin) === Base

let
    # Test example from docstring to ExtensibleUnions.addmethod!
    g(x) = x + 13
    ci = code_lowered(g)[1]
    function f end
    ExtensibleUnions.addmethod!(f, (Any,), ci)
    @test f(1) === 14

    # Alternative syntax
    function f2 end
    @test ExtensibleUnions.makesig(f2, (Any,)) === Tuple{typeof(f2), Any}
    ExtensibleUnions.addmethod!(Tuple{typeof(f2), Any}, ci)
    @test f2(1) === 14

end

let
    # Test example from docstring to codetransform!
    g(x) = x + 13
    function e end
    ExtensibleUnions.codetransform!(g => e) do ci
        for ex in ci.code
            if ex isa Expr
                map!(x -> x === 13 ? 7 : x, ex.args, ex.args)
            end
        end
        ci
    end
    @test e(1) === 8
    @test g(1) === 14
end

let
    a = Vector{T} where T
    b = ExtensibleUnions.typevars(a)
    @test b isa Tuple
    @test length(b) == 1
    @test b[1] isa TypeVar
    @test b[1].name == :T
    @test b[1].lb === Union{}
    @test b[1].ub === Any
end


text:
abstract type EllipticalCopula{d,MT} <: Copula{d} end
# N(::Type{EllipticalCopula) = @error "Not Implemented"
# U(::EllipticalCopula) = @error "Not Implemented"
Base.eltype(C::CT) where CT<:EllipticalCopula = Base.eltype(N(CT)(C.Σ))
function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT <: EllipticalCopula}
    Random.rand!(rng,N(CT)(C.Σ),x)
    x .= cdf.(U(CT)(),x)
end
function Base.rand(rng::Distributions.AbstractRNG,C::CT) where {CT <: EllipticalCopula}
    x = Random.rand(rng,N(CT)(C.Σ))
    x .= cdf.(U(CT)(),x)
end
function Distributions._logpdf(C::CT, u) where {CT <: EllipticalCopula}
    x = quantile.(U(CT)(),u)
    return Distributions.logpdf(N(CT)(C.Σ),x) - sum(Distributions.logpdf.(U(CT)(),x))
end
function make_cor!(Σ)
    # Verify that Σ is a correlation matrix, otherwise make it so :
    d = size(Σ,1)
    σ = [1/sqrt(Σ[i,i]) for i in 1:d]
    for i in 1:d
        for j in 1:d
            Σ[i,j] *= σ[i] .* σ[j]
        end
    end
end