text:
# Functions for estimating information measures for continuous or discrete datasets.
# These are intended to maximise ease of use for setting different options. For
# maximum speed, it may be better to use the functions from Formulae.jl directly.

export discretize_values,
	get_probabilities,
	get_entropy,
	get_conditional_entropy,
	get_mutual_information,
	get_conditional_mutual_information,
	get_total_correlation,
	get_interaction_information,
	get_partial_information_decomposition,
	get_redundancy,
	get_cross_entropy

"""
    discretize_values(values_x; <keyword arguments>)

Assign continuous measurements to discrete bins.

# Arguments
* `values_x`: the data values.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
"""
function discretize_values(values_x...; mode = "uniform_width", number_of_bins = 0, get_number_of_bins = get_root_n)

	if number_of_bins == 0
		number_of_bins = get_number_of_bins(values_x...)
	end

	return get_frequencies(mode, number_of_bins, values_x...)
end

"""
    get_probabilities(estimator, frequencies; <keyword arguments>)

Estimate probabilities from a set of discrete values.

# Arguments:
* `estimator`: the entropy estimator.
* `frequencies`: the bin frequencies for the discretized data values.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_probabilities(estimator, frequencies; lambda = nothing, prior = 1)

	if estimator == "maximum_likelihood" || estimator == "miller_madow"
		probabilities = get_probabilities_maximum_likelihood(frequencies)
	elseif estimator == "shrinkage"
		probabilities = get_probabilities_shrinkage(frequencies, lambda)
	elseif estimator == "dirichlet"
		probabilities = get_probabilities_dirichlet(frequencies, prior)
	end

	return probabilities
end

"""
    get_entropy(values_x; <keyword arguments>)

Estimate entropy (or joint entropy, if more than one set of data values is given).

# Arguments:
* `values_x`: the data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `discretized=false`: whether the data values are already discretized.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_entropy(values_x...; estimator = "maximum_likelihood", base = 2, mode = "uniform_width", number_of_bins = 0,
	get_number_of_bins = get_root_n, discretized = false, lambda = nothing, prior = 1)

	# If discretized, values_x should be one array of frequencies
	# (because higher dimensional frequencies cannot be reconstructed
	# from multiple lower dimensional frequencies), so values_x... will
	# be ([f1, f2, f3, ...],)
	frequencies = discretized ?
		values_x[1] :
		discretize_values(values_x..., mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	probabilities = get_probabilities(estimator, frequencies, lambda = lambda, prior = prior)

	entropy = apply_entropy_formula(probabilities, base)

	if estimator == "miller_madow"
		entropy += (countnz(probabilities) - 1) / (2 * length(values_x[1]))
	end

	# Make sure we don't return -0.0
	return entropy + 0.0
end

"""
    get_conditional_entropy(values_x, values_y; <keyword arguments>)

Estimate conditional entropy of one set of values conditioned on another set of values.

# Arguments:
* `values_x`: the data values.
* `values_y`: the data values to be conditioned on.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_conditional_entropy(values_x, values_y; estimator = "maximum_likelihood", base = 2, mode = "uniform_width",
	number_of_bins = 0, get_number_of_bins = get_root_n, lambda = nothing, prior = 1)

	frequencies_xy = discretize_values(values_x, values_y, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_conditional_entropy(frequencies_xy; estimator = estimator, base = base, lambda = lambda, prior = prior)
end
"""
    get_conditional_entropy(xy; <keyword arguments>)

# Arguments:
* `xy`: the 2D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xy` is probabilities. If `false`, `xy` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_conditional_entropy(xy; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1)

	probabilities_xy = probabilities ? xy : get_probabilities(estimator, xy, lambda = lambda, prior = prior)
	probabilities_y = sum(probabilities_xy, 2)

	entropy_xy = apply_entropy_formula(probabilities_xy, base)
	entropy_y = apply_entropy_formula(probabilities_y, base)

	# TODO: either deprecate or add warning about inappropriate use
	if estimator == "miller_madow"
		entropy_xy += (countnz(probabilities_xy) - 1) / (2 * length(probabilities_xy))
		entropy_y += (countnz(probabilities_y) - 1) / (2 * length(probabilities_y))
	end

	return apply_conditional_entropy_formula(entropy_xy, entropy_y)
end

"""
    get_mutual_information(values_x, values_y; <keyword arguments>)

Estimate the mutual information between two sets of values.

# Arguments:
* `values_x`: the data values.
* `values_y`: a second set of data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_mutual_information(values_x, values_y; estimator = "maximum_likelihood", base = 2, mode = "uniform_width",
	number_of_bins = 0, get_number_of_bins = get_root_n, lambda = nothing, prior = 1)

	frequencies_xy = discretize_values(values_x, values_y, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_mutual_information(frequencies_xy; estimator = estimator, base = base, lambda = lambda, prior = prior)
end
"""
    get_mutual_information(xy; <keyword arguments>)

# Arguments:
* `xy`: the 2D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xy` is probabilities. If `false`, `xy` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_mutual_information(xy; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1)

	probabilities_xy = probabilities ? xy : get_probabilities(estimator, xy, lambda = lambda, prior = prior)
	probabilities_x = sum(probabilities_xy, 1)
	probabilities_y = sum(probabilities_xy, 2)

	# TODO: either deprecate or add warning about inappropriate use
	if estimator == "miller_madow"
		entropy_xy = apply_entropy_formula(probabilities_xy, base)
		entropy_x = apply_entropy_formula(probabilities_x, base)
		entropy_y = apply_entropy_formula(probabilities_y, base)
		entropy_xy += (countnz(probabilities_xy) - 1) / (2 * length(probabilities_xy))
		entropy_x += (countnz(probabilities_x) - 1) / (2 * length(probabilities_x))
		entropy_y += (countnz(probabilities_y) - 1) / (2 * length(probabilities_y))
		return apply_mutual_information_formula(entropy_x, entropy_y, entropy_xy)
	end

	return apply_mutual_information_formula(probabilities_xy, probabilities_x, probabilities_y, base)
end

"""
    get_conditional_mutual_information(values_x, values_y, values_z; <keyword arguments>)

Estimate the conditional mutual information between two sets of values, conditioned on a third.

# Arguments:
* `values_x`: the data values.
* `values_y`: a second set of data values.
* `values_z`: the data values to be conditioned on.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_conditional_mutual_information(values_x, values_y, values_z; estimator = "maximum_likelihood", base = 2, mode = "uniform_width",
	number_of_bins = 0, get_number_of_bins = get_root_n, lambda = nothing, prior = 1)

	frequencies_xyz = discretize_values(values_x, values_y, values_z, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_conditional_mutual_information(frequencies_xyz; estimator = estimator, base = base, lambda = lambda, prior = prior)
end
"""
    get_conditional_mutual_information(xyz; <keyword arguments>)

# Arguments:
* `xyz`: the 3D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xyz` is probabilities. If `false`, `xyz` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_conditional_mutual_information(xyz; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1)

	probabilities_xyz = probabilities ? xyz : get_probabilities(estimator, xyz, lambda = lambda, prior = prior)
	probabilities_xz = sum(probabilities_xyz, 2)
	probabilities_yz = sum(probabilities_xyz, 1)
	probabilities_z = sum(probabilities_xz, 1) # xz not a typo

	entropy_xyz = apply_entropy_formula(probabilities_xyz, base)
	entropy_xz = apply_entropy_formula(probabilities_xz, base)
	entropy_yz = apply_entropy_formula(probabilities_yz, base)
	entropy_z = apply_entropy_formula(probabilities_z, base)

	# TODO: either deprecate or add warning about inappropriate use
	if estimator == "miller_madow"
		entropy_xyz += (countnz(probabilities_xyz) - 1) / (2 * length(probabilities_xyz))
		entropy_xz += (countnz(probabilities_xz) - 1) / (2 * length(probabilities_xz))
		entropy_yz += (countnz(probabilities_yz) - 1) / (2 * length(probabilities_yz))
		entropy_z += (countnz(probabilities_z) - 1) / (2 * length(probabilities_z))
	end

	return apply_conditional_mutual_information_formula(entropy_xz, entropy_yz, entropy_xyz, entropy_z)
end

"""
    get_interaction_information(values_x, values_y, values_z; <keyword arguments>)

Estimate the interaction information between three sets of values.

# Arguments:
* `values_x`: the data values.
* `values_y`: a second set of data values.
* `values_z`: a third set of data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_interaction_information(values_x, values_y, values_z; estimator = "maximum_likelihood", base = 2, mode = "uniform_width",
	number_of_bins = 0, get_number_of_bins = get_root_n, lambda = nothing, prior = 1)

	frequencies_xyz = discretize_values(values_x, values_y, values_z, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_interaction_information(frequencies_xyz; estimator = estimator, base = base, lambda = lambda, prior = prior)
end
"""
    get_interaction_information(xyz; <keyword arguments>)

# Arguments:
* `xyz`: the 3D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xyz` is probabilities. If `false`, `xyz` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_interaction_information(xyz; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1)

	probabilities_xyz = probabilities ? xyz : get_probabilities(estimator, xyz, lambda = lambda, prior = prior)
	probabilities_xy = probabilities_xz = sum(probabilities_xyz, 3)

	conditional_mutual_information = get_conditional_mutual_information(probabilities_xyz, estimator = estimator, base = base, probabilities = true, lambda = lambda, prior = prior)
	mutual_information = get_mutual_information(probabilities_xy, estimator = estimator, base = base, probabilities = true, lambda = lambda, prior = prior)

	return apply_interaction_information_formula(conditional_mutual_information, mutual_information)
end

"""
    get_total_correlation(values_x, values_y, values_z; <keyword arguments>)

Estimate the total correlation between three sets of values.

# Arguments:
* `values_x`: the data values.
* `values_y`: a second set of data values.
* `values_z`: a third set of data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_total_correlation(values_x, values_y, values_z; estimator = "maximum_likelihood", base = 2, mode = "uniform_width", number_of_bins = 0,
	get_number_of_bins = get_root_n, lambda = nothing, prior = 1)

	frequencies_xyz = discretize_values(values_x, values_y, values_z, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_total_correlation(frequencies_xyz; estimator = estimator, base = base, lambda = lambda, prior = prior)
end
"""
    get_total_correlation(xyz; <keyword arguments>)

# Arguments:
* `xyz`: the 3D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xyz` is probabilities. If `false`, `xyz` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_total_correlation(xyz; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1)

	probabilities_xyz = probabilities ? xyz : get_probabilities(estimator, xyz, lambda = lambda, prior = prior)
	probabilities_xy = sum(probabilities_xyz, 3)
	probabilities_x = sum(probabilities_xy, 2)
	probabilities_y = sum(probabilities_xy, 1)
	probabilities_z = sum(sum(probabilities_xyz, 1), 2)

	entropy_xyz = apply_entropy_formula(probabilities_xyz, base)
	entropy_x = apply_entropy_formula(probabilities_x, base)
	entropy_y = apply_entropy_formula(probabilities_y, base)
	entropy_z = apply_entropy_formula(probabilities_z, base)

	return apply_total_correlation_formula(entropy_x, entropy_y, entropy_z, entropy_xyz)
end

"""
    get_partial_information_decomposition(values_x, values_y, values_z; <keyword arguments>)

Estimate the partial information decomposition between three sets of values.

Performance can be improved by setting `all_orientations`, `include_unique` and `include_synergy` according to the use case.

# Arguments:
* `values_x`: the data values.
* `values_y`: a second set of data values.
* `values_z`: a third set of data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equqivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
* `all_orientations=false`: whether to use each set of values as the target. If `false`, only `values_z` is the target.
* `include_unique=true`: whether to get the unique information.
* `include_synergy=true`: whether to get the synergy.
"""
function get_partial_information_decomposition(values_x, values_y, values_z; estimator = "maximum_likelihood", base = 2, mode = "uniform_width", number_of_bins = 0,
	get_number_of_bins = get_root_n, lambda = nothing, prior = 1, all_orientations = false, include_unique = true, include_synergy = true)

	frequencies_xyz = discretize_values(values_x, values_y, values_z, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_partial_information_decomposition(frequencies_xyz; estimator = estimator, base = base, lambda = lambda, prior = prior, all_orientations = all_orientations,
		include_unique = include_unique, include_synergy = include_synergy)
end
"""
    get_partial_information_decomposition(xyz; <keyword arguments>)

# Arguments:
* `xyz`: the 3D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xyz` is probabilities. If `false`, `xyz` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
* `all_orientations=false`: whether to use each set of values as the target. If `false`, only `values_z` is the target.
* `include_unique=true`: whether to get the unique information.
* `include_synergy=true`: whether to get the synergy.
"""
function get_partial_information_decomposition(xyz; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1, all_orientations = false,
	include_unique = true, include_synergy = true)

	probabilities_xyz = probabilities ? xyz : get_probabilities(estimator, xyz, lambda = lambda, prior = prior)

	probabilities_xz = sum(probabilities_xyz, 2)
	probabilities_yz = sum(probabilities_xyz, 1)
	probabilities_xy = sum(probabilities_xyz, 3)
	probabilities_x = sum(probabilities_xz, 3)
	probabilities_y = sum(probabilities_yz, 3)
	probabilities_z = sum(probabilities_xz, 1)

	# Needed for unique information and synergy
	if include_unique || include_synergy
		entropy_x = apply_entropy_formula(probabilities_x, base)
		entropy_y = apply_entropy_formula(probabilities_y, base)
		entropy_z = apply_entropy_formula(probabilities_z, base)
		entropy_xy = apply_entropy_formula(probabilities_xy, base)
		entropy_xz = apply_entropy_formula(probabilities_xz, base)
		entropy_yz = apply_entropy_formula(probabilities_yz, base)
		mutual_information_xy = apply_mutual_information_formula(entropy_x, entropy_y, entropy_xy)
	end

	# Needed for unique information only
	if include_unique
		mutual_information_xz = apply_mutual_information_formula(entropy_x, entropy_z, entropy_xz)
		mutual_information_yz = apply_mutual_information_formula(entropy_y, entropy_z, entropy_yz)
	end

	# Needed for synergy only
	if include_synergy
		entropy_xyz = apply_entropy_formula(probabilities_xyz, base)
		interaction_information = apply_interaction_information_formula(
			apply_conditional_mutual_information_formula(entropy_xz, entropy_yz, entropy_xyz, entropy_z), mutual_information_xy
		)
	end

	pid = Dict()

	if all_orientations
		redundancy_xy_z = apply_redundancy_formula(probabilities_xz, probabilities_yz, probabilities_x, probabilities_y, probabilities_z, (1, 2), base)
		redundancy_xz_y = apply_redundancy_formula(probabilities_xy, probabilities_yz, probabilities_x, probabilities_z, probabilities_y, (1, 3), base)
		redundancy_yz_x = apply_redundancy_formula(probabilities_xy, probabilities_xz, probabilities_y, probabilities_z, probabilities_x, (2, 3), base)

		pid["z"] = Dict(
			"redundancy" => redundancy_xy_z
		)
		pid["y"] = Dict(
			"redundancy" => redundancy_xz_y
		)
		pid["x"] = Dict(
			"redundancy" => redundancy_yz_x
		)

		if include_unique
			pid["z"]["unique_1"] = mutual_information_xz - redundancy_xy_z
			pid["z"]["unique_2"] = mutual_information_yz - redundancy_xy_z
			pid["y"]["unique_1"] = mutual_information_xy - redundancy_xz_y
			pid["y"]["unique_2"] = mutual_information_yz - redundancy_xz_y
			pid["x"]["unique_1"] = mutual_information_xy - redundancy_yz_x
			pid["x"]["unique_2"] = mutual_information_xz - redundancy_yz_x
		end

		if include_synergy
			pid["z"]["synergy"] = interaction_information + redundancy_xy_z
			pid["y"]["synergy"] = interaction_information + redundancy_xz_y
			pid["x"]["synergy"] = interaction_information + redundancy_yz_x
		end

		# Rounding errors may lead to slightly negative results
		for orientation in ["x", "y", "z"]
			for measure in ["redundancy", "unique_1", "unique_2", "synergy"]
				if haskey(pid[orientation], measure)
					pid[orientation][measure] = pid[orientation][measure] < 0 ? 0 : pid[orientation][measure]
				end
			end
		end

	else

		redundancy = apply_redundancy_formula(probabilities_xz, probabilities_yz, probabilities_x, probabilities_y, probabilities_z, (1, 2), base)
		pid["redundancy"] = redundancy

		if include_unique
			unique_x = mutual_information_xz - redundancy
			unique_y = mutual_information_yz - redundancy
			# Rounding errors may lead to slightly negative results
			pid["unique_1"] = unique_x < 0 ? 0 : unique_x
			pid["unique_2"] = unique_y < 0 ? 0 : unique_y
		end

		if include_synergy
			synergy = interaction_information + redundancy
			# Rounding errors may lead to slightly negative results
			pid["synergy"] = synergy < 0 ? 0 : synergy
		end

	end

	return pid
end

"""
    get_redundancy(values_x, values_y, values_z; <keyword arguments>)

Estimate the redundancy between three sets of values.

# Arguments:
* `values_x`: the data values.
* `values_y`: a second set of data values.
* `values_z`: a third set of data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equqivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
* `all_orientations=false`: whether to use each set of values as the target. If `false`, only `values_z` is the target.
"""
function get_redundancy(values_x, values_y, values_z; estimator = "maximum_likelihood", base = 2, mode = "uniform_width", number_of_bins = 0,
	get_number_of_bins = get_root_n, lambda = nothing, prior = 1, all_orientations = false)

	frequencies_xyz = discretize_values(values_x, values_y, values_z, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

	return get_redundancy(frequencies_xyz; estimator = estimator, base = base, lambda = lambda, prior = prior, all_orientations = all_orientations)
end
"""
    get_redundancy(xyz; <keyword arguments>)

# Arguments:
* `xyz`: the 3D frequencies or probabilities.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `probabilities=false`: whether `xyz` is probabilities. If `false`, `xyz` is frequencies.
* `lambda=nothing`: the shrinkage instensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
* `all_orientations=false`: whether to use each set of values as the target. If `false`, only `values_z` is the target.
"""
function get_redundancy(xyz; estimator = "maximum_likelihood", base = 2, probabilities = false, lambda = nothing, prior = 1, all_orientations = false)
	probabilities_xyz = probabilities ? xyz : get_probabilities(estimator, xyz, lambda = lambda, prior = prior)

	probabilities_xz = sum(probabilities_xyz, 2)
	probabilities_yz = sum(probabilities_xyz, 1)
	probabilities_xy = sum(probabilities_xyz, 3)
	probabilities_x = sum(probabilities_xz, 3)
	probabilities_y = sum(probabilities_yz, 3)
	probabilities_z = sum(probabilities_xz, 1)

	if all_orientations
		redundancy_xy_z = apply_redundancy_formula(probabilities_xz, probabilities_yz, probabilities_x, probabilities_y, probabilities_z, (1, 2), base)
		redundancy_xz_y = apply_redundancy_formula(probabilities_xy, probabilities_yz, probabilities_x, probabilities_z, probabilities_y, (1, 3), base)
		redundancy_yz_x = apply_redundancy_formula(probabilities_xy, probabilities_xz, probabilities_y, probabilities_z, probabilities_x, (2, 3), base)

		return (redundancy_xy_z, redundancy_xz_y, redundancy_yz_x)
	else
		return apply_redundancy_formula(probabilities_xz, probabilities_yz, probabilities_x, probabilities_y, probabilities_z, (1, 2), base)
	end
end

"""
    get_cross_entropy(values_x, values_y; <keyword arguments>)
Estimate the cross-entropy between two sets of values.
# Arguments:
* `values_x`: the first set of data values.
* `values_y`: the second set of data values.
* `estimator="maximum_likelihood"`: the entropy estimator.
* `base=2`: the base of the logarithm, equivalent to the units of information.
* `mode="uniform_width"`: the discretization algorithm.
* `number_of_bins=0`: the number of bins (will be overridden if `mode` is `"bayesian_blocks"`).
* `get_number_of_bins=get_root_n`: the method for calculating the number of bins (only called if `number_of_bins` is `0`).
* `discretized=false`: whether the data values are already discretized.
* `lambda=nothing`: the shrinkage intensity, only used if `estimator` is `"shrinkage"`.
* `prior=1`: the Dirichlet prior, only used if `estimator` is `"dirichlet"`.
"""
function get_cross_entropy(values_x, values_y; estimator = "maximum_likelihood", base = 2, mode = "uniform_width", number_of_bins = 0, get_number_of_bins = get_root_n, discretized = false, lambda = nothing, prior = 1)
    frequency_x = discretized ? values_x : discretize_values(values_x, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)
    frequency_y = discretized ? values_y : discretize_values(values_y, mode = mode, number_of_bins = number_of_bins, get_number_of_bins = get_number_of_bins)

    probability_x = get_probabilities(estimator, frequency_x, lambda = lambda, prior = prior)
    probability_y = get_probabilities(estimator, frequency_y, lambda = lambda, prior = prior)

    cross_entropy = apply_cross_entropy_formula(probability_x, probability_y, base)

    if estimator == "miller_madow"
        println("WARNING: Miller-Madow correction not implemented for the cross-entropy. ")
        println("The calculation was performed without applying the correction.")
    end

    return cross_entropy
end


text:
##Testing the Generation of Kronecker Graph
include("KronGraph500NoPerm.jl")
using PyPlot, SparseArrays

# Set scale, average degree, and determine number of vertices and edges
SCALE = 12
EdgesPerVertex = 16
Nmax = 2 .^ SCALE
M = EdgesPerVertex .* Nmax

println("Maximum number of vertices: "*string(Nmax))
println("Number of edges: "*string(M))

# Generate graph and create adjacency matrix
v1,v2 = KronGraph500NoPerm(SCALE,EdgesPerVertex)
A = sparse(v1,v2,1)

# Visualize graph with spy plot
figure()
spy(convert(Array,A))

# Plot degree distribution
figure()
loglog(convert(Array,OutDegree(A)'),"o")
xlabel("degree")
ylabel("count")

ext: jl


text:
"""
    NonlinearContact{T,N} <: Contact{T,N}

    contact object for impact and friction with a nonlinear friction cone

    friction_coefficient: value of friction coefficient
    contact_tangent: mapping from world frame to surface tangent frame
    contact_normal: inverse/complement of contact_tangent
    contact_origin: position of contact on Body relative to center of mass
    contact radius: radius of contact
"""
mutable struct NonlinearContact{T,N} <: Contact{T,N}
    friction_coefficient::T
    friction_parameterization::SMatrix{2,2,T,4}
    collision::Collision{T,2,3,6}
end

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, contact::NonlinearContact)
#     summary(io, contact)
#     println(io, "")
#     println(io, "friction_coefficient:      "*string(contact.friction_coefficient))
#     println(io, "friction_parameterization: "*string(contact.friction_parameterization))
#     println(io, "collision:                 "*string(contact.collision))
# end

function NonlinearContact(body::Body{T}, normal::AbstractVector, friction_coefficient;
    contact_origin=szeros(T, 3),
    contact_radius=0.0) where T

    # contact directions
    V1, V2, V3 = orthogonal_columns(normal)
    A = [V1 V2 V3]
    Ainv = inv(A)
    contact_normal = Ainv[3, SA[1; 2; 3]]'
    contact_tangent = Ainv[SA[1; 2], SA[1; 2; 3]]

    # friction parameterization
    parameterization = SA{T}[
         1.0  0.0
         0.0  1.0
    ]

    # collision
    collision = SphereHalfSpaceCollision(contact_tangent, contact_normal, SVector{3}(contact_origin), contact_radius)

    NonlinearContact{T,8}(friction_coefficient, parameterization, collision)
end

function constraint(mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}
    # contact model
    model = contact.model

    # parent
    pbody = get_body(mechanism, contact.parent_id)
    xp, vp, qp, ωp = next_configuration_velocity(pbody.state, mechanism.timestep)

    # child
    cbody = get_body(mechanism, contact.child_id)
    xc, vc, qc, ωc = next_configuration_velocity(cbody.state, mechanism.timestep)

    # distance
    d = distance(model.collision, xp, qp, xc, qc)

    # relative tangential velocity
    vt = relative_tangential_velocity(model, xp, qp, vp, ωp, xc, qc, vc, ωc)

    # unpack contact variables
    γ = contact.impulses[2]
    s = contact.impulses_dual[2]

    SVector{N½,T}(
        d - s[1],
        model.friction_coefficient * γ[1] - γ[2],
        (model.friction_parameterization * vt - s[@SVector [3,4]])...)
end

function constraint_jacobian(contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}
    friction_coefficient = contact.model.friction_coefficient
    γ = contact.impulses[2] + REG * neutral_vector(contact.model) # TODO need to check this is legit
    s = contact.impulses_dual[2] + REG * neutral_vector(contact.model) # TODO need to check this is legit

    ∇s1 = [γ[SA[1]]; szeros(T,3)]'
    ∇s2 = [szeros(T,3,1) cone_product_jacobian(γ[SA[2,3,4]])]
    ∇s3 = Diagonal(SVector{4,T}(-1.0, 0.0, -1.0, -1.0))
    ∇s = [∇s1; ∇s2; ∇s3]

    ∇γ1 = [s[SA[1]]; szeros(T, 3)]'
    ∇γ2 = [szeros(T,3,1) cone_product_jacobian(s[SA[2,3,4]])]
    ∇γ3 = SA[0.0                   0.0 0.0 0.0;
             friction_coefficient -1.0 0.0 0.0;
             0.0                   0.0 0.0 0.0;
             0.0                   0.0 0.0 0.0;]
    ∇γ = [∇γ1; ∇γ2; ∇γ3]

    return [∇s ∇γ]
end

neutral_vector(::NonlinearContact{T,N}) where {T,N} = [sones(T, 2); szeros(T, Int(N/2) - 2)]

cone_degree(::NonlinearContact) = 2


text:
include("api.jl")

function WMATA_auth(SubscriptionKey::String)
    url = "https://api.wmata.com/Misc/Validate"
    subscription_key = Dict("api_key" => SubscriptionKey)

    try
        r = request("GET", url, subscription_key)
    catch
        @error "Invalid Subscription Key. Please refer to your WMATA account."
    end

    global wmata = wmata_API(
        SubscriptionKey,
        "https://api.wmata.com/Rail.svc/json/jStations",
        "https://api.wmata.com/Rail.svc/json/jStationTimes",
        "https://api.wmata.com/StationPrediction.svc/json/GetPrediction/",
        "https://api.wmata.com/Rail.svc/json/jPath?",
        "https://api.wmata.com/Rail.svc/json/jSrcStationToDstStationInfo",
        "https://api.wmata.com/TrainPositions/TrainPositions?contentType=json",
        "https://api.wmata.com/Incidents.svc/json/Incidents"
    )

    println("Authentication complete.")
end


text:
using GasDispersion
using Test

include("test_scenarios.jl")

@testset "GasDispersion.jl tests" begin
    @test_throws MethodError Scenario()

    @test_throws MethodError plume()
    @test_throws ErrorException plume(test_scenario, model=:error)
    @test_throws ErrorException plume(bad_class)

    @test_throws MethodError puff()
    @test_throws ErrorException puff(test_scenario, model=:error)
    @test_throws ErrorException puff(bad_class)

end

# testing utilities
include("util_tests.jl")

# testing specific models
include("gaussian_tests.jl")
include("britter_tests.jl")


text:
module JuMPDEA

import JuMP
import MathProgBase

#type DEAModel
#  model::JuMP.Model

#  limitSize::Int#variable size of smaller LP problem
#  #option::Symbol
#  #increase

#  subset::Vector
#  pns::Vector#samples not drawing to subsets

#  originConsMatrix::SparseMatrixCSC#constraint matrix of orignal LP problem (the bogger one)
#  fullOriginConMatrix::Matrix(Float64)
#  transposeOriginMatrix::Matrix(Float64)

#  #data of smaller LP problem
#  smallConsMatrix::SparseMatrixCSC
#  smallVarLB::Vector(Float64)
#  smallVarUB::Vector(Float64)
#  smallConsLB::Vector(Float64)
#  smallConsUB::Vector(Float64)

#  smallObjCoeff::Vector(Float64)
#  smallSense::Symbol

#  optStatus::Bool
#  iterations::Int

#  # DEA data
#  inputNum::Int
#  outputNum::Int
#end

#function DEAModel(m::JuMP.Model; size=300)
#  objCoeff, varLB, varUB = JuMP.prepProblemBounds(m)

#  DEAModel(m, size, )
#end
export solveDEA

function solveDEA(m::JuMP.Model; limitSize=300)
  ##############  get original model data  ##############
  originalConstrMatrix = JuMP.prepConstrMatrix(m)
  originalObjCoeff, originalConstrLB, originalConstrUB = JuMP.prepProblemBounds(m)
  originalVarLB = m.colLower
  originalVarUB = m.colUpper
  originalSense = m.objSense

  originalVarNum = m.numCols

  ###############  Standardize the model  ###############
  thetaPosition = 1
  for i=1:originalVarNum
    if originalVarLB[i] == -Inf
      thetaPosition = i
    end
  end

  tempMatrix = originalConstrMatrix[:, thetaPosition]
  originalConstrMatrix[:, thetaPosition] = originalConstrMatrix[:, originalVarNum]
  originalConstrMatrix[:, originalVarNum] = tempMatrix

  temp = originalVarLB[thetaPosition]
  originalVarLB[thetaPosition] = originalVarLB[originalVarNum]
  originalVarLB[originalVarNum] = temp

  temp = originalVarUB[thetaPosition]
  originalVarUB[thetaPosition] = originalVarUB[originalVarNum]
  originalVarUB[originalVarNum] = temp

  temp = originalObjCoeff[thetaPosition]
  originalObjCoeff[thetaPosition] = originalObjCoeff[originalVarNum]
  originalObjCoeff[originalVarNum] = temp

  ##################  random sampling  ##################
  fullOriginalConstrMatrix = full(originalConstrMatrix)
  transOriginalConstrMatrix = transpose(fullOriginalConstrMatrix)
  set, pns, smallConstrMatrix = sampling(limitSize, originalConstrMatrix)

  ################  build our model data  ################
  smallVarLB = reCreatV(originalVarLB, limitSize)
  smallVarUB = reCreatV(originalVarUB, limitSize)
  smallConstrLB = originalConstrLB
  smallConstrUB = originalConstrUB
  smallObjCoeff = reCreatV(originalObjCoeff, limitSize)
  smallSense = m.objSense

  ########################  loop  ########################
  optStatus = false
  iterations = 0
  overNum = 0
  lm = MathProgBase.model(m.solver)

  while optStatus == false
    iterations = iterations + 1
    optStatus = true

    lm = MathProgBase.model(m.solver)
    MathProgBase.loadproblem!(lm, smallConstrMatrix, smallVarLB, smallVarUB, smallObjCoeff, smallConstrLB, smallConstrUB, smallSense)
    MathProgBase.updatemodel!(lm)
    MathProgBase.optimize!(lm)

    if MathProgBase.status(lm) == :Infeasible
      optStatus, overNum = checkOpt(MathProgBase.getinfeasibilityray(lm), pns, transOriginalConstrMatrix, smallSense)
    else
      optStatus, overNum = checkOpt(MathProgBase.getconstrduals(lm), pns, transOriginalConstrMatrix, smallSense)
    end

    if optStatus == false
      reSampling(overNum, iterations, set, pns, smallConstrMatrix, originalConstrMatrix)
    end
  end

  m.objVal = MathProgBase.getobjval(lm)
  m.objVal += m.obj.aff.constant

end

#################  function: sampling  #################
function sampling(subsetSize::Int, originConsMatrix::SparseMatrixCSC)
  rowNum = size(originConsMatrix)[1]
  colNum = size(originConsMatrix)[2]
  pns = [1:colNum-1]
  set = Array(Int, subsetSize-1)
  randNum = 0
  smallConsMatrix = originConsMatrix[1:rowNum, 1:subsetSize]
  temp = 0.0

  for i = 1:rowNum
    smallConsMatrix[i,subsetSize] = originConsMatrix[i,colNum]
  end

  for i = 1:subsetSize-1
    randNum = rand(1:colNum-i)
    set[i] = pns[randNum]
    temp = pns[randNum]
    pns[randNum] = pns[colNum-i]
    pns[colNum-i] = temp

    for j = 1:rowNum
      smallConsMatrix[j, i] = originConsMatrix[j, set[i]]
    end
  end

  pns = pns[1:colNum-subsetSize]

  return set, pns, smallConsMatrix
end

######  function: recreat varLB, varUB, objCoeff  ######
function reCreatV(v::Vector{Float64}, subsetSize::Int)
  newV = Array(Float64, subsetSize)
  newV[1:subsetSize-1] = v[1:subsetSize-1]
  newV[subsetSize] = v[length(v)]
  return newV
end

#############  function: check optimality  #############
function checkOpt(uv::Vector{Float64}, pns::Vector{Int}, data::Matrix{Float64}, sense::Symbol)
  optimal = true
  overNum = 0
  checkValue = 0.0

  if sense == :Min
    checkValue, overNum = findmax(data[pns, :]*uv)
    optimal = checkValue <= 0
  else
    checkValue, overNum = findmin(data[pns, :]*uv)
    optimal = checkValue >= 0
  end

  optimal, overNum
end

################  function: reSampling  ################
function reSampling(overNum::Int, iteration::Int, set::Vector{Int}, pns::Vector{Int},
                    smallConsMatrix::SparseMatrixCSC, originConsMatrix::SparseMatrixCSC)
  iteration = mod(iteration, length(set))
  if iteration == 0
    iteration = length(set)
  end

  for i = 1:size(smallConsMatrix)[1]
    smallConsMatrix[i, iteration] = originConsMatrix[i, pns[overNum]]
  end

  temp = set[iteration]
  set[iteration] = pns[overNum]
  pns[overNum] = temp
end

end


text:
using Dates: Millisecond
using Random: randstring
import Base: getindex, convert
using .LazyTrees: LazyTree, data, children

mutable struct GeometryData{G <: Union{AbstractGeometry, AbstractMesh}, C <: Colorant, T <: Transformation}
    geometry::G
    color::C
    transform::T
end

GeometryData(g, c::Colorant) = GeometryData(g, c, IdentityTransformation())
GeometryData(g, t::Transformation) = GeometryData(g, RGBA(1, 0, 0, 0.5), t)

GeometryData(g::Union{AbstractGeometry, AbstractMesh}) =
    GeometryData(g, RGBA(1, 0, 0, 0.5))
# convert(::Type{G}, g::Union{AbstractGeometry, AbstractMesh}) where {G <: GeometryData} = GeometryData(g, RGBA(1, 0, 0, 0.5))

mutable struct VisData
    transform::Transformation
    geometries::Vector{GeometryData}
end

VisData() = VisData(IdentityTransformation(), GeometryData[])

const Path = Vector{Symbol}

struct CommandQueue
    delete::Set{Path}
    setgeometry::Set{Path}
    settransform::Set{Path}

    CommandQueue() = new(Set{Path}(), Set{Path}(), Set{Path}())
end

isempty(queue::CommandQueue) = isempty(queue.delete) && isempty(queue.setgeometry) && isempty(queue.settransform)

function empty!(queue::CommandQueue)
    empty!(queue.delete)
    empty!(queue.setgeometry)
    empty!(queue.settransform)
end

const VisTree = LazyTree{Symbol, VisData}

mutable struct CoreVisualizer
    lcm::LCM
    client_id::String
    tree::VisTree
    queue::CommandQueue
    publish_immediately::Bool

    function CoreVisualizer(lcm::LCM=LCM())
        client_id = "jl_$(randstring())"  # 10^14 possibilities
        vis = new(lcm, client_id, VisTree(), CommandQueue(), true)
        function handle_msg(channel, msg)
            try
                onresponse(vis, msg)
            catch e
                warn("""
An error ocurred while handling the viewer response:
    error: $e
    response: $msg
""")
            end
        end
        sub = subscribe(lcm, response_channel(vis), handle_msg, Comms.CommsT)
        @async while true
            handle(lcm)
        end
        vis
    end
end

request_channel(vis::CoreVisualizer) = "DIRECTOR_TREE_VIEWER_REQUEST_<$(vis.client_id)>"
response_channel(vis::CoreVisualizer) = "DIRECTOR_TREE_VIEWER_RESPONSE_<$(vis.client_id)>"

function setgeometry!(vis::CoreVisualizer, path::AbstractVector)
    push!(vis.queue.setgeometry, path)
    settransform!(vis, path)
end

function setgeometry!(vis::CoreVisualizer, path::AbstractVector, geoms::AbstractVector)
    vis.tree[path].data.geometries = geoms
    setgeometry!(vis, path)
end

function addgeometry!(vis::CoreVisualizer, path::AbstractVector, geoms::AbstractVector)
    append!(vis.tree[path].data.geometries, geoms)
    setgeometry!(vis, path)
end

setgeometry!(vis::CoreVisualizer, path::AbstractVector, geom) = setgeometry!(vis, path, [geom])
addgeometry!(vis::CoreVisualizer, path::AbstractVector, geom) = addgeometry!(vis, path, [geom])

setgeometry!(vis::CoreVisualizer, path::AbstractVector, geom::AbstractGeometry) =
    setgeometry!(vis, path, GeometryData(geom))
addgeometry!(vis::CoreVisualizer, path::AbstractVector, geom::AbstractGeometry) =
    addgeometry!(vis, path, GeometryData(geom))

function settransform!(vis::CoreVisualizer, path::AbstractVector)
    push!(vis.queue.settransform, path)
    if vis.publish_immediately
        publish!(vis)
    end
end

function settransform!(vis::CoreVisualizer, path::AbstractVector, tform)
    vis.tree[path].data.transform = tform
    settransform!(vis, path)
end

function delete!(vis::CoreVisualizer, path::AbstractVector)
    delete!(vis.tree, path)
    push!(vis.queue.delete, path)
    if vis.publish_immediately
        publish!(vis)
    end
end

function publish!(vis::CoreVisualizer)
    if !isempty(vis.queue)
        data = serialize(vis, vis.queue)
        msg = to_lcm(data)
        publish(vis.lcm, request_channel(vis), msg)
        empty!(vis.queue)
    end
end

function onresponse(vis::CoreVisualizer, msg)
    data = JSON.parse(IOBuffer(msg.data))
    if data["status"] == 0
        empty!(vis.queue)
    elseif data["status"] == 1
        for path in LazyTrees.descendants(vis.tree)
            push!(vis.queue.setgeometry, path)
            push!(vis.queue.settransform, path)
        end
        publish!(vis)
    else
        error("unhandled: $data")
    end
end

function to_lcm(data::AbstractDict)
    jsondata = JSON.json(data)
    Comms.CommsT(
        data["utime"],
        "treeviewer_json",
        1,
        0,
        length(jsondata),
        codeunits(jsondata))
end

struct Visualizer
    core::CoreVisualizer
    path::Vector{Symbol}

    Visualizer(lcm::LCM=LCM()) = new(CoreVisualizer(lcm), Symbol[])
    Visualizer(core::CoreVisualizer, path::AbstractVector) = new(core, path)
end

show(io::IO, vis::Visualizer) = print(io, "Visualizer with path prefix $(vis.path) using LCM $(vis.core.lcm)")

function setgeometry!(vis::Visualizer)
    setgeometry!(vis.core, vis.path)
    vis
end

function setgeometry!(vis::Visualizer, geom)
    setgeometry!(vis.core, vis.path, geom)
    vis
end

function addgeometry!(vis::Visualizer, geom)
    addgeometry!(vis.core, vis.path, geom)
    vis
end

settransform!(vis::Visualizer) = settransform!(vis.core, vis.path)
settransform!(vis::Visualizer, transform) = settransform!(vis.core, vis.path, transform)
delete!(vis::Visualizer) = delete!(vis.core, vis.path)
publish!(vis::Visualizer) = publish!(vis.core)

getindex(vis::Visualizer, path::Symbol) = Visualizer(vis.core, vcat(vis.path, path))
getindex(vis::Visualizer, path::AbstractVector) = Visualizer(vis.core, vcat(vis.path, path))

function batch(func, vis::Visualizer)
    old_publish_flag = vis.core.publish_immediately
    try
        vis.core.publish_immediately = false
        func(vis)
        if old_publish_flag
            publish!(vis)
        end
    finally
        vis.core.publish_immediately = old_publish_flag
    end
end

# Old-style visualizer interface
function Visualizer(geom::GeometryData)
    vis = Visualizer()[:body1]
    setgeometry!(vis, geom)
    vis
end

Visualizer(geom::Union{AbstractGeometry, AbstractMesh}) = Visualizer(GeometryData(geom))


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 8cf4d8ca-84eb-11eb-22d2-255ce7237090
begin
	using MITgcmTools, PlutoUI, Printf
	exps=verification_experiments()
end

# ╔═╡ f588eaba-84ef-11eb-0755-bf1b85b2b561
begin
md"""# Standard MITgcm configurations

This notebook scans configuration folders of [MITgcm](https://mitgcm.readthedocs.io/en/latest/?badge=latest) within `MITgcm/verification` using [MITgcmTools.jl](https://gaelforget.github.io/MITgcmTools.jl/dev/). It then let's user inspect parameters interactively.

!!! tip
	For more on compiling and running a model configuration, please refer to the [examples deck](https://gaelforget.github.io/MITgcmTools.jl/dev/examples/) for more on that topic.

"""
end

# ╔═╡ 98b6621c-85ab-11eb-29d1-af0433598c6a
md"""## Select Model Configuration

$(TableOfContents())

!!! note
	If you use a live version of this notebook, selecting a different configuration from the list below will make the other notebook cells react (e.g. displayed contents). If you visualize an html version of this notebook, then cells wont react.
"""

# ╔═╡ a28f7354-84eb-11eb-1830-1f401bf2db97
@bind myexp Select([exps[i].configuration for i in 1:length(exps)],default="advect_cs")

# ╔═╡ f051e094-85ab-11eb-22d4-5bd61ac572a1
md"""## Select Parameter Group

Each configuration has various groups of parameters, often called namelists, depending on compiled packages. The dropdown menus below let you navigate the various parameter groups -- the result is show in the next section.

The chosen defaults, `data` and `PARM01`, are expected to be found in any MITgcm run directory. However, there can be many more parameters.
"""

# ╔═╡ f40e76c4-86d5-11eb-15b0-cd55d6cd1e65
md"""### Appendices

The following cells select Julia packages and perform basic operations.
"""

# ╔═╡ 168e178c-dd09-4e27-8cb6-fc0479a55f75
begin
	🏁 = "🏁"
	imgA="https://user-images.githubusercontent.com/20276764/111042787-12377e00-840d-11eb-8ddb-64cc1cfd57fd.png"
	imgB="https://user-images.githubusercontent.com/20276764/97648227-970b9780-1a2a-11eb-81c4-65ec2c87efc6.png"
	md"""$(Resource(imgB, :width => 120))"""
end

# ╔═╡ 4965715d-93ca-496b-8ab1-238e9c6e34b4
begin
	iexp=findall([exps[i].configuration==myexp for i in 1:length(exps)])[1]
	builddir=joinpath(MITgcm_path[1],"verification",myexp)
	rundir=joinpath(exps[iexp].folder,string(exps[iexp].ID),"run")
	!isdir(rundir) ? setup(exps[iexp]) : nothing
	🏁
end

# ╔═╡ d7f2c656-8512-11eb-2fdf-47a3e57a55e6
begin
	function list_namelist_files(pth)
		tmpA=readdir(pth)
		tmpA=tmpA[findall([length(tmpA[i])>3 for i in 1:length(tmpA)])]
		tmpA=tmpA[findall([tmpA[i][1:4]=="data" for i in 1:length(tmpA)])]
		tmpA=[tmpA[:];"eedata"]
	end

	dats=list_namelist_files(rundir)
	try
		@bind mydats Select([dats[i] for i in 1:length(dats)])
	catch e
		"Error: could not find any namelist in $(pth)"
	end
end

# ╔═╡ 348c692e-84fe-11eb-3288-dd0a1dedce90
begin
	fil=joinpath(rundir,mydats)
	nml=read(fil,MITgcm_namelist())
	🏁
end

# ╔═╡ ca7bb004-8510-11eb-379f-632c3b40723d
try
	@bind nmlgroup Select(String.(nml.groups))
catch e
	"Error: could not find any namelist in $(rundir)"
end

# ╔═╡ e73fda3a-f05a-49b4-a83d-e7b535467106
md"""## Browse Parameters

#

Now displaying 👉 **$myexp / $mydats : $nmlgroup**"""

# ╔═╡ 9bdb94da-8510-11eb-01a6-c9a1519baa68
begin
	inml=findall(nml.groups.==Symbol(nmlgroup))[1]
	🏁
end

# ╔═╡ e50726aa-86d3-11eb-0418-fff8fb79ef95
nml.params[inml]

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
MITgcmTools = "62725fbc-3a66-4df3-9000-e33e85b3a198"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[compat]
MITgcmTools = "~0.1.32"
PlutoUI = "~0.7.25"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.AWS]]
deps = ["Base64", "Compat", "Dates", "Downloads", "GitHub", "HTTP", "IniFile", "JSON", "MbedTLS", "Mocking", "OrderedCollections", "Retry", "Sockets", "URIs", "UUIDs", "XMLDict"]
git-tree-sha1 = "07d944e4d9946c2061f97c1564d1b7ae8ea8f189"
uuid = "fbe9abb3-538b-5e4e-ba9e-bc94f4f92ebc"
version = "1.61.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "abb72771fd8895a7ebd83d5632dc4b989b022b5b"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.2"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Blosc]]
deps = ["Blosc_jll"]
git-tree-sha1 = "575bdd70552dd9a7eaeba08ef2533226cdc50779"
uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
version = "0.7.2"

[[deps.Blosc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "91d6baa911283650df649d0aea7c28639273ae7b"
uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
version = "1.21.1+0"

[[deps.CFTime]]
deps = ["Dates", "Printf"]
git-tree-sha1 = "bca6cb6ee746e6485ca4535f6cc29cf3579a0f20"
uuid = "179af706-886a-5703-950a-314cd64e0468"
version = "0.1.1"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "49f14b6c56a2da47608fe30aed711b5882264d7a"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.9.11"

[[deps.CatViews]]
deps = ["Random", "Test"]
git-tree-sha1 = "23d1f1e10d4e24374112fcf800ac981d14a54b24"
uuid = "81a5f4ea-a946-549a-aa7e-2a7f63a27d31"
version = "1.0.0"

[[deps.ClimateModels]]
deps = ["AWS", "CFTime", "CSV", "DataFrames", "Dates", "Downloads", "Git", "NetCDF", "OrderedCollections", "Pkg", "Statistics", "Suppressor", "TOML", "Test", "UUIDs", "Zarr"]
git-tree-sha1 = "0e5e942b23049bca23fd66a08d549c234373dd1c"
uuid = "f6adb021-9183-4f40-84dc-8cea6f651bb0"
version = "0.1.21"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "cfdfef912b7f93e4b848e80b9befdf9e331bc05a"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.1"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DiskArrays]]
git-tree-sha1 = "6a50d800025a1664c99a8e819e0568c75e3ac0c7"
uuid = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
version = "0.2.12"

[[deps.Distances]]
deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.7"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[deps.ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[deps.EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "04d13bfa8ef11720c24e4d840c0033d145537df7"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.17"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Gettext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "8c14294a079216000a0bdca5ec5a447f073ddc9d"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.20.1+7"

[[deps.Git]]
deps = ["Git_jll"]
git-tree-sha1 = "d7bffc3fe097e9589145493c08c41297b457e5d0"
uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
version = "1.2.1"

[[deps.GitHub]]
deps = ["Base64", "Dates", "HTTP", "JSON", "MbedTLS", "Sockets", "SodiumSeal"]
git-tree-sha1 = "c8594dff1ed76e232d8063b2a2555335900af6f3"
uuid = "bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26"
version = "5.7.0"

[[deps.Git_jll]]
deps = ["Artifacts", "Expat_jll", "Gettext_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "33be385f3432a5a5b7f6965af9592d4407f3167f"
uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
version = "2.31.0+0"

[[deps.HDF5_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "fd83fa0bde42e01952757f01149dd968c06c4dba"
uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
version = "1.12.0+1"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "8d70835a3759cdd75881426fced1508bb7b7e1b6"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.1"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LRUCache]]
git-tree-sha1 = "d64a0aff6691612ab9fb0117b0995270871c5dfc"
uuid = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
version = "1.3.0"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Lz4_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
version = "1.9.3+0"

[[deps.MITgcmTools]]
deps = ["Artifacts", "ClimateModels", "DataFrames", "Dates", "LazyArtifacts", "MeshArrays", "NetCDF", "OrderedCollections", "Printf", "SparseArrays", "Suppressor", "UUIDs"]
git-tree-sha1 = "85fd18c07803b16af3e317e9568d526f5ca61853"
uuid = "62725fbc-3a66-4df3-9000-e33e85b3a198"
version = "0.1.32"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MeshArrays]]
deps = ["CatViews", "Dates", "Downloads", "LazyArtifacts", "NearestNeighbors", "Pkg", "Printf", "SparseArrays", "Statistics", "Unitful"]
git-tree-sha1 = "c5b9b98540a900934d9b531f8a153ce49016cec7"
uuid = "cb8c808f-1acf-59a3-9d2b-6e38d009f683"
version = "0.2.31"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NearestNeighbors]]
deps = ["Distances", "StaticArrays"]
git-tree-sha1 = "16baacfdc8758bc374882566c9187e785e85c2f0"
uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
version = "0.4.9"

[[deps.NetCDF]]
deps = ["DiskArrays", "Formatting", "NetCDF_jll"]
git-tree-sha1 = "23b0e32fde256a4e2e497e678abcf956ed26204b"
uuid = "30363a11-5582-574a-97bb-aa9a979735b9"
version = "0.11.3"

[[deps.NetCDF_jll]]
deps = ["Artifacts", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Pkg", "Zlib_jll", "nghttp2_jll"]
git-tree-sha1 = "0cf4d1bf2ef45156aed85c9ac5f8c7e697d9288c"
uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3"
version = "400.702.400+0"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PCRE2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "d7fa6237da8004be601e19bd6666083056649918"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "93cf0910f09a9607add290a3a2585aa376b4feb6"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.25"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "b7ff9f9ce50eab241e978cd975ad4ae113f5a41f"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.0"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Retry]]
git-tree-sha1 = "41ac127cd281bb33e42aba46a9d3b25cd35fc6d5"
uuid = "20febd7b-183b-5ae2-ac4a-720e7ce64774"
version = "0.4.1"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "244586bc07462d22aed0113af9c731f2a518c93e"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.10"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SodiumSeal]]
deps = ["Base64", "Libdl", "libsodium_jll"]
git-tree-sha1 = "80cef67d2953e33935b41c6ab0a178b9987b1c99"
uuid = "2133526b-2bfb-4018-ac12-889fb3908a75"
version = "0.1.1"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.1.0"

[[deps.Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Unitful]]
deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "0992ed0c3ef66b0390e5752fe60054e5ff93b908"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.9.2"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

[[deps.XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[deps.XMLDict]]
deps = ["EzXML", "IterTools", "OrderedCollections"]
git-tree-sha1 = "d9a3faf078210e477b291c79117676fca54da9dd"
uuid = "228000da-037f-5747-90a9-8195ccbf91a5"
version = "0.4.1"

[[deps.Zarr]]
deps = ["AWS", "Blosc", "CodecZlib", "DataStructures", "Dates", "DiskArrays", "HTTP", "JSON", "LRUCache", "OffsetArrays", "Pkg", "URIs"]
git-tree-sha1 = "7238cf588d2def313a65b63ea9bba07aa762f26b"
uuid = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"
version = "0.7.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.libsodium_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "848ab3d00fe39d6fbc2a8641048f8f272af1c51e"
uuid = "a9144af2-ca23-56d9-984f-0d03f7b5ccf8"
version = "1.0.20+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─f588eaba-84ef-11eb-0755-bf1b85b2b561
# ╟─98b6621c-85ab-11eb-29d1-af0433598c6a
# ╟─a28f7354-84eb-11eb-1830-1f401bf2db97
# ╟─f051e094-85ab-11eb-22d4-5bd61ac572a1
# ╟─d7f2c656-8512-11eb-2fdf-47a3e57a55e6
# ╟─ca7bb004-8510-11eb-379f-632c3b40723d
# ╟─e73fda3a-f05a-49b4-a83d-e7b535467106
# ╟─e50726aa-86d3-11eb-0418-fff8fb79ef95
# ╟─f40e76c4-86d5-11eb-15b0-cd55d6cd1e65
# ╟─8cf4d8ca-84eb-11eb-22d2-255ce7237090
# ╟─9bdb94da-8510-11eb-01a6-c9a1519baa68
# ╟─348c692e-84fe-11eb-3288-dd0a1dedce90
# ╟─4965715d-93ca-496b-8ab1-238e9c6e34b4
# ╟─168e178c-dd09-4e27-8cb6-fc0479a55f75
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    tricontourf(x, y, z, t, levels)

Draw a filled contour plot of unstructured triangular data. `x`, `y`, and `z` are
one-dimensional input arrays of the same length (the number of points). `t` is an integer
array of size `(3,nt)`, where `nt` is the number of triangles. The coordinates of the `i`th
vertex of the `j`th triangle are given by `(x[t[i,j]], y[t[i,j]])`. Function values at
vertices are specified by the `z` array, and contours levels are defined by `levels.`

`tricontourf` results in regions *between* the contour levels, such that specifying `n`
levels will result in `n-1` regions. If `levels` is an array, it will be augmented to
include the minimum and maximum values of the `z` data.
"""
function tricontourf(x, y, z, t, levels)
    m = TriMesh(x, y, t)
    epsmin = eps(minimum(z))
    epsmax = eps(maximum(z))
    if levels isa Integer
        levels = range(minimum(z)-epsmin, maximum(z)+epsmax, length=levels)
    else
        levels = copy(levels)
        minz,maxz = extrema(z)
        minimum(levels) >= minz && pushfirst!(levels, minz-epsmin)
        maximum(levels) <= maxz && push!(levels, maxz+epsmax)
    end
    tricontourf(m, z, levels)
end

function tricontourf(m::TriMesh, z, levels)
    @assert issorted(levels)
    nlevels = length(levels)
    lmin = minimum(levels)
    lmax = (nlevels > 1) ? maximum(levels[1:end-1]) : lmin
    GR.setfillintstyle(1)
    for i=1:nlevels-1
        lower = levels[i]
        upper = levels[i+1]
        contours = generate_filled_contours(m, z, lower, upper)
        x = Float64[]
        y = Float64[]
        for c=contours
            append_with_nan!(x, first.(c))
            append_with_nan!(y, last.(c))
        end
        GR.setfillcolorind(getcolorind(levels[i],lmin,lmax))
        GR.fillarea(x,y)
    end
end

function append_with_nan!(a, b)
    append!(a, b)
    push!(a, NaN)
end

function generate_filled_contours(m::TriMesh, z, lower, upper)
    reset_visited!(m)
    contours = Vector{Contour}()
    add_filled_bdr_lines!(contours, m, z, lower, upper)
    add_interior_lines!(contours, m, z, lower; on_upper=false, filled=true);
    add_interior_lines!(contours, m, z, upper; on_upper=true,  filled=true);
    contours
end

function follow_boundary!(contour::Contour, m::TriMesh, edge::TriEdge, z, lower, upper, on_upper)
    ib,ie = m.bdr_map[edge]
    m.bdr_used[ib] = true

    stop = false
    first_edge = true
    z_end = 0.0
    while !stop
        @assert !m.bdr_visited[ib][ie]
        m.bdr_visited[ib][ie] = true
        if first_edge
            z_start = z[m.t[edge.ie,edge.it]]
        else
            z_start = z_end
        end
        z_end = z[m.t[mod1(edge.ie+1,3),edge.it]]

        if z_end > z_start
            if (on_upper || !first_edge) && (z_end >= lower) && (z_start < lower)
                stop = true
                on_upper = false
            elseif (z_end >= upper) && (z_start < upper)
                stop = true
                on_upper = true
            end
        else
            if (!on_upper || !first_edge) && (z_start >= upper) && (z_end < upper)
                stop = true
                on_upper = true
            elseif (z_start >= lower) && (z_end < lower)
                stop = true
                on_upper = false
            end
        end
        first_edge = false
        if !stop
            ie = mod1(ie+1, length(m.bdr_loops[ib]))
            edge = m.bdr_loops[ib][ie]
            point = m.t[edge.ie,edge.it]
            push!(contour, (m.x[point], m.y[point]))
        end
    end
    return on_upper,edge
end

function add_filled_bdr_lines!(contours, m::TriMesh, z, lower, upper)
    for i=1:length(m.bdr_loops)
        bdr = m.bdr_loops[i]
        for j=1:length(bdr)
            m.bdr_visited[i][j] && continue
            edge = bdr[j]

            z_start = z[m.t[edge.ie,edge.it]]
            z_end = z[m.t[mod1(edge.ie+1,3),edge.it]]

            # Does this boundary edge's z increase through upper level and/or decrease
            # through lower level?
            incr_upper = (z_start < upper && z_end >= upper)
            decr_lower = (z_start >= lower && z_end < lower)

            if decr_lower || incr_upper
                start_edge = edge
                contour = Contour()
                on_upper = incr_upper
                while true
                    level_int = on_upper ? upper : lower
                    edge = follow_interior!(contour, m, edge, z, level_int; end_on_bdr=true, on_upper)
                    on_upper,edge = follow_boundary!(contour, m, edge, z, lower, upper, on_upper)
                    edge == start_edge && break
                end
                # Ensure closed loop
                push!(contour, first(contour))
                push!(contours, contour)
            end
        end
    end

    # Add full boundaries that lie between the lower and upper levels.  These are boundaries
    # that have not been touched by an internal contour line.
    for i=1:length(m.bdr_loops)
        m.bdr_used[i] && continue
        bdr = m.bdr_loops[i]
        zval = z[m.t[bdr[1].ie,bdr[1].it]]
        if (zval >= lower) && (zval < upper)
            contour = Contour()
            for j=1:length(bdr)
                edge = bdr[j]
                point = m.t[edge.ie,edge.it]
                push!(contour, (m.x[point], m.y[point]))
            end
            push!(contour, first(contour))
            push!(contours,contour)
        end
    end
end


text:
libap = dlopen("libargyris_pack.so")

__ap_ref_values         = dlsym(libap, :ap_ref_values)
__ap_ref_gradients      = dlsym(libap, :ap_ref_gradients)
__ap_ref_hessians       = dlsym(libap, :ap_ref_hessians)
__ap_physical_maps      = dlsym(libap, :ap_physical_maps)
__ap_physical_values    = dlsym(libap, :ap_physical_values)
__ap_physical_gradients = dlsym(libap, :ap_physical_gradients)
__ap_physical_hessians  = dlsym(libap, :ap_physical_hessians)
__ap_matrix_mass        = dlsym(libap, :ap_matrix_mass)
__ap_matrix_stiffness   = dlsym(libap, :ap_matrix_stiffness)
__ap_matrix_biharmonic  = dlsym(libap, :ap_matrix_biharmonic)

# ------------------------------------------------------------------------------
# Julia interfaces to the .so file.
# ------------------------------------------------------------------------------
function ap_ref_values{T}(x::Vector{T}, y::Vector{T})
# evaluate the 21 Argyris basis functions at quadrature points specified by
# 'points'.
    check_size(x,y)
    function_values = zeros(21,size(x)[1])
    ccall(__ap_ref_values, Void,
          (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{Float64}),
          x, y, length(x), function_values)
    return function_values
end

function ap_ref_gradients{T}(x::Vector{T}, y::Vector{T})
# evaluate the derivatives of the 21 Argyris basis functions at quadrature
# quadrature points specified by 'points'.
    check_size(x,y)
    dx = zeros(21,size(x)[1])
    dy = zeros(21,size(x)[1])
    ccall(__ap_ref_gradients, Void,
          (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{Float64}, Ptr{Float64}),
          x, y, size(x)[1], dx, dy)
    return dx, dy
end

function ap_ref_hessians{T}(x::Vector{T}, y::Vector{T})
# evaluate the second derivatives of the 21 Argyris basis functions at quadrature
# quadrature points specified by 'points'.
    check_size(x,y)

    dxx = zeros(21,size(x)[1])
    dxy = zeros(21,size(x)[1])
    dyy = zeros(21,size(x)[1])
    ccall(__ap_ref_hessians, Void,
          (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}),
          x, y, length(x), dxx, dxy, dyy)
    return dxx, dxy, dyy
end

function ap_physical_maps{T}(x::Vector{T}, y::Vector{T})
# Use the change-of-basis matrix C to convert reference values to physical values.
    check_size(x,y)
    if length(x) != 3
        error("There should be three corners to a triangle.")
    end
    C = zeros(21,21)
    B = zeros(2,2)
    b = zeros(2)
    ccall(__ap_physical_maps, Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
          x, y, C, B, b)
    return C, B, b
end

function ap_physical_values(C, ref_values)
# Use the change-of-basis matrix C to convert reference values to physical values.
    check_transformations(C)
    check_ref_values(ref_values)

    function_values = zeros(21,size(ref_values)[2])
    ccall(__ap_physical_values, Void, (Ptr{Float64}, Ptr{Float64}, Int32,
          Ptr{Float64}), C, ref_values, size(ref_values)[2],
          function_values)
    return function_values
end

function ap_physical_gradients(C, B, ref_dx, ref_dy)
# evaluate the derivatives of the 21 Argyris basis functions given values of
# the derivatives on the reference triangle.
    check_transformations(C,B)
    check_ref_values(ref_dx, ref_dy)

    dx = zeros(size(ref_dx))
    dy = zeros(size(ref_dx))
    ccall(__ap_physical_gradients, Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32,
           Ptr{Float64}, Ptr{Float64}),
          C, B, ref_dx, ref_dy, size(ref_dx)[2], dx, dy)
    return dx, dy
end

function ap_physical_hessians(C, B, ref_dxx, ref_dxy, ref_dyy)
# Evaluate second derivatives of the basis functions given reference values.
    check_transformations(C,B)
    check_ref_values(ref_dxx, ref_dxy, ref_dyy)

    dxx = zeros(size(ref_dxx))
    dxy = zeros(size(ref_dxx))
    dyy = zeros(size(ref_dxx))
    ccall(__ap_physical_hessians, Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
          C, B, ref_dxx, ref_dxy, ref_dyy, size(ref_dxx)[2], dxx, dxy, dyy)
    return dxx, dxy, dyy
end

function ap_matrix_mass(C, B, ref_functions, weights)
# Evaluate the mass matrix given reference values.
    check_transformations(C,B)
    check_ref_values(ref_values)
    check_weights(ref_values, weights)
    mass = zeros(21,21)
    ccall(__ap_matrix_mass, Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32,
           Ptr{Float64}),
          C, B, ref_values, weights, size(weights)[1], mass)
    return mass
end

function ap_matrix_stiffness(C, B, ref_dx, ref_dy, weights)
# Evaluate the stiffness matrix given reference values.
    check_transformations(C,B)
    check_ref_values(ref_dx, ref_dy)
    check_weights(ref_dx, weights)

    stiffness = zeros(21,21)
    ccall(__ap_matrix_stiffness, Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Int32, Ptr{Float64}),
          C, B, ref_dx, ref_dy, weights, size(weights)[1], stiffness)
    return stiffness
end

function ap_matrix_biharmonic(C, B, ref_dxx, ref_dxy, ref_dyy, weights)
# Evaluate the biharmonic matrix given reference values.
    check_transformations(C,B)
    check_ref_values(ref_dxx, ref_dxy, ref_dyy)
    check_weights(ref_dxx, weights)

    biharmonic = zeros(21,21)
    ccall(__ap_matrix_biharmonic, Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
           Ptr{Float64}, Int32, Ptr{Float64}),
          C, B, ref_dxx, ref_dxy, ref_dyy, weights, size(weights)[1],
          biharmonic)
    return biharmonic
end

# ------------------------------------------------------------------------------
# These functions pertain to validating input.
# ------------------------------------------------------------------------------
function check_size(arrays...)
# raise an exception if the arrays are not equal in size.
    if length(arrays) == 1
        return
    else
        for i=2:length(arrays)
            if  size(arrays[i-1]) != size(arrays[i])
                error("Mismatch in array sizes.")
            end
        end
    end
end

function check_transformations(transformations...)
# check the sizes of the matrix transformations C and B.
    if size(transformations[1]) != (21,21)
        error("Incorrect dimensions in the 21x21 'C' matrix.")
    end
    if length(transformations) > 1
        if size(transformations[2]) != (2,2)
            error("Incorrect dimensions in the 2x2 'B' matrix.")
        end
    end
end

function check_ref_values(values...)
# Raise an exception if the input arrays do not all have 21 rows and the same
# number of columns.
    if size(values[1])[1] != 21
        error("There should be 21 rows corresponding to 21 basis functions.")
    end
    if length(values) == 1
        return
    else
        for i=2:length(values)
            if  size(values[i-1]) != size(values[i])
                error("Mismatch in array sizes.")
            end
        end
    end
end

function check_weights{T}(values, weights::Vector{T})
    if size(weights)[1] != size(values)[2]
        error("Mismatch in number of weights and number of quadrature points.")
    end
end

max_stars_repo_name: VT-ICAM/ArgyrisPack

max_issues_repo_name: VT-ICAM/ArgyrisPack

max_forks_repo_name: VT-ICAM/ArgyrisPack


text:
import RDates: build

build()


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    dcc_dropdown(id = "input",
        options = [
            (label="regular", value="regular"),
            (label="prevent", value="prevent"),
            (label="no_update", value="no_update"),
        ]
    ),
    html_div(id="output"),
    html_div(id="regular_output")
end
callback!(app, callid"input.value => regular_output.children") do input
    return input
end
callback!(app, callid"input.value => output.children") do input
    input == "prevent" && throw(PreventUpdate())
    input == "no_update" && return no_update()
    return input
end

run_server(app)


text:
using Test
using BenchmarkTools

using PersistentCollections

const ENV_DIR = "env.lmdb"
const UNSAFE_ENV_DIR = "unsafe_env.lmdb"
const TEST_DIRS = [ENV_DIR, UNSAFE_ENV_DIR]
# clean previous runs
for d in TEST_DIRS
    rm(d, force=true, recursive=true)
    mkdir(d)
end

@test string(LMDB.LMDBError(-30798)) == "MDB_NOTFOUND: No matching key/data pair found (-30798)"

env = LMDB.Environment(ENV_DIR)
@test isopen(env)
@test env == LMDB.Environment(ENV_DIR)
@test close(env)
@test !isopen(env)
@test !close(env) # nothing bad should happen

env = LMDB.Environment(ENV_DIR, maxdbs=1)

struct Immutable
    a::Int
    b::Float64
    c::Tuple{Int,Float64}
end
Base.:(==)(a::Immutable, b::Immutable) = a.a == b.a && a.b == b.b && a.c == b.c

mutable struct Mutable
    a::Int
    b::Float64
    c::Tuple{Int,Float64}
end
Base.:(==)(a::Mutable, b::Mutable) = a.a == b.a && a.b == b.b && a.c == b.c

ro = Immutable(1, 2.5, (2, 3.4))
@assert isimmutable(ro)

rw = Mutable(1, 2.5, (2, 3.4))
@assert !isimmutable(rw)

fastkey = LMDB.MDBValue("fastkey")
fastval = LMDB.MDBValue("fastval")

# Open up default database as a Dict
d = PersistentDict{Any,Any}(env, id="foo")

@test_throws BoundsError first(d)
@test_throws BoundsError last(d)

d["stringkey"] = "stringval"
d["byteskey"] = Vector{UInt8}("bytesval")
d["intkey"] = 1234
d["floatkey"] = 2.5
d["tuplekey"] = (1, 2.5)
d["immutable_struct_key"] = ro
d["mutable_struct_key"] = rw
d[fastkey] = fastval
# Set use-case
d["nothing_key"] = nothing

# sorted lexically
@test first(d) == Pair(LMDB.MDBValue("byteskey"), LMDB.MDBValue(Vector{UInt8}("bytesval")))
@test last(d) == Pair(LMDB.MDBValue("tuplekey"), LMDB.MDBValue((1, 2.5)))

@test get(d, "stringkey", "") == "stringval"
@test get(d, "byteskey", UInt8[]) == Vector{UInt8}("bytesval")
@test get(d, "intkey", 0) == 1234
@test get(d, "floatkey", 0.0) == 2.5
@test get(d, "tuplekey", (0, 0.0)) == (1, 2.5)
@test get(d, "immutable_struct_key", Immutable(0, 0.0, (0, 0.0))) == ro
@test get(d, "immutable_struct_key", Mutable(0, 0.0, (0, 0.0))) == rw
@test get(d, fastkey, "") == "fastval"
# Set use-case
notfound = UInt8[0x1]
@test get(d, "nothing_key", notfound) != notfound

# delete
@test delete!(d, fastkey)
@test delete!(d, "stringkey")
# idempotency
@test !delete!(d, fastkey)

# check deleted
@test_throws KeyError d[fastkey]
@test_throws KeyError d["stringkey"]

@test convert(Int, d["intkey"]) == 1234
@test length(collect(keys(d))) > 0
@test length(collect(values(d))) == length(collect(keys(d)))
@test !isempty(collect(d))
@test (d["dictkey"] = "dictval") == "dictval"

# manual sync
unsafe_env = LMDB.Environment(UNSAFE_ENV_DIR, maxdbs=1, flags=LMDB.MDB_WRITEMAP | LMDB.MDB_MAPASYNC)
unsafe_dict = PersistentDict{String,Vector{UInt8}}(unsafe_env, id="foo")
flush(unsafe_env) do
    unsafe_dict["unsafe_key1"] = Vector{UInt8}("unsafe_val1")
    unsafe_dict["unsafe_key2"] = Vector{UInt8}("unsafe_val2")
end
# flushed when the flush block exits
unsafe_dict["unsafe_key1"] == Vector{UInt8}("unsafe_val1")
unsafe_dict["unsafe_key2"] == Vector{UInt8}("unsafe_val1")

# multi-threading - don't crash or deadlock while re-writing key/val
if Threads.nthreads() > 1
    @info "doing $(Threads.nthreads()) parallel writes for ~10s (should not crash or deadlock) ..."
    randvals = [rand(Int) for _ in 1:Threads.nthreads()]
    deadline = time() + 10
    writes = [zero(Int) for _ in 1:Threads.nthreads()]
    Threads.@threads for i in 1:Threads.nthreads()
        while time() < deadline
            d[fastkey] = randvals[i]
            writes[i] += 1
        end
    end
    @info "... $(sum(writes)) writes completed!"
else
    @warn "skipping multi-threading tests (re-run with JULIA_NUM_THREADS=999999 to enable)"
end

# == Benchmarks ==

if get(ENV, "BENCH", "") == "y"
    longkey = "012345678901234567890123456789012345678901234567890123456789"

    @info "Benchmarking setindex!(::PersistentDict) ..."
    @btime setindex!(d, v, longkey) setup=(v=rand(UInt8, 500))

    @info "Benchmarking getindex(::PersistentDict) ..."
    @btime getindex(d, longkey)

    n = length(collect(keys(d)))
    @info "Benchmarking itertion: keys(::PersistentDict)) ($n entries) ..."
    @btime for _ in keys(d) end

    @info "Benchmarking itertion: iterated(::PersistentDict)) ($n entries) ..."
    @btime for _ in d end

    @info "Benchmarking last(::PersistentDict) ..."
    @btime last(d)

    @info "Benchmarking (MDB_WRITEMAP | MDB_MAPASYNC) with setindex!(::PersistentDict) ..."
    @btime setindex!(unsafe_dict, v, longkey) setup=(v=rand(UInt8, 500))
end


text:
# -----------------------------------------------------------------------------
# Implementation of Atomistic MolecularDynamicsResult API
# -----------------------------------------------------------------------------

"""
    NBSResult <: MolecularDynamicsResult

The result generated from running an `NBSimulator`.

**Field descriptions**
- `result::SimulationResult` the standard simulation result from `NBodySimulator`
- `energy_cache::Vector{Float64}` cached potential energy values from the simulation
"""
struct NBSResult <: MolecularDynamicsResult
    result::SimulationResult
    energy_cache::Vector{ENERGY_TYPE}
end

get_time_range(result::NBSResult) = result.result.solution.t * TIME_UNIT
get_num_bodies(result::NBSResult) = length(result.result.simulation.system.bodies)
get_bounding_box(result::NBSResult) = get_bounding_box(result.result.simulation.boundary_conditions)
get_boundary_conditions(result::NBSResult) = get_boundary_conditions(result.result.simulation.boundary_conditions)

reference_temperature(result::NBSResult) = reference_temperature(result.result.simulation.thermostat)
reference_temperature(thermostat::Thermostat) = thermostat.T * TEMPERATURE_UNIT
reference_temperature(::NullThermostat) = missing

function get_positions(result::NBSResult, t::Integer)
    positions = get_position(result.result, result.result.solution.t[t])
    [bound_position(SVector{3}(p), result.result.simulation.boundary_conditions) * LENGTH_UNIT for p ∈ eachcol(positions)]
end
function get_velocities(result::NBSResult, t::Integer)
    velocities = get_velocity(result.result, result.result.solution.t[t])
    [SVector{3}(v) * VELOCITY_UNIT for v ∈ eachcol(velocities)]
end
get_particles(result::NBSResult, t::Integer) = AtomsBase.Atom.(result.result.simulation.system.bodies, get_positions(result, t), get_velocities(result, t))

temperature(result::NBSResult, t::Integer) = NBodySimulator.temperature(result.result, result.result.solution.t[t]) * TEMPERATURE_UNIT
kinetic_energy(result::NBSResult, t::Integer) = NBodySimulator.kinetic_energy(result.result, result.result.solution.t[t]) * ENERGY_UNIT
potential_energy(result::NBSResult, t::Integer) = result.energy_cache[t]

"""
    animate(result::NBSResult, filename::String; kwargs...)

Animate an `NBSResult` using the interface provided by NBodySimulator.jl and store the result in a file.
"""
Plots.animate(result::NBSResult, filename::String; kwargs...) = animate(result.result, filename; kwargs...)


text:
@testset "Solver" begin
    srand(42)
    N = 100
    x = sort(10 .* rand(N))
    y = sin.(x)
    yerr = 0.01 .+ rand(N) ./ 100

    kernel = celerite.RealTerm(0.5, 1.0) + celerite.SHOTerm(0.1, 2.0, -0.5)
    gp = celerite.Celerite(kernel)

    K = celerite.get_matrix(gp, x)
    for n in 1:N
        K[n, n] = K[n, n] + yerr[n]^2
    end

    # Compute using celerite
#    celerite.compute!(gp, x, yerr)
    celerite.compute_ldlt!(gp, x, yerr)
#    ll = celerite.log_likelihood(gp, y)
    ll = celerite.log_likelihood_ldlt(gp, y)

    # Compute directly
    ll0 = -0.5*(sum(y .* (K \ y)) + logdet(K) + N*log(2.0*pi))

    @test isapprox(ll, ll0)
end


text:
export mfactor

const smsievelib =  "libsmsieve.so"

type Msieveopts
    n::AbstractString
    deadline::Int
    logfile::AbstractString
    deepecm::Bool
    info::Bool
end

# Send the string to msieve and return c struct msieve_obj
function runmsieve(opts::Msieveopts)
    n = opts.n  # n is a string
    d = opts.deadline
    logfile = opts.logfile
    deepecm = opts.deepecm
    numcores = 1
    ecmflag = deepecm ? 1 : 0
    infoflag = opts.info ? 1 : 0
    res = try
        if (logfile == "")
            ccall((:factor_from_string,smsievelib), Ptr{Void}, (Ptr{UInt8},Int,Int,Ptr{UInt8},Int,Int),
                  n, numcores, d, C_NULL, ecmflag, infoflag)
        else
            ccall((:factor_from_string,smsievelib), Ptr{Void}, (Ptr{UInt8},Int,Int,Ptr{UInt8},Int, Int),
                  n, numcores, d, logfile, ecmflag, infoflag)
        end
    catch
        error("factor_from_string failed")
    end
    res == C_NULL && throw(InterruptException())
    res
end

# Send ptr to msieve_obj and get ptr to struct factors.
getfactors(obj) = ccall((:get_factors_from_obj,smsievelib), Ptr{Void}, (Ptr{Void},), obj)
# Sent ptr to struct factors and get number of factors.
get_num_factors(factors) = ccall((:get_num_factors,smsievelib), Int, (Ptr{Void},), factors)
msieve_free(obj) =  ccall((:msieve_obj_free_2,smsievelib), Void, (Ptr{Void},), obj)

# Send ptr to struct factor and get string rep of one factor.
# A ptr to next struct factor, correponding to next factor, is returned.
function get_one_factor_value(factor)
    factorstring = Array{UInt8}(500) # max num digits input to msieve is 300
    nextfactor = ccall((:get_one_factor_value,smsievelib), Ptr{Void}, (Ptr{Void},Ptr{UInt8},Int),
                       factor,factorstring,length(factorstring))
    factorstring[end] = 0
    return(nextfactor,unsafe_string(pointer(factorstring)))
end

# Send ptr to first struct factor. Return all factors as array of strings
function get_all_factor_values(factor)
    allf = Array{AbstractString}(0)
    nfactor = factor
    n = get_num_factors(factor)
    for i in 1:n
        (nfactor,sfact) = get_one_factor_value(nfactor)
        push!(allf,sfact)
    end
    return allf
end

# Send n as string to msieve, return all factors as array of strings
function runallmsieve(opts::Msieveopts)
    obj = runmsieve(opts)
    sfactors = get_all_factor_values(getfactors(obj))
    msieve_free(obj)
    sfactors
end

# input factors as Array of strings. Output Array of Integers (Usually of type Int)
function factor_strings_to_integers(sfactors::Array{AbstractString})
    m = length(sfactors)
    n1 = eval(parse(BigInt,sfactors[m]))  ## why eval ?
    T = typeof(n1)
    arr = Array{T}(m)
    arr[m] = n1
    for i in 1:m-1
        arr[i] = parse(T,sfactors[i])
    end
    arr
end

# Send string to msieve. Return factors as list of Integers.
mfactorl(opts::Msieveopts) = factor_strings_to_integers(runallmsieve(opts))

# Send string to msieve. Return factors in Dict, like Base.factor
function mfactor(opts::Msieveopts)
    arr = mfactorl(opts)
    T = eltype(arr)
    d = Dict{T,Int}()
    for i in arr d[i] = get(d,i,0) + 1 end
    d
end

function mfactor(x::Union{AbstractString,Integer}; deadline::Integer = 0, logfile::AbstractString = "", ecm::Bool = false,
                 info::Bool = false)
    mfactor(Msieveopts(string(x),deadline,logfile,ecm,info))
end

for (thetype) in ( :AbstractString, :Integer )
    @eval begin
        function mfactor{T<:$thetype}(a::AbstractArray{T,1}; dl::Integer=0, logfile::AbstractString = "",
                                      ecm::Bool = false, info::Bool = false)
            outa = Array{Any}(0)
            for x in a
                res = mfactor(x; deadline=dl, logfile=logfile, ecm=ecm, info=info)
                push!(outa,res)
            end
            outa
        end
    end
end