text:
"""
This struct contains the reserve up market parameters
for price forcasts using CEM and endogeneous Economic Dispatch.
The Reserve Up Market is implemented using a piece-wise linear Operating Reserve Demand Curve (ORDC),
with break-points and price-points parameterized within this struct.
"""

struct ReserveUpMarket

    break_points::AxisArrays.AxisArray{Vector{Float64}}   # ORDC break points
    price_points::AxisArrays.AxisArray{Vector{Float64}}    # $/MWh
    demand::AxisArrays.AxisArray{Float64, 2}

    function ReserveUpMarket(b, p, d)
        @assert length(p) == length(b)
        @assert all(d .>= 0)
        for i = 1:length(b)
            @assert all(b[i] .>= 0)
            @assert all(p[i] .>= 0)

            @assert length(b[i]) > 1
            @assert length(b[i]) == length(p[i])
        end
        new(b, p, d)
    end

end


text:
import Test: @testset, @test

include("src.jl")

@testset "tvl_or truth table" begin
    @test @_tvl_or(true, true) === true
    @test @_tvl_or(true, false) === true
    @test @_tvl_or(true, missing) === true
    @test @_tvl_or(false, true) === true
    @test @_tvl_or(false, false) === false
    @test @_tvl_or(false, missing) === missing
    @test @_tvl_or(missing, true) === true
    @test @_tvl_or(missing, false) === missing
    @test @_tvl_or(missing, missing) === missing
end

@testset "tvl_and truth table" begin
    @test @_tvl_and(true, true) === true
    @test @_tvl_and(true, false) === false
    @test @_tvl_and(true, missing) === missing
    @test @_tvl_and(false, true) === false
    @test @_tvl_and(false, false) === false
    @test @_tvl_and(false, missing) === false
    @test @_tvl_and(missing, true) === missing
    @test @_tvl_and(missing, false) === false
    @test @_tvl_and(missing, missing) === missing
end

# We define a function that prints out a unique ID for each argument
# to a Boolean operator. By wrapping all Boolean values in calls to this
# function, we're able to check that the order of evaluation and
# side-effects of the short-circuiting operators are retained by our
# macro rewrites.

function f(io, i, x)
    print(io, i)
    x
end

@testset "Order of evaluation for tvl" begin
    for x in (true, false)
        for y in (true, false)
            for z in (true, false)
                io = IOBuffer()

                a = f(io, 1, x) && f(io, 2, y) || f(io, 3, z)
                order_a = String(take!(io))
                b = @tvl f(io, 1, x) && f(io, 2, y) || f(io, 3, z)
                order_b = String(take!(io))

                @test a === b
                @test order_a === order_b

                a = f(io, 1, x) || f(io, 2, y) && f(io, 3, z)
                order_a = String(take!(io))
                b = @tvl f(io, 1, x) || f(io, 2, y) && f(io, 3, z)
                order_b = String(take!(io))

                @test a === b
                @test order_a === order_b
            end
        end
    end
end

ext: jl


text:
# ----------------------------------------------------
# --- Run time UHD enumeration
# ----------------------------------------------------
# Direclry inherited from C file tune_request.h
@enum uhd_tune_request_policy_t begin
	UHD_TUNE_REQUEST_POLICY_NONE=78;
	UHD_TUNE_REQUEST_POLICY_AUTO=65;
	UHD_TUNE_REQUEST_POLICY_MANUAL=77
end
# Directly inherited from C file usrp.h
@enum uhd_stream_mode_t  begin
	UHD_STREAM_MODE_START_CONTINUOUS   = 97;
	UHD_STREAM_MODE_STOP_CONTINUOUS    = 111;
	UHD_STREAM_MODE_NUM_SAMPS_AND_DONE = 100;
	UHD_STREAM_MODE_NUM_SAMPS_AND_MORE = 109
end
@enum error_code_t begin
	ERROR_CODE_NONE = 0x0;
	ERROR_CODE_TIMEOUT = 0x1;
	ERROR_CODE_LATE_COMMAND = 0x2;
	ERROR_CODE_BROKEN_CHAIN = 0x4;
	ERROR_CODE_OVERFLOW = 0x8;
	ERROR_CODE_ALIGNMENT = 0xc;
	ERROR_CODE_BAD_PACKET = 0xf;
	BIG_PROBLEM;
end
@enum uhd_error begin
	UHD_ERROR_NONE = 0;
	UHD_ERROR_INVALID_DEVICE = 1;
	UHD_ERROR_INDEX = 10;
	UHD_ERROR_KEY = 11;
	UHD_ERROR_NOT_IMPLEMENTED = 20;
	UHD_ERROR_USB = 21;
	UHD_ERROR_IO = 30;
	UHD_ERROR_OS = 31;
	UHD_ERROR_ASSERTION = 40;
	UHD_ERROR_LOOKUP = 41;
	UHD_ERROR_TYPE = 42;
	UHD_ERROR_VALUE = 43;
	UHD_ERROR_RUNTIME = 44;
	UHD_ERROR_ENVIRONMENT = 45;
	UHD_ERROR_SYSTEM = 46;
	UHD_ERROR_EXCEPT = 47;
	UHD_ERROR_BOOSTEXCEPT = 60;
	UHD_ERROR_STDEXCEPT = 70;
	UHD_ERROR_UNKNOWN = 100
end

# ----------------------------------------------------
# --- Runtime structures
# ----------------------------------------------------
# These structures are necessary to run the wrapper
# They are used for parameter identification, Rx Tx stage config
# and to configure the LO
struct uhd_stream_args_t
	cpu_format::Cstring
	otw_format::Cstring;
	args::Cstring;
	channel_list::Ref{Csize_t};
	n_channels::Cint;
end
struct uhd_tune_request_t
	target_freq::Cdouble;
	rf_freq_policy::uhd_tune_request_policy_t;
	rf_freq::Cdouble;
	dsp_freq_policy::uhd_tune_request_policy_t;
	dsp_freq::Cdouble;
	args::Cstring;
end
struct uhd_tune_result
	clipped_rf_freq::Cdouble;
	target_rf_freq::Cdouble;
	actual_rf_freq::Cdouble;
	target_dsp_freq::Cdouble;
	actual_dsp_freq::Cdouble;
end
struct stream_cmd
	stream_mode::uhd_stream_mode_t;
	num_samps::Csize_t;
	stream_now::Cint;
	time_spec_full_secs::Cintmax_t;
	time_spec_frac_secs::Cdouble;
end
# --- Structure to manipulate UHD time stamps
struct Timestamp
	intPart::FORMAT_LONG;
	fracPart::Cdouble;
end

"""
" @assert_uhd macro
# Get the current UHD flag and raise an error if necessary
"""
macro assert_uhd(ex)
	quote
		local flag = $(esc(ex));
		if flag == UHD_ERROR_KEY
			# --- Specific
			error("Unable to create the UHD device. No attached UHD device found.");
		elseif flag != UHD_ERROR_NONE
			error("Unable to create or instantiate the UHD device. The return error flag is $flag");
		end
	end
end

# --- Global UHD structure from UHD
mutable struct uhd_usrp
end

# ----------------------------------------------------
# --- Rx structures
# ----------------------------------------------------
# --- C UHD streamer
mutable struct uhd_rx_streamer
end
# --- Metadata structure
struct uhd_rx_metadata
	has_time_spec::Cuchar;
	time_spec::Clonglong;
	time_spec_frac::Cdouble;
	more_fragments::Cuchar;
	fragment_offset::Csize_t;
	start_of_burst::Cuchar;
	end_of_burst::Cuchar;
	eov_positions::Ref{Csize_t};
	eov_positions_size::Csize_t;
	eov_positions_count::Csize_t;
	error_code::error_code_t;
	out_of_sequence::Cuchar;
end
# --- Structure with pointer reference
struct UHDRxWrapper
	flag::Bool;
	pointerUSRP::Ptr{uhd_usrp};
	pointerStreamer::Ptr{uhd_rx_streamer};
	pointerMD::Ptr{uhd_rx_metadata};
	addressStream::Ref{Ptr{uhd_rx_streamer}};
	addressMD::Ref{Ptr{uhd_rx_metadata}};
	pointerSamples::Ref{Csize_t}
end
# --- Main Rx structure
mutable struct UHDRx
	uhd::UHDRxWrapper;
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Int,Float64};
	antenna::String;
	packetSize::Csize_t;
	released::Int;
end

# ----------------------------------------------------
# --- Tx structures
# ----------------------------------------------------
# --- C UHD streamer
mutable struct uhd_tx_streamer
end
# --- Metadata structure
struct uhd_tx_metadata
	has_time_spec::Cuchar;
	time_spec::Clonglong;
	time_spec_frac::Cdouble;
	start_of_burst::Cuchar;
	end_of_burst::Cuchar;
	eov_positions::Ref{Csize_t};
	eov_positions_size::Csize_t;
end

# --- Structure with pointer reference
struct UHDTxWrapper
	flag::Bool;
	pointerUSRP::Ptr{uhd_usrp};
	pointerStreamer::Ptr{uhd_tx_streamer};
	pointerMD::Ptr{uhd_tx_metadata};
	addressStream::Ref{Ptr{uhd_tx_streamer}};
	addressMD::Ref{Ptr{uhd_tx_metadata}};
end
# --- Main Tx structure
mutable struct UHDTx
	uhd::UHDTxWrapper;
	carrierFreq::Float64;
	samplingRate::Float64;
	gain::Union{Int,Float64};
	antenna::String;
	packetSize::Csize_t;
	released::Int;
end

# ----------------------------------------------------
# --- UHDBinding structure
# ----------------------------------------------------
mutable struct UHDBinding
	addressUSRP::Ref{Ptr{uhd_usrp}};
	rx::UHDRx;
	tx::UHDTx;
end

# ----------------------------------------------------
# --- Print function
# ----------------------------------------------------
# To print fancy message with different colors with Tx and Rx
function customPrint(str,handler;style...)
    msglines = split(chomp(str), '\n')
    printstyled("┌",handler,": ";style...)
    println(msglines[1])
    for i in 2:length(msglines)
        (i == length(msglines)) ? symb="└ " : symb = "|";
        printstyled(symb;style...);
        println(msglines[i]);
    end
end
# define macro for printing Rx info
macro inforx(str)
    quote
        customPrint($(esc(str)),"Rx";bold=true,color=:light_green)
    end
end
# define macro for printing Rx warning
macro warnrx(str)
    quote
        customPrint($(esc(str)),"Rx Warning";bold=true,color=:light_yellow)
    end
end
# define macro for printing Tx info
macro infotx(str)
    quote
        customPrint($(esc(str)),"Tx";bold=true,color=:light_blue)
    end
end
# define macro for printing Tx warning
macro warntx(str)
    quote
        customPrint($(esc(str)),"Tx Warning";bold=true,color=:light_yellow)
    end
end


text:
#### Domain && Support

struct RealInterval
    lb::Float64
    ub::Float64

    RealInterval(lb::Real, ub::Real) = new(Float64(lb), Float64(ub))
end

minimum(r::RealInterval) = r.lb
maximum(r::RealInterval) = r.ub
extrema(r::RealInterval) = (r.lb, r.ub)
in(x::Real, r::RealInterval) = (r.lb <= Float64(x) <= r.ub)

isbounded(d::Union{D,Type{D}}) where {D<:UnivariateDistribution} = isupperbounded(d) && islowerbounded(d)

islowerbounded(d::Union{D,Type{D}}) where {D<:UnivariateDistribution} = minimum(d) > -Inf
isupperbounded(d::Union{D,Type{D}}) where {D<:UnivariateDistribution} = maximum(d) < +Inf

hasfinitesupport(d::Union{D,Type{D}}) where {D<:DiscreteUnivariateDistribution} = isbounded(d)
hasfinitesupport(d::Union{D,Type{D}}) where {D<:ContinuousUnivariateDistribution} = false

"""
    params(d::UnivariateDistribution)

Return a tuple of parameters. Let `d` be a distribution of type `D`, then `D(params(d)...)`
will construct exactly the same distribution as ``d``.
"""
params(d::UnivariateDistribution)

"""
    scale(d::UnivariateDistribution)

Get the scale parameter.
"""
scale(d::UnivariateDistribution)

"""
    location(d::UnivariateDistribution)

Get the location parameter.
"""
location(d::UnivariateDistribution)

"""
    shape(d::UnivariateDistribution)

Get the shape parameter.
"""
shape(d::UnivariateDistribution)

"""
    rate(d::UnivariateDistribution)

Get the rate parameter.
"""
rate(d::UnivariateDistribution)

"""
    ncategories(d::UnivariateDistribution)

Get the number of categories.
"""
ncategories(d::UnivariateDistribution)

"""
    ntrials(d::UnivariateDistribution)

Get the number of trials.
"""
ntrials(d::UnivariateDistribution)

"""
    dof(d::UnivariateDistribution)

Get the degrees of freedom.
"""
dof(d::UnivariateDistribution)

"""
    minimum(d::UnivariateDistribution)

Return the minimum of the support of `d`.
"""
minimum(d::UnivariateDistribution)

"""
    maximum(d::UnivariateDistribution)

Return the maximum of the support of `d`.
"""
maximum(d::UnivariateDistribution)

"""
    extrema(d::UnivariateDistribution)

Return the minimum and maximum of the support of `d` as a 2-tuple.
"""
extrema(d::UnivariateDistribution) = (minimum(d), maximum(d))

"""
    insupport(d::UnivariateDistribution, x::Any)

When `x` is a scalar, it returns whether x is within the support of `d`
(e.g., `insupport(d, x) = minimum(d) <= x <= maximum(d)`).
When `x` is an array, it returns whether every element in x is within the support of `d`.

Generic fallback methods are provided, but it is often the case that `insupport` can be
done more efficiently, and a specialized `insupport` is thus desirable.
You should also override this function if the support is composed of multiple disjoint intervals.
"""
insupport{D<:UnivariateDistribution}(d::Union{D, Type{D}}, x::Any)

function insupport!(r::AbstractArray, d::Union{D,Type{D}}, X::AbstractArray) where D<:UnivariateDistribution
    length(r) == length(X) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    for i in 1 : length(X)
        @inbounds r[i] = insupport(d, X[i])
    end
    return r
end

insupport(d::Union{D,Type{D}}, X::AbstractArray) where {D<:UnivariateDistribution} =
     insupport!(BitArray(undef, size(X)), d, X)

insupport(d::Union{D,Type{D}},x::Real) where {D<:ContinuousUnivariateDistribution} = minimum(d) <= x <= maximum(d)
insupport(d::Union{D,Type{D}},x::Real) where {D<:DiscreteUnivariateDistribution} = isinteger(x) && minimum(d) <= x <= maximum(d)

support(d::Union{D,Type{D}}) where {D<:ContinuousUnivariateDistribution} = RealInterval(minimum(d), maximum(d))
support(d::Union{D,Type{D}}) where {D<:DiscreteUnivariateDistribution} = round(Int, minimum(d)):round(Int, maximum(d))

# Type used for dispatch on finite support
# T = true or false
struct FiniteSupport{T} end

## macros to declare support

macro distr_support(D, lb, ub)
    D_has_constantbounds = (isa(ub, Number) || ub == :Inf) &&
                           (isa(lb, Number) || lb == :(-Inf))

    paramdecl = D_has_constantbounds ? :(d::Union{$D, Type{$D}}) : :(d::$D)

    # overall
    esc(quote
        minimum($(paramdecl)) = $lb
        maximum($(paramdecl)) = $ub
    end)
end

##### generic methods (fallback) #####

## sampling

# multiple univariate, must allocate array
rand(rng::AbstractRNG, s::Sampleable{Univariate}, dims::Dims) =
    rand!(rng, sampler(s), Array{eltype(s)}(undef, dims))

# multiple univariate with pre-allocated array
function rand!(rng::AbstractRNG, s::Sampleable{Univariate}, A::AbstractArray)
    smp = sampler(s)
    for i in eachindex(A)
        @inbounds A[i] = rand(rng, smp)
    end
    return A
end

"""
    rand(rng::AbstractRNG, d::UnivariateDistribution)

Generate a scalar sample from `d`. The general fallback is `quantile(d, rand())`.
"""
rand(rng::AbstractRNG, d::UnivariateDistribution) = quantile(d, rand(rng))

"""
    rand!(rng::AbstractRNG, ::UnivariateDistribution, ::AbstractArray)

Sample a univariate distribution and store the results in the provided array.
"""
rand!(rng::AbstractRNG, ::UnivariateDistribution, ::AbstractArray)

## statistics

"""
    mean(d::UnivariateDistribution)

Compute the expectation.
"""
mean(d::UnivariateDistribution)

"""
    var(d::UnivariateDistribution)

Compute the variance. (A generic std is provided as `std(d) = sqrt(var(d))`)
"""
var(d::UnivariateDistribution)

"""
    std(d::UnivariateDistribution)

Return the standard deviation of distribution `d`, i.e. `sqrt(var(d))`.
"""
std(d::UnivariateDistribution) = sqrt(var(d))

"""
    median(d::UnivariateDistribution)

Return the median value of distribution `d`.
"""
median(d::UnivariateDistribution) = quantile(d, 0.5)

"""
    modes(d::UnivariateDistribution)

Get all modes (if this makes sense).
"""
modes(d::UnivariateDistribution) = [mode(d)]

"""
    mode(d::UnivariateDistribution)

Returns the first mode.
"""
mode(d::UnivariateDistribution)

"""
    skewness(d::UnivariateDistribution)

Compute the skewness.
"""
skewness(d::UnivariateDistribution)

"""
    entropy(d::UnivariateDistribution)

Compute the entropy value of distribution `d`.
"""
entropy(d::UnivariateDistribution)

"""
    entropy(d::UnivariateDistribution, b::Real)

Compute the entropy value of distribution `d`, w.r.t. a given base.
"""
entropy(d::UnivariateDistribution, b::Real) = entropy(d) / log(b)

"""
    isplatykurtic(d)

Return whether `d` is platykurtic (*i.e* `kurtosis(d) < 0`).
"""
isplatykurtic(d::UnivariateDistribution) = kurtosis(d) < 0.0

"""
    isleptokurtic(d)

Return whether `d` is leptokurtic (*i.e* `kurtosis(d) > 0`).
"""
isleptokurtic(d::UnivariateDistribution) = kurtosis(d) > 0.0

"""
    ismesokurtic(d)

Return whether `d` is mesokurtic (*i.e* `kurtosis(d) == 0`).
"""
ismesokurtic(d::UnivariateDistribution) = kurtosis(d) ≈ 0.0

"""
    kurtosis(d::UnivariateDistribution)

Compute the excessive kurtosis.
"""
kurtosis(d::UnivariateDistribution)

"""
    kurtosis(d::Distribution, correction::Bool)

Computes excess kurtosis by default. Proper kurtosis can be returned with correction=false
"""
function kurtosis(d::Distribution, correction::Bool)
    if correction
        return kurtosis(d)
    else
        return kurtosis(d) + 3.0
    end
end

excess(d::Distribution) = kurtosis(d)
excess_kurtosis(d::Distribution) = kurtosis(d)
proper_kurtosis(d::Distribution) = kurtosis(d, false)

"""
    mgf(d::UnivariateDistribution, t)

Evaluate the moment generating function of distribution `d`.
"""
mgf(d::UnivariateDistribution, t)

"""
    cf(d::UnivariateDistribution, t)

Evaluate the characteristic function of distribution `d`.
"""
cf(d::UnivariateDistribution, t)

#### pdf, cdf, and friends

# pdf

"""
    pdf(d::UnivariateDistribution, x::Real)

Evaluate the probability density (mass) at `x`.

See also: [`logpdf`](@ref).
"""
pdf(d::UnivariateDistribution, x::Real)

"""
    pdf(d::DiscreteUnivariateDistribution, x::T) where {T<:Real}

Evaluate the probability density (mass) at `x`. If `T` is not an `Integer`
type but `x` is integer, the value is converted to `Int`.

The version with `x::Integer` must be implemented by
discrete distributions.

See also: [`logpdf`](@ref).
"""
pdf(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) ? pdf(d, round(Int, x)) : zero(partype(d))

"""
    logpdf(d::UnivariateDistribution, x::Real)

Evaluate the logarithm of probability density (mass) at `x`.
Whereas there is a fallback implemented `logpdf(d, x) = log(pdf(d, x))`.
Relying on this fallback is not recommended in general, as it is prone to overflow or underflow.
"""
logpdf(d::UnivariateDistribution, x::Real) = log(pdf(d, x))
logpdf(d::DiscreteUnivariateDistribution, x::Integer) = log(pdf(d, x))

"""
    cdf(d::UnivariateDistribution, x::Real)

Evaluate the cumulative probability at `x`.

See also [`ccdf`](@ref), [`logcdf`](@ref), and [`logccdf`](@ref).
"""
cdf(d::UnivariateDistribution, x::Real)
cdf(d::DiscreteUnivariateDistribution, x::Integer) = cdf(d, x, FiniteSupport{hasfinitesupport(d)})

# Discrete univariate with infinite support
function cdf(d::DiscreteUnivariateDistribution, x::Integer, ::Type{FiniteSupport{false}})
    c = 0.0
    for y = minimum(d):x
        c += pdf(d, y)
    end
    return c
end

# Discrete univariate with finite support
function cdf(d::DiscreteUnivariateDistribution, x::Integer, ::Type{FiniteSupport{true}})
    # calculate from left if x < (min + max)/2
    # (same as infinite support version)
    x <= div(minimum(d) + maximum(d),2) && return cdf(d, x, FiniteSupport{false})

    # otherwise, calculate from the right
    c = 1.0
    for y = x+1:maximum(d)
        c -= pdf(d, y)
    end
    return c
end

cdf(d::DiscreteUnivariateDistribution, x::Real) = cdf(d, floor(Int,x))
cdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(cdf, (d, x)))

"""
    ccdf(d::UnivariateDistribution, x::Real)

The complementary cumulative function evaluated at `x`, i.e. `1 - cdf(d, x)`.
"""
ccdf(d::UnivariateDistribution, x::Real) = 1.0 - cdf(d, x)
ccdf(d::DiscreteUnivariateDistribution, x::Integer) = 1.0 - cdf(d, x)
ccdf(d::DiscreteUnivariateDistribution, x::Real) = ccdf(d, floor(Int,x))

"""
    logcdf(d::UnivariateDistribution, x::Real)

The logarithm of the cumulative function value(s) evaluated at `x`, i.e. `log(cdf(x))`.
"""
logcdf(d::UnivariateDistribution, x::Real) = log(cdf(d, x))
logcdf(d::DiscreteUnivariateDistribution, x::Integer) = log(cdf(d, x))
logcdf(d::DiscreteUnivariateDistribution, x::Real) = logcdf(d, floor(Int,x))

"""
    logccdf(d::UnivariateDistribution, x::Real)

The logarithm of the complementary cumulative function values evaluated at x, i.e. `log(ccdf(x))`.
"""
logccdf(d::UnivariateDistribution, x::Real) = log(ccdf(d, x))
logccdf(d::DiscreteUnivariateDistribution, x::Integer) = log(ccdf(d, x))
logccdf(d::DiscreteUnivariateDistribution, x::Real) = logccdf(d, floor(Int,x))

"""
    quantile(d::UnivariateDistribution, q::Real)

Evaluate the inverse cumulative distribution function at `q`.

See also: [`cquantile`](@ref), [`invlogcdf`](@ref), and [`invlogccdf`](@ref).
"""
quantile(d::UnivariateDistribution, p::Real)

"""
    cquantile(d::UnivariateDistribution, q::Real)

The complementary quantile value, i.e. `quantile(d, 1-q)`.
"""
cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p)

"""
    invlogcdf(d::UnivariateDistribution, lp::Real)

The inverse function of logcdf.
"""
invlogcdf(d::UnivariateDistribution, lp::Real) = quantile(d, exp(lp))

"""
    invlogcdf(d::UnivariateDistribution, lp::Real)

The inverse function of logcdf.
"""
invlogccdf(d::UnivariateDistribution, lp::Real) = quantile(d, -expm1(lp))

# gradlogpdf

gradlogpdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(gradlogpdf, (d, x)))

function _pdf_fill_outside!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
    vl = vfirst = first(X)
    vr = vlast = last(X)
    n = vlast - vfirst + 1
    if islowerbounded(d)
        lb = minimum(d)
        if vl < lb
            vl = lb
        end
    end
    if isupperbounded(d)
        ub = maximum(d)
        if vr > ub
            vr = ub
        end
    end

    # fill left part
    if vl > vfirst
        for i = 1:(vl - vfirst)
            r[i] = 0.0
        end
    end

    # fill central part: with non-zero pdf
    fm1 = vfirst - 1
    for v = vl:vr
        r[v - fm1] = pdf(d, v)
    end

    # fill right part
    if vr < vlast
        for i = (vr-vfirst+2):n
            r[i] = 0.0
        end
    end
    return vl, vr, vfirst, vlast
end

function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)

    # fill central part: with non-zero pdf
    fm1 = vfirst - 1
    for v = vl:vr
        r[v - fm1] = pdf(d, v)
    end
    return r
end

abstract type RecursiveProbabilityEvaluator end

function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange, rpe::RecursiveProbabilityEvaluator)
    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)

    # fill central part: with non-zero pdf
    if vl <= vr
        fm1 = vfirst - 1
        r[vl - fm1] = pv = pdf(d, vl)
        for v = (vl+1):vr
            r[v - fm1] = pv = nextpdf(rpe, pv, v)
        end
    end

    return r
end

## loglikelihood
"""
    loglikelihood(d::UnivariateDistribution, X::AbstractArray)

The log-likelihood of distribution `d` w.r.t. all samples contained in array `x`.
"""
loglikelihood(d::UnivariateDistribution, X::AbstractArray) = sum(x -> logpdf(d, x), X)

### macros to use StatsFuns for method implementation

macro _delegate_statsfuns(D, fpre, psyms...)
    dt = eval(D)
    T = dt <: DiscreteUnivariateDistribution ? :Int : :Real

    # function names from StatsFuns
    fpdf = Symbol(fpre, "pdf")
    flogpdf = Symbol(fpre, "logpdf")
    fcdf = Symbol(fpre, "cdf")
    fccdf = Symbol(fpre, "ccdf")
    flogcdf = Symbol(fpre, "logcdf")
    flogccdf = Symbol(fpre, "logccdf")
    finvcdf = Symbol(fpre, "invcdf")
    finvccdf = Symbol(fpre, "invccdf")
    finvlogcdf = Symbol(fpre, "invlogcdf")
    finvlogccdf = Symbol(fpre, "invlogccdf")

    # parameter fields
    pargs = [Expr(:(.), :d, Expr(:quote, s)) for s in psyms]

    esc(quote
        pdf(d::$D, x::$T) = $(fpdf)($(pargs...), x)
        logpdf(d::$D, x::$T) = $(flogpdf)($(pargs...), x)

        cdf(d::$D, x::$T) = $(fcdf)($(pargs...), x)
        ccdf(d::$D, x::$T) = $(fccdf)($(pargs...), x)
        logcdf(d::$D, x::$T) = $(flogcdf)($(pargs...), x)
        logccdf(d::$D, x::$T) = $(flogccdf)($(pargs...), x)

        quantile(d::$D, q::Real) = convert($T, $(finvcdf)($(pargs...), q))
        cquantile(d::$D, q::Real) = convert($T, $(finvccdf)($(pargs...), q))
        invlogcdf(d::$D, lq::Real) = convert($T, $(finvlogcdf)($(pargs...), lq))
        invlogccdf(d::$D, lq::Real) = convert($T, $(finvlogccdf)($(pargs...), lq))
    end)
end

##### specific distributions #####

const discrete_distributions = [
    "bernoulli",
    "betabinomial",
    "binomial",
    "discreteuniform",
    "discretenonparametric",
    "categorical",
    "geometric",
    "hypergeometric",
    "negativebinomial",
    "noncentralhypergeometric",
    "poisson",
    "skellam",
    "poissonbinomial"
]

const continuous_distributions = [
    "arcsine",
    "beta",
    "betaprime",
    "biweight",
    "cauchy",
    "chernoff",
    "chisq",    # Chi depends on Chisq
    "chi",
    "cosine",
    "epanechnikov",
    "exponential",
    "fdist",
    "frechet",
    "gamma", "erlang",
    "pgeneralizedgaussian", # GeneralizedGaussian depends on Gamma
    "generalizedpareto",
    "generalizedextremevalue",
    "gumbel",
    "inversegamma",
    "inversegaussian",
    "kolmogorov",
    "ksdist",
    "ksonesided",
    "laplace",
    "levy",
    "locationscale",
    "logistic",
    "noncentralbeta",
    "noncentralchisq",
    "noncentralf",
    "noncentralt",
    "normal",
    "normalcanon",
    "normalinversegaussian",
    "lognormal",    # LogNormal depends on Normal
    "logitnormal",    # LogitNormal depends on Normal
    "pareto",
    "rayleigh",
    "semicircle",
    "studentizedrange",
    "symtriangular",
    "tdist",
    "triangular",
    "triweight",
    "uniform",
    "vonmises",
    "weibull"
]

for dname in discrete_distributions
    include(joinpath("univariate", "discrete", "$(dname).jl"))
end

for dname in continuous_distributions
    include(joinpath("univariate", "continuous", "$(dname).jl"))
end


text:
#using Modia
@static if VERSION < v"0.7.0-DEV.2005"
  using Base.Test
else
  using Test
end

@testset "RunExamples" begin

Modia.ModiaLogging.resetTestStatus()

global figure = 1
include("CurrentController.jl")
include("Rectifier.jl")
@static if ! (VERSION < v"0.7.0-DEV.2005")
  include("CauerLowPassFilter.jl")
end
include("LinearSystems.jl")
include("SynchronousExamples.jl")
include("ElectricalVehicleAndCharger.jl")
include("CollidingBalls.jl")
include("HeatTransfer2D.jl")

Modia.ModiaLogging.printTestStatus()

end # testset
nothing


text:
"""
    annotations(strings::Vector{String}, positions::Vector{Point})

Plots an array of texts at each position in `positions`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Annotations, text, position) do scene
    default_theme(scene, Text)
end

function convert_arguments(::Type{<: Annotations},
                           strings::AbstractVector{<: AbstractString},
                           text_positions::AbstractVector{<: Point{N}}) where N
    return (map(strings, text_positions) do str, pos
        (String(str), Point{N, Float32}(pos))
    end,)
end

function plot!(plot::Annotations)
    sargs = (
        plot.model, plot.font,
        plot[1],
        getindex.(plot, (:color, :textsize, :align, :rotation, :justification, :lineheight))...,
    )
    atlas = get_texture_atlas()
    combinedpos = [Point3f0(0)]
    colors = RGBAf0[RGBAf0(0,0,0,0)]
    textsize = Float32[0]
    fonts = [defaultfont()]
    rotations = Quaternionf0[Quaternionf0(0,0,0,0)]

    tplot = text!(plot, " ",
        align = Vec2f0(0), model = Mat4f0(I),
        position = combinedpos, color = colors, visible = plot.visible,
        textsize = textsize, font = fonts, rotation = rotations
    )

    onany(sargs...) do model, pfonts, text_pos, args...
        io = IOBuffer();
        empty!(combinedpos); empty!(colors); empty!(textsize); empty!(fonts); empty!(rotations)
        broadcast_foreach(1:length(text_pos), to_font(pfonts), text_pos, args...) do idx, f,
                (text, startpos), color, tsize, alignment, rotation, justification, lineheight
            c = to_color(color)
            rot = to_rotation(rotation)
            pos = layout_text(text, startpos, tsize, f, alignment, rot, model, justification, lineheight)
            print(io, text)
            n = length(pos)
            append!(combinedpos, pos)
            append!(textsize, repeated(tsize, n))
            append!(colors, repeated(c, n))
            append!(fonts, one_attribute_per_char(f, text))
            append!(rotations, repeated(rot, n))
        end
        str = String(take!(io))
        # update string the signals
        tplot[1] = str
        return
    end
    # update one time in the beginning, since otherwise the above won't run
    notify!(sargs[1])
    plot
end


text:
# ---
# title: 1467. Probability of a Two Boxes Having The Same Number of Distinct Balls
# id: problem1467
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Backtracking
# link: <https://leetcode.com/problems/probability-of-a-two-boxes-having-the-same-number-of-distinct-balls/description/>
# hidden: true
# ---
#
# Given `2n` balls of `k` distinct colors. You will be given an integer array
# `balls` of size `k` where `balls[i]` is the number of balls of color `i`.
#
# All the balls will be **shuffled uniformly at random** , then we will
# distribute the first `n` balls to the first box and the remaining `n` balls to
# the other box (Please read the explanation of the second example carefully).
#
# Please note that the two boxes are considered different. For example, if we
# have two balls of colors `a` and `b`, and two boxes `[]` and `()`, then the
# distribution `[a] (b)` is considered different than the distribution `[b] (a)
# `(Please read the explanation of the first example carefully).
#
# We want to _calculate the probability_ that the two boxes have the same number
# of distinct balls.
#
#
#
# **Example 1:**
#
#
#
#     Input: balls = [1,1]
#     Output: 1.00000
#     Explanation: Only 2 ways to divide the balls equally:
#     - A ball of color 1 to box 1 and a ball of color 2 to box 2
#     - A ball of color 2 to box 1 and a ball of color 1 to box 2
#     In both ways, the number of distinct colors in each box is equal. The probability is 2/2 = 1
#
#
# **Example 2:**
#
#
#
#     Input: balls = [2,1,1]
#     Output: 0.66667
#     Explanation: We have the set of balls [1, 1, 2, 3]
#     This set of balls will be shuffled randomly and we may have one of the 12 distinct shuffles with equale probability (i.e. 1/12):
#     [1,1 / 2,3], [1,1 / 3,2], [1,2 / 1,3], [1,2 / 3,1], [1,3 / 1,2], [1,3 / 2,1], [2,1 / 1,3], [2,1 / 3,1], [2,3 / 1,1], [3,1 / 1,2], [3,1 / 2,1], [3,2 / 1,1]
#     After that we add the first two balls to the first box and the second two balls to the second box.
#     We can see that 8 of these 12 possible random distributions have the same number of distinct colors of balls in each box.
#     Probability is 8/12 = 0.66667
#
#
# **Example 3:**
#
#
#
#     Input: balls = [1,2,1,2]
#     Output: 0.60000
#     Explanation: The set of balls is [1, 2, 2, 3, 4, 4]. It is hard to display all the 180 possible random shuffles of this set but it is easy to check that 108 of them will have the same number of distinct colors in each box.
#     Probability = 108 / 180 = 0.6
#
#
# **Example 4:**
#
#
#
#     Input: balls = [3,2,1]
#     Output: 0.30000
#     Explanation: The set of balls is [1, 1, 1, 2, 2, 3]. It is hard to display all the 60 possible random shuffles of this set but it is easy to check that 18 of them will have the same number of distinct colors in each box.
#     Probability = 18 / 60 = 0.3
#
#
# **Example 5:**
#
#
#
#     Input: balls = [6,6,6,6,6,6]
#     Output: 0.90327
#
#
#
#
# **Constraints:**
#
#   * `1 <= balls.length <= 8`
#   * `1 <= balls[i] <= 6`
#   * `sum(balls)` is even.
#   * Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

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


text:
# ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ #
#  Verification   #
# ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ #

function test_GKGT_random( ; tol=1e-1)

    t1, t2    = randn(2)*3
    ℓ         = rand()*(4.95 - 0.05) + 0.05
    w1        = rand()*(10_000 - 3_000) + 3_000.0
    mass      = exp(rand()*(log(1e9) - log(1e6)) + log(1e6))
    accretion = rand()*0.9 + 0.05

    test_GKGT(; t1=t1, t2=t2, ℓ=ℓ, w1=w1, mass=mass,
                accretion=accretion, tol=tol)

end

function test_GKGT(;t1=t1,t2=t2,ℓ=ℓ,w1=w1,mass=mass,accretion=accretion, tol=tol)

    @printf("times: t1 = %.4f, t2 = %.4f\n", t1, t2)
    @printf("ℓ = %.4f\n", ℓ)
    @printf("wavelength = %.4f\n", w1)
    @printf("mass = %.4e\n", mass)
    @printf("accretion = %.4f\n", accretion)

    k1 = symGKGᵀ(mass,accretion, w1, true)
    k2 = symGKGᵀ(mass,accretion, w1, false)

    # all results obtained below must agree with each other

    a1 = k1(t1, t2, ℓ)

    a11 = k2(t2, t1, ℓ)

    a2 = test_approx_physical(t1, t2, ℓ; mass = mass, accretion = accretion, wavelengths=[w1; w1])

    a3 = test_nonsymmetricconvolution(t1, t2, ℓ; mass = mass, accretion = accretion, wavelengths=[w1; w1])

    a4 = test_symmetricconvolution_cuba(t1, t2, ℓ; mass = mass, accretion = accretion, wavelength = w1)

    display([a1; a11; a2; a3; a4])

    (abs(a1 - a11) < tol) &&
    (abs(a1 - a2)  < tol) &&
    (abs(a1 - a3)  < tol) &&
    (abs(a1 - a4)  < tol) &&
    (abs(a2 - a11) < tol) &&
    (abs(a2 - a3)  < tol) &&
    (abs(a2 - a4)  < tol) &&
    (abs(a3 - a11) < tol) &&
    (abs(a3 -  a4) < tol)

end


text:
export PseudoHilbert,Hilbert,SingularIntegral

#############
# Hilbert implements the Hilbert operator as a contour integral:
#
#        H f(z) := 1/π\int_Γ f(t)/(t-z) dt,  z ∈ Γ
#
# SingularIntegral implements the Hilbert operator as a line integral:
#
#       SI f(z) := 1/π\int_Γ f(t)/(t-z) ds(t),  z ∉ Γ
#
# PseudoHilbert corresponds to pseudocauchy, which may not be normalized at
# infinity.
#############

@calculus_operator(PseudoHilbert)
@calculus_operator(Hilbert)
@calculus_operator(SingularIntegral)

for Op in (:PseudoHilbert,:Hilbert,:SingularIntegral)
    ConcOp=Meta.parse("Concrete"*string(Op))
    OpWrap=Meta.parse(string(Op)*"Wrapper")
    OffOp=Meta.parse("Off"*string(Op))
    @eval begin
        ## Convenience routines
        $Op(d::IntervalOrSegmentDomain, n::Int) =
            $Op(JacobiWeight(-.5,-.5,Chebyshev(d)),n)
        $Op(d::IntervalOrSegmentDomain) =
            $Op(JacobiWeight(-.5,-.5,Chebyshev(d)))
        $Op(d::PeriodicDomain,n::Int) = $Op(Laurent(d),n)
        $Op(d::PeriodicDomain) = $Op(Laurent(d))

        ## Modifiers for SumSpace, ArraySpace, and PiecewiseSpace
        function $Op(S::PiecewiseSpace,n::Int)
            sp = components(S)
            m = length(sp)
            diag = Any[$Op(sp[k],n) for k=1:m]
            C = Any[k==j ? diag[k] : $OffOp(sp[k],rangespace(diag[j]),n) for j=1:m,k=1:m]
            D = Operator{mapreduce(i->eltype(C[i]),promote_type,1:m^2)}[C[j,k] for j=1:m,k=1:m]
            D = promotespaces(D)
            $OpWrap(InterlaceOperator(D,S,PiecewiseSpace(map(rangespace,D[:,1]))),n)
        end

        bandwidths(::$ConcOp{Hardy{s,DD,RR}}) where {s,DD,RR} = 0,0
        domainspace(H::$ConcOp{Hardy{s,DD,RR}}) where {s,DD,RR} = H.space
        rangespace(H::$ConcOp{Hardy{s,DD,RR}}) where {s,DD,RR} = H.space

        bandwidths(::$ConcOp{Laurent{DD,RR}}) where {DD,RR} = 0,0
        domainspace(H::$ConcOp{Laurent{DD,RR}}) where {DD,RR} = H.space
        rangespace(H::$ConcOp{Laurent{DD,RR}}) where {DD,RR} = H.space

        bandwidths(H::$ConcOp{Fourier{DD,RR}}) where {DD,RR} = H.order,H.order
        domainspace(H::$ConcOp{Fourier{DD,RR}}) where {DD,RR} = H.space
        rangespace(H::$ConcOp{Fourier{DD,RR}}) where {DD,RR} = H.space

        function rangespace(H::$ConcOp{<:JacobiWeight{<:Chebyshev}})
            @assert domainspace(H).α == domainspace(H).β == -0.5
            H.order==0 ? Chebyshev(domain(H)) : Ultraspherical(H.order,domain(H))
        end
        function rangespace(H::$ConcOp{<:JacobiWeight{<:Ultraspherical{Int}}})
            @assert order(domainspace(H).space) == 1
            @assert domainspace(H).α==domainspace(H).β==0.5
            (H.order==1||H.order==0) ? Chebyshev(domain(H)) : Ultraspherical(H.order-1,domain(H))
        end
        # bandwidths{λ,DD,RR}(H::$ConcOp{JacobiWeight{Ultraspherical{λ,DD,RR},DD,RR}})=-λ,H.order-λ
        bandwidths(H::$ConcOp{<:JacobiWeight{<:Chebyshev}}) = 0,H.order
        bandwidths(H::$ConcOp{<:JacobiWeight{<:Ultraspherical{Int}}}) =
            H.order > 0 ? (1,H.order-1) : (2,0)

        choosedomainspace(H::$Op{UnsetSpace}, sp::Ultraspherical) =
            ChebyshevWeight(ChebyshevDirichlet{1,1}(domain(sp)))
        choosedomainspace(H::$Op{UnsetSpace}, sp::Chebyshev) =
            ChebyshevWeight(ChebyshevDirichlet{1,1}(domain(sp)))
        choosedomainspace(H::$Op{UnsetSpace}, sp::PiecewiseSpace) =
            PiecewiseSpace(map(s->choosedomainspace(H,s),sp.spaces))
    end
end

for TYP in (:SumSpace,:PiecewiseSpace),(Op,OpWrap) in ((:PseudoHilbert,:PseudoHilbertWrapper),
                          (:Hilbert,:HilbertWrapper),
                          (:SingularIntegral,:SingularIntegralWrapper))
    @eval $Op(S::$TYP,k) = $OpWrap(InterlaceOperator(Diagonal([map(s->$Op(s,k),S.spaces)...]),$TYP),k)
end

for (Op,OpWrap) in ((:PseudoHilbert,:PseudoHilbertWrapper),
                          (:Hilbert,:HilbertWrapper),
                          (:SingularIntegral,:SingularIntegralWrapper))
    @eval function $Op(S::ArraySpace,k)
        ops = map(s->$Op(s,k),vec(S.spaces))
        $OpWrap(InterlaceOperator(Diagonal(ops),S,ArraySpace(reshape(map(rangespace,ops),size(S)))),k)
    end
end

# Length catch

ConcreteHilbert(sp::Space{D,R},n) where {D,R<:Complex} =
    ConcreteHilbert{typeof(sp),typeof(n),prectype(sp)}(sp,n)
ConcreteSingularIntegral(sp::Space{D,R},n) where {D,R<:Complex} =
    ConcreteSingularIntegral{typeof(sp),typeof(n),prectype(sp)}(sp,n)

# Override sumspace

for TYP in (:Hilbert,:SingularIntegral)
    ConcOp=Meta.parse("Concrete"*string(TYP))
    @eval function $TYP(F::Fourier{DD,RR},n) where {DD<:Circle,RR}
        if !domain(F).orientation
            R=reverseorientation(F)
            Conversion(R,F)*(-$TYP(R,n))*Conversion(F,R)
        else
            $ConcOp(F,n)
        end
    end
end

### Operator Entries

## Circle

for Typ in (:Hilbert,:SingularIntegral)
    ConcTyp = Meta.parse("Concrete"*string(Typ))
    WrapTyp = Meta.parse(string(Typ)*"Wrapper")
    @eval begin
        $Typ(S::Hardy{s,DD,RR},m::Integer) where {s,DD<:Circle,RR} =
            m ≤ 1 ? $ConcTyp(S,m) : $WrapTyp(Derivative(m-1)*Hilbert(S),m)
        $Typ(S::Laurent{DD,RR},m::Integer) where {DD<:Circle,RR} =
            m ≤ 1 ? $ConcTyp(S,m) : $WrapTyp(Derivative(m-1)*Hilbert(S),m)
        $Typ(S::Fourier{DD,RR},m::Integer) where {DD<:Circle,RR} =
            m ≤ 1 ? $ConcTyp(S,m) : $WrapTyp(Derivative(m-1)*Hilbert(S),m)
    end
end

function getindex(H::ConcreteHilbert{Hardy{true,DD,RR},OT,T},
                  k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
    ##TODO: Add scale for different radii.
    m = H.order
    d = domain(H)
    r = d.radius

    if k == j
        if m==0
            k == 1 ? -T(2r*log(r)) : T(r/(k-1))
        else
            im*(im*(k-1)/r)^(m-1)
        end
    else
        zero(T)
    end
end

function getindex(H::ConcreteHilbert{Hardy{false,DD,RR},OT,T},
                  k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
    ##TODO: Add scale for different radii.
    m = H.order
    d = domain(H)
    r = d.radius

    if k == j
        -im*(im*k/r)^(m-1)
    else
        zero(T)
    end
end

function getindex(H::ConcreteHilbert{Laurent{DD,RR},OT,T},
                  k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
    ##TODO: Add scale for different radii.
    m = H.order
    d = domain(H)
    r = d.radius

    if k == j
        if m == 0
            if k == 1
                -T(2r*log(r))
            elseif isodd(k)
                T(r/((k-1)÷2))
            else
                T(-r/(k÷2))
            end
        else
            if isodd(k)
                T(im*(im*((k-1)÷2)/r)^(m-1))
            else
                -T(im*(im*(k÷2)/r)^(m-1))
            end
        end
    else
        zero(T)
    end
end

function getindex(H::ConcreteHilbert{Fourier{DD,RR},OT,T},k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
    m = H.order
    d = domain(H)
    r = d.radius
    o = d.orientation

    if m == 0 && k==j
        T(k==1 ? 2r*log(r) : (-r/(k÷2)))
    elseif m == 1
        if k==j==1
            T(o ? im : -im)
        elseif iseven(k) && j==k+1
            -one(T)
        elseif isodd(k) && j==k-1
            one(T)
        else
            zero(T)
        end
    else
        zero(T)
    end
end

function getindex(H::ConcreteSingularIntegral{Hardy{true,DD,RR},OT,T},k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
##TODO: Add scale for different radii.
    m = H.order
    d = domain(H)
    r = d.radius

    if k == j
        if m == 0
            k == 1 ? T(2r*log(r)) : T(-r/(k-1))
        else
            k == 1 ? zero(T) : im*(im*(k-1)/r)^(m-1)
        end
    else
        zero(T)
    end
end

function getindex(H::ConcreteSingularIntegral{Hardy{false,DD,RR},OT,T},k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
##TODO: Add scale for different radii.
    m = H.order
    d = domain(H)
    r = d.radius

    if k == j
        -im*(im*k/r)^(m-1)
    else
        zero(T)
    end
end

function getindex(H::ConcreteSingularIntegral{Laurent{DD,RR},OT,T},k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
    ##TODO: Add scale for different radii.
    m = H.order
    d = domain(H)
    r = d.radius

    if k == j
        if m == 0
            if k == 1
                T(2r*log(r))
            else
                T(-r/(k÷2))
            end
        else
            if isodd(k)
                k == 1 ? zero(T) : im*(im*((k-1)÷2)/r)^(m-1)
            else
                -im*(im*(k÷2)/r)^(m-1)
            end
        end
    else
        zero(T)
    end
end

function getindex(H::ConcreteSingularIntegral{Fourier{DD,RR},OT,T},k::Integer,j::Integer) where {DD<:Circle,RR,OT,T}
    m = H.order
    d = domain(H)
    r = d.radius

    if m == 0
        if k == j
            k == 1 ? T(2r*log(r)) : T(-r/(k÷2))
        end
    elseif m == 1
        if k+1 == j
            -one(T)
        elseif k-1 == j
            one(T)
        else
            zero(T)
        end
    else
        error("SingularIntegral order $(H.order) not implemented for Fourier")
    end
end

## JacobiWeight

for (Op,Len) in ((:Hilbert,:complexlength),
                        (:SingularIntegral,:arclength))
    ConcOp=Meta.parse("Concrete"*string(Op))
    OpWrap=Meta.parse(string(Op)*"Wrapper")

    @eval begin
        function $Op(S::JacobiWeight{Chebyshev{DD,RR},DD},m::Int) where {DD<:IntervalOrSegment,RR}
            if S.α==S.β==-0.5
                if m==0
                    $ConcOp(S,m)
                else
                    d=domain(S)
                    C=(4/$Len(d))^(m-1)
                    $OpWrap(SpaceOperator(
                        ToeplitzOperator([0.],[zeros(m);C]),S,Ultraspherical(m,d)),m)
                end
            elseif S.α==S.β==0.5
                d=domain(S)
                if m==1
                    J=JacobiWeight(0.5,0.5,Ultraspherical(1,d))
                    $OpWrap($Op(J,m)*Conversion(S,J),m)
                else
                    J=JacobiWeight(-0.5,-0.5,Chebyshev(d))
                    $OpWrap($Op(J,m)*Conversion(S,J),m)
                end
            else
                error(string($Op)*" not implemented for parameters $(S.α),$(S.β)")
            end
        end

        function getindex(H::$ConcOp{<:JacobiWeight{<:Chebyshev{<:IntervalOrSegment}},OT,T},k::Integer,j::Integer) where {OT,T}
            sp=domainspace(H)
            @assert H.order == 0
            @assert sp.α==sp.β==-0.5

            if k==j
                d=domain(H)
                C=$Len(d)/4
                -T(k==1 ? -2C*log(C) : 2C/(k-1))
            else
                zero(T)
            end
        end

        # we always have real for n==1
        function $Op(S::JacobiWeight{<:Ultraspherical{Int,<:IntervalOrSegment}},m)
            @assert order(S.space) == 1
            if S.α==S.β==0.5
                if m==1
                    d=domain(S)
                    $OpWrap(SpaceOperator(
                        ToeplitzOperator([-1.0],[0.]),S,m==1 ? Chebyshev(d) : Ultraspherical(m-1,d)),m)
                else
                    $ConcOp(S,m)
                end
            else
                error(string($Op)*" not implemented for parameters $(S.α),$(S.β)")
            end
        end
        function getindex(H::$ConcOp{<:JacobiWeight{<:Ultraspherical{Int,<:IntervalOrSegment}},OT,T},k::Integer,j::Integer) where {OT,T}
            # order(domainspace(H))==1
            m=H.order
            d=domain(H)
            sp=domainspace(H)
            @assert sp.α==sp.β==0.5
            @assert m ≠ 1

            C=(4/$Len(d))^(m-1)
            if m == 0
                if k==j==1
                    T(C*log(C))
                elseif k==j
                    -T(C/(k-1))
                elseif j==k-2
                    T(C/(k-1))
                else
                    zero(T)
                end
            elseif j==k+m-2
                -T(.5C*k/(m-1))
            else
                zero(T)
            end
        end
    end
end

## PseudoHilbert
# The default is Hilbert

# getindex(H::ConcretePseudoHilbert,k::Integer,j::Integer)=Hilbert(H.space,H.order)[k,j]
# bandwidths(H::ConcretePseudoHilbert)=bandwidths(Hilbert(H.space,H.order))


text:
@doc """
	newcorner!(corners::Tessellation, generator::IndexablePoint2D, corner::Point2D)

Update `corners` with a new `corner` of the cell belonging to a particular `generator`.
If `generator` is already in `corners`, the entry in `corners` is updated with `corner` and otherwise a new cell is added.
"""->
function newcorner!(corners::Tessellation, generator::IndexablePoint2D, corner::AbstractPoint2D)
	index = getindex(generator)

	if haskey( corners, index )
		if !contains(corners[index], corner)
			push!( corners[index], corner )
		end
	else
		corners[index] = [ corner ]
	end
end

@doc """
	newedge!(corners::Tessellation, edge::VoronoiEdge)

Update `corners` with the corners of `edge`.
See also `newcorner!`.
"""->
function newedge!(corners::Tessellation, edge::VoronoiDelaunay.VoronoiEdge{IndexablePoint2D})
	# TODO: Import edge type?

	# Clip edge to bounding box
	A = small2large(geta(edge))
	B = small2large(getb(edge))
	if !isinside(A) || !isinside(B)
		A, B = clip(A, B)
		if isa(A,Void) || isa(B,Void)
			return nothing
		end
	end

	generator = getgena(edge)
	newcorner!(corners, generator, A)
	newcorner!(corners, generator, B)

	generator = getgenb(edge)
	newcorner!(corners, generator, A)
	newcorner!(corners, generator, B)
end

@doc """
	voronoicells(generators::IndexablePoints2D) -> Tessellation

Collect the Voronoi cells from a set of `generators`.
"""->
function voronoicells(generators::IndexablePoints2D)
	# Transform points to the middle square
	tgen = large2small(generators)

	# VoronoiDelaunay data structure
	Ngen = length(tgen)
	tess = DelaunayTessellation2D{IndexablePoint2D}(Ngen)
	push!(tess, tgen)

	# Initialize output
	corn = Tessellation()
	sizehint!(corn, Ngen)

	# Dict for indices of quadrant neighbors
	Q = Dict{Int64, Vector{Int64}}(1=>[], 2=>[], 3=>[], 4=>[])

	for edge in voronoiedges(tess)
		newedge!(corn, edge)
		quadrant!(Q, edge)
	end

	# Add corners of bounding box to corn to cor
	add_bounding_corners!(corn, generators, Q)

	return corn
end

@doc """
`Q` is a `Dict` with keys 1 through 4 representing the four quadrants of the bounding box.
The points with cells that border quadrant `q` are `Q[q]`.

This function may include too many points, but not so many that it is important for the performance.
"""->
function quadrant!{T<:AbstractPoint2D}(Q::Dict{Int64, Vector{Int64}}, edge::VoronoiDelaunay.VoronoiEdge{T})
	# TODO: Should test if edge is outside the middle square/map edge to
	# full square. Mapping takes time and the clip functions needs to be
	# modified to test in the middle square.
	if isoutside(edge)
		return nothing
	end

	gena = getgena(edge)
	genb = getgenb(edge)
	if getindex(gena) == -1
		q = quadrant(gena)
		push!(Q[q], getindex(genb))
	end

	if getindex(genb) == -1
		q = quadrant(genb)
		push!(Q[q], getindex(gena))
	end
end

@doc """
	add_bounding_corners!(corn, generators, Q)

Adds each of the quadrants to the cell in the tesselation `corn` whose generator is closest.
`Q` holds the indices of the cells that border the bounding box corners (see `quadrant!`).
"""->
function add_bounding_corners!(corn::Tessellation, generators::IndexablePoints2D, Q)
	sort!(generators, by=getindex)

	for quad in keys(Q)
		gen_dist = Inf
		neighbor_index = -1
		for idx in Q[quad]
			D = dist_squared(BoxCorners[quad], generators[idx])
			if D < gen_dist
				gen_dist = D
				neighbor_index = idx
			end
		end

		newcorner!(corn, generators[neighbor_index], BoxCorners[quad])
	end
end

function large2small(p::IndexablePoint2D)
	IndexablePoint2D( 0.5*getx(p)+0.75, 0.5*gety(p)+0.75, getindex(p) )
end

function large2small(pts::Vector{IndexablePoint2D})
	[large2small(p) for p in pts]
end

function small2large(p::Point2D)
	Point2D( 2.0*getx(p)-1.5, 2.0*gety(p)-1.5 )
end


text:
# follows Conjugate gradient algorithm for optimization under unitary Traian Abrudan  ,1,2 , Jan Eriksson 2 , Visa Koivunen
# matrix constraint Signal Processing 89 (2009) 1704–1714

function step_size(w_k, h_k; q, euclidian_gradient, p = 5)

    """

        step_size(w_k, h_k; q, euclidian_gradient, p = 5)

    implements the geodesic search algorithm of table 1 in Conjugate gradient algorithm for optimization under unitary Traian Abrudan  ,1,2 , Jan Eriksson 2 , Visa Koivunen matrix constraint Signal Processing 89 (2009) 1704–1714"""

    n = size(w_k)[1]

    omega_max = maximum(abs.(eigvals(h_k)))

    t_mu = 2pi/(q*omega_max)

    mu = [i * t_mu/p for i in 0:p]

    r = Matrix{ComplexF64}[]

    push!(r, Matrix{ComplexF64}(I,n,n))
    push!(r, exp(-t_mu/p * h_k))
    for i in 2:p
        push!(r, r[2]*r[end])
    end

    derivatives = zeros(ComplexF64, p+1)

    for i = 0:p
        derivatives[i+1] = -2 * real(tr(euclidian_gradient(r[i+1]*w_k) * w_k'*r[i+1]'*h_k'))
    end

    a = zeros(ComplexF64, p+1)

    a[1] = derivatives[1]

    mat_mu = zeros(ComplexF64, p,p)

    for i = 1:p
        mat_mu[:, i] = mu[2:end].^i
    end

    a[2:end] = inv(mat_mu) * (derivatives[2:end].-a[1])

    positive_roots = []

    for root in roots(reverse(a))
        if imag(root) ≈ 0. atol = 1e-7
            if real(root) > 0.
                push!(positive_roots, real(root))
            end
        end
    end

    #step_size
    if positive_roots == []
        return 0.
    else
        return minimum(positive_roots)
    end

end

"""

    minimize_over_unitary_matrices(;euclidian_gradient, q, n, p=5, tol = 1e-6, max_iter = 100)

returns the (optimized matrix, optimization_success)

implements the minimization algorithm of table 3 in  Conjugate gradient algorithm for optimization under unitary Traian Abrudan  ,1,2 , Jan Eriksson 2 , Visa Koivunen matrix constraint Signal Processing 89 (2009) 1704–1714"""
function minimize_over_unitary_matrices(;euclidian_gradient, q, n, p=5, tol = 1e-6, max_iter = 100)

    optimization_success = false

    function inner_product(a,b)
        0.5 * real(tr(a'*b))
    end

    k = 0
    w_k = Matrix{ComplexF64}(I,n,n)
    g_k = similar(w_k)
    h_k = similar(w_k)

    while k < max_iter

        #println("iter $(k+1)/$max_iter")
        if k % n^2 == 0
            gamma_k = euclidian_gradient(w_k)
            g_k = gamma_k * w_k' - w_k * gamma_k'
            h_k = g_k
        end

        if inner_product(g_k,g_k) < tol
            #println("optimisation achieved")
            optimization_success = true
            break

        end

        w_k_plus_1 = exp(-step_size(w_k, h_k; q = q, euclidian_gradient = euclidian_gradient, p = p) * h_k)*w_k

        gamma_k_plus_1 = euclidian_gradient(w_k_plus_1)

        g_k_plus_1 = gamma_k_plus_1 * w_k_plus_1' - w_k_plus_1*gamma_k_plus_1'
        gamma_small_k = inner_product((g_k_plus_1 - g_k) , g_k_plus_1) / inner_product(g_k , g_k)

        h_k_plus_1 = g_k_plus_1 + gamma_small_k * h_k

        if inner_product(h_k_plus_1, g_k_plus_1) < 0.
            h_k_plus_1 = g_k_plus_1
        end

        k+=1
        gamma_k = gamma_k_plus_1
        w_k = w_k_plus_1
        h_k = h_k_plus_1

    end

    if optimization_success == false
        #println("optimisation failed")
    end

    return (w_k, optimization_success)

end

function run_brockett_optimization_test(;maxiter = 100)

    achieved_an_optimization = false
    iter = 0
    this_test = nothing

    while achieved_an_optimization == false && iter < maxiter

        iter += 1
        n = 3
        N = diagm(1:n)
        sigma = hermitian_test_matrix(n)
        q_brockett = 2

        function objective_function(w, sigma, N)
            tr(w' * sigma * w * N)
        end

        function euclidian_gradient_brockett(w,sigma,N)
            sigma*w*N
        end

        function euclidian_gradient(w)
            euclidian_gradient_brockett(w,sigma,N)
        end

        result= minimize_over_unitary_matrices(;euclidian_gradient = euclidian_gradient, q = q_brockett,n=n, tol = 1e-7, max_iter = 100000)

        achieved_an_optimization = result[2]

        if achieved_an_optimization
            optimized_matrix = result[1]
            this_test = @test real.(optimized_matrix' * sigma * optimized_matrix) ≈ diagm(sort(eigvals(sigma), rev = true)) atol = 1e-2

            break
        end

    end

    if achieved_an_optimization == false
        println("no optimization succeeded during the test, maybe be statistical but looks like an anomaly")
    end

    this_test

end


text:
## Diversity Measures

function shannon(v::AbstractVector{T}) where T<:Real
    total = sum(v)
    relab = map(x-> x/total, v)
    return -sum([log(x^x) for x in relab])
end

function shannon(abt::AbstractComMatrix)
    return shannon.(eachcol(occurrences(abt)))
end

function ginisimpson(v::AbstractVector{T}) where T<:Real
    total = sum(v)
    relab = map(x-> x/total, v)
    return 1 - sum([x^2 for x in relab])
end

function ginisimpson(abt::AbstractComMatrix)
    return ginisimpson.(eachcol(occurrences(abt)))
end


text:
using CallByName
VERSION < v"0.7-DEV" ? using Base.Test : using Test

a = 42
x = @stub a + 1 Int
@test apply(x) == 43
@test a == 42

function f(x::CallArg{T}, y) where T<:Number
    (apply(x), apply(x), y)
end

@test f(a, 1) == (42, 42, 1)

@test f((@stub a Int), 1) == (42, 42, 1)

@test f((@stub 13 Int), 1) == (13, 13, 1)

@test f((@stub global a = a + 1 Int), 1) == (43, 44, 1)

@test_throws UndefVarError f((@stub a = a + 1 Int), 2) == (43, 44, 2)

function g()
    a = 402 + 2.0im
    f((@stub a = a + 1 Complex{Int}), 10)
end

@test g() == (403 + 2im, 404 + 2im, 10)


text:
using CryptoTools
using Test
using SafeTestsets

@testset "CryptoTools.jl" begin
    @time @safetestset "Block Cipher" begin include("BlockCipher.jl") end
    @time @safetestset "Stream Cipher" begin include("StreamCipher.jl") end
    @time @safetestset "Secret Sharing" begin include("SecretSharing.jl") end
end


text:
"""
以下のプログラムを Julia に翻訳してみる。

レーダーチャート
http://aoki2.si.gunma-u.ac.jp/R/radar.html

使用例

using RDatasets
swiss = dataset("datasets", "swiss")
test = swiss[1:5, 2:7]
radar(test)
radar(test, col=[:black, :red, :green4, :blue, :brown])
radar(test, maxx=[100, 50, 20, 20, 100, 25], normalize=false,
        col=[:black, :red, :green4, :blue, :brown])
radar(test, maxx=[100, 50, 20, 20, 100, 25],
        minx=[70, 10, 0, 0, 0, 20], normalize=false,
        col=[:black, :red, :green4, :blue, :brown])
radar(test, maxx=100, minx=0, normalize=false,
        col=[:black, :red, :green4, :blue, :brown])

"""

using Statistics, Plots

function radar(df; maxx=Inf, minx=Inf, normalize=true, col="", lty="", title="")
    function drawnet(x, border, lty)
        scale = (x .- minx) ./ (maxx .- minx)
        scale = vcat(scale, scale[1])
        if all(0 .< scale .<= 1.2)
            plot!(vcat(cosine, cosine[1]).*scale, vcat(sine, sine[1]).*scale,
                  seriestype=:path, label="", linecolor=border, linestyle=lty)
        end
    end
    pyplot()
    axisname = names(df)
    df = Matrix(df)
    n, m = size(df)
    m >= 3 || error("3変数以上であるべし！")
    p = plot(xlimit=(-1.4, 1.4), ylimit=(-1.2, 1.2), aspect_ratio=1, axis=false,
             ticks=false, xlabel="", ylabel="", title=title, size=(600, 400))
    θ = π .*(0.5 .- collect(0:m-1) .* 2 ./ m)
    cosine = cos.(θ)
    sine = sin.(θ)
    if normalize
        means = mean(df, dims=1)
        stds = std(df, dims=1)
        for i = 1:m
            df[:, i] = (df[:, i] .- means[i]) ./ stds[i]
        end
        maxx = maximum(df)
        minx = minimum(df)
        w = (maxx - minx) * 0.1
        maxx = maxx + w
        minx = minx - w
        for i = -3:3
            x = fill(i, m)
            drawnet(x, :gray, i==0 ? :solid : :dot)
        end
    else
        maxx != Inf || (maxx = vec(maximum(df, dims=1)))
        minx != Inf || (minx = vec(minimum(df, dims=1)))
        w = (maxx - minx) * 0.1
        maxx = maxx + w
        minx = minx - w
        for i = 1:5
            x = cosine * i / 5
            y = sine * i / 5
            plot!(vcat(x, x[1]), vcat(y, y[1]), seriestype=:path,
                    linecolor=:gray, linestyle=:dot, label="")
        end
    end
    for i = 1:m
        plot!([0, cosine[i]*1.05], [0, sine[i]*1.05], seriestype=:path,
                color=:gray, label="")
        annotate!(cosine[i]*1.1, sine[i]*1.1, text(axisname[i], 10,
                i > m/2 ? :right : :left))
    end
    col != "" || (col = fill(:red, m))
    lty != "" || (lty = fill(:solid, m))
    for i = 1:n
        drawnet(df[i, :], col[i], lty[i])
    end
    display(p)
end

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
__precompile__()
module CasADi

using PyCall

import Base: convert, getproperty, hcat, length, promote_rule, show, size, vcat
import Base: +, -, *, /, \, ^
import Base: >, >=, <, <=, ==
import LinearAlgebra: ×

export CasadiSymbolicObject, SX, MX
export casadi, to_julia, substitute

include("types.jl")
include("constructors.jl")
include("numbers.jl")
include("mathops.jl")
include("mathfuns.jl")
include("generic.jl")
include("utils.jl")

##################################################

const casadi = PyNULL()

function __init__()
    ## Define casadi
    copy!(casadi, pyimport_conda("casadi", "casadi", "conda-forge"))

    pytype_mapping(casadi.SX, SX)
    pytype_mapping(casadi.MX, MX)
end

## Add generic and new methods
include("importexport.jl")

end # module


text:
using JuMP
using Gurobi
using DataFrames
using CSV

##### Problem 2 code skeleton

### Solution

# Note: @__DIR__ points to the directory which contains this script file
# Read data
cd(string(@__DIR__, "/Data"))
sessions = convert(Matrix,CSV.File("Pb2_sessions.csv"; header=true) |> DataFrame!);
courses = convert(Matrix,CSV.File("Pb2_courses.csv"; header=true) |> DataFrame!);
rooms = convert(Matrix,CSV.File("Pb2_rooms.csv"; header=true) |> DataFrame!);
preferences = convert(Matrix,CSV.File("Pb2_preferences.csv"; header=true) |> DataFrame!);
cohorts = convert(Matrix,CSV.File("Pb2_cohorts.csv"; header=true) |> DataFrame!);

cd(@__DIR__)

# Initialize sets and parameters
J = size(sessions,1); #number of sessions
C = size(courses,1); #number of courses
R = size(rooms,1); #number of classrooms
S = size(cohorts,1); #number of student cohorts
D = 2; #number of days
T = 6; #number of time blocks
Q = rooms[:,2]; #classroom capacities
U = courses[:,2]; #course units
N = cohorts[:,2]; #number of students per cohort
ms = cohorts[:,3]; #minimum number of units per cohort
Ms = cohorts[:,4]; #maximum number of units per cohort

mandatory = zeros(S,C); #mandatory courses per cohort
mandatory[1:6,1:5] .= 1;
mandatory[7:13,6:9] .= 1;
mandatory[14:15,10:11] .= 1;
mandatory[16:17,12:13] .= 1;

# Start building model
timelimit = 10*60;
model = Model(Gurobi.Optimizer)
set_optimizer_attribute(model, "TimeLimit", timelimit)

# VARIABLES
@variable(model,z[1:S,1:J,1:D,1:T],Bin);
@variable(model,x[1:J,1:D,1:T,1:R],Bin);

# CONSTRAINTS
@constraint(
    model, [s in 1:S, j in 1:J, d in 1:D, t in 1:T],
    z[s,j,d,t] <= sum(x[j,d,t,r] for r in 1:R)
)

@constraint(
    model,[s in 1:S, d in 1:D, t in 1:T],
    sum(z[s,j,d,t] for j in 1:J) <= 1
)

for c in 1:C
    @constraint(model,
        [s in 1:S],
        sum(z[s,j,d,t]
            for d in 1:D,
                t in 1:T,
                j in sessions[findall(x -> x==c,sessions[:,2]),1])
            <= 1);
end

for s in 1:S
    mandatory_s = findall(x -> x==1,mandatory[s,:]);
    for c in mandatory_s
        @constraint(model,
            sum(z[s,j,d,t]
            for d in 1:D,
                t in 1:T,
                j in sessions[findall(x -> x==c,sessions[:,2]),1])
            >= 1);
    end
end

@constraint(model,[r in 1:R, d in 1:D, t in 1:T],sum(x[j,d,t,r] for j in 1:J) <= 1);

@constraint(model,[j in 1:J],sum(x[j,d,t,r] for r in 1:R, d in 1:D, t in 1:T) == 1);

@constraint(model,[j in 1:J, d in 1:D, t in 1:T],
            sum(N[s]*z[s,j,d,t] for s in 1:S) <= sum(Q[r]*x[j,d,t,r] for r in 1:R))


text:
#using Distributions
#include("constants.jl")
#include("orbit.jl")
#include("planet.jl")

if !isdefined(:PlanetarySystemAbstract)
  @compat abstract type PlanetarySystemAbstract end

  immutable PlanetarySystem{StarT<:StarAbstract} <: PlanetarySystemAbstract
    star::StarT
    planet::Vector{Planet}         # TODO OPT: If want to optimize, consider using  something like ImmmutableArrays?
    orbit::Vector{Orbit}           # TODO OPT: If want to optimize, consider using  something like ImmmutableArrays?
    # TODO DETAIL: Setup inner constructor to enforce equal number of planets & orbits
    #PlanetarySystem{StarT}(s::StarAbstract, p::Vector{Planet}, o::Vector{Orbit})
    #  @assert(length(p)==length(o)) # else error(string("Number of planets must match number of orbits: Np= ",length(p)," No= ",length(o)))
    #  new(s,p,o)
    #end
  end
  #typealias PlanetarySystemSingleStar  PlanetarySystem{SingleStar}
  PlanetarySystemSingleStar = PlanetarySystem{SingleStar}

end
PlanetarySystem{StarT<:StarAbstract}(s::StarT) = PlanetarySystem(s,Array{Planet}(0),Array{Orbit}(0))  # Constructor for a Planetary System with no planets

PlanetarySystem{StarT<:StarAbstract}(s::StarT, p::Planet, o::Orbit) = PlanetarySystem(s,[p],[o])  # Constructor for a single Planet System

function PlanetarySystem{StarT}(ps::PlanetarySystem{StarT}, keep::Vector{Int64}) # Why doesn't this work?
   PlanetarySystem{StarT}(ps.star,ps.planet[keep],ps.orbit[keep])
end

#function PlanetarySystemSingleStar(ps::PlanetarySystemSingleStar, keep::Vector{Int64})
#function PlanetarySystem{StarT<:StarAbstract}(ps::PlanetarySystem{StarT}, keep::Vector{Int64})
#   PlanetarySystem(ps.star,ps.planet[keep],ps.orbit[keep])
#end

flux(ps::PlanetarySystem{StarAbstract}) = flux(ps.star)
flux(ps::PlanetarySystem{Star}) = flux(ps.star)
flux(ps::PlanetarySystem{BinaryStar}) = flux(ps.star)
flux(ps::PlanetarySystem{MultipleStar}) = flux(ps.star)

function num_planets(s::PlanetarySystem{Star})
  @assert( length(s.planet) == length(s.orbit) )    # TODO OPT: Deactivate inner assert's like this for speed once tested
  return length(s.planet)
end

function generate_planet_mass_from_radius_powerlaw(r::Float64, s::Star, o::Orbit, sim_param::SimParam) # TODO SCI: IMPORTANT once have stability criteria, replace w/ better M-R relationship
  const mr_power_index::Float64 = get_real(sim_param,"mr_power_index")
  const mr_const::Float64 = get_real(sim_param,"mr_const")
  m = mr_const*earth_mass*(r/earth_radius)^mr_power_index
  return m
end

function generate_num_planets_poisson(lambda::Real, max_planets_in_sys::Integer; min_planets_in_sys::Integer = 0)
  if lambda < min_planets_in_sys*1e-3
      return min_planets_in_sys
  end
  d = Distributions.Truncated(Distributions.Poisson(lambda),min_planets_in_sys,max_planets_in_sys)
  n = rand(d)
  #=
  n = -1
  while !(min_planets_in_sys<=n<=max_tranets_in_sys)
     n = rand(Distributions.Poisson(lambda))
  end
  =#
 return n
end

function generate_num_planets_poisson(s::Star, sim_param::SimParam)
  const lambda::Float64 = exp(get_real(sim_param,"log_eta_pl"))
  const max_tranets_in_sys::Int64 = get_int(sim_param,"max_tranets_in_sys")
  generate_num_planets_poisson(lambda,max_tranets_in_sys)
end

function generate_period_and_sizes_log_normal(s::Star, sim_param::SimParam; num_pl::Integer = 1)  # TODO SCI:  IMPORTANT: Improve how periods and sizes are drawn
    const mu_log_r::Float64 = get_real(sim_param,"mean_log_planet_radius")
    const sigma_log_r::Float64 = get_real(sim_param,"sigma_log_planet_radius")
    const mu_log_P::Float64 = get_real(sim_param,"mean_log_planet_period")
    const sigma_log_P::Float64 = get_real(sim_param,"sigma_log_planet_period")
    const min_period::Float64 = get_real(sim_param,"min_period")
    const max_period::Float64 = get_real(sim_param,"max_period")
    const min_radius::Float64 = get_real(sim_param,"min_radius")
    const max_radius::Float64 = get_real(sim_param,"max_radius")
    const max_draws::Int64 = 100

    if   sigma_log_r <= 0. || sigma_log_P<=0.
     println("# mu_log_r= ", mu_log_r, " sigma_log_r= ", sigma_log_r, " mu_log_P= ", mu_log_P, " sigma_log_P= ", sigma_log_P)
    end
    const rdist = LogNormal(mu_log_r,sigma_log_r)
    const Pdist = LogNormal(mu_log_P,sigma_log_P)
    #Rlist = rand(rdist,num_pl)
    #Plist = rand(Pdist,num_pl)
    #idx_keep = find(i->(min_radius<=Rlist[i]<=max_radius) && (min_period<=Plist[i]<=max_period), 1:num_pl )
    #return Plist[idx_keep], Rlist[idx_keep]  # replaced because want to return exactly num_pl
    Rlist = zeros(num_pl)
    Plist = zeros(num_pl)
    for i in 1:num_pl
      j = 0
      while ! (min_radius<Rlist[i]<=max_radius) && j<max_draws
            Rlist[i] = rand(rdist)
            j+=1
      end
      if j>=max_draws
         println("# Struggled to draw size for: ",mu_log_r, " ", sigma_log_r)
      end
      j = 0
      while ! (min_period<Plist[i]<=max_period) && j<max_draws
            Plist[i] = rand(Pdist)
            j+=1
      end
      if j>=max_draws
         println("# Struggled to draw period for: ",mu_log_P, " ", sigma_log_P)
      end
    end
    return Plist, Rlist
end

function draw_power_law(n::Real, x0::Real, x1::Real, num_pl::Integer)
     ((x1^(n+1) - x0^(n+1))*rand(num_pl) + x0^(n+1)).^(1/(n+1))
end

function generate_period_and_sizes_power_law(s::Star, sim_param::SimParam; num_pl::Integer = 1)
    const power_law_r::Float64 = get_real(sim_param,"power_law_r")
    const power_law_P::Float64 = get_real(sim_param,"power_law_P")
    const min_period::Float64 = get_real(sim_param,"min_period")
    const max_period::Float64 = get_real(sim_param,"max_period")
    const min_radius::Float64 = get_real(sim_param,"min_radius")
    const max_radius::Float64 = get_real(sim_param,"max_radius")
    Rlist = power_law_r!=-1.0 ? draw_power_law(power_law_r,min_radius,max_radius, num_pl) : exp(log(min_radius).+rand()*log(max_radius/min_radius))
    Plist = power_law_P!=-1.0 ? draw_power_law(power_law_P,min_period,max_period, num_pl) : exp(log(min_period).+rand()*log(max_period/min_period))
    return Plist, Rlist
end

function generate_e_omega_rayleigh(sigma_hk::Float64)
  h = k = 1.0
  while h*h+k*k >= 1.0
    h = sigma_hk*randn()
    k = sigma_hk*randn()
  end
  ecc = sqrt(h*h+k*k)
  w = atan2(k,h)
  return ecc, w
end

function generate_e_omega_rayleigh(sim_param::SimParam)
  sigma_hk::Float64 = get_real(sim_param,"sigma_hk")
  generate_e_omega_rayleigh(sigma_hk)
end

function generate_planetary_system_hardcoded_example(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # in this version we specify fixed functions that are known at compile time, allowing for additional optimizations (~0.6 second faster per Kepler catalog out of ~3.6 sec on my laptop w/ 1 core)
  const  generate_planet_mass_from_radius = generate_planet_mass_from_radius_powerlaw
  const  generate_num_planets = generate_num_planets_poisson
  const  generate_period_and_sizes = generate_period_and_sizes_log_normal

  const  generate_e_omega =  generate_e_omega_rayleigh

  # const  generate_star = get_function(sim_param,"generate_star")
  # const star::StarAbstract = generate_star(sim_param)
  const num_pl::Int64 = generate_num_planets(star, sim_param)

  if( num_pl==0 )
    return PlanetarySystem(star)
  else
     pl = Array{Planet}(num_pl)
     orbit = Array{Orbit}(num_pl)
     (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
     idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred

    for i in 1:num_pl
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      (ecc::Float64,  omega::Float64) = generate_e_omega(sim_param)
      incl::Float64 = acos(rand())
      orbit[idx[i]] = Orbit(P,ecc,incl,omega,2pi*rand(),2pi*rand())
      # set!(orbit[idx[i]],P,ecc,incl,omega,2pi*rand(),2pi*rand())
      mass::Float64 = generate_planet_mass_from_radius(Rlist[idx[i]], star, orbit[idx[i]], sim_param)
      pl[i] = Planet( Rlist[idx[i]],  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

function generate_planetary_system_uncorrelated_incl(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # load functions to use for drawing parameters
  const  generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
  const  generate_num_planets = get_function(sim_param,"generate_num_planets")
  const  generate_period_and_sizes = get_function(sim_param,"generate_period_and_sizes")
  const  generate_e_omega = get_function(sim_param,"generate_e_omega")

  # const  generate_star = get_function(sim_param,"generate_star")
  # const star::StarAbstract = generate_star(sim_param)
  const num_pl::Int64 = generate_num_planets(star, sim_param)::Int64
  const sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

  if( num_pl==0 )
    return PlanetarySystem(star)::PlanetarySystem
  else
     pl = Array{Planet}(num_pl)
     orbit = Array{Orbit}(num_pl)
     (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
     idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred

    for i in 1:num_pl
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      (ecc::Float64,  omega::Float64) = generate_e_omega(sim_param)
      incl::Float64 = acos(rand())
      orbit[idx[i]] = Orbit(P,ecc,incl,omega,2pi*rand(),2pi*rand())
      # set!(orbit[idx[i]],P,ecc,incl,omega,2pi*rand(),2pi*rand())
      mass::Float64 = generate_planet_mass_from_radius(Rlist[idx[i]], star, orbit[idx[i]], sim_param)
      pl[i] = Planet( Rlist[idx[i]],  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

# This version generates more systems roughly near a common plane, but until incorporate CORBITS data, ABC can't match input param
function generate_planetary_system_simple(star::StarAbstract, sim_param::SimParam; verbose::Bool = false)
  # load functions to use for drawing parameters
  const  generate_planet_mass_from_radius = get_function(sim_param,"generate_planet_mass_from_radius")
  const  generate_num_planets = get_function(sim_param,"generate_num_planets")
  #const  generate_num_planets = generate_num_planets_christiansen
  #const  generate_period_and_sizes = generate_period_and_sizes_christiansen
  const  generate_period_and_sizes = get_function(sim_param,"generate_period_and_sizes")
  const  generate_e_omega = get_function(sim_param,"generate_e_omega")
  const  sigma_incl = deg2rad(get_real(sim_param,"sigma_incl"))

  # const  generate_star = get_function(sim_param,"generate_star")
  # const star::StarAbstract = generate_star(sim_param)
  const num_pl = generate_num_planets(star, sim_param)::Int64
  const sigma_ecc::Float64 = haskey(sim_param,"sigma_hk") ? get_real(sim_param,"sigma_hk") : 0.0

  if( num_pl==0 )
    return PlanetarySystem(star)
  else
     pl = Array{Planet}(num_pl)
     orbit = Array{Orbit}(num_pl)
     (Plist::Vector{Float64}, Rlist::Vector{Float64}) = generate_period_and_sizes(star, sim_param, num_pl=num_pl)
     idx = sortperm(Plist)                   # TODO OPT: Check to see if sorting is significant time sink.  If so, it might could be deferred
     incl_sys = acos(rand())

    for i in 1:num_pl
      # if verbose   println("i=",i," idx=",idx," Plist=",Plist[idx] );     end
      P = Plist[idx[i]]
      if haskey(sim_param,"sigma_hk_one") && haskey(sim_param,"sigma_hk_multi")
         sigma_ecc = num_pl == 1 ? get_real(sim_param,"sigma_hk_one") : get_real(sim_param,"sigma_hk_multi")
      end
      (ecc,  omega) = generate_e_omega(sim_param)::Tuple{Float64,Float64}
      incl_mut = sigma_incl*sqrt(randn()^2+randn()^2) # rand(Distributions.Rayleigh(sigma_incl)) # sigma_incl*randn()
      asc_node = 2pi*rand()
      mean_anom = 2pi*rand()
      #incl = incl_sys + sigma_incl*randn()
      incl =  incl_mut!=zero(incl_mut) ? acos( cos(incl_sys)*cos(incl_mut) + sin(incl_sys)*sin(incl_mut)*cos(asc_node) ) : incl_sys
      orbit[idx[i]] = Orbit(P,ecc,incl,omega,asc_node,mean_anom)
      # set!(orbit[idx[i]], P,ecc,incl,omega,asc_node,mean_anom) # if Orbit were mutable
      mass = generate_planet_mass_from_radius(Rlist[idx[i]], star, orbit[idx[i]], sim_param)::Float64
      pl[i] = Planet( Rlist[idx[i]],  mass )
    end
  return PlanetarySystem(star,pl,orbit)
  end
end

function test_planetary_system_constructors(sim_param::SimParam)
  generate_star = get_function(sim_param,"generate_star")
  star = generate_star(sim_param)
  empty_sys = PlanetarySystem(star)
  earth = Planet(earth_radius,earth_mass)
  earth_orbit = Orbit(365.2425,0.0167,0.5*pi,0.0,0.0,0.0)
  solar_sys = PlanetarySystem(star, earth,earth_orbit)
  m = generate_planet_mass_from_radius_powerlaw(0.02,star,earth_orbit,sim_param)/earth_mass
  generate_planetary_system_simple(star,sim_param,verbose=true)
end


text:
function produce_mode_query(modes::Vector{String}) :: String

    mode_query = ""
    if length(modes) > 0
        mode_query = "AND ("
        for j in 1:length(modes)
            mode_query *= "mode like '%$(modes[j])%' "

            if j < length(modes)
                mode_query *= "OR "
            end

        end
        mode_query *= ")"
    end

    return mode_query

end

function extract_OCO_sounding_id(sounding_id::Int)
    return parse(Int, string(sounding_id)[end])
end

function convert_OCO_sounding_id_to_fp(sounding_id::Int)
    return(parse(Int, string(sounding_id)[end]))
end

function convert_OCO_sounding_id_to_date(sounding_id::Int)

    # Conert to String
    snid = string(sounding_id)

    # OCO-type sounding ID must be 16 digits long
    if length(snid) != 16
        @error "Sounding ID $(sonding_id) must be 16 characters/digits long!"
        return nothingx
    end

    id_year = parse(Int, snid[1:4])
    id_month = parse(Int, snid[5:6])
    id_day = parse(Int, snid[7:8])
    id_hour = parse(Int, snid[9:10])
    id_minute = parse(Int, snid[11:12])
    id_second = parse(Int, snid[13:14])
    id_milli = parse(Int, snid[15:16])

    return DateTime(id_year, id_month, id_day,
                    id_hour, id_minute, id_second, id_milli)
end

# These dictionaries contain the footprint-dependent coefficients
# to produce the SIF single-sounding uncertainty of the form
# sqrt(A + B * continuum_radiance)

OCO2_noise_dict = Dict(
    1 => (0.00243, 0.00192),
    2 => (0.00204, 0.00179),
    3 => (0.00365, 0.00174),
    4 => (0.00331, 0.0018) ,
    5 => (0.00329, 0.00188),
    6 => (0.00301, 0.00203),
    7 => (0.00241, 0.00211),
    8 => (0.00373, 0.00242)
)

OCO3_noise_dict = Dict(
    1 => (0.00016, 0.00161),
    2 => (0.00011, 0.00156),
    3 => (3.0e-5, 0.00144),
    4 => (3.0e-5, 0.00135),
    5 => (0.00011, 0.00131),
    6 => (0.00021, 0.0014),
    7 => (0.00037, 0.0015),
    8 => (0.00034, 0.00164)
)

function calculate_OCO_SIF_uncertainty(continuum::Real, instrument::String, fp::Int)
    if instrument == "OCO-3"
        noise_dict = OCO3_noise_dict
    elseif instrument == "OCO-2"
        noise_dict = OCO2_noise_dict
    else
        @error "Need to be either OCO-2 or OCO-3, but I got $(instrument)."
    end

    return sqrt(noise_dict[fp][1] + noise_dict[fp][2] * continuum)
end

function convert_OCO_df_to_scenes(df::DataFrame)

    # At this point, we expect all the instrument labels here
    # to be a unique

    unique_instrument = unique(df.instrument)

    if length(unique_instrument) != 1
        @error "We have multiple instruments via the SQLITE query: $(unique_instrument)"
        return nothing
    end

    instrument = unique_instrument[1]

    # Here we pre-load the surface data that we then sample
    # in the scene loop below. The strategy is clear: accessing
    # the surface data HDF file point-by-point is slow,
    # so we read in some intermediate objects here by only
    # loading a small subset that bounds the ROI. Once this is
    # in memory, we can very quickly sample and calculate reflectivity etc.

    lon_bound_min = minimum(df.lon)
    lon_bound_max = maximum(df.lon)
    lat_bound_min = minimum(df.lat)
    lat_bound_max = maximum(df.lat)

    vnp_sd = create_VNP_SD_from_locbounds(
        lon_bound_min,
        lat_bound_min,
        lon_bound_max,
        lat_bound_max
    )

    #=
    tropomi_data = create_TROPOMI_SIF_from_locbounds(
        lon_bound_min,
        lat_bound_min,
        lon_bound_max,
        lat_bound_max
    )
    =#

    locarray = Geolocation[]
    scenearray = Scene[]

    # Go through scenes and create scene and location objects
    for row in eachrow(df)

        this_footprint = convert_OCO_sounding_id_to_fp(row.sounding_id)

        this_loc = Geolocation(row.lon, row.lat)

        this_loctime = GeolocationTime(
            this_loc,
            convert_OCO_sounding_id_to_date(row.sounding_id)
        )

        # New unction for this particular footprint
        this_ucert_function(x) = calculate_OCO_SIF_uncertainty(
            x,
            instrument,
            this_footprint
        )

        this_scene = create_SIF_scene(
            instrument,
            row.mode,
            this_loctime,
            row.vza,
            this_ucert_function,
            vnp_sd
        )

        push!(locarray, this_loc)
        push!(scenearray, this_scene)
    end

    @info "Populated $(length(locarray)) scene objects."

    # We must time-order scenes (measurements)
    time_sort = sortperm(get_time.(scenearray))
    # Not sure what order we want for locations?

    return instrument, locarray[time_sort], scenearray[time_sort]

end

"""
    OCOSampling(target_lon, target_lat, radius, oco_location_file)

Constructor for an OCOSampling type object. Reads files from "datapath"
and subsets those scenes by only retain those locations which are
within the radius of the user-supplied target lon/lat location.

# Arguments
- 'target_lon::Number': Longitude of target location
- 'target_lat::Number': Latitude of target location
- 'radius::Number': Radius [km] within target location where scenes will be kept
- 'oco_location_file::String': Path to location of OCO-x locations file

# Returns
- 'OCOSampling': OCO sampling object

"""
function OCOSampling(target_lon::Real, target_lat::Real,
                     radius::Real, oco_db_file::String;
                     modes::Vector{String}=String[])

    # First - construct a sensible bounding box which will definitely contain the
    # locations in the circular search.

    wgs84 = Projection(Proj4.epsg[4326])
    xy = lonlat2xy([Float64(target_lon), Float64(target_lat)], wgs84)
    r = radius * 1000.0 * 1.02 # Add x percent as buffer, radius is given in [km] but we need [m]

    north = xy2lonlat(geod_destination(xy, 0.0, r, wgs84), wgs84)
    east = xy2lonlat(geod_destination(xy, 90.0, r, wgs84), wgs84)
    south = xy2lonlat(geod_destination(xy, 180.0, r, wgs84), wgs84)
    west = xy2lonlat(geod_destination(xy, 270.0, r, wgs84), wgs84)

    lon_min = minimum([west[1], south[1]])
    lat_min = minimum([west[2], south[2]])
    lon_max = maximum([east[1], north[1]])
    lat_max = maximum([east[2], north[2]])

    # Check for bounding box extent
    if (lon_min > lon_max)
        @error "Longitude maximum [$(lon_max)] is smaller than longitude minimum [$(lon_min)]"
        return nothing
    end

    if (lat_min > lat_max)
        @error "Latitude maximum [$(lon_max)] is smaller than latitude minimum [$(lat_min)]"
        return nothing
    end

    # If the user requests certain modes only, we inject the following additional query
    # into the main query.
    mode_query = produce_mode_query(modes)

    # Construct SQL query
    query = """
    SELECT sounding_id, vza, ST_X(geometry) AS lon, ST_Y(geometry) AS lat, mode
    FROM locations WHERE
    PtDistWithin (MakePoint($(target_lon), $(target_lat), 4326), geometry, $(radius * 1000.0)) = 1
    $(mode_query)
    AND rowid IN
    (
        SELECT rowid FROM SpatialIndex WHERE f_table_name = 'locations' AND
        search_frame = BuildMBR($(lon_min), $(lat_min), $(lon_max), $(lat_max))
    );
    """

    # Load up DB and execute query
    db = SQLite.DB(oco_db_file)
    SQLite.enable_load_extension(db, true)
    DBInterface.execute(db, "SELECT load_extension('/home/psomkuti/miniconda3/lib/mod_spatialite.so');")
    df = DBInterface.execute(db, query) |> DataFrame

    # Obtain instrument label from database and attach to dataframe
    this_instrument = SQLite.getvalue(DBInterface.execute(db, "SELECT * FROM instrument;"), 1, String)
    df[!, "instrument"] .= this_instrument

    # Turn DataFrame into location and scene arrays
    instrument, locarray, scenearray = convert_OCO_df_to_scenes(df)

    # Create OCOSampling object and return
    info = "OCO-type sampling with center point $(target_lon), $(target_lat) and radius $(radius) m."

    return OCOSampling(
        info,
        [instrument],
        locarray,
        scenearray
    )

end

"""
    OCOSampling(bbox, oco_location_file)

Constructor for an OCOSampling type object. Reads files from "datapath"
and subsets those scenes by only retain those locations which are
within the bounding box defined as [lon_min, lat_min, lon_max, lat_max]

# Arguments
- 'bbox::Array{Number, 1}': Bounding box [lon_min, lat_min, lon_max, lat_max]
- 'oco_db_file::String': Path to location of OCO-x locations SQLite file

# Returns
- 'OCOSampling': OCO sampling object

"""
function OCOSampling(bbox::Vector{<:Real}, oco_db_file::String)

    # Check for validity of bounding box length
    if length(bbox) != 4
        @error "Length of bounding box vector must be 4!"
        return nothing
    end

    lon_min, lat_min, lon_max, lat_max = bbox

    # Check for bounding box extent
    if (lon_min > lon_max)
        @error "Longitude maximum [$(lon_max)] is smaller than longitude minimum [$(lon_min)]"
        return nothing
    end

    if (lat_min > lat_max)
        @error "Latitude maximum [$(lon_max)] is smaller than latitude minimum [$(lat_min)]"
        return nothing
    end

    # Construct SQL query
    query = """
    SELECT sounding_id, vza, ST_X(geometry) AS lon, ST_Y(geometry) AS lat, mode
    FROM locations WHERE
    MBRWithin(geometry, BuildMBR($(lon_min), $(lat_min), $(lon_max), $(lat_max))) = 1
    AND rowid IN
    (
        SELECT rowid FROM SpatialIndex WHERE f_table_name = 'locations' AND
        search_frame = BuildMBR($(lon_min), $(lat_min), $(lon_max), $(lat_max))
    );
    """

    # Load up DB and execute query
    db = SQLite.DB(oco_db_file)
    SQLite.enable_load_extension(db, true)
    DBInterface.execute(db, "SELECT load_extension('/home/psomkuti/miniconda3/lib/mod_spatialite.so');")
    df = DBInterface.execute(db, query) |> DataFrame

    # Turn DataFrame into location and scene arrays

    instrument, locarray, scenearray = convert_OCO_df_to_scenes(df)

    # Create OCOSampling object and return
    info = "OCO-type sampling with bounding box $(bbox) with N=$(size(df)[1])."

    return OCOSampling(
        info,
        [instrument],
        locarray,
        scenearray
    )

end


text:
JSON.lower(ds::DimensionSpec) = non_nothing_dict(ds)
JSON.lower(ef::ExtractionFunction) = non_nothing_dict(ef)

"""
    DefaultDS(dimension::String;
        outputName=nothing, outputType=nothing)

The "default" type dimension spec.

outputName and outputType should be `String`s if provided.
"""
struct DefaultDS <: DimensionSpec
    type::String
    dimension::String
    outputName
    outputType
    function DefaultDS(dimension; outputName=nothing, outputType=nothing)
        nothing_or_type(outputName, String)
        outputType === nothing || nothing_or_type(outputType, String) && (outputType = uppercase(outputType)) ∈ ["STRING", "LONG", "FLOAT"] || error("Invalid outputType")
        new("default", dimension, outputName, outputType)
    end
end

"""
    ExtractionDS(dimension::String, extractionFn::ExtractionFunction;
        outputName=nothing, outputType=nothing)

Dimension spec to use with an extraction function.

outputName and outputType should be `String`s if provided.
"""
struct ExtractionDS <: DimensionSpec
    type::String
    dimension::String
    extractionFn::ExtractionFunction
    outputName
    outputType
    function ExtractionDS(dimension, extractionFn; outputName=nothing, outputType=nothing)
        nothing_or_type(outputName, String)
        outputType === nothing || nothing_or_type(outputType, String) && (outputType = uppercase(outputType)) ∈ ["STRING", "LONG", "FLOAT"] || error("Invalid outputType")
        new("extraction", dimension, extractionFn, outputName, outputType)
    end
end

"""
    ListFiltered(delegate::DimensionSpec, values::Vector{String}; isWhitelist=nothing)

Filter the `delegate` dimension spec with the `values`.

isWhitelist should be a Bool if provided.
"""
struct ListFiltered <: DimensionSpec
    type::String
    delegate::DimensionSpec
    values::Vector{String}
    isWhitelist
    function ListFiltered(delegate, values; isWhitelist=nothing)
        nothing_or_type(isWhitelist, Bool)
        new("listFiltered", delegate, values, isWhitelist)
    end
end

"""
    RegexFiltered(delegate::DimensionSpec, pattern::String)

Filter the `delegate` dimension spec with the regex `pattern`.
"""
struct RegexFiltered <: DimensionSpec
    type::String
    delegate::DimensionSpec
    pattern::String
    RegexFiltered(delegate, pattern) = new("regexFiltered", delegate, pattern)
end

"""
    PrefixFiltered(delegate::DimensionSpec, prefix::String)

Filter the `delegate` dimension spec with the `prefix`.
"""
struct PrefixFiltered <: DimensionSpec
    type::String
    delegate::DimensionSpec
    prefix::String
    PrefixFiltered(delegate, prefix) = new("prefixFiltered", delegate, prefix)
end

"""
    LookupDS(dimension::String, name::String; outputName=nothing)

A dimension spec to use with the `name` lookup.

outputName should be a String if provided.
"""
struct LookupDS <: DimensionSpec
    type::String
    dimension::String
    name::String
    outputName
    function LookupDS(dimension, name; outputName=nothing)
        nothing_or_type(outputName, String)
        new("lookup", dimension, name, outputName)
    end
end

"""
    Map(isOneToOne::Bool, dict::Dict)

The lookup map used in the map lookup dimension spec.
"""
struct Map
    isOneToOne::Bool
    dict::Dict
end

JSON.lower(m::Map) = Dict("type" => "map", "map" => m.dict, "isOneToOne" => m.isOneToOne)

"""
    MapLookupDS(dimension::String, lookup::Map;
        outputName=nothing, retainMissingValue=nothing, replaceMissingValueWith=nothing, optimize=nothing)

Map lookup dimension spec.

outputName and replaceMissingValueWith should be `String`s, retainMissingValue
and optimize should be `Bool`s if provided.
"""
struct MapLookupDS <: DimensionSpec
    type::String
    dimension::String
    lookup::Map
    outputName
    retainMissingValue
    replaceMissingValueWith
    optimize
    function MapLookupDS(dimension, lookup; outputName=nothing, retainMissingValue=nothing, replaceMissingValueWith=nothing, optimize=nothing)
        nothing_or_type(outputName, String)
        nothing_or_type(retainMissingValue, Bool)
        nothing_or_type(replaceMissingValueWith, String)
        nothing_or_type(optimize, Bool)
        retainMissingValue != true || !isa(replaceMissingValueWith, String) || error("Cannon specify replaceMissingValueWith when retainMissingValue == true")
        new("lookup", dimension, lookup, outputName, retainMissingValue, replaceMissingValueWith, optimize)
    end
end

"""
    RegexEF(expr::String; index=nothing, replaceMissingValue=nothing, replaceMissingValueWith=nothing)

Regex based extraction function.

index should be a non negative integer, replaceMissingValue a Bool and
replaceMissingValueWith a String if provided.
"""
struct RegexEF <: ExtractionFunction
    type::String
    expr::String
    index
    replaceMissingValue
    replaceMissingValueWith
    function RegexEF(expr; index=nothing, replaceMissingValue=nothing, replaceMissingValueWith=nothing)
        index === nothing || (isa(index, Integer) && index >= 0) || error("index must be a non-negative integer")
        nothing_or_type(replaceMissingValue, Bool)
        new("regex", expr, index, replaceMissingValue, replaceMissingValueWith)
    end
end

"""
    PartialEF(expr::String)

Regex based extraction function that return nulls for non-matches.
"""
struct PartialEF <: ExtractionFunction
    type::String
    expr::String
    PartialEF(expr) = new("partial", expr)
end

"""
    SearchQueryEF(query::SearchQuerySpec)

Extraction using a search query spec.
"""
struct SearchQueryEF <: ExtractionFunction
    type::String
    query::SearchQuerySpec
    SearchQueryEF(query) = new("searchQuery", query)
end

"""
    SubstringEF(index::UInt64; length=nothing)

Substring extraction based on index and length.

length should be a non negative integer if provided.
"""
struct SubstringEF <: ExtractionFunction
    type::String
    index::UInt64
    length
    function SubstringEF(index; length=nothing)
        length === nothing || (isa(length, Integer) && length >= 0) || error("length must be a non-negative integer")
        new("substring", index, length)
    end
end

"""
    StrlenEF()

Extract string length.
"""
struct StrlenEF <: ExtractionFunction
    type::String
    StrlenEF() = new("strlen")
end

"""
    TimeFormatEF(; format=nothing, timeZone=nothing, locale=nothing, granularity=nothing, asMillis=nothing)

Extract timestamps.

format, timeZone, locale should be `String`s, granularity should be a
Granularity and asMillis should be a Bool if provided.
"""
struct TimeFormatEF <: ExtractionFunction
    type::String
    format
    timeZone
    locale
    granularity
    asMillis
    function TimeFormatEF(; format=nothing, timeZone=nothing, locale=nothing, granularity=nothing, asMillis=nothing)
        nothing_or_type(format, String)
        nothing_or_type(timeZone, String)
        nothing_or_type(locale, String)
        nothing_or_type(granularity, Granularity)
        nothing_or_type(asMillis, Bool)
        new("timeFormat", format, timeZone, locale, granularity, asMillis)
    end
end

"""
    TimeParseEF(timeFormat::String, resultFormat::String; joda=nothing)

Parse and extract time.

joda should be a Bool if provided.
"""
struct TimeParseEF <: ExtractionFunction
    type::String
    timeFormat::String
    resultFormat::String
    joda
    function TimeParseEF(timeFormat, resultFormat; joda=nothing)
        nothing_or_type(joda, Bool)
        new("time", timeFormat, resultFormat, joda)
    end
end

"""
    JavaScriptEF(jsfunction::String; injective=nothing)

Extract using a custom JavaScript function.

injective should be a Bool if provided.
"""
struct JavaScriptEF <: ExtractionFunction
    type::String
    jsfunction::String
    injective
    function JavaScriptEF(jsfunction; injective=nothing)
        nothing_or_type(injective, Bool)
        new("javascript", jsfunction, injective)
    end
end

function JSON.lower(ef::JavaScriptEF)
    d = Dict()
    d["type"] = ef.type
    d["function"] = ef.jsfunction
    if !(ef.injective === nothing)
        d["injective"] = ef.injective
    end
    d
end

"""
    RegisteredLookupEF(lookup::String;
        retainMissingValue=nothing, replaceMissingValueWith=nothing, injective=nothing, optimize=nothing)

Extract using a registered lookup.

retainMissingValue, injective, optimize should be `Bool`s and
replaceMissingValueWith should be a String if provided.
"""
struct RegisteredLookupEF <: ExtractionFunction
    type::String
    lookup::String
    retainMissingValue
    replaceMissingValueWith
    injective
    optimize
    function RegisteredLookupEF(lookup; retainMissingValue=nothing, replaceMissingValueWith=nothing, injective=nothing, optimize=nothing)
        nothing_or_type(retainMissingValue, Bool)
        nothing_or_type(replaceMissingValueWith, String)
        nothing_or_type(injective, Bool)
        nothing_or_type(optimize, Bool)
        retainMissingValue != true || !isa(replaceMissingValueWith, String) || error("Cannon specify replaceMissingValueWith when retainMissingValue == true")
        new("registeredLookup", lookup, retainMissingValue, replaceMissingValueWith, injective, optimize)
    end
end

"""
    InlineLookupEF(lookup::Map;
        retainMissingValue=nothing, replaceMissingValueWith=nothing, injective=nothing, optimize=nothing)

Extract using an inline lookup.

retainMissingValue, injective, optimize should be `Bool`s and
replaceMissingValueWith should be a String if provided.
"""
struct InlineLookupEF <: ExtractionFunction
    type::String
    lookup::Map
    retainMissingValue
    replaceMissingValueWith
    injective
    optimize
    function InlineLookupEF(lookup; retainMissingValue=nothing, replaceMissingValueWith=nothing, injective=nothing, optimize=nothing)
        nothing_or_type(retainMissingValue, Bool)
        nothing_or_type(replaceMissingValueWith, String)
        nothing_or_type(injective, Bool)
        nothing_or_type(optimize, Bool)
        retainMissingValue != true || !isa(replaceMissingValueWith, String) || error("Cannon specify replaceMissingValueWith when retainMissingValue == true")
        new("lookup", lookup, retainMissingValue, replaceMissingValueWith, injective, optimize)
    end
end

"""
    CascadeEF(extractionFns::Vector{<:ExtractionFunction})

Apply multiple extraction functions in sequence.
"""
struct CascadeEF <: ExtractionFunction
    type::String
    extractionFns::Vector{<:ExtractionFunction}
    CascadeEF(extractionFns) = new("cascade", extractionFns)
end

"""
    StringFormatEF(format::String; nullHandling=nothing)

Extract strings using a format.

nullHandling should be a String if provided.
"""
struct StringFormatEF <: ExtractionFunction
    type::String
    format::String
    nullHandling
    function StringFormatEF(format; nullHandling=nothing)
        nullHandling === nothing || nullHandling ∈ ["nullString", "emptyString", "returnNull"] || error("Invalid nullHandling")
        new("stringFormat", format, nullHandling)
    end
end

"""
    UpperEF(; locale=nothing)

Extract strings as upper case.

locale should be a String if provided.
"""
struct UpperEF <: ExtractionFunction
    type::String
    locale
    function UpperEF(;locale=nothing)
        nothing_or_type(locale, String)
        new("upper", locale)
    end
end

"""
    LowerEF(; locale=nothing)

Extract strings as lower case.

locale should be a String if provided.
"""
struct LowerEF <: ExtractionFunction
    type::String
    locale
    function LowerEF(;locale=nothing)
        nothing_or_type(locale, String)
        new("lower", locale)
    end
end

"""
    BucketEF(; size=nothing, offset=nothing)

Bucket numeric values into `size` buckets starting at `offset`.

size and offset should be non negative integers if provided.
"""
struct BucketEF <: ExtractionFunction
    type::String
    size
    offset
    function BucketEF(;size=nothing, offset=nothing)
        size === nothing || (isa(size, Integer) && size >= 0) || error("size must be a non-negative integer")
        offset === nothing || (isa(offset, Integer) && offset >= 0) || error("offset must be a non-negative integer")
        new("bucket", size, offset)
    end
end