text:
"""
    censored(d0::UnivariateDistribution; [lower::Real], [upper::Real])
    censored(d0::UnivariateDistribution, lower::Real, upper::Real)

A _censored distribution_ `d` of a distribution `d0` to the interval
``[l, u]=```[lower, upper]` has the probability density (mass) function:

```math
f(x; d_0, l, u) = \\begin{cases}
    P_{Z \\sim d_0}(Z \\le l), & x = l \\\\
    f_{d_0}(x),              & l < x < u \\\\
    P_{Z \\sim d_0}(Z \\ge u), & x = u \\\\
  \\end{cases}, \\quad x \\in [l, u]
```
where ``f_{d_0}(x)`` is the probability density (mass) function of ``d_0``.

If ``Z \\sim d_0``, and `X = clamp(Z, l, u)`, then ``X \\sim d``. Note that this implies
that even if ``d_0`` is continuous, its censored form assigns positive probability to the
bounds ``l`` and ``u``. Therefore, a censored continuous distribution has atoms and is a
mixture of discrete and continuous components.

The function falls back to constructing a [`Distributions.Censored`](@ref) wrapper.

# Usage

```julia
censored(d0; lower=l)           # d0 left-censored to the interval [l, Inf)
censored(d0; upper=u)           # d0 right-censored to the interval (-Inf, u]
censored(d0; lower=l, upper=u)  # d0 interval-censored to the interval [l, u]
censored(d0, l, u)              # d0 interval-censored to the interval [l, u]
```

# Implementation

To implement a specialized censored form for distributions of type `D`, instead of
overloading a method with one of the above signatures, one or more of the following methods
should be implemented:
- `censored(d0::D, l::T, u::T) where {T <: Real}`
- `censored(d0::D, ::Nothing, u::Real)`
- `censored(d0::D, l::Real, ::Nothing)`
"""
function censored end
function censored(d0::UnivariateDistribution, l::T, u::T) where {T<:Real}
    return Censored(d0, l, u)
end
function censored(d0::UnivariateDistribution, ::Nothing, u::Real)
    return Censored(d0, nothing, u)
end
function censored(d0::UnivariateDistribution, l::Real, ::Nothing)
    return Censored(d0, l, nothing)
end
censored(d0::UnivariateDistribution, l::Real, u::Real) = censored(d0, promote(l, u)...)
censored(d0::UnivariateDistribution, ::Nothing, ::Nothing) = d0
function censored(
    d0::UnivariateDistribution;
    lower::Union{Real,Nothing} = nothing,
    upper::Union{Real,Nothing} = nothing,
)
    return censored(d0, lower, upper)
end

"""
    Censored

Generic wrapper for a [`censored`](@ref) distribution.
"""
struct Censored{
    D<:UnivariateDistribution,
    S<:ValueSupport,
    T<:Real,
    TL<:Union{T,Nothing},
    TU<:Union{T,Nothing},
} <: UnivariateDistribution{S}
    uncensored::D      # the original distribution (uncensored)
    lower::TL      # lower bound
    upper::TU      # upper bound
    function Censored(d0::UnivariateDistribution, lower::T, upper::T; check_args::Bool=true) where {T<:Real}
        @check_args(Censored, lower ≤ upper)
        new{typeof(d0), value_support(typeof(d0)), T, T, T}(d0, lower, upper)
    end
    function Censored(d0::UnivariateDistribution, l::Nothing, u::Real; check_args::Bool=true)
        new{typeof(d0), value_support(typeof(d0)), typeof(u), Nothing, typeof(u)}(d0, l, u)
    end
    function Censored(d0::UnivariateDistribution, l::Real, u::Nothing; check_args::Bool=true)
        new{typeof(d0), value_support(typeof(d0)), typeof(l), typeof(l), Nothing}(d0, l, u)
    end
end

const LeftCensored{D<:UnivariateDistribution,S<:ValueSupport,T<:Real} = Censored{D,S,T,T,Nothing}
const RightCensored{D<:UnivariateDistribution,S<:ValueSupport,T<:Real} = Censored{D,S,T,Nothing,T}

function censored(d::Censored, l::T, u::T) where {T<:Real}
    return censored(
        d.uncensored,
        d.lower === nothing ? l : max(l, d.lower),
        d.upper === nothing ? u : min(u, d.upper),
    )
end
function censored(d::Censored, ::Nothing, u::Real)
    return censored(d.uncensored, d.lower, d.upper === nothing ? u : min(u, d.upper))
end
function censored(d::Censored, l::Real, ::Nothing)
    return censored(d.uncensored, d.lower === nothing ? l : max(l, d.lower), d.upper)
end

function params(d::Censored)
    d0params = params(d.uncensored)
    return (d0params..., d.lower, d.upper)
end

function partype(d::Censored{<:UnivariateDistribution,<:ValueSupport,T}) where {T}
    return promote_type(partype(d.uncensored), T)
end

Base.eltype(::Type{<:Censored{D,S,T}}) where {D,S,T} = promote_type(T, eltype(D))

#### Range and Support

isupperbounded(d::LeftCensored) = isupperbounded(d.uncensored)
isupperbounded(d::Censored) = isupperbounded(d.uncensored) || _ccdf_inclusive(d.uncensored, d.upper) > 0

islowerbounded(d::RightCensored) = islowerbounded(d.uncensored)
islowerbounded(d::Censored) = islowerbounded(d.uncensored) || cdf(d.uncensored, d.lower) > 0

maximum(d::LeftCensored) = max(maximum(d.uncensored), d.lower)
maximum(d::Censored) = min(maximum(d.uncensored), d.upper)

minimum(d::RightCensored) = min(minimum(d.uncensored), d.upper)
minimum(d::Censored) = max(minimum(d.uncensored), d.lower)

function insupport(d::Censored, x::Real)
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    return (
        (_in_open_interval(x, lower, upper) && insupport(d0, x)) ||
        (x == lower && cdf(d0, lower) > 0) ||
        (x == upper && _ccdf_inclusive(d0, upper) > 0)
    )
end

#### Show

function show(io::IO, ::MIME"text/plain", d::Censored)
    print(io, "Censored(")
    d0 = d.uncensored
    uml, namevals = _use_multline_show(d0)
    uml ? show_multline(io, d0, namevals; newline=false) : show_oneline(io, d0, namevals)
    if d.lower === nothing
        print(io, "; upper=$(d.upper))")
    elseif d.upper === nothing
        print(io, "; lower=$(d.lower))")
    else
        print(io, "; lower=$(d.lower), upper=$(d.upper))")
    end
end

_use_multline_show(d::Censored) = _use_multline_show(d.uncensored)

#### Statistics

quantile(d::Censored, p::Real) = _clamp(quantile(d.uncensored, p), d.lower, d.upper)

median(d::Censored) = _clamp(median(d.uncensored), d.lower, d.upper)

# the expectations use the following relation:
# 𝔼_{x ~ d}[h(x)] = P_{x ~ d₀}(x < l) h(l) + P_{x ~ d₀}(x > u) h(u)
#                 + P_{x ~ d₀}(l ≤ x ≤ u) 𝔼_{x ~ τ}[h(x)],
# where d₀ is the uncensored distribution, d is d₀ censored to [l, u],
# and τ is d₀ truncated to [l, u]

function mean(d::LeftCensored)
    lower = d.lower
    log_prob_lower = _logcdf_noninclusive(d.uncensored, lower)
    log_prob_interval = log1mexp(log_prob_lower)
    μ = xexpy(lower, log_prob_lower) + xexpy(mean(_to_truncated(d)), log_prob_interval)
    return μ
end
function mean(d::RightCensored)
    upper = d.upper
    log_prob_upper = logccdf(d.uncensored, upper)
    log_prob_interval = log1mexp(log_prob_upper)
    μ = xexpy(upper, log_prob_upper) + xexpy(mean(_to_truncated(d)), log_prob_interval)
    return μ
end
function mean(d::Censored)
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    log_prob_lower = _logcdf_noninclusive(d0, lower)
    log_prob_upper = logccdf(d0, upper)
    log_prob_interval = log1mexp(logaddexp(log_prob_lower, log_prob_upper))
    μ = (xexpy(lower, log_prob_lower) + xexpy(upper, log_prob_upper) +
         xexpy(mean(_to_truncated(d)), log_prob_interval))
    return μ
end

function var(d::LeftCensored)
    lower = d.lower
    log_prob_lower = _logcdf_noninclusive(d.uncensored, lower)
    log_prob_interval = log1mexp(log_prob_lower)
    dtrunc = _to_truncated(d)
    μ_interval = mean(dtrunc)
    μ = xexpy(lower, log_prob_lower) + xexpy(μ_interval, log_prob_interval)
    v_interval = var(dtrunc) + abs2(μ_interval - μ)
    v = xexpy(abs2(lower - μ), log_prob_lower) + xexpy(v_interval, log_prob_interval)
    return v
end
function var(d::RightCensored)
    upper = d.upper
    log_prob_upper = logccdf(d.uncensored, upper)
    log_prob_interval = log1mexp(log_prob_upper)
    dtrunc = _to_truncated(d)
    μ_interval = mean(dtrunc)
    μ = xexpy(upper, log_prob_upper) + xexpy(μ_interval, log_prob_interval)
    v_interval = var(dtrunc) + abs2(μ_interval - μ)
    v = xexpy(abs2(upper - μ), log_prob_upper) + xexpy(v_interval, log_prob_interval)
    return v
end
function var(d::Censored)
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    log_prob_lower = _logcdf_noninclusive(d0, lower)
    log_prob_upper = logccdf(d0, upper)
    log_prob_interval = log1mexp(logaddexp(log_prob_lower, log_prob_upper))
    dtrunc = _to_truncated(d)
    μ_interval = mean(dtrunc)
    μ = (xexpy(lower, log_prob_lower) + xexpy(upper, log_prob_upper) +
         xexpy(μ_interval, log_prob_interval))
    v_interval = var(dtrunc) + abs2(μ_interval - μ)
    v = (xexpy(abs2(lower - μ), log_prob_lower) + xexpy(abs2(upper - μ), log_prob_upper) +
         xexpy(v_interval, log_prob_interval))
    return v
end

# this expectation also uses the following relation:
# 𝔼_{x ~ τ}[-log d(x)] = H[τ] - log P_{x ~ d₀}(l ≤ x ≤ u)
#   + (P_{x ~ d₀}(x = l) (log P_{x ~ d₀}(x = l) - log P_{x ~ d₀}(x ≤ l)) +
#      P_{x ~ d₀}(x = u) (log P_{x ~ d₀}(x = u) - log P_{x ~ d₀}(x ≥ u))
#   ) / P_{x ~ d₀}(l ≤ x ≤ u),
# where H[τ] is the entropy of τ.

function entropy(d::LeftCensored)
    d0 = d.uncensored
    lower = d.lower
    log_prob_lower_inc = logcdf(d0, lower)
    if value_support(typeof(d0)) === Discrete
        logpl = logpdf(d0, lower)
        log_prob_lower = logsubexp(log_prob_lower_inc, logpl)
        xlogx_pl = xexpx(logpl)
    else
        log_prob_lower = log_prob_lower_inc
        xlogx_pl = 0
    end
    log_prob_interval = log1mexp(log_prob_lower)
    entropy_bound = -xexpx(log_prob_lower_inc)
    dtrunc = _to_truncated(d)
    entropy_interval = xexpy(entropy(dtrunc), log_prob_interval) - xexpx(log_prob_interval) + xlogx_pl
    return entropy_interval + entropy_bound
end
function entropy(d::RightCensored)
    d0 = d.uncensored
    upper = d.upper
    log_prob_upper = logccdf(d0, upper)
    if value_support(typeof(d0)) === Discrete
        logpu = logpdf(d0, upper)
        log_prob_upper_inc = logaddexp(log_prob_upper, logpu)
        xlogx_pu = xexpx(logpu)
    else
        log_prob_upper_inc = log_prob_upper
        xlogx_pu = 0
    end
    log_prob_interval = log1mexp(log_prob_upper)
    entropy_bound = -xexpx(log_prob_upper_inc)
    dtrunc = _to_truncated(d)
    entropy_interval = xexpy(entropy(dtrunc), log_prob_interval) - xexpx(log_prob_interval) + xlogx_pu
    return entropy_interval + entropy_bound
end
function entropy(d::Censored)
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    log_prob_lower_inc = logcdf(d0, lower)
    log_prob_upper = logccdf(d0, upper)
    if value_support(typeof(d0)) === Discrete
        logpl = logpdf(d0, lower)
        logpu = logpdf(d0, upper)
        log_prob_lower = logsubexp(log_prob_lower_inc, logpl)
        log_prob_upper_inc = logaddexp(log_prob_upper, logpu)
        xlogx_pl = xexpx(logpl)
        xlogx_pu = xexpx(logpu)
    else
        log_prob_lower = log_prob_lower_inc
        log_prob_upper_inc = log_prob_upper
        xlogx_pl = xlogx_pu = 0
    end
    log_prob_interval = log1mexp(logaddexp(log_prob_lower, log_prob_upper))
    entropy_bound = -(xexpx(log_prob_lower_inc) + xexpx(log_prob_upper_inc))
    dtrunc = _to_truncated(d)
    entropy_interval = xexpy(entropy(dtrunc), log_prob_interval) - xexpx(log_prob_interval) + xlogx_pl + xlogx_pu
    return entropy_interval + entropy_bound
end

#### Evaluation

function pdf(d::Censored, x::Real)
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    px = float(pdf(d0, x))
    return if _in_open_interval(x, lower, upper)
        px
    elseif x == lower
        x == upper ? one(px) : oftype(px, cdf(d0, x))
    elseif x == upper
        if value_support(typeof(d0)) === Discrete
            oftype(px, ccdf(d0, x) + px)
        else
            oftype(px, ccdf(d0, x))
        end
    else  # not in support
        zero(px)
    end
end

function logpdf(d::Censored, x::Real)
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    logpx = logpdf(d0, x)
    return if _in_open_interval(x, lower, upper)
        logpx
    elseif x == lower
        x == upper ? zero(logpx) : oftype(logpx, logcdf(d0, x))
    elseif x == upper
        if value_support(typeof(d0)) === Discrete
            oftype(logpx, logaddexp(logccdf(d0, x), logpx))
        else
            oftype(logpx, logccdf(d0, x))
        end
    else  # not in support
        oftype(logpx, -Inf)
    end
end

function loglikelihood(d::Censored, x::AbstractArray{<:Real})
    d0 = d.uncensored
    lower = d.lower
    upper = d.upper
    logpx = logpdf(d0, first(x))
    log_prob_lower = lower === nothing ? zero(logpx) : oftype(logpx, logcdf(d0, lower))
    log_prob_upper = upper === nothing ? zero(logpx) : oftype(logpx, _logccdf_inclusive(d0, upper))
    logzero = oftype(logpx, -Inf)
    return sum(x) do xi
        _in_open_interval(xi, lower, upper) && return logpdf(d0, xi)
        xi == lower && return log_prob_lower
        xi == upper && return log_prob_upper
        return logzero
    end
end

function cdf(d::Censored, x::Real)
    lower = d.lower
    upper = d.upper
    result = cdf(d.uncensored, x)
    return if lower !== nothing && x < lower
        zero(result)
    elseif upper === nothing || x < upper
        result
    else
        one(result)
    end
end

function logcdf(d::Censored, x::Real)
    lower = d.lower
    upper = d.upper
    result = logcdf(d.uncensored, x)
    return if d.lower !== nothing && x < d.lower
        oftype(result, -Inf)
    elseif d.upper === nothing || x < d.upper
        result
    else
        zero(result)
    end
end

function ccdf(d::Censored, x::Real)
    lower = d.lower
    upper = d.upper
    result = ccdf(d.uncensored, x)
    return if lower !== nothing && x < lower
        one(result)
    elseif upper === nothing || x < upper
        result
    else
        zero(result)
    end
end

function logccdf(d::Censored{<:Any,<:Any,T}, x::Real) where {T}
    lower = d.lower
    upper = d.upper
    result = logccdf(d.uncensored, x)
    return if lower !== nothing && x < lower
        zero(result)
    elseif upper === nothing || x < upper
        result
    else
        oftype(result, -Inf)
    end
end

#### Sampling

rand(rng::AbstractRNG, d::Censored) = _clamp(rand(rng, d.uncensored), d.lower, d.upper)

#### Utilities

# utilities to handle intervals represented with possibly `nothing` bounds

_in_open_interval(x::Real, l::Real, u::Real) = l < x < u
_in_open_interval(x::Real, ::Nothing, u::Real) = x < u
_in_open_interval(x::Real, l::Real, ::Nothing) = x > l

_clamp(x, l, u) = clamp(x, l, u)
_clamp(x, ::Nothing, u) = min(x, u)
_clamp(x, l, ::Nothing) = max(x, l)

_to_truncated(d::Censored) = truncated(d.uncensored, d.lower, d.upper)

# utilities for non-inclusive CDF p(x < u) and inclusive CCDF (p ≥ u)

_logcdf_noninclusive(d::UnivariateDistribution, x) = logcdf(d, x)
function _logcdf_noninclusive(d::DiscreteUnivariateDistribution, x)
    return logsubexp(logcdf(d, x), logpdf(d, x))
end

_ccdf_inclusive(d::UnivariateDistribution, x) = ccdf(d, x)
_ccdf_inclusive(d::DiscreteUnivariateDistribution, x) = ccdf(d, x) + pdf(d, x)

_logccdf_inclusive(d::UnivariateDistribution, x) = logccdf(d, x)
function _logccdf_inclusive(d::DiscreteUnivariateDistribution, x)
    return logaddexp(logccdf(d, x), logpdf(d, x))
end

# like xlogx but for input on log scale, safe when x == -Inf
function xexpx(x::Real)
    result = x * exp(x)
    return x == -Inf ? zero(result) : result
end

# x * exp(y) with correct limit for y == -Inf
function xexpy(x::Real, y::Real)
    result = x * exp(y)
    return y == -Inf && !isnan(x) ? zero(result) : result
end


text:
const rchsVersion    = "CRINEX VERS   / TYPE"
const rchsProg       = "CRINEX PROG / DATE"
# Nav
const nhsVersion     = "RINEX VERSION / TYPE"
const nhsRunBy       = "PGM / RUN BY / DATE"
const nhsComment     = "COMMENT"
const nhsIonoCorr    = "IONOSPHERIC CORR"
const nhsTimeSysCorr = "TIME SYSTEM CORR"
const nhsLeapSeconds = "LEAP SECONDS"
    # R2.10GLO
const nhsCorrSysTime = "CORR TO SYSTEM TIME"
    # R2.11GPS
const nhsDeltaUTC    = "DELTA-UTC: A0,A1,T,W"
    # R2.11GEO
const nhsDUTC        = "D-UTC A0,A1,T,W,S,U"
    # R2.11
const nhsIonAlpha    = "ION ALPHA"
    # R2.11
const nhsIonBeta     = "ION BETA"
const nhsEoH         = "END OF HEADER"

# Obs
const ohsVersion           = nhsVersion
const ohsRunBy             = nhsRunBy
const ohsComment           = nhsComment
const ohsMarkerName        = "MARKER NAME"
const ohsMarkerNumber      = "MARKER NUMBER"
const ohsMarkerType        = "MARKER TYPE"
const ohsObserver          = "OBSERVER / AGENCY"
const ohsReceiver          = "REC # / TYPE / VERS"
const ohsAntennaType       = "ANT # / TYPE"
const ohsAntennaPosition   = "APPROX POSITION XYZ"
const ohsAntennaDeltaHEN   = "ANTENNA: DELTA H/E/N"
const ohsAntennaDeltaXYZ   = "ANTENNA: DELTA X/Y/Z"
const ohsAntennaPhaseCtr   = "ANTENNA: PHASECENTER"
const ohsAntennaBsightXYZ  = "ANTENNA: B.SIGHT XYZ"
const ohsAntennaZeroDirAzi = "ANTENNA: ZERODIR AZI"
const ohsAntennaZeroDirXYZ = "ANTENNA: ZERODIR XYZ"
const ohsCenterOfMass      = "CENTER OF MASS: XYZ"
const ohsNumObs            = "# / TYPES OF OBSERV"
const ohsSystemNumObs      = "SYS / # / OBS TYPES"
const ohsWaveFact          = "WAVELENGTH FACT L1/2"
const ohsSigStrengthUnit   = "SIGNAL STRENGTH UNIT"
const ohsInterval          = "INTERVAL"
const ohsFirstTime         = "TIME OF FIRST OBS"
const ohsLastTime          = "TIME OF LAST OBS"
const ohsReceiverOffset    = "RCV CLOCK OFFS APPL"
const ohsSystemDCBSapplied = "SYS / DCBS APPLIED"
const ohsSystemPCVSapplied = "SYS / PCVS APPLIED"
const ohsSystemScaleFac    = "SYS / SCALE FACTOR"
const ohsSystemPhaseShift  = "SYS / PHASE SHIFT"
const ohsGlonassSlotFreqNo = "GLONASS SLOT / FRQ #"
const ohsGlonassCodPhsBias = "GLONASS COD/PHS/BIS"
const ohsLeapSeconds       = nhsLeapSeconds
const ohsNumSats           = "# OF SATELLITES"
const ohsPrnObs            = "PRN / # OF OBS"
const ohsEoH               = nhsEoH

function get_header!(header_hs::Dict{String,String},
    f::IOStream,
    hstrings_set::Set{String},
    end_label::String=ohsEoH)::Bool

    while !eof(f)
        line = rstrip(readline(f))
        label = line[61:end]
        if !in(label, hstrings_set)
            @warn "don't have such label in rinex |$label|"
        end
        if haskey(header_hs, label)
            header_hs[label] = string(header_hs[label], "\n", line[1:60], )
        else
            header_hs[label] = line[1:60]
        end
        if isequal(label, end_label)
            return true
        end
    end
    throw(error("file no $end_label"))
    false
end

function returninput(input::String)::String
    return input
end

#####################################
function parse_rchsVersion(input::String)::String
    return input
end
function parse_rchsProg(input::String)::String
    return input
end
#####################################
# nav header parse start
# begin

function setFileSystem(c::Char, v::Float64)
    if c == "M" && v < 3.0
        throw(error("rinex 2 'Mixed' Nav files do not existed"))
    end
    char2SatelliteSystemType(c)
end

function parse_nhsVersion(input::String)
    version  = parse(Float64, input[ 1:20])
    fileType = input[21:40]
    ftc  = uppercase(fileType[1])
    fileSys  = input[41:60]

    if version > 2.99
        if ftc != 'N'
            throw(error("File type is not Nav $fileType"))
        end
        ssc = uppercase(fileSys[1])
    else
        if ftc == 'N'     ssc = 'G'
        elseif ftc == 'G' ssc = 'R'
            elseif ftc == 'H' ssc = 'S'
            else
                throw(error("version 2 file type is invalid $fileType"))
        end
    end
    ss = setFileSystem(ssc, version)

        # fileType = strip(input[21:40])
        # fileSys  = strip(input[41:60])
    return version, fileType, fileSys, ss
end
    # Nav
function parse_nhsRunBy(input::String)
    return strip(input[1:20]), strip(input[21:40]), strip(input[41:60])
end
parse_nhsComment(input::String)     = returninput(input)
parse_nhsIonoCorr(input::String)    = returninput(input)
function parse_nhsTimeSysCorr(input::String)
    len = length(input) ÷ 60
    arr = Dict{TimeCorrectionType,TimeSystemCorrection}()
    for i in 1:len
        start_index = (i - 1) * 61 + 1
        line = SubString(input, start_index, start_index + 59)
        t = string(line[1:4])
        tp = str2TimeCorrectionType(t)
        A0 = parse(Float64, line[6:22])
        A1 = parse(Float64, line[23:38])
        refSOW = parse(Int64, line[39:45])
        refWeek = parse(Int64, line[46:51])
        geoProvider = string(line[52:57])
        geoUTCid = string(line[58:60])
        arr[tp]  = getTimeSystemCorrection(tp, A0, A1,
                                        refSOW, refWeek,
                                        geoProvider, geoUTCid)
    end
    arr
end

function parse_lpn(s::SubString)
    a = strip(s)
    if a == ""
        return 0
    else
        return parse(Int64, a)
    end
end
function parse_nhsLeapSeconds(line::String)
    if line == "" return 0, 0, 0, 0 end
    leapSeconds = parse(Int64, SubString(line, 1, 6))
    leapDelta = parse_lpn(SubString(line, 7, 12))
    leapWeek  = parse_lpn(SubString(line, 13, 18))
    leapDay   = parse_lpn(SubString(line, 19, 20))
    leapSeconds, leapDelta, leapWeek, leapDay
end
    # R2.10GLO
function parse_nhsCorrSysTime(line::String)
    arr = Dict{TimeCorrectionType,TimeSystemCorrection}()
        # line =  "  2011    09    08    1.578591763973e-07                    CORR TO SYSTEM TIME"
    tp = str2TimeCorrectionType("GLUT")
    A0 = parse(Float64, replace(line[22:40], 'D' => 'E'))
    A1 = 0.0
    refYr  = parse(Int64, line[1:6])
    refMon = parse(Int64, line[7:12])
    refDay = parse(Int64, line[13:18])
    gldt = DateTime(refYr, refMon, refDay, 0, 0, 0)

    utdt = gldt + diff2UTC(gldt, tstGLO)

    refWeek, refSOW = HTgnss.datetime2gpsws(utdt)
        # @show gldt, utdt, refWeek, refSOW
    geoProvider = "    "
    geoUTCid = " 3"
    frTST, toTST = tct2TimeSystemTypes(tp)
    arr[tp] = TimeSystemCorrection(tp, frTST, toTST,
                        A0, A1, refWeek, refSOW,
                        refYr, refMon, refDay,
                        geoProvider, geoUTCid)
    arr
end
    # R2.11GPS
function parse_nhsDeltaUTC(line::String)
    arr = Dict{TimeCorrectionType,TimeSystemCorrection}()
        # line =  "   -0.186264514923D-08-0.106581410364D-13   503808     1652 DELTA-UTC: A0,A1,T,W"
    tp = str2TimeCorrectionType("GPUT")
    A0 = parse(Float64, replace(line[4:22], 'D' => 'E'))
    A1 = parse(Float64, replace(line[23:41], 'D' => 'E'))

    refSOW = parse(Int64, line[42:50])
    refWeek = parse(Int64, line[51:59])

    geoProvider = "    "
    geoUTCid = " 0"

    arr[tp]  = getTimeSystemCorrection(tp, A0, A1,
                            refSOW, refWeek,
                            geoProvider, geoUTCid)
    arr
end
    # R2.11GEO
function parse_nhsDUTC(line::String)
    arr = Dict{TimeCorrectionType,TimeSystemCorrection}()
        # line =
    tp = str2TimeCorrectionType("SBUT")
    A0 = parse(Float64, line[1:19])
    A1 = parse(Float64, line[20:38])

    refSOW = parse(Int64, line[38:45])
    refWeek = parse(Int64, line[46:51])

    geoProvider = line[52:57]
    geoUTCid = line[58:59]

    arr[tp]  = getTimeSystemCorrection(tp, A0, A1,
                            refSOW, refWeek,
                            geoProvider, geoUTCid)
    arr
end
    # R2.11
parse_nhsIonAlpha(input::String)    = returninput(input)
    # R2.11
parse_nhsIonBeta(input::String)     = returninput(input)
parse_nhsEoH(input::String)         = returninput(input)

# end

const nav_hs_function = Dict{String,Function}(
    rchsVersion    => parse_rchsVersion,
    rchsProg       => parse_rchsProg,
    nhsVersion     => parse_nhsVersion,
    nhsRunBy       => parse_nhsRunBy,
    nhsComment     => parse_nhsComment,
    nhsIonoCorr    => parse_nhsIonoCorr,
    nhsTimeSysCorr => parse_nhsTimeSysCorr,
    nhsLeapSeconds => parse_nhsLeapSeconds, # R2.10GLO
    nhsCorrSysTime => parse_nhsCorrSysTime, # R2.11GPS
    nhsDeltaUTC    => parse_nhsDeltaUTC,    # R2.11GEO
    nhsDUTC        => parse_nhsDUTC,        # R2.11
    nhsIonAlpha    => parse_nhsIonAlpha,    # R2.11
    nhsIonBeta     => parse_nhsIonBeta,
    nhsEoH         => parse_nhsEoH,

)

const nav_hstrings_set = Set(keys(nav_hs_function))
# nav header parse end

#####################################
# obs header parse begin
begin
    function parse_ohsVersion(input::String)
    version  = parse(Float64, input[ 1:20])
    fileType = uppercase(input[21])

    fileSys  = char2SatelliteSystemType(input[41])
        # fileType = strip(input[21:40])
        # fileSys  = strip(input[41:60])
    return version, fileType, fileSys
end

    parse_ohsRunBy(input::String) = parse_nhsRunBy(input)
    parse_ohsComment(input::String) = parse_nhsComment(input)
    parse_ohsMarkerName(input::String) = returninput(input)
    parse_ohsMarkerNumber(input::String) = returninput(input)
    parse_ohsMarkerType(input::String) = returninput(input)
    function parse_ohsObserver(input::String)
    return strip(input[1:20]), strip(input[21:60])
end
    parse_ohsReceiver(input::String) = returninput(input)
    parse_ohsAntennaType(input::String) = returninput(input)
    function parse_ohsAntennaPosition(input::String)
    x = parse(Float64, input[1:14])
    y = parse(Float64, input[15:28])
    z = parse(Float64, input[29:42])
    ECEF{Float64}(x, y, z)
end
    parse_ohsAntennaDeltaHEN(input::String) = returninput(input)
    parse_ohsAntennaDeltaXYZ(input::String) = returninput(input)
    parse_ohsAntennaPhaseCtr(input::String) = returninput(input)
    parse_ohsAntennaBsightXYZ(input::String) = returninput(input)
    parse_ohsAntennaZeroDirAzi(input::String) = returninput(input)
    parse_ohsAntennaZeroDirXYZ(input::String) = returninput(input)
    parse_ohsCenterOfMass(input::String) = returninput(input)
    parse_ohsNumObs(input::String) = returninput(input)
    function NumObs2SystemNumObs(obs_codes::Array{String,1},
            obsn::Int64, sys::SatelliteSystemType,
            obst::D_SystemNumObs)
    this_obst = Array{ObsID,1}(undef, obsn)
    for iot in 1:obsn
            # @show obs_codes[iot]
        this_obst[iot] = ObsID(obs_codes[iot], sys)
    end
    obst[sys] = obsn, this_obst
    return
end
    function parse_ohsSystemNumObs2(input::String, sys::SatelliteSystemType=sstMixed)

    obsn = parse(Int64, input[1:6])
    maxObsPerLine = 9
    len_max = obsn ÷ maxObsPerLine
    res     = obsn % maxObsPerLine
    len_max += res == 0 ? 0 : 1
        lines = split(input, "\n")
        length(lines) == len_max ||
                throw(error("line number error for $obsn!\n|$lines|"))
        i_lines = 1

        obs_codes = Array{String,1}()
        this_line = lines[i_lines]
        start_index = 5
        for iot in  1:obsn
            start_index += 6
            c = this_line[start_index:start_index + 1]
            if c[1] == 'C'     cc = string('C', c[2], 'C')
            elseif c[1] == 'P' cc = string('C', c[2], 'P')
            else cc = string(c, " ")
            end
            push!(obs_codes, cc)
            if iot % maxObsPerLine == 0 && i_lines < len_max
                i_lines += 1
                this_line = lines[i_lines]
                start_index = 5
            end
        end
        # @show obs_codes
        obst = D_SystemNumObs()
        if sys == sstMixed
            for isys in [sstGPS,sstGlonass,sstGalileo,
                            sstGeosync,sstBeiDou,sstQZSS,
                            sstIRNSS]
                NumObs2SystemNumObs(obs_codes, obsn, isys, obst)
            end
        else
            NumObs2SystemNumObs(obs_codes, obsn, sys, obst)
        end
        obst
    end
    function parse_ohsSystemNumObs3(input::String)
        lines = split(input, "\n")
        lines_n = length(lines) + 1
        i_lines = 1
        obst = D_SystemNumObs()
        while i_lines < lines_n
            this_line = lines[i_lines]
            i_lines += 1
            satSysC = this_line[1]
            satSys = char2SatelliteSystemType(satSysC)
            obsn = parse(Int64, this_line[3:6])
            this_obst = Array{ObsID,1}(undef, obsn)
            start_index = 4
            for iot in 1:obsn
                start_index += 4
                if start_index > 58
                    this_line = lines[i_lines]
                    i_lines +=  1
                    start_index = 4
                end
                this_obst[iot] =
                    ObsID(this_line[start_index:start_index + 2],
                            satSys)
            end
            obst[satSys] = (obsn, this_obst)
        end
        obst
    end

    parse_ohsWaveFact(input::String) = returninput(input)
    parse_ohsSigStrengthUnit(input::String) = returninput(input)
    function parse_ohsInterval(input::String)
        if input == "" return NaN end
        parse(Float64, input[1:10])
    end
    # todo get time system

    function parse_ohsFirstTime(input::String)
        # @show SubString(input, 3, 43)
        dtstr = SubString(input, 3, 43)
        y = parse(Int64, dtstr[1:4])
        m = parse(Int64, dtstr[9:10])
        d = parse(Int64, dtstr[15:16])
        H = parse(Int64, dtstr[21:22])
        M = parse(Int64, dtstr[27:28])
        S = parse(Int64, dtstr[32:33])
        s = parse(Int64, dtstr[35:41])
        #  2011     9     8     0     0    0.1330000     GPS         TIME OF FIRST OBS
        dt = DateTime(y, m, d, H, M, S, s ÷ 1e4)
        tststr = input[49:50]
        haskey(str2tst, tststr) || throw(error("no such time system:\n$input"))
        d = diff2UTC(dt, str2tst[tststr])
        dt + d,  d
    end
    parse_ohsLastTime(input::String) = returninput(input)
    parse_ohsReceiverOffset(input::String) = returninput(input)
    parse_ohsSystemDCBSapplied(input::String) = returninput(input)
    parse_ohsSystemPCVSapplied(input::String) = returninput(input)
    parse_ohsSystemScaleFac(input::String) = returninput(input)
    parse_ohsSystemPhaseShift(input::String) = returninput(input)
    parse_ohsGlonassSlotFreqNo(input::String) = returninput(input)
    parse_ohsGlonassCodPhsBias(input::String) = returninput(input)
    parse_ohsLeapSeconds(input::String) = returninput(input)
    parse_ohsNumSats(input::String) = returninput(input)
    parse_ohsPrnObs(input::String) = returninput(input)
    parse_ohsEoH(input::String) = returninput(input)
    function parse_ohsFile(intput::String)
        return basename(intput)
    end
end

const obs_hs_function = Dict{String,Function}(
    rchsVersion          => parse_rchsVersion,
    rchsProg             => parse_rchsProg,
    ohsVersion           => parse_ohsVersion,
    ohsRunBy             => parse_ohsRunBy,
    ohsComment           => parse_ohsComment,
    ohsMarkerName        => parse_ohsMarkerName,
    ohsMarkerNumber      => parse_ohsMarkerNumber,
    ohsMarkerType        => parse_ohsMarkerType,
    ohsObserver          => parse_ohsObserver,
    ohsReceiver          => parse_ohsReceiver,
    ohsAntennaType       => parse_ohsAntennaType,
    ohsAntennaPosition   => parse_ohsAntennaPosition,
    ohsAntennaDeltaHEN   => parse_ohsAntennaDeltaHEN,
    ohsAntennaDeltaXYZ   => parse_ohsAntennaDeltaXYZ,
    ohsAntennaPhaseCtr   => parse_ohsAntennaPhaseCtr,
    ohsAntennaBsightXYZ  => parse_ohsAntennaBsightXYZ,
    ohsAntennaZeroDirAzi => parse_ohsAntennaZeroDirAzi,
    ohsAntennaZeroDirXYZ => parse_ohsAntennaZeroDirXYZ,
    ohsCenterOfMass      => parse_ohsCenterOfMass,
    ohsNumObs            => parse_ohsNumObs,
    ohsSystemNumObs      => parse_ohsSystemNumObs3,
    ohsWaveFact          => parse_ohsWaveFact,
    ohsSigStrengthUnit   => parse_ohsSigStrengthUnit,
    ohsInterval          => parse_ohsInterval,
    ohsFirstTime         => parse_ohsFirstTime,
    ohsLastTime          => parse_ohsLastTime,
    ohsReceiverOffset    => parse_ohsReceiverOffset,
    ohsSystemDCBSapplied => parse_ohsSystemDCBSapplied,
    ohsSystemPCVSapplied => parse_ohsSystemPCVSapplied,
    ohsSystemScaleFac    => parse_ohsSystemScaleFac,
    ohsSystemPhaseShift  => parse_ohsSystemPhaseShift,
    ohsGlonassSlotFreqNo => parse_ohsGlonassSlotFreqNo,
    ohsGlonassCodPhsBias => parse_ohsGlonassCodPhsBias,
    ohsLeapSeconds       => parse_ohsLeapSeconds,
    ohsNumSats           => parse_ohsNumSats,
    ohsPrnObs            => parse_ohsPrnObs,
    ohsEoH               => parse_ohsEoH,
)

const obs_hstrings_set = Set(keys(obs_hs_function))

# obs header parse end

ext: jl


text:
# Not everything has an SONAME
get_soname(oh::ObjectHandle) = nothing

# Auto-open a path into an ObjectHandle
function get_soname(path::AbstractString)
    try
        readmeta(get_soname, path)
    catch e
        @warn "Could not probe $(path) for an SONAME!" exception=(e, catch_backtrace())
        return nothing
    end
end

function get_soname(oh::ELFHandle)
    # Get the dynamic entries, see if it contains a DT_SONAME
    es = ELFDynEntries(oh)
    soname_idx = findfirst(e -> e.entry.d_tag == ELF.DT_SONAME, es)
    if soname_idx === nothing
        # If all else fails, just return the filename.
        return nothing
    end

    # Look up the SONAME from the string table
    return strtab_lookup(es[soname_idx])
end

function get_soname(oh::MachOHandle)
    # Get the dynamic entries, see if it contains an ID_DYLIB_CMD
    lcs = MachOLoadCmds(oh)
    id_idx = findfirst(lc -> typeof(lc) <: MachOIdDylibCmd, lcs)
    if id_idx === nothing
        # If all else fails, just return the filename.
        return nothing
    end

    # Return the Dylib ID
    return dylib_name(lcs[id_idx])
end

function ensure_soname(prefix::Prefix, path::AbstractString, platform::Platform;
                       verbose::Bool = false, autofix::Bool = false)
    # Skip any kind of Windows platforms
    if platform isa Windows
        return true
    end

    # Skip if this file already contains an SONAME
    rel_path = relpath(realpath(path), realpath(prefix.path))
    soname = get_soname(path)
    if soname != nothing
        if verbose
            @info("$(rel_path) already has SONAME \"$(soname)\"")
        end
        return true
    else
        soname = basename(path)
    end

    # If we're not allowed to fix it, fail out
    if !autofix
        return false
    end

    # Otherwise, set the SONAME
    ur = preferred_runner()(prefix.path; cwd="/workspace/", platform=platform)
    set_soname_cmd = ``

    if Sys.isapple(platform)
        install_name_tool = "/opt/x86_64-apple-darwin14/bin/install_name_tool"
        set_soname_cmd = `$install_name_tool -id $(soname) $(rel_path)`
    elseif Sys.islinux(platform) || Sys.isbsd(platform)
        patchelf = "/usr/bin/patchelf"
        set_soname_cmd = `$patchelf $(patchelf_flags(platform)) --set-soname $(soname) $(rel_path)`
    end

    # Create a new linkage that looks like @rpath/$lib on OSX,
    retval = with_logfile(prefix, "set_soname_$(basename(rel_path))_$(soname).log") do io
        run(ur, set_soname_cmd, io; verbose=verbose)
    end

    if !retval
        @warn("Unable to set SONAME on $(rel_path)")
        return false
    end

    # Read the SONAME back in and ensure it's set properly
    new_soname = get_soname(path)
    if new_soname != soname
        @warn("Set SONAME on $(rel_path) to $(soname), but read back $(string(new_soname))!")
        return false
    end

    if verbose
        @info("Set SONAME of $(rel_path) to \"$(soname)\"")
    end

    return true
end

"""
    symlink_soname_lib(path::AbstractString)

We require that all shared libraries are accessible on disk through their
SONAME (if it exists).  While this is almost always true in practice, it
doesn't hurt to make doubly sure.
"""
function symlink_soname_lib(path::AbstractString; verbose::Bool = false,
                                                  autofix::Bool = false)
    # If this library doesn't have an SONAME, then just quit out immediately
    soname = get_soname(path)
    if soname === nothing
        return true
    end

    # Absolute path to where the SONAME-named file should be
    soname_path = joinpath(dirname(path), basename(soname))
    if !isfile(soname_path)
        if autofix
            if verbose
                @info("Library $(soname) does not exist, creating...")
            end
            symlink(path, soname_path)
        else
            if verbose
                @info("Library $(soname) does not exist, failing out...")
            end
            return false
        end
    end
    return true
end


text:
import RepoSnapshots

const GIT_DESTINATION_HOST = lowercase(strip(get(ENV, "GIT_DESTINATION_HOST", "",)))

if GIT_DESTINATION_HOST == "github"
    const GITHUB_SRC_ORGANIZATION = "UnofficialJuliaMirror"
    const GITHUB_DST_ORGANIZATION = "UnofficialJuliaMirrorSnapshots"
    const GITHUB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const GITHUB_BOT_PERSONAL_ACCESS_TOKEN =
        ENV["GITHUB_BOT_PERSONAL_ACCESS_TOKEN"]
    const src_provider =
        RepoSnapshots.Hosts.GitHubSrcHost.new_github_session(
            ;
            github_organization = GITHUB_SRC_ORGANIZATION,
            github_bot_username = GITHUB_BOT_USERNAME,
            github_bot_personal_access_token =
                GITHUB_BOT_PERSONAL_ACCESS_TOKEN,
            )
    const dst_provider =
        RepoSnapshots.Hosts.GitHubDstHost.new_github_session(
            ;
            github_organization = GITHUB_DST_ORGANIZATION,
            github_bot_username = GITHUB_BOT_USERNAME,
            github_bot_personal_access_token =
                GITHUB_BOT_PERSONAL_ACCESS_TOKEN,
            )
elseif GIT_DESTINATION_HOST == "gitlab"
    const GITLAB_SRC_GROUP = "UnofficialJuliaMirror"
    const GITLAB_DST_GROUP = "UnofficialJuliaMirrorSnapshots"
    const GITLAB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const GITLAB_BOT_PERSONAL_ACCESS_TOKEN =
        ENV["GITLAB_BOT_PERSONAL_ACCESS_TOKEN"]
    const src_provider =
        RepoSnapshots.Hosts.GitLabSrcHost.new_gitlab_session(
            ;
            gitlab_group = GITLAB_SRC_GROUP,
            gitlab_bot_username = GITLAB_BOT_USERNAME,
            gitlab_bot_personal_access_token =
                GITLAB_BOT_PERSONAL_ACCESS_TOKEN,
            )
    const dst_provider =
        RepoSnapshots.Hosts.GitLabDstHost.new_gitlab_session(
            ;
            gitlab_group = GITLAB_DST_GROUP,
            gitlab_bot_username = GITLAB_BOT_USERNAME,
            gitlab_bot_personal_access_token =
                GITLAB_BOT_PERSONAL_ACCESS_TOKEN,
            )
elseif GIT_DESTINATION_HOST == "bitbucket"
    const GITHUB_SRC_ORGANIZATION = "UnofficialJuliaMirror"
    const GITHUB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const GITHUB_BOT_PERSONAL_ACCESS_TOKEN =
        ENV["GITHUB_BOT_PERSONAL_ACCESS_TOKEN"]
    const BITBUCKET_DST_TEAM = "UnofficialJuliaMirrorSnapshots"
    const BITBUCKET_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const BITBUCKET_BOT_APP_PASSWORD =
        ENV["BITBUCKET_BOT_APP_PASSWORD"]
    const src_provider =
        RepoSnapshots.Hosts.GitHubSrcHost.new_github_session(
            ;
            github_organization = GITHUB_SRC_ORGANIZATION,
            github_bot_username = GITHUB_BOT_USERNAME,
            github_bot_personal_access_token =
                GITHUB_BOT_PERSONAL_ACCESS_TOKEN,
            )
    const dst_provider =
        RepoSnapshots.Hosts.BitbucketDstHost.new_bitbucket_session(
            ;
            bitbucket_team = BITBUCKET_DST_TEAM,
            bitbucket_bot_username = BITBUCKET_BOT_USERNAME,
            bitbucket_bot_app_password = BITBUCKET_BOT_APP_PASSWORD,
            )
else
    error(
        string(
            "You must set the \"GIT_DESTINATION_HOST\" environment variable, ",
            "and it must be set to one of the following values: ",
            "GITHUB GITLAB BITBUCKET",
            )
        )
end


text:
if ! haskey(ENV, "BULLET3")
  run(`./install-bullet3.sh`)
  println("Installed bullet3 at ", joinpath(ENV["HOME"], "bullet3"))
else
  println("Found bullet3 at $(ENV["BULLET3"])")
end


text:
using ..BasisDefinition, ..DynamicDefinition, ..Contractors, ..PwDynamicDefinition
using ValidatedNumerics, LinearAlgebra
import Base: iterate
import ..BasisDefinition: one_vector, integral_covector, is_integral_preserving, strong_norm, weak_norm, aux_norm

"""
Equispaced Ulam basis on [0,1] of size n
"""
struct Ulam <:Basis
	n::Integer #TODO, change to partition
end

Base.length(B::Ulam) = B.n

Base.iterate(B::Ulam, state = 1) = state < length(B)+1 ? (B[state-1], state+1) : nothing

"""
Returns the left endpoint of the i-th element of the Ulam basis
"""
Base.getindex(B::Ulam, i) = Float64(i)/B.n

"""
This iterator returns the preimages of the endpoints
of the intervals defining the Ulam basis through the dynamic
"""
function Base.iterate(S::DualComposedWithDynamic{Ulam, <:Dynamic}, state = (1, 1))
	i, k = state

	if i == length(S.basis)+1
			return nothing
	end

	# remark that this version supposes that for each i there exists a preimage
	# another more specific version should be implemented for maps with
	# incomplete branches

	x₁ = preim(S.dynamic, k, getindex(S.basis, i-1), S.ϵ)
	x₂ = preim(S.dynamic, k, getindex(S.basis, i), S.ϵ)

	lower, upper = x₁, x₂

	if k == nbranches(S.dynamic)
		return ((i, (lower, upper)), (i+1, 1))
	else
		return ((i, (lower, upper)), (i, k+1))
	end
end

function Base.iterate(S::DualComposedWithDynamic{Ulam, PwMap}, state = (1, 1))
	i, k = state

	if i == length(S.basis)+1
			return nothing
	end

	# remark that this version supposes that for each i there exists a preimage
	# another more specific version should be implemented for maps with
	# incomplete branches

	x₁ = preim(S.dynamic, k, getindex(S.basis, i-1), S.ϵ)
	#@info "x₁" x₁
	x₂ = preim(S.dynamic, k, getindex(S.basis, i), S.ϵ)
	#@info "x₂" x₂

	if S.dynamic.orientations[k]>0
			if isempty(x₁) && !isempty(x₁)
				x₁ = S.dynamic.endpoints[k]
			end
			if isempty(x₂) && !isempty(x₁)
				x₂ = S.dynamic.endpoints[k+1]
			end
			lower, upper = x₁, x₂
	elseif	S.dynamic.orientations[k]<0
			if isempty(x₂) && !isempty(x₁)
				x₂ = S.dynamic.endpoints[k]
			end
			if isempty(x₁) && !isempty(x₂)
				x₁ = S.dynamic.endpoints[k+1]
			end
			lower, upper = x₂, x₁
	end

	if k == nbranches(S.dynamic)
		if isempty(lower) && isempty(upper)
			return ((i, :∅), (i+1, 1))
		end
		return ((i, (lower, upper)), (i+1, 1))
	else
		if isempty(lower) && isempty(upper)
			return ((i, :∅), (i, k+1))
		end
		return ((i, (lower, upper)), (i, k+1))
	end
end

"""
Returns the indices of the elements of the Ulam basis that intersect with the interval y
"""
BasisDefinition.nonzero_on(B::Ulam, y) = max(floor(Int64, y[1].lo*length(B)), 1), min(ceil(Int64, y[2].hi*length(B)), length(B))

function relative_measure(S::Ulam, y, a, b)
	lower = max(y[1], a)
	upper = min(y[2], b)
	return max(upper-lower, 0)/(b-a)
end

"""
Given a preimage of an interval ```I_i```, this iterator returns
its relative intersection with all the elements of the Ulam basis that
have nonzero intersection with him
"""
function Base.iterate(S::ProjectDualElement{Ulam}, state = S.j_min)
	if state == S.j_max+1
		return nothing
	end

	return ((state, relative_measure(S.basis, S.dual_element,
			getindex(S.basis, state-1),
			getindex(S.basis, state))),
		    state+1)
end

BasisDefinition.evaluate(B::Ulam, i, x) = (x>(i-1)/n) && (x<i/n) ? 1 : 0

BasisDefinition.evaluate_integral(B::Ulam, i, T::Type)  = T(i)/length(B)

function Base.iterate(S::AverageZero{Ulam}, state = 1)
	n = length(S.basis)
	if state == n
		return nothing
	end
	v = zeros(Float64, n)
	v[1] = 1
	v[state+1]=-1
	return (v, state+1)
end

BasisDefinition.is_refinement(Bf::Ulam, Bc::Ulam) = rem(Bf.n, Bc.n) == 0

function integral_covector(B::Ulam)
	n = length(B)
	return 1/n * ones(Interval{Float64}, n)'
end

is_integral_preserving(B::Ulam) = true

function one_vector(B::Ulam)
	return ones(length(B))
end

BasisDefinition.strong_norm(B::Ulam) = TotalVariation
BasisDefinition.weak_norm(B::Ulam) = L1
BasisDefinition.aux_norm(B::Ulam) = L1

# See BasisDefinition for docs on these constants
# These must be rounded up correctly!

BasisDefinition.weak_projection_error(B::Ulam) = 0.5 ⊘₊ Float64(length(B), RoundDown)
BasisDefinition.aux_normalized_projection_error(B::Ulam) = 0.
BasisDefinition.strong_weak_bound(B::Ulam) = Float64(length(B), RoundUp)
BasisDefinition.aux_weak_bound(B::Ulam) = 1.
BasisDefinition.weak_by_strong_and_aux_bound(B::Ulam) = (0., 1.)
BasisDefinition.bound_weak_norm_from_linalg_norm(B::Ulam) = (1., 0.)
BasisDefinition.bound_linalg_norm_L1_from_weak(B::Ulam) = 1.
BasisDefinition.bound_linalg_norm_L∞_from_weak(B::Ulam) = Float64(length(B), RoundUp)

function BasisDefinition.invariant_measure_strong_norm_bound(B::Ulam, D::Dynamic)
	A, B = dfly(strong_norm(B), aux_norm(B), D)
	@assert A < 1.
	return B ⊘₊ (1. ⊖₋ A)
end

using RecipesBase
using LaTeXStrings

@userplot PlotUlam
@recipe function f(h::PlotUlam)
	if length(h.args)<2 || (typeof(h.args[1])!= Ulam) || !(typeof(h.args[2])<:AbstractVector)
		error("Plot Ulam needs as an input a Ulam Basis and a vector")
	end
	B = h.args[1]
	w = h.args[2]
	D = h.args[3]

	layout := (3, 1)
	link := :both

	if eltype(w) <: Interval
		w = mid.(w)
	end

	# bar plot
	@series begin
		linecolor := :blue
		linealpha := 0.8
		fillalpha := 0.8
		seriestype := :bar
		label := L"f_{\delta}"
		collect(B), w
	end

	#dynamic plot
	@series begin
		linecolor := :green
	    seriestype := :path
	    label := L"T(x)"
	    G = x-> plottable(D, x)
	    collect(B), G.(collect(B))
	end

	@series begin
		linecolor := :red
	    seriestype := :path
	    label := L"T'(x)"
	    G = x-> derivative(D, x)
	    collect(B), G.(collect(B))
	end
end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
"""
Randomly selects the genetic operator from the vector
according to its weight and applies it.
"""
struct FixedGeneticOperatorsMixture <: GeneticOperatorsMixture
    operators::Vector{GeneticOperator} # available operations
    weights::Weights{Float64, Float64, Vector{Float64}}  # fixed weights

    function FixedGeneticOperatorsMixture(
        operators::AbstractVector{<:GeneticOperator},
        rates::AbstractVector{Float64} = fill(1.0/length(operators), length(operators)) # defaults to uniform distribution of rates
    )
        length(operators) == length(rates) ||
            throw(DimensionMismatch("Number of mutators does not match the number of their rates"))
        new(GeneticOperator[op for op in operators], weights(rates))
    end
end

"""
Default implementation of `apply!()` for operators mixture.
"""
function apply!(opmix::GeneticOperatorsMixture,
                v::AbstractVector{<:Real}, target_index::Int)
    op, tag = next(opmix)
    apply!(op, v, target_index)
end

"""
    next(fa::FixedGeneticOperatorsMixture)

Gets the random genetic operator from the mixture.

The probability to select an operator is proportional to its weight.

Returns a tuple of the genetic operator and its index in the mix.
"""
function next(opmix::FixedGeneticOperatorsMixture)
    i = sample(1:length(opmix.operators), opmix.weights)
    return opmix.operators[i], i
end

"""
Frequency-adapting genetic operators mixture.
"""
struct FAGeneticOperatorsMixture <: GeneticOperatorsMixture
    operators::Vector{GeneticOperator} # available operators
    fa::FrequencyAdapter               # adapter of operator frequencies

    FAGeneticOperatorsMixture(
        operators::AbstractVector{<:GeneticOperator};
        c = 1E-2, eta = 1E-2, pmin = 0.01, pmax = 1.0
    ) = new(GeneticOperator[op for op in operators],
            FrequencyAdapter(length(operators);
                             c = c, eta = eta, pmin = pmin, pmax = pmax))
end

frequencies(opmix::FAGeneticOperatorsMixture) = frequencies(opmix.fa)

function next(opmix::FAGeneticOperatorsMixture)
    i = next(opmix.fa)
    return opmix.operators[i], i
end

function adjust!(opmix::FAGeneticOperatorsMixture, op_index::Int, candi_index::Int,
                 new_fitness::F, old_fitness::F, is_improved::Bool) where F
    # KLUDGE we don't know the fitness scheme, but if there is no improvement in fitness,
    # new_fitness==old_fitness, so it's ok to take the abs()
    update!(opmix.fa, op_index, is_improved ? abs(new_fitness-old_fitness) : 0.0)

    # also adjust the actual operator
    adjust!(opmix.operators[op_index], 0, candi_index, new_fitness, old_fitness, is_improved)
end

# FIXME use logging
function trace_state(io::IO, fa::FAGeneticOperatorsMixture, mode::Symbol)
    if mode == :verbose
        println(io, "FrequencyAdapting mixture rates:")
        for i in eachindex(fa.operators)
            @printf(io, " %s=%.2f", typeof(fa.operators[i]), fa.fa.p[i])
        end
    end
    println(io)
end


text:
export Dirichlet

"""
Description:

    Dirichlet factor node

    Multivariate:
    f(out, a) = Dir(out|a)
              = Γ(Σ_i a_i)/(Π_i Γ(a_i)) Π_i out_i^{a_i}
    where 'a' is a vector with every a_i > 0

    Matrix variate:
    f(out, a) = Π_k Dir(out|a_*k)
    where 'a' represents a left-stochastic matrix with every a_jk > 0

Interfaces:

    1. out
    2. a

Construction:

    Dirichlet(id=:some_id)
"""
mutable struct Dirichlet <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Dirichlet(out, a; id=generateId(Dirichlet))
        @ensureVariables(out, a)
        self = new(id, Array{Interface}(undef, 2), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:a] = self.interfaces[2] = associate!(Interface(self), a)

        return self
    end
end

slug(::Type{Dirichlet}) = "Dir"

format(dist::ProbabilityDistribution{V, Dirichlet}) where V<:VariateType = "$(slug(Dirichlet))(a=$(format(dist.params[:a])))"

ProbabilityDistribution(::Type{MatrixVariate}, ::Type{Dirichlet}; a=ones(3,3)) = ProbabilityDistribution{MatrixVariate, Dirichlet}(Dict(:a=>a))
ProbabilityDistribution(::Type{Multivariate}, ::Type{Dirichlet}; a=ones(3)) = ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>a))
ProbabilityDistribution(::Type{Dirichlet}; a=ones(3)) = ProbabilityDistribution{Multivariate, Dirichlet}(Dict(:a=>a))

dims(dist::ProbabilityDistribution{<:VariateType, Dirichlet}) = size(dist.params[:a])

vague(::Type{Dirichlet}, dims::Int64) = ProbabilityDistribution(Multivariate, Dirichlet, a=ones(dims))
vague(::Type{Dirichlet}, dims::Tuple{Int64}) = ProbabilityDistribution(Multivariate, Dirichlet, a=ones(dims))
vague(::Type{Dirichlet}, dims::Tuple{Int64, Int64}) = ProbabilityDistribution(MatrixVariate, Dirichlet, a=ones(dims))

isProper(dist::ProbabilityDistribution{V, Dirichlet}) where V<:VariateType = all(dist.params[:a] .> 0.0)

unsafeMean(dist::ProbabilityDistribution{Multivariate, Dirichlet}) = dist.params[:a]./sum(dist.params[:a])
unsafeMean(dist::ProbabilityDistribution{MatrixVariate, Dirichlet}) = dist.params[:a]./sum(dist.params[:a],dims=1) # Normalize columns

function unsafeLogMean(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    a = clamp.(dist.params[:a], tiny, huge)
    return digamma.(a) .- digamma.(sum(a))
end
function unsafeLogMean(dist::ProbabilityDistribution{MatrixVariate, Dirichlet})
    a = clamp.(dist.params[:a], tiny, huge)
    return digamma.(a) .- digamma.(sum(a,dims=1)) # Normalize columns
end

logPdf(dist::ProbabilityDistribution{Multivariate, Dirichlet}, x) = sum((dist.params[:a].-1).*log.(x)) - sum(loggamma.(dist.params[:a])) + loggamma(sum(dist.params[:a]))
logPdf(dist::ProbabilityDistribution{MatrixVariate, Dirichlet}, x) = sum(sum((dist.params[:a].-1).*log.(x),dims=1) - sum(loggamma.(dist.params[:a]), dims=1) + loggamma.(sum(dist.params[:a],dims=1)))

function unsafeVar(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    a_sum = sum(dist.params[:a])
    return dist.params[:a].*(a_sum .- dist.params[:a])./(a_sum^2*(a_sum + 1.0))
end

function prod!( x::ProbabilityDistribution{V, Dirichlet},
                y::ProbabilityDistribution{V, Dirichlet},
                z::ProbabilityDistribution{V, Dirichlet}=ProbabilityDistribution(V, Dirichlet, a=ones(dims(x)))) where V<:VariateType

    z.params[:a] = x.params[:a] + y.params[:a] .- 1.0

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{Multivariate, Dirichlet},
                            y::ProbabilityDistribution{Multivariate, PointMass},
                            z::ProbabilityDistribution{Multivariate, PointMass}=ProbabilityDistribution(Multivariate, PointMass, m=[NaN]))

    a_y = clamp.(y.params[:m], 0.0, 1.0)
    z.params[:m] = deepcopy(a_y)

    return z
end

@symmetrical function prod!(x::ProbabilityDistribution{MatrixVariate, Dirichlet},
                            y::ProbabilityDistribution{MatrixVariate, PointMass},
                            z::ProbabilityDistribution{MatrixVariate, PointMass}=ProbabilityDistribution(MatrixVariate, PointMass, m=mat(NaN)))

    A_y = clamp.(y.params[:m], 0.0, 1.0)
    z.params[:m] = deepcopy(A_y)

    return z
end

function sample(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    smpl = Vector{Float64}(undef, length(dist.params[:a]))
    for (i, alpha_i) in enumerate(dist.params[:a])
        smpl[i] = sample(ProbabilityDistribution(Univariate,Gamma,a=alpha_i,b=1.0))
    end
    smpl = smpl./sum(smpl)
    return smpl
end

naturalParams(dist::ProbabilityDistribution{Multivariate, Dirichlet}) = dist.params[:a] .- 1.0 # Variant 2 of https://en.wikipedia.org/wiki/Exponential_family
naturalParams(dist::ProbabilityDistribution{MatrixVariate, Dirichlet}) = vec(dist.params[:a]) .- 1.0

standardDistribution(::Type{Multivariate}, ::Type{Dirichlet}; η::Vector) = ProbabilityDistribution(Multivariate, Dirichlet, a=η.+1.0)
standardDistribution(::Type{MatrixVariate}, ::Type{Dirichlet}; η::Vector, dims::Tuple) = ProbabilityDistribution(MatrixVariate, Dirichlet, a=reshape(η, dims).+1.0) # Include dimensionality argument for rectangular case

logNormalizer(::Type{Multivariate}, ::Type{Dirichlet}; η::Vector) = sum(loggamma.(η.+1.0)) - loggamma(sum(η.+1.0))
logNormalizer(::Type{MatrixVariate}, ::Type{Dirichlet}; η::Vector, dims::Tuple) = sum(loggamma.(η.+1.0)) - sum(loggamma.(sum(reshape(η,dims).+1.0, dims=1)))

logPdf(V::Type{Multivariate}, F::Type{Dirichlet}, x::Vector; η::Vector) = log.(x)'*η - logNormalizer(V, F; η=η)
logPdf(V::Type{MatrixVariate}, F::Type{Dirichlet}, x::Matrix; η::Vector) = log.(vec(x))'*η - logNormalizer(V, F; η=η, dims=size(x))

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Multivariate, Dirichlet})
    a_sum = sum(dist.params[:a])

    -sum( (dist.params[:a] .- 1.0).*(digamma.(dist.params[:a]) .- digamma.(a_sum)) ) -
    labsgamma(a_sum) +
    sum(labsgamma.(dist.params[:a]))
end

function differentialEntropy(dist::ProbabilityDistribution{MatrixVariate, Dirichlet})
    H = 0.0
    for k = 1:dims(dist)[2] # For all columns
        a_sum = sum(dist.params[:a][:,k])

        H += -sum( (dist.params[:a][:,k] .- 1.0).*(digamma.(dist.params[:a][:,k]) .- digamma.(a_sum)) ) -
        labsgamma(a_sum) +
        sum(labsgamma.(dist.params[:a][:,k]))
    end

    return H
end

# Average energy functional
function averageEnergy(::Type{Dirichlet}, marg_out::ProbabilityDistribution{Multivariate}, marg_a::ProbabilityDistribution{Multivariate, PointMass})
    a_sum = sum(marg_a.params[:m])

    -labsgamma(a_sum) +
    sum(labsgamma.(marg_a.params[:m])) -
    sum( (marg_a.params[:m] .- 1.0).*unsafeLogMean(marg_out) )
end

function averageEnergy(::Type{Dirichlet}, marg_out::ProbabilityDistribution{Multivariate}, marg_a::ProbabilityDistribution{Multivariate, SampleList})
    samples, weights = marg_a.params[:s], marg_a.params[:w]
    S = length(weights) #number of samples
    log_gamma_of_sum, sum_of_log_gamma = 0.0, 0.0

    for s=1:S
        log_gamma_of_sum += weights[s]*loggamma(sum(samples[s]))
        sum_of_log_gamma += weights[s]*sum(loggamma.(samples[s]))
    end

    return -log_gamma_of_sum + sum_of_log_gamma -
            sum((unsafeMeanVector(marg_a) .- 1.0).*unsafeLogMean(marg_out))
end

function averageEnergy(::Type{Dirichlet}, marg_out::ProbabilityDistribution{MatrixVariate}, marg_a::ProbabilityDistribution{MatrixVariate, PointMass})
    (dims(marg_out) == dims(marg_a)) || error("Distribution dimensions must agree")

    log_mean_marg_out = unsafeLogMean(marg_out)

    H = 0.0
    for k = 1:dims(marg_out)[2] # For all columns
        a_sum = sum(marg_a.params[:m][:,k])

        H += -labsgamma(a_sum) +
        sum(labsgamma.(marg_a.params[:m][:,k])) -
        sum( (marg_a.params[:m][:,k] .- 1.0).*log_mean_marg_out[:,k] )
    end

    return H
end


text:
module Runtests_withPlot

using ModiaLang
using Test

const  test_title = "Test ModiaLang (version=$(ModiaLang.Version) with " * currentPlotPackage() * ")"

println("\n... $test_title")

@testset "$test_title" begin
    include("include_all.jl")
end

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
# fall back. An alternative definition should be provided
# if an array can be accessed via linear indexing in order (ie, 1,2,3..),
# and all elements are reached before length(A). Eg, if
# a triangular matrix is defined and does not store the off triangle.
# type_length(::NTuple{N}) where N= N
# type_length(::SArray{S,T,N,L}) where {S,T,N,L} = L
type_length(::Type{NTuple{N}}) where {N} = N
# type_length(::Type{SArray{S,T,N,L}}) where {S,T,N,L} = L
#type_length(::T) where {T} = type_length(T)
type_length(::T) where {T} = (@show T; throw("$T fallback type_length"))
type_length(::Type{A}) where {T,A<:AbstractArray{T}} = sizeof(A) ÷ sizeof(T)

# type_length(::Type{T}) where T = throw("type_length(::Type{$T}) has not been defined.")
# function type_length(::T) where T
#     @show T
#     throw("Stop!")
# end

"""
LinearStorage allows one to define a linear storing behavior different from standard linear indexing,
in case the two may differ. For example, consider a `L = LowerTrianglar{MMatrix{3,3,Float64,9}}`.
One is likely to want the standard linear indexing behavior for operations like `L .* A`.
However, for the construction of ScatteredArrays, one would prefer to only access the lower triangular elements.
Thus, ScatteredArrays accesses elements via LinearStorage, which default to the standard LinearIndex,
so that one may optionally provide other methods if one wants different behavior.
"""
struct LinearStorage end
@inline Base.getindex(A::AbstractArray, ::LinearStorage, i::Vararg{Any,N} where N) = getindex(A, i...)
@inline Base.setindex!(A::AbstractArray, v, ::LinearStorage, i::Vararg{Any,N} where N) = setindex!(A, v, i...)

"""
AbstractScatteredArray{E,M,A,N,Np1}
Parameters are:
E: Element type of the constituent arrays.
M: Dimensionality of the constituent arrays.
A: Type of the constituent arrays; element type of the scattered array.
N: Dimensionality of the scattered arrays.
Np1: N + 1; dimensionality of the stored data.
"""
abstract type AbstractScatteredArray{E,M,A<:AbstractArray{E,M},N,Np1} <: AbstractArray{A,N} end

@generated function Base.getindex(A::PaddedMatrices.AbstractFixedSizeArray{S,T,N,X,L}, ::LinearStorage, i::Vararg{<:Integer,NI}) where {S,T,N,X,L,NI}
    if N == 1 || ( ((S.parameters[1])::Int == (X.parameters[2])::Int) && (NI == 1) )
        @assert NI == 1
        return quote
            $(Expr(:meta,:inline))
            ind = i[1]
            @boundscheck ind > $L && PaddedMatrices.ThrowBoundsError()
            @inbounds A[ind]
        end
    elseif NI == 1
        return quote
            $(Expr(:meta,:inline))
            nrep, nrem = divrem(i[1], $(S.parameters[1]))
            ind = nrep * $P + nrem
            @boundscheck ind > $L && PaddedMatrices.ThrowBoundsError()
            @inbounds A[ind]
        end
    else
        @assert N == NI
        return quote
            $(Expr(:meta,:inline))
            $(Expr(:meta,:propagate_inbounds))
            A[i...]
        end
    end
end

## Lifted from metaprogramming documentation
# @generated function sub2ind_quote(dims::NTuple{N}, I::NTuple{N}) where N
#    ex = :(I[$N] - 1)
#    for i = (N - 1):-1:1
#        ex = :(I[$i] - 1 + dims[$i] * $ex)
#    end
#    :($ex + 1)
# end
@generated function sub2ind(dims::NTuple{N}, I::NTuple{N}) where N
   ex = :(I[$N] - 1)
   for i = (N - 1):-1:1
       ex = :(I[$i] - 1 + dims[$i] * $ex)
   end
   quote
       $(Expr(:meta,:inline))
       $ex + 1
   end
end

@noinline thrownotisbitserror(T) = throw("The type is required to be isbits, but $T is not.")
@noinline throwdimensionmismatcherror(dimA, dimB) = throw("Dimensions $dimA != $dimB.")


text:
module FooTests

using Test
import Foo

@test 2 == 2
Foo.greet()

end


text:
# ---
# layout: post
# title: "Pentachoron"
# date: 2019-04-03 00:00:00 +0000
# categories: blog
# mathjax: true
# ---

# > This page contains quite a few animated GIFs. I apologise in advance if your browser or network connection doesn't like them. It's much easier to add a GIF than it is to embed a video into a GitHub pages blog, unfortunately.

# ### Pentachoron

# ![image label](IMAGEFOLDER/header-image.svg){: .center-image}

# I encountered the word "pentachoron" recently for the first time, and I thought it had to be something cool. Turns out it is, so I thought I'd have a go at drawing one...

# Unfortunately, it's a 4-dimensional object, which is going to tax my limited 2D skills. But I like being out of my depth.

# ## First steps

# The first job is to define a couple of types to help organize all the various points.

# Here's an immutable structure to hold the coordinates of a point in four dimensions. In four dimensions, the axes are usually called `x`, `y`, `z`, and `w`.

struct Point4D <: AbstractArray{Float64, 1}
    x::Float64
    y::Float64
    z::Float64
    w::Float64
end

Point4D(a::Array{Float64, 1}) = Point4D(a...)

Base.size(pt::Point4D) = (4, )

Base.getindex(pt::Point4D, i) = [pt.x, pt.y, pt.z, pt.w][i]

# By adopting the AbstractArray as our stepmother type, we can pick up many useful behaviors. `size()` and `getindex()` can be taught how to handle Point4Ds. If you browse the [official documentation](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array-1) for a while you'll find out how to effectively tap in to many desirable pre-defined abilities that our new type can inherit.

# While we're here, let's do three-dimensional points as well:

struct Point3D <: AbstractArray{Float64, 1}
    x::Float64
    y::Float64
    z::Float64
end

Base.size(pt::Point3D) = (3, )

# ## Conversion

# The primary task we have to address is how to convert a 4D point into a 2D point. Let's start with the easier task: how to convert a 3D point into a 2D point, ie how can we draw a 3D shape on a flat surface?

# Consider a simple cube. The front face and the back face could have the same X and Y coordinates, and vary only by their Z values.

# (This is just a diagram, it's not really 3D...)

# ![image label](IMAGEFOLDER/cube.svg){: .center-image}

using Luxor

@svg begin
    fontface("Menlo")
    fontsize(8)
    setblend(blend(
        boxtopcenter(BoundingBox()),
        boxmiddlecenter(BoundingBox()),
        "skyblue",
        "white"))
    box(boxtopleft(BoundingBox()),
        boxmiddleright(BoundingBox()), :fill)

    setblend(blend(
        boxmiddlecenter(BoundingBox()),
        boxbottomcenter(BoundingBox()),
        "grey95",
        "grey45"
        ))
    box(boxmiddleleft(BoundingBox()),
        boxbottomright(BoundingBox()), :fill)

    sethue("black")
    setline(2)
    bx1 = box(O, 250, 250, vertices=true)
    poly(bx1, :stroke, close=true)
    label.(["-1/1/1", "-1/-1/1", "1/-1/1", "1/1/1"],
        slope.(O, bx1), bx1)

    setline(1)
    bx2 = box(O, 150, 150, vertices=true)
    poly(bx2, :stroke, close=true)
    label.(["-1/1/0", "-1/-1/0", "1/-1/0", "1/1/0"],
        slope.(O, bx2), bx2, offset=-45)

    map((x, y) -> line(x, y, :stroke), bx1, bx2)
end 500 500 "images/pentachoron/cube.svg"

# So the idea is to project the cube from 3D to 2D by keeping the first two values, and multiplying or modifying them by the third value. We'll make a `convert()` function to do this:

function convert(Point, pt3::Point3D)
    k = 1/(K - pt3.z)
    return Point(pt3.x * k, pt3.y * k)
end

# `K` is just a constant which provides a consistent value for depth:

const K = 4.0

# Testing this quickly in Luxor.jl gives promising results:

@svg begin
    cube = Point3D[
        Point3D(-1, -1,  1),
        Point3D(-1,  1,  1),
        Point3D( 1, -1,  1),
        Point3D( 1,  1,  1),
        Point3D(-1, -1, -1),
        Point3D(-1,  1, -1),
        Point3D( 1, -1, -1),
        Point3D( 1,  1, -1),
    ]
    circle.(convert.(Point, cube) * 300, 5, :fill)
end 500 500 "images/pentachoron/basic3_to_2.png"

# ![image label](IMAGEFOLDER/basic3_to_2.png){: .center-image}

# It's a simple type of perspective projection.

# Using the same principle, let's make a method for converting a 4D point.

function convert(Point3D, pt4::Point4D)
    k = 1/(K - pt4.w)
    return Point3D(pt4.x * k, pt4.y * k, pt4.z * k)
end

# Let's hope it works.

# We can combine these into a single utility function called `flatten()` that takes a list of 4D points and double-maps them into a list of 2D points suitable for drawing.

function flatten(shape4)
    return map(pt3 -> convert(Point, pt3), map(pt4 -> convert(Point3D, pt4), shape4))
end

# ## Going up a level

# To test this, we'll define the vertices of a unit pentachoron:

const n = -1/√5

const pentachoron = [Point4D(vertex...) for vertex in [
    [ 1.0,  1.0,  1.0, n],
    [ 1.0, -1.0, -1.0, n],
    [-1.0,  1.0, -1.0, n],
    [-1.0, -1.0,  1.0, n],
    [ 0.0,  0.0,  0.0, n + √5]]]

#

#-

# (According to [Wikipedia](https://en.wikipedia.org/wiki/5-cell), other names for the pentachoron include:

# * Regular 5-cell
# * C5
# * pentatope
# * pentahedroid
# * tetrahedral pyramid
# * 4-simplex

# I think *pentachoron* is the coolest, though.)

#-

#

# Here's a list of "faces", with each face defined by three of the vertices:

const pentachoronfaces = [
    [1, 2, 3],
    [1, 2, 4],
    [1, 2, 5],
    [1, 3, 4],
    [1, 3, 5],
    [1, 4, 5],
    [2, 3, 4],
    [2, 3, 5],
    [2, 4, 5],
    [3, 4, 5]]

# A quick test:

@svg begin
    setopacity(0.2)
    pentachoron2D = flatten(pentachoron)
    for (n, face) in enumerate(pentachoronfaces)
        randomhue()
        poly(1500 * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end 600 250 "images/pentachoron/firstdraw.svg"

# ![image label](IMAGEFOLDER/firstdraw.svg){: .center-image}

# This isn't very interesting, although I *think* it's correct. To see a more appealing display, let's make it dance...

# ## Swings and roundabouts

# You usually rotate 2D points about a 1D point. You usually rotate 3D points about a 2D line (often one of the XYZ axes). So logically you'd rotate 4D points with reference to a 3D plane. We'll define some matrices that do 4D rotations relative to a plane defined by two of the X, Y, Z, and W axes.

# As I'm thinking about it, the XY plane is typically the plane of the drawing surface. If you think of the XY plane as a computer screen straight in front of you, the XZ plane is parallel to your desk or floor, and the YZ plane is like a wall beside your desk on your right or left side.

# But what about the XW plane? And the YW and ZW planes, for that matter? This is the mystery of 4D shapes: we can't see these planes, we can only imagine their existence by watching shapes moving through and around them.

# Here's a function that generates a suitable rotation matrix for rotating a 4D point in the XY plane:

function XY(θ)
    [cos(θ)   -sin(θ)  0        0;
    sin(θ)    cos(θ)   0        0;
    0         0        1        0;
    0         0        0        1]
end

# and here's another one that rotates a 4D point in the XW plane:

function XW(θ)
    [cos(θ)   0        0        -sin(θ);
    0         1        0        0;
    0         0        1        0;
    sin(θ)    0        0        cos(θ)]
end

# Oh, it looks like we're having a matrix party, so we might as well do the others while we're here:

function XZ(θ)
    [cos(θ)   0        -sin(θ)  0;
    0         1        0        0;
    sin(θ)    0        cos(θ)   0;
    0         0        0        1]
end

function YZ(θ)
    [1        0        0        0;
    0         cos(θ)   -sin(θ)  0;
    0         sin(θ)   cos(θ)   0;
    0         0        0        1]
end

function YW(θ)
    [1        0        0        0;
    0         cos(θ)   0        -sin(θ);
    0         0        1        0;
    0         sin(θ)   0        cos(θ)]
end

function ZW(θ)
    [1        0        0        0;
    0         1        0        0;
    0         0        cos(θ)   -sin(θ);
    0         0        sin(θ)   cos(θ)];
end

# We'll make a handy utility function that uses one of these matrix functions to rotate an array of 4D points `A` in the plane defined by the matrix function:

function rotate4(A, matrixfunction)
    return map(A) do pt4
        Point4D(matrixfunction * pt4)
    end
end

# So we can write `rotate4(A, XW(π/2))` to rotate the array of points in `A` in the XW plane by π/2. We can also write `rotate4(A, XW(π/2) * XY(π/2))` to rotate the points first in the XW plane and then in the XY plane.

# Notice that the `rotate4()` function returns a new array, rather than modifying the original one. This probably has some disadvantages, but for now it means that we don't have to keep track of mutating objects so carefully.

# ## It moves

# We're now able to make an animation that rotates our pentachoron in one or more planes.

# The `frame()` function generates a single frame, using the framenumber (`eased_n` provides a normalized version of the framenumber modulated by the easing function) to control the rotation between 0 and 2π. The scalefactor makes the unit shape big enough to fill the frame.

# It's helpful to color the faces very slightly, using `poly(..., :fill)`-er.

using ColorSchemes

function frame(scene, framenumber, scalefactor=1000)
    background("antiquewhite")
    setlinejoin("bevel")
    setline(1.0)
    sethue("black")
    eased_n       = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    pentachoron′  = rotate4(pentachoron, XZ(eased_n * 2π))
    pentachoron2D = flatten(pentachoron′)

    setopacity(0.2)
    for (n, face) in enumerate(pentachoronfaces)
        sethue(get(ColorSchemes.diverging_rainbow_bgymr_45_85_c67_n256,
            n/length(pentachoronfaces)))
        poly(scalefactor * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

# And this function generates a series of frames and saves them to an animation:

function makemovie(w, h, fname;
        scalefactor=1000)
    movie1 = Movie(w, h, "4D movie")
    animate(movie1,
        Scene(movie1, (s, f)  -> frame(s, f, scalefactor),
            1:300,
            easingfunction=easeinoutsine),
        creategif=true,
        pathname="images/pentachoron/$(fname)")
end

makemovie(600, 600, "pentachoron-xz.gif", scalefactor=2000)

# This is a rotation in the "parallel to the desk my computer's resting on" XZ plane.

# ![image label](IMAGEFOLDER/pentachoron-xz.gif){: .center-image}

# It's not a realistic rendering with hidden-surface removal—I wouldn't know how to do that in 3D, let alone 4D. But it does make it slightly easier to follow some of the faces. Something to do one day might be to sort the "faces" to determine which are drawn first, but given the maths involved this is a challenge I'm dubious about tackling...

# We can modify `frame()` to allow various permutations of XY, XW, XZ, YZ, YW, and ZW rotations. Here's a version that rotates in XZ, followed by YW:

function frame(scene, framenumber, scalefactor=1000)
    background("antiquewhite")
    setlinejoin("bevel")
    setline(1.0)
    setopacity(0.2)

    eased_n       = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    pentachoron2D = flatten(
        rotate4(
            pentachoron,
                XZ(eased_n * 2π) *
                YW(eased_n * 2π)))
    for (n, face) in enumerate(pentachoronfaces)
        sethue(get(ColorSchemes.diverging_rainbow_bgymr_45_85_c67_n256,
            n/length(pentachoronfaces)))
        poly(scalefactor * pentachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

makemovie(800, 800, "pentachoron-xz-yw.gif", scalefactor=2000)

# ![image label](IMAGEFOLDER/pentachoron-xz-yw.gif){: .center-image}

# These are starting to make me think of kneading machines for making pizza dough.

# You can draw a few different models at once, placing each one in a separate cell:

# ![image label](IMAGEFOLDER/showall-pentachoronmovie.gif){: .center-image}

# This is showing 2D shadows of 3D shadows of 4D objects. Do these objects exist? I don't know...

# ### Other shapes

# While the word "pentachoron" was unfamiliar to me until recently, the word "tesseract" is much more familiar. After all, didn't the tesseract contain the Space Stone, one of the six Infinity Stones that apparently predate the universe and possess unlimited energy? (According to the scientists at Marvel Comics, at least...) It's also the name of a 4-dimensional hypercube, one step up from the pentachoron.

# Here are the vital numbers for the tesseract:

const tesseract = [Point4D(vertex...) for vertex in [
    [-1, -1, -1,  1],
    [ 1, -1, -1,  1],
    [ 1,  1, -1,  1],
    [-1,  1, -1,  1],
    [-1, -1,  1,  1],
    [ 1, -1,  1,  1],
    [ 1,  1,  1,  1],
    [-1,  1,  1,  1],
    [-1, -1, -1, -1],
    [ 1, -1, -1, -1],
    [ 1,  1, -1, -1],
    [-1,  1, -1, -1],
    [-1, -1,  1, -1],
    [ 1, -1,  1, -1],
    [ 1,  1,  1, -1],
    [-1,  1,  1, -1]]]

const tesseractfaces = [
    [1,   2,   3,   4],
    [1,   2,  10,   9],
    [1,   4,   8,   5],
    [1,   5,   6,   2],
    [1,   9,  12,   4],
    [2,   3,  11,  10],
    [2,   3,   7,   6],
    [3,   4,   8,   7],
    [5,   6,  14,  13],
    [5,   6,   7,   8],
    [5,   8,  16,  13],
    [6,   7,  15,  14],
    [7,   8,  16,  15],
    [9,  10,  11,  12],
    [9,  10,  14,  13],
    [9,  13,  16,  12],
    [10, 11,  15,  14],
    [13, 14,  15,  16]]

#-
# Replace `pentachoron` everywhere by `tesseract` in the above code, choosing the rotations to taste.
#-

function frame(scene, framenumber, scalefactor=1000)
    background("black")
    setlinejoin("bevel")
    setline(10.0)
    setopacity(0.7)

    eased_n       = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    tesseract2D = flatten(
        rotate4(
            tesseract,
                XZ(eased_n * 2π) *
                YW(eased_n * 2π)))
    for (n, face) in enumerate(tesseractfaces)
        sethue([Luxor.lighter_blue, Luxor.lighter_green,
            Luxor.lighter_purple, Luxor.lighter_red][mod1(n, 4)]...)
        poly(scalefactor * tesseract2D[face], :fillpreserve, close=true)
        sethue([Luxor.darker_blue, Luxor.darker_green,
            Luxor.darker_purple, Luxor.darker_red][mod1(n, 4)]...)
        strokepath()
    end
end

makemovie(800, 800, "tesseract-xz-yw.gif", scalefactor=2000)

# ![image label](IMAGEFOLDER/tesseract-xz-yw.gif){: .center-image}

# (Just for fun I temporarily switched the color scheme over to the Julia logo colors. You can see the 3D rendering problems more clearly, as well!)

# Here's a composite of four pairs of rotations:

# ![image label](IMAGEFOLDER/showall-tesseractmovie.gif){: .center-image}

# ## Enough is enough

# It's tempting to go further, but these web pages do get very big with all these GIFs. But just one more... I quite like this one, the *hexadecachoron*:

const hexadecachoron = [Point4D(vertex...) for vertex in [
        [ 1,  0,  0,  0],
        [ 0,  1,  0,  0],
        [ 0,  0,  1,  0],
        [ 0,  0,  0,  1],
        [-1,  0,  0,  0],
        [ 0, -1,  0,  0],
        [ 0,  0, -1,  0],
        [ 0,  0,  0, -1]]]

using Combinatorics

const hexadecachoronfaces = combinations(1:8, 3)

function frame(scene, framenumber, scalefactor=1000)
    background("antiquewhite")
    setlinejoin("bevel")
    setline(1.0)
    setopacity(0.15)
    eased_n          =
        scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    hexadecachoron2D = flatten(
        rotate4(
            hexadecachoron,
                YZ(eased_n * 2π) *
                XZ(eased_n * 2π) *
                YW(eased_n * 2π)))
    for (n, face) in enumerate(hexadecachoronfaces)
        sethue(get(ColorSchemes.magma,
            n/length(hexadecachoronfaces)))
        poly(scalefactor * hexadecachoron2D[face], :fillpreserve, close=true)
        sethue("black")
        strokepath()
    end
end

makemovie(650, 650, "movie-hexadecachoron.gif", scalefactor=5000)

# ![image label](IMAGEFOLDER/movie-hexadecachoron.gif){: .center-image}

# There's lots more fun to be had (combining two or more different shapes is fun), but I'm worried about your network and my brain cells, so that's enough animated GIFs for today.

# [2018-04-03]

# ![cormullion signing off](http://steampiano.net/cormullionknot.gif?pentachoron){: .center-image}

                                                                                #src
using Literate                                                                           #src
# preprocess for notebooks                                                               #src
function setimagefolder(content)                                                         #src
    content = replace(content, "IMAGEFOLDER" => "$IMAGEFOLDER")                          #src
    return content                                                                       #src
end                                                                                #src
# for Jupyter notebook, put images in subfolder                                          #src
#IMAGEFOLDER = "images/pentachoron"                                                       #src
#Literate.notebook("source/pentachoron.jl", "notebooks", preprocess = setimagefolder)     #src
# for Markdown/Jekyll notebook, put images in "/images"                                  #src
IMAGEFOLDER = "/images/pentachoron"                                                      #src
Literate.markdown("source/pentachoron.jl", ".", name="_posts/2019-04-03-pentachoron",    #src
 preprocess = setimagefolder,                                                            #src
 codefence = "{% highlight julia %}" => "{% endhighlight julia %}",                      #src
 documenter=false)                                                                       #src
                                                                                #src


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ResourceProviderOperationDisplayProperties <: SwaggerModel
    publisher::Any # spec type: Union{ Nothing, String } # spec name: publisher
    provider::Any # spec type: Union{ Nothing, String } # spec name: provider
    resource::Any # spec type: Union{ Nothing, String } # spec name: resource
    operation::Any # spec type: Union{ Nothing, String } # spec name: operation
    description::Any # spec type: Union{ Nothing, String } # spec name: description

    function ResourceProviderOperationDisplayProperties(;publisher=nothing, provider=nothing, resource=nothing, operation=nothing, description=nothing)
        o = new()
        validate_property(ResourceProviderOperationDisplayProperties, Symbol("publisher"), publisher)
        setfield!(o, Symbol("publisher"), publisher)
        validate_property(ResourceProviderOperationDisplayProperties, Symbol("provider"), provider)
        setfield!(o, Symbol("provider"), provider)
        validate_property(ResourceProviderOperationDisplayProperties, Symbol("resource"), resource)
        setfield!(o, Symbol("resource"), resource)
        validate_property(ResourceProviderOperationDisplayProperties, Symbol("operation"), operation)
        setfield!(o, Symbol("operation"), operation)
        validate_property(ResourceProviderOperationDisplayProperties, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        o
    end
end # type ResourceProviderOperationDisplayProperties

const _property_map_ResourceProviderOperationDisplayProperties = Dict{Symbol,Symbol}(Symbol("publisher")=>Symbol("publisher"), Symbol("provider")=>Symbol("provider"), Symbol("resource")=>Symbol("resource"), Symbol("operation")=>Symbol("operation"), Symbol("description")=>Symbol("description"))
const _property_types_ResourceProviderOperationDisplayProperties = Dict{Symbol,String}(Symbol("publisher")=>"String", Symbol("provider")=>"String", Symbol("resource")=>"String", Symbol("operation")=>"String", Symbol("description")=>"String")
Base.propertynames(::Type{ ResourceProviderOperationDisplayProperties }) = collect(keys(_property_map_ResourceProviderOperationDisplayProperties))
Swagger.property_type(::Type{ ResourceProviderOperationDisplayProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ResourceProviderOperationDisplayProperties[name]))}
Swagger.field_name(::Type{ ResourceProviderOperationDisplayProperties }, property_name::Symbol) =  _property_map_ResourceProviderOperationDisplayProperties[property_name]

function check_required(o::ResourceProviderOperationDisplayProperties)
    true
end

function validate_property(::Type{ ResourceProviderOperationDisplayProperties }, name::Symbol, val)
end


text:
function handle_step(obj::MultiSatStepHandler, arg0::List, arg1::jboolean)
    return jcall(obj, "handleStep", void, (List, jboolean), arg0, arg1)
end

function init(obj::MultiSatStepHandler, arg0::List, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (List, AbsoluteDate), arg0, arg1)
end


text:
using Sundials, Test

println("Test error handling")

f_error(u, p, t) = u / t
u0 = 1.0
prob = ODEProblem(f_error, u0, (0.0, 1.0))
sol = solve(prob, CVODE_BDF())
sol = solve(prob, CVODE_BDF(), verbose = false)

f_error2(du, u, p, t) = u / t - 1
u0 = 1.0;
du0 = 1.0;
prob = DAEProblem(f_error2, u0, du0, (0.0, 1.0), differential_vars = [1])
sol = solve(prob, IDA())
sol = solve(prob, IDA(), verbose = false)

@test sol.retcode == :InitialFailure


text:
"""
    streamgraph(x, y, group)

Creates an `EChart` where region below plotted line filled with color.

## Methods
```julia
streamgraph(x::AbstractVector{String},y::AbstractVector{<:Union{Missing, Real}},
            group::AbstractVector)
streamgraph(x::AbstractVector{<:Union{Missing, Real}},
            y::AbstractVector{<:Union{Missing, Real}},
            group::AbstractVector)
streamgraph(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol)
streamgraph(x::AbstractVector{<:Dates.TimeType},y::AbstractVector{<:Union{Missing, Real}},             group::AbstractVector)
```

## Arguments
* `legend::Bool = false` : display legend?
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

## Examples
```julia
s_df = readtable(Pkg.dir("ECharts", "exampledata/streamdata.csv"))
sg = streamgraph(s_df[:date], s_df[:value], s_df[:key], legend = true)
```
"""
function streamgraph(x::AbstractVector{String},
                    y::AbstractVector{<:Union{Missing, Real}},
                    group::AbstractVector;
          legend::Bool = false,
          kwargs...)

    ec = ECharts.newplot(kwargs, ec_charttype = "streamgraph")

    data = ECharts.arrayofarray(x,y,group)

    #Character value would represent dates
    ec.singleAxis = ECharts.SingleAxis(_type = "time")

    ec.series = [ThemeRiverSeries(data = data)]

    #need to modify legend!
    legend ? legend!(ec) : nothing

    return ec

end

function streamgraph(x::AbstractVector{<:Union{Missing, Real}},
                     y::AbstractVector{<:Union{Missing, Real}},
                     group::AbstractVector;
          legend::Bool = false,
          kwargs...)

    ec = ECharts.newplot(kwargs, ec_charttype = "streamgraph")

    data = ECharts.arrayofarray(x,y,group)

    #Character value would represent dates
    ec.singleAxis = ECharts.SingleAxis()

    ec.series = [ThemeRiverSeries(data = data)]

    #need to modify legend!
    legend ? legend!(ec) : nothing

    return ec

end

streamgraph(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol;
            legend::Bool = false,
            kwargs...) = streamgraph(df[x], df[y], df[group]; legend = legend, kwargs...)

streamgraph(x::AbstractVector{<:Dates.TimeType},
            y::AbstractVector{<:Union{Missing, Real}},
            group::AbstractVector;
            legend::Bool = false,
            kwargs...) = streamgraph(x, y, group; legend = legend, kwargs...)


text:
export BlockScan

mutable struct BlockScan
    file::String
    startbyte::Int
    endbyte::Int
    summary::Dict{String, Array{Int32, 1}}
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

for f in (:exp, :expm1)
    @eval begin
        ($f)(a::Interval{T}) where T = @round( ($f)(a.lo), ($f)(a.hi) )
        ($f)(x::EmptyInterval) = x
    end
end

ln_interval(::Type{T}, ::Val{N}) where {T, N} = log(unsafe_atomic(T, N))
ln2_interval(::Type{T}) where T = ln_interval(T, Val(2))
ln10_interval(::Type{T}) where T = ln_interval(T, Val(10))

# no CRlibm version
exp2(x::Interval{T}) where T = exp(ln2_interval(T) * x)
exp10(x::Interval{T}) where T = exp(ln10_interval(T) * x)

for f in (:exp2, :exp10)
  # add BigFloat functions with rounding:
    @eval ($f)(a::Interval{BigFloat}) = @round( ($f)(a.lo), ($f)(a.hi) )
    @eval ($f)(x::EmptyInterval) = x
end

for f in (:log, :log2, :log10)
    @eval ($(Symbol(f, "_inner")))(x::EmptyInterval) = x
    @eval ($(Symbol(f, "_inner")))(a::Interval) =
        a == zero(a) ? emptyinterval(a) : @round(($f)(a.lo), ($f)(a.hi))
    @eval function ($f)(a::Interval{T}) where T
            domain = nonnegatives(T)
            return ($(Symbol(f, "_inner")))(a ∩ domain)
        end
    @eval ($f)(x::EmptyInterval) = x
end

log1p_inner(a::Interval{T}) where T =
    a == -one(a) ? emptyinterval(T) : @round( log1p(a.lo), log1p(a.hi) )
log1p_inner(x::EmptyInterval) = x

function log1p(a::Interval{T}) where T
    domain = unsafe_interval(-one(T), T(Inf))
    return log1p_inner(a ∩ domain)
end
log1p(x::EmptyInterval) = x


text:
module MyExample

using ForwardDiff

include("extra_file.jl")

#= You should export only functions for public use.
	Internal function should NEVER be export but imported (import, not using)
=#
export my_f, derivative_of_my_f

end # module