text:
# Autogenerated wrapper script for Libical_jll for x86_64-apple-darwin
export libical, libical_cxx, libicalss, libicalss_cxx, libicalvcal

using BerkeleyDB_jll
using Glib_jll
using ICU_jll
using XML2_jll
JLLWrappers.@generate_wrapper_header("Libical")
JLLWrappers.@declare_library_product(libical, "@rpath/libical.3.dylib")
JLLWrappers.@declare_library_product(libical_cxx, "@rpath/libical_cxx.3.dylib")
JLLWrappers.@declare_library_product(libicalss, "@rpath/libicalss.3.dylib")
JLLWrappers.@declare_library_product(libicalss_cxx, "@rpath/libicalss_cxx.3.dylib")
JLLWrappers.@declare_library_product(libicalvcal, "@rpath/libicalvcal.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header(BerkeleyDB_jll, Glib_jll, ICU_jll, XML2_jll)
    JLLWrappers.@init_library_product(
        libical,
        "lib/libical.3.0.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libical_cxx,
        "lib/libical_cxx.3.0.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicalss,
        "lib/libicalss.3.0.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicalss_cxx,
        "lib/libicalss_cxx.3.0.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicalvcal,
        "lib/libicalvcal.3.0.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module IntervalLapper

export Interval, overlap, intersect
export Lapper, find, seek, lower_bound, merge_overlaps!, coverage, union_and_intersect
export Bits, count

"""
Represents an interval of (start, stop] with a label val.
"""
struct Interval{T}
	start::Int
	stop::Int
	val::T
end

"""
Check if two intervals overlap.

# Example
```jldoctest
julia> overlap(Interval(0, 10, "cats"), Interaval(5, 15, "Dogs"))
true
```
"""
function overlap(a::Interval, b::Interval)
	a.start < b.stop && a.stop > b.start
end

function overlap(a::Interval, start::Int, stop::Int)
	a.start < stop && a.stop > start
end

"""
Determine the number of positions shared by two intervals

# Example
```jldoctest
julia> intersect(Interval(0, 10, "cats"), Interaval(5, 15, "Dogs"))
5
```
"""
function intersectlen(a::Interval, b::Interval)
	diff = min(a.stop, b.stop) - max(a.start, b.start)
	diff >= 0 ? diff : 0
end

"""
Primary object of the library. The intervals can be used for iterating / pulling values out of the tree
"""
mutable struct Lapper{T}
	intervals::Vector{Interval{T}}
	max_len::Int
	cursor::Int
	cov::Union{Nothing, Int}
	overlaps_merged::Bool
end

function Lapper(intervals::Vector{Interval{T}}) where T
	sort!(intervals, by = x -> (x.start))
	max_len = 0
	for interval in intervals
		iv_len = interval.stop - interval.start # add an check for intervals where this could be negative?
		if iv_len > max_len
			max_len = iv_len
		end
	end
	return Lapper(intervals, max_len, 0, nothing, false)
end

function lower_bound(start::Int, intervals::Vector{Interval{T}}) where T
	size = length(intervals)
	low = 1

	@inbounds while size > 1
		half = div(size, 2)
		other_half = size - half
		probe = low + half
		other_low = low + other_half
		v = intervals[probe]
		size = half
		low = v.start < start ? other_low : low
	end

	low
end

"""
Merge any intervals that overlap with eachother within the Lapper. This is an easy
way to speed up queries.
"""
function merge_overlaps!(lapper::Lapper{T}) where T
	if length(lapper.intervals) <= 1
		lapper.overlaps_merged = true
		return
	elseif lapper.overlaps_merged
		return
	end

	stack = Vector{Interval{T}}()
	first = lapper.intervals[1]
	push!(stack, first)
	for interval in lapper.intervals[2:end]
		top = pop!(stack)
		if top.stop < interval.start
			push!(stack, top)
			push!(stack, interval)
		elseif top.stop < interval.stop
			top = Interval{T}(top.start, interval.stop, interval.val)
			push!(stack, top)
		else
			# they were equal
			push!(stack, top)
		end
	end
	lapper.overlaps_merged = true
	lapper.intervals = stack
end

"""
Calculate the nuber of positions covered by the intervals in Lapper.
"""
function coverage(lapper::Lapper{T}) where T
	moving_start = 0
	moving_stop = 0
	cov = 0

	for interval in lapper.intervals
		if overlap(interval, moving_start, moving_stop)
			moving_start = min(moving_start, interval.start)
			moving_stop = max(moving_stop, interval.stop)
		else
			cov += moving_stop - moving_start
			moving_start = interval.start
			moving_stop = interval.stop
		end
	end
	cov += moving_stop - moving_start
	cov
end

"""
Find the union and the intersect of two lapper objects.
Union: The set of positions found in both lappers
Intersect: The number of positions where both lappers intersect. Note that a position only
ounts one time, multiple Intervals covering the same position don't add up.
TODO: make this and other funcitons more generic an not depend on T if they don't have to
"""
function union_and_intersect(self::Lapper{T}, other::Lapper{T}, self_cov::Union{Nothing, Int}=nothing, other_cov::Union{Nothing, Int}=nothing) where T
	cursor = Ref(1)
	if !self.overlaps_merged || !other.overlaps_merged
		intersections = Vector{Interval{Bool}}()
		for self_iv in self.intervals
			for other_iv in seek(other, self_iv.start, self_iv.stop, cursor)
				start = max(self_iv.start, other_iv.start)
				stop = min(self_iv.stop, other_iv.stop)
				push!(intersections, Interval(start, stop, true))
			end
		end
		temp_lapper = Lapper(intersections)
		merge_overlaps!(temp_lapper)
		temp_cov = coverage(temp_lapper)
		other_cov = isnothing(other_cov) ? coverage(other) : other_cov
		self_cov = isnothing(self_cov) ? coverage(self) : self_cov
		union = self_cov + other_cov - temp_cov
		return union, temp_cov
	else
		intersect = 0
		for c1_iv in self.intervals
			for c2_iv in seek(other, c1_iv.start, c1_iv.stop, cursor)
				local_intersect = intersectlen(c1_iv, c2_iv)
				intersect += local_intersect
			end
		end
		other_cov = isnothing(other_cov) ? coverage(other) : other_cov
		self_cov = isnothing(self_cov) ? coverage(self) : self_cov
		union =  self_cov + other_cov - intersect
		return (union, intersect)
	end
end

#=
#
# Find Iterator / Seek Iterator
#
=#

@inline function checked_sub(a::Int, b::Int, or=1)
	maybe = a - b
	maybe >= 1 ? maybe : or
end

struct FindIter{T}
	inner::Lapper{T}
	start::Int
	stop::Int
end

struct SeekIter{T}
	inner::Lapper{T}
	start::Int
	stop::Int
	cursor::Ref{Int}
end

@inline function _find(iter::Union{FindIter, SeekIter}, offset::Int)
    while offset <= length(iter.inner.intervals)
	interval = iter.inner.intervals[offset]
	offset += 1
        if overlap(interval, iter.start, iter.stop)
	    return (interval, offset)
        elseif interval.start >= iter.stop
            break
        end
    end
    nothing
end

find(lapper::Lapper{T}, start::Int, stop::Int) where T = FindIter(lapper, start, stop)
Base.iterate(iter::FindIter, offset=lower_bound(checked_sub(iter.start, iter.inner.max_len), iter.inner.intervals)) = _find(iter, offset)
Base.IteratorSize(::FindIter) = Base.SizeUnknown()

function seek(lapper::Lapper{T}, start::Int, stop::Int, cursor::Ref{Int}) where T
	if cursor[] <= 1 || (cursor[] <= length(lapper.intervals) && lapper.intervals[cursor[]].start > start)
		cursor[] = lower_bound(checked_sub(start, lapper.max_len), lapper.intervals)
	end

	while cursor[] + 1 <= length(lapper.intervals) && lapper.intervals[cursor[] + 1].start < checked_sub(start, lapper.max_len)
		cursor[] += 1
	end
	SeekIter(lapper, start, stop, cursor)
end

Base.iterate(iter::SeekIter, offset=iter.cursor[]) = _find(iter, offset)
Base.IteratorSize(::SeekIter) = Base.SizeUnknown()

#=
# Depth Iterator
=#

struct DepthIter
	inner::Lapper
	# Lapper that is merged lapper of the inner
	merged::Lapper{Bool}
	merged_len::Int
end

"""
Return the contiguous intervals of coverage, `val` represents the number of intervals
covering the returned interval.
"""
function depth(lapper::Lapper)
	merged_lapper = Lapper(collect(map(x -> Interval(x.start, x.stop, true), lapper.intervals)))
	merge_overlaps!(merged_lapper)
	merged_len = length(merged_lapper.intervals)
	DepthIter(lapper, merged_lapper, merged_len)
end
Base.IteratorSize(::DepthIter) = Base.SizeUnknown()

function Base.iterate(iter::DepthIter, (curr_merged_pos, curr_pos, cursor)=(1, 1, Ref(1)))
	interval = iter.merged.intervals[curr_pos]
	if curr_merged_pos == 1
		curr_merged_pos = interval.start
	end
	if interval.stop == curr_merged_pos
		if curr_pos + 1 <= iter.merged_len
			curr_pos += 1
			interval = iter.merged.intervals[curr_pos]
			curr_merged_pos = interval.start
		else
			return nothing
		end
	end
	start = curr_merged_pos
	depth_at_point = 0
	for _ in seek(iter.inner, curr_merged_pos, curr_merged_pos + 1, cursor)
		depth_at_point += 1
	end
	new_depth_at_point = depth_at_point
	while new_depth_at_point == depth_at_point && curr_merged_pos < interval.stop
		curr_merged_pos += 1

		tmp = 0
		for _ in seek(iter.inner, curr_merged_pos, curr_merged_pos + 1, cursor)
			tmp += 1
		end
		new_depth_at_point = tmp
	end
	return (Interval(start, curr_merged_pos, depth_at_point), (curr_merged_pos, curr_pos, cursor))
end

"""
A data structure for counting all intervals that overlap start .. stop. It is very fast.
Two binary searches are performed to fina all the excluded elements, then the intersections
can be deduced from there. See [BITS](https://arxiv.org/pdf/1208.3407.pdf) for more info.
"""
struct Bits
	starts::Vector{Int}
	stops::Vector{Int}
end

@inline unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))

function Bits(intervals::Vector{Interval{T}}) where T
	starts, stops = unzip(map( x -> (x.start, x.stop), intervals))
	Bits(sort!(starts), sort!(stops))
end

@inline function bsearch_seq(key::Int, elems::Vector{Int})
    if elems[1] > key
        return 1
    end

    high = length(elems) + 1
    low = 1

    while high - low > 1
        mid = div(high + low, 2)
        if elems[mid] < key
            low = mid
        else
            high = mid
        end
    end
    high
end

function count(bits::Bits, start::Int, stop::Int)
	len = length(bits.starts)
	first = bsearch_seq(start, bits.stops)
	last = bsearch_seq(stop, bits.starts)
	while first <= len && bits.stops[first] == start
		first += 1
	end

	num_cant_after = len - last
	len - first - num_cant_after
end

end # module

ext: jl


text:
"""
TrackedFiles

Convenience type to keep track of files to watch.
"""
const TrackedFiles = Dict{Pair{String, String}, Float64}

"""
$(SIGNATURES)

Prepare the output directory `PATHS[:pub]`.

* `clear=true` removes the content of the output directory if it exists to start from a blank
slate
"""
function prepare_output_dir(clear::Bool=true)::Nothing
    # if required to start from a blank slate -> remove the output dir
    (clear & isdir(PATHS[:pub])) && rm(PATHS[:pub], recursive=true)
    # create the output dir and the css dir if necessary
    !isdir(PATHS[:pub]) && mkdir(PATHS[:pub])
    !isdir(PATHS[:css]) && mkdir(PATHS[:css])
    return nothing
end

"""
$(SIGNATURES)

Take a `root` path to an input file and convert to output path. If the output path does not exist,
create it.
"""
function out_path(root::String)::String
    len_in = lastindex(joinpath(PATHS[:src], ""))
    length(root) <= len_in && return PATHS[:folder]
    dpath = root[nextind(root, len_in):end]
    # construct the out path
    f_out_path = joinpath(PATHS[:folder], dpath)
    f_out_path = replace(f_out_path, r"([^a-zA-Z\d\s_:])pages" => s"\1pub")
    # if it doesn't exist, make the path
    !ispath(f_out_path) && mkpath(f_out_path)
    return f_out_path
end

"""
$(SIGNATURES)

Update the dictionaries referring to input files and their time of last change. The variable `verb`
propagates verbosity.
"""
function scan_input_dir!(md_files::TrackedFiles, html_files::TrackedFiles,
                         other_files::TrackedFiles, infra_files::TrackedFiles,
                         literate_files::TrackedFiles, verb::Bool=false)::Nothing
    # top level files (src/*)
    for file ∈ readdir(PATHS[:src])
        isfile(joinpath(PATHS[:src], file)) || continue
        # skip if it has to be ignored
        file ∈ IGNORE_FILES && continue
        fname, fext = splitext(file)
        fpair = (PATHS[:src] => file)
        if file == "config.md"
            add_if_new_file!(infra_files, fpair, verb)
        elseif fext == ".md"
            add_if_new_file!(md_files, fpair, verb)
        else
            add_if_new_file!(html_files, fpair, verb)
        end
    end
    # pages files (src/pages/*)
    for (root, _, files) ∈ walkdir(PATHS[:src_pages])
        for file ∈ files
            isfile(joinpath(root, file)) || continue
            # skip if it has to be ignored
            file ∈ IGNORE_FILES && continue
            fname, fext = splitext(file)
            fpair = (root => file)
            if fext == ".md"
                add_if_new_file!(md_files, fpair, verb)
            elseif fext == ".html"
                add_if_new_file!(html_files, fpair, verb)
            else
                add_if_new_file!(other_files, fpair, verb)
            end
        end
    end
    # infastructure files (src/_css/* and src/_html_parts/*)
    for d ∈ (:src_css, :src_html), (root, _, files) ∈ walkdir(PATHS[d])
        for file ∈ files
            isfile(joinpath(root, file)) || continue
            fname, fext = splitext(file)
            # skipping files that are not of the type INFRA_FILES
            fext ∉ INFRA_FILES && continue
            add_if_new_file!(infra_files, root=>file, verb)
        end
    end
    # literate script files if any, note that the folder may not exist
    if isdir(PATHS[:literate])
        for (root, _, files) ∈ walkdir(PATHS[:literate])
            for file ∈ files
                isfile(joinpath(root, file)) || continue
                fname, fext = splitext(file)
                # skipping files that are not script file
                fext != ".jl" && continue
                add_if_new_file!(literate_files, root=>file, verb)
            end
        end
    end
    return nothing
end

"""
$(SIGNATURES)

Helper function, if `fpair` is not referenced in the dictionary (new file) add the entry to the
dictionary with the time of last modification as val.
"""
function add_if_new_file!(dict::TrackedFiles, fpair::Pair{String,String}, verb::Bool)::Nothing
    haskey(dict, fpair) && return nothing
    # it's a new file
    verb && println("tracking new file '$(fpair.second)'.")
    dict[fpair] = mtime(joinpath(fpair...))
    return nothing
end


text:
# --- User Parameters

# Number of samples to use to estimate probability distributions
sample_size = 1000;

# --- Imports

import Distributions
using Distributions: MvNormal, ContinuousUnivariateDistribution, cdf
using HypothesisTests
using LinearAlgebra
using Plots
using Random
using StatsBase: Histogram, fit

# --- Internal Parameters

# ------ Uniformity test parameters

# p-value cutoff
p_value_cutoff = 0.1;

# maximum number of attempts to generate a sample of vectors
max_attempts = 10;

# --- Probability distribution for θ (angle between vector and x_1-axis)

struct ThetaDistribution <: ContinuousUnivariateDistribution
    n::Integer
end

function Distributions.cdf(dist::ThetaDistribution, x::Real)
    # Preparations
    k = dist.n-2

    # Compute pdf
    value = 0
    if k % 2 == 0
        # Precompute re-used coefficients
        sgn = (-1)^(k÷2)
        k_choose_k_over_2 = binomial(big(k), k÷2)

        for j = 0:(k÷2 - 1)
            value += (-1)^j * binomial(big(k), j) * sin((k - 2*j) * x) / (k - 2*j)
        end
        value *= 2
        value += sgn * k_choose_k_over_2 * x  # contribution from i = k÷2 term
        value *= sgn / π / k_choose_k_over_2

    else
        normalization = 0
        for j = 0:(k-1)÷2

            coef = (-1)^j * binomial(big(k), j) * (1 / (k - 2*j) - 1 / (2*j - k))
            value += coef * (1 - cos((k - 2*j) * x))
            normalization += coef
        end
        normalization = 0.5 / normalization
        value *= normalization

    end

    return value
end

# --- Tests

# 2D: cdf(x) = x / π
theta_dist_2d = ThetaDistribution(2)
for x = range(0, π; length=10)
    computed_cdf = cdf(theta_dist_2d, x)
    expected_cdf = x / π
    @assert (computed_cdf ≈ expected_cdf) "Expected: $(expected_cdf). Got: $(computed_cdf)"
end

# 3D: cdf(x) = 0.5 * (1 - cos(x))
theta_dist_3d = ThetaDistribution(3)
for x = range(0, π; length=10)
    computed_cdf = cdf(theta_dist_3d, x)
    expected_cdf = 0.5 * (1 - cos(x))
    @assert (computed_cdf ≈ expected_cdf) "Expected: $(expected_cdf). Got: $(computed_cdf)"
end

function generate_vectors(n::Integer, sample_size::Integer)::Tuple
    """
    Generate a sample of unit vectors in an n-dimensional space drawn from a uniform
    distribution over the surface of an unit n-sphere.

    Parameters
    ----------
    n: dimension of space

    sample_size: number of blades to include in sample

    Return values
    -------------
    vectors: sample of vectors

    uniformity_stats: results of tests for uniformity of sample on unit n-sphere
    """
    # --- Set up uniformity tests

    theta_dist = ThetaDistribution(n)
    AD_p_value = 0
    KS_p_value = 0

    # --- Generate sample of vectors

    num_attempts = 0
    vectors = undef
    while ( ((AD_p_value < p_value_cutoff) || (KS_p_value < p_value_cutoff)) &&
            (num_attempts < max_attempts) )

        # Generate new sample
        dist = MvNormal(zeros(n), ones(n))
        num_vectors = sample_size
        vectors = rand(dist, num_vectors)
        for i = 1:num_vectors
            vectors[:, i] /= norm(vectors[:, i])
        end

        # Perform test for uniformity on unit n-sphere
        thetas = map(i -> atan(norm(vectors[:, i][2:end]), vectors[:, i][1]), 1:num_vectors)

        # Anderson-Darling Test
        AD_p_value = pvalue(OneSampleADTest(thetas, theta_dist))

        # Kolmogorov-Smirnov Test
        KS_p_value = pvalue(ExactOneSampleKSTest(thetas, theta_dist))

        # Update attempt count
        num_attempts += 1
    end

    # --- Save uniformity test results

    uniformity_stats = (
        success=(num_attempts < max_attempts),
        num_attempts=num_attempts,
        AD_p_value=AD_p_value,
        KS_p_value=KS_p_value,
    )

    return vectors, uniformity_stats
end

# --- Generate sample of vectors

# Dimension of space
n_2d = 2;

# Generate sample of vectors
vectors, uniformity_stats = generate_vectors(n_2d, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit circle.")
    display(uniformity_stats)
end

# --- Compute inner products

inner_products = Vector{Float64}()
for i = 1:sample_size
    vector_i = vectors[:, i]

    for j = 1:sample_size
        # Skip inner product of vector with itself
        if i == j
            continue
        end

        push!(inner_products, dot(vector_i, vectors[:, j]))
    end
end

# --- Compute distributions

# Analytical PDF
pdf_2d(x) = 1 / π / sqrt(1 - x^2)

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins = range(-1, 1; length=num_hist_bins)
plt = plot(histogram(inner_products; bins=hist_bins, normalize=true))
plot!(pdf_2d, range(-1, 1; length=100))

hist_bins_abs = range(0, 1; length=num_hist_bins)
abs_inner_products = abs.(inner_products)
plt_abs = plot(histogram(abs_inner_products; bins=hist_bins_abs, normalize=true))

display(plt)
display(plt_abs)

# Compute empirical volume distributions
hist = fit(Histogram, inner_products; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

hist_abs = fit(Histogram, abs_inner_products; nbins=num_hist_bins)
edf_abs = normalize(hist; mode=:pdf);

# --- Generate sample of vectors
# Dimension of space
n_3d = 3;

# Generate sample of vectors
vectors, uniformity_stats = generate_vectors(n_3d, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit sphere.")
    display(uniformity_stats)
end

# --- Compute inner products

inner_products = Vector{Float64}()
for i = 1:sample_size
    vector_i = vectors[:, i]

    for j = 1:sample_size
        # Skip inner product of vector with itself
        if i == j
            continue
        end

        push!(inner_products, dot(vector_i, vectors[:, j]))
    end
end

# --- Compute distributions

# Analytical PDF
pdf_3d(x) = 0.5

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins = range(-1, 1; length=num_hist_bins)
plt = plot(histogram(inner_products; bins=hist_bins, normalize=true))
plot!(pdf_3d, range(-1, 1; length=100))

hist_bins_abs = range(0, 1; length=num_hist_bins)
abs_inner_products = abs.(inner_products)
plt_abs = plot(histogram(abs_inner_products; bins=hist_bins_abs, normalize=true))

display(plt)
display(plt_abs)

# Compute empirical volume distributions
hist = fit(Histogram, inner_products; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

hist_abs = fit(Histogram, abs_inner_products; nbins=num_hist_bins)
edf_abs = normalize(hist; mode=:pdf);

# --- Generate sample of vectors

# Dimension of space
n = 100;

# Generate sample of vectors
vectors, uniformity_stats = generate_vectors(n, sample_size)
if !uniformity_stats[:success]
    println("FAILED to generated sample with sufficient uniformity on unit $(n)-sphere.")
    display(uniformity_stats)
end

# --- Compute inner products

inner_products = Vector{Float64}()
for i = 1:sample_size
    vector_i = vectors[:, i]

    for j = 1:sample_size
        # Skip inner product of vector with itself
        if i == j
            continue
        end

        push!(inner_products, dot(vector_i, vectors[:, j]))
    end
end

# --- Compute distributions

# ------ Analytical PDF: proportional to (sin(θ))^(n-3)

# Compute normalization factor
k = n - 3
if k % 2 == 0
    norm_ = 1 / π
else
    norm_ = 0
    for j = 0:(k-1)÷2
        norm_ += (-1)^j * binomial(big(k), j) / (k - 2*j)
    end
    norm_ *= 4
    norm_ = abs(norm_)
end
pdf_nd(x) = (2^big(k) / norm_) * (1 - x^2)^(k/2)

# Histogram parameters
num_hist_bins = 50

# Plot histograms
hist_bins = range(-1, 1; length=num_hist_bins)
plt = plot(histogram(inner_products; bins=hist_bins, normalize=true))
plot!(pdf_nd, range(-1, 1; length=100))

hist_bins_abs = range(0, 1; length=num_hist_bins)
abs_inner_products = abs.(inner_products)
plt_abs = plot(histogram(abs_inner_products; bins=hist_bins_abs, normalize=true))

display(plt)
display(plt_abs)

# Compute empirical volume distributions
hist = fit(Histogram, inner_products; nbins=num_hist_bins)
edf = normalize(hist; mode=:pdf);

hist_abs = fit(Histogram, abs_inner_products; nbins=num_hist_bins)
edf_abs = normalize(hist; mode=:pdf);

max_stars_repo_name: ktchu/RESEARCH-math-exploration

max_issues_repo_name: ktchu/RESEARCH-math-exploration

max_forks_repo_name: ktchu/RESEARCH-math-exploration


text:
using DiffEqSensitivity, OrdinaryDiffEq, ParameterizedFunctions, Test

A = [1,0,2,3]
function f_sobol(p)
    x = A[1]*p[1] + A[2]*p[2]
    y = A[3]*p[1] + A[4]*p[2]
    [x,y]
end
p_range = [[0.5,1.5],[0.5,1.5]]
N = 1000000
total = sobol_sensitivity(f_sobol,p_range,N,0)
first_order = sobol_sensitivity(f_sobol,p_range,N,1)
second_order = sobol_sensitivity(f_sobol,p_range,N,2)

@test [total[1][1],total[2][1]] ≈ [first_order[1][1],first_order[2][1]] .+ second_order[1][1] atol=1e-1
@test [total[1][2],total[2][2]] ≈ [first_order[1][2],first_order[2][2]] .+ second_order[1][2] atol=1e-1


text:
#=

functions to estimate area and lineage trait averages
and lineage specific differences

Ignacio Quintero Mächler

t(-_-t)

May 15 2017

=#

"""
    deltaX!(δX   ::Array{Float64,3},
            X    ::Array{Float64,2},
            m    ::Int64,
            ntip ::Int64,
            narea::Int64)

Estimate pairwise trait and range distances between all lineages
conditional on a given `δY`.
"""
function deltaX!(δX   ::Array{Float64,3},
                 X    ::Array{Float64,2},
                 wcol ::Array{Array{Int64,1},1},
                 m    ::Int64,
                 ntip ::Int64,
                 narea::Int64)

  @inbounds begin

    for i = Base.OneTo(m), l = wcol[i], j = wcol[i]
      l == j && continue
      δX[j,l,i] = X[i,j] - X[i,l]
    end

  end

  return nothing
end

"""
    deltaY!(δY   ::Array{Float64,3},
            Y    ::Array{Int64,3},
            m    ::Int64,
            ntip ::Int64,
            narea::Int64)

Estimate pairwise trait and range distances between all lineages
conditional on a given `δY`.
"""
function deltaY!(δY   ::Array{Float64,3},
                 Y    ::Array{Int64,3},
                 wcol ::Array{Array{Int64,1},1},
                 m    ::Int64,
                 ntip ::Int64,
                 narea::Int64)

  @inbounds begin

    for i = Base.OneTo(m), l = wcol[i], j = wcol[i]

      j == l && continue

      sl        = 0.0
      δY[l,j,i] = 0.0
      @simd for k = Base.OneTo(narea)
        if Y[i,j,k] == 1
          sl        += 1.0
          δY[l,j,i] += Float64(Y[i,l,k])
        end
      end
      δY[l,j,i] /= sl
    end

  end

  return nothing
end

"""
    deltaXY!(δX   ::Array{Float64,3},
             δY   ::Array{Float64,3},
             X    ::Array{Float64,2},
             Y    ::Array{Int64,3},
             m    ::Int64,
             ntip ::Int64,
             narea::Int64)

Estimate pairwise trait and range distances between all lineages.
"""
function deltaXY!(δX   ::Array{Float64,3},
                  δY   ::Array{Float64,3},
                  X    ::Array{Float64,2},
                  Y    ::Array{Int64,3},
                  wcol ::Array{Array{Int64,1},1},
                  m    ::Int64,
                  ntip ::Int64,
                  narea::Int64)

  @inbounds begin

    for i = Base.OneTo(m), l = wcol[i], j = wcol[i]

      j == l && continue

      # X differences
      δX[j,l,i] = X[i,j] - X[i,l]

      # Y area overlap
      sl        = 0.0
      δY[l,j,i] = 0.0
      @simd for k = Base.OneTo(narea)
        if Y[i,j,k] == 1
          sl        += 1.0
          δY[l,j,i] += Float64(Y[i,l,k])
        end
      end
      δY[l,j,i] /= sl
    end

  end

  return nothing
end

"""
    sde!(LA   ::Array{Float64,2},
         δX   ::Array{Float64,3},
         δY   ::Array{Float64,3},
         m    ::Int64,
         ntip ::Int64)

Estimate the lineage averages for the SDE of trait evolution
"""
function sde!(LAp   ::Array{Float64,2},
              LAn   ::Array{Float64,2},
              δX   ::Array{Float64,3},
              δY   ::Array{Float64,3},
              wcol ::Array{Array{Int64,1},1},
              m    ::Int64,
              ntip ::Int64)

  @inbounds begin

    for i = Base.OneTo(m), j = wcol[i]
      LAp[i,j] = 0.0
      LAn[i,j] = 0.0
      for l = wcol[i]
        l == j && continue
        y = δY[l,j,i]
        iszero(y) && continue
        x = δX[l,j,i]
        LAp[i,j] += x*y
        LAn[i,j] += sign(x) * y * exp(-abs(x))
      end
    end
  end

  return nothing
end

"""
    lindiff!(LD   ::Array{Float64,3},
             δX   ::Array{Float64,3},
             Y    ::Array{Int64,3},
             m    ::Int64,
             ntip ::Int64,
             narea::Int64)

Estimate area-specific distance averages
for colonization and extinction rates.
"""
function lindiff!(LD   ::Array{Float64,3},
                  δX   ::Array{Float64,3},
                  Y    ::Array{Int64,3},
                  wcol ::Array{Array{Int64,1},1},
                  m    ::Int64,
                  ntip ::Int64,
                  narea::Int64)
  @inbounds begin

    for k = Base.OneTo(narea), i = Base.OneTo(m), l = wcol[i]

      xmin = Inf
      for j = wcol[i]
        j == l && continue
        if Y[i,j,k] == 1
          x = abs(δX[j,l,i])
          iszero(x) && continue
          xmin = x < xmin ? x : xmin
        end
      end
      LD[i,l,k] = isinf(xmin) ? 0.0 : xmin
    end
  end

  return nothing
end

"""
    Xupd_linavg!(δxi  ::Array{Float64,2},
                 lapi ::Array{Float64,1},
                 lani ::Array{Float64,1},
                 ldi  ::Array{Float64,2},
                 wci  ::Array{Int64,1},
                 xpi  ::Array{Float64,1},
                 xi   ::Int64,
                 xj   ::Int64,
                 Y    ::Array{Int64,3},
                 δyi  ::Array{Float64,2},
                 narea::Int64)

Re-estimate lineage specific means
for a node update when `ωx >= 0.0`
"""
function Xupd_linavg!(δxi  ::Array{Float64,2},
                      lapi ::Array{Float64,1},
                      lani ::Array{Float64,1},
                      ldi  ::Array{Float64,2},
                      wcol ::Array{Array{Int64,1},1},
                      xpi  ::Array{Float64,1},
                      xi   ::Int64,
                      xj   ::Int64,
                      Y    ::Array{Int64,3},
                      δY   ::Array{Float64,3},
                      narea::Int64)

  @inbounds begin

    # estimate pairwise distances
    wci = wcol[xi]
    for l = wci, j = wci
      l == j && continue
      δxi[j,l] = xpi[j] - xpi[l]
    end

    # estimate lineage averages
    for l = wci
      lapi[l] = 0.0
      lani[l] = 0.0
      for j = wci
        j == l && continue
        y = δY[j,l,xi]
        iszero(y) && continue
        x = δxi[j,l]
        lapi[l] += x * y
        lani[l] += sign(x) * y * exp(-abs(x))
      end
    end

    # estimate lineage sum of distances
    for k = Base.OneTo(narea), l = wci
      xmin = Inf
      for j = wci
        j == l && continue
        if Y[xi,j,k] == 1
          x = abs(δxi[j,l])
          iszero(x) && continue
          xmin = x < xmin ? x : xmin
        end
      end
      ldi[l,k] = xmin == Inf ? 0.0 : xmin
    end
  end

  return nothing
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# Interval
# ========
#
# Base interval types and utilities.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Note, just to be clear: this shadows IntervalTrees.Interval
"A genomic interval specifies interval with some associated metadata"
type Interval{T} <: AbstractInterval{Int64}
    seqname::StringField
    first::Int64
    last::Int64
    strand::Strand
    metadata::T

    function Interval(seqname, first, last, strand, metadata)
        return new(seqname, first, last, strand, metadata)
    end

    function Interval()
        return new(StringField(), 0, 0, STRAND_NA, T())
    end
end

function Interval{T}(seqname::AbstractString, first::Integer, last::Integer,
                    strand::Union{Strand,Char}, metadata::T)
    return Interval{T}(convert(StringField, seqname), first, last, strand, metadata)
end

function Interval(seqname::AbstractString, first::Integer, last::Integer,
                  strand::Union{Strand,Char}=STRAND_BOTH)
    return Interval{Void}(seqname, first, last, strand, nothing)
end

function Base.copy{T}(interval::Interval{T})
    return Interval{T}(copy(interval.seqname), interval.first, interval.last,
                       interval.strand, copy(interval.metadata))
end

function Bio.Seq.seqname(i::Interval)
    return i.seqname
end

function Bio.Seq.metadata(i::Interval)
    return i.metadata
end

function strand(i::Interval)
    return i.strand
end

"""
    leftposition(i::Interval)

Return the leftmost position of `i`.
"""
function leftposition(i::Interval)
    return i.first
end

"""
    rightposition(i::Interval)

Return the rightmost position of `i`.
"""
function rightposition(i::Interval)
    return i.last
end

IntervalTrees.first(i::Interval) = i.first
IntervalTrees.last(i::Interval) = i.last

function Base.isless{T}(a::Interval{T}, b::Interval{T},
                        seqname_isless::Function=alphanum_isless)
    if a.seqname != b.seqname
        return seqname_isless(a.seqname, b.seqname)::Bool
    elseif a.first != b.first
        return a.first < b.first
    elseif a.last != b.last
        return a.last < b.last
    elseif a.strand != b.strand
        return a.strand < b.strand
    else
        return false
    end
end

"""
Check if two intervals are well ordered.

Intervals are considered well ordered if a.seqname <= b.seqnamend and
a.first <= b.first.
"""
function isordered{T}(a::Interval{T}, b::Interval{T},
                      seqname_isless::Function=alphanum_isless)
    if a.seqname != b.seqname
        return seqname_isless(a.seqname, b.seqname)::Bool
    elseif a.first != b.first
        return a.first < b.first
    else
        return true
    end
end

"""
Return true if interval `a` entirely precedes `b`.
"""
function precedes{T}(a::Interval{T}, b::Interval{T},
                     seqname_isless::Function=alphanum_isless)
    return (a.last < b.first && a.seqname == b.seqname) ||
        seqname_isless(a.seqname, b.seqname)::Bool
end

function Base.:(==){T}(a::Interval{T}, b::Interval{T})
    return a.seqname  == b.seqname &&
           a.first    == b.first &&
           a.last     == b.last &&
           a.strand   == b.strand &&
           a.metadata == b.metadata
end

"Return true if interval `a` overlaps interval `b`, with no consideration to strand"
function isoverlapping{S, T}(a::Interval{S}, b::Interval{T})
    return a.first <= b.last && b.first <= a.last && a.seqname == b.seqname
end

function Base.show(io::IO, i::Interval)
    if get(io, :compact, false)
        print(io, i.seqname, ":", i.first, "-", i.last, "  ", i.strand, "  ", i.metadata)
    else
        println(io, summary(i), ':')
        println(io, "  sequence name: ", i.seqname)
        println(io, "  leftmost position: ", i.first)
        println(io, "  rightmost position: ", i.last)
        println(io, "  strand: ", i.strand)
          print(io, "  metadata: ", i.metadata)
    end
end

"""
A comparison function used to sort on numbers within text.

This is useful since sequences are often named things like "chr12" or
"read1234". Treating the numbers as numbers and not text gives a more natural
ordering.

This is similar to the '--version-sort' option in GNU coreutils sort.
"""
function alphanum_isless(a::AbstractString, b::AbstractString)
    i = 1
    j = 1

    # match up to the first digit
    k0 = 0 # position of first digit
    while i <= length(a) && j <= length(b)
        if isdigit(a[i]) && isdigit(b[j])
            k0 = i
            break
        else
            if a[i] != b[j]
                return a[i] < b[j]
            end
        end
        i = nextind(a, i)
        j = nextind(b, j)
    end

    # match numbers
    ka1, kb1 = 0, 0
    while i <= length(a) && isdigit(a[i])
        ka1 = i
        i = nextind(a, i)
    end
    while j <= length(b) && isdigit(b[j])
        kb1 = j
        j = nextind(b, j)
    end

    if ka1 == 0 && kb1 != 0
        return true
    elseif ka1 != 0 && kb1 == 0
        return false
    elseif ka1 != 0 && kb1 != 0
        aval = parse(Int, a[k0:ka1])
        bval = parse(Int, b[k0:kb1])
        if aval != bval
            return aval < bval
        end
    end

    # match suffixes
    while i <= length(a) && j <= length(b)
        if a[i] != b[j]
            return a[i] < b[j]
        end
        i = nextind(a, i)
        j = nextind(b, j)
    end

    return j <= length(b)
end

"""
A type deriving `IntervalStream{T}` must be iterable and produce
Interval{T} objects in sorted order.
"""
abstract IntervalStream{T}

typealias IntervalStreamOrArray{T} Union{Vector{Interval{T}},IntervalStream{T},Bio.IO.AbstractReader}

function metadatatype{T}(::IntervalStream{T})
    return T
end


text:
#list of added packages

#routines for adding the packages
add Atom Juno
add DataFrames Revise Distributions StatsBase GLM Dates NLopt ForwardDiff Formatting
add DataStructures LinearAlgebra Gadfly Measures Formatting CodecZlib CodecLz4
add CSV GZip NLopt CategoricalArrays Serialization Compose Distributed Random Cairo
add Fontconfig Compose StaticArrays VegaLite CuArrays
add https://github.com/clintonTE/Finometrics

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# This file contains utility functions used throughout the module
# which are not to be exported to the user

function check_pa_error(err::Int32)
    if err < 0
        errtext = ccall((:Pa_GetErrorText, portaudio), Ptr{Uint8}, (Int32,), err)
        error(bytestring(errtext))
    end
end


text:
import ImageCore: clamp01, gamutmax
import Statistics: mean

"""
    ColorJitter <: ColorOperation

Description
--------------

Adjusts the brightness and contrast of an image according to the formula
`α * image[i] + β * M`, where `M` is either `mean(image)` or the maximum
intensity value.

Usage
--------------

    ColorJitter()
    ColorJitter(α, β; [usemax])

Arguments
--------------

- **`α`** : `Real` or `AbstractVector` of `Real` that denote the coefficient(s)
    for contrast adjustment. Defaults to `0.8:0.1:1.2`.
- **`β`** : `Real` or `AbstractVector` of `Real` that denote the coefficient(s)
    for brightness adjustment. Defaults to `-0.2:0.1:0.2`.
- **`usemax::Bool`**: Optional. If `true`, the brightness will be adjusted by
    the maximum intensity value; otherwise, the image mean will be used.
    Defaults to `true`.

Examples
--------------

```
using Augmentor
img = testpattern()

# use exactly 1.2 for contrast, and one of 0.5 and 0.8 for brightness
augment(img, ColorJitter(1.2, [0.5, 0.8]))

# pick the coefficients randomly from the specified ranges
augment(img, ColorJitter(0.8:0.1:2.0, 0.5:0.1:1.1))
```
"""
struct ColorJitter{A<:AbstractVector, B<:AbstractVector} <: ColorOperation
    α::A
    β::B
    usemax::Bool

    function ColorJitter(α::A, β::B, usemax) where {A<:AbstractVector{<:Real},
                                                    B<:AbstractVector{<:Real}}
        length(α) > 0 || throw(ArgumentError("Range $(α) is empty"))
        length(β) > 0 || throw(ArgumentError("Range $(β) is empty"))
        new{A, B}(α, β, usemax)
    end
end

# The default values for α and β are taken from Albumentations
ColorJitter() = ColorJitter(0.8:0.1:1.2, -0.2:0.1:0.2, true)
ColorJitter(α, β; usemax=true) = ColorJitter(vectorize(α), vectorize(β), usemax)

randparam(op::ColorJitter, img) = (safe_rand(op.α), safe_rand(op.β))

@inline supports_eager(::Type{<:ColorJitter}) = true
@inline supports_lazy(::Type{<:ColorJitter}) = true

function applyeager(op::ColorJitter, img::AbstractArray, (α, β))
    M = _get_M(op, img)
    return _map_pix.(α, β, M, img)
end

function applylazy(op::ColorJitter, img::AbstractArray, (α, β))
    M = _get_M(op, img)
    return mappedarray(pix -> _map_pix(α, β, M, pix), img)
end

function showconstruction(io::IO, op::ColorJitter)
    print(io, typeof(op).name.name, '(', op.α, ", ", op.β, ", ", op.usemax, ')')
end

function Base.show(io::IO, op::ColorJitter)
    if get(io, :compact, false)
        maxmsg = op.usemax ? "max. intensity" : "mean value";
        print(io, "Color jitter with coffecients α=$(op.α) and β=$(op.β) (w.r.t. $(maxmsg))")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end

_map_pix(α, β, M, pix) = clamp01(α * pix + β * M)

function _get_M(op::ColorJitter, img)
    T = eltype(img)
    if op.usemax
        return T(gamutmax(T)...)
    else
        return convert(T, mean(img))
    end
end


text:
# TODO: Add grid based visualization

"""
    write_vtk(name, surfaces, [surface_properties]; kwargs...)
    write_vtk(name, wakes; kwargs...)
    write_vtk(name, surfaces, wakes, [surface_properties]; kwargs...)

Write geometry from surfaces and/or wakes to Paraview files for visualization.

# Arguments
 - `name`: Base name for the generated files
 - `surfaces`:
   - Vector of grids of shape (3, nc+1, ns+1) which represent lifting surfaces
   or
   - Vector of matrices of shape (nc, ns) containing surface panels (see
    [`SurfacePanel`](@ref))
   where `nc` is the number of chordwise panels and `ns` is the number of
   spanwise panels
 - `wakes`: (optional) Vector of wakes corresponding to each surface, represented
    by matrices of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns) where
    `nw` is the number of chordwise wake panels and `ns` is the number of
    spanwise panels.
 - `surface_properties`: (optional) Vector of surface panel properties for each
    surface, stored as matrices of panel properties (see [`PanelProperties`](@ref))
    of shape (nc, ns) where `nc` is the number of chordwise panels and `ns` is
    the number of spanwise panels

# Keyword Arguments:
 - `symmetric`: (required if `surface_properties` is provided) Flags indicating whether a
    mirror image (across the X-Z plane) was used when calculating induced velocities
    for each surface.
 - `trailing_vortices`: Flag indicating whether the model uses trailing vortices.
    Defaults to `true` when wake panels are absent, `false` otherwise
 - `xhat`: Direction in which trailing vortices extend if used. Defaults to [1, 0, 0].
 - `wake_length`: Distance to extend trailing vortices. Defaults to 10
 - `metadata`: Dictionary of metadata to include in generated files
"""
write_vtk(name, surfaces::AbstractVector{<:AbstractMatrix}, args...; kwargs...)

function write_vtk(name, surfaces::AbstractVector{<:AbstractMatrix{<:SurfacePanel}},
    properties=nothing; symmetric=fill(nothing, length(surfaces)), kwargs...) where T

    # create paraview multiblock file
    vtk_multiblock(name) do vtmfile
        # loop through all surfaces
        for i = 1:length(surfaces)
            # add paraview files corresponding to the surface to the multiblock file
            write_vtk!(vtmfile, surfaces[i], properties[i]; symmetric=symmetric[i], kwargs...)
        end
    end

    return nothing
end

function write_vtk(name, wakes::AbstractVector{<:AbstractMatrix{<:WakePanel}}; kwargs...) where T

    # create paraview multiblock file
    vtk_multiblock(name) do vtmfile
        # loop through all wakes
        for i = 1:length(wakes)
            # add paraview files corresponding to the wake to the multiblock file
            write_vtk!(vtmfile, wakes[i]; kwargs..., surface_index = i)
        end
    end

    return nothing
end

function write_vtk(name, surfaces::AbstractVector{<:AbstractMatrix{<:SurfacePanel}},
    wakes::AbstractVector{<:AbstractMatrix{<:WakePanel}}, properties=nothing; kwargs...)

    # create multiblock file
    vtk_multiblock(name) do vtmfile

        # loop through all surfaces
        for i = 1:length(surfaces)

            # extract circulation at the trailing edge of `surface`
            if isnothing(properties)
                surface_circulation = zeros(size(surfaces[i], 2))
            else
                surface_circulation = getproperty.(properties[i][end,:], :gamma)
            end

            # extract circulation at the leading edge of `wake`
            wake_circulation = getproperty.(wakes[i][1,:], :gamma)

            # add paraview files corresponding to the surface to the multiblock file
            write_vtk!(vtmfile, surfaces[i], properties[i]; wake_circulation,
                trailing_edge = isempty(wakes[i]), kwargs..., trailing_vortices = false)

            # add paraview files corresponding to the wake to the multiblock file
            write_vtk!(vtmfile, wakes[i]; surface_circulation, kwargs...)
        end
    end

    return nothing
end

"""
    write_vtk(name, surface_history, property_history, wake_history; kwargs...)

Writes unsteady simulation geometry to Paraview files for visualization.

# Arguments
 - `name`: Base name for the generated files
 - `surface_history`: Vector of surfaces at each time step, where each surface is
    represented by a matrix of surface panels (see [`SurfacePanel`](@ref)) of shape
    (nc, ns) where `nc` is the number of chordwise panels and `ns` is the number
    of spanwise panels
 - `property_history`: Vector of surface properties for each surface at each
    time step, where surface properties are represented by a matrix of panel
    properties (see [`PanelProperties`](@ref)) of shape (nc, ns) where `nc` is
    the number of chordwise panels and `ns` is the number of spanwise panels
 - `wake_history`: Vector of wakes corresponding to each surface at each time step,
    where each wake is represented by a matrix of wake panels (see [`WakePanel`](@ref))
    of shape (nw, ns) where `nw` is the number of chordwise wake panels and
    `ns` is the number of spanwise panels.
 - `dt`: Time step vector

# Keyword Arguments:
 - `symmetric`: (required if `properties` is provided) Flags indicating whether a
    mirror image (across the X-Z plane) was used when calculating induced velocities
    for each surface.
 - `wake_length`: Distance to extend trailing vortices. Defaults to 10
 - `metadata`: Dictionary of metadata to include in generated files
"""
function write_vtk(name, surface_history::AbstractVector{<:AbstractVector{<:AbstractMatrix}},
    property_history::AbstractVector{<:AbstractVector{<:AbstractMatrix}},
    wake_history::AbstractVector{<:AbstractVector{<:AbstractMatrix}}, dt;
    symmetric = fill(nothing, length(surface_history[1])), kwargs...)

    symmetric = isa(symmetric, Number) ? fill(symmetric, length(surface_history[1])) : symmetric

    # create paraview collection file
    paraview_collection(name) do pvdfile

        # construct time vector
        time = cumsum(dt)

        # loop through each time step
        for it = 1:length(time)

            # construct multiblock file for each time step
            vtk_multiblock(name*"-step$it") do vtmfile

                # loop through all surfaces
                for i = 1:length(surface_history[it])

                    # extract circulation at the trailing edge of `surface`
                    surface_circulation = getproperty.(property_history[it][i][end,:], :gamma)

                    # extract circulation at the leading edge of `wake`
                    if isempty(wake_history[it][i])
                        wake_circulation = zeros(size(surface_history[it][i], 2))
                    else
                        wake_circulation = getproperty.(wake_history[it][i][1,:], :gamma)
                    end

                    # add paraview files corresponding to the surface to the multiblock file
                    write_vtk!(vtmfile, surface_history[it][i], property_history[it][i];
                        wake_circulation,
                        trailing_edge = isempty(wake_history[it][i]),
                        symmetric = symmetric[i],
                        kwargs...,
                        trailing_vortices = false)

                    # add paraview files corresponding to the wake to the multiblock file
                    write_vtk!(vtmfile, wake_history[it][i];
                        surface_circulation,
                        symmetric = symmetric[i],
                        kwargs...)
                end

                # add multiblock file to the paraview collection file
                pvdfile[time[it]] = vtmfile
            end
        end
    end

    return nothing
end

"""
    write_vtk!(vtmfile, surface, [surface_properties]; kwargs...)

Writes geometry to Paraview files for visualization.

# Arguments
 - `vtmfile`: Multiblock file handle
 - `surface`: Matrix of surface panels (see [`SurfacePanel`](@ref)) of shape
    (nc, ns) where `nc` is the number of chordwise panels and `ns` is the number
    of spanwise panels
 - `surface_properties`: (optional) Matrix of panel properties for each non-wake panel
    where each element of the matrix is of type [`PanelProperties`](@ref).

# Keyword Arguments:
 - `symmetric`: (required if `properties` is provided) Flag indicating whether a
    mirror image (across the X-Z plane) was used when calculating induced velocities.
 - `trailing_vortices = true`: Flag indicating whether the model uses trailing vortices
 - `xhat = [1, 0, 0]`: Direction in which trailing vortices extend if used
 - `wake_length = 10`: Distance to extend trailing vortices
 - `wake_circulation = zeros(size(surfaces, 2))`: Contribution to the trailing
    edge circulation from the wake attached to this surface
 - `metadata = Dict()`: Dictionary of metadata to include in generated files
"""
function write_vtk!(vtmfile, surface::AbstractMatrix{<:SurfacePanel}, properties=nothing;
    symmetric = nothing,
    trailing_vortices = true,
    trailing_edge = true,
    xhat = SVector(1, 0, 0),
    wake_length = 10,
    wake_circulation = zeros(size(surface, 2)),
    metadata = Dict())

    # get float type
    TF = eltype(eltype(surface))

    # check to make sure `symmetric` is provided if `properties` is provided
    @assert !(!isnothing(properties) && isnothing(symmetric)) "Keyword argument `symmetric` is required when optional argument `properties` is provided"

    # grid dimensions
    nc, ns = size(surface)
    N = length(surface)

    # extract geometry as a grid
    xyz = Array{TF, 4}(undef, 3, nc+1, ns+1, 1)
    for (i, I) in enumerate(CartesianIndices((nc+1, ns+1)))
        if I[1] <= nc && I[2] <= ns
            xyz[:, I, 1] = top_left(surface[I[1], I[2]])
        elseif I[1] <= nc && I[2] == ns + 1
            xyz[:, I, 1] = top_right(surface[I[1], I[2]-1])
        elseif I[1] == nc + 1 && I[2] <= ns
            xyz[:, I, 1] = bottom_left(surface[I[1]-1, I[2]])
        else # I[1] == nc + 1 && I[2] == ns + 1
            xyz[:, I, 1] = bottom_right(surface[I[1]-1, I[2]-1])
        end
    end

    # convert grid to points
    points = reshape(xyz, 3, :)

    # now extract bound vortex geometries
    li = LinearIndices((size(xyz, 2), size(xyz, 3)))
    # horizontal bound vortices
    lines_h = [MeshCell(PolyData.Lines(), [
        li[i,j],
        li[i,j+1]
        ]) for j = 1:ns for i = 1:nc]
    # vertical bound vortices
    lines_v = [MeshCell(PolyData.Lines(), [
        li[i,j],
        li[i+1,j]
        ]) for j = 1:ns+1 for i = 1:nc]

    # now extract data (if applicable)
    if !isnothing(properties)
        # horizontal bound vortex circulation strength
        gamma_h = Matrix{TF}(undef, nc, ns)
        for i = 1:nc, j = 1:ns
            previous_gamma = i == 1 ? 0.0 : properties[i-1, j].gamma
            current_gamma = properties[i,j].gamma
            # check if we need to account for symmetry
            if symmetric && on_symmetry_plane(top_left(surface[i,j]), top_right(surface[i,j]))
                gamma_h[i,j] = 0.0
            else
                gamma_h[i,j] = current_gamma - previous_gamma
            end
        end

        # horizontal bound vortex force coefficient
        cf_h = Array{TF}(undef, 3, nc, ns)
        for i = 1:nc, j = 1:ns

            # check if we need to account for symmetry
            if symmetric && on_symmetry_plane(top_left(surface[i,j]), top_right(surface[i,j]))
                cf_h[:,i,j] .= 0.0
            else
                cf_h[:,i,j] = properties[i,j].cfb
            end
        end

        # horizontal bound vortex local velocity
        v_h = Array{TF}(undef, 3, nc, ns)
        for i = 1:nc, j = 1:ns
            v_h[:,i,j] = properties[i,j].velocity
        end

        # vertical bound vortex circulation strength
        gamma_v = Matrix{TF}(undef, nc, ns+1)
        for i = 1:nc
            current_gamma = properties[i,1].gamma
            # check if we need to account for symmetry
            if symmetric && on_symmetry_plane(bottom_left(surface[i,1]), top_left(surface[i,1]))
                gamma_v[i,1] = 0.0
            else
                gamma_v[i,1] = -current_gamma
            end
            for j = 2:ns
                previous_gamma = current_gamma
                current_gamma = properties[i,j].gamma
                # check if we need to account for symmetry
                if symmetric && on_symmetry_plane(bottom_left(surface[i,j]), top_left(surface[i,j]))
                    gamma_v[i,j] = 0.0
                else
                    gamma_v[i,j] = previous_gamma - current_gamma
                end
            end
            previous_gamma = current_gamma
            # check if we need to account for symmetry
            if symmetric && on_symmetry_plane(bottom_right(surface[i,end]), top_right(surface[i,end]))
                gamma_v[i,end] = 0.0
            else
                gamma_v[i,end] = previous_gamma
            end
        end

        # vertical bound vortex force coefficient
        cf_v = Array{TF}(undef, 3, nc, ns+1)
        for i = 1:nc
            #  check if we need to account for symmetry
            if symmetric && on_symmetry_plane(bottom_left(surface[i,1]), top_left(surface[i,1]))
                cf_v[:,i,1] .= 0.0
            else
                cf_v[:,i,1] = properties[i,1].cfl
            end
            for j = 2:ns
                previous_cf = properties[i,j-1].cfr
                current_cf = properties[i,j].cfl
                #  check if we need to account for symmetry
                if symmetric && on_symmetry_plane(bottom_left(surface[i,j]), top_left(surface[i,j]))
                    cf_v[i,j] = 0.0
                else
                    cf_v[:,i,j] = previous_cf + current_cf
                end
            end
            #  check if we need to account for symmetry
            if symmetric && on_symmetry_plane(bottom_right(surface[i,end]), top_right(surface[i,end]))
                cf_v[:,i,end] .= 0.0
            else
                cf_v[:,i,end] = properties[i,end].cfr
            end
        end
    end

    # horizontal bound vortices
    vtk_grid(vtmfile, points, lines_h) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        if !isnothing(properties)
            # circulation strength
            vtkfile["circulation"] = reshape(gamma_h, :)

            # local velocity
            vtkfile["velocity"] = reshape(v_h, 3, :)

            # force coefficient
            vtkfile["force"] = reshape(cf_h, 3, :)
        end
    end

    # vertical bound vortices
    vtk_grid(vtmfile, points, lines_v) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        if !isnothing(properties)
            # circulation strength
            vtkfile["circulation"] = reshape(gamma_v, :)

            # force coefficient
            vtkfile["force"] = reshape(cf_v, 3, :)
        end
    end

    # --- control points ---

    points_cp = Matrix{TF}(undef, 3, N)
    for i = 1:N
        ipoint = i
        points_cp[:,ipoint] = controlpoint(surface[i])
    end

    cells_cp = [MeshCell(PolyData.Verts(), [i]) for i = 1:N]

    vtk_grid(vtmfile, points_cp, cells_cp) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # add normal
        data = Matrix{TF}(undef, 3, N)
        for i = 1:length(surface)
            data[:,i] = normal(surface[i])
        end
        vtkfile["normal"] = data

    end

    # trailing edge and/or trailing vortices
    if trailing_vortices

        # trailing vortices points as a grid
        xyz_t = Array{TF}(undef, 3, 2, ns+1)
        for j = 1:ns
            xyz_t[:,1,j] = bottom_left(surface[end,j])
            xyz_t[:,2,j] = xyz_t[:,1,j] + wake_length*xhat
        end
        xyz_t[:,1,end] = bottom_right(surface[end,end])
        xyz_t[:,2,end] = xyz_t[:,1,end] + wake_length*xhat

        # generate points
        points_t = reshape(xyz_t, 3, :)

        li = LinearIndices((2, ns+1))

        # horizontal bound vortices
        lines_h = [MeshCell(PolyData.Lines(), [li[1,j], li[1,j+1]]) for j = 1:ns]

        # vertical bound vortices
        lines_v = [MeshCell(PolyData.Lines(), [li[1,j], li[2,j]]) for j = 1:ns+1]

        # combine
        lines_t = vcat(lines_h, lines_v)

        # now extract data (if applicable)
        if !isnothing(properties)
            # horizontal bound vortex circulation strength
            gamma_h = zeros(ns)

            # vertical bound vortex circulation strength
            gamma_v = Vector{TF}(undef, ns+1)
            current_gamma = properties[end,1].gamma

            # check if we need to account for symmetry
            if symmetric && on_symmetry_plane(bottom_left(surface[end,1]))
                gamma_v[1] = 0.0
            else
                gamma_v[1] = -current_gamma
            end

            for j = 2:ns
                previous_gamma = current_gamma
                current_gamma = properties[end,j].gamma

                # check if we need to account for symmetry
                if symmetric && on_symmetry_plane(bottom_left(surface[end,j]))
                    gamma_v[j] = 0.0
                else
                    gamma_v[j] = previous_gamma - current_gamma
                end

            end
            previous_gamma = current_gamma

            # check if we need to account for symmetry
            if symmetric && on_symmetry_plane(bottom_right(surface[end,end]))
                gamma_v[end] = 0.0
            else
                gamma_v[end] = previous_gamma
            end

            # combine
            gamma_t = vcat(gamma_h, gamma_v)
        end

    else
        # only generate if we plan to include the trailing edge
        if trailing_edge

            # generate points
            points_t = Matrix{TF}(undef, 3, ns+1)
            for j = 1:ns
                points_t[:,j] = bottom_left(surface[end,j])
            end
            points_t[:,end] = bottom_right(surface[end,end])

            # horizontal bound vortices
            lines_t = [MeshCell(PolyData.Lines(), j:j+1) for j = 1:ns]

            # now extract data (if applicable)
            if !isnothing(properties)
                # horizontal bound vortex circulation strength
                gamma_t = Vector{TF}(undef, ns)
                for j = 1:ns
                    # check if we need to account for symmetry
                    if symmetric && on_symmetry_plane(bottom_left(surface[end,j]), bottom_right(surface[end,j]))
                        gamma_t[j] = 0.0
                    else
                        gamma_t[j] = wake_circulation[j] - properties[end,j].gamma
                    end
                end
            end
        end

    end

    if trailing_vortices || trailing_edge

        # trailing edge and trailing vortices
        vtk_grid(vtmfile, points_t, lines_t) do vtkfile

            # add metadata
            for (key, value) in pairs(metadata)
                vtkfile[string(key)] = value
            end

            if !isnothing(properties)
                # circulation strength
                vtkfile["circulation"] = gamma_t
            end
        end
    end

    return nothing
end

"""
    write_vtk!(vtmfile, wake; kwargs...)

Writes geometry to Paraview files for visualization.

# Arguments
 - `vtmfile`: Paraview file handle
 - `wake`: Matrix of wake panels (see [`WakePanel`](@ref)) of shape (nw, ns)
    where `nw` is the number of chordwise wake panels and `ns` is the number of
    spanwise panels

# Keyword Arguments:
 - `symmetric`: (required) Flag indicating whether a mirror image (across the
    X-Z plane) was used when calculating induced velocities.
 - `trailing_vortices = false`: Flag indicating whether the model uses trailing vortices
 - `xhat = [1, 0, 0]`: Direction in which trailing vortices extend if used
 - `wake_length = 10`: Distance to extend trailing vortices
 - `surface_circulation = zeros(size(wake, 2))`: Contribution to the leading edge
    circulation from the surface attached to this wake.
 - `metadata = Dict()`: Dictionary of metadata to include in generated files
"""
function write_vtk!(vtmfile, wake::AbstractMatrix{<:WakePanel};
    symmetric,
    trailing_vortices = false,
    xhat = SVector(1, 0, 0),
    wake_length = 10,
    surface_circulation = zeros(size(wake, 2)),
    metadata=Dict())

    # do nothing if no wake panels are present
    if isempty(wake)
        return vtmfile
    end

    # extract float type
    TF = eltype(eltype(wake))

    # get wake dimensions
    nc, ns = size(wake)
    N = length(wake)

    # extract geometry as a grid
    xyz = Array{TF, 4}(undef, 3, nc+1, ns+1, 1)
    for (i, I) in enumerate(CartesianIndices((nc+1, ns+1)))
        if I[1] <= nc && I[2] <= ns
            xyz[:, I, 1] = top_left(wake[I[1], I[2]])
        elseif I[1] <= nc && I[2] == ns + 1
            xyz[:, I, 1] = top_right(wake[I[1], I[2]-1])
        elseif I[1] == nc + 1 && I[2] <= ns
            xyz[:, I, 1] = bottom_left(wake[I[1]-1, I[2]])
        else # I[1] == nc + 1 && I[2] == ns + 1
            xyz[:, I, 1] = bottom_right(wake[I[1]-1, I[2]-1])
        end
    end

    # convert to points
    points = reshape(xyz, 3, :)

    # now extract bound vortex geometries
    li = LinearIndices(size(xyz)[2:3])
    # horizontal bound vortices
    lines_h = [MeshCell(PolyData.Lines(), [
        li[i,j],
        li[i,j+1]
        ]) for j = 1:ns for i = 1:nc]
    # vertical bound vortices
    lines_v = [MeshCell(PolyData.Lines(), [
        li[i,j],
        li[i+1,j]
        ]) for j = 1:ns+1 for i = 1:nc]

    # horizontal bound vortex circulation strength
    gamma_h = Matrix{TF}(undef, nc, ns)
    for i in 1:nc, j = 1:ns
        previous_gamma = i == 1 ? surface_circulation[j] : wake[i-1, j].gamma
        current_gamma = wake[i,j].gamma
        if symmetric && on_symmetry_plane(top_left(wake[nc, ns]), top_right(wake[nc, ns]))
            gamma_h[i,j] = 0.0
        else
            gamma_h[i,j] = current_gamma - previous_gamma
        end
    end

    # vertical bound vortex circulation strength
    gamma_v = Matrix{TF}(undef, nc, ns+1)
    for i = 1:nc
        current_gamma = wake[i,1].gamma
        if symmetric && on_symmetry_plane(bottom_left(wake[i,1]), top_left(wake[i,1]))
            gamma_v[i,1] = 0.0
        else
            gamma_v[i,1] = -current_gamma
        end
        for j = 2:ns
            previous_gamma = current_gamma
            current_gamma = wake[i,j].gamma
            if symmetric && on_symmetry_plane(bottom_left(wake[i,j]), top_left(wake[i,j]))
                gamma_v[i,j] = 0.0
            else
                gamma_v[i,j] = previous_gamma - current_gamma
            end
        end
        previous_gamma = current_gamma
        if symmetric && on_symmetry_plane(bottom_right(wake[i,end]), top_right(wake[i,end]))
            gamma_v[i,end] = 0.0
        else
            gamma_v[i,end] = previous_gamma
        end
    end

    # horizontal bound vortices
    vtk_grid(vtmfile, points, lines_h) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # circulation strength
        vtkfile["circulation"] = reshape(gamma_h, :)
    end

    # vertical bound vortices
    vtk_grid(vtmfile, points, lines_v) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # circulation strength
        vtkfile["circulation"] = reshape(gamma_v, :)
    end

    if trailing_vortices

        # trailing vortices points as a grid
        xyz_t = Array{TF}(undef, 3, 2, ns+1)
        for j = 1:ns
            xyz_t[:,1,j] = bottom_left(wake[end,j])
            xyz_t[:,2,j] = xyz_t[:,1,j] + wake_length*xhat
        end
        xyz_t[:,1,end] = bottom_right(wake[end,end])
        xyz_t[:,2,end] = xyz_t[:,1,end] + wake_length*xhat

        # generate points
        points_t = reshape(xyz_t, 3, :)

        li = LinearIndices((2, ns+1))

        # horizontal bound vortices
        lines_h = [MeshCell(PolyData.Lines(), [li[1,j], li[1,j+1]]) for j = 1:ns]

        # vertical bound vortices
        lines_v = [MeshCell(PolyData.Lines(), [li[1,j], li[2,j]]) for j = 1:ns+1]

        # combine
        lines_t = vcat(lines_h, lines_v)

        # horizontal bound vortex circulation strength
        gamma_h = zeros(ns)

        # vertical bound vortex circulation strength
        gamma_v = Vector{TF}(undef, ns+1)
        current_gamma = wake[end,1].gamma
        if symmetric && on_symmetry_plane(bottom_left(wake[end,1]))
            gamma_v[1] = 0.0
        else
            gamma_v[1] = -current_gamma
        end
        for j = 2:ns
            previous_gamma = current_gamma
            current_gamma = wake[end,j].gamma
            if symmetric && on_symmetry_plane(bottom_left(wake[end,j]))
                gamma_v[j] = 0.0
            else
                gamma_v[j] = previous_gamma - current_gamma
            end
        end
        previous_gamma = current_gamma
        if symmetric && on_symmetry_plane(bottom_right(wake[end,end]))
            gamma_v[j] = 0.0
        else
            gamma_v[end] = previous_gamma
        end

        # combine
        gamma_t = vcat(gamma_h, gamma_v)

    else

        # generate points
        points_t = Matrix{TF}(undef, 3, ns+1)
        for j = 1:ns
            points_t[:,j] = bottom_left(wake[end,j])
        end
        points_t[:,end] = bottom_right(wake[end,end])

        # horizontal bound vortices
        lines_t = [MeshCell(PolyData.Lines(), j:j+1) for j = 1:ns]

        # horizontal bound vortex circulation strength
        gamma_t = Vector{TF}(undef, ns)
        for j = 1:ns
            if symmetric && on_symmetry_plane(bottom_left(wake[end,j]), bottom_right(wake[end,j]))
                gamma_t[j] = 0.0
            else
                gamma_t[j] = -wake[end,j].gamma
            end
        end

    end

    # trailing edge and trailing vortices
    vtk_grid(vtmfile, points_t, lines_t) do vtkfile

        # add metadata
        for (key, value) in pairs(metadata)
            vtkfile[string(key)] = value
        end

        # circulation strength
        vtkfile["circulation"] = reshape(gamma_t, :)
    end

    return nothing
end


text:
""" Tests of limit and colimit data structures.

More extensive tests are provided by (co)limits for specific categories such as
FinSet.
"""
module TestLimits
using Test

using Catlab.Theories, Catlab.CategoricalAlgebra

A, B, C = Ob(FreeCategory, :A, :B, :C)

# Limits
########

# Products.
f, g = Hom(:f, C, A), Hom(:g, C, B)
lim = Limit(ObjectPair(A,B), Span(f,g))
@test lim isa BinaryProduct
@test ob(lim) == C
@test apex(lim) == C
@test legs(lim) == [f,g]

lim = Limit(DiscreteDiagram([A,B]), Span(f,g))
@test lim isa Product

# Colimits
##########

# Coproducts.
f, g = Hom(:f, A, C), Hom(:g, B, C)
colim = Colimit(ObjectPair(A,B), Cospan(f,g))
@test colim isa BinaryCoproduct
@test ob(colim) == C
@test apex(colim) == C
@test legs(colim) == [f,g]

colim = Colimit(DiscreteDiagram([A,B]), Cospan(f,g))
@test colim isa Coproduct

end


text:
using NeuralVerification, LazySets, Test, LinearAlgebra
import NeuralVerification: ReLU, Id

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

net_path = joinpath(@__DIR__, "../networks/")

function test_runtime(problem, solver)
	solvername = string(typeof(solver))
	println("Solver: "*solvername)
	timed_result =@timed solve(solver, problem)
	println(" - Time: " * string(timed_result[2]) * " s")
	println(" - Output: ")
	println(timed_result[1])
end

# Problem: Small NN, input: HPolytope, output: HPolytope

small_nnet_file = net_path * "small_nnet.nnet"
small_nnet_id_file = net_path * "small_nnet_id.nnet"
small_nnet  = read_nnet(small_nnet_file, last_layer_activation = ReLU())
small_nnet_id  = read_nnet(small_nnet_id_file, last_layer_activation = Id())

in_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])
in_hpoly  = convert(HPolytope, in_hyper)
out_superset    = Hyperrectangle(low = [30.0], high = [80.0])

problem_small_HH    = Problem(small_nnet, in_hpoly, convert(HPolytope, out_superset))

# Problem: Small NN, input: Hyperrectangle, output: HPolytope

in_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])
problem_small_RH = Problem(small_nnet, in_hyper, HPolytope([HalfSpace([1.], 100.)]))

# Problem: Small NN, input: Hyperrectangle, output: Hyperrectangle

out_hyper    = Hyperrectangle(low = [30.0], high = [80.0])
problem_small_RR = Problem(small_nnet, in_hyper, out_hyper)

# Problem: ACAS, input: HPolytope, output: HPolytope

acas_file = net_path * "ACASXU_run2a_4_5_batch_2000.nnet"
acas_nnet = read_nnet(acas_file, last_layer_activation = Id())

A0 = Matrix{Float64}(I, 5, 5)
A1 = -Matrix{Float64}(I, 5, 5)
A = vcat(A0, A1)
b_lower = [ 0.21466922,  0.11140846, -0.4999999 ,  0.3920202 ,  0.4      ]
b_upper = [ 0.58819589,  0.4999999 , -0.49840835,  0.66474747,  0.4      ]
in_hyper  = Hyperrectangle(low = b_lower, high = b_upper)
inputSet = convert(HPolytope, in_hyper)

A = Matrix(undef, 2, 1)
A = [1.0, 0.0, 0.0, 0.0, -1.0]'
b = [0.0]
outputSet = HPolytope(A, b)

problem_acas_HH = Problem(acas_nnet, inputSet, outputSet)

# Problem: ACAS, input: Hyperrectangle, output: HPolytope

in_hyper  = Hyperrectangle(low = b_lower, high = b_upper)
problem_acas_HrH = Problem(acas_nnet, in_hyper, outputSet)

# Problem: ACAS, input: Hyperrectangle, output: HPolytope

in_hyper  = Hyperrectangle(low = b_lower, high = b_upper)
out_hyper  = Hyperrectangle(low = [-0.1], high = [3.0])
problem_acas_RR = Problem(acas_nnet, in_hyper, out_hyper)

# Problem: ACAS1 (one output), input: Hyperrectangle, output: Hyperrectangle
acas_file = net_path * "ACASXU_run2a_4_5_batch_2000_out1.nnet"
#ACASXU_run2a_4_5_batch_2000.nnet
acas_nnet1 = read_nnet(acas_file, last_layer_activation = Id())
in_hyper  = Hyperrectangle(low = b_lower, high = b_upper)
out_hyper  = Hyperrectangle(low = [-0.1], high = [3.0])

problem_acas1_RR = Problem(acas_nnet1, in_hyper, out_hyper)


text:
include("status-info.jl")
include("file-system.jl")
include("parsing.jl")
include("multiset.jl")
include("equality.jl")

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
include("decompose_cell.jl")
include("decompose_unstruct.jl")
include("decompose_struct.jl")

decompose(dataset::AbstractVTKMultiblockData, target::String="Faces") = decompose(VTKUnstructuredData(dataset), target)

function decompose(dataset::AbstractTimeSeriesVTKData, target::String="Faces")
    return VTKTimeSeriesData(dataset.timemarkers, [decompose(dataset[i], target) for
        i in 1:length(dataset)])
end


text:
"""
    group(by, iter)

Group the elements `x` of the iterable `iter` into groups labeled by `by(x)`.

# Example

```jldoctest
julia> group(iseven, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [1, 3, 5, 7, 9]
  true  => [2, 4, 6, 8, 10]
```
"""
group(by, iter) = group(by, identity, iter)

"""
    group(by, f, iter...)

Sorts the elements `x` of the iterable `iter` into groups labeled by `by(x)`,
transforming each element to `f(x)`. If multiple collections (of the same length)
are provided, the transformations occur elementwise.

# Example

```jldoctest
julia> group(iseven, x -> x ÷ 2, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [0, 1, 2, 3, 4]
  true  => [1, 2, 3, 4, 5]
```
"""
function group(by, f, iter)
    # TODO Do this inference-free, like comprehensions...
    T = eltype(iter)
    K = promote_op(by, T)
    V = promote_op(f, T)

    out = Dict{K, Vector{V}}()
    for x ∈ iter
        key = by(x)
        push!(get!(()->Vector{V}(), out, key), f(x))
    end
    return out
end

group(by, f, iter1, iter2, iters...) = group((x -> by(x...)), (x -> f(x...)), zip(iter1, iter2, iters...))

"""
    groupinds(by, iter...)

Sorts the indices `i` of `iter` into groups labeled by `by(iter[i])`. If multiple
collections (with matching indices) are provided, the groups are formed elementwise.

# Example

```jldoctest
julia> groupinds(iseven, [3,4,2,6,5,8])
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [1, 5]
  true  => [2, 3, 4, 6]
```
"""
function groupinds(by, iter)
    T = eltype(iter)
    K = promote_op(by, T)
    inds = keys(iter)
    V = eltype(inds)

    out = Dict{K, Vector{V}}()
    for i ∈ inds
        @inbounds x = iter[i]
        key = by(x)
        push!(get!(()->Vector{V}(), out, key), i)
    end
    return out
end

groupinds(by, iter1, iter2, iters...) = groupinds((x -> by(x...)), zip(iter1, iter2, iters...))

struct Groups{K, V, T, Inds} <: AbstractDict{K, V}
    parent::T
    inds::Inds
end

Base.keys(g::Groups) = keys(g.inds)
Base.@propagate_inbounds function Base.getindex(g::Groups, k)
    inds = g.inds[k]
    return view(g.parent, inds)
end
Base.length(g::Groups) = length(keys(g))
Base.start(g::Groups) = start(g.inds)
function Base.next(g::Groups, i)
    ((key, inds), i2) = next(g.inds, i)
    return (key => g.parent[inds], i2)
end
Base.done(g::Groups, i) = done(g.inds, i)

"""
    groupview(by, iter)

Like `group`, but each grouping is a view of the input collection `iter`.

# Example

```jldoctest
julia> v = [3,4,2,6,5,8]
6-element Array{Int64,1}:
 3
 4
 2
 6
 5
 8

julia> groups = groupview(iseven, v)
SAC.Groups{Bool,Any,Array{Int64,1},Dict{Bool,Array{Int64,1}}} with 2 entries:
  false => [3, 5]
  true  => [4, 2, 6, 8]

julia> groups[false][:] = 99
99

julia> v
6-element Array{Int64,1}:
 99
  4
  2
  6
 99
  8
```
"""
function groupview(by, iter)
    inds = groupinds(by, iter)
    V = promote_op(view, iter, eltype(inds))
    return Groups{keytype(inds), V, typeof(iter), typeof(inds)}(iter, inds)
end

groupview(by, iter1, iter2, iters...) = groupview((x -> by(x...)), zip(iter1, iter2, iters...))

"""
    groupreduce(by, op, iter)

Like `groupreduce(by, identity, op, iter)`.

# Example

```jldoctest
julia> groupreduce(iseven, +, 1:10)
Dict{Bool,Int64} with 2 entries:
false => 25
true  => 30
```
"""
groupreduce(by, op, iter) = groupreduce(by, identity, op, iter)

"""
    groupreduce(by, f, op, iter)

Like `groupreduce(by, f, op, v0, iter)`, except the reductions are started with the first
element in the group rather than `v0`.
"""
function groupreduce(by, f, op, iter)
    # TODO Do this inference-free, like comprehensions...
    T = eltype(iter)
    K = promote_op(by, T)
    T2 = promote_op(f, T)
    V = promote_op(op, T2, T2)

    out = Dict{K, V}()
    for x ∈ iter
        key = by(x)
        dict_index = Base.ht_keyindex2!(out, key)
        if dict_index > 0
            @inbounds out.vals[dict_index] = op(out.vals[dict_index], f(x))
        else
            Base._setindex!(out, convert(V, f(x)), key, -dict_index)
        end
    end
    return out
end

"""
    groupreduce(by, f, op, v0, iter...)

Applies a `mapreduce`-like operation on the groupings labeled by passing the elements of
`iter` through `by`. Mostly equivalent to `map(g -> reduce(op, v0, g), group(by, f, iter))`,
but designed to be more efficient. If multiple collections (of the same length) are
provided, the transformations `by` and `f` occur elementwise.
"""
function groupreduce(by, f, op, v0, iter)
    # TODO Do this inference-free, like comprehensions...
    T = eltype(iter)
    K = promote_op(by, T)
    T2 = promote_op(f, T)
    V = promote_op(op, T2, T2)

    out = Dict{K, V}()
    for x ∈ iter
        key = by(x)
        dict_index = Base.ht_keyindex2!(out, key)
        if dict_index > 0
            @inbounds out.vals[dict_index] = op(out.vals[dict_index], f(x))
        else
            Base._setindex!(out, convert(V, op(v0, f(x))), key, -dict_index)
        end
    end
    return out
end

groupreduce(by, f, op, v0, iter1, iter2, iters...) = groupreduce((x -> by(x...)), (x -> f(x...)), op, v0, zip(iter1, iter2, iters...))


text:
# deprecated formulation of the problem
module OldJuBasicMD

using LinearAlgebra

"""
    Boltzmann_likelihood(x, V, β)
Compute the unnormalized Boltzmann density, exp(-β V(x))
for potential V at inverse temperature β
"""
function Boltzmann_likelihood(x, V, β)
    w = exp(-β * V(x));
    return w
end

"""
    MALA_likelihood(X₀, X₁, gradV0, β, Δt)
Compute the likelihood of the MALA proposal X₀→X₁ for
potential with gradient gradV0 at inverse temperature β
and time step Δt
"""
function MALA_likelihood(X₀, X₁, gradV0, β, Δt)
    #norm of increment squared
    norminc2 = 0.0;
    for i = 1:length(X₀)
        norminc2 += (X₁[i]-X₀[i] + Δt * gradV0[i])^2;
    end
    return exp( -β * (norminc2) / (4*Δt) );
end

"""
    RWM(x₀, V, β, Δt, n_iters, return_trajectory=true)
Run an RWM sampler for distribution exp(-β V).  Proposal variance is
sqrt(2*Δt/β).  If `return_trajectory=true`, then the entire time series is
returned.  The acceptance rates are also returned.
"""
function RWM(x₀, V, β, Δt, n_iters; return_trajectory=true)

    # preallocate data structures
    X₀ = copy(x₀);
    Xp = similar(x₀);
    d = length(x₀)

    naccept = 0;

    if(return_trajectory)
        Xvals =zeros(d,n_iters);
        avals =zeros(n_iters);
    end

    V0 = V(X₀);
    gaussian_coef =  sqrt(2 * Δt /β);
    for j = 1:n_iters

        @. Xp = X₀ + gaussian_coef * randn();

        Vp = V(Xp);
        a = min(1, exp(β*(V0-Vp)))

        if rand()<a
            naccept = naccept+1;
            @. X₀ = Xp;
            V0 = Vp;
        end

        if(return_trajectory)
            @. Xvals[:,j] = X₀;
            avals[j] = naccept/j;
        end

    end
    if return_trajectory
        return Xvals, avals
    else
        return X₀, naccept/n_iters
    end

end # end RWM

"""
    RWM!(X₀, V, β, Δt, n_iters)
Run an in place RWM sampler for distribution exp(-β V).  Proposal variance is
sqrt(2*Δt/β).
"""
function RWM!(X₀, V, β, Δt, n_iters)

    # preallocate data structures
    Xp = similar(X₀);

    V0 = V(X₀);
    gaussian_coef =  sqrt(2 * Δt/β);
    for j = 1:n_iters

        @. Xp = X₀ + gaussian_coef * randn();

        Vp = V(Xp);
        a = min(1, exp(β*(V0-Vp)))

        if rand()<a
            @. X₀ = Xp;
            V0 = Vp;
        end
    end
    X₀
end # end RWM!

"""
    MALA(x₀, V, gradV!, β, Δt, n_iters, return_trajectory=true)
Run a MALA sampler for distribution exp(-β V).  Proposal variance is
sqrt(2*Δt/β).  If `return_trajectory=true`, then the entire time series is
returned.  The acceptance rates are also returned.
"""
function MALA(x₀, V, gradV!, β, Δt, n_iters; return_trajectory=true)

    # preallocate data structures
    X₀ = copy(x₀);
    Xp = similar(x₀);
    gradV0 = similar(x₀);
    gradVp = similar(x₀);
    d = length(x₀)

    naccept = 0;

    if(return_trajectory)
        Xvals =zeros(d,n_iters);
        avals =zeros(n_iters);
    end

    V0 = V(X₀);
    gaussian_coef =  sqrt(2 * Δt/β);
    gradV!(gradV0,X₀);
    for j = 1:n_iters

        @. Xp = X₀ - Δt * gradV0 + gaussian_coef * randn();

        Vp = V(Xp);

        gradV!(gradVp,Xp);

        g0 = MALA_likelihood(X₀, Xp, gradV0, β, Δt);
        gp = MALA_likelihood(Xp, X₀, gradVp, β, Δt);

        a = min(1, gp/g0 * exp(β*(V0-Vp)));

        if rand()<a
            naccept = naccept+1;
            @. X₀ = Xp;
            @. gradV0 = gradVp;
            V0 = Vp;
        end
        if(return_trajectory)
            @. Xvals[:,j] = X₀;
            avals[j] = naccept/j;
        end

    end
    if return_trajectory
        return Xvals, avals
    else
        return X₀, naccept/n_iters
    end

end # end MALA

"""
    MALA!(x₀, V, gradV!, β, Δt, n_iters, return_trajectory=true)
Run an in place MALA sampler for distribution exp(-β V).  Proposal variance is
sqrt(2*Δt/β).
"""
function MALA!(X₀, V, gradV!, β, Δt, n_iters)

    # preallocate data structures
    Xp = similar(X₀);
    gradV0 = similar(X₀);
    gradVp = similar(X₀);

    V0 = V(X₀);
    gaussian_coef =  sqrt(2 * Δt/β);
    gradV!(gradV0,X₀);
    for j = 1:n_iters

        @. Xp = X₀ - Δt * gradV0 + gaussian_coef * randn();

        Vp = V(Xp);

        gradV!(gradVp,Xp);

        g0 = MALA_likelihood(X₀, Xp, gradV0, β, Δt);
        gp = MALA_likelihood(Xp, X₀, gradVp, β, Δt);

        a = min(1, gp/g0 * exp(β*(V0-Vp)));

        if rand()<a
            @. X₀ = Xp;
            @. gradV0 = gradVp;
            V0 = Vp;
        end
    end
    X₀
end # end MALA!

"""
    EM(x₀, V, gradV!, β, Δt, n_iters, return_trajectory=true)
Run the Euler-Maruyama integrator for dX = -∇V(X)dt + sqrt(2/β)dW.   If
If `return_trajectory=true`, then the entire time series is returned.
"""
function EM(x₀, gradV!, β, Δt, n_iters; return_trajectory=true)

    # preallocate data structures
    X₀ = copy(x₀);
    gradV0 = similar(x₀);
    d = length(x₀)

    if(return_trajectory)
        Xvals =zeros(d, n_iters);
    end

    gaussian_coef = sqrt(2 * Δt / β);

    for j = 1:n_iters
        gradV!(gradV0,X₀);

        @. X₀ = X₀ - gradV0 * Δt +  gaussian_coef * randn();

        if(return_trajectory)
            @. Xvals[:,j] = X₀;
        end
    end
    if(return_trajectory)
        return Xvals
    else
        return X₀
    end

end

"""
    EM!(X₀, V, gradV!, β, Δt, n_iters, return_trajectory=true)
Run an in place Euler-Maruyama integrator for dX = -∇V(X)dt + sqrt(2/β)dW.
"""
function EM!(X₀, gradV!, β, Δt, n_iters)
    # preallocate data structures
    gradV0 = similar(X₀);

    gaussian_coef = sqrt(2 * Δt / β);

    for j = 1:n_iters
        gradV!(gradV0,X₀);

        @. X₀ = X₀ - gradV0 * Δt +  gaussian_coef * randn();
    end
    X₀
end # end EM!

"""
    BBK(q₀,p₀, gradV!, β, γ, M, Δt, n_iters, return_trajectory=true)
Run the BBK integrator for interial Langevin   If `return_trajectory=true`,
then the entire time series is returned.
"""
function BBK(q₀, p₀, gradV!, β, γ, M, Δt, n_iters; return_trajectory=true)

    # diffusion coefficient
    σ=sqrt(2*γ/β);
    d = length(q₀);

    c₀ = sqrt(0.5*Δt);

    Q₀ = copy(q₀);
    P₀ = copy(p₀);
    Phalf = similar(p₀);
    gradV0 = similar(q₀);

    if return_trajectory
        Qvals = zeros(d,n_iters);
        Pvals = zeros(d,n_iters);
    end

    gradV!(gradV0,Q₀);
    for j = 1:n_iters
        @. Phalf = P₀ - 0.5 * Δt * gradV0 - 0.5 * Δt * γ/M * P₀ + c₀ * σ * randn();
        @. Q₀ = Q₀ + Δt/M * Phalf;
        gradV!(gradV0, Q₀);
        @. P₀ = (Phalf -0.5 * Δt * gradV0 + c₀ * σ * randn())/(1 + 0.5 * Δt * γ/M);

        if return_trajectory
            @. Qvals[:,j] = Q₀;
            @. Pvals[:,j] = P₀;
        end

    end
    if return_trajectory
        return Qvals, Pvals
    else
        return Q₀, P₀
    end
end # end BBK

"""
    BBK!(q₀,p₀, gradV!, β, Δt, n_iters)
Run an in place BBK integrator for interial Langevin.
"""
function BBK!(Q₀, P₀,gradV!, β, γ, M, Δt, n_iters)

    # diffusion coefficient
    σ=sqrt(2*γ/β);
    d = length(Q₀);

    Phalf = similar(P₀);
    gradV0 = similar(Q₀);

    gradV!(gradV0,Q₀);
    for j = 1:n_iters
        @. Phalf = P₀ - 0.5 * Δt * gradV0 - 0.5 * Δt * γ/M * P₀ + sqrt(0.5*Δt) * σ * randn();
        @. Q₀ = Q₀ + Δt/M * Phalf;
        gradV!(gradV0, Q₀);
        @. P₀ = (Phalf -0.5 * Δt * gradV0 + sqrt(0.5 * Δt) * σ * randn())/(1 + 0.5 * Δt * γ/M);

    end
    Q₀, P₀

end # end BBK!

"""
    ABOBA(q₀,p₀, gradV!, β, γ, M, Δt, n_iters, return_trajectory=true)
Run the ABOBA integrator for interial Langevin.  If `return_trajectory=true`,
then the entire time series is returned.
"""
function ABOBA(q₀, p₀, gradV!, β, γ, M, Δt, n_iters; return_trajectory=true)

    c₀ = exp(-Δt * γ);
    c₁ = sqrt((1 - exp(-2*γ*Δt))/β);
    d = length(q₀);

    Q₀ = copy(q₀);
    P₀ = copy(p₀);
    Qhalf = similar(q₀);
    Phalf = similar(p₀);
    P̂half = similar(p₀);
    gradVhalf = similar(q₀);
    sqrtM = sqrt.(M);

    if return_trajectory
        Qvals = zeros(d,n_iters);
        Pvals = zeros(d,n_iters);
    end

    for j in 1:n_iters
        @. Qhalf = Q₀ + 0.5 * Δt * P₀/M;
        gradV!(gradVhalf, Qhalf);
        @. Phalf = P₀ - 0.5 * Δt * gradVhalf;
        @. P̂half = c₀ * Phalf + c₁ * sqrtM * randn();
        @. P₀ = P̂half - 0.5 * Δt * gradVhalf;
        @. Q₀ = Qhalf + 0.5 * Δt * P₀/M;

        if return_trajectory
            @. Qvals[:,j] = Q₀;
            @. Pvals[:,j] = P₀;
        end

    end
    if return_trajectory
        return Qvals, Pvals
    else
        return Q₀, P₀
    end
end

"""
    ABOBA!(Q₀,P₀, gradV!, β, γ, M, Δt, n_iters, return_trajectory=true)
Run an in place ABOBA integrator for interial Langevin.
"""
function ABOBA!(Q₀, P₀, gradV!, β, γ, M, Δt, n_iters)

    c₀ = exp(-Δt * γ);
    c₁ = sqrt((1 - exp(-2*γ*Δt))/β);
    d = length(Q₀);

    Qhalf = similar(Q₀);
    Phalf = similar(P₀);
    P̂half = similar(P₀);
    gradVhalf = similar(Q₀);
    sqrtM = sqrt.(M);

    for i in 1:n_iters
        @. Qhalf = Q₀ + 0.5 * Δt * P₀/M;
        gradV!(gradVhalf, Qhalf);
        @. Phalf = P₀ - 0.5 * Δt * gradVhalf;
        @. P̂half = c₀ * Phalf + c₁ * sqrtM * randn();
        @. P₀ = P̂half - 0.5 * Δt * gradVhalf;
        @. Q₀ = Qhalf + 0.5 * Δt * P₀/M;

    end
    Q₀, P₀

end

"""
    BAOAB(q₀,p₀, gradV!, β, γ, M, Δt, n_iters, return_trajectory=true)
Run the BAOAB integrator for interial Langevin.  If `return_trajectory=true`,
then the entire time series is returned.
"""
function BAOAB(q₀, p₀, gradV!, β, γ, M, Δt, n_iters; return_trajectory=true)

    c₀ = exp(-Δt * γ);
    c₁ = sqrt((1 - exp(-2*γ*Δt))/β);
    d = length(q₀);

    Q₀ = copy(q₀);
    P₀ = copy(p₀);
    Qhalf = similar(q₀);
    Phalf = similar(p₀);
    P̂half = similar(p₀);
    gradV0 = similar(q₀);
    sqrtM = sqrt.(M);

    if return_trajectory
        Qvals = zeros(d,n_iters);
        Pvals = zeros(d,n_iters);
    end

    gradV!(gradV0, Q₀);
    for j in 1:n_iters
        @. Phalf = P₀ - 0.5 * Δt * gradV0;
        @. Qhalf = Q₀ + 0.5 * Δt * Phalf/M;
        @. P̂half = c₀ * Phalf + c₁ * sqrtM * randn();
        @. Q₀ = Qhalf + 0.5 * Δt * P̂half/M;
        gradV!(gradV0, Q₀);
        @. P₀ = P̂half - 0.5 * Δt * gradV0;

        if return_trajectory
            @. Qvals[:,j] = Q₀;
            @. Pvals[:,j] = P₀;
        end

    end
    if return_trajectory
        return Qvals, Pvals
    else
        return Q₀, P₀
    end
end

"""
    BAOAB!(Q₀,P₀, gradV!, β, γ, M, Δt, n_iters, return_trajectory=true)
Run an in place BAOAB integrator for interial Langevin.
"""
function BAOAB!(Q₀, P₀, gradV!, β, γ, M, Δt, n_iters)

    c₀ = exp(-Δt * γ);
    c₁ = sqrt((1 - exp(-2*γ*Δt))/β);
    d = length(Q₀);

    Qhalf = similar(Q₀);
    Phalf = similar(P₀);
    P̂half = similar(P₀);
    gradV0 = similar(Q₀);
    sqrtM = sqrt.(M);

    gradV!(gradV0, Q₀);
    for j in 1:n_iters
        @. Phalf = P₀ - 0.5 * Δt * gradV0;
        @. Qhalf = Q₀ + 0.5 * Δt * Phalf/M;
        @. P̂half = c₀ * Phalf + c₁ * sqrtM * randn();
        @. Q₀ = Qhalf + 0.5 * Δt * P̂half/M;
        gradV!(gradV0, Q₀);
        @. P₀ = P̂half - 0.5 * Δt * gradV0;

    end
    Q₀, P₀

end

"""
    Verlet(q₀,p₀, gradV!, M, Δt, n_iters, return_trajectory=true)
Run the Verlet integrator.  If `return_trajectory=true`, then the entire time
series is returned.
"""
function Verlet(q₀, p₀, gradV!, M, Δt, n_iters; return_trajectory=true)

    d = length(q₀);

    Q₀ = copy(q₀);
    P₀ = copy(p₀);
    gradV0 = similar(Q₀);
    Phalf = similar(P₀);

    if return_trajectory
        Qvals = zeros(d,n_iters);
        Pvals = zeros(d,n_iters);
    end

    gradV!(gradV0, Q₀);
    for j in 1:n_iters
        @. Phalf = P₀ - 0.5 * Δt * gradV0;
        @. Q₀ = Q₀ + Δt * Phalf/M;
        gradV!(gradV0,Q₀);
        @. P₀ = Phalf - 0.5 * Δt * gradV0;

        if return_trajectory
            @. Qvals[:,j] = Q₀;
            @. Pvals[:,j] = P₀;
        end

    end

    if return_trajectory
        return Qvals, Pvals
    else
        return Q₀, P₀
    end

end

"""
    Verlet!(q₀,p₀, gradV!, M, Δt, n_iters, return_trajectory=true)
Run an in place Verlet integrator.
"""
function Verlet!(Q₀,P₀, gradV!, M, Δt, n_iters)

    gradV0 = similar(Q₀);
    Phalf = similar(P₀);

    gradV!(gradV0, Q₀);
    for j in 1:n_iters
        @. Phalf = P₀ - 0.5 * Δt * gradV0;
        @. Q₀ = Q₀ + Δt * Phalf/M;
        gradV!(gradV0,Q₀);
        @. P₀ = Phalf - 0.5 * Δt * gradV0;

    end
    Q₀, P₀
end

"""
    HMC(x₀, V, gradV!, β, M, Δt, nΔt, n_iters)
Run the HMC sampler
"""

function HMC(x₀, V, gradV!, β, M, Δt, nΔt, n_iters; return_trajectory = true)

    Q₀ = copy(x₀);
    P₀ = similar(x₀);
    Qp = similar(x₀);
    Pp = similar(x₀);
    gradVp = similar(x₀);
    Phalf = similar(P₀);
    d = length(x₀);

    naccept = 0;

    if(return_trajectory)
        Xvals =zeros(d,n_iters);
        avals =zeros(n_iters);
    end

    V₀ = V(Q₀);
    for j in 1:n_iters
        @. P₀ = sqrt(M/β) * randn();
        E₀ = V₀ + 0.5 * P₀⋅ (P₀ ./M);

        # run Verlet integrator
        Qp = copy(Q₀);
        Pp = copy(P₀);
        gradV!(gradVp, Qp);
        for k in 1:nΔt
            @. Phalf = Pp - 0.5 * Δt * gradVp;
            @. Qp = Qp + Δt * Phalf/M;
            gradV!(gradVp,Qp);
            @. Pp = Phalf - 0.5 * Δt * gradVp;
        end

        # compute energy
        Vp = V(Qp);
        Ep = Vp + 0.5 * Pp⋅ (Pp ./M);

        # accept/reject
        a = min(1, exp( β *(E₀-Ep)));
        if rand()<a
            naccept = naccept+1;
            @. Q₀ = Qp;
            V₀ = Vp;
        end

        if(return_trajectory)
            @. Xvals[:,j] = Q₀;
            avals[j] = naccept/j;
        end
    end

    if return_trajectory
        return Xvals, avals
    else
        return Q₀, naccept/n_iters
    end
end

"""
    HMC!(Q₀, V, gradV!, β, M, Δt, nΔt, n_iters)
In place HMC sampler
"""

function HMC!(Q₀, V, gradV!, β, M, Δt, nΔt, n_iters)

    P₀ = similar(x₀);
    Qp = similar(x₀);
    Pp = similar(x₀);
    gradVp = similar(x₀);
    Phalf = similar(P₀);
    d = length(x₀);

    V₀ = V(Q₀);
    for j in 1:n_iters
        @. P₀ = sqrt(M/β) * randn();
        E₀ = V₀ + 0.5 * P₀⋅ (P₀ ./M);

        # run Verlet integrator
        Qp = copy(Q₀);
        Pp = copy(P₀);
        gradV!(gradVp, Qp);
        for k in 1:nΔt
            @. Phalf = Pp - 0.5 * Δt * gradVp;
            @. Qp = Qp + Δt * Phalf/M;
            gradV!(gradVp,Qp);
            @. Pp = Phalf - 0.5 * Δt * gradVp;
        end

        # compute energy
        Vp = V(Qp);
        Ep = Vp + 0.5 * Pp⋅ (Pp ./M);

        # accept/reject
        a = min(1, exp( β *(E₀-Ep)));
        if rand()<a
            naccept = naccept+1;
            @. Q₀ = Qp;
            V₀ = Vp;
        end
    end
    Q₀
end

end # end module

max_stars_repo_name: cortner/JuBasicMD

max_issues_repo_name: cortner/JuBasicMD

max_forks_repo_name: cortner/JuBasicMD


text:
function MassCG(CG,Param)
OrdPoly=CG.OrdPoly;
nz=Param.Grid.nz;
M=zeros(CG.NumG,nz);
J = Param.cache.J
for iF=1:Param.Grid.NumFaces
  M[CG.Faces[iF].Glob,:]=M[CG.Faces[iF].Glob,:]+
    reshape(0.5*(J[:,:,1,iF,:]+J[:,:,2,iF,:]),
    (OrdPoly+1)*(OrdPoly+1),nz);
end
MW=0.5*(M[:,1:end-1]+M[:,2:end]);
return (M,MW)
end


text:
module SymSparseMatrixCSRTests

using Test
using SparseMatricesCSR
using SparseArrays
using LinearAlgebra

function test_csr(Bi,Tv,Ti)
  maxrows=5
  maxcols=6

  I = Ti[1,2,1,3,2,5,2,5]
  J = Ti[1,1,2,2,3,2,5,5]
  V = Tv[2,3,3,4,4,7,7,8]

  I_up = Ti[1,1,2,2,5]
  J_up = Ti[1,2,3,5,5]
  V_up = Tv[2,3,4,7,8]

  CSC = sparse(I,J,V)
  if Bi == 1
    CSR = symsparsecsr(I_up,J_up,V_up)
    @test CSR == CSC
    CSR = symsparsecsr(copy(I),copy(J),copy(V);symmetrize=true)
    @test CSR == CSC
  end

  CSR = symsparsecsr(Val(Bi),I_up,J_up,V_up)
  @test CSR == CSC
  CSR = symsparsecsr(Val(Bi),copy(I),copy(J),copy(V);symmetrize=true)
  @test CSR == CSC
  @test eltype(CSR) == Tv
  @test isa(CSR,SymSparseMatrixCSR{Bi,Tv,Ti})

  for i=1:size(CSR,1)
    for j=1:size(CSR,2)
      if (i,j) in zip(I,J)
         CSR[i,j] = eltype(V)(i+j)
         @test CSR[i,j] ≈ eltype(V)(i+j)
      else
         try
          CSR[i,j] = eltype(V)(i+j)
         catch e
          @test isa(e,ArgumentError)
         end
      end
    end
  end

  CSC = sparse(I,J,V,maxrows,maxcols)
  if Bi == 1
    CSR = symsparsecsr(I_up,J_up,V_up,maxrows,maxcols)
    @test CSR == CSC
    CSR = symsparsecsr(copy(I),copy(J),copy(V),maxrows,maxcols;symmetrize=true)
    @test CSR == CSC
  end

  CSR = symsparsecsr(Val(Bi),I_up,J_up,V_up,maxrows,maxcols)
  @test CSR == CSC
  CSR = symsparsecsr(Val(Bi),copy(I),copy(J),copy(V),maxrows,maxcols;symmetrize=true)
  @test CSR == CSC

  CSR2 = convert(typeof(CSR),CSR)
  @test CSR2 === CSR
  CSR2 = convert(SymSparseMatrixCSR{0,Float64,Int32},CSR)
  @test CSR2 == CSR

  @test size(CSR) == size(CSC)
  @test eltype(CSR) == Tv
  @test isa(CSR,SymSparseMatrixCSR{Bi,Tv,Ti})
  @test issparse(CSR)
  @test getBi(CSR) == Bi
  @test getoffset(CSR) == 1-Bi
  @test nnz(CSR) == nnz(CSR.uppertrian)
  @test length(nonzeros(CSR)) == nnz(CSR)
  @test nonzeros(CSR) === CSR.uppertrian.nzval
  @test colvals(CSR) === CSR.uppertrian.colval
  i,j,v = findnz(CSR)
  csr = symsparsecsr(Val(Bi),i,j,v,maxrows,maxcols)
  @test csr == CSR
  @test count(v->v>0,CSR) == count(v->v>0,nonzeros(CSR))

  x = rand(Tv,maxcols)
  y = Vector{Tv}(undef,maxrows)
  z = Vector{Tv}(undef,maxrows)
  mul!(y,CSR,x)
  mul!(z,CSC,x)
  @test y ≈ z

  mul!(y,CSR,x,1,2)
  mul!(z,CSC,x,1,2)
  @test y ≈ z

  @test CSR*x ≈ CSC*x

  out = LinearAlgebra.fillstored!(CSR,3.33)
  @test out === CSR
  LinearAlgebra.fillstored!(CSC,3.33)
  mul!(y,CSR,x)
  mul!(z,CSC,x)
  @test y ≈ z

  _CSR = copy(CSR)
  out = LinearAlgebra.rmul!(CSR,-1)
  @test out === CSR
  @test _CSR ≈ -1*CSR

end

for Bi in (0,1)
  for Tv in (Float32,Float64)
    for Ti in (Int32,Int64)
      test_csr(Bi,Tv,Ti)
    end
  end
end

end # module


text:
module Helper
using SpecialFunctions
using Distributions
import Random: AbstractRNG
import StatsFuns: logtwo, normlogpdf, normlogcdf

std_t_lcdf(x, df) = log(std_t_cdf(x, df))
function std_t_cdf(x, df)
  x_t = df / (x^2 + df)
  neg_cdf = 0.5 * beta_inc(0.5 * df, 0.5, x_t)[1]
  return x < 0 ? neg_cdf : 1 - neg_cdf
end

struct SkewT{Tdf<:Real, Tloc<:Real, Tscale<:Real, Tskew<:Real} <: ContinuousUnivariateDistribution
  df::Tdf
  loc::Tloc
  scale::Tscale
  skew::Tskew
end

function skewt_logpdf(df, loc, scale, skew, x)
  z = (x - loc) / scale
  u = skew * z * sqrt((df + 1) / (df + z ^ 2))
  kernel = logpdf(TDist(df), z) + std_t_lcdf(u, df + 1)
  # kernel = logpdf(TDist(df), z) + logcdf(TDist(df + 1), u)
  return kernel + log(2 / scale)
end

function Distributions.logpdf(d::SkewT{Tdf, Tloc, Tscale, Tskew}, x::Real) where {Tdf<:Real,
                                                                                Tloc<:Real,
                                                                                Tscale<:Real,
                                                                                Tskew<:Real}
  return skewt_logpdf(d.df, d.loc, d.scale, d.skew, x)
end

function Distributions.rand(rng::AbstractRNG,
                            d::SkewT{Tdf, Tloc, Tscale, Tskew}) where {Tdf<:Real,
                                                                       Tloc<:Real,
                                                                       Tscale<:Real,
                                                                       Tskew<:Real}
  w = rand(rng, Gamma(d.df/2, 2/d.df))
  z = rand(rng, truncated(Normal(0, sqrt(1/w)), 0, Inf))
  delta = d.skew / sqrt(1 + d.skew ^ 2)
  loc = d.loc + d.scale * z * delta
  scale = d.scale * sqrt(1 - delta ^ 2)

  return randn(rng) * scale + loc
end

function skewnormal_logpdf(loc, scale, skew, y)
  z = (y - loc) / scale
  return logtwo - log(scale) + normlogpdf(z) + normlogcdf(z * skew)
end

end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using Test, Astro

print("13.a Equatorial to ecliptical coordinates")
# convert RA/Dec to celestial longitude and latitude
# equ_to_ecl(ra, declin, obliquity)
# ra is usually in hours, so convert to radians
# Calculate the ecliptical coordinate the star Pollux whose equatorial coordinates are
# ra = a2000 = 7h 45m 12.946s = 116.328942deg, dec = d2000 = +28d01m34.26
# L should be 113.21563 degrees2radians
# B should be +6.68417 degrees

L, B = equ_to_ecl(hms_to_fday(7,45,12.946) * 2 * pi, deg2rad(dms_to_d(28,1,34.26)), deg2rad(23.4392911))
# L = longitude
@test isapprox(rad2deg(L), 113.21563, atol=0.1)
# B = latitude
@test isapprox(rad2deg(B), 6.68417, atol=0.1)
println(" passed")

print("13.a1  Ecliptic to Equatorial")
ra, decl = ecl_to_equ(L, B, deg2rad(23.4392911))
# ra is right ascension in radians
ra_deg, decl_deg = map(rad2deg, (ra,decl))
@test isapprox(ra_deg, 116.328942, atol=0.1)
@test isapprox(decl_deg, 28.026183, atol=0.1)
ra_day = fday_to_hms(ra_deg)
ra_hms =  fday_to_hms(rad2deg(ra)/360)
@test ra_hms[1:2] == (7.0,45.0)
@test isapprox(ra_hms[3], 12.946, atol=0.1)
println(" passed")

print("13.a2 Ecliptical to equatorial coordinates")
# what are the right ascension and declination of a planet whose ecliptic coordinates are
# ecl longitude 139-41-10 and ecl latitude 4-52-31"
# RA should be 09h34m53.6s
# decl should be 10d32m14,2s
ra, decl = ecl_to_equ(deg2rad(dms_to_d(139, 41, 10)), deg2rad(dms_to_d(4, 52, 31)), deg2rad(23.4392911))
ra_deg, decl_deg = map(rad2deg, (ra, decl))
@test isapprox(ra_deg/15, 9.581551, atol=0.01)
@test isapprox(decl_deg, dms_to_d(19,32,14.2), atol=0.01)
println(" passed")

print("13 supplemental Ecliptical to equatorial coordinates (RA, decl)")
ra, decl = ecl_to_equ(deg2rad(dms_to_d(139, 41, 10)), deg2rad(dms_to_d(4, 52, 31)), deg2rad(23.4392911))
ra_deg, decl_deg = map(rad2deg, (ra, decl))
@test isapprox(ra_deg/15, 9.581534330766397, atol=0.01)
@test isapprox(decl_deg, 19.536745492338458, atol=0.01)
println(" passed")

print("Equatorial (HA/Dec) to horizontal coordinates (Alt/Azi)")
# need to set latitude
set_latitude(deg2rad(52.0))
println(" using latitude $(rad2deg(latitude))")

# what is the altitude and azimuth of a star whose hour-angle is 05h51m44s and declination is 23d13m10s,
# latitude 52N
# altitude should be 19d20m04s
# azimuth should be 283d16m16s

set_latitude(deg2rad(52))
h = hangle_to_dec_deg(5,51,44) # returns degrees
d = dms_to_d(23,13,10) # returns degrees
azi,alt = equ_to_horiz(deg2rad(h), deg2rad(d))
azi_deg = rad2deg(azi)
alt_deg = rad2deg(alt)
@test(d_to_dms(azi_deg)[1] == 283)
@test isapprox(d_to_dms(azi_deg)[2], 16, atol=0.1)
@test d_to_dms(azi_deg)[3] ≈  16 atol= 0.4# only accurate to 0.4
@test(d_to_dms(alt_deg)[1] == 19.0)
@test(d_to_dms(alt_deg)[2] == 20.0)
@test d_to_dms(alt_deg)[3] ≈  4 atol= 0.4# only accurate to 0.4 seconds???
println(" passed")

print("40.a topocentric right ascension and declination")
# calculate the topocentric RA and Dec of Mars on 2003 Aug 28 at 3h17m00 UT at Palomar,
# longitude 7h47m27s, altitude 1706
# Mars' equatorial coordinates: RA 339.530208, decl -15.771083
"""
      phi: geocentric radius in radians
        H : observer's altitude (above sea level) in meters
        L ; observer longitude (in radians) (hangle_to_dec_deg(7,47,27))
        ra : body right ascension
        decl : body declination in radians
        d : body distance (in AU)
        jd : observer's hour angle in radian

"""
ra, decl = geocentric_to_topocentric(deg2rad(33.356111), 1706, deg2rad(116.8625), deg2rad(339.530208), deg2rad(-15.771083), 0.37276, cal_to_jd(2003, 8, 28.136806))
# topocentric right ascension
@test isapprox(rad2deg(ra), 339.5, atol=.1)
# topocentric declination
@test isapprox(rad2deg(decl), -15.775, atol=.1)
println(" passed")