text:
module NvMon

using ..LIKWID:
    LibLikwid, gputopo_initialized, nvmon_initialized, init_topology_gpu, GroupInfoCompact, get_gpu_topology

function init(gpus::AbstractVector{Int32}=Int32[i-1 for i in 1:(get_gpu_topology().numDevices)])
    nvmon_initialized[] && finalize()

    if !gputopo_initialized[]
        init_topology_gpu() || error("Couldn't init gpu topology.")
    end

    ngpus = length(gpus)
    ret = LibLikwid.nvmon_init(ngpus, gpus)
    if ret == 0
        nvmon_initialized[] = true
        return true
    end
    return false
end

init(gpus::AbstractVector{<:Integer}) = init(convert(Vector{Int32}, gpus))

function finalize()
    LibLikwid.nvmon_finalize()
    nvmon_initialized[] = false
    return nothing
end

_check_groupid(gid) = 0 ≤ gid < get_number_of_groups()
_check_eventidx(gid, eidx) = 0 ≤ eidx < get_number_of_events(gid)
_check_metricidx(gid, eidx) = 0 ≤ eidx < get_number_of_metrics(gid)
_check_gpuid(gpuid) = 0 ≤ gpuid < get_number_of_gpus()

"""
Return the number of GPUs initialized in the nvmon module.
"""
get_number_of_gpus() = LibLikwid.nvmon_getNumberOfGPUs()

"""
Return the number of groups currently registered in the nvmon module.
"""
get_number_of_groups() = LibLikwid.nvmon_getNumberOfGroups()

"""
Return the number of events in the group.
"""
get_number_of_events(groupid::Integer) = LibLikwid.nvmon_getNumberOfEvents(groupid)

"""
Return the number of metrics in the group.
Always zero for custom event sets.
"""
get_number_of_metrics(groupid::Integer) = LibLikwid.nvmon_getNumberOfMetrics(groupid)

"""
Return a list of all available nvmon groups for the GPU identified by `gpuid`.

# Examples
```jldoctest
julia> LIKWID.NvMon.get_groups()
4-element Vector{LIKWID.GroupInfoCompact}:
 DATA => Load to store ratio
 FLOPS_SP => Single-precision floating point
 FLOPS_HP => Half-precision floating point
 FLOPS_DP => Double-precision floating point
```
"""
function get_groups(gpuid::Integer=0)
    if !gputopo_initialized[]
        init_topology_gpu() || error("Couldn't init gpu topology.")
    end
    if !nvmon_initialized[]
        init(Int32[gpuid]) || error("Couldn't init nvmon.")
    end
    _check_gpuid(gpuid) || return nothing

    # refs to char**
    groups_ref = Ref{Ptr{Ptr{Cchar}}}()
    shorts_ref = Ref{Ptr{Ptr{Cchar}}}()
    longs_ref = Ref{Ptr{Ptr{Cchar}}}()

    ret = LibLikwid.nvmon_getGroups(gpuid, groups_ref, shorts_ref, longs_ref)
    ret <= 0 && return nothing

    groups_vec = unsafe_wrap(Array, groups_ref[], ret)
    shorts_vec = unsafe_wrap(Array, shorts_ref[], ret)
    longs_vec = unsafe_wrap(Array, longs_ref[], ret)
    res = Vector{GroupInfoCompact}(undef, ret)
    for i in 1:ret
        res[i] = GroupInfoCompact(
            unsafe_string(groups_vec[i]),
            unsafe_string(shorts_vec[i]),
            unsafe_string(longs_vec[i]),
        )
    end
    LibLikwid.nvmon_returnGroups(ret, groups_ref[], shorts_ref[], longs_ref[])
    return res
end

"""
    add_event_set(estr) -> groupid
Add a performance group or a custom event set to the nvmon module.
Returns a `groupid` which is required to later specify the event set.
"""
function add_event_set(estr::AbstractString)
    nvmon_initialized[] || return nothing
    groupid = LibLikwid.nvmon_addEventSet(estr)
    return Int(groupid)
end

"""
Return the name of the group identified by `groupid`.
If it is a custom event set, the name is set to `Custom`.
"""
function get_name_of_group(groupid::Integer)
    nvmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    name = unsafe_string(LibLikwid.nvmon_getGroupName(groupid))
    return name
end

"""
Return the short information about a performance group.
"""
function get_shortinfo_of_group(groupid::Integer)
    nvmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    sinfo = unsafe_string(LibLikwid.nvmon_getGroupInfoShort(groupid))
    return sinfo
end

"""
Return the (long) description of a performance group.
"""
function get_longinfo_of_group(groupid::Integer)
    nvmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    linfo = unsafe_string(LibLikwid.nvmon_getGroupInfoLong(groupid))
    return linfo
end

"""
Return the name of the event identified by `groupid` and `eventidx`.
"""
function get_name_of_event(groupid::Integer, eventidx::Integer)
    if !nvmon_initialized[] ||
       !_check_groupid(groupid) ||
       !_check_eventidx(groupid, eventidx)
        return nothing
    end
    name = unsafe_string(LibLikwid.nvmon_getEventName(groupid, eventidx))
    return name
end

"""
Return the name of the counter register identified by `groupid` and `eventidx`.
"""
function get_name_of_counter(groupid::Integer, eventidx::Integer)
    if !nvmon_initialized[] ||
       !_check_groupid(groupid) ||
       !_check_eventidx(groupid, eventidx)
        return nothing
    end
    name = unsafe_string(LibLikwid.nvmon_getCounterName(groupid, eventidx))
    return name
end

"""
Return the name of a derived metric identified by `groupid` and `metricidx`.
"""
function get_name_of_metric(groupid::Integer, metricidx::Integer)
    if !nvmon_initialized[] ||
       !_check_groupid(groupid) ||
       !_check_metricidx(groupid, metricidx)
        return nothing
    end
    name = unsafe_string(LibLikwid.nvmon_getMetricName(groupid, metricidx))
    return name
end

"""
Program the counter registers to measure all events in group `groupid`. Returns `true` on success.
"""
function setup_counters(groupid::Integer)
    nvmon_initialized[] || return nothing
    _check_groupid(groupid) || return false
    ret = LibLikwid.nvmon_setupCounters(groupid)
    return ret == 0
end

"""
Start the counter registers. Returns `true` on success.
"""
function start_counters()
    nvmon_initialized[] || return nothing
    ret = LibLikwid.nvmon_startCounters()
    return ret == 0
end

"""
Stop the counter registers. Returns `true` on success.
"""
function stop_counters()
    nvmon_initialized[] || return nothing
    ret = LibLikwid.nvmon_stopCounters()
    return ret == 0
end

"""
Read the counter registers.
To be executed after `start_counters` and before `stop_counters`.
Returns `true` on success.
"""
function read_counters()
    nvmon_initialized[] || return nothing
    ret = LibLikwid.nvmon_readCounters()
    return ret == 0
end

"""
Return the `groupid` of the currently activate group.
"""
function get_id_of_active_group()
    nvmon_initialized[] || return nothing
    return LibLikwid.nvmon_getIdOfActiveGroup()
end

"""
Switch currently active group to `groupid`. Returns `true` on success.
"""
function switch_group(groupid::Integer)
    nvmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    groupid == get_id_of_active_group() && return true
    ret = LibLikwid.nvmon_switchActiveGroup(groupid)
    return ret == 0
end

"""
Return the raw counter register result of the last measurement cycle identified by group `groupid` and the indices for event `eventidx` and thread `threadidx`.
"""
function get_last_result(groupid::Integer, eventidx::Integer, threadidx::Integer)
    nvmon_initialized[] || return nothing
    _check_eventidx(groupid, eventidx) || return nothing
    _check_gpuid(threadidx) || return nothing
    res = LibLikwid.nvmon_getLastResult(groupid, eventidx, threadidx)
    return res
end

"""
Return the raw counter register result of all measurements identified by group `groupid` and the indices for event `eventidx` and thread `threadidx`.
"""
function get_result(groupid::Integer, eventidx::Integer, threadidx::Integer)
    nvmon_initialized[] || return nothing
    _check_eventidx(groupid, eventidx) || return nothing
    _check_gpuid(threadidx) || return nothing
    res = LibLikwid.nvmon_getResult(groupid, eventidx, threadidx)
    return res
end

"""
Return the derived metric result of all measurements identified by group `groupid` and the indices for metric `metricidx` and thread `threadidx`.
"""
function get_metric(groupid::Integer, metricidx::Integer, threadidx::Integer)
    nvmon_initialized[] || return nothing
    _check_metricidx(groupid, metricidx) || return nothing
    _check_gpuid(threadidx) || return nothing
    res = LibLikwid.nvmon_getMetric(groupid, metricidx, threadidx)
    return res
end

"""
Return the derived metric result of the last measurement cycle identified by group `groupid` and the indices for metric `metricidx` and thread `threadidx`.
"""
function get_last_metric(groupid::Integer, metricidx::Integer, threadidx::Integer)
    nvmon_initialized[] || return nothing
    _check_metricidx(groupid, metricidx) || return nothing
    _check_gpuid(threadidx) || return nothing
    res = LibLikwid.nvmon_getLastMetric(groupid, metricidx, threadidx)
    return res
end

"""
Return the measurement time for group identified by `groupid`.
"""
function get_time_of_group(groupid::Integer)
    nvmon_initialized[] || return nothing
    _check_groupid(groupid) || return nothing
    time = LibLikwid.nvmon_getTimeOfGroup(groupid)
    return time
end

end # module


text:
using RunMetrics

tcx_dir = joinpath(@__DIR__, "..", "sample_data", "multi")

files = readdir(tcx_dir)
run_sums = Vector{RunSummary}(undef, 0)

for file in files
    if endswith(file, ".tcx")
        path = joinpath(tcx_dir, file)
        if isrun(path)
            rs = parse_run(path)
            push!(run_sums, rs)

        end
    end
end

β = RunMetrics.multi_regress(run_sums, 116.0)
# β_speed = β[1]
# β_alt = β[2]
# β_time = β[3:end]
# dates = [x.start_time for x in run_sums]

ext: jl


text:
abstract type Policy end

struct LinearFeedback <: Policy
	input
	output
	idx_input
	idx_input_nom
	idx_output
	p
end

function linear_feedback(input, output;
		idx_input = (1:input),
		idx_input_nom = (1:input),
		idx_output = (1:output))

	return LinearFeedback(input, output,
		idx_input, idx_input_nom, idx_output,
		input * output)
end

# linear state-feedback policy
function eval_policy(policy::LinearFeedback, θ, x, x̄, ū)
	view(ū, policy.idx_output) - reshape(θ, policy.output, policy.input) * (view(x, policy.idx_input) - view(x̄, policy.idx_input_nom))
end

struct PolicyConstraint <: Constraints
	n
	ineq
	control_dim
end

function policy_constraint(model, T; control_dim = model.m)
	n = control_dim * (T - 1)
	ineq = (1:0)

	return PolicyConstraint(n, ineq, control_dim)
end

struct PolicyConstraints <: Constraints
	con
	n
	ineq
end

function policy_constraints(prob, N;
		control_dim = prob.mean.model.m)

	T = prob.nom.T

	con = [policy_constraint(prob.mean.model, T,
		control_dim = control_dim)]
	n = con[end].n

	for i = 1:N
		push!(con, policy_constraint(prob.sample[i].model, T,
			control_dim = control_dim))
		n += con[i + 1].n
	end

	ineq = (1:0)

	return PolicyConstraints(con, n, ineq)
end

function constraints!(c, Θ, τ_nom, τ_sample, con::PolicyConstraint,
		policy::Policy,
		model_nom, model_sample,
		idx_τ_nom, idx_τ_sample, idx_θ,
		h, T)

	p = con.control_dim

	for t = 1:T-1
		θ = view(Θ, idx_θ[t])

		x̄ = view(τ_nom, idx_τ_nom.x[t])
		ū = view(τ_nom, idx_τ_nom.u[t])

		x = view(τ_sample, idx_τ_sample.x[t])
		u = view(τ_sample, idx_τ_sample.u[t])

		c[(t - 1) * p .+ (1:p)] = eval_policy(policy, θ, x, x̄, ū) - view(u, 1:p)
	end
    return nothing
end

function constraints_jacobian!(∇c, Θ, τ_nom, τ_sample, con::PolicyConstraint,
		policy::Policy,
		model_nom, model_sample,
		idx_τ_nom, idx_τ_sample, idx_θ,
		h, T)

	p = con.control_dim

	shift = 0

	for t = 1:T-1
		θ = view(Θ, idx_θ[t])

		x̄ = view(τ_nom, idx_τ_nom.x[t])
		ū = view(τ_nom, idx_τ_nom.u[t])

		x = view(τ_sample, idx_τ_sample.x[t])
		u = view(τ_sample, idx_τ_sample.u[t])

		pθ(y) = eval_policy(policy, y, x, x̄, ū) - view(u, 1:p)
		px(y) = eval_policy(policy, θ, y, x̄, ū) - view(u, 1:p)
		px̄(y) = eval_policy(policy, θ, x, y, ū) - view(u, 1:p)
		pū(y) = eval_policy(policy, θ, x, x̄, y) - view(u, 1:p)
		pu(y) = eval_policy(policy, θ, x, x̄, ū) - view(y, 1:p)

		# c[(t - 1) * p .+ (1:p)] = policy(model_sample, θ, x, x̄, ū) - view(u, 1:p)

		r_idx = (t - 1) * p .+ (1:p)
		c_idx = idx_θ[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(pθ, θ))
		shift += len

		c_idx = idx_τ_sample.x[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(px, x))
		shift += len

		c_idx = idx_τ_nom.x[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(px̄, x̄))
		shift += len

		c_idx = idx_τ_nom.u[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(pū, ū))
		shift += len

		c_idx = idx_τ_sample.u[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(pu, u))
		shift += len
	end

	return nothing
end

function constraints_sparsity(con::PolicyConstraint,
		policy,
		model_nom, model_sample,
		idx_τ_nom, idx_τ_sample, idx_θ,
		idx_z_nom, idx_z_sample, idx_z_policy,
		T;
		shift_row = 0, shift_col = 0)

	row = []
    col = []

	p = con.control_dim

	for t = 1:T-1
		r_idx = shift_row + (t - 1) * p .+ (1:p)

		c_idx = shift_col .+ idx_z_policy[idx_θ[t]]
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ idx_z_sample[idx_τ_sample.x[t]]
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ idx_z_nom[idx_τ_nom.x[t]]
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ idx_z_nom[idx_τ_nom.u[t]]
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ idx_z_sample[idx_τ_sample.u[t]]
		row_col!(row, col, r_idx, c_idx)
	end

    return collect(zip(row, col))
end

function constraints!(c, Z, con::PolicyConstraints,
	policy::Policy,
	prob::DPOProblems, idx::DPOIndices, N)

	shift = 0

	# mean
	Θ = view(Z, idx.policy)
	τ_nom = view(Z, idx.nom)
	τ_mean = view(Z, idx.mean)

	constraints!(view(c, shift .+ (1:con.con[1].n)),
			Θ, τ_nom, τ_mean, con.con[1], policy,
			prob.nom.model, prob.mean.model,
			prob.nom.idx, prob.mean.idx, idx.θ,
			prob.nom.h, prob.nom.T)

	shift += con.con[1].n

	for i = 1:N
		τ_sample = view(Z, idx.sample[i])
		constraints!(view(c, shift .+ (1:con.con[i + 1].n)),
				Θ, τ_nom, τ_sample,	con.con[i + 1], policy,
				prob.nom.model, prob.sample[i].model,
				prob.nom.idx, prob.sample[i].idx, idx.θ,
				prob.nom.h, prob.nom.T)
		shift += con.con[i + 1].n
	end
	nothing
end

function constraints_jacobian!(∇c, Z, con::PolicyConstraints,
		policy::Policy,
		prob::DPOProblems, idx::DPOIndices, N)

	shift = 0

	Θ = view(Z, idx.policy)
	τ_nom = view(Z, idx.nom)
	τ_mean = view(Z, idx.mean)

	len = length(constraints_sparsity(con.con[1],
			policy,
	 		prob.nom.model, prob.mean.model,
			prob.nom.idx, prob.mean.idx, idx.θ,
			idx.nom, idx.mean, idx.policy,
			prob.nom.T))

	constraints_jacobian!(view(∇c, shift .+ (1:len)),
			Θ, τ_nom, τ_mean, con.con[1], policy,
			prob.nom.model, prob.mean.model,
			prob.nom.idx, prob.mean.idx, idx.θ,
			prob.nom.h, prob.nom.T)

	shift += len

	for i = 1:N
		τ_sample = view(Z, idx.sample[i])

		len = length(constraints_sparsity(con.con[i + 1],
				policy,
		 		prob.nom.model, prob.sample[i].model,
				prob.nom.idx, prob.sample[i].idx, idx.θ,
				idx.nom, idx.mean, idx.policy,
				prob.nom.T))

		constraints_jacobian!(view(∇c, shift .+ (1:len)),
				Θ, τ_nom, τ_sample, con.con[i + 1], policy,
				prob.nom.model, prob.sample[i].model,
				prob.nom.idx, prob.sample[i].idx, idx.θ,
				prob.nom.h, prob.nom.T)

		shift += len
	end
	nothing
end

function constraints_sparsity(con::PolicyConstraints,
		policy::Policy,
		prob::DPOProblems, idx::DPOIndices, N;
		shift_row = 0, shift_col = 0)

	spar = constraints_sparsity(con.con[1],
			policy,
			prob.nom.model, prob.mean.model,
			prob.nom.idx, prob.mean.idx, idx.θ,
			idx.nom, idx.mean, idx.policy,
			prob.nom.T,
			shift_row = shift_row, shift_col = shift_col)
	shift_row += con.con[1].n

	for i = 1:N
		spar = collect([spar...,
				constraints_sparsity(con.con[i + 1],
				policy,
				prob.nom.model, prob.sample[i].model,
				prob.nom.idx, prob.sample[i].idx, idx.θ,
				idx.nom, idx.sample[i], idx.policy,
				prob.nom.T,
				shift_row = shift_row, shift_col = shift_col)...])
		shift_row += con.con[i + 1].n
	end

	return spar
end

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
#######################################
#
# Plot spectrum of signal
#
#######################################
"""
    plot_spectrum(signal::Vector, fs::Real; kwargs...)

Plot the spectrum of a signal.
"""
function plot_spectrum(
    signal::Vector,
    fs::Real;
    titletext::S = "",
    Fmin::Number = 0,
    Fmax::Number = 90,
    targetFreq::F = 0.0,
    dBPlot::Bool = true,
    noise_level::Number = 0,
    signal_level::Number = 0,
) where {S<:AbstractString,F<:AbstractFloat}

    # Determine fft frequencies
    signal_length = length(signal)
    frequencies = range(0, stop = 1, length = Int(signal_length / 2 + 1)) * fs / 2

    # Calculate fft and convert to power
    fftSweep = 2 / signal_length * fft(signal)
    spectrum = abs.(fftSweep[1:div(signal_length, 2)+1])  # Amplitude
    spectrum = spectrum .^ 2

    valid_idx = (frequencies .<= Fmax) .& (frequencies .>= Fmin)
    spectrum = spectrum[valid_idx]
    frequencies = frequencies[valid_idx]

    # Want a log plot?
    if dBPlot
        spectrum = 10 * log10.(spectrum)
        ylabel = "Response Power (dB)"
    else
        ylabel = "Response Power (uV^2)"
    end

    # Plot signal
    Plots.plot(frequencies, spectrum, lab = "Spectrum", color = :black)
    xlims!(Fmin, Fmax)
    xlabel!("Frequency (Hz)")
    ylabel!(ylabel)
    p = title!(titletext)

    # Plot the noise level if requested
    if noise_level != 0
        if dBPlot
            noise_level = 10 * log10.(noise_level)
        end
        p = plot!(
            [Fmin, targetFreq + 2],
            [noise_level, noise_level],
            lab = "Noise",
            color = :red,
        )
    end

    # Plot the signal level if requested
    if signal_level != 0
        if dBPlot
            signal_level = 10 * log10.(signal_level)
        end
        plot!(
            [Fmin, targetFreq],
            [signal_level, signal_level],
            lab = "Signal",
            color = :green,
        )

        targetFreqIdx = something(
            findfirst(
                isequal(minimum(abs.(frequencies .- targetFreq))),
                abs.(frequencies .- targetFreq),
            ),
            0,
        )
        targetFreq = frequencies[targetFreqIdx]
        targetResults = spectrum[targetFreqIdx]
        #TODO remove label for circle rather than make it empty
        p = plot!(
            [targetFreq],
            [targetResults],
            marker = (:circle, 5, 0.1, :green),
            markerstrokecolor = :green,
            lab = "",
        )
    end
    return p
end

#######################################
#
# Plot time series
#
# Automatic plotting for single channel (with units) or multichannel (with channel names) time series
#
#######################################

"""
    plot_single_channel_timeseries(signal::AbstractVector{T}, fs::Real; kwargs...)

Plot a single channel time series

# Input

* signal: Vector of data
* fs: Sample rate
* channels: Name of channel to plot
* plot_points: Number of points to plot, they will be equally spread. Used to speed up plotting
* Other optional arguements are passed to gadfly plot function

# Output

Returns a figure
"""
function plot_single_channel_timeseries(
    signal::AbstractVector{T},
    fs::Real;
    xlabel::S = "Time (s)",
    ylabel::S = "Amplitude (uV)",
    lab::S = "",
    kwargs...,
) where {T<:Number,S<:AbstractString}

    @debug("Plotting single channel waveform of size $(size(signal))")

    time_s = collect(1:size(signal, 1)) / fs   # Create time axis

    Plots.plot(
        time_s,
        signal,
        t = :line,
        c = :black,
        lab = lab,
        xlabel = xlabel,
        ylabel = ylabel,
    )
end

"""
    plot_multi_channel_timeseries(signals::Array{T,2}, fs::Number, channels::Array{S}; kwargs...)

Plot a multi channel time series

# Input

* signals: Array of data
* fs: Sample rate
* channels: Name of channels
* plot_points: Number of points to plot, they will be equally spread. Used to speed up plotting
* Other optional arguements are passed to gadfly plot function

#### Output

Returns a figure
"""
function plot_multi_channel_timeseries(
    signals::Array{T,2},
    fs::Number,
    channels::Array{S};
    xlabel::S = "Time (s)",
    ylabel::S = "Amplitude (uV)",
    kwargs...,
) where {T<:Number,S<:AbstractString}

    @debug("Plotting multi channel waveform of size $(size(signals))")

    time_s = collect(1:size(signals, 1)) / fs                        # Create time axis

    variances = var(signals, dims = 1)          # Variance of each figure for rescaling
    mean_variance = Statistics.mean(variances)     # Calculate for rescaling figures

    p = Plots.plot(
        t = :line,
        c = :black,
        xlabel = xlabel,
        ylabel = ylabel,
        ylim = (-0.5, size(signals, 2) - 0.5),
    )

    for c = 1:size(signals, 2)                                  # Plot each channel
        signals[:, c] = signals[:, c] .- Statistics.mean(signals[:, c])      # Remove mean
        signals[:, c] = signals[:, c] ./ (mean_variance ./ 4) .+ (c - 1)     # Rescale and shift so all chanels are visible
        p = plot!(time_s, signals[:, c], c = :black, lab = "")
    end
    p = plot!(yticks = (0:length(channels)-1, channels))

    return p
end

#######################################
#
# Filter response
#
#######################################

# Plot filter response
function plot_filter_response(
    zpk_filter::FilterCoefficients,
    fs::Integer;
    lower::Number = 1,
    upper::Number = 30,
    sample_points::Int = 1024,
)

    frequencies = range(lower, stop = upper, length = 1024)
    h = freqresp(zpk_filter, frequencies * ((2pi) / fs))
    magnitude_dB = 20 * log10.(convert(Array{Float64}, abs.(h)))
    phase_response = (360 / (2 * pi)) * unwrap(convert(Array{Float64}, angle.(h)))

    p1 = plot(frequencies, magnitude_dB, lab = "")
    p2 = plot(frequencies, phase_response, lab = "")

    p = plot(
        p1,
        p2,
        ylabel = ["Magnitude (dB)" "Phase (degrees)"],
        xlabel = "Frequency (Hz)",
        layout = @layout([a; b]),
    )
end


text:
module CUDNN

export
    conv2d!,
    conv2d,
    conv2d_grad_x!,
    conv2d_grad_x,
    conv2d_grad_w!,
    conv2d_grad_w,
    deconv2d,
    deconv2d_grad_x,
    deconv2d_grad_w,
    pool,
    pool_grad,
    unpool,
    unpool_grad,
    batchnorm_train!,
    batchnorm_infer!,
    batchnorm_grad!,
    batchnorm_train,
    batchnorm_infer,
    batchnorm_grad,
    BatchNormState,
    softmax4d!,
    softmax4d,
    softmax4d_grad!,
    softmax4d_grad

include("core.jl")

end # module


text:
using CollisionDetection
using CompScienceMeshes
using JLD

m = meshsphere(1.0, 0.2)
@show numcells(m)

centroid(ch) = cartesian(meshpoint(ch, [1,1]/3))

ctrs = [centroid(simplex(cellvertices(m,i))) for i in 1:numcells(m)]
rads = [maximum(norm(v-ctrs[i]) for v in cellvertices(m,i)) for i in eachindex(ctrs)]
tree = Octree(ctrs, rads)

fn = joinpath(@__FILE__,"..","center_sizes.jld")
save(fn,
    "ctrs", ctrs,
    "rads", rads)

max_stars_repo_name: shakisparki/CollisionDetection

max_issues_repo_name: shakisparki/CollisionDetection

max_forks_repo_name: shakisparki/CollisionDetection


text:
#===== Zygote =====#

if VERSION < v"1.6.9" # Zygote is failing on nightly

t5 = time()
using Zygote
# patch for https://github.com/FluxML/Zygote.jl/issues/897
@eval Zygote begin
   function _pullback(cx::AContext, ::typeof(sum), f, xs::AbstractArray)
      y, back = pullback(cx, ((f, xs) -> sum(f.(xs))), f, xs)
      y, ȳ -> (nothing, back(ȳ)...)
   end
end

GRAD = :Zygote
_gradient(x...) = Zygote.gradient(x...)

@tullio grad=Base
@testset "gradients: Zygote + DiffRules" begin include("gradients.jl") end

@tullio grad=Dual
@testset "gradients: Zygote + ForwardDiff" begin include("gradients.jl") end

@tullio grad=Base
@testset "complex gradients with Zygote" begin

    x0 = [1,2,3] .+ [5im, 0, -11im]
    # y0 = rand(Int8,3) .+ im .* rand(Int8,3) .+ 0.0
    @testset "analytic" begin

        g1 = _gradient(x -> real(sum(x)), x0)[1]
        g1i = _gradient(x -> imag(sum(x)), x0)[1]
        @test g1 ≈ _gradient(x -> real(@tullio y := x[i]), x0)[1]
        @test g1i ≈ _gradient(x -> imag(@tullio y := x[i]), x0)[1]

        g2 = _gradient(x -> real(sum(exp, x)), x0)[1]
        g2i = _gradient(x -> imag(sum(exp, x)), x0)[1]
        @test g2 ≈ _gradient(x -> real(@tullio y := exp(x[i])), x0)[1]
        @test g2i ≈ _gradient(x -> imag(@tullio y := exp(x[i])), x0)[1]

        g3 = _gradient(x -> real(sum(1 ./ (x.+im).^2)), x0)[1]
        g3i = _gradient(x -> imag(sum(1 ./ (x.+im).^2)), x0)[1]
        @test g3 ≈ _gradient(x -> real(@tullio y := 1/(x[i] + im)^2), x0)[1]
        @test g3 ≈ _gradient(x -> real(@tullio y := inv(x[i] + im)^2), x0)[1]
        @test g3i ≈ _gradient(x -> imag(@tullio y := 1/(x[i] + im)^2), x0)[1]
        @test g3i ≈ _gradient(x -> imag(@tullio y := inv(x[i] + im)^2), x0)[1]

        # with finaliser
        g7 = _gradient(x -> real(sum(sqrt.(sum(exp.(x), dims=2)))), x0 .+ x0')[1]
        g7i = _gradient(x -> imag(sum(sqrt.(sum(exp.(x), dims=2)))), x0 .+ x0')[1]
        @test g7 ≈ _gradient(x -> real(sum(@tullio y[i] := sqrt <| exp(x[i,j]) )), x0 .+ x0')[1]
        @test g7i ≈ _gradient(x -> imag(sum(@tullio y[i] := sqrt <| exp(x[i,j]) )), x0 .+ x0')[1]

    end
    @testset "non-analytic" begin

        g4 = _gradient(x -> real(sum(x * x')), x0)[1]
        g4i = _gradient(x -> imag(sum(x * x')), x0)[1] # zero!
        @test_broken g4 ≈ _gradient(x -> real(@tullio y := x[i] * conj(x[j])), x0)[1]
        @test_broken g4i ≈ _gradient(x -> imag(@tullio y := x[i] * conj(x[j])), x0)[1]
        @test_broken g4 ≈ _gradient(x -> real(@tullio y := x[i] * adjoint(x[j])), x0)[1]
        @test_broken g4i ≈ _gradient(x -> imag(@tullio y := x[i] * adjoint(x[j])), x0)[1]

        g5 = _gradient(x -> real(sum(abs2.(x .+ 2 .+ im))), x0)[1]
        g5i = _gradient(x -> imag(sum(abs2.(x .+ 2 .+ im))), x0)[1] # zero!
        @test_broken g5 ≈ _gradient(x -> real(@tullio y := abs2(x[i] + 2 + im)), x0)[1]
        @test_broken g5i ≈ _gradient(x -> real(@tullio y := abs2(x[i] + 2 + im)), x0)[1]

        g6 = _gradient(x -> real(sum(abs.(x.^3))), x0)[1]
        g6i = _gradient(x -> imag(sum(abs.(x.^3))), x0)[1] # zero!
        @test_broken g6 ≈ _gradient(x -> real(@tullio y := abs(x[i]^3)), x0)[1]
        @test_broken g6i ≈ _gradient(x -> real(@tullio y := abs(x[i]^3)), x0)[1]

    end
end

@info @sprintf("Zygote tests took %.1f seconds", time()-t5)

end # VERSION

#===== ReverseDiff =====#
#=
t6 = time()
using ReverseDiff

GRAD = :ReverseDiff
_gradient(x...) = ReverseDiff.gradient(x...) # ??

@tullio grad=Base
@testset "gradients: ReverseDiff + DiffRules" begin include("gradients.jl") end

@tullio grad=Dual
@testset "gradients: ReverseDiff + ForwardDiff" begin include("gradients.jl") end

@info @sprintf("ReverseDiff tests took %.1f seconds", time()-t6)
=#

#===== Yota =====#
#=
t7 = time()
using Yota

GRAD = :Yota
_gradient(x...) = Yota.grad(x...)[2]

@tullio grad=Base
@testset "gradients: Yota + DiffRules" begin include("gradients.jl") end

@tullio grad=Dual
@testset "gradients: Yota + ForwardDiff" begin include("gradients.jl") end

@info @sprintf("Yota tests took %.1f seconds", time()-t7)
=#

#===== LoopVectorization =====#

t8 = time()
using LoopVectorization
using VectorizationBase

#=
if isdefined(VectorizationBase, :SVec) # LV version 0.8, for Julia <= 1.5
    using VectorizationBase: SVec, Mask
else # LV version 0.9, supports Julia >= 1.5
    using VectorizationBase: Vec, Mask
    SVec{N,T} = Vec{N,T}
end

@testset "LoopVectorization onlyone" begin
    ms = Mask{4,UInt8}(0x03) # Mask{4,Bool}<1, 1, 0, 0>
    sv = SVec{4,Int}(1,2,3,4) # SVec{4,Int64}<1, 2, 3, 4>

    # preliminaries:
    @test Tullio.allzero(sv) === false
    @test Tullio.allzero(zero(sv)) === true

    @test Tullio.anyone(ms) === true

    # the main function:
    @test Tullio.onlyone(false, 0) === false
    @test Tullio.onlyone(true, 0) === true
    @test Tullio.onlyone(true, 1) === false

    @test Tullio.onlyone(ms, 0) === Mask{4}(0x02)
    # @test Tullio.onlyone(ms, 0).u == 0x02
    @test Tullio.onlyone(ms, sv) === Mask{4}(0x00)
    # @test Tullio.onlyone(ms, sv).u == 0x00
    @test Tullio.onlyone(ms, zero(sv)) === Mask{4}(0x02)
    # @test Tullio.onlyone(ms, zero(sv)).u == 0x02
end
=#

@testset "parsing + LoopVectorization" begin include("parsing.jl") end

if test_group != "3" # Github CI fails, on some runs, "ERROR: Package Tullio errored during testing (received signal: KILL)"
    # https://github.com/mcabbott/Tullio.jl/pull/57/checks?check_run_id=1753332805

    using Tracker
    GRAD = :Tracker
    _gradient(x...) = Tracker.gradient(x...)

    @tullio grad=Base
    @testset "gradients: Tracker + DiffRules + LoopVectorization" begin include("gradients.jl") end

    @tullio grad=Dual
    @testset "gradients: Tracker + ForwardDiff + LoopVectorization" begin include("gradients.jl") end

end

@info @sprintf("LoopVectorization tests took %.1f seconds", time()-t8)

@tullio avx=false

#===== TensorOperations =====#

t9 = time()
using TensorOperations

using Tracker
GRAD = :Tracker
_gradient(x...) = Tracker.gradient(x...)

@tullio grad=Base
@testset "gradients: Tracker + TensorOperations" begin include("gradients.jl") end

if VERSION < v"1.6-" # Zygote isn't working on 1.6

using Zygote
GRAD = :Zygote
_gradient(x...) = Zygote.gradient(x...)

@tullio grad=Base
@testset "gradients: Zygote + TensorOperations" begin include("gradients.jl") end

@testset "complex gradients with TensorOperations" begin

    x0 = [1 2; 3 4] .+ [5im 0; 7im -8im]

    @testset "analytic" begin

        g1 = _gradient(x -> real(sum(x * x)), x0)[1]
        g1i = _gradient(x -> imag(sum(x * x)), x0)[1]
        @test g1 ≈ _gradient(x -> real(sum(@tullio y[i,j] := x[i,k] * x[k,j])), x0)[1]
        @test g1i ≈ _gradient(x -> imag(sum(@tullio y[i,j] := x[i,k] * x[k,j])), x0)[1]

    end
    @testset "non-analytic" begin

        g2 = _gradient(x -> real(sum(x * x')), x0)[1]
        g2i = _gradient(x -> imag(sum(x * x')), x0)[1] # zero
        @test_broken g2 ≈ _gradient(x -> real(sum(@tullio y[i,j] := x[i,k] * conj(x[j,k]))), x0)[1]
        @test_broken g2i ≈ _gradient(x -> imag(sum(@tullio y[i,j] := x[i,k] * conj(x[j,k]))), x0)[1]

    end
end

end # VERSION

@testset "parsing + TensorOperations" begin include("parsing.jl") end # testing correct fallback

@info @sprintf("TensorOperations tests took %.1f seconds", time()-t9)

#===== done! =====#


text:
using SweepContractor

L=6; d=2

LTN = LabelledTensorNetwork{Tuple{Int,Int}}()
for i∈1:L, j∈1:L
    adj=[
        (mod1(i-1,L),mod1(j,L)),
        (mod1(i+1,L),mod1(j,L)),
        (mod1(i,L),mod1(j-1,L)),
        (mod1(i,L),mod1(j+1,L))
    ]
    LTN[i,j] = Tensor(adj, randn(d,d,d,d), i+0.01*rand(), j+0.01*rand())
end

println("example_2d_periodic")
for χ ∈ d.^(4:9)
    sweep = sweep_contract(LTN,χ, 2*χ)
    println("χ=$χ:\t",ldexp(sweep...))
end


text:
using Luxor

include("../tools/play.jl")
include("../test/environments/collect_points_cpu.jl")

config_environment = OrderedDict()
config_environment["maze_columns"] = 7
config_environment["maze_rows"] = 10
config_environment["maze_cell_size"] = 80
config_environment["agent_radius"] = 12
config_environment["point_radius"] = 8
config_environment["agent_movement_range"] = 10.0
config_environment["use_sensors"] = true
config_environment["reward_per_collected_positive_point"] = 500.0
config_environment["reward_per_collected_negative_point"] = -700.0
config_environment["number_time_steps"] = 500

configuration = CollectPointsCpuCfg(config_environment)

width = configuration.maze_cell_size * configuration.maze_columns
height = configuration.maze_cell_size * configuration.maze_rows
number_iterations = configuration.number_time_steps

env = CollectPointsCpu(configuration)

@play width height number_iterations begin

    step(env)

    # Draw white background
    background("white")

    # Draw agent
    sethue("green")
    circle(env.agent_position_x, env.agent_position_y, configuration.agent_radius, :fill)

    # Draw positive point
    sethue("blue")
    circle(env.point_x, env.point_y, configuration.point_radius, :fill)

    # Draw negative point
    sethue("red")
    circle(env.negative_point_x, env.negative_point_y, configuration.point_radius, :fill)

    # Draw maze
    sethue("black")
    setline(2)
    for cell_x = 1:configuration.maze_columns
        for cell_y = 1:configuration.maze_rows

            x_left, x_right, y_top, y_bottom = get_coordinates_maze_cell(cell_x, cell_y, configuration.maze_cell_size)

            cell = env.maze[cell_x, cell_y]

            # Draw walls
            if cell.walls[North] == true
                line(Point(x_left, y_top), Point(x_right, y_top), :stroke)
            end

            if cell.walls[South] == true
                line(Point(x_left, y_bottom), Point(x_right, y_bottom), :stroke)
            end

            if cell.walls[East] == true
                line(Point(x_right, y_top), Point(x_right, y_bottom), :stroke)
            end

            if cell.walls[West] == true
                line(Point(x_left, y_top), Point(x_left, y_bottom), :stroke)
            end

        end
    end

    sleep(0.05)

end

println("Finished")


text:
function ybus!(Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, b::Line)

    Y_l = (1 / (b.r + b.x * 1im))

    Y11 = Y_l + (1im * b.b.from);
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.from.number] += Y11;

    Y12 = -Y_l;
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.to.number] += Y12;
    #Y21 = Y12
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.from.number] += Y12;

    Y22 = Y_l + (1im * b.b.to);
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.to.number] += Y22;

end

function ybus!(Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, b::Transformer2W)

    Y_t = 1 / (b.r + b.x * 1im)

    Y11 = Y_t
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.from.number] += Y11;
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.to.number] += -Y_t;
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.from.number] += -Y_t;
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.to.number] += Y_t + (1im * b.primaryshunt);

end

function ybus!(Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, b::TapTransformer)

    Y_t = 1 / (b.r + b.x * 1im)
    c = 1 / b.tap

    Y11 = (Y_t * c^2);
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.from.number] += Y11;
    Y12 = (-Y_t*c) ;
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.to.number] += Y12;
    #Y21 = Y12
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.from.number] += Y12;
    Y22 = Y_t;
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.to.number] += Y22 + (1im * b.primaryshunt);

end

# TODO: Add testing for Ybus of a system with a PS Transformer
function ybus!(Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, b::PhaseShiftingTransformer)

    Y_t = 1 / (b.r + b.x * 1im)
    tap =  (b.tap * exp(b.α * 1im))
    c_tap =  (b.tap * exp(-1*b.α * 1im))

    Y11 = (Y_t/abs(tap)^2);
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.from.number] += Y11;
    Y12 = (-Y_t/c_tap);
    Ybus[b.connectionpoints.from.number,
        b.connectionpoints.to.number] += Y12;
    Y21 = (-Y_t/tap);
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.from.number] += Y21;
    Y22 = Y_t;
    Ybus[b.connectionpoints.to.number,
        b.connectionpoints.to.number] += Y22 + (1im * b.primaryshunt);

end

#=
function ybus!(Ybus::SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}, b::Transformer3W)

    @warn "Data contains a 3W transformer"

    Y11 = (1 / (b.line.r + b.line.x * 1im) + (1im * b.line.b) / 2);
    Ybus[b.line.connectionpoints.from.number,
        b.line.connectionpoints.from.number] += Y11;
    Y12 = (-1 ./ (b.line.r + b.line.x * 1im));
    Ybus[b.line.connectionpoints.from.number,
        b.line.connectionpoints.to.number] += Y12;
    #Y21 = Y12
    Ybus[b.line.connectionpoints.to.number,
        b.line.onnectionpoints[1].number] += Y12;
    #Y22 = Y11;
    Ybus[b.line.connectionpoints.to.number,
        b.line.connectionpoints.to.number] += Y11;

    y = 1 / (b.transformer.r + b.transformer.x * 1im)
    y_a = y / (b.transformer.tap * exp(b.transformer.α * 1im * (π / 180)))
    c = 1 / b.transformer.tap

    Y11 = (y_a + y * c * (c - 1) + (b.transformer.zb));
    Ybus[b.transformer.connectionpoints.from.number,
        b.transformer.connectionpoints.from.number] += Y11;
    Y12 = (-y_a) ;
    Ybus[b.transformer.connectionpoints.from.number,
        b.transformer.connectionpoints.to.number] += Y12;
    #Y21 = Y12
    Ybus[b.transformer.connectionpoints.to.number,
        b.transformer.connectionpoints.from.number] += Y12;
    Y22 = (y_a + y * (1 - c)) ;
    Ybus[b.transformer.connectionpoints.to.number,
        b.transformer.connectionpoints.to.number] += Y22;

end
=#

function build_ybus(buscount::Int64, branches::Array{T}) where {T <: Branch}

    Ybus = SparseArrays.spzeros(Complex{Float64}, buscount, buscount)

    for b in branches

        if b.name == "init"
            @error "The data in Branch is incomplete" # TODO: raise error here?
        end

        ybus!(Ybus, b)

    end

    return Ybus

end


text:
@testset "felix output" begin
    data = ComplexF64[1 4im;  # will be written as ComplexF32
                      2 5im;
                      3 6im]

    io = IOBuffer()
    bandwidth = (123u"Hz", 456u"Hz")
    frequency = 500.5u"MHz"
    save_felix(io, data; bandwidth=bandwidth, frequency=frequency)

    words = reinterpret(UInt32, take!(io))

    # Just some basic tests here, to check that the metadata gets in the
    # correct place and the data is appended.
    @test words[1] == 0x04030201
    @test words[2] == 256 # Size of header

    header = words[3:3+256-1]
    @test header[121] == 1
    @test header[122] == size(data,1)
    @test header[123] == size(data,2)

    @test all(header[128:133] .== 1)

    fheader = reinterpret(Float32, header)
    @test fheader[134] == frequency/u"MHz"
    @test fheader[135] == frequency/u"MHz"

    @test fheader[140] == bandwidth[1]/u"Hz"
    @test fheader[141] == bandwidth[2]/u"Hz"

    @test words[3+256] == size(data,1)*2
    @test reinterpret(ComplexF32, words[3+256 .+ (1:6)]) == conj.(data[:,1])
    @test words[3+256 + 6 + 1] == size(data,1)*2
    @test reinterpret(ComplexF32, words[3+256 + 6 + 1 .+ (1:6)]) == conj.(data[:,2])
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
# DEMAND GENERATOR
libs = Libraries("/home/kreiton/.julia/dev/DispatchOps/data/laptri", false)
libs.demand_realization = noisify_fixed(libs.demand_forecast, 0)

# INITIALIZATION
sim = Simulation(
    terminating_timestep = 2,
    libs = libs,
    cons = Constants(1, 0.001)
)

# RUN
run!(sim)


text:
"""
    arrowsplit(line: String) -> String[]

Return a line as a vector, splitting on an arrow " -> ".

This is used for **@piping.**

**Example**
```julia
arrowsplit("a -> b")  # returns ["a", "b"]
```
"""
arrowsplit(line) = split(line, " -> ")


text:
const version = "0.1.0"
const HGS_CVRP_SRC = "https://github.com/chkwon/HGS-CVRP/archive/v$(version).tar.gz"

const HGS_CVRP_WIN = "https://github.com/chkwon/HGS_CVRP_jll.jl/releases/download/HGS_CVRP-v0.1.0%2B0/libhgscvrp.v0.1.0.x86_64-w64-mingw32-cxx11.tar.gz"

const SRC_DIR = "HGS-CVRP-$version"

import BinDeps

function build_HGS()
    cd(SRC_DIR)
    mkdir("build")
    cd("build")
    run(`cmake -DCMAKE_BUILD_TYPE=Release ../`)
    run(`make hgscvrp`)

    return joinpath(@__DIR__, SRC_DIR, "build", "libhgscvrp")
end

function download_HGS_LIB_WIN()
    rm("win", recursive=true, force=true)
    mkdir("win")
    lib_tarball = joinpath(@__DIR__, "win", "hgs_win.tar.gz")
    download(HGS_CVRP_WIN, lib_tarball)
    cd("win")
    run(BinDeps.unpack_cmd(lib_tarball, ".", ".gz", ".tar"))

    return joinpath(@__DIR__, "win", "bin", "libhgscvrp")
end

function download_HGS()
    src_tarball = joinpath(@__DIR__, "hgs_cvrp.tar.gz")
    download(HGS_CVRP_SRC, src_tarball)
    rm(SRC_DIR, recursive=true, force=true)
    run(BinDeps.unpack_cmd(src_tarball, ".", ".gz", ".tar"))
end

function install_HGS()
    lib = get(ENV, "HGS_CVRP_SHARED_LIBRARY", nothing)
    if !haskey(ENV, "HGS_CVRP_SHARED_LIBRARY")
        if Sys.iswindows()
            lib = download_HGS_LIB_WIN()
        else
            download_HGS()
            lib = build_HGS()
        end
        ENV["HGS_CVRP_SHARED_LIBRARY"] = lib
    end

    if lib === nothing
        error("Environment variable `HGS_CVRP_SHARED_LIBRARY` not found.")
    else
        # pr2392 = joinpath(@__DIR__, LKH_VERSION, "pr2392.par")
        # run(`$(executable) $(pr2392)`)
    end

    open(joinpath(@__DIR__, "deps.jl"), "w") do io
        write(io, "const LIBHGSCVRP = \"$(escape_string(lib))\"\n")
    end
end

install_HGS()


text:
module ClimaCorePlots

import RecipesBase
import TriplotBase

import ClimaCore:
    ClimaCore,
    DataLayouts,
    Fields,
    Geometry,
    Meshes,
    Operators,
    Topologies,
    Spaces

RecipesBase.@recipe function f(space::Spaces.SpectralElementSpace2D;)
    quad = Spaces.quadrature_style(space)
    quad_name = Base.typename(typeof(quad)).name
    dof = Spaces.Quadratures.degrees_of_freedom(quad)

    topology = Spaces.topology(space)
    mesh = topology.mesh
    @assert mesh isa Meshes.RectilinearMesh "plotting only defined for RectilinearMesh"

    n1 = Meshes.nelements(mesh.intervalmesh1)
    n2 = Meshes.nelements(mesh.intervalmesh2)

    coord_field = Fields.coordinate_field(space)
    x1coord = vec(parent(coord_field)[:, :, 1, :])
    x2coord = vec(parent(coord_field)[:, :, 2, :])

    coord_symbols = propertynames(coord_field)

    seriestype := :scatter
    title --> "$n1 × $n2 $quad_name{$dof} element space"
    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    marker --> :cross
    markersize --> 1
    seriescolor --> :blue
    legend --> false

    (x1coord, x2coord)
end

RecipesBase.@recipe function f(space::Spaces.ExtrudedFiniteDifferenceSpace)
    coord_field = Fields.coordinate_field(space)
    data = Fields.field_values(coord_field)
    Ni, Nj, _, Nv, Nh = size(data)

    #TODO: assumes VIFH layout
    @assert Nj == 1 "plotting only defined for 1D extruded fields"

    hspace = space.horizontal_space

    quad = Spaces.quadrature_style(hspace)
    quad_name = Base.typename(typeof(quad)).name
    dof = Spaces.Quadratures.degrees_of_freedom(quad)

    coord_symbols = propertynames(coord_field)
    hcoord = vec(parent(coord_field)[:, :, 1, :])
    vcoord = vec(parent(coord_field)[:, :, 2, :])

    stagger = space.staggering isa Spaces.CellCenter ? :center : :face

    seriestype := :scatter
    title --> "$Nh $quad_name{$dof} element × $Nv level $stagger space"
    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    marker --> :cross
    markersize --> 1
    seriescolor --> :blue
    legend --> false

    (hcoord, vcoord)
end

RecipesBase.@recipe function f(field::Fields.FiniteDifferenceField)
    # unwrap the data to plot
    space = axes(field)
    coord_field = Fields.coordinate_field(space)

    xdata = parent(field)[:, 1]
    ydata = parent(Spaces.coordinates_data(space))[:, 1]

    coord_symbols = propertynames(coord_field)

    # set the plot attributes
    xguide --> "value"
    yguide --> (
        if field isa Spaces.FaceFiniteDifferenceSpace
            "$(coord_symbols[1]) faces"
        else
            "$(coord_symbols[1]) centers"
        end
    )

    # fix the ylim to the column space (domain)
    ylims --> extrema(ydata)
    linewidth --> 2
    legend --> false

    (xdata, ydata)
end

RecipesBase.@recipe function f(
    field::Fields.SpectralElementField2D;
    interpolate = 10,
)
    @assert interpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    # compute the interpolated data to plot
    space = axes(field)
    topology = Spaces.topology(space)
    mesh = topology.mesh

    Nu = interpolate
    coord_field = Fields.coordinate_field(space)

    M_coords = Operators.matrix_interpolate(coord_field, Nu)
    M = Operators.matrix_interpolate(field, Nu)

    domain = Meshes.domain(mesh)
    x1min = Geometry.component(domain.interval1.coord_min, 1)
    x2min = Geometry.component(domain.interval2.coord_min, 1)
    x1max = Geometry.component(domain.interval1.coord_max, 1)
    x2max = Geometry.component(domain.interval2.coord_max, 1)

    # our interpolated field is transposed
    x1coord = [Geometry.component(pt, 1) for pt in M_coords[:, 1]]
    x2coord = [Geometry.component(pt, 2) for pt in M_coords[1, :]]

    coord_symbols = propertynames(coord_field)

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    (x1coord, x2coord, M')
end

function _unfolded_pannel_matrix(field, interpolate)
    space = axes(field)
    topology = Spaces.topology(space)
    mesh = topology.mesh
    nelem = Meshes.nelements(mesh)
    panel_size = mesh.ne

    quad_from = Spaces.quadrature_style(space)
    quad_to = Spaces.Quadratures.Uniform{interpolate}()
    Imat = Spaces.Quadratures.interpolation_matrix(Float64, quad_to, quad_from)

    dof = interpolate

    pannel_range(i) =
        ((panel_size * dof) * (i - 1) + 1):((panel_size * dof) * i)

    # construct a matrix to fill in the rotated / flipped pannel data
    unfolded_panels =
        fill(NaN, ((panel_size * dof) * 3, (panel_size * dof) * 4))

    # temporary pannels as we have to rotate / flip some and not all operators are in place
    # TODO: inefficient memory wise, but good enough for now
    panels = [fill(NaN, (panel_size * dof, panel_size * dof)) for _ in 1:6]

    interpolated_data =
        DataLayouts.IJFH{Float64, interpolate}(Array{Float64}, nelem)
    field_data = Fields.field_values(field)

    Operators.tensor_product!(interpolated_data, field_data, Imat)

    # element index ordering defined by a specific layout
    eidx = 1
    for panel_idx in 1:6
        panel_data = panels[panel_idx]
        # elements are ordered along fastest axis defined in sphere box mesh
        for ex2 in 1:panel_size, ex1 in 1:panel_size
            # compute the nodal extent index range for this element
            x1_nodal_range = (dof * (ex1 - 1) + 1):(dof * ex1)
            x2_nodal_range = (dof * (ex2 - 1) + 1):(dof * ex2)
            # transpose the data as our plotting axis order is
            # reverse nodal element order (x1 axis varies fastest)
            data_element = permutedims(parent(interpolated_data)[:, :, 1, eidx])
            panel_data[x2_nodal_range, x1_nodal_range] = data_element
            eidx += 1
        end
    end

    # (px, py, rot) for each panel
    # px, py are the locations of the panel
    # rot is number of clockwise rotations (0:3)
    # https://extranet.gfdl.noaa.gov/~atw/ferret/cubed_sphere/
    # Equatorial strip, poles connected through the Americas
    panel_locations =
        [(4, 2, 0), (1, 2, 0), (3, 3, 2), (2, 2, 1), (3, 2, 1), (3, 1, 1)]
    for (i, (px, py, rot)) in enumerate(panel_locations)
        unfolded_panels[pannel_range(py), pannel_range(px)] .= if rot == 0
            panels[i]
        elseif rot == 1
            reverse(transpose(panels[i]), dims = 1)
        elseif rot == 2
            reverse(panels[i], dims = (1, 2))
        else
            reverse(transpose(panels[i]), dims = 2)
        end
    end
    #=
    unfolded_panels[pannel_range(1), pannel_range(2)] =
        reverse(panels[5], dims = 1)
    unfolded_panels[pannel_range(2), pannel_range(1)] =
        reverse(panels[4], dims = 2)
    unfolded_panels[pannel_range(2), pannel_range(2)] = transpose(panels[1])
    unfolded_panels[pannel_range(2), pannel_range(3)] = transpose(panels[2])
    unfolded_panels[pannel_range(2), pannel_range(4)] =
        reverse(panels[6], dims = 2)
    unfolded_panels[pannel_range(3), pannel_range(2)] =
        reverse(panels[3], dims = 2)
    =#
    return unfolded_panels
end

RecipesBase.@recipe function f(
    field::Fields.CubedSphereSpectralElementField2D;
    interpolate = 10,
)
    @assert interpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    unfolded_panels = _unfolded_pannel_matrix(field, interpolate)

    # construct the title for info about the field space
    space = axes(field)
    nelem = Topologies.nlocalelems(Spaces.topology(space))
    quad_from = Spaces.quadrature_style(space)
    dof_in = Spaces.Quadratures.degrees_of_freedom(quad_from)
    quad_from_name = Base.typename(typeof(quad_from)).name

    # set the plot attributes
    seriestype := :heatmap
    title --> "$nelem $quad_from_name{$dof_in} element space"
    xguide --> "panel x1"
    yguide --> "panel x2"
    seriescolor --> :balance

    (unfolded_panels)
end

RecipesBase.@recipe function f(
    field::Fields.Field{<:Any, S};
    level = nothing,
    hinterpolate = 10,
) where {
    S <: Spaces.ExtrudedFiniteDifferenceSpace{
        <:Any,
        <:Spaces.CubedSphereSpectralElementSpace2D,
    },
}
    @assert hinterpolate ≥ 1 "number of element quadrature points for uniform interpolation must be ≥ 1"

    space = axes(field)
    if level === nothing
        if space isa Spaces.FaceExtrudedFiniteDifferenceSpace
            level = Utilities.PlusHalf(0)
        else
            level = 1
        end
    end
    level_field = Fields.level(field, level)
    unfolded_panels = _unfolded_pannel_matrix(level_field, hinterpolate)

    # construct the title for info about the field space
    nlevel = Spaces.nlevels(space)
    nelem = Topologies.nlocalelems(Spaces.topology(space))
    quad_from = Spaces.quadrature_style(space)
    dof_in = Spaces.Quadratures.degrees_of_freedom(quad_from)
    quad_from_name = Base.typename(typeof(quad_from)).name

    # set the plot attributes
    seriestype := :heatmap
    title -->
    "level $level of $nlevel × $nelem $quad_from_name{$dof_in} element space"
    xguide --> "panel x1"
    yguide --> "panel x2"
    seriescolor --> :balance

    (unfolded_panels)
end

RecipesBase.@recipe function f(
    field::Fields.Field{<:Any, S};
    hinterpolate = 0,
    ncolors = 256,
) where {
    S <: Spaces.ExtrudedFiniteDifferenceSpace{
        <:Any,
        <:Spaces.SpectralElementSpace1D,
    },
}
    data = Fields.field_values(field)
    Ni, Nj, _, Nv, Nh = size(data)

    space = axes(field)
    @assert Nj == 1

    coord_symbols = propertynames(Fields.coordinate_field(space))
    hcoord_field = getproperty(Fields.coordinate_field(space), 1)
    vcoord_field = getproperty(Fields.coordinate_field(space), 2)

    if hinterpolate ≥ 1
        Nu = hinterpolate
        uquad = Spaces.Quadratures.ClosedUniform{Nu}()
        M_hcoord = Operators.matrix_interpolate(hcoord_field, uquad)
        M_vcoord = Operators.matrix_interpolate(vcoord_field, uquad)
        M_data = Operators.matrix_interpolate(field, uquad)

        hcoord_data = vec(M_hcoord)
        vcoord_data = vec(M_vcoord)
        data = vec(M_data)
        triangles = hcat(Spaces.triangles(Nv, Nu, Nh)...)'
    else
        hcoord_data = vec(parent(hcoord_field))
        vcoord_data = vec(parent(vcoord_field))
        data = vec(parent(data))
        triangles = hcat(Spaces.triangles(Nv, Ni, Nh)...)'
    end
    cmap = range(extrema(data)..., length = ncolors)

    z = TriplotBase.tripcolor(
        hcoord_data,
        vcoord_data,
        data,
        triangles,
        cmap;
        bg = NaN,
        # unique number of nodal element coords
        # (number of nodal values minus number of shared faces)
        px = (Ni * Nh) - (Nh - 1),
        py = Nv,
    )

    # set the plot attributes
    seriestype := :heatmap

    xguide --> "$(coord_symbols[1])"
    yguide --> "$(coord_symbols[2])"
    seriescolor --> :balance

    # some plots backends need coordinates in sorted order
    (sort(unique(hcoord_data)), sort(unique(vcoord_data)), z')
end

end # module


text:
__precompile__()

module BoxAtlasControl

include("qp_input.jl")
include("state.jl")

end # module


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# deprecated require
module OldRequire

function find_in_path(name::AbstractString)
    isabspath(name) && return name
    isfile(name) && return abspath(name)
    base = name
    if endswith(name,".jl")
        base = name[1:end-3]
    else
        name = string(base,".jl")
        isfile(name) && return abspath(name)
    end
    for prefix in [Pkg.dir(); LOAD_PATH]
        path = joinpath(prefix, name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, base, "src", name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, name, "src", name)
        isfile(path) && return abspath(path)
    end
    return nothing
end

find_in_node1_path(name) = myid()==1 ?
    find_in_path(name) : remotecall_fetch(find_in_path, 1, name)

# Store list of files and their load time
package_list = Dict{String,Float64}()
# to synchronize multiple tasks trying to require something
package_locks = Dict{String,Any}()

# only broadcast top-level (not nested) requires and reloads
toplevel_load = true

require(fname::AbstractString) = require(String(fname))
require(f::AbstractString, fs::AbstractString...) = (require(f); for x in fs require(x); end)

function require(name::String)
    path = find_in_node1_path(name)
    path === nothing && error("$name not found")

    if myid() == 1 && toplevel_load
        refs = Any[ @spawnat p _require(path) for p in filter(x->x!=1, procs()) ]
        _require(path)
        for r in refs; wait(r); end
    else
        _require(path)
    end
    nothing
end

function _require(path)
    global toplevel_load
    if haskey(package_list,path)
        wait(package_locks[path])
    else
        last = toplevel_load
        toplevel_load = false
        try
            reload_path(path)
        finally
            toplevel_load = last
        end
    end
end

# remote/parallel load

function source_path(default::Union{AbstractString,Void}="")
    t = current_task()
    while true
        s = t.storage
        if s !== nothing && haskey(s, :SOURCE_PATH)
            return s[:SOURCE_PATH]
        end
        if t === t.parent
            return default
        end
        t = t.parent
    end
end

function include_from_node1(path::AbstractString)
    prev = source_path(nothing)
    path = (prev === nothing) ? abspath(path) : joinpath(dirname(prev),path)
    tls = task_local_storage()
    tls[:SOURCE_PATH] = path
    local result
    try
        if myid()==1
            # sleep a bit to process file requests from other nodes
            nprocs()>1 && sleep(0.005)
            result = Core.include(path)
            nprocs()>1 && sleep(0.005)
        else
            result = include_string(remotecall_fetch(readstring, 1, path), path)
        end
    finally
        if prev === nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
    result
end

function reload_path(path::AbstractString)
    had = haskey(package_list, path)
    if !had
        package_locks[path] = RemoteChannel()
    end
    package_list[path] = time()
    tls = task_local_storage()
    prev = pop!(tls, :SOURCE_PATH, nothing)
    try
        eval(Main, :(Base.include_from_node1($path)))
    catch e
        had || delete!(package_list, path)
        rethrow(e)
    finally
        if prev !== nothing
            tls[:SOURCE_PATH] = prev
        end
    end
    if !isready(package_locks[path])
        put!(package_locks[path],nothing)
    end
    nothing
end

end # module

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
## Use the astropy python package to read a votable.

module GASS

using PyCall
using DataFrames
using Printf
using LsqFit , FFTW
using Statistics , Distributions ,Random
using StatsBase

import CSV , JLD
import PyPlot

include("types.jl")
export inputCfg , observation , subarrayParameters , weight , GA , cfg , synthbeam,  population

include("inputcfg.jl")
export input_parameters , parse_input , parse_parameters , read_cfg

include("utils.jl")
export init_pop

include("simobserve.jl")
export calc_baselines , calc_uv , calc_dirtybeam , fit_beam , calc_mrs

include("evolution.jl")
export fitness_subarray , create_population , get_elitism , get_parents , get_crossover1 ,
  get_mutation , get_evolution , gass_optimization

include("output.jl")
export save_CASAarr , sort_population , full_run

end

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
"""
Store pre-allocated memory in this data structures, which contains
intermediate values used in the ELBO calculation.
"""
struct HessianSubmatrices{NumType <: Number}
    u_u::Matrix{NumType}
    shape_shape::Matrix{NumType}
end

"""
Pre-allocated memory for efficiently accumulating certain sub-matrices
of the E_G_s and E_G2_s Hessian.

Args:
    NumType: The numeric type of the hessian.
    i: The type of celestial source, from 1:NUM_SOURCE_TYPES
"""
function HessianSubmatrices(NumType::DataType, i::Int)
    @assert 1 <= i <= NUM_SOURCE_TYPES
    shape_p = length(shape_standard_alignment[i])

    u_u = zeros(NumType, 2, 2)
    shape_shape = zeros(NumType, shape_p, shape_p)
    HessianSubmatrices{NumType}(u_u, shape_shape)
end

struct ElboIntermediateVariables{NumType <: Number}
    # Vectors of star and galaxy bvn quantities from all sources for a pixel.
    # The vector has one element for each active source, in the same order
    # as ea.active_sources.

    fs0m::SensitiveFloat{NumType}
    fs1m::SensitiveFloat{NumType}

    # Brightness values for a single source
    E_G_s::SensitiveFloat{NumType}
    E_G2_s::SensitiveFloat{NumType}
    var_G_s::SensitiveFloat{NumType}

    # Subsets of the Hessian of E_G_s and E_G2_s that allow us to use BLAS
    # functions to accumulate Hessian terms. There is one submatrix for
    # each celestial object type in 1:NUM_SOURCE_TYPES
    E_G_s_hsub_vec::Vector{HessianSubmatrices{NumType}}
    E_G2_s_hsub_vec::Vector{HessianSubmatrices{NumType}}

    # Expected pixel intensity and variance for a pixel from all sources.
    E_G::SensitiveFloat{NumType}
    var_G::SensitiveFloat{NumType}

    # Pre-allocated memory for the gradient and Hessian of combine functions.
    combine_grad::Vector{NumType}
    combine_hess::Matrix{NumType}

    # A placeholder for the log term in the ELBO.
    elbo_log_term::SensitiveFloat{NumType}

    # The ELBO itself.
    elbo::SensitiveFloat{NumType}

    active_pixel_counter::Ref{Int64}
    inactive_pixel_counter::Ref{Int64}
end

"""
Args:
    - num_active_sources: The number of actives sources (with derivatives)
    - calculate_gradient: If false, only calculate values
    - calculate_hessian: If false, only calculate gradients. Note that if
                calculate_gradient = false, then hessians will not be
                calculated irrespective of the value of calculate_hessian.
"""
function ElboIntermediateVariables(NumType::DataType,
                                   num_active_sources::Int,
                                   calculate_gradient::Bool=true,
                                   calculate_hessian::Bool=true)
    @assert NumType <: Number

    # fs0m and fs1m accumulate contributions from all bvn components
    # for a given source.
    fs0m = SensitiveFloat{NumType}(length(StarPosParams), 1,
                                calculate_gradient, calculate_hessian)
    fs1m = SensitiveFloat{NumType}(length(GalaxyPosParams), 1,
                                calculate_gradient, calculate_hessian)

    E_G_s = SensitiveFloat{NumType}(length(CanonicalParams), 1,
                                    calculate_gradient, calculate_hessian)
    E_G2_s = SensitiveFloat(E_G_s)
    var_G_s = SensitiveFloat(E_G_s)

    E_G_s_hsub_vec =
        HessianSubmatrices{NumType}[ HessianSubmatrices(NumType, i) for i=1:NUM_SOURCE_TYPES ]
    E_G2_s_hsub_vec =
        HessianSubmatrices{NumType}[ HessianSubmatrices(NumType, i) for i=1:NUM_SOURCE_TYPES ]

    E_G = SensitiveFloat{NumType}(length(CanonicalParams), num_active_sources,
                                  calculate_gradient, calculate_hessian)
    var_G = SensitiveFloat(E_G)

    combine_grad = zeros(NumType, 2)
    combine_hess = zeros(NumType, 2, 2)

    elbo_log_term = SensitiveFloat(E_G)
    elbo = SensitiveFloat(E_G)

    ElboIntermediateVariables{NumType}(
        fs0m, fs1m,
        E_G_s, E_G2_s, var_G_s, E_G_s_hsub_vec, E_G2_s_hsub_vec,
        E_G, var_G, combine_grad, combine_hess,
        elbo_log_term, elbo, 0, 0)
end

function zero!(elbo_vars::ElboIntermediateVariables{T}) where {T<:Number}
    SensitiveFloats.zero!(elbo_vars.fs0m)
    SensitiveFloats.zero!(elbo_vars.fs1m)
    SensitiveFloats.zero!(elbo_vars.E_G_s)
    SensitiveFloats.zero!(elbo_vars.E_G2_s)
    SensitiveFloats.zero!(elbo_vars.var_G_s)

    for i in 1:NUM_SOURCE_TYPES
        fill!(elbo_vars.E_G_s_hsub_vec[i].u_u, zero(T))
        fill!(elbo_vars.E_G_s_hsub_vec[i].shape_shape, zero(T))
        fill!(elbo_vars.E_G2_s_hsub_vec[i].u_u, zero(T))
        fill!(elbo_vars.E_G2_s_hsub_vec[i].shape_shape, zero(T))
    end

    SensitiveFloats.zero!(elbo_vars.E_G)
    SensitiveFloats.zero!(elbo_vars.var_G)

    fill!(elbo_vars.combine_grad, zero(T))
    fill!(elbo_vars.combine_hess, zero(T))

    SensitiveFloats.zero!(elbo_vars.elbo_log_term)
    SensitiveFloats.zero!(elbo_vars.elbo)
end

"""
If Infs/NaNs have crept into the ELBO evaluation (a symptom of poorly conditioned optimization),
this helps catch them immediately.
"""
function assert_all_finite{NumType <: Number}(sf::SensitiveFloat{NumType})
    @assert isfinite(sf.v[]) "Value is Inf/NaNs"
    @assert all(isfinite, sf.d) "Gradient contains Inf/NaNs"
    @assert all(isfinite, sf.h) "Hessian contains Inf/NaNs"
end

"""
Some parameter to a function has invalid values. The message should explain what parameter is
invalid and why.
"""
struct InvalidInputError <: Exception
    message::String
end

"""
ElboArgs stores the arguments needed to evaluate the variational objective
function.
"""
struct ElboArgs
    # the overall number of sources: we don't necessarily visit them
    # all or optimize them all, but if we do visit a pixel where any
    # of these are active, we use it in the elbo calculation
    S::Int64

    # number of active sources (see below)
    Sa::Int64

    # the number of images
    N::Int64

    # The number of components in the point spread function.
    psf_K::Int64
    images::Vector{Image}

    # subimages is a better name for patches: regions of an image
    # around a particular light source
    patches::Matrix{SkyPatch}

    # the sources to optimize
    active_sources::Vector{Int}

    # If false, elbo = elbo_likelihood
    include_kl::Bool
end

function ElboArgs(
            images::Vector{Image},
            patches::Matrix{SkyPatch},
            active_sources::Vector{Int};
            psf_K::Int=2,
            include_kl::Bool=true)
    S = size(patches, 1)
    Sa = length(active_sources)
    N = length(images)

    @assert size(patches, 2) == N
    @assert psf_K > 0

    @assert(length(active_sources) <= 5 || !calculate_hessian,
            "too many active_sources to store a hessian")

    ElboArgs(S, Sa, N, psf_K, images, patches,
             active_sources, include_kl)
end


text:
module SSLTestRescale
using SimpleSDMLayers
using Test

M1 = rand(Float64, (5, 10)).*4
S1 = SimpleSDMPredictor(M1, 0.0, 1.0, 0.0, 1.0)

M2 = (rand(Float64, (5, 10)).+1).*5
S2 = SimpleSDMPredictor(M2, 0.0, 1.0, 0.0, 1.0)

@test extrema(rescale(S1, S2)) == extrema(S2)
@test extrema(rescale(S2, S1)) == extrema(S1)

rescale!(S1, [0.0, 0.5, 1.0])
@test sort(unique(collect(S1))) == [0.0, 0.5, 1.0]

end


text:
# sample from weights
# import IncrementalInference: AliasingScalarSampler

"""
$(TYPEDEF)

Sampler from intensity map given Euclidean domain `x` and probability weights `p_x`.

Example

`AliasingScalarSampler(x::Vector{<:Real}, p_x::Vector{<:Real}; SNRfloor::Float64=0.0)`
"""
struct AliasingScalarSampler
  domain::Vector{Float64}
  weights::StatsBase.ProbabilityWeights
  AliasingScalarSampler() = new()
  AliasingScalarSampler(x::Vector{<:Real}, p_x::Vector{<:Real}; SNRfloor::Float64=0.0) = begin

    # pxf = Float64.(p_x)
    # pxf .-= quantile(pxf, SNRfloor)
    # pxf[pxf .< 0.0] = 0.0
    # pxf ./= norm(pxf)
    # new implementation : pxf should be an empirical pmf before use in statsbase
    pxf = Float64.(p_x)
    pxf[pxf.<0.0] .= 0.0 # no negative values!
    pxf ./=sum(pxf)  # must sum to 1
        pxf2 = pxf .- quantile(pxf,SNRfloor) # remove lowest quantile
        pxf2[pxf2.<0.0] .= 0.0
        pxf2s = sum(pxf2)
        pxf[:] = 1e-10 < pxf2s ? pxf2 : pxf
        pxf ./= sum(pxf)
        sum(isnan.(pxf)) == 0 ? nothing : error("AliasingScalarSampler got NaN because of particular values in p_x")
        # pxf .-= quantile(pxf,SNRfloor) # remove lowest quantile
        # pxf[pxf.<0.0] .= 0.0
        # pxf ./=sum(pxf)
    wim = StatsBase.ProbabilityWeights(pxf)
    new(x, wim)
  end
end

function rand!(ass::AliasingScalarSampler, smpls::Array{Float64})
    StatsBase.alias_sample!(ass.domain, ass.weights, smpls)
    nothing
end

function rand(ass::AliasingScalarSampler, N::Int=1)
  smpls = zeros(N)
  rand!(ass,smpls)
  smpls
end

function convert(::Type{AliasingScalarSampler}, str::AS) where {AS <: AbstractString}
  sstr = split(split(str, "AliasingScalarSampler")[end], '[')
  sstr = sstr[length.(sstr) .> 2]
  ssstr = split.(sstr, ']')
  domain = parse.(Float64, strip.(split(ssstr[1][1], ',')) )
  weight = parse.(Float64, strip.(split(ssstr[2][1], ',')) )
  AliasingScalarSampler(domain, weight)
end
# str = "IncrementalInference.AliasingScalarSampler([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0,24.0, 25.0, 26.0, 27.0, 28.0, 29.0, 30.0], [0.156102, 0.117163, 0.221049, 0.275905, 0.0488494, 0.0731541, 0.107584,0.313848, 0.0309002, 0.0, 0.0384554, 0.155308, 0.276917, 0.0271168, 0.293263, 0.171316, 0.27459, 0.175323, 0.0535772, 0.181663, 0.295042, 0.104593, 0.0472137, 0.326016, 0.055283, 0.0737767, 0.302647, 0.0291257, 0.0206642, 0.223375])"


text:
using Neighborhood
using Distances

export garcia_almeida_embedding
export n_statistic

"""
    garcia_almeida_embedding(s; kwargs...) → Y, τ_vals, ts_vals, FNNs ,NS
A unified approach to properly embed a time series (`Vector` type) or a
set of time series (`Dataset` type) based on the papers of Garcia & Almeida
[^Garcia2005a],[^Garcia2005b].

## Keyword arguments

* `τs= 0:50`: Possible delay values `τs` (in sampling time units). For each of
  the `τs`'s the N-statistic gets computed.
* `w::Int = 1`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `r1 = 10`: The threshold, which defines the factor of tolerable stretching for
  the d_E1-statistic.
* `r2 = 2`: The threshold for the tolerable relative increase of the distance
  between the nearest neighbors, when increasing the embedding dimension.
* `fnn_thres= 0.05`: A threshold value defining a sufficiently small fraction
  of false nearest neighbors, in order to the let algorithm terminate and stop
  the embedding procedure (`0 ≤ fnn_thres < 1).
* `T::Int = 1`: The forward time step (in sampling units) in order to compute the
  `d_E2`-statistic (see algorithm description). Note that in the paper this is
  not a free parameter and always set to `T=1`.
* `metric = Euclidean()`: metric used for finding nearest neigbhors in the input
  phase space trajectory `Y`.
* `max_num_of_cycles = 50`: The algorithm will stop after that many cycles no matter what.

## Description
The method works iteratively and gradually builds the final embedding vectors
`Y`. Based on the `N`-statistic the algorithm
picks an optimal delay value `τ` for each embedding cycle as the first local
minimum of `N`. In case of multivariate embedding, i.e. when embedding a set of
time series (`s::Dataset`), the optimal delay value `τ` is chosen as the first
minimum from all minimum's of all considered `N`-statistics for each embedding
cycle. The range of considered delay values is determined in `τs` and for the
nearest neighbor search we respect the Theiler window `w`. After each embedding
cycle the FNN-statistic `FNNs` [^Hegger1999][^Kennel1992] is being checked and
as soon as this statistic drops below the threshold `fnn_thres`, the algorithm
breaks. In order to increase the  practability of the method the algorithm also
breaks, when the FNN-statistic `FNNs` increases . The final embedding vector is
stored in `Y` (`Dataset`). The chosen delay values for each embedding cycle are
stored in the `τ_vals` and the according time series number chosen for the
according delay value in `τ_vals` is stored in `ts_vals`. For univariate
embedding (`s::Vector`) `ts_vals` is a vector of ones of length `τ_vals`,
because there is simply just one time series to choose from. The function also
returns the `N`-statistic `NS` for each embedding cycle as an `Array` of
`Vector`s.

Notice that we were *not* able to reproduce the figures from the papers with our
implementation (which nevertheless we believe is the correct one).

[^Garcia2005a]: Garcia, S. P., Almeida, J. S. (2005). [Nearest neighbor embedding with different time delays. Physical Review E 71, 037204](https://doi.org/10.1103/PhysRevE.71.037204).
[^Garcia2005b]: Garcia, S. P., Almeida, J. S. (2005). [Multivariate phase space reconstruction by nearest neighbor embedding with different time delays. Physical Review E 72, 027205](https://doi.org/10.1103/PhysRevE.72.027205).
"""
function garcia_almeida_embedding(s::Vector{F}; τs = 0:50 , w::Int = 1,
    r1::Real = 10, r2::Real = 2, fnn_thres::Real = 0.05,
    T::Int = 1, metric = Euclidean(), max_num_of_cycles = 50) where {F<:Real}

    @assert 0 ≤ fnn_thres < 1 "Please select a valid breaking criterion, i.e. a threshold value `fnn_thres` ∈ [0 1)"
    @assert all(x -> x ≥ 0, τs)
    s_orig = s
    s = standardize(s) # especially important for fnn-computation
    # define actual phase space trajectory
    Y_act = Dataset(s)

    # compute nearest neighbor distances in the first embedding dimension for
    # FNN statistic
    vtree = KDTree(Y_act, metric)
    allNNidxs, NNdist_old = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    # preallocate output variables
    τ_vals = Int64[0]
    ts_vals = Int64[1]
    FNNs = Float64[]
    NS = Array{F}(undef, length(τs), max_num_of_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        Y_act, NNdist_new = garcia_embedding_cycle!(
            Y_act, flag, s, τs, w, counter, NS, τ_vals, metric, r1, r2, T,
            FNNs, fnn_thres, ts_vals, NNdist_old
        )

        flag = garcia_almeida_break_criterion(FNNs, counter, fnn_thres, max_num_of_cycles)
        counter += 1
        # for the next embedding cycle save the distances of NN for this dimension
        NNdist_old = NNdist_new
    end

    # construct final reconstruction vector
    Y_final = s_orig
    for i = 2:length(τ_vals[1:counter-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final,s_orig,τ_vals[i])
    end

    return Y_final, τ_vals[1:counter-1], ts_vals[1:counter-1], FNNs, NS[:,1:counter-1]

end

function garcia_almeida_embedding(Y::AbstractDataset{D, F}; τs = 0:50 , w::Int = 1,
    r1::Real = 10, r2::Real = 2, fnn_thres::Real = 0.05,
    T::Int = 1, metric = Euclidean(), max_num_of_cycles = 50) where {D, F<:Real}

    @assert 0 ≤ fnn_thres < 1 "Please select a valid breaking criterion, i.e. a threshold value `fnn_thres` ∈ [0 1)"
    @assert all(x -> x ≥ 0, τs)
    Y_orig = Y
    Y = standardize(Y) # especially important for fnn-computation
    # define actual phase space trajectory and NN-distances
    Y_act = []
    NNdist_old = []

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    # preallocate output variables
    τ_vals = Int64[0]
    ts_vals = Int64[1]
    FNNs = Float64[]
    NS = fill(zeros(F, length(τs), size(Y,2)), 1, max_num_of_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        Y_act, NNdist_new = garcia_multivariate_embedding_cycle!(
            Y_act, flag, Y, τs, w, counter, NS, τ_vals, metric, r1, r2, T,
            FNNs, fnn_thres, ts_vals, NNdist_old
        )

        flag = garcia_almeida_break_criterion(FNNs, counter, fnn_thres, max_num_of_cycles)
        counter += 1
        # for the next embedding cycle save the distances of NN for this dimension
        NNdist_old = NNdist_new
    end

    # construct final reconstruction vector
    Y_final = Y_orig[:,ts_vals[1]]
    for i = 2:length(τ_vals[1:counter-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final,Y_orig[:,ts_vals[i]],τ_vals[i])
    end

    return Y_final, τ_vals[1:counter-1], ts_vals[1:counter-1], FNNs, NS[:,1:counter-1]

end

## core: Garcia-Almeida-method for one arbitrary embedding cycle

"""
    n_statistic(Y, s; kwargs...) → N, d_E1
Perform one embedding cycle according to the method proposed in [^Garcia2005a]
for a given phase space trajectory `Y` (of type `Dataset`) and a time series `s
(of type `Vector`). Return the proposed N-Statistic `N` and all nearest
neighbor distances `d_E1` for each point of the input phase space trajectory
`Y`. Note that `Y` is a single time series in case of the first embedding cycle.

## Keyword arguments
* `τs= 0:50`: Considered delay values `τs` (in sampling time units). For each of
  the `τs`'s the N-statistic gets computed.
* `r = 10`: The threshold, which defines the factor of tolerable stretching for
  the d_E1-statistic (see algorithm description).
* `T::Int = 1`: The forward time step (in sampling units) in order to compute the
  `d_E2`-statistic (see algorithm description). Note that in the paper this is
  not a free parameter and always set to `T=1`.
* `w::Int = 0`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors. Note that in the paper this is not
  a free parameter and always `w=0`.
* `metric = Euclidean()`: metric used for finding nearest neigbhors in the input
  phase space trajectory `Y`.

## Description
For a range of possible delay values `τs` one constructs a temporary embedding
matrix. That is, one concatenates the input phase space trajectory `Y` with the
`τ`-lagged input time series `s`. For each point on the temporary trajectory one
computes its nearest neighbor, which is denoted as the `d_E1`-statistic for a
specific `τ`. Now one considers the distance between the reference point and its
nearest neighbor `T` sampling units ahead and calls this statistic `d_E2`.
[^Garcia2005a] strictly use `T=1`, so they forward each reference point and its
corresponding nearest neighbor just by one (!) sampling unit. Here it is a free
parameter.

The `N`-statistic is then the fraction of `d_E2`/`d_E1`-pairs which exceed a
threshold `r`.

Plotted vs. the considered `τs`-values it is proposed to pick the `τ`-value for
this embedding cycle as the value, where `N` has its first local minimum.

[^Garcia2005a]: Garcia, S. P., Almeida, J. S. (2005). [Nearest neighbor embedding with different time delays. Physical Review E 71, 037204](https://doi.org/10.1103/PhysRevE.71.037204).
"""
function n_statistic(Y::AbstractDataset{D, F}, s::Vector{F}; τs = 0:50 , r::Real = 10,
    T::Int = 1, w::Int = 1, metric = Euclidean()) where {D, F<:Real}

    # assert a minimum length of the input time series
    @assert length(s) ≥ length(Y) "The length of the input time series `s` must be at least the length of the input trajectory `Y` "
    @assert all(x -> x ≥ 0, τs)
    τ_max = maximum(τs)

    # preallocation of output
    N_stat = zeros(length(τs))
    NN_distances = [AbstractArray[] for j in τs]

    for (i,τ) in enumerate(τs)
        # build temporary embedding matrix Y_temp
        Y_temp = hcat_lagged_values(Y,s,τ)
        NN = length(Y_temp)     # length of the temporary phase space vector
        N = NN-T               # accepted length w.r.t. the time horizon `T

        vtree = KDTree(Y_temp[1:N], metric) # tree for input data
        # nearest neighbors (d_E1)
        allNNidxs, d_E1 = all_neighbors(vtree, Y_temp[1:N], 1:N, 1, w)
        NN_distances[i] = d_E1

        # for each point, consider its next iterate and compute the distance to
        # its former(!) nearest neighbour iterate
        newfididxs = [j .+ T for j = 1:N]
        newNNidxs = [j .+ T for j in allNNidxs]

        # compute d_E2
        d_E2 = [evaluate(metric,Y_temp[newfididxs[i]],Y_temp[newNNidxs[i][1]]) for i = 1:N]

        # look at ratio d_E2/d_E1
        distance_ratios = d_E2./vcat(d_E1...)

        # number of distance ratio larger than r
        N_stat[i] = sum(distance_ratios.>r)/NN

    end

    return N_stat, NN_distances
end

# Here we separate the inner core loop of the mdop_embedding into another function
# not only for clarity of source code but also to introduce a function barrier
# for the type instability of `Y → Y_act`
function garcia_embedding_cycle!(
    Y_act, flag, s, τs, w, counter, NS, τ_vals, metric, r1, r2, T,
    FNNs, fnn_thres, ts_vals, NNdist_old)

    N, NNdistances = n_statistic(Y_act, s; τs = τs , r = r1,
        T = T, w = w, metric = metric)
    NS[:, counter] = N

    # determine the optimal tau value from the N-statistic
    min_idx = findlocalminima(N)

    # store chosen delay (and chosen time series)
    push!(τ_vals, τs[min_idx[1]])
    push!(ts_vals, 1)

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act,s,τ_vals[counter+1])
    vtree = KDTree(Y_act, metric)
    allNNidxs, NNdist_new = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # compute FNN-statistic and store vals, while also
    # truncating distance-vector to the "new" length of the actual embedding vector
    fnns = fnn_embedding_cycle(view(NNdist_old, 1:length(Y_act)), NNdist_new, r2)
    push!(FNNs,fnns)
    return Y_act, NNdist_new
end

function garcia_multivariate_embedding_cycle!(
    Y_act, flag, Ys, τs, w, counter, NS, τ_vals, metric, r1, r2, T,
    FNNs, fnn_thres, ts_vals, NNdist_old
    )

    M = size(Ys,2)
    # in the 1st cycle we have to check all (size(Y,2)^2 combinations
    if counter == 1
        Y_act, NNdist_new = first_embedding_cycle!(M, Ys, τs, r1, r2, T, w,
                                                metric, τ_vals, ts_vals,
                                                NS, FNNs)

    # in all other cycles we just have to check (size(Y,2)) combinations
    else
        Y_act, NNdist_new = embedding_cycle!(Y_act, M, counter, Ys, τs, r1, r2, T, w,
                                            metric, NNdist_old, τ_vals, ts_vals,
                                            NS, FNNs)
    end
    return Y_act, NNdist_new
end

"""
Perform the first embedding cycle of the multivariate embedding. Return the
actual reconstruction vector `Y_act` and all nearest neighbor distances
`NNdist_new` of `Y_act`.
"""
function first_embedding_cycle!(M, Ys, τs, r1, r2, T, w, metric, τ_vals, ts_vals,
                                NS, FNNs)
    counter = 1
    Ns = zeros(length(τs), M*M)
    n_counter = 1
    for ts = 1:M
        for ts2 = 1:M
            Ns[:,n_counter], _ = n_statistic(Dataset(Ys[:,ts]), Ys[:,ts2];
                                            τs = τs , r = r1, T = T, w = w,
                                            metric = metric)
            n_counter += 1
        end
    end
    # determine the optimal tau value from the N-statistic
    τ_opt, ts_opt1, ts_opt2 = choose_optimal_tau1_garcia(Ns, M)
    # store chosen delay, time series and according N-statistic
    push!(τ_vals, τs[τ_opt])
    push!(ts_vals, ts_opt1)
    push!(ts_vals, ts_opt2)
    NS[counter] = Ns[:,(ts_vals[counter]-1)*M+1:(ts_vals[counter]-1)*M+M]

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Ys[:,ts_vals[counter]],
                                Ys[:,ts_vals[counter+1]],τ_vals[counter+1])

    # compute nearest neighbor distances in the first embedding dimension for
    # FNN statistic
    Y_old = Dataset(Ys[:,ts_vals[counter]])
    vtree = KDTree(Y_old, metric)
    _, NNdist_old = all_neighbors(vtree, Y_old, 1:length(Y_old), 1, w)

    vtree = KDTree(Y_act, metric)
    _, NNdist_new = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # compute FNN-statistic and store vals, while also
    # truncating distance-vector to the "new" length of the actual embedding vector
    fnns = fnn_embedding_cycle(view(NNdist_old, 1:length(Y_act)), NNdist_new, r2)
    push!(FNNs,fnns)

    return Y_act, NNdist_new
end

"""
Perform an embedding cycle of the multivariate embedding, but the 1st. Return the
actual reconstruction vector `Y_act` and all nearest neighbor distances
`NNdist_new` of `Y_act`.
"""
function embedding_cycle!(Y_act, M, counter, Ys, τs, r1, r2, T, w, metric, NNdist_old,
                                            τ_vals, ts_vals, NS, FNNs)
    Ns = zeros(length(τs), M)
    for ts = 1:M
        Ns[:,ts], _ = n_statistic(Y_act, Ys[:,ts]; τs = τs , r = r1, T = T,
                                    w = w, metric = metric)
    end

    # determine the optimal tau value from the N-statistic
    τ_opt, ts_opt = choose_optimal_tau2_garcia(Ns)
    # store chosen delay, time series and according n-statistic
    push!(τ_vals, τs[τ_opt])
    push!(ts_vals, ts_opt)
    NS[counter] = Ns

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act,
                                Ys[:,ts_vals[counter+1]],τ_vals[counter+1])
    vtree = KDTree(Y_act, metric)
    _, NNdist_new = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # compute FNN-statistic and store vals, while also
    # truncating distance-vector to the "new" length of the actual embedding vector
    fnns = fnn_embedding_cycle(view(NNdist_old, 1:length(Y_act)), NNdist_new, r2)
    push!(FNNs,fnns)

    return Y_act, NNdist_new
end

"""
This function chooses the optimal τ value as the 1st minimum of all the
N-statistics 1st minima. It returns the index of this 1st minimum as well as the
time series number corresponding to that 1st minimum of all 1st minima. In this
first embedding cycle it also returns the time series number, which act as Y_act.
"""
function choose_optimal_tau1_garcia(Ns::Array{T, 2}, M::Int) where {T}
    NN = size(Ns,2)
    min_idx = zeros(Int, NN)
    for i = 1:NN
        # determine optimal tau value and minimum vals from all N-statistic's
        minimas = findlocalminima(Ns[:,i])
        min_idx[i] = minimas[1]
    end
    # take the 1st of all 1st minimas
    _, idx = findmin(min_idx)
    # determine the first time series and the according second one
    a = idx/M
    ts_1 = Int(ceil(a))
    ts_2 = idx-M*(ts_1-1)
    return min_idx[idx], ts_1, ts_2
end

"""
This function chooses the optimal τ value as the 1st minimum of the 1st minima
of the N-statistics. It returns the index of this 1st minimum as well as the
time series number corresponding to that 1st minimum of the 1st minima.
"""
function choose_optimal_tau2_garcia(Ns::Array{T, 2}) where {T}
    NN = size(Ns,2)
    min_idx = zeros(Int, NN)
    for i = 1:NN
        # determine optimal tau value and minimum vals from all N-statistic's
        minimas = findlocalminima(Ns[:,i])
        min_idx[i] = minimas[1]
    end
    # take the 1st of all 1st minimas
    _, idx = findmin(min_idx)
    return min_idx[idx], idx
end

function garcia_almeida_break_criterion(FNNs, counter, fnn_thres, max_num_of_cycles)
    flag = true
    if FNNs[counter] ≤ fnn_thres
        flag = false
        println("Algorithm stopped due to sufficiently small FNNs. "*
                "Valid embedding achieved ✓.")
    end
    if counter > 1 && FNNs[counter] > FNNs[counter-1]
        flag = false
        println("Algorithm stopped due to rising FNNs. "*
                "Valid embedding achieved ✓.")
    end
    if max_num_of_cycles == counter
        println("Algorithm stopped due to hitting max cycle number. "*
                "Valid embedding NOT achieved ⨉.")
        flag = false
    end
    return flag
end


text:
using Distributions, NumericExtensions

module GLM

    using Distributions, NumericExtensions, NumericFuns
    using Base.LinAlg.LAPACK: potrf!, potrs!
    using Base.LinAlg.BLAS: gemm!, gemv!
    using Base.LinAlg: QRCompactWY, Cholesky
    using StatsBase: StatsBase, CoefTable, StatisticalModel, RegressionModel
    using Distributions: sqrt2, sqrt2π

    import Base: (\), cholfact, cor, scale, show, size
    import StatsBase: coef, coeftable, confint, deviance, loglikelihood, nobs, stderr,
                      vcov, residuals, predict, fit
    import NumericExtensions: evaluate, result_type

    export                              # types
        CauchitLink,
        CloglogLink,
        DensePred,
        DensePredQR,
        DensePredChol,
        GeneralizedLinearModel,
        GlmResp,
        IdentityLink,
        InverseLink,
        LinearModel,
        Link,
        LinPred,
        LinPredModel,
        LogitLink,
        LogLink,
        LmResp,
        ProbitLink,
        SqrtLink,

        # Deprecated
        LmMod,
        GlmMod,

                                        # functions
        canonicallink,  # canonical link function for a distribution
        contr_treatment,# treatment contrasts
        delbeta!,       # evaluate the increment in the coefficient vector
        deviance,       # deviance of fitted and observed responses
        devresid,       # vector of squared deviance residuals
        df_residual,    # degrees of freedom for residuals
        drsum,          # sum of squared deviance residuals
        fit,            # function to fit models, from StatsBase
        formula,        # extract the formula from a model
        glm,            # general interface
        linkfun!,       # mutating link function
        linkfun,        # link function mapping mu to eta, the linear predictor
        linkinv!,       # mutating inverse link
        linkinv,        # inverse link mapping eta to mu
        linpred,        # linear predictor
        linpred!,       # update the linear predictor
        lm,             # linear model (QR factorization)
        lmc,            # linear model (Cholesky factorization)
        mueta!,         # mutating derivative of inverse link
        mueta,          # derivative of inverse link
        mustart,        # derive starting values for the mu vector
        nobs,           # total number of observations
        objective,      # the objective function in fitting a model
        predict,        # make predictions
        sqrtwrkwt,      # square root of the working weights
        stderr,         # standard errors of the coefficients
        updatemu!,      # update the response type from the linear predictor
        var!,           # mutating variance function
        wrkresid!,      # mutating working residuals function
        wrkresid,       # extract the working residuals
        wrkresp         # working response

    typealias FP FloatingPoint
    typealias FPVector{T<:FloatingPoint} DenseArray{T,1}

    abstract ModResp                   # model response

    abstract LinPred                   # linear predictor in statistical models
    abstract DensePred <: LinPred      # linear predictor with dense X
    abstract LinPredModel <: RegressionModel # model based on a linear predictor

    include("linpred.jl")
    include("lm.jl")
    include("glmtools.jl")
    include("glmfit.jl")
    include("deprecated.jl")

end # module


text:
for (A,F) in ((:gamma, :arb_gamma),
              (:lgamma, :arb_lgamma),(:rgamma, :arb_rgamma), (:digamma, :arb_digamma),
              (:zeta, :arb_zeta),

              (:erf, :arb_hypgeom_erf), (:erfc, :arb_hypgeom_erfc), (:erfi, :arb_hypgeom_erfi),
              (:ei, :arb_hypgeom_ei), (:si, :arb_hypgeom_si), (:ci, :arb_hypgeom_ci),
              (:shi, :arb_hypgeom_shi), (:chi, :arb_hypgeom_chi),
             )
    @eval begin
        function ($A)(x::ArbReal{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, prec)
            return z
         end
    end
end

const Cint0 = zero(Cint)
const FmpzZero = Fmpz(0)

function lambertw(x::ArbReal{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_lambertw), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Cint, Clong), z, x, Cint0, P)
    return z
end

lambertw(x::ArbFloat{P}) where {P} = ArbFloat{P}(lambertw(ArbReal{P}(x)))

function lambertw(x::ArbComplex{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_lambertw), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{Fmpz}, Cint, Clong), z, x, FmpzZero, 0, 128)
    return z
end

for (A,F) in ((:gamma, :acb_gamma),
              (:lgamma, :acb_lgamma),(:rgamma, :acb_rgamma), (:digamma, :acb_digamma),

              (:logsinpi, :acb_log_sin_pi), (:barnesg, :acb_barnes_g), (:logbarnesg, :acb_log_barnes_g),
              (:agm1, :acb_agm1),

              (:erf, :acb_hypgeom_erf), (:erfc, :acb_hypgeom_erfc), (:erfi, :acb_hypgeom_erfi),
              (:ei, :acb_hypgeom_ei), (:si, :acb_hypgeom_si), (:ci, :acb_hypgeom_ci),
              (:shi, :acb_hypgeom_shi), (:chi, :acb_hypgeom_chi),
             )
    @eval begin
        function ($A)(x::ArbComplex{P}, prec::Int=P) where P
            z = ArbComplex{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, x, prec)
            return z
         end
    end
end

for (A,F) in ((:elliptick, :acb_elliptic_k), (:elliptice, :acb_elliptic_e),
              (:eta, :acb_dirichlet_eta), (:xi, :acb_dirichlet_xi)
             )
    @eval begin
        function ($A)(x::ArbReal{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            xc = ArbComplex{P}(x)
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, xc, prec)
            rea = real(z)
            return rea
         end
    end
end
#=
for (A,F) in ((:elliptick, :elliptic_k), (:elliptice, :elliptic_e))
    @eval begin
        function ($A)(k::ArbReal{P}, prec::Int=P) where {P}
           m = square(k)
           return ($F)(m, prec)
        end
    end
end
=#
for (A,F) in ((:agm, :arb_agm), )
    @eval begin
        function ($A)(x::ArbReal{P}, y::ArbReal{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, y, prec)
            return z
         end
    end
end

for (A,F) in ((:ellipticp, :acb_elliptic_p), (:ellipticpi, :acb_elliptic_pi),
              (:ellipticzeta, :acb_elliptic_zeta), (:ellipticsigma, :acb_elliptic_sigma),
             )
    @eval begin
        function ($A)(x::ArbComplex{P}, y::ArbComplex{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Ref{ArbReal}, Clong), z, x, y, prec)
            return z
         end
    end
end

for (A,F) in ((:ellipticpi, :acb_elliptic_pi), (:ellipticp, :acb_elliptic_p),
              (:ellipticzeta, :acb_elliptic_zeta), (:ellipticsigma, :acb_elliptic_sigma),
             )
    @eval begin
        function ($A)(x::ArbReal{P}, y::ArbReal{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            xc = ArbComplex{P}(x)
            yc = ArbComplex{P}(y)
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, xc, yc, prec)
            rea = real(z)
            return rea
         end
    end
end
#=
function ellipticpi(n::ArbComplex{P}, k::ArbComplex{P}, prec::Int=P) where {P}
    m = square(k)
    return elliptic_pi(n, m, prec)
end
=#

for (A,F) in ((:gamma, :arb_gamma),
              (:lgamma, :arb_lgamma),(:rgamma, :arb_rgamma), (:digamma, :arb_digamma),
              (:zeta, :arb_zeta),

              (:erf, :arb_hypgeom_erf), (:erfc, :arb_hypgeom_erfc), (:erfi, :arb_hypgeom_erfi),
              (:ei, :arb_hypgeom_ei), (:si, :arb_hypgeom_si), (:ci, :arb_hypgeom_ci),
              (:shi, :arb_hypgeom_shi), (:chi, :arb_hypgeom_chi),
             )
    @eval begin
        function ($A)(x::ArbFloat{P}, prec::Int=P) where {P}
            z = ArbReal{P}()
            xb = ArbReal{P}(x)
            ccall(@libarb($F), Cvoid, (Ref{ArbReal}, Ref{ArbReal}, Clong), z, xb, prec)
            return midpoint_byref(z)
         end
    end
end

for (A,F) in ((:elliptick, :acb_elliptic_k), (:elliptice, :acb_elliptic_e),
              (:eta, :acb_dirichlet_eta), (:xi, :acb_dirichlet_xi))
    @eval begin
        function ($A)(x::ArbFloat{P}, prec::Int=P) where {P}
            z = ArbComplex{P}()
            xc = ArbComplex{P}(x)
            ccall(@libarb($F), Cvoid, (Ref{ArbComplex}, Ref{ArbComplex}, Clong), z, xc, prec)
            return midpoint_byref(real(z))
         end
    end
end

function agm1(x::ArbReal{P}) where {P}
    x1 = ArbComplex{P}(x)
    w  = agm1(x1)
    z = real(w)
    return z
end

function agm1(x::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    w  = agm1(x1)
    z = midpoint_byref(w)
    return z
end

function agm(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
    x1 = ArbReal{P}(x)
    y1 = ArbReal{P}(y)
    w  = agm(x1, y1)
    z = midpoint_byref(w)
    return z
end


text:
using Test, OrderedCollections
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

@testset "Average updating method" begin
    x_new = rand(3)
    x_old = rand(3)
    learning_rates = rand(10)
    for i in 1:10
        x_avg = average_update(x_new, x_old, learning_rates[i])
        @test x_avg == (x_new * learning_rates[i] + x_old * (1 - learning_rates[i]))
    end
end

@testset "Methods for calculating errors during functional iteration" begin
    funcvar1 = OrderedDict(:a => rand(3), :b => rand(3))
    funcvar2 = OrderedDict(:a => rand(3), :b => rand(3))

    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :total_error) == sum(abs.(funcvar1[:a] - funcvar2[:a]))
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :L∞) == maximum(abs.(funcvar1[:a] - funcvar2[:a]))
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :Linf) == maximum(abs.(funcvar1[:a] - funcvar2[:a]))
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :max_abs_error) == maximum(abs.(funcvar1[:a] - funcvar2[:a]))
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :maximum_absolute_error) == maximum(abs.(funcvar1[:a] - funcvar2[:a]))
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :L²) == sum((funcvar1[:a] - funcvar2[:a]) .^ 2)
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :L2) == sum((funcvar1[:a] - funcvar2[:a]) .^ 2)
    @test calculate_func_error(funcvar1[:a], funcvar2[:a], :squared_error) == sum((funcvar1[:a] - funcvar2[:a]) .^ 2)

    @test calculate_func_error(funcvar1, funcvar2, :total_error) == sum(map(x -> sum(abs.(funcvar1[x] - funcvar2[x])), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :L∞) ==
        maximum(map(x -> maximum(abs.(funcvar1[x] - funcvar2[x])), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :Linf) ==
        maximum(map(x -> maximum(abs.(funcvar1[x] - funcvar2[x])), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :max_abs_error) ==
        maximum(map(x -> maximum(abs.(funcvar1[x] - funcvar2[x])), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :maximum_absolute_error) ==
        maximum(map(x -> maximum(abs.(funcvar1[x] - funcvar2[x])), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :L²) == sum(map(x -> sum((funcvar1[x] - funcvar2[x]) .^ 2), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :L2) == sum(map(x -> sum((funcvar1[x] - funcvar2[x]) .^ 2), [:a, :b]))
    @test calculate_func_error(funcvar1, funcvar2, :squared_error) == sum(map(x -> sum((funcvar1[x] - funcvar2[x]) .^ 2), [:a, :b]))

    @test calculate_func_error(funcvar1, funcvar2, :total_error; vars = [:a]) == sum(map(x -> sum(abs.(funcvar1[x] - funcvar2[x])), [:a]))
    @test calculate_func_error(funcvar1, funcvar2, :L∞; vars = [:a]) ==
        maximum(map(x -> maximum(abs.(funcvar1[x] - funcvar2[x])), [:a]))
    @test calculate_func_error(funcvar1, funcvar2, :L²; vars = [:a]) == sum(map(x -> sum((funcvar1[x] - funcvar2[x]) .^ 2), [:a]))
end

@testset "Construction of non-uniform grid spacing"  begin
    stategrid = StateGrid(OrderedDict(:x => sort(rand(10))))
    bc = [(0., 1.)]
    @test_throws AssertionError nonuniform_grid_spacing(stategrid, [(0., 1.), (0., 1.)])
    @test nonuniform_grid_spacing(stategrid) ≈ vcat(diff(stategrid[:x][1:2]), diff(stategrid[:x]), diff(stategrid[:x][end - 1:end]))
    @test nonuniform_grid_spacing(stategrid) ≈ vcat(stategrid[:x][1], diff(stategrid[:x]), 1. - stategrid[:x][end])
end


text:
using TimeSeriesClustering

ts_input_data = load_timeseries_data(:CEP_GER1)

# when using kmedoids-exact, one needs to supply the optimizer. Make sure the optimizer is added through Pkg.add()

using Cbc
optimizer = Cbc.Optimizer
out = run_clust(ts_input_data;method="kmedoids_exact",representation="medoid",n_clust=ts_input_data.5,n_init=1,kmexact_optimizer=optimizer)

using Gurobi
optimizer = Gurobi.Optimizer
out = run_clust(ts_input_data;method="kmedoids_exact",representation="medoid",n_init=1,kmexact_optimizer=optimizer)