text:
# ---
# title: 1499. Max Value of Equation
# id: problem1499
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Sliding Window
# link: <https://leetcode.com/problems/max-value-of-equation/description/>
# hidden: true
# ---
#
# Given an array `points` containing the coordinates of points on a 2D plane,
# sorted by the x-values, where `points[i] = [xi, yi]` such that `xi < xj` for
# all `1 <= i < j <= points.length`. You are also given an integer `k`.
#
# Find the _maximum value of the equation_`yi + yj + |xi - xj|` where `|xi - xj|
# <= k` and `1 <= i < j <= points.length`. It is guaranteed that there exists at
# least one pair of points that satisfy the constraint `|xi - xj| <= k`.
#
#
#
# **Example 1:**
#
#
#
#     Input: points = [[1,3],[2,0],[5,10],[6,-10]], k = 1
#     Output: 4
#     Explanation: The first two points satisfy the condition |xi - xj| <= 1 and if we calculate the equation we get 3 + 0 + |1 - 2| = 4. Third and fourth points also satisfy the condition and give a value of 10 + -10 + |5 - 6| = 1.
#     No other pairs satisfy the condition, so we return the max of 4 and 1.
#
# **Example 2:**
#
#
#
#     Input: points = [[0,0],[3,0],[9,2]], k = 3
#     Output: 3
#     Explanation: Only the first two points have an absolute difference of 3 or less in the x-values, and give the value of 0 + 0 + |0 - 3| = 3.
#
#
#
#
# **Constraints:**
#
#   * `2 <= points.length <= 10^5`
#   * `points[i].length == 2`
#   * `-10^8 <= points[i][0], points[i][1] <= 10^8`
#   * `0 <= k <= 2 * 10^8`
#   * `points[i][0] < points[j][0]` for all `1 <= i < j <= points.length`
#   * `xi` form a strictly increasing sequence.
#
#
## @lc code=start
using LeetCode

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


text:
# Inverse CDF method
#
# copied from Distributions.jl, Feb 14, 2015
function randtn_invcdf(d::Truncated)
    d0 = d.untruncated
    if d.tp > 0.25
        while true
            r = rand(d0)
            if d.lower <= r <= d.upper
                return r
            end
        end
    else
        return quantile(d0, d.lcdf + rand() * d.tp)
    end
end

# Rejection sampler based on algorithm from Robert (1995)
#
#  - Available at http://arxiv.org/abs/0907.4010
#
# copied from Distributions.jl, Feb 14, 2015
function randtn_robert(d::Truncated{Normal})
    d0 = d.untruncated
    μ = mean(d0)
    σ = std(d0)
    a = (d.lower - μ) / σ
    b = (d.upper - μ) / σ
    z = randtn_robert_base(a, b, d.tp)
    return μ + σ * z
end

function randtn_robert_base(lb::Float64, ub::Float64, tp::Float64)
    r::Float64
    if tp > 0.3   # has considerable chance of falling in [lb, ub]
        r = randn()
        while r < lb || r > ub
            r = randn()
        end
        return r

    else
        span = ub - lb
        if lb > 0 && span > 2.0 / (lb + sqrt(lb^2 + 4.0)) * exp((lb^2 - lb * sqrt(lb^2 + 4.0)) / 4.0)
            a = (lb + sqrt(lb^2 + 4.0))/2.0
            while true
                r = rand(Exponential(1.0 / a)) + lb
                u = rand()
                if u < exp(-0.5 * (r - a)^2) && r < ub
                    return r
                end
            end
        elseif ub < 0 && ub - lb > 2.0 / (-ub + sqrt(ub^2 + 4.0)) * exp((ub^2 + ub * sqrt(ub^2 + 4.0)) / 4.0)
            a = (-ub + sqrt(ub^2 + 4.0)) / 2.0
            while true
                r = rand(Exponential(1.0 / a)) - ub
                u = rand()
                if u < exp(-0.5 * (r - a)^2) && r < -lb
                    return -r
                end
            end
        else
            while true
                r = lb + rand() * (ub - lb)
                u = rand()
                if lb > 0
                    rho = exp((lb^2 - r^2) * 0.5)
                elseif ub < 0
                    rho = exp((ub^2 - r^2) * 0.5)
                else
                    rho = exp(-r^2 * 0.5)
                end
                if u < rho
                    return r
                end
            end
        end
    end
end

ext: jl

max_stars_repo_name: jdrugo/fasttn

max_issues_repo_name: jdrugo/fasttn

max_forks_repo_name: jdrugo/fasttn


text:
export Posterior, asflat, ascube, flatten, logdensity, transform, inverse, dimension

import DensityInterface
import ParameterHandling
import ParameterHandling: flatten
using HypercubeTransform
using TransformVariables
using ValueShapes: NamedTupleDist

"""
    Posterior(lklhd, prior, model)
Posterior density that follows obeys the [DensityInferface](https://github.com/JuliaMath/DensityInferface.jl)

The expected arguments are:

- lklhd: Which should be an intance of RadioLikelihood with whatever data products you want to fit
- prior: This should be a `NamedTuple` with the priors for each model ParameterHandling
- model: Function that takes a `NamedTuple` of parameters and constructs the `Comrade` model.

To evaluate the logdensity of the posterior use can either use `DensityInterface.logdensityof`
or `logdensity`.
"""
struct Posterior{L,P,F}
    lklhd::L
    prior::P
    model::F
end

function Posterior(lklhd, prior::NamedTuple, model)
    return Posterior(lklhd, NamedTupleDist(prior), model)
end

@inline DensityInterface.DensityKind(::Posterior) = IsDensity()

function DensityInterface.logdensityof(post::Posterior, x)
    pr = logdensity(post.prior, x)
    !isfinite(pr) && return -Inf
    return logdensity(post.lklhd, post.model(x)) + pr
end

"""
    `$(TYPEDEF)`
A transformed version of the posterior `lpost`. This is an internal type
that an end user shouldn't have to directly construct.

To construct a transformed posterior see the `asflat`, `ascube`, and `flatten`
functions.

When evaluating the `logdensity` of the `TransformedPosterior`, `logdenisty`
will include any jacobian terms automatically, and expects the argument to be the
vector of transformed quantities. This is usually easier to incorporate with
samplers, and optimizers, which expect homogenous vectors as arguments.
"""
struct TransformedPosterior{P<:Posterior,T}
    lpost::P
    transform::T
end

"""
    transform(posterior::TransformedPosterior, x)
Transforms the value `x` into parameter space, i.e. usually a NamedTuple.
"""
HypercubeTransform.transform(p::TransformedPosterior, x) = transform(p.transform, x)

"""
    inverse(posterior::TransformedPosterior, x)
Transforms the value model parameters `x` into the flattened transformed space.
"""
HypercubeTransform.inverse(p::TransformedPosterior, y) = inverse(p.transform, y)

@inline DensityInterface.DensityKind(::TransformedPosterior) = IsDensity()

MeasureBase.logdensity(tpost::Union{Posterior,TransformedPosterior}, x) = DensityInterface.logdensityof(tpost, x)

"""
    asflat(post::Posterior)
Construct a flattened version of the posterior, where the parameters are transformed so that
their support is from (-∞, ∞). This uses [TransformVariables](https://github.com/tpapp/TransformVariables.jl)

The transformed posterior can then be evaluated by the `logdensityof(transformed_posterior,x)`
method following the `DensityInterface`, where `x` is a flattened vector of the infinite support
variables. **Note** this already includes the jacobian of the transformation so this does
not need to be added.

This is useful for optimization and sampling algorithms such as HMC that will use gradients
to explore the posterior surface.
"""
function HypercubeTransform.asflat(post::Posterior)
    pr = getfield(post.prior, :_internal_distributions)
    tr = asflat(pr)
    return TransformedPosterior(post, tr)
end

function DensityInterface.logdensityof(post::TransformedPosterior{P, T}, x) where {P, T<:TransformVariables.AbstractTransform}
    p, logjac = transform_and_logjac(post.transform, x)
    return DensityInterface.logdensityof(post.lpost, p) + logjac
end

HypercubeTransform.dimension(post::TransformedPosterior) = dimension(post.transform)
HypercubeTransform.dimension(post::Posterior) = length(first(ParameterHandling.flatten(rand(post.prior))))

"""
    ascube(post::Posterior)
Construct a flattened version of the posterior, where the parameters are transformed to live
in the unit hypercube. In astronomy parlance, we are transforming the variables to the unit
hypercube. This is done using the [HypercubeTransform](https://github.com/ptiede/HypercubeTransform.jl)
package.

The transformed posterior can then be evaluated by the `logdensityof(transformed_posterior,x)`
method following the `DensityInterface`, where `x` vector that lives in the unit hypercube.
**Note** this already includes the jacobian of the transformation so this does
not need to be added.

This transform is useful for NestedSampling methods that often assume that the model is written
to live in the unit hypercube.
"""
function HypercubeTransform.ascube(post::Posterior)
    pr = getfield(post.prior, :_internal_distributions)
    tr = ascube(pr)
    return TransformedPosterior(post, tr)
end

function DensityInterface.logdensityof(tpost::TransformedPosterior{P, T}, x) where {P, T<:HypercubeTransform.AbstractHypercubeTransform}
    # Check that x really is in the unit hypercube. If not return -Inf
    for xx in x
        (xx > 1 || xx < 0) && return -Inf
    end
    p = transform(tpost.transform, x)
    post = tpost.lpost
    return logdensity(post.lklhd, post.model(p))
end

struct FlatTransform{T}
    transform::T
end

HypercubeTransform.transform(t::FlatTransform, x) = t.transform(x)
HypercubeTransform.inverse(::FlatTransform, x) = first(ParameterHandling.flatten(x))
HypercubeTransform.dimension(t::FlatTransform) = t.transform.unflatten.sz[end]
"""
    flatten(post::Posterior)
Flatten the representation of the posterior. Internally this uses ParameterHandling to
construct a flattened version of the posterior.

Note this is distinct from `asflat` that transforms the variables to live in (-∞,∞).
Instead this method just flattens the repsentation of the model from a NamedTuple to a vector.
This allows the easier integration to optimization and sampling algorithms.
"""
function ParameterHandling.flatten(post::Posterior)
    x0 = rand(post.prior)
    _, unflatten = ParameterHandling.flatten(x0)
    return TransformedPosterior(post, FlatTransform(unflatten))
end

function DensityInterface.logdensityof(post::TransformedPosterior{P,T}, x) where {P, T<: FlatTransform}
    return logdensity(post.lpost, transform(post.transform, x))
end


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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a565c4ee-d425-11eb-0a5c-d7d322fb95c6
using IndexFunArrays, Plots, Colors, ImageShow, PlutoUI, MappedArrays

# ╔═╡ eb26d519-61f9-49c8-bf13-fb0bd5b01be6
md"#### IndexFunArrays.jl | Fun with indices (and functions on them)
* _Rainer Heintzmann, Leibniz Institute of Photonic Technology in Jena, Germany_
* _Felix Wechsler, Leibniz Institute of Photonic Technology in Jena, Germany_

Main features:
* Arrays without memory allocations
* entries are calculated once accessed
* feels like an arrays but is some kind of generator

"

# ╔═╡ 840c51ce-1300-41b0-be75-4b298c6c3390
md"### Motivation

* Calculation with indices works efficiently with `CartesianIndices`
* However, as you can see it can get complicated
"

# ╔═╡ 1d880f4b-71fb-4b83-8d23-4a46a1481423
function rr_cartesian_indices(s)
	# entry wise function
	rr_fun(index) = √(sum(abs2.(Tuple(index) .- (s .÷ 2 .+ 1))));
	# array of tuple indices
	arr_indices = CartesianIndices(s)
	# apply rr_fun elementwise
	rr_fun.(arr_indices)
end

# ╔═╡ ab66c67d-832b-40ff-bbc2-f52761912b2d
rr_cartesian_indices((5,6))

# ╔═╡ f09a33fa-dd83-4dbd-948b-78303d3fb76d
md"### MappedArrays
* in principle offers similar functionality
* but doesn't provide the explicit functions (like `rr2`, `rr`)
"

# ╔═╡ e70c0bdb-aa02-4017-ae32-ba3b9fb50906
function rr_mapped_arrays(s)
	# entry wise function
	rr_fun(index) = √(sum(abs2.(Tuple(index) .- (s .÷ 2 .+ 1))));
	# array of tuple indices
	mappedarray(rr_fun, CartesianIndices(s))
end

# ╔═╡ bee02ceb-24a4-4d48-bfe0-c24fdfde0472
rr_mapped_arrays((5,6))

# ╔═╡ 01651a42-dce0-4015-a9d0-a92dc02cbd89
rr((5, 6))

# ╔═╡ 3faba71d-ffa3-4bad-96b8-76f488dfc600
md"### Data Structure Behind It"

# ╔═╡ da4c455d-f657-4fc8-b8d5-ac1497e29071
struct IndexFunArray_demo{T, N, F} <: AbstractArray{T, N} where {F}
    # stores the generator function to be applied to the indices.
    generator::F
    # output size of the array
    size::NTuple{N, Int}
end

# ╔═╡ 5b7b9e0b-9f93-43ed-aba5-56cbf9f1bea8
Base.getindex(A::IndexFunArray{T,N}, I::Vararg{Int, N}) where {T,N} =
    return A.generator(I)

# ╔═╡ 8604d201-c4ba-4314-bc0b-df73f8d3e24a
md"### How To Use It"

# ╔═╡ 44af166f-a01f-43c2-a774-9a09f5d42712
hello = IndexFunArray(x -> "hello", (3, 3))

# ╔═╡ 3a2ddf05-0d2a-401e-bbad-adccf7ab21bf
hello[2, 3]

# ╔═╡ 5d5835da-4331-4f06-a24c-e60177f38cba
img = IndexFunArray(x -> x[1] + x[2], (3, 3))

# ╔═╡ df51963f-ceea-4279-aa51-6d8c8afaa9a0
distance_to_corner = IndexFunArray(x -> sqrt((x[1]-1)^2 + (x[2]-1)^2), (3, 3))

# ╔═╡ dc58b7a5-8c64-4c3f-aa7b-8261647d2513
md"### What is already defined?

#### Window Functions
* Gaussian
* Hamming
* Hanning
* Linear
* ...

#### Distance Functions
* radius
* radius squared
* 1D ramps
* ...
"

# ╔═╡ 455d650b-cb51-4dc3-a02f-e6c832e3e5e0

# ╔═╡ d48829c1-c04f-47d3-aa0d-af656b6a5276

# ╔═╡ fbfbc63b-4188-4c03-9bc3-b89b4c6c6298

# ╔═╡ c0eb6cff-a47c-485e-8983-2bf4070d5223
md"## Defining apertures in 1 LOC"

# ╔═╡ d4e01cd2-3eff-404f-86c6-fcd4b772f478
@bind radius Slider(0:1:200, show_value=true)

# ╔═╡ a9ab1f95-42d2-4d64-b015-423ec4fcfc8a
@bind c_ind Select(["1" => "CtrFT", "2" => "CtrEnd", "3" => "CtrFFT"])

# ╔═╡ 04b7c231-0868-493e-99ef-3eb4d89b7ad4
ctr = [CtrFT, CtrEnd, CtrFFT][parse(Int, c_ind)];

# ╔═╡ 4512c50f-e094-4ef9-a177-27a9fdc58f09
[Gray.(rr((300, 300), offset=ctr) .<= radius) Gray.(disc((300, 300), radius, offset=ctr))]

# ╔═╡ d7878c97-068d-4040-9f73-2400e05636bf

# ╔═╡ 2d23b2d8-5dac-422a-804d-e2ccd27de12e
Gray.(0.5 .+ xx((300, 300), offset=ctr, scale=ScaNorm))

# ╔═╡ 465522c9-d209-4476-867f-f3ba1f1563ca
md"## Different kinds of Window functions"

# ╔═╡ fabeb838-fb8e-46b9-ad97-9271d4ada93b
md"
$(@bind border_out Slider(0:0.01:1, show_value=true))

$(@bind border_in Slider(0:0.01:1, show_value=true))

"

# ╔═╡ 54023a37-4843-48a2-bb94-262590f310c3
begin
	w_gauss = window_gaussian((256, 256); border_in, border_out)
	w_linear = window_linear((256, 256); border_in, border_out)
end;

# ╔═╡ fedff1f5-5595-4ac9-b53f-a11bd864554c
[Gray.(w_gauss) Gray.(w_linear) ]

# ╔═╡ f0c86a16-f255-4b85-b32a-e0adba958504
md"# Gauss"

# ╔═╡ 5c691197-0516-4ab6-9804-0bb6d0df5a66
md"
$(@bind σ1 Slider(0:0.01:200, show_value=true))
$(@bind σ2 Slider(0:0.01:200, show_value=true))
"

# ╔═╡ 7f72a169-5c33-49f4-982e-e85181425374
w_gauss2 = gaussian((256, 256), sigma=(σ1, σ2));

# ╔═╡ 0fa4b07a-5fe2-4294-8753-9adb1ed7c6f7
Gray.(w_gauss2)

# ╔═╡ Cell order:
# ╠═a565c4ee-d425-11eb-0a5c-d7d322fb95c6
# ╟─eb26d519-61f9-49c8-bf13-fb0bd5b01be6
# ╟─840c51ce-1300-41b0-be75-4b298c6c3390
# ╠═1d880f4b-71fb-4b83-8d23-4a46a1481423
# ╠═ab66c67d-832b-40ff-bbc2-f52761912b2d
# ╟─f09a33fa-dd83-4dbd-948b-78303d3fb76d
# ╠═e70c0bdb-aa02-4017-ae32-ba3b9fb50906
# ╠═bee02ceb-24a4-4d48-bfe0-c24fdfde0472
# ╠═01651a42-dce0-4015-a9d0-a92dc02cbd89
# ╟─3faba71d-ffa3-4bad-96b8-76f488dfc600
# ╠═da4c455d-f657-4fc8-b8d5-ac1497e29071
# ╠═5b7b9e0b-9f93-43ed-aba5-56cbf9f1bea8
# ╟─8604d201-c4ba-4314-bc0b-df73f8d3e24a
# ╠═44af166f-a01f-43c2-a774-9a09f5d42712
# ╠═3a2ddf05-0d2a-401e-bbad-adccf7ab21bf
# ╠═5d5835da-4331-4f06-a24c-e60177f38cba
# ╠═df51963f-ceea-4279-aa51-6d8c8afaa9a0
# ╟─dc58b7a5-8c64-4c3f-aa7b-8261647d2513
# ╟─455d650b-cb51-4dc3-a02f-e6c832e3e5e0
# ╟─d48829c1-c04f-47d3-aa0d-af656b6a5276
# ╟─fbfbc63b-4188-4c03-9bc3-b89b4c6c6298
# ╟─c0eb6cff-a47c-485e-8983-2bf4070d5223
# ╟─d4e01cd2-3eff-404f-86c6-fcd4b772f478
# ╟─a9ab1f95-42d2-4d64-b015-423ec4fcfc8a
# ╟─04b7c231-0868-493e-99ef-3eb4d89b7ad4
# ╠═4512c50f-e094-4ef9-a177-27a9fdc58f09
# ╠═d7878c97-068d-4040-9f73-2400e05636bf
# ╠═2d23b2d8-5dac-422a-804d-e2ccd27de12e
# ╟─465522c9-d209-4476-867f-f3ba1f1563ca
# ╟─fabeb838-fb8e-46b9-ad97-9271d4ada93b
# ╠═54023a37-4843-48a2-bb94-262590f310c3
# ╠═fedff1f5-5595-4ac9-b53f-a11bd864554c
# ╟─f0c86a16-f255-4b85-b32a-e0adba958504
# ╟─5c691197-0516-4ab6-9804-0bb6d0df5a66
# ╠═7f72a169-5c33-49f4-982e-e85181425374
# ╠═0fa4b07a-5fe2-4294-8753-9adb1ed7c6f7


text:
x = randn(10)
y = randn(10)
r = x .- y

@testset "Robust Loss" begin
    δ = 0.5
    rlδ = RobustLoss(Huber(δ))
    @test rlδ isa RobustLoss{HuberRho{δ}}
    @test rlδ(r) == rlδ(x, y) == sum(ifelse(abs(rᵢ)≤δ, rᵢ^2/2, δ*(abs(rᵢ)-δ/2)) for rᵢ in r)

    rlδ = RobustLoss(Andrews(δ))
    p1(z) = -cos(π * z/δ)/(π/δ)^2
    p2(z) = p1(δ)
    @test rlδ(r) ≈ sum(ifelse(abs(rᵢ)≤δ, p1.(rᵢ), p2.(rᵢ)) for rᵢ in r)

    rlδ = RobustLoss(Bisquare(δ))
    p1(z) = δ^2/6 * (1-(1-(z/δ)^2)^3)
    p2(z) = p1(δ)
    @test rlδ(r) ≈ sum(ifelse(abs(rᵢ)≤δ, p1.(rᵢ), p2.(rᵢ)) for rᵢ in r)

    rlδ = RobustLoss(Logistic(δ))
    @test rlδ(r) ≈ sum(δ^2 * log(cosh(rᵢ/δ)) for rᵢ in r)

    rlδ = RobustLoss(FairRho(δ))
    @test rlδ(r) ≈ sum(δ^2 * (abs(rᵢ)/δ - log(1+abs(rᵢ/δ))) for rᵢ in r)

    rlδ = RobustLoss(Talwar(δ))
    p1(z) = z^2/2
    p2(z) = p1(δ)
    @test rlδ(r) ≈ sum(ifelse(abs(rᵢ)≤δ, p1.(rᵢ), p2.(rᵢ)) for rᵢ in r)
end


text:
###############################################################################
#
# Linear regression of Y~X, test slope
#
###############################################################################
"""
    test_slope(xs::Array,
               ys::Array;
               slope::Number = 0,
               intercept::Bool = true
    )

Make linear regression and return the p value of whether the regression slope
    differs from the given slope, given
- `xs` Array of x, can be NaN
- `ys` Array of y, can be NaN
- `slope` Slope to test
- `intercept` Optional: if true use intercept in the fitting
"""
function test_slope(
            xs::Array,
            ys::Array;
            slope::Number = 0,
            intercept::Bool = true
)
    # filter out NaN from the lists
    new_x = Float64[];
    new_y = Float64[];
    for (x,y) in zip(xs, ys)
        if ~isnan(x) && ~isnan(y)
            push!(new_x, x);
            push!(new_y, y - slope*x);
        end
    end

    # put new_x and new_y into a DataFrame
    df = DataFrame(X=new_x, Y=new_y);

    # run the fitting
    lr = linear_df_xy(df, intercept);

    # calculate the p value
    if intercept
        _slope_p = coeftable(lr).cols[4][2];
    else
        _slope_p = coeftable(lr).cols[4][1];
    end

    return _slope_p
end


text:
# Implicit types
i = 2
v = 5.
x = .05
println("i=$i, v=$v, x=$x")
# Type annotation
let
    z::Int8    = 10
    y::Float16 = 5.5
    print("z=$z, y=$y")
end

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
module TestPopfirst

include("preamble.jl")

@testset begin
    @test popfirst!!([0, 1]) == ([1], 0)
    @test popfirst!!((0, 1)) === ((1,), 0)
    @test popfirst!!((a=0, b=1)) === ((b=1,), 0)
    @test popfirst!!(SVector(0, 1)) === (SVector(1), 0)
end

@testset "mutation" begin
    @testset for xs in [
        [0, 1],
    ]
        @test popfirst!!(xs)[1] === xs
    end
end

end  # module


text:
using Test
using Verificatum: unmarshal_publickey, ElGamal, Enc, Tree, decode #, encode

PUBLIC_KEY = "$(@__DIR__)/publicKey"
CIPHERTEXTS = "$(@__DIR__)/ciphertexts"

pk, g = let
    bytes = read(PUBLIC_KEY)
    tree = decode(bytes)
    unmarshal_publickey(tree)
end

enc = Enc(pk, g)

𝐦 = g .^ (2:11)
𝐫 = 102:111

𝐞 = enc(𝐦, 𝐫)

tree = Tree(𝐞)
write(CIPHERTEXTS, tree)

### Checking that writing was succesful

bytes = read(CIPHERTEXTS)
tree′ = decode(bytes)

𝓖 = group(g)
𝐞′ = convert(ElGamal{Generator{𝓖}}, tree′)

@test 𝐞 == 𝐞′

#dec = Dec(sk)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Box and Ball System (BBS)
#
# ## References
#
# - [箱玉系の数理](https://mathsoc.jp/meeting/kikaku/2008aki/2008_aki_tokihiro.pdf)
# - https://ci.nii.ac.jp/naid/110009579038

# ## Package のインポート

using Colors

# +
function case1()
    n=10
    river = zeros(Bool, n)
    river[2]=1
    river[5]=1
    river
end

function case2()
    n=10
    river = zeros(Bool, n)
    river[2]=1
    river[3]=1
    river[5]=1
    river
end

function case3()
    n = 30
    river = zeros(Bool, n)
    river[2]=1
    river[3]=1
    river[4]=1

    river[9]=1
    river[10]=1

    river[15]=1
    river
end

# -

river = case1()

findfirst(river)

function evolve!(river::AbstractVector{Bool})
    T = eltype(river)
    n = count(river)
    bs = findall(river)
    negative = .! river
    offset=1
    for i in 1:n
        b = bs[i]
        p = findfirst(negative[b+1:end])
        if !isnothing(p)
            ∅ = p + b
            river[∅] = one(T)
        end
        river[b] = zero(T)
        negative = .! river
    end
    return river
end

river = case1()
expected = Bool[0,0,1,0,0,1,0,0,0,0]
next = evolve!(river)
@assert expected == next "え？ちゃうで？ expected $expected, actual $river"

river = case2()

river = case2()
expected = Bool[0,0,0,1,0,1,1,0,0,0]
next = evolve!(river)
@assert expected == next "え？ちゃうで？ expected $expected, actual $river"

river = case2()
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,1,0,0,1,1,0]
@assert expected == next

river = case2()
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,0,1,0,0,0,1]
@assert expected == next

river = case2()
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,0,0,1,0,0,0]
@assert expected == next

river = case2()
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,0,0,0,1,0,0]
@assert expected == next

river = case2()
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,0,0,0,0,1,0]
@assert expected == next

river = case2()
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,0,0,0,0,0,1]
@assert expected == next

river = case2()
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
next = evolve!(river)
expected = Bool[0,0,0,0,0,0,0,0,0,0]
@assert expected == next

function case4()
    n = 100
    river = zeros(Bool, n)
    river[2]=1
    river[3]=1
    river[4]=1

    river[25]=1
    river[26]=1

    river[50]=1
    river
end

# +
function play()
    river = case4()
    n = length(river)
    grid=zeros(Bool,n,n)
    for i ∈ 1:n
        IJulia.clear_output(true)
        grid[i,:] = evolve!(river)
        sleep(0.5)
        display(Gray.(.!grid))
    end
end

play()


text:
include("softmax.jl")
include("egreedy.jl")
include("normal.jl")
export LinearSoftmaxPolicy, LinearNormalPolicy
export EpsilonGreedy

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
@info "Running tests"
using MonteCarloMeasurements, Distributions
using Test, LinearAlgebra, Statistics, Random, GenericLinearAlgebra
import MonteCarloMeasurements: ⊗, gradient, optimize, DEFAULT_NUM_PARTICLES
@info "import Plots"
import Plots
@info "import Plots done"

Random.seed!(0)

@testset "MonteCarloMeasurements.jl" begin
    @info "Testing MonteCarloMeasurements"

    # σ/√N = σm
    @time @testset "sampling" begin
        @info "Testing sampling"
        for _ = 1:10
            @test -3 < mean(systematic_sample(100))*sqrt(100) < 3
            @test -3 < mean(systematic_sample(10000))*sqrt(10000) < 3
            @test -0.9 < std(systematic_sample(100)) < 1.1
            @test -0.9 < std(systematic_sample(10000)) < 1.1
        end
        params = Distributions.params
        @test systematic_sample(10000, Normal(1,1)) |> Base.Fix1(fit, Normal) |> params |> x-> all(isapprox.(x,(1,1), atol=0.1))
        systematic_sample(10000, Gamma(1,1)) #|> Base.Fix1(fit, Gamma)
        systematic_sample(10000, TDist(1)) #|> Base.Fix1(fit, TDist)
        @test systematic_sample(10000, Beta(1,1)) |> Base.Fix1(fit, Beta) |> params |> x-> all(isapprox.(x,(1,1), atol=0.1))

        for i = 1:100
            @test MonteCarloMeasurements.ess(Particles(10000)) > 7000
            x = randn(5000)
            v = vec([x';x'])
            @test 3000 < MonteCarloMeasurements.ess(Particles(v)) < 5300
        end
    end

    @time @testset "Basic operations" begin
        @info "Creating the first StaticParticles"
        @test 0 ∓ 1 isa StaticParticles
        @test [0,0] ∓ 1. isa MonteCarloMeasurements.MvParticles
        @test [0,0] ∓ [1.,1.] isa MonteCarloMeasurements.MvParticles
        @test -50 ⊞ Normal(0,1) ≈ -50 ± 1
        @test 10 ⊠ Normal(0,1) ≈ 10*Particles(Normal(0,1))

        @info "Done"
        PT = Particles
        for PT = (Particles, StaticParticles)
            @testset "$(repr(PT))" begin
                @info "Running tests for $PT"
                p = PT(100)
                @test_nowarn MonteCarloMeasurements.shortform(p)
                @test_nowarn println(p)
                @test (p+p+p).particles ≈ 3p.particles # Test 3arg operator
                @test (p+p+1).particles ≈ 1 .+ 2p.particles # Test 3arg operator
                @test (1+p+1).particles ≈ 2 .+ p.particles # Test 3arg operator
                @test (p+1+p).particles ≈ 1 .+ 2p.particles # Test 3arg operator
                @test 0 ± 1 ≈ p
                @test 0 ∓ 1 ≈ p
                @test sum(p) ≈ 0
                @test pcov(p) ≈ 1 atol=0.2
                @test pstd(p) ≈ 1 atol=0.2
                @test pvar(p) ≈ 1 atol=0.2
                @test meanvar(p) ≈ 1/(nparticles(p)) atol=5e-3
                @test meanstd(p) ≈ 1/sqrt(nparticles(p)) atol=5e-3
                @test minmax(1+p,p) == (p, 1+p)
                b = PT(100)
                uvec = unique([p, p, b, p, b, b]) # tests hash
                @test length(uvec) == 2
                @test p ∈ uvec
                @test b ∈ uvec

                @test !(p ≲ p)
                @test !(p ≳ p)
                @test (p ≲ 2.1)
                @test !(p ≲ 1.9)
                @test (p ≳ -2.1)
                @test !(p ≳ -1.9)
                @test (-2.1 ≲ p)
                @test !(-1.9 ≲ p)
                @test (2.1 ≳ p)
                @test !(1.9 ≳ p)
                @test p ≈ p
                @test p ≈ 0
                @test 0 ≈ p
                @test p != 0
                @test p != 2p
                @test p ≈ 1.9pstd(p)
                @test !(p ≈ 2.1pstd(p))
                @test !(p ≉ p)
                @test !(pmean(p) ≉ p)
                @test p ≉ 2.1pstd(p)
                @test !(p ≉ 1.9pstd(p))
                @test (5 ± 0.1) ≳ (1 ± 0.1)
                @test (1 ± 0.1) ≲ (5 ± 0.1)

                a = rand()
                pa = Particles([a])
                @test a == pa
                @test a ≈ pa
                @test pa ≈ a

                v = randn(5)
                @test Vector(PT(v)) == v
                @test Array(PT(v)) == v

                @testset "unsafe comparisons" begin
                    unsafe_comparisons(false)
                    @test_throws ErrorException p<p
                    @test_throws ErrorException p>p
                    @test_throws ErrorException p>=p
                    @test_throws ErrorException p<=p

                    for mode in (:montecarlo, :reduction, :safe)
                        @show mode
                        unsafe_comparisons(mode, verbose=false)
                        @test p<100+p
                        @test p+100>p
                        @test p+100>=p
                        @test p<=100+p

                        @test p<100
                        @test 100>p
                        @test 100>=p
                        @test p<=100
                        @unsafe begin
                            @test -10 < p
                            @test p <= p
                            @test p >= p
                            @test !(p < p)
                            @test !(p > p)
                            @test (p < 1+p)
                            @test (p+1 > p)
                        end
                    end

                    @test_throws ErrorException p<p
                    @test_throws ErrorException p>p
                    @test_throws ErrorException @unsafe error("") # Should still be safe after error

                    @test_throws ErrorException p>=p
                    @test_throws ErrorException p<=p
                    unsafe_comparisons(:montecarlo, verbose=false)
                    @test p>=p
                    @test p<=p
                    @test !(p<p)
                    @test !(p>p)
                    @test_throws ErrorException p < Particles(p.particles[randperm(nparticles(p))])
                    unsafe_comparisons(false)

                    @unsafe tv = 2
                    @test tv == 2
                    @unsafe tv1,tv2 = 1,2
                    @test (tv1,tv2) == (1,2)
                    @unsafe tv3,tv4 = range(1, stop=3, length=5), range(1, stop=3, length=5)
                    @test (tv3,tv4) == (range(1, stop=3, length=5),range(1, stop=3, length=5))
                    @test MonteCarloMeasurements.COMPARISON_FUNCTION[] == pmean
                    set_comparison_function(pmedian)
                    @test MonteCarloMeasurements.COMPARISON_FUNCTION[] == pmedian
                    cp = PT(10)
                    cmp = @unsafe p < cp
                    @test cmp == (pmedian(p) < pmedian(cp))
                    set_comparison_function(pmean)
                end

                f = x -> 2x + 10
                @test 9.6 < pmean(f(p)) < 10.4
                # @test 9.6 < f(p) < 10.4
                @test f(p) ≈ 10
                @test !(f(p) ≲ 11)
                @test f(p) ≲ 15
                @test 5 ≲ f(p)
                @test 1 ≲ 2

                @test Normal(f(p)).μ ≈ pmean(f(p))
                @test fit(Normal, f(p)).μ ≈ pmean(f(p))

                f = x -> x^2
                p = PT(100)
                @test 0.9 < pmean(f(p)) < 1.1
                @test 0.9 < pmean(f(p)) < 1.1
                @test f(p) ≈ 1
                @test !(f(p) ≲ 1)
                @test f(p) ≲ 5
                @test -3 ≲ f(p)
                @test MvNormal([f(p),p]) isa MvNormal

                A = randn(3,3) .+ [PT(100) for i = 1:3, j = 1:3]
                a = [PT(100) for i = 1:3]
                b = [PT(100) for i = 1:3]
                @test sum(a.*b) ≈ 0
                @test all(A*b .≈ [0,0,0])
                @test A*b .+ 1 ≈ [1,1,1]
                @test [1,1,1] ≈ A*b .+ 1

                @test all(A\b .≈ zeros(3))
                @test_nowarn @unsafe qr(A)
                @test_nowarn Particles(100, MvNormal(2,1)) ./ Particles(100, Normal(2,1))
                pn = Particles(100, Normal(2,1), systematic=false)
                @test pn ≈ 2
                @test !issorted(pn.particles)
                @test !issorted(p.particles)

                pn = Particles(100, Normal(2,1), systematic=true, permute=false)
                @test pn ≈ 2
                @test issorted(pn.particles)

                rng = MersenneTwister(657)
                pn1 = Particles(rng, 100, Normal(2,1), systematic=true, permute=true)
                rng = MersenneTwister(657)
                pn2 = Particles(rng, 100, Normal(2,1), systematic=true, permute=true)
                @test pn1 == pn2

                rng = MersenneTwister(27)
                pn1 = Particles(rng, Normal(2,1), systematic=true, permute=true)
                rng = MersenneTwister(27)
                pn2 = Particles(rng, Normal(2,1), systematic=true, permute=true)
                @test pn1 == pn2

                rng = MersenneTwister(932)
                pn1 = Particles(rng, 100, systematic=true, permute=true)
                rng = MersenneTwister(932)
                pn2 = Particles(rng, 100, systematic=true, permute=true)
                @test pn1 == pn2

                @test PT(Float32) isa PT{Float32}
                @test PT(Float64) isa PT{Float64}
                @test PT(Float32, 250) isa PT{Float32, 250}
                @test PT(Float32, 250, Normal(0.1f0)) isa PT{Float32, 250}
                @test_throws ArgumentError PT(Float32, 250, Gamma(0.1))

                @info "Tests for $PT done"

                p = PT{Float64,10}(2)
                @test p isa PT{Float64,10}
                @test all(p.particles .== 2)

                @test Particles(100) + Particles(randn(Float32, 100)) ≈ 0
                @test_throws MethodError p + Particles(randn(Float32, 200)) # Npart and Float type differ
                @test_throws MethodError p + Particles(200) # Npart differ

                mp = MvParticles([randn(10) for _ in 1:3])
                @test length(mp) == 10
                @test nparticles(mp) == 3

                v = [(1,2), (3,4)]
                pv = MvParticles(v)
                @test length(pv) == 2
                @test pv[1].particles[1] == v[1][1]
                @test pv[1].particles[2] == v[2][1]
                @test pv[2].particles[1] == v[1][2]
                @test pv[2].particles[2] == v[2][2]
                @test length(pv) == 2

                v = [(a=1,b=randn(2)), (a=3,b=randn(2))]
                pv = MvParticles(v)
                @test length(pv) == 2
                @test pv.a.particles[1] == v[1].a
                @test pv.a.particles[2] == v[2].a
                @test pv.b == Particles([v[1].b v[2].b]')
                @test length(pv) == 2

                @testset "discrete distributions" begin
                    p = PT(Poisson(50))
                    @test p isa PT{Int}
                    @test (p^2 - 1) isa PT{Int}
                    @test exp(p) isa PT{Float64}

                    # mainly just test that printing PT{Int} doesn't error
                    io = IOBuffer()
                    show(io, p)
                    s = String(take!(io))
                    @test occursin('±', s)

                    # issue #50
                    @test 2.5 * p isa PT{Float64}
                    @test p / 3 isa PT{Float64}
                    @test sqrt(p) isa PT{Float64}
                end
            end
        end
    end

    @time @testset "Multivariate Particles" begin
        for PT = (Particles, StaticParticles)
            @testset "$(repr(PT))" begin
                @info "Running tests for multivariate $PT"
                p = PT(100, MvNormal(2,1))
                @test_nowarn sum(p)
                @test pcov(p) ≈ I atol=0.6
                @test pcor(p) ≈ I atol=0.6
                @test pmean(p) ≈ [0,0] atol=0.2
                m = Matrix(p)
                @test size(m) == (100,2)
                # @test m[1,2] == p[1,2]

                p = PT(100, MvNormal(2,2))
                @test pcov(p) ≈ 4I atol=2
                @test [0,0] ≈ pmean(p) atol=1
                @test size(Matrix(p)) == (100,2)

                p = PT(100, MvNormal(2,2))
                @test fit(MvNormal, p).μ ≈ pmean(p)
                @test MvNormal(p).μ ≈ pmean(p)
                @test cov(MvNormal(p)) ≈ pcov(p)
                @info "Tests for multivariate $PT done"
            end
        end
    end

    @testset "sigmapoints" begin
        @info "Testing sigmapoints"
        m = 1
        Σ = 3
        s = sigmapoints(m,Σ)
        @test var(s) ≈ Σ
        @test mean(s) == m
        @test sigmapoints(Normal(m,√(Σ))) == s

        m = [1,2]
        Σ = [3. 1; 1 4]
        s = sigmapoints(m,Σ)
        @test cov(s) ≈ Σ
        @test mean(s, dims=1)' ≈ m
        @test sigmapoints(MvNormal(m,Σ)) == s
    end

    @testset "transform_moments" begin
        m, Σ   = [1,2], [2 1; 1 4] # Desired mean and covariance
        C = randn(2,2)
        C = cholesky(C'C + 5I).L
        particles = transform_moments((C*randn(2,DEFAULT_NUM_PARTICLES))', m, Σ)
        @test mean(particles, dims=1)[:] ≈ m
        @test cov(particles) ≈ Σ
        particles = transform_moments((C*randn(2,DEFAULT_NUM_PARTICLES))', m, Σ, preserve_latin=true)
        @test mean(particles, dims=1)[:] ≈ m
        @test Diagonal(cov(particles)) ≈ Diagonal(Σ) atol=2
    end

    @time @testset "gradient" begin
        @info "Testing gradient"
        e = 0.001
        p = 3 ± e
        f = x -> x^2
        fp = f(p)
        @test gradient(f,p)[1] ≈ 6 atol=1e-4
        @test gradient(f,p)[2] ≈ 2e atol=1e-4
        # @test gradient(f,3) > 6 # Convex function
        @test gradient(f,3) ≈ 6

        A = randn(3,3)
        H = A'A
        h = randn(3)
        c = randn()
        @assert isposdef(H)
        f = x -> (x'H*x + h'x) + c
        j = x -> H*x + h

        e = 0.001
        x = randn(3)
        xp = x ± e
        g = 2H*x + h
        @test MonteCarloMeasurements.gradient(f,xp) ≈ g atol = 0.1
        @test MonteCarloMeasurements.jacobian(j,xp) ≈ H
    end

    @time @testset "leastsquares" begin
        @info "Testing leastsquares"
        n, m = 10000, 3
        A = randn(n,m)
        x = randn(m)
        y = A*x
        σ = 0.1
        yn = y .+ σ.*randn()
        # xh = A\y
        C1 = σ^2*inv(A'A)

        yp = yn .+ σ.*Particles.(2000)
        xhp = (A'A)\(A'yp)
        @test xhp ≈ A\yp
        @test sum(abs, tr((pcov(xhp) .- C1) ./ abs.(C1))) < 0.2

        @test norm(pcov(xhp) .- C1) < 1e-7
        @test xhp ≈ x
        @test pmean(xhp) ≈ x atol=3sum(sqrt.(diag(C1)))
        yp = nothing; GC.gc(true) # make sure this big matrix is deallocated
    end

    @time @testset "misc" begin
        @info "Testing misc"
        p = 0 ± 1
        @test MonteCarloMeasurements.particletypetuple(p) == (Float64, DEFAULT_NUM_PARTICLES, Particles)
        @test MonteCarloMeasurements.particletypetuple(typeof(p)) == (Float64, DEFAULT_NUM_PARTICLES, Particles)
        @test_nowarn display(p)
        @test_nowarn println(p)
        @test_nowarn show(stdout, MIME"text/x-latex"(), p); println()
        @test_nowarn println(0p)
        @test_nowarn show(stdout, MIME"text/x-latex"(), 0p); println()
        @test_nowarn show(stdout, MIME"text/x-latex"(), p + im*p); println()
        @test_nowarn show(stdout, MIME"text/x-latex"(), p - im*p); println()
        @test_nowarn show(stdout, MIME"text/x-latex"(), im*p); println()
        @test_nowarn show(stdout, MIME"text/x-latex"(), -im*p); println()

        @test_nowarn show(stdout, MIME"text/plain"(), p); println()
        @test_nowarn println(0p)
        @test_nowarn show(stdout, MIME"text/plain"(), 0p); println()
        @test_nowarn show(stdout, MIME"text/plain"(), p + im*p); println()
        @test_nowarn show(stdout, MIME"text/plain"(), p - im*p); println()
        @test_nowarn show(stdout, MIME"text/plain"(), im*p); println()
        @test_nowarn show(stdout, MIME"text/plain"(), -im*p -10im); println()

        @test_nowarn display([p, p])
        @test_nowarn println([p, p])
        @test_nowarn println([p, 0p])

        @test Particles{Float64,DEFAULT_NUM_PARTICLES}(p) == p
        @test Particles{Float64,5}(0) == 0*Particles(5)
        @test length(Particles(100, MvNormal(2,1))) == 2
        @test nparticles(p) == DEFAULT_NUM_PARTICLES
        @test ndims(p) == 0
        @test particleeltype(p) == Float64
        @test eltype(typeof(p)) == typeof(p)
        @test eltype(p) == typeof(p)
        @test convert(Int, 0p) == 0
        @test promote_type(Particles{Float64,10}, Float64) == Particles{Float64,10}
        @test promote_type(Particles{Float64,10}, Int64) == Particles{Float64,10}
        @test promote_type(Particles{Float64,10}, ComplexF64) == Complex{Particles{Float64,10}}
        @test promote_type(Particles{Float64,10}, Missing) == Union{Particles{Float64,10},Missing}
        @testset "promotion of $PT" for PT in (Particles, StaticParticles)
            @test promote_type(PT{Float64,10}, PT{Float64,10}) == PT{Float64,10}
            @test promote_type(PT{Float64,10}, PT{Int,10}) == PT{Float64,10}
            @test promote_type(PT{Int,5}, PT{Float64,10}) == PT

            @test promote_type(PT{Float64,10}, PT{Float32,10}) == PT{Float64,10}
            @test promote_type(StaticParticles{Float64,10}, PT{Float32,10}) == StaticParticles{Float64,10}
        end
        @test promote_type(Particles{Float64,10}, StaticParticles{Float64,10}) == StaticParticles{Float64,10}
        @test promote_type(Particles{Int,10}, StaticParticles{Float64,10}) == StaticParticles{Float64,10}
        @test promote_type(Particles{Float64,10}, StaticParticles{Float32,10}) == StaticParticles{Float64,10}
        @test promote_type(Particles{Int,10}, StaticParticles{Float32,10}) == StaticParticles{Float32,10}
        @test convert(Float64, 0p) isa Float64
        @test convert(Float64, 0p) == 0
        @test convert(Int, 0p) isa Int
        @test convert(Int, 0p) == 0
        @test convert(Particles{Float64,100}, Particles(randn(Float32, 100))) isa Particles{Float64,100}
        @test_throws ArgumentError convert(Int, p)
        @test_throws ArgumentError AbstractFloat(p)
        @test AbstractFloat(0p) == 0.0
        @test Particles(500) + Particles(randn(Float32, 500)) isa typeof(Particles(500))
        @test isfinite(p)
        @test iszero(0p)
        @test iszero(p, 0.1)
        @test !iszero(p)
        @test !(!p)
        @test !(0p)
        @test round(p) ≈ 0 atol=0.1
        @test round(Int,p) == 0
        @test round(Int,p) isa Int
        @test sincos(p) == (sin(p), cos(p))
        @test norm(p) == abs(p)
        @test pmean(norm([p,p]) - sqrt(2p^2)) < sqrt(eps()) # ≈ with atol fails on mac
        @test pmean(LinearAlgebra.norm2([p,p]) - sqrt(2p^2)) < sqrt(eps()) # ≈ with atol fails on mac
        @test MvNormal(Particles(500, MvNormal(2,1))) isa MvNormal
        @test eps(typeof(p)) == eps(Float64)
        @test eps(p) == eps(Float64)
        A = randn(2,2)
        B = A .± 0
        @test sum(pmean, exp(A) .- exp(B)) < 1e-9
        @test sum(pmean, abs.(log(A)) .- abs.(log(B))) < 1e-9
        @test sum(pmean, abs.(eigvals(A)) .- abs.(eigvals(B))) < 1e-9

        @test @unsafe pmean(sum(abs, sort(eigvals([0 1 ± 0.001; -1. 0]), by=imag) - [ 0.0 - (1.0 ± 0.0005)*im
                                                0.0 + (1.0 ± 0.0005)*im])) < 0.002

        e = eigvals([1 ± 0.001 0; 0 1.])
        @test e isa Vector{Complex{Particles{Float64, DEFAULT_NUM_PARTICLES}}}
        @test all(isapprox.(e, [1.0 ± 0.00058, 1.0 ± 0.00058], atol=1e-2))

        @test (1 .. 2) isa Particles
        @test std(diff(sort((1 .. 2).particles))) < sqrt(eps())
        @test pmaximum((1 .. 2)) <= 2
        @test pminimum((1 .. 2)) >= 1

        pp = [1. 0; 0 1] .± 0.0
        @test lyap(pp,pp) == lyap([1. 0; 0 1],[1. 0; 0 1])

        @test intersect(p,p) == union(p,p)
        @test nparticles(intersect(p, 1+p)) < 2nparticles(p)
        @test nparticles(union(p, 1+p)) == 2nparticles(p)

        p = 2 ± 0
        q = 3 ± 0
        @test sqrt(complex(p,p)) == sqrt(complex(2,2))
        @test exp(complex(p,p)) == exp(complex(2,2))
        @test sqrt!(fill(complex(1.,1.), DEFAULT_NUM_PARTICLES), complex(p,p)) == sqrt(complex(2,2))
        @test exp!(fill(complex(1.,1.), DEFAULT_NUM_PARTICLES), complex(p,p)) == exp(complex(2,2))
        y = Particles(100)
        @test exp(im*y) ≈ cos(y) + im*sin(y)
        @test complex(p,p)/complex(q,q) == complex(2,2)/complex(3,3)
        @test p/complex(q,q) == 2/complex(3,3)
        @test 2/complex(q,q) == 2/complex(3,3)
        @test !isinf(complex(p,p))
        @test isfinite(complex(p,p))

        z = complex(1 ± 0.1, 1 ± 0.1)
        @unsafe @test abs(sqrt(z ^ 2) - z) < eps()
        @unsafe @test abs(sqrt(z ^ 2.0) - z) < eps()
        @test z/z ≈ 1

        z = complex(1 ± 0.1, 0 ± 0.1)
        @test real(2 ^ z) ≈ 2 ^ real(z)
        @test real(2.0 ^ z) ≈ 2.0 ^ real(z)
        @test real(z ^ z) ≈ real(z) ^ real(z)
        p = 2 ± 0.1
        q = 3 ± 0.1
        @test wasserstein(p,p,1) == 0
        @test wasserstein(p,q,1) >= 0
        @test bootstrap(p) ≈ p
        rng = MersenneTwister(453)
        p1 = bootstrap(rng,p)
        rng = MersenneTwister(453)
        p2 = bootstrap(rng,p)
        @test p1 == p2
        @test nparticles(bootstrap(p, 10)) == 10
        @test_nowarn bootstrap([p; p])
    end

    @time @testset "mutation" begin
        @info "Testing mutation"
        function adder!(x)
            for i = eachindex(x)
                x[i] += 1
            end
            x
        end
        x = (1:5) .± 1
        adder!(x)
        @test x ≈ ((2:6) .± 1)
    end

    @time @testset "outer_product" begin
        @info "Testing outer product"
        d = 2
        μ = zeros(d)
        σ = ones(d)
        p = μ ⊗ σ
        @test length(p) == 2
        @test nparticles(p[1]) <= 100_000
        @test pcov(p) ≈ I atol=1e-1
        p = μ ⊗ 1
        @test length(p) == 2
        @test nparticles(p[1]) <= 100_000
        @test pcov(p) ≈ I atol=1e-1
        p = 0 ⊗ σ
        @test length(p) == 2
        @test nparticles(p[1]) <= 100_000
        @test pcov(p) ≈ I atol=1e-1

        rng = MersenneTwister(38)
        p1 = outer_product(rng, Normal.(μ,σ))
        rng = MersenneTwister(38)
        p2 = outer_product(rng, Normal.(μ,σ))
        @test p1 == p2
    end

    @time @testset "plotting" begin
        @info "Testing plotting"
        p = Particles(100)
        v = randn(3) .+ Particles.(10)
        M = randn(3,2) .+ Particles.(10)
        @test_nowarn Plots.plot(p)
        @test_nowarn Plots.plot(v)
        @test_nowarn Plots.plot(v, N=8)
        @test_nowarn Plots.plot(M, N=10)
        @test_nowarn Plots.plot(x->x^2,v)
        @test_nowarn Plots.plot(v,v)
        @test_nowarn Plots.plot(v,v, N=10)
        @test_nowarn Plots.plot(v,v; points=true)
        @test_nowarn Plots.plot(v,ones(3))
        @test_nowarn Plots.plot(1:3,v)
        @test_nowarn Plots.plot(1:3, v, N=10)
        @test_nowarn Plots.plot(1:3, M, N=10)

        @test_nowarn errorbarplot(1:3,v)
        @test_nowarn errorbarplot(1:3,[v v])
        @test_nowarn mcplot(1:3,v)
        @test_nowarn mcplot(1:3,v, 10)
        @test_nowarn mcplot(1:3,[v v])
        @test_nowarn ribbonplot(1:3,v)
        @test_nowarn ribbonplot(1:3,v,(0.1,0.9))

        @test_nowarn errorbarplot(v)
        @test_nowarn errorbarplot([v v])
        @test_nowarn mcplot(v)
        @test_nowarn mcplot(v, 10)
        @test_nowarn mcplot([v v])
        @test_nowarn ribbonplot(v)
        @test_nowarn ribbonplot(v,(0.1,0.9))
        @test_nowarn ribbonplot(v, N = 2)

        @test_nowarn errorbarplot(v, 0.1)
        @test_nowarn errorbarplot([v v], 0.1)
        @test_nowarn mcplot(v, 0.1)
        @test_nowarn ribbonplot(v, 0.1)

        @test_throws ArgumentError errorbarplot(1:3, (1:3) .± 0.1, 1,1)

        @test_nowarn MonteCarloMeasurements.print_functions_to_extend()
    end

    @time @testset "optimize" begin
        @info "Testing optimization"
        function rosenbrock2d(x)
            return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
        end
        @test any(1:10) do i
            p = -1ones(2) .+ 2 .*Particles.(200) # Optimum is in [1,1]
            popt = optimize(rosenbrock2d, deepcopy(p))
            popt ≈ [1,1]
        end

        p = -1ones(2) .+ 2 .*Particles.(200) # Optimum is in [1,1]
        rng = MersenneTwister(876)
        popt1 = optimize(rng, rosenbrock2d, deepcopy(p))
        rng = MersenneTwister(876)
        popt2 = optimize(rng, rosenbrock2d, deepcopy(p))
        @test popt1 == popt2
    end

    @testset "Particle BLAS" begin
        @info "Testing Particle BLAS"
        p = ones(10) .∓ 1
        A = randn(20,10)
        @test pmean(sum(abs, A*p - MonteCarloMeasurements._pgemv(A,p))) < 1e-12
        @test pmean(sum(abs, A*Particles.(p) - A*p)) < 1e-12

        v = randn(10)
        @test pmean(sum(abs, v'p - MonteCarloMeasurements._pdot(v,p))) < 1e-12
        @test pmean(sum(abs, p'v - MonteCarloMeasurements._pdot(v,p))) < 1e-12
        @test pmean(sum(abs, v'*Particles.(p) - v'p)) < 1e-12

        @test pmean(sum(abs, axpy!(2.0,Matrix(p),copy(Matrix(p))) - Matrix(copy(axpy!(2.0,p,copy(p)))))) < 1e-12
        @test pmean(sum(abs, axpy!(2.0,Matrix(p),copy(Matrix(p))) - Matrix(copy(axpy!(2.0,p,copy(p)))))) < 1e-12

        y = randn(20) .∓ 1
        @test pmean(sum(abs, mul!(y,A,p) - mul!(Particles.(y),A,Particles.(p)))) < 1e-12

        for PT in (Particles, StaticParticles)
            for x in (1.0, 1 + PT()), y in (1.0, 1 + PT()), z in (1.0, 1 + PT())
                x == y == z == 1.0 && continue
                @test (x*y+z).particles ≈ muladd(x,y,z).particles
            end
        end

        #
        # @btime $A*$p
        # @btime _pgemv($A,$p)
        #
        # @btime sum($A*$p)
        # @btime sum(_pgemv($A,$p))
        #
        # @btime $v'*$p
        # @btime _pdot($v,$p)
        #
        # @btime sum($v'*$p)
        # @btime sum(_pdot($v,$p))

        # @btime mul!($y,$A,$p)
        # @btime MonteCarloMeasurements.pmul!($y,$A,$p)
        # 178.373 μs (6 allocations: 336 bytes)
        # 22.320 μs (0 allocations: 0 bytes)
        # 3.705 μs (0 allocations: 0 bytes)
    end

    @time @testset "bymap" begin
        @info "Testing bymap"

        p = 0 ± 1
        ps = 0 ∓ 1

        f(x) = 2x
        f(x,y) = 2x + y

        @test f(p) ≈ bymap(f,p)
        @test f(p) ≈ @bymap f(p)
        @test typeof(f(p)) == typeof(@bymap(f(p)))
        @test typeof(f(ps)) == typeof(@bymap(f(ps)))
        @test f(p,p) ≈ @bymap f(p,p)
        @test f(p,10) ≈ @bymap f(p,10)
        @test !(f(p,10) ≈ @bymap f(p,-10))

        @test f(p) ≈ @bypmap f(p)
        @test f(p,p) ≈ @bypmap f(p,p)
        @test f(p,10) ≈ @bypmap f(p,10)
        @test !(f(p,10) ≈ @bypmap f(p,-10))

        g(x,y) = sum(x) + sum(y)
        @test g([p,p], [p,p]) ≈ @bymap g([p,p], [p,p])
        @test g([p,p], p) ≈ @bymap g([p,p], p)
        @test g([p p], [p,p]) ≈ @bymap g([p p], [p,p])

        @test g([p,p], [p,p]) ≈ @bypmap g([p,p], [p,p])
        @test g([p,p], p) ≈ @bypmap g([p,p], p)
        @test g([p p], [p,p]) ≈ @bypmap g([p p], [p,p])

        h(x,y) = x .* y'
        Base.Cartesian.@nextract 4 p d-> 0±1
        @test all(h([p_1,p_2], [p_3,p_4]) .≈ bymap(h, [p_1,p_2], [p_3,p_4]))
        @test all(h([p_1,p_2], [p_3,p_4]) .≈ bypmap(h, [p_1,p_2], [p_3,p_4]))

        h2(x,y) = x .* y
        Base.Cartesian.@nextract 4 p d-> 0±1
        @test h2([p_1,p_2], [p_3,p_4]) ≈ @bymap  h2([p_1,p_2], [p_3,p_4])
        @test h2([p_1,p_2], [p_3,p_4]) ≈ @bypmap h2([p_1,p_2], [p_3,p_4])

        g(nt::NamedTuple) = nt.x^2 + nt.y^2
        @test g((x=p_1, y=p_2)) == p_1^2 + p_2^2

        g2(a,nt::NamedTuple) = a + nt.x^2 + nt.y^2
        @test g2(p_3, (x=p_1, y=p_2)) == p_3 + p_1^2 + p_2^2

        @test_throws ErrorException bymap(x->ones(3,3,3,3), p)
        @test_throws ErrorException bypmap(x->ones(3,3,3,3), p)
        @test MonteCarloMeasurements.arggetter(1,1) == 1

        @test MonteCarloMeasurements.particletype(p) == Particles{Float64,DEFAULT_NUM_PARTICLES}
        @test MonteCarloMeasurements.particletype(Particles{Float64,DEFAULT_NUM_PARTICLES}) == Particles{Float64,DEFAULT_NUM_PARTICLES}
        @test MonteCarloMeasurements.particletype([p,p]) == Particles{Float64,DEFAULT_NUM_PARTICLES}
        @test nparticles(p) == DEFAULT_NUM_PARTICLES

        @testset "@prob" begin
            @info "Testing @prob"
            p = Particles()
            q = Particles()

            @test mean((p).particles .< 1) == @prob p < 1
            @test mean((p+2).particles .< 1) == @prob p+2 < 1
            @test mean((2+p).particles .< 1) == @prob 2+p < 1
            @test mean((p+p).particles .< 1) == @prob p+p < 1
            @test mean((p).particles .< p.particles) == @prob p < p
            @test mean((p+1).particles .< p.particles) == @prob p+1 < p

            @test mean((p+q).particles .< 1) == @prob p+q < 1
            @test mean((p).particles .< q.particles) == @prob p < q
            @test mean((p+1).particles .< q.particles) == @prob p+1 < q

            @test mean((p+q).particles .> 1) == @prob p+q > 1
            @test mean((p).particles .> q.particles) == @prob p > q
            @test mean((p+1).particles .> q.particles) == @prob p+1 > q

            @test mean((p+q).particles .>= 1) == @prob p+q >= 1
            @test mean((p).particles .>= q.particles) == @prob p >= q
            @test mean((p+1).particles .>= q.particles) == @prob p+1 >= q

            @test mean((abs(p)).particles .> sin(q).particles) == @prob abs(p) > sin(q)

        end

    end

    @testset "inference" begin
        @inferred zero(Particles{Float64,1})
        @inferred zeros(Particles{Float64,1}, 5)
        @inferred bymap(sin, 1 ± 2)
    end

    @testset "nominal values" begin
        @info "Testing nominal values"

        p = 1 ± 0.1
        n = 0
        pn = with_nominal(p, n)
        @test nominal(pn) == pn.particles[1] == n
        @test nominal(p) != n

        p = [p, p]
        n = [0, 1]
        pn = with_nominal(p, n)
        @test nominal(pn) == MonteCarloMeasurements.vecindex.(pn, 1) == n
        @test nominal(p) != n

        p = 1 ∓ 0.1
        n = 0
        pn = with_nominal(p, n)
        @test nominal(pn) == pn.particles[1] == n
        @test nominal(p) != n

        p = [p, p]
        n = [0, 1]
        pn = with_nominal(p, n)
        @test nominal(pn) == MonteCarloMeasurements.vecindex.(pn, 1) == n
        @test nominal(p) != n

        P = complex(1 ± 0.1, 2 ± 0.1)
        @test nominal(P) == complex(real(P).particles[1], imag(P).particles[1])

    end

    include("test_unitful.jl")
    include("test_forwarddiff.jl")
    include("test_deconstruct.jl")
    include("test_sleefpirates.jl")
    include("test_measurements.jl")

end

# These can not be inside a testset, causes "testf not defined"
testf(x,y) = sum(x+y)
@test_nowarn register_primitive(testf)
p = 1 ± 0.1
@test testf(p,p) == sum(p+p)

# Integration tests and bechmarks

# using BenchmarkTools
# A = [StaticParticles(100) for i = 1:3, j = 1:3]
# B = similar(A, Float64)
# @btime qr($(copy(A)))
# @btime map(_->qr($B), 1:100);

#
# # Benchmark and comparison to Measurements.jl
# using BenchmarkTools, Printf, ControlSystems
# using MonteCarloMeasurements, Measurements
# using Measurements: ±
# using MonteCarloMeasurements: ∓
# w = exp10.(LinRange(-0.7,0.3,50))
#
# p = 1 ± 0.1
# ζ = 0.3 ± 0.1
# ω = 1 ± 0.1
# Gm = tf([p*ω], [1, 2ζ*ω, ω^2])
# # tm = @belapsed bode($Gm,$w)
#
# p = 1 ∓ 0.1
# ζ = 0.3 ∓ 0.1
# ω = 1 ∓ 0.1
# Gmm = tf([p*ω], [1, 2ζ*ω, ω^2])
# # tmm = @belapsed bode($Gmm,$w)
#
# σquant = 1-(cdf(Normal(0,1), 1)-cdf(Normal(0,1), -1))
#
# magm = bode(Gm,w)[1][:]
# magmm = bode(Gmm,w)[1][:]
# errorbarplot(w,magmm, σquant/2, xscale=:log10, yscale=:log10, lab="Particles", linewidth=2)
# plot!(w,magm, lab="Measurements")


text:
using MyterialColors
using Test

# ------------------------------- utility code ------------------------------- #

"""
Lenient comparison operator for `struct`, both mutable and immutable (type with \\eqsim).

From: https://discourse.julialang.org/t/what-general-purpose-commands-do-you-usually-end-up-adding-to-your-projects/4889/2
"""
@generated function ≂(x, y)
    if !isempty(fieldnames(x)) && x == y
        mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), fieldnames(x))
    else
        :(x == y)
    end
end

"""
    Checks if two palettes have the same colors
"""
function palette_equal(p1::Palette, p2::Palette)
    same_colors = map((x) -> x[1] ≂ x[2], zip(p1.colors, p2.colors))
    return sum(same_colors) == length(same_colors)
end

# ---------------------------------------------------------------------------- #
#                                 BASE IMPORTS                                 #
# ---------------------------------------------------------------------------- #
@testset "MyterialColors.jl" begin
    # test we can import color codes
    @test typeof(MyterialColors.salmon) == String

end

# ---------------------------------------------------------------------------- #
#                                  color types                                 #
# ---------------------------------------------------------------------------- #
@testset "Colors Types" begin
    # ------------------------------- constructors ------------------------------- #
    rgb1 = RGB(255, 255, 255)
    rgb2 = RGB(1.0, 1.0, 1.0)
    rgb3 = RGB((255, 255, 255))

    @test rgb1  ≂ rgb2
    @test rgb2  ≂ rgb3

    h = Hex("#ffffff")
    @test h.code == "ffffff"
    @test h.string == "#ffffff"

    @test_throws AssertionError Hex("ff")
    @test_throws AssertionError Hex("#ff")

    # -------------------------------- coinversion ------------------------------- #
    rgb1_as_hex = rgb2hex(rgb1)
    rgb1_v2_as_hex = rgb2hex(255, 255, 255)

    @test rgb1_as_hex ≂ h
    @test rgb1_as_hex ≂ rgb1_v2_as_hex

    h_as_rgb = hex2rgb(h)

    @test h_as_rgb ≂ rgb1

    h_as_rgb_as_h = rgb2hex(h_as_rgb)
    @test  h_as_rgb_as_h.code == h.code
    @test  h_as_rgb_as_h ≂ h

end

# ---------------------------------------------------------------------------- #
#                                   palettte                                   #
# ---------------------------------------------------------------------------- #
# @testset "Palette" begin
#     # create a palette from two colors
#     rgb1 = RGB(0, 0, 0)
#     rgb2 = RGB(255, 255, 255)
#     palette = Palette(rgb1, rgb2, 100)

#     @test length(palette.colors) == 100
#     @test palette.colors[1] ≂ rgb1
#     @test palette.colors[end] ≂ rgb2

#     # check that palette equals is working
#     @test palette_equal(palette, palette)
# end


text:
"""

```
distance(x,y)
```

Computes the minimum distance between two sets of atoms, between an atom and a set of atoms, or simply the distance between two atoms. The input may be a vector of `Atom`s, or the matrices that are output of the `coor` function. Using the matrices may result in greater speed, particularly if the greatest dimension of the matrices is the first one (`xyz_in_cols=true` option, which is the default option).

### Examples

```julia-repl
julia> model = wget("1BSX");

julia> protein = select(model,"protein");

julia> ligand = select(model,"resname T3");

julia> distance(protein,ligand)
2.7775834820937417

julia> distance(protein[1],ligand[3])
36.453551075306784

julia> using BenchmarkTools

julia> @btime distance(\$ligand[1],\$protein)
  10.092 μs (0 allocations: 0 bytes)
4.186781819010877

julia> xprot = coor(protein);

julia> @btime distance(\$ligand[1],\$xprot)
  5.590 μs (0 allocations: 0 bytes)
4.186781819010877

```

"""
@inline distance(x₁,y₁,z₁,x₂,y₂,z₂) = sqrt((x₂-x₁)^2+(y₂-y₁)^2+(z₂-z₁)^2)
@inline distance(x,y,z,at::Atom) = distance(x,y,z,at.x,at.y,at.z)
@inline distance(at::Atom,x,y,z) = distance(x,y,z,at)
@inline distance(x::Atom,y::Atom) = distance(x.x,x.y,x.z,y.x,y.y,y.z)

@inline distance_sq(x₁,y₁,z₁,x₂,y₂,z₂) = (x₂-x₁)^2+(y₂-y₁)^2+(z₂-z₁)^2
@inline distance_sq(x,y,z,at::Atom) = distance_sq(x,y,z,at.x,at.y,at.z)
@inline distance_sq(at::Atom,x,y,z) = distance_sq(x,y,z,at)
@inline distance_sq(x::Atom,y::Atom) = distance_sq(x.x,x.y,x.z,y.x,y.y,y.z)

"""

```
closest(x,y)
```

Computes the minimum distance between two sets of atoms and returns the index(es) of the atoms and their distance. The input may be a vector of `Atom`s, or the matrices that are output of the `coor` function. Using the matrices results in greater speed, particularly if the greatest dimension of the matrices is the first one (`xyz_in_cols=true` option, which is the default option).

### Examples

```julia-repl
julia> model = wget("1BSX");

julia> protein = select(model,"protein");

julia> ligand = select(model,"resname T3");

julia> closest(ligand,protein)
((43, 3684), 2.7775834820937417)

julia> ligand[43]
    4037   O1      T3     B        2      512  -22.568   81.625    3.159 36.59  1.00     1       -      4041

julia> closest(ligand[1],protein)
(728, 4.186781819010877)

julia> using BenchmarkTools

julia> @btime closest(\$ligand[1],\$protein)
  9.940 μs (0 allocations: 0 bytes)
(728, 4.186781819010877)

julia> xprot = coor(protein);

julia> @btime closest(\$ligand[1],\$xprot)
  5.451 μs (0 allocations: 0 bytes)
(728, 4.186781819010877)

```

"""
function closest(x,y,z,atoms::AbstractVector{Atom})
  d = +Inf
  imin = -1
  for (i,at) in pairs(atoms)
    dᵢ = distance_sq(x,y,z,at)
    if dᵢ < d
      d = dᵢ
      imin = i
    end
  end
  imin, sqrt(d)
end
closest(at::Atom, y::AbstractVector{Atom}) = closest(at.x,at.y,at.z,y)
distance(at::Atom, y::AbstractVector{Atom}) = closest(at,y)[2]

# Two vectors of atoms
function closest(x::AbstractVector{Atom}, y::AbstractVector{Atom})
  lx = length(x)
  ly = length(y)
  if lx >= ly
    outer = x
    inner = y
  else
    outer = y
    inner = x
  end
  d = +Inf
  pair = (-1, -1)
  for (i,atx) in pairs(outer)
    j, dᵢⱼ = closest(atx,inner)
    if dᵢⱼ < d
      d = dᵢⱼ
      pair = (i,j)
    end
  end
  if lx >= ly
    pair, d
  else
    (pair[2], pair[1]), d
  end
end
distance(x::AbstractVector{Atom}, y::AbstractVector{Atom}) = closest(x,y)[2]

# Two residues
closest(r1::Residue,r2::Residue) =
  closest(@view(r1.atoms[r1.range]),@view(r2.atoms[r2.range]))
distance(r1::Residue,r2::Residue) =
  distance(@view(r1.atoms[r1.range]),@view(r2.atoms[r2.range]))

#
# The following versions are thought to be faster because the coordinates
# are provided as matrices
#

# Three coordinates and a matrix of coordinates
function closest(x,y,z,A::Matrix; xyz_in_cols=true)
  n, m = size(A)
  imin = 1
  if xyz_in_cols
    @assert m == 3 "Number of columns of coordinates matrix must be 3 when xyz_in_cols=true"
    d = distance_sq(x,y,z,A[1,1],A[1,2],A[1,3])
    for i in 2:n
      dᵢ = distance_sq(x,y,z,A[i,1],A[i,2],A[i,3])
      if dᵢ < d
        d = dᵢ
        imin = i
      end
    end
  else
    @assert n == 3 "Number of rows of coordinates matrix must be 3 when xyz_in_cols=false"
    d = distance_sq(x,y,z,A[1,1],A[2,1],A[3,1])
    for j in 2:m
      dⱼ = distance_sq(x,y,z,A[1,j],A[2,j],A[3,j])
      if dⱼ < d
        d = dⱼ
        imin = j
      end
    end
  end
  imin, sqrt(d)
end
distance(x,y,z,A::Matrix; xyz_in_cols=true) = closest(x,y,z,A, xyz_in_cols=xyz_in_cols)[2]

# One Atom and a matrix of coordinates
closest(at::Atom,A::Matrix; xyz_in_cols=true) =
  closest(at.x,at.y,at.z,A,xyz_in_cols=xyz_in_cols)
distance(at::Atom,A::Matrix; xyz_in_cols=true) =
  distance(at.x,at.y,at.z,A,xyz_in_cols=xyz_in_cols)

#
# Two matrices of coordinates. The smallest matrix is iterated atom by atom,
# and the largest matrix is iterated on the above function, which is loop-vectorized
#
function closest(A::Matrix,B::Matrix; xyz_in_cols=true)
  nA, mA = size(A)
  nB, mB = size(B)
  pair = (1,1)
  if xyz_in_cols
    @assert (mA == 3 && mB == 3) "Matrices must have 3 columns if xyz_in_cols=true"
    d = distance(B[1,1],B[1,2],B[1,3],A[1,1],A[1,2],A[1,3])
    if nA > nB
      for i in 2:nB
        j, dᵢⱼ = closest(B[i,1],B[i,2],B[i,3],A,xyz_in_cols=true)
        if dᵢⱼ < d
          d = dᵢⱼ
          pair = (i,j)
        end
      end
    else
      for i in 2:nA
        j, dᵢⱼ = closest(A[i,1],A[i,2],A[i,3],B,xyz_in_cols=true)
        if dᵢⱼ < d
          d = dᵢⱼ
          pair = (i,j)
        end
      end
    end
  else
    @assert (nA == 3 && nB == 3) "Matrices must have 3 rows if xyz_in_cols=false"
    d = distance(B[1,1],B[2,1],B[3,1],A[1,1],A[2,1],A[3,1])
    if mA > mB
      for j in 2:mB
        i, dᵢⱼ = closest(B[1,j],B[2,j],B[3,j],A,xyz_in_cols=false)
        if dᵢⱼ < d
          d = dᵢⱼ
          pair = (i,j)
        end
      end
    else
      for j in 2:mA
        i, dᵢⱼ = closest(A[1,j],A[2,j],A[3,j],B,xyz_in_cols=false)
        if dᵢⱼ < d
          d = dᵢⱼ
          pair = (i,j)
        end
      end
    end
  end
  pair, d
end
distance(A::Matrix,B::Matrix; xyz_in_cols=true) = closest(A,B,xyz_in_cols=xyz_in_cols)[2]


text:
#=
  From a BLOG model in
  Yi Wu, Lei Li, Stuart Russell, Rastislav Bodik
  "Swift: Compiled Inference for Probabilistic Programming Languages"
  Page 2
  """
  [The model] demonstrates the open-universe urn-ball model. In this version
  the query asks for the color of the next random pick from an urn given the
  colors of balls drawn previously.
  """

  Summary Statistics
  parameters      mean       std   naive_se      mcse        ess       rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64   Float64?   Float64?      Float64?

    numBalls    2.6863    3.6429     0.0364    0.2942   107.7806     1.0000        4.7524
    ...
     drawnD5    1.7703    1.8844     0.0188    0.1436   120.3713     1.0076        5.3076
     colorD5    1.8228    0.3819     0.0038    0.0296   107.6107     0.9999        4.7450

   Distributions of variable numBalls (num:0)
   1.00000 =>    5969  (0.596900)
   2.00000 =>    1409  (0.140900)
   3.00000 =>     843  (0.084300)
   4.00000 =>     489  (0.048900)
   9.00000 =>     319  (0.031900)
   5.00000 =>     242  (0.024200)
   19.00000 =>     209  (0.020900)
   7.00000 =>     161  (0.016100)
   6.00000 =>     109  (0.010900)
   15.00000 =>      72  (0.007200)
   16.00000 =>      54  (0.005400)
   20.00000 =>      46  (0.004600)
   11.00000 =>      42  (0.004200)
   12.00000 =>      36  (0.003600)

   Distributions of variable drawnD5 (num:0)
   1.00000 =>    7070  (0.707000)
   2.00000 =>    1149  (0.114900)
   3.00000 =>     905  (0.090500)
   4.00000 =>     449  (0.044900)
   10.00000 =>     109  (0.010900)
   9.00000 =>      72  (0.007200)
   5.00000 =>      68  (0.006800)
   8.00000 =>      54  (0.005400)
   6.00000 =>      42  (0.004200)
   19.00000 =>      38  (0.003800)
   7.00000 =>      36  (0.003600)
   11.00000 =>       8  (0.000800)

   Distributions of variable colorD5
   blue       =>    8228  (0.822800)
   green      =>    1772  (0.177200)

   I.e. it's quite probable that there is just one ball (or perhaps two).

  Cf ~/blog/urn_ball.blog
     ~/psi/urn_ball.psi
     ~/webppl/urn_ball.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

# Blue: 1 Green: 2
# We observe 4 Green balls
@model function urn_ball(obs=[2,2,2,2])
    n = 5

    blue = 1
    green = 2

    numBalls ~ DiscreteUniform(1,20)

    # PG(5) throws
    # ERROR: LoadError: CTaskException:
    # DimensionMismatch("tried to assign 7 elements to 16 destinations")
    # for this filldist() expression
    # color ~ filldist(Categorical([0.9,0.1]),numBalls) # [Blue,Green]

    # PG(5) don't have any problems with this, though
    color = tzeros(numBalls)
    for i in 1:numBalls
        color[i] ~ Categorical([0.9,0.1]) # [Blue,Green]
    end

    drawn ~ filldist(DiscreteUniform(1,numBalls),n)

    for i in 1:length(obs)
        obs[i] ~ Dirac(color[drawn[i]])
    end

    drawnD5 ~ Dirac(drawn[5])
    colorD5 ~ Dirac(color[drawn[5]])

end

obs=[2,2,2,2] # Observed only green balls
model = urn_ball(obs)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:numBalls)
show_var_dist_pct(chns,:drawnD5)
show_var_dist_pct(chns,:colorD5,["green","blue"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["pa_ringbuffer"], :pa_ringbuffer),
    # FileProduct(prefix, "include/pa_ringbuffer.h", :pa_ringbuffer_h),
]

# Download binaries from hosted location
bin_prefix = "https://s3.us-east-2.amazonaws.com/ringbuffersbuilder"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/pa_ringbuffer.v19.6.0.aarch64-linux-gnu.tar.gz", "8724038e5c3bec6159e5f5816d9b2e4aea14f1006c019dcba3d2fe1f345b7639"),
    Linux(:aarch64, :musl) => ("$bin_prefix/pa_ringbuffer.v19.6.0.aarch64-linux-musl.tar.gz", "61818246bac6925543906319648eb50d8c2deac1380cb892f4ee0ecbd6fc8389"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/pa_ringbuffer.v19.6.0.arm-linux-gnueabihf.tar.gz", "c24a9abc0a07e48d9cc2ebd332479e2affac505e295593d37585339680ad0a37"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/pa_ringbuffer.v19.6.0.arm-linux-musleabihf.tar.gz", "8301c7b87e75bd049a41d72c3ee1f902912fd73315fcc73bbb3705598aa32407"),
    Linux(:i686, :glibc) => ("$bin_prefix/pa_ringbuffer.v19.6.0.i686-linux-gnu.tar.gz", "f00ab4e1727ba60e4305363241ce834e32898448f1eec8c6664f3a4a729886f6"),
    Linux(:i686, :musl) => ("$bin_prefix/pa_ringbuffer.v19.6.0.i686-linux-musl.tar.gz", "a7c5bcac2f24448ced2b2b6ef350fc8b8124be3e834b410733a6f0b67a2f0df3"),
    Windows(:i686) => ("$bin_prefix/pa_ringbuffer.v19.6.0.i686-w64-mingw32.tar.gz", "134a1c5e27ed634183527aefab5e62715ebaf23b462ef006fc72d4e9eaac2b7b"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/pa_ringbuffer.v19.6.0.powerpc64le-linux-gnu.tar.gz", "fb69d38b4221c6163072b86b7a11efd78d91c5279d89e005bc166bc74d805af7"),
    MacOS(:x86_64) => ("$bin_prefix/pa_ringbuffer.v19.6.0.x86_64-apple-darwin14.tar.gz", "d0441ce81f7dedd42e8fef353d946e3e45da8d3b4537d76ad0086d263e14da26"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/pa_ringbuffer.v19.6.0.x86_64-linux-gnu.tar.gz", "7a7dfae0106f31e62f5f4a3de1860c84dfafb6395dd6a507027d9a0c405c42bb"),
    Linux(:x86_64, :musl) => ("$bin_prefix/pa_ringbuffer.v19.6.0.x86_64-linux-musl.tar.gz", "e6de813624b9023a638bce804de5df012093f90e2fecf7b9aad91cecc1a4dc7e"),
    FreeBSD(:x86_64) => ("$bin_prefix/pa_ringbuffer.v19.6.0.x86_64-unknown-freebsd11.1.tar.gz", "a5abb7de8881c59aadc6c3ff028ed360cdb6ad457f432808af4a4934aa11b44e"),
    Windows(:x86_64) => ("$bin_prefix/pa_ringbuffer.v19.6.0.x86_64-w64-mingw32.tar.gz", "2908d3cca37e8af3d307b290f8c0ccafd0c835f61a0f11707af9d90ccc055079"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error(string("$(triplet(platform_key())) is not a supported platform. ",
                 "Please see https://github.com/JuliaAudio/RingBuffersBuilder if ",
                 "you'd like to build your own"))
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
@testset "Show recorder events in EmulationModel" begin
    template = get_thermal_standard_uc_template()
    c_sys5_uc_re = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc_re";
        add_single_time_series=true,
        force_build=true,
    )
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    model = EmulationModel(template, c_sys5_uc_re; optimizer=GLPK_optimizer)

    @test build!(model; executions=10, output_dir=mktempdir(cleanup=true)) ==
          BuildStatus.BUILT
    @test run!(model) == RunStatus.SUCCESSFUL

    recorder_log = joinpath(PSI.get_recorder_dir(model), "execution.log")
    events = list_recorder_events(PSI.ParameterUpdateEvent, recorder_log)
    @test !isempty(events)
    events = list_recorder_events(PSI.InitialConditionUpdateEvent, recorder_log)
    @test !isempty(events)
    for wall_time in (true, false)
        show_recorder_events(
            devnull,
            PSI.InitialConditionUpdateEvent,
            recorder_log,
            wall_time=wall_time,
        )
    end
end


text:
using Profile, ChromeProfileFormat

function foo()
    A = randn(1000, 1000)

    return inv(A)
end

foo()
Profile.clear()
@profile foo()

data, lines = Profile.retrieve()
data2, lines2 = Profile.flatten(data, lines)

ChromeProfileFormat.save_cpuprofile("test.cpuprofile", data2, lines2)


text:
using FactorGraph
using Test

@testset "fgraph" begin
    include("fgraph.jl")
end

@testset "contract" begin
    include("contract.jl")
end


text:
using GoogleDrive
using Test

@testset "GoogleDrive.jl" begin
    # Write your own tests here.
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using SparseArrays
Random.seed!(233)

################## Load Operator ##################
if Sys.islinux()
py"""
import tensorflow as tf
libSparseSolver = tf.load_op_library('build/libSparseSolver.so')
@tf.custom_gradient
def sparse_solver(ii,jj,vv,kk,ff,d):
    u = libSparseSolver.sparse_solver(ii,jj,vv,kk,ff,d)
    def grad(dy):
        return libSparseSolver.sparse_solver_grad(dy, u, ii,jj,vv,kk,ff,d)
    return u, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libSparseSolver = tf.load_op_library('build/libSparseSolver.dylib')
@tf.custom_gradient
def sparse_solver(ii,jj,vv,kk,ff,d):
    u = libSparseSolver.sparse_solver(ii,jj,vv,kk,ff,d)
    def grad(dy):
        return libSparseSolver.sparse_solver_grad(dy, u, ii,jj,vv,kk,ff,d)
    return u, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libSparseSolver = tf.load_op_library('build/libSparseSolver.dll')
@tf.custom_gradient
def sparse_solver(ii,jj,vv,kk,ff,d):
    u = libSparseSolver.sparse_solver(ii,jj,vv,kk,ff,d)
    def grad(dy):
        return libSparseSolver.sparse_solver_grad(dy, u, ii,jj,vv,kk,ff,d)
    return u, grad
"""
end

sparse_solver = py"sparse_solver"
################## End Load Operator ##################

# TODO:
d0 = 30
nv = 100
nf = 100
aug = Array(1:d0)
ii = [rand(1:d0,nv);aug...]
jj = [rand(1:d0,nv);aug...]
vv = [rand(nv);rand(d0)]
A = sparse(ii,jj,vv,d0,d0)
kk = rand(1:d0,nf)
ff = rand(nf)
rhs = (sparse(kk,ones(Int64,nf),ff,d0,1)|>Array)[:,1]
u_ = A\rhs
# @show Array(A), rhs
ii = constant(ii,dtype=Int32)
jj = constant(jj,dtype=Int32)
vv = constant(vv)
kk = constant(kk,dtype=Int32)
ff = constant(ff)
d = constant(d0,dtype=Int32)

# ii = constant([1;2;3;2],dtype=Int32)
# jj = constant([1;2;3;3],dtype=Int32)
# vv = constant([1.0;2.0;2.0;1.0])
# kk = constant([1;2;3],dtype=Int32)
# ff = constant([1.0;2.0;3.0])
# d = constant(3,dtype=Int32)
u = sparse_solver(ii,jj,vv,kk,ff,d)
sess = Session()
init(sess)
@show norm(run(sess, u)-u_)
# error("")
# TODO:

# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(m)
    return sum(tanh(sparse_solver(ii,jj,m,kk,ff,d)))
end

m_ = vv
v_ = rand(nv+d0)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

sess = Session()
@show run(sess, dy_)


text:
# PostGIS
"""
    dbconnect()::Connection
Returns a connection to the `postgis_1` OSS [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_1&db=oss&ns=universe) (i.e., your UVA computing ID).

# Example

```
conn = dbconnect()
```
"""
function dbconnect()
    (@isdefined(db_user) && @isdefined(db_pwd)) ||
        throw(ArgumentError("Run sdad_setup! before making a connection."))
    Connection("""
               host = $db_host
               port = $db_port
               dbname = $dbname
               user = $db_user
               password = $db_pwd
               """)
end
"""
    get_licenses()::Vector{String}
List of SPDX ID for every OSI-approved license. (source: [SPDX](https://github.com/spdx/license-list-data))

# Example
```
licences = get_licenses()
```
"""
function get_licenses()
    conn = dbconnect()
    output = execute(conn, "select spdx from universe.licenses") |>
        (licenses -> getproperty.(licenses, :spdx))
    close(conn)
    output
end
# execute(conn,
#         """
#         create table universe.github_repos(
#           id int not null,
#           slug varchar not null,
#           created_at timestamptz not null,
#           is_archived bool not null,
#           is_fork bool not null,
#           is_mirror bool not null,
#           spdx varchar not null,
#           query char(52) not null,
#           as_of timestamptz not null,
#           primary key (id)
#         )
#         """)
# execute(conn,
#         """
#         create table universe.github_repos_tracker(
#           spdx varchar(36) not null,
#           query char(52) not null,
#           total int not null,
#           primary key (spdx, query)
#         )
#         """)
# execute(conn, "drop table universe.github_commits")
# execute(conn, "truncate table universe.github_commits")
# execute(conn,
#         """
#         create table universe.github_commits(
#           slug text not null,
#           author text,
#           sha char(40) not null,
#           additions int not null,
#           deletions int not null,
#           datetime timestamptz not null,
#           primary key (slug, sha)
#         )
#         """)
# execute(conn,
#         """
#         create table universe.bad_repos(
#           slug text not null,
#           primary key (slug)
#         )
#         """)
"""
    insert_record_repos_by_license!(conn::Connection,
                                    license::AbstractString,
                                    created_at::AbstractString,
                                    as_of::AbstractString,
                                    node)
It inserts each record to the `universe.github_repos` table.
"""
function insert_record_repos_by_license!(conn::Connection,
                                         license::AbstractString,
                                         created_at::AbstractString,
                                         as_of::AbstractString,
                                         node)
    @unpack isPrivate, databaseId, nameWithOwner, createdAt,
            isArchived, isFork, isMirror = node
    node.isPrivate && return
    execute(conn,
            """insert into universe.github_repos values(
               '$databaseId', '$nameWithOwner', '$createdAt',
               $isArchived, $isFork, $isMirror,
               '$license', '$created_at', '$as_of'
               )
               on conflict (id) do update set
               slug = excluded.slug,
               is_archived = excluded.is_archived,
               is_fork = excluded.is_fork,
               is_mirror = excluded.is_mirror,
               spdx = excluded.spdx,
               created_at = excluded.created_at,
               as_of = excluded.as_of
            """)
end
"""
    insert_commit!(conn::Connection, slug::AbstractString, node)
It inserts each record to the `universe.github_commits` table.
"""
function insert_commit!(conn::Connection, slug::AbstractString, node)
    @unpack author, oid, committedDate, additions, deletions = node
    user = author.user
    login = isnothing(user) ? "null" : "'$(user.login)'"
    execute(conn,
            """insert into universe.github_commits values(
               '$slug', $login, '$oid',
               $additions, $deletions, '$committedDate'
               )
               on conflict (slug, sha) do nothing
            """)
end


text:
const RandSeed     = 1050868099
const SplitMixSeed = 0x3ff9e3779b97f4a8

#=
    __SplitMix__ generates 64 pseudorandom bits (UInt64)
	               intended to be used initializing values
				         (internal state) for more robust RNGs
=#

splitmix_state = OneInternalState(SplitMixSeed)

function splitmix(internal::OneInternalState{T}) where {T<:IUI64}
	s1  = state(internal)
	s1 += 0x9e3779b97f4a7c15
	z   = s1
	z  ⊻= z >> 30
	z  ⊻= z >> 27
	z  ⊻= z >> 31

	state(internal, s1)
	return z
end

splitmix() = splitmix(splitmix_state)

function seed(::Type{OneInternalState{T}}, n::Int=RandSeed) where {T<:IUI32}
	srand(n)
	s1 = rand(T)
	while iszero(s1) s1 = rand(T) end
	return OneInternalState{T}(s1)
end

function seed(::Type{TwoInternalStates{T}}, n::Int=RandSeed) where {T<:IUI64}
	srand(n)
    s1, s2 = rand(T, 2)
	while iszero(&(s1,s2)) s1, s2 = rand(T, 2) end
	return TwoInternalStates{T}(s1, s2)
end

function seed(::Type{FourInternalStates{T}}, n::Int=RandSeed) where {T<:IUI64}
	srand(n)
	s1, s2, s3, s4 = rand(T, 4)
	while iszero(&(s1,s2,s3,s4)) s1,s2,s3,s4 = rand(T, 4) end
	return FourInternalStates{T}(s1, s2, s3, s4)
end

function seed(::Type{EightInternalStates{T}}, n::Int=RandSeed) where {T<:IUI64}
	srand(n)
	s1, s2, s3, s4, s5, s6, s7, s8 = rand(UInt64, 8)
	while iszero(&(s1,s2,s3,s4,s5,s6,s7,s8)) s1,s2,s3,s4,s5,s6,s7,s8 = rand(T, 8) end
	return EightInternalStates{T}(s1, s2, s3, s4, s5, s6, s7, s8)
end


text:
@inline max_min(x,y) = ifelse(abs(y) < abs(x), (x,y) , (y,x))

function maxtomin(a::T, b::T, c::T, d::T) where {T}
    a, b = max_min(a, b)
    c, d = max_min(c, d)

    a, c = max_min(a, c)
    b, d = max_min(b, d)

    b, c = max_min(b, c)

    return a, b, c, d
end

function maxtomin(a::T, b::T, c::T, d::T) where {T}
    a, b = max_min(a, b)
    c, d = max_min(c, d)
    a, c = max_min(a, c)
    b, d = max_min(b, d)
    b, c = max_min(b, c)
    return a, b, c, d
end

function vec_sum(x0::T, x1::T, x2::T, x3::T) where {T}
    s3 = x3
    s2, e3 = two_sum(x2, s3)
    s1, e2 = two_sum(x1, s2)
    s0, e1 = two_sum(x0, s1)
    return s0,e1,e2,e3
end

function vecsum_errbranch(x::NTuple{4,T}) where {T}
    y = r = e = zeros(T, 4)
    j = 1
    e[1] = x[1]
    for i = 1:2
        r[i], t = two_sum(e[i], x[i+1])
        if t !== zero(T)
            y[j] = r[i]
            e[i+1] = t
            j += 1
        else
            e[i+1] = r[i]
        end
    end
    y[j], y[j+1] = two_sum(e[3], x[4])
    return y
end

function fast_vecsum_errbranch(x::NTuple{4,T}) where {T}
    y = zeros(T, 4)
    j = 1
    # e[1] = x1
    # i = 1
    r, t = two_sum(x[1], x[2])
    if t !== zero(T)
       y[j] = r
       e = t
       j += 1
    else
       e = r
    end
    # i = 2
    r, t = two_sum(e, x[3])
    if t !== zero(T)
       y[j] = r
       e = t
       j += 1
    else
       e = r
    end

    y[j], y[j+1] = two_sum(e, x[4])
    return y
end

function fast_vecsum_errbranch(x1::T,x2::T,x3::T,x4::T) where {T}
    y = zeros(T, 4)
    j = 1
    # e[1] = x1
    # i = 1
    r, t = two_sum(x1, x2)
    if t !== zero(T)
       y[j] = r
       e = t
       j += 1
    else
       e = r
    end
    # i = 2
    r, t = two_sum(e, x3)
    if t !== zero(T)
       y[j] = r
       e = t
       j += 1
    else
       e = r
    end

    y[j], y[j+1] = two_sum(e, x4)
    return y
end

function quadword(x1::T, x2::T, x3::T, x4::T) where {T}
    a1, a2 = two_sum(x1, x2)
    b1, b2 = two_sum(x3, x4)
    c1, c2 = two_sum(a1, b1)
    d1, d2 = two_sum(a2, b2)
    e1to4 = vec_sum(c1,c2,d1,d2)
    y = vecsum_errbranch(e1to4)
    return (y...,)
end

@inline function fast_quadword(x1::T, x2::T, x3::T, x4::T) where {T}
    a,b,c,d = maxtomin(x1,x2,x3,x4)
    return fast_vecsum_errbranch(a,b,c,d)
end


text:
if Sys.isunix()
    cd(joinpath(dirname(@__FILE__), "src", "fdasrsf"))
    suffix = Sys.isapple() ? "dylib" : "so"
    run(`make clean`)
    run(`make SUFFIX=$suffix`)
    cd(joinpath(dirname(@__FILE__), "src", "gropt"))
    run(`make clean`)
    run(`make SUFFIX=$suffix`)
    cd(joinpath(dirname(@__FILE__), "src", "fdaqmap"))
    run(`make clean`)
    run(`make SUFFIX=$suffix`)
else # Windows
    # these binaries were cross-compiled from Cygwin for x86_64 only using
    # the Makefile_win in the corresponding src directories and the windows bin
    # release of openblas
    url1 = "https://github.com/jdtuck/ElasticFDA.jl/releases/download/v0.4.0/fdaqmap.7z"
    url2 = "https://github.com/jdtuck/ElasticFDA.jl/releases/download/v0.4.0/fdasrsf.7z"
    url3 = "https://github.com/jdtuck/ElasticFDA.jl/releases/download/v0.4.0/gropt.7z"
    try
        run(`curl -LO $url1`)
        run(`curl -LO $url2`)
        run(`curl -LO $url3`)
    catch
        run(`powershell -Command "(new-object net.webclient).DownloadFile(\"$url1\", \"fdaqmap.7z\")"`)
        run(`powershell -Command "(new-object net.webclient).DownloadFile(\"$url2\", \"fdasrsf.7z\")"`)
        run(`powershell -Command "(new-object net.webclient).DownloadFile(\"$url3\", \"gropt.7z\")"`)
    end
    run(`7z x -y fdaqmap.7z`)
    run(`7z x -y fdasrsf.7z`)
    run(`7z x -y gropt.7z`)
end


text:
# This file is a part of UltraNest.jl, licensed under the MIT License (MIT).

using UltraNest
using Test

using Random, StatsBase, Distributions
using HypothesisTests

@testset "test_ultranest" begin
    dist = product_distribution([
        MixtureModel([truncated(Normal(-1, 0.1), -2, 0), truncated(Normal(1, 0.1), 0, 2)], [0.5, 0.5]),
        MixtureModel([truncated(Normal(-2, 0.25), -3, -1), truncated(Normal(2, 0.25), 1, 3)], [0.3, 0.7]),
        MixtureModel([truncated(Normal(-5, 0.25), -6, -4), truncated(Normal(5, 0.25), 4, 6)], [0.2, 0.8]),
    ])

    prior_transform = let dist=dist
        (u::AbstractVector{<:Real}) -> begin
            x0 = minimum.(dist.v)
            Δx = maximum.(dist.v) .- x0
            u .* Δx .+ x0
        end
    end

    prior_transform_vectorized = let trafo = prior_transform
        (U::AbstractMatrix{<:Real}) -> reduce(vcat, (u -> trafo(u)').(eachrow(U)))
    end

    loglikelihood = let dist = dist
        function (x::AbstractVector{<:Real})
            ll = logpdf(dist, x)
            # lofpdf on MixtureModel returns NaN in gaps between distributions, and UltraNest
            # doesn't like -Inf, so return -1E10
            T = promote_type(Float32, typeof(ll))
            isnan(ll) ? T(-1E10) : T(ll)
        end
    end

    loglikelihood_vectorized = let loglikelihood = loglikelihood
        # UltraNest has variate in rows:
        (X::AbstractMatrix{<:Real}) -> loglikelihood.(eachrow(X))
    end

    paramnames = ["a", "b", "c"]
    smplr = ultranest.ReactiveNestedSampler(paramnames, loglikelihood_vectorized, transform = prior_transform_vectorized, vectorized = true)
    result = smplr.run(min_num_live_points = 4000, cluster_num_live_points = 400)
    @test result isa Dict

    X_un = result["samples"]::AbstractMatrix{<:Real}
    @test X_un isa AbstractMatrix{<:Real}

    X_iid = Array(rand(dist, 10^4)')
    pvalues_ad = [pvalue(KSampleADTest(X_un[:,i], X_iid[:,i])) for i in axes(X_un, 2)]
    @test pvalues_ad[1] > 1E-4 && pvalues_ad[2] > 1E-4 && pvalues_ad[3] > 1E-4

    logz = result["logz"]
    logzerr = result["logzerr"]
    logz_expected = -log(prod(maximum.(dist.v) .- minimum.(dist.v)))
    @test isapprox(logz, logz_expected, atol = 10 * logzerr)
end


text:
"""
Constructor for the POMCPSolver.

Use keyword arguments to specify any field values.

For documentation of each field, see the docstring for the POMCPSolver type.
"""
function POMCPSolver(;eps=0.01,
                      max_depth=typemax(Int),
                      c=1.,
                      tree_queries=100,
                      rng=Base.GLOBAL_RNG,
                      node_belief_updater=DefaultReinvigoratorStub(),
                      estimate_value=RolloutEstimator(POMDPToolbox.RandomSolver()),
                      init_V=0.0,
                      init_N=0,
                      num_sparse_actions=0,
                      default_action=ExceptionRethrow())

    return POMCPSolver(eps,
                       max_depth,
                       c,
                       tree_queries,
                       rng,
                       node_belief_updater,
                       estimate_value,
                       init_V,
                       init_N,
                       num_sparse_actions,
                       default_action)
end

"""
Constructor for the POMCPDPWSolver.

Use keyword arguments to specify any field values.

For documentation of each field, see the docstring for POMCPDPWSolver type.
"""
function POMCPDPWSolver(;eps=0.01,
                      max_depth=typemax(Int),
                      c=1.,
                      tree_queries=100,
                      rng=Base.GLOBAL_RNG,
                      node_belief_updater=DefaultReinvigoratorStub(),
                      estimate_value=RolloutEstimator(POMDPToolbox.RandomSolver()),
                      enable_action_pw=true,
                      alpha_observation::Float64=0.5,
                      k_observation::Float64=10.,
                      alpha_action::Float64=0.5,
                      k_action::Float64=10.,
                      init_V=0.0,
                      init_N=0,
                      next_action=RandomActionGenerator(),
                      default_action=ExceptionRethrow())

    return POMCPDPWSolver(eps,
                       max_depth,
                       c,
                       tree_queries,
                       rng,
                       node_belief_updater,
                       estimate_value,
                       enable_action_pw,
                       alpha_observation,
                       k_observation,
                       alpha_action,
                       k_action,
                       init_V,
                       init_N,
                       next_action,
                       default_action)
end


text:
using Test

pair_legal_monads = (0, 1, -1, 4, -4, 500000000, -500000000, 847658787.834, -847658787.834, 2.568368, -2.568368)
pair_illegal_monads = (NaN16, NaN32, NaN64)

pair_legal_dyads = [(x, y) for x = pair_legal_monads, y = pair_legal_monads]
pair_illegal_dyads = [(x, y) for x = pair_illegal_monads, y = pair_legal_monads]

function test_one_pair_constructor(type::Type, args=[0.0])
    secondindex = length(args) < 2 ? 1 : 2
    @eval begin
        testpair = $type($args...)
        @test testpair.x ≈ $args[1] ≈ testpair[1]
        @test testpair.y ≈ $args[$secondindex] ≈ testpair[$secondindex]
        if $args[1] ≈ $args[$secondindex]
            @test testpair[1] ≈ testpair[$secondindex]
        end
    end
end

function test_pair_constructor_cases_with(type, legal_args, illegal_args)
    @testset "with legal input: $arg" for arg ∈ legal_args
        if eltype(type) == Int && !all(isa.(arg, Int))
            @eval begin
                @test_throws InexactError $type($(arg)...)
            end
        else
            test_one_pair_constructor(type, arg)
        end
    end

    @testset "with illegal input: $arg" for arg ∈ illegal_args
        if any(isnan.(arg))
            if eltype(type) == Int
                @eval begin
                    @test_throws InexactError $type($(arg)...)
                end
            else
                @eval begin
                    @test_throws AssertionError $type($(arg)...)
                end
            end
        else
            @eval begin
                @test_throws MethodError $type($(arg)...)
            end
        end
    end

    @testset "with irrational input" begin
        if eltype(type) == Int
            @eval begin
                @test_throws MethodError $type(π)
            end
        else
            test_one_pair_constructor(type, π)
        end
    end
end

function test_all_pair_constructors()
    @testset "Pair constructor tests" begin
        log_test_batch_begin("Pair constructor")
        @testset "Constructor test set with type: $name" for (type, name) ∈ pair_type_iterable
            @testset "Niladic constructor for $name" begin
                test_one_pair_constructor(type)
                increment_progress()
            end #nilad

            @testset "Monadic constructors for $name" begin
                test_pair_constructor_cases_with(type, pair_legal_monads, pair_illegal_monads)
                increment_progress()
            end #monad

            @testset "Dyadic constructors for $name" begin
                test_pair_constructor_cases_with(type, pair_legal_dyads, pair_illegal_dyads)
                increment_progress()
            end
        end
        log_test_batch_end()
    end
end


text:
addprocs(3)
using Primes.factor
ref = @spawn factor(21883298135690819)
factors = fetch(ref)
@assert factors == Dict(234711901=>1,93234719=>1)


text:
## low-level API

const nvtxRangeId_t = UInt64
function start_range(msg::String)
    active[] && ccall((:nvtxRangeStartA, libnvtx), nvtxRangeId_t, (Cstring,), msg)
end
function end_range(id::nvtxRangeId_t)
    active[] && ccall((:nvtxRangeEnd, libnvtx), Cvoid, (nvtxRangeId_t,), id)
end

function push_range(msg::String)
    if active[]
        rv = ccall((:nvtxRangePushA, libnvtx), Cint, (Cstring,), msg)
        if rv < 0
            # FIXME: this fails when not running under `nvprof`
            #warn("Could not push a new range")
        end
    end
end
function pop_range()
    if active[]
        rv = ccall((:nvtxRangePop, libnvtx), Cint, ())
        if rv < 0
            # FIXME: this fails when not running under `nvprof`
            #warn("Could not pop a new range")
        end
    end
end

## high-level API

struct Range
    id::nvtxRangeId_t
end

"""
    range(msg)

Create and start a new range. The range is not automatically stopped, use
[`stop(::Range)`](@ref) for that.

Use this API if you need overlapping ranges, for scope-based use [`@range`](@ref) instead.
"""
function range(msg)
    Range(start_range(msg))
end

stop(r::Range) = end_range(r.id)

"""
    @range "msg" ex

Create a new range and execute `ex`. The range is popped automatically afterwards.

See also: [`range`](@ref)
"""
macro range(msg, ex)
    quote
        push_range($(esc(msg)))
        local ret = $(esc(ex))
        pop_range()
        ret
    end
end


text:
import HierarchicalUtils: NodeType, children, InnerNode, LeafNode, printtree, nodeshow, nodecommshow

# for schema structures
NodeType(::Type{<:Union{Nothing, Entry}}) = LeafNode()  # because sometimes we have empty array extractor
NodeType(::Type{Pair{Symbol,Any}}) = InnerNode()  # because sometimes we have empty array extractor
NodeType(::Type{<:Union{ArrayEntry, DictEntry, MultiEntry}}) = InnerNode()

nodeshow(io::IO, ::Nothing) = print(io, "Nothing")
nodeshow(io::IO, n::Entry) = print(io, "[Scalar - $(join(sort(string.(types(n))), ","))], $(length(keys(n.counts))) unique values")
nodeshow(io::IO, @nospecialize(n::ArrayEntry)) = print(io, "[" * (isnothing(n.items) ? "Empty " : "") * "List]")
nodeshow(io::IO, @nospecialize(n::DictEntry)) = print(io, "[" * (isnothing(n.childs) ? "Empty " : "") * "Dict]")
nodeshow(io::IO, @nospecialize(n::MultiEntry)) = print(io, "[" * (isempty(n.childs) ? "Empty " : "") * "MultiEntry]")

nodecommshow(io::IO, @nospecialize(n::JSONEntry)) = print(io, "\t# updated = $(n.updated)")

children(n::ArrayEntry) = (n.items,)
# using vector of pairs because splatting to named tuple is not good for compiler
children(n::DictEntry) = collect(n.childs)
children(n::MultiEntry) = [Symbol(k) => v for (k,v) in enumerate(n.childs)]

# for extractor structures
# default extractor
NodeType(::Type{<:AbstractExtractor}) = LeafNode()
NodeType(::Type{<:Union{ExtractArray, ExtractDict, MultipleRepresentation, ExtractKeyAsField, AuxiliaryExtractor}}) = InnerNode()

nodeshow(io::IO, ::T) where T <: AbstractExtractor = print(io, "$(nameof(T))")
nodeshow(io::IO, @nospecialize(::ExtractArray)) = print(io, "Array of")
nodeshow(io::IO, @nospecialize(::ExtractDict)) = print(io, "Dict")
nodeshow(io::IO, @nospecialize(n::ExtractCategorical)) = print(io, "Categorical d = $(n.n)")
nodeshow(io::IO, ::ExtractScalar{T}) where {T} = print(io, "$(T)")
nodeshow(io::IO, @nospecialize(::ExtractString)) = print(io, "String")
nodeshow(io::IO, @nospecialize(n::ExtractVector)) = print(io, "FeatureVector with $(n.n) items")
nodeshow(io::IO, @nospecialize(::MultipleRepresentation)) = print(io, "MultiRepresentation")
nodeshow(io::IO, @nospecialize(::ExtractKeyAsField)) = print(io, "KeyAsField")
nodeshow(io::IO, @nospecialize(::AuxiliaryExtractor)) = print(io, "Auxiliary extractor with")

children(n::ExtractArray) = (n.item,)
children(n::Pair) = (n.second,)
children(n::MultipleRepresentation) = n.extractors
children(e::ExtractKeyAsField) = (e.key, e.item)
children(n::AuxiliaryExtractor) = (n.extractor,)
children(n::ExtractDict) = (; Dict(Symbol(k)=>v for (k,v) in n.dict)...)


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function get_par_unb(θ_bound::Array{Float64,1}, MIN::Array{Float64,1}, MAX::Array{Float64,1}, opt_transf::Array{Int64, 1})

# ----------------------------------------------------------------------------------------------------------------------
# Change domain: from bounded to unbounded domain
# - Transformation: 1 natural logarithm, 2 no transformations, 3 generalized logit
# ----------------------------------------------------------------------------------------------------------------------

     i1 = opt_transf .== 1;
     i2 = opt_transf .== 2;
     i3 = opt_transf .== 3;

     θ_unb     = zeros(size(θ_bound));
     θ_unb[i1] = log.(θ_bound[i1] - MIN[i1]);
     θ_unb[i2] = θ_bound[i2];
     θ_unb[i3] = log.((θ_bound[i3] - MIN[i3]) ./ (MAX[i3] - θ_bound[i3]));

     return θ_unb;
end