text:
using TrussMorph
tm = TrussMorph
using Test
using Makie

fp = "/Users/yijiangh/.julia/dev/TrussMorph/test/truss_json/2D_truss_1.json"
t = tm.parse_truss_json(fp)
load_fp = "/Users/yijiangh/.julia/dev/TrussMorph/test/truss_json/2D_truss_load_case.json"

# truss
node_dof = 3 # 2 for truss, 3 for frame
full_node_dof = 3

load = tm.parse_load_json(load_fp, node_dof)

n_v = size(t.X,1)
n_e = size(t.T,1)

# TODO: perturb this and see resulting weight variance
# initial cross secs radius
r = ones(size(t.T,1)) * 0.01

# pertubing the truss randomly within a box
iter = 100

F = tm.assemble_load_vector(load, n_v, node_dof)
perm_vec, perm, n_dof_free = tm.dof_permutation(t.S, n_v, node_dof)

F_perm = perm * F
F_m = Array(F_perm[1:n_dof_free])

weight_fn = tm.get_weight_calculation_fn(t.X, r, t.T, F_m, perm, n_dof_free, node_dof, full_node_dof, t.mp, collect(1:n_v*2))

st_t = time()
for i=1:iter
    box_range = 6
    X = copy(t.X)
    X[2:3,:] += hcat([-1, -1] .+ 2 * rand(Float64,2), [1, 1]*(-box_range/2) .+ box_range * rand(Float64,2))
    # @show X
    weight = weight_fn(reshape(X, 2*n_v))
    @show weight
end
println("avg time: ", (time() - st_t) / iter)

# scene = Scene()
# tm.draw_truss!(scene, t, r, supp_scale=0.2)
# tm.draw_load!(scene, t, load, load_scale=0.05)
# # tm.draw_deformed!(scene, t, U, node_dof)


text:
## Laplacian

function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, L::Laplacian, border::NoPad, inds)
    TT = eltype(out) # accumtype(eltype(out), eltype(A))
    n = 2*length(L.offsets)
    R = CartesianIndices(inds)
    @inbounds for I in R
        tmp = convert(TT, - n * A[I])
        for J in L.offsets
            tmp += A[I+J]
            tmp += A[I-J]
        end
        out[I] = tmp
    end
    out
end

## imgradients
"""
```julia
    imgradients(img, kernelfun=KernelFactors.ando3, border="replicate") -> gimg1, gimg2, ...
```
Estimate the gradient of `img` in the direction of the first and second dimension
at all points of the image, using a kernel specified by `kernelfun`.

# Output

The gradient is returned as a tuple-of-arrays, one for each dimension of the
input; `gimg1` corresponds to the derivative with respect to the first
dimension, `gimg2` to the second, and so on.

# Details

To appreciate the difference between various gradient estimation methods
it is helpful to distinguish between: (1) a continuous scalar-valued
*analogue* image ``f_\\textrm{A}(x_1,x_2)``, where ``x_1,x_2 \\in
\\mathbb{R}``, and (2) its discrete *digital* realization
``f_\\textrm{D}(x_1',x_2')``, where ``x_1',x_2' \\in \\mathbb{N}``, ``1
\\le x_1' \\le M`` and ``1 \\le x_2' \\le N``.

## Analogue image

The gradient of a continuous analogue image ``f_{\\textrm{A}}(x_1,x_2)`` at
location ``(x_1,x_2)`` is defined as the vector
```math
\\nabla \\mathbf{f}_{\\textrm{A}}(x_1,x_2) = \\frac{\\partial
f_{\\textrm{A}}(x_1,x_2)}{\\partial x_1} \\mathbf{e}_{1} +
\\frac{\\partial f_{\\textrm{A}}(x_1,x_2)}{\\partial x_2} \\mathbf{e}_{2},
```
where ``\\mathbf{e}_{d}`` ``(d = 1,2)`` is the unit
vector in the ``x_d``-direction. The gradient points in the direction of
maximum rate of change of ``f_{\\textrm{A}}`` at the coordinates
``(x_1,x_2)``. The gradient can be used to compute the derivative of a
function in an arbitrary direction. In particular, the derivative of
``f_{\\textrm{A}}`` in the direction of a unit vector ``\\mathbf{u}`` is
given by ``\\nabla_{\\mathbf{u}}f_\\textrm{A}(x_1,x_2) = \\nabla
\\mathbf{f}_{\\textrm{A}}(x_1,x_2) \\cdot \\mathbf{u}``, where
``\\cdot`` denotes the dot product.

## Digital image

In practice, we acquire a digital image ``f_\\textrm{D}(x_1',x_2')`` where
the light intensity is known only at a discrete set of locations. This
means that the required partial derivatives are undefined and need to be
approximated using discrete difference formulae [1].

A straightforward way to approximate the partial derivatives is to use
central-difference formulae
```math
 \\frac{\\partial f_{\\textrm{D}}(x_1',x_2')}{\\partial x_1'}  \\approx
        \\frac{f_{\\textrm{D}}(x_1'+1,x_2') - f_{\\textrm{D}}(x_1'-1,x_2') }{2}
```
and
```math
 \\frac{\\partial f_{\\textrm{D}}(x_1',x_2')}{\\partial x_2'}   \\approx
         \\frac{f_{\\textrm{D}}(x_1',x_2'+1) - f_{\\textrm{D}}(x_1',x_2'+1)}{2}.
```
However, the central-difference formulae are very sensitive to noise.
When working with noisy image data,
one can obtain a better approximation of the partial
derivatives by using a suitable weighted combination of the neighboring
image intensities. The weighted combination can be represented as a
*discrete convolution* operation between the image and a
*kernel* which characterizes the requisite weights. In particular,
if ``h_{x_d}`` (``d = 1,2)`` represents a ``2r+1 \\times 2r+1`` kernel, then
```math
 \\frac{\\partial f_{\\textrm{D}}(x_1',x_2')}{\\partial x_d'}  \\approx
\\sum_{i = -r}^r \\sum_{j = -r}^r
f_\\textrm{D}(x_1'-i,x_2'-j)
  h_{x_d}(i,j).
```
The kernel is frequently also called a *mask* or *convolution matrix*.

### Weighting schemes and approximation error

The choice of weights determines the magnitude of the approximation
error and whether the finite-difference scheme is *isotropic*. A
finite-difference scheme is isotropic if the approximation error does
not depend on the orientation of the coordinate system and
*anisotropic* if the approximation error has a directional bias [2].
With a continuous analogue image the magnitude of the gradient would be
invariant upon rotation of the coordinate system, but in practice one
cannot obtain perfect isotropy with a finite set of discrete points.
Hence a finite-difference scheme is typically considered isotropic if
the leading error term in the approximation does not have preferred
directions.

Most finite-difference schemes that are used in image processing are
based on ``3 \\times 3`` kernels, and as noted by [7], many can also be
parametrized by a single parameter ``\\alpha`` as follows:

```math
\\mathbf{H}_{x_{1}} =
\\frac{1}{4 + 2\\alpha}
\\begin{bmatrix}
-1 & -\\alpha & -1 \\\\
0 & 0 & 0 \\\\
 1 & \\alpha & 1
\\end{bmatrix}
\\quad
\\text{and}
\\quad
\\mathbf{H}_{x_{2}} =
\\frac{1}{2 + 4\\alpha}
\\begin{bmatrix}
-1 & 0 & 1 \\\\
-\\alpha & 0 & \\alpha \\\\
 -1 & 0 & 1
\\end{bmatrix},
```
where
```math
\\alpha =
\\begin{cases}
0,  & \\text{Simple Finite Difference}; \\\\
1, &  \\text{Prewitt}; \\\\
2, &  \\text{Sobel}; \\\\
2.4351, &  \\text{Ando}; \\\\
\\frac{10}{3}, &  \\text{Scharr}; \\\\
4, &  \\text{Bickley}.
\\end{cases}
```

## Separable kernel

A kernel is called *separable* if it can be expressed as the
convolution of two one-dimensional filters. With a matrix representation
of the kernel, separability means that the kernel matrix can be written
as an outer product of two vectors. Separable kernels offer
computational advantages since instead of performing a two-dimensional
convolution one can perform a sequence of one-dimensional convolutions.

# Options
You can specify your choice of the finite-difference scheme via the `kernelfun`
parameter. You can also indicate how to deal with the pixels on the border
of the image with the `border` parameter.

## Choices for `kernelfun`
In general `kernelfun` can be any function which satisfies the following
interface:
```julia
    kernelfun(extended::NTuple{N,Bool}, d) -> kern_d,
```
where `kern_d` is the kernel for producing the derivative with respect to
the ``d``th dimension of an ``N``-dimensional array. The parameter `extended[i]` is true
if the image is of size > 1 along dimension ``i``. The parameter `kern_d` may be
provided as a dense or factored kernel, with factored representations
recommended when the kernel is separable.

Some valid `kernelfun` options are described below.

### `KernelFactors.prewitt`

With the *prewit* option [3] the computation of the gradient is based on
the kernels
```math
\\begin{aligned}
\\mathbf{H}_{x_1} & = \\frac{1}{6}
    \\begin{bmatrix}
    -1 & -1 & -1 \\\\
    0 & 0 & 0 \\\\
    1 & 1 & 1
    \\end{bmatrix}
&
\\mathbf{H}_{x_2} & =  \\frac{1}{6}
    \\begin{bmatrix}
    -1 & 0 & 1 \\\\
    -1 & 0 & 1 \\\\
    -1 & 0 & 1
    \\end{bmatrix} \\\\
& = \\frac{1}{6}
    \\begin{bmatrix}
    1 \\\\
    1  \\\\
    1
    \\end{bmatrix}
    \\begin{bmatrix}
    -1 & 0 & 1
    \\end{bmatrix}
&
& = \\frac{1}{6}
    \\begin{bmatrix}
    -1 \\\\
    0  \\\\
    1
    \\end{bmatrix}
    \\begin{bmatrix}
    1 & 1 & 1
    \\end{bmatrix}.
\\end{aligned}
```
See also: [`KernelFactors.prewitt`](@ref) and [`Kernel.prewitt`](@ref)

### `KernelFactors.sobel`

The *sobel* option [4] designates the kernels
```math
\\begin{aligned}
\\mathbf{H}_{x_1} & = \\frac{1}{8}
    \\begin{bmatrix}
    -1 & -2 & -1 \\\\
     0 & 0 & 0 \\\\
     1 & 2 & 1
    \\end{bmatrix}
&
\\mathbf{H}_{x_2} & = \\frac{1}{8}
    \\begin{bmatrix}
    -1 & 0 & 1 \\\\
    -2 & 0 & 2 \\\\
    -1 & 0 & 1
    \\end{bmatrix} \\\\
& = \\frac{1}{8}
    \\begin{bmatrix}
    -1 \\\\
    0  \\\\
    1
    \\end{bmatrix}
    \\begin{bmatrix}
    1 & 2 & 1
    \\end{bmatrix}
&
& = \\frac{1}{8}
    \\begin{bmatrix}
    1 \\\\
    2  \\\\
    1
    \\end{bmatrix}
    \\begin{bmatrix}
    -1 & 0 & 1
    \\end{bmatrix}.
\\end{aligned}
```
See also:  [`KernelFactors.sobel`](@ref) and [`Kernel.sobel`](@ref)

### `KernelFactors.ando3`
The *ando3* option [5] specifies the kernels
```math
\\begin{aligned}
\\mathbf{H}_{x_1} &  =
    \\begin{bmatrix}
    -0.112737 & -0.274526 & -0.112737 \\\\
     0 & 0 & 0 \\\\
     0.112737 & 0.274526 & 0.112737
    \\end{bmatrix}
&
\\mathbf{H}_{x_2}  & =
    \\begin{bmatrix}
    -0.112737 & 0 & 0.112737 \\\\
    -0.274526 & 0 & 0.274526 \\\\
    -0.112737 & 0 & 0.112737
    \\end{bmatrix} \\\\
&  = \\begin{bmatrix}
    -1 \\\\
    0  \\\\
    1
    \\end{bmatrix}
    \\begin{bmatrix}
    0.112737 & 0.274526 & 0.112737
    \\end{bmatrix}
&
&  = \\begin{bmatrix}
    0.112737 \\\\
    0.274526  \\\\
    0.112737
    \\end{bmatrix}
    \\begin{bmatrix}
    -1 & 0 & 1
    \\end{bmatrix}.
\\end{aligned}
```
See also:  [`KernelFactors.ando3`](@ref), and
[`Kernel.ando3`](@ref);  [`KernelFactors.ando4`](@ref), and
[`Kernel.ando4`](@ref); [`KernelFactors.ando5`](@ref), and
[`Kernel.ando5`](@ref)

### `KernelFactors.scharr`

The *scharr* option [6] designates the kernels
```math
\\begin{aligned}
\\mathbf{H}_{x_{1}} & =
\\frac{1}{32}
\\begin{bmatrix}
-3 & -10 & -3 \\\\
0 & 0 & 0 \\\\
 3 & 10 & 3
\\end{bmatrix}
&
\\mathbf{H}_{x_{2}} & =
\\frac{1}{32}
\\begin{bmatrix}
-3 & 0 & 3 \\\\
-10 & 0 & 10\\\\
-3 & 0 & 3
\\end{bmatrix} \\\\
& = \\frac{1}{32}
\\begin{bmatrix}
    -1 \\\\
    0  \\\\
    1
\\end{bmatrix}
\\begin{bmatrix}
    3 & 10 & 3
\\end{bmatrix}
&
& = \\frac{1}{32}
\\begin{bmatrix}
    3 \\\\
    10  \\\\
    3
\\end{bmatrix}
\\begin{bmatrix}
    -1 & 0 & 1
\\end{bmatrix}.
\\end{aligned}
```
See also:  [`KernelFactors.scharr`](@ref) and [`Kernel.scharr`](@ref)

### `KernelFactors.bickley`

The *bickley* option [7,8] designates the kernels
```math
\\begin{aligned}
\\mathbf{H}_{x_1} & = \\frac{1}{12}
    \\begin{bmatrix}
        -1 & -4 & -1 \\\\
         0 & 0 & 0 \\\\
         1 & 4 & 1
    \\end{bmatrix}
&
\\mathbf{H}_{x_2} & = \\frac{1}{12}
    \\begin{bmatrix}
        -1 & 0 & 1 \\\\
        -4 & 0 & 4 \\\\
        -1 & 0 & 1
    \\end{bmatrix} \\\\
& = \\frac{1}{12}
    \\begin{bmatrix}
        -1 \\\\
        0  \\\\
        1
    \\end{bmatrix}
    \\begin{bmatrix}
        1 & 4 & 1
    \\end{bmatrix}
&
&  = \\frac{1}{12}
   \\begin{bmatrix}
        1 \\\\
        4  \\\\
        1
   \\end{bmatrix}
   \\begin{bmatrix}
        -1 & 0 & 1
   \\end{bmatrix}.
\\end{aligned}
```
See also:  [`KernelFactors.bickley`](@ref) and [`Kernel.bickley`](@ref)

## Choices for `border`
At the image edge, `border` is used to specify the padding which will be used
to extrapolate the image beyond its original bounds. As an indicative example
of each option the results of the padding are illustrated on an image consisting of
a row of six pixels which are specified alphabetically: ``\\boxed{a \\, b \\, c \\, d \\, e \\, f}``.
We show the effects of padding only on the left and right border, but analogous
consequences hold for the top and bottom border.

### `"replicate"`

The border pixels extend beyond the image boundaries.
```math
\\boxed{
\\begin{array}{l|c|r}
  a\\, a\\, a\\, a  &  a \\, b \\, c \\, d \\, e \\, f & f \\, f \\, f \\, f
\\end{array}
}
```
See also: [`Pad`](@ref), [`padarray`](@ref), [`Inner`](@ref) and
[`NoPad`](@ref)

### `"circular"`

The border pixels wrap around. For instance, indexing beyond the left border
returns values starting from the right border.
```math
\\boxed{
\\begin{array}{l|c|r}
  c\\, d\\, e\\, f  &  a \\, b \\, c \\, d \\, e \\, f & a \\, b \\, c \\, d
\\end{array}
}
```
See also: [`Pad`](@ref), [`padarray`](@ref), [`Inner`](@ref) and
[`NoPad`](@ref)

### `"symmetric"`
The border pixels reflect relative to a position between pixels. That is, the
border pixel is omitted when mirroring.
```math
\\boxed{
\\begin{array}{l|c|r}
  e\\, d\\, c\\, b  &  a \\, b \\, c \\, d \\, e \\, f & e \\, d \\, c \\, b
\\end{array}
}
```
See also: [`Pad`](@ref), [`padarray`](@ref), [`Inner`](@ref) and
[`NoPad`](@ref)

### `"reflect"`
The border pixels reflect relative to the edge itself.
```math
\\boxed{
\\begin{array}{l|c|r}
  d\\, c\\, b\\, a  &  a \\, b \\, c \\, d \\, e \\, f & f \\, e \\, d \\, c
\\end{array}
}
```
See also: [`Pad`](@ref), [`padarray`](@ref), [`Inner`](@ref) and
[`NoPad`](@ref)

# Example

This example compares the quality of the gradient estimation methods in terms of
the accuracy with which the orientation of the gradient is estimated.

```julia
using Images

values = LinRange(-1,1,128);
w = 1.6*pi;

# Define a function of a sinusoidal grating, f(x,y) = sin( (w*x)^2 + (w*y)^2 ),
# together with its exact partial derivatives.
I = [sin( (w*x)^2 + (w*y)^2 ) for y in values, x in values];
Ix = [2*w*x*cos( (w*x)^2 + (w*y)^2 ) for y in values, x in values];
Iy = [2*w*y*cos( (w*x)^2 + (w*y)^2 ) for y in values, x in values];

# Determine the exact orientation of the gradients.
direction_true = atan.(Iy./Ix);

for kernelfunc in (KernelFactors.prewitt, KernelFactors.sobel,
                   KernelFactors.ando3, KernelFactors.scharr,
                   KernelFactors.bickley)

    # Estimate the gradients and their orientations.
    Gy, Gx = imgradients(I,kernelfunc, "replicate");
    direction_estimated = atan.(Gy./Gx);

    # Determine the mean absolute deviation between the estimated and true
    # orientation. Ignore the values at the border since we expect them to be
    # erroneous.
    error = mean(abs.(direction_true[2:end-1,2:end-1] -
                     direction_estimated[2:end-1,2:end-1]));

    error = round(error, digits=5);
    println("Using \$kernelfunc results in a mean absolute deviation of \$error")
end

# output

Using ImageFiltering.KernelFactors.prewitt results in a mean absolute deviation of 0.01069
Using ImageFiltering.KernelFactors.sobel results in a mean absolute deviation of 0.00522
Using ImageFiltering.KernelFactors.ando3 results in a mean absolute deviation of 0.00365
Using ImageFiltering.KernelFactors.scharr results in a mean absolute deviation of 0.00126
Using ImageFiltering.KernelFactors.bickley results in a mean absolute deviation of 0.00038
```
# References
  1. B. Jahne, *Digital Image Processing* (5th ed.). Springer Publishing Company, Incorporated, 2005. [10.1007/3-540-27563-0](https://doi.org/10.1007/3-540-27563-0)
  2. M. Patra  and  M. Karttunen, "Stencils with isotropic discretization error for differential operators," *Numer. Methods Partial Differential Eq.*, vol. 22, pp. 936–953, 2006. [doi:10.1002/num.20129](https://doi.org/doi:10.1002/num.20129)
  3. J. M. Prewitt, "Object enhancement and extraction," *Picture processing and Psychopictorics*, vol. 10, no. 1, pp. 15–19, 1970.
  4. P.-E. Danielsson and O. Seger, "Generalized and separable sobel operators," in  *Machine Vision for Three-Dimensional Scenes*,  H. Freeman, Ed.  Academic Press, 1990,  pp. 347–379. [doi:10.1016/b978-0-12-266722-0.50016-6](https://doi.org/doi:10.1016/b978-0-12-266722-0.50016-6)
  5. S. Ando, "Consistent gradient operators," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 22, no.3, pp. 252–265, 2000. [doi:10.1109/34.841757](https://doi.org/doi:10.1109/34.841757)
  6. H. Scharr and  J. Weickert, "An anisotropic diffusion algorithm with optimized rotation invariance," *Mustererkennung 2000*, pp. 460–467, 2000. [doi:10.1007/978-3-642-59802-9_58](https://doi.org/doi:10.1007/978-3-642-59802-9_58)
  7. A. Belyaev, "Implicit image differentiation and filtering with applications to image sharpening," *SIAM Journal on Imaging Sciences*, vol. 6, no. 1, pp. 660–679, 2013. [doi:10.1137/12087092x](https://doi.org/doi:10.1137/12087092x)
  8. W. G. Bickley, "Finite difference formulae for the square lattice," *The Quarterly Journal of Mechanics and Applied Mathematics*, vol. 1, no. 1, pp. 35–42, 1948.  [doi:10.1093/qjmam/1.1.35](https://doi.org/doi:10.1093/qjmam/1.1.35)

***

"""
function imgradients(img::AbstractArray, kernelfun::Function, border="replicate")
    extended = map(isextended, axes(img))
    _imgradients(extended, img, kernelfun, extended, border)
end

isextended(ind) = length(ind) > 1

# Add the next dimension to G
function _imgradients(donewhenempty::NTuple{M}, img::AbstractArray{T,N}, kernelfun::Function, extended, border) where {T,M,N}
    d = N-M+1  # the dimension we're working on now
    kern = kernelfun(extended, d)
    return (imfilter(img, kern, border), _imgradients(Base.tail(donewhenempty), img, kernelfun, extended, border)...)
end
# When all N gradients have been calculated, return the result
_imgradients(::Tuple{}, img::AbstractArray, kernelfun::Function, extent, border) = ()

ext: jl


text:
include("xyz_to_json.jl")

directory = joinpath(@__DIR__) * "/S22_3/"
all_files = readdir(directory)
all_files .= directory .* all_files
all_files = filter(x->occursin(".xyz",x), all_files)
#display(all_files)

create_input_rhf.(all_files)


text:
# (a<1) & (a>0)
# (β<1) & (β>0)
# -a+β < 0

using NLopt
using Revise

function func(x, p)
    # simple linear func
    x*p[1] + p[2]
end

function cost(p, xvec, actual)
    # predicted
    pred = func.(xvec, (p,))

    # cost
    sum((pred .- actual).^2)
end

function constraint1(p)
    p[1] - p[2]
end

function main()
    xvec = Float64.(0.0:0.1:10)
    actualp = [6.0, 5.0]
    actualresponse = func.(xvec, (actualp,))

    opt = Opt(:LN_COBYLA, 2)
    opt.lower_bounds = [-Inf, -Inf]
    opt.upper_bounds = [Inf, Inf]
    opt.xtol_rel = 1e-8

    opt.min_objective = (params, grad) -> cost(params, xvec, actualresponse)

    inequality_constraint!(opt, (p, g) -> constraint1(p))

    (minf,minx,ret) = optimize(opt, [1.234, 3.678])
    numevals = opt.numevals # the number of function evaluations
    println("got $minf at $minx after $numevals iterations (returned $ret)")
end

main()

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
"""
    GatesHollowayElbow()

Simple two state elbow potential from Gates and Holloway:
[Journal of Electron Spectroscopy and Related Phenomena, 64/65 (1993) 633-639](https://doi.org/10.1016/0368-2048(93)80131-5)

Has two diabatic states each comprised of the sum of a Morse and a repulsive potential.
The coupling between them is an exponential function of `z` (distance from the surface).
"""
Parameters.@with_kw struct GatesHollowayElbow <: DiabaticModel
    λ₁::Float64 = 3.5
    λ₂::Float64 = 3.5
    z₀::Float64 = 1.4
    x₀::Float64 = 0.6
    α::Float64 = 1.028
    d::Float64 = austrip(4.67u"eV")
    z12::Float64 = 0.5
    c::Float64 = austrip(0.5u"eV")
    γ::Float64 = 1.0
end

NQCModels.nstates(::GatesHollowayElbow) = 2
NQCModels.ndofs(::GatesHollowayElbow) = 1

function NQCModels.potential(model::GatesHollowayElbow, R::AbstractVector)
    Parameters.@unpack λ₁, λ₂, z₀, x₀, α, d, z12, c, γ = model

    repel(x, λ, d) = exp(-λ*(x+d))
    morse(x, d, α) = d*(1-exp(-α*x))^2

    x = R[1]
    z = R[2]

    V11 = morse(x, d, α) + repel(z, λ₁, z₀)
    V22 = morse(z, d, α) + repel(x, λ₂, x₀)
    V12 = c * repel(z, γ,-z12)

    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative!(model::GatesHollowayElbow, D, R::AbstractVector)
    Parameters.@unpack λ₁, λ₂, z₀, x₀, α, d, z12, c, γ = model

    drepel(x, λ, d) = -λ*exp(-λ*(x+d))
    dmorse(x, d, α) = 2α*d*(1-exp(-α*x))*exp(-α*x)

    x = R[1]
    z = R[2]

    # x derivative
    D11 = dmorse(x, d, α)
    D22 = drepel(x, λ₂, x₀)
    D[1] = Hermitian(SMatrix{2,2}(D11, 0, 0, D22))

    # z derivative
    D11 = drepel(z, λ₁, z₀)
    D22 = dmorse(z, d, α)
    D12 = c * drepel(z, γ, -z12)
    D[2] = Hermitian(SMatrix{2,2}(D11, D12, D12, D22))

    return D
end


text:
@testset "SymSparseMatrixCSR" begin
    maxnz=5
    maxrows=5
    maxcols=5
    int_types=(Int32,Int64)
    float_types=(Float32,Float64)
    Bi_types=(0,1)

    for Ti in int_types
        for Tv in float_types
            for Bi in Bi_types
                I = Vector{Ti}()
                J = Vector{Ti}()
                V = Vector{Tv}()
                for (ik, jk, vk) in zip(rand(1:maxrows, maxnz), rand(1:maxcols, maxnz), rand(1:Tv(maxnz), maxnz) )
                    push_coo!(SymSparseMatrixCSR{Bi,Tv,Ti},I,J,V,ik,jk,vk)
                end
                finalize_coo!(SymSparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows, maxcols)
                SYMCSC = Symmetric(sparse(I, J, V, maxrows, maxcols),:U)
                SYMCSR = symsparsecsr(SymSparseMatrixCSR{Bi,Tv,Ti},I, J, V, maxrows, maxcols)

                @test size(SYMCSC)==size(SYMCSR)
                @test SYMCSC == SYMCSR

                @test convert(SymSparseMatrixCSR{Bi}, SYMCSR) === SYMCSR

                @test hasrowmajororder(SYMCSR) == true
                @test hascolmajororder(SYMCSR) == false
                @test getptr(SYMCSR)           == SYMCSR.uppertrian.rowptr
                @test getindices(SYMCSR)       == colvals(SYMCSR)

                @test nnz(SYMCSC.data) == nnz(SYMCSR.uppertrian) <= nnz(SYMCSR)
                @test count(!iszero, SYMCSC.data) == count(!iszero, SYMCSR.uppertrian)

                ICSC,JCSC,VCSC= findnz(SYMCSC.data)
                ICSR,JCSR,VCSR= findnz(SYMCSR)

                @test sort(ICSC)==sort(JCSR) && sort(JCSC)==sort(ICSR) && sort(VCSC)==sort(VCSR)

                v = rand(size(SYMCSC)[2])
                @test SYMCSC*v == SYMCSR*v

                for cBi in Bi_types
                    if Bi == cBi
                        @test convert(SymSparseMatrixCSR{Bi}, SYMCSR) === SYMCSR
                    else
                        SYMCSRC = convert(SymSparseMatrixCSR{cBi}, SYMCSR)
                        @test SYMCSRC == SYMCSR
                        @test SYMCSRC !== SYMCSR
                    end
                end

                for cTi in int_types
                    for cTv in float_types
                        if (Ti,Tv) == (cTi,cTv)
                            @test convert(SymSparseMatrixCSR{Bi,Tv,Ti}, SYMCSR) === SYMCSR
                        else
                            SYMCSRC = convert(SymSparseMatrixCSR{Bi,cTv,cTi}, SYMCSR)
                            @test SYMCSRC == SYMCSR
                            @test SYMCSRC !== SYMCSR
                        end
                    end
                end
            end
        end
    end

end


text:
## Julia program for Polar Simulation
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-20

using Statistics

function polarsim(x::Array)
    while true
        phi1 = rand()*2*pi
        phi2 = rand()*pi - pi/2
        u = rand()
        xdotxi = x[1]*cos(phi1) + x[2]*sin(phi1)*cos(phi2) + x[3]*sin(phi1)*sin(phi2)
        if u <= exp(xdotxi - sum(x.^2)/2)
            return(randn() + xdotxi)
        end
    end
end

# approximate E^pi
function Epi(m, x = [0.1, 1.2, -0.7])
    rhos = ones(m)
    for i = 1:m
        rhos[i] = 1/(2*polarsim(x)^2+3)
    end
    return(mean(rhos))
end

# example
Epi(10)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
@system Context begin
    context ~ ::Nothing
    config ~ ::Config(override)
    clock(config) ~ ::Clock
end

#TODO: remove once Context and Clock get merged
@system DailyContext{Clock => DailyClock}(Context) <: Context

timeunit(c::Context) = timeunit(c.clock)

#HACK: explicitly set up timeunit for default Context
#TODO: merge Context and Clock to remove boilerplates
timeunit(::Type{Context}) = timeunit(Clock)
timeunit(::Type{typefor(Context)}) = timeunit(Context)

#HACK: fallback when timeunit() not available for custom Context
timeunit(C::Type{<:Context}) = only(filter(v -> v.name == :clock, geninfos(C))).type |> scopeof(C).eval |> timeunit

export Context


text:
using Synchrony, Base.Test, CrossDecomposition

# Tests for statistics determined by cross spectrum
angles = [
    2.439564585801219,
    3.0190627944596296,
    1.9558533611578697,
    2.9276497655265747,
    -2.668278147116921,
    2.9670038532462017,
    -2.489899551868197,
    -2.7098277543134612,
    -2.0699115143373175,
    -1.6574845096744018,
    -1.7483177390187856,
    3.0954095279265252,
    2.8604385096072225,
    -3.04237881355147,
    2.6904076043411562
]
r = [
    0.5222318110603368,
    0.29671462136718185,
    0.0766100484932224,
    1.3353190416626122,
    1.2002648555961346,
    1.0636420608553432,
    1.745582800792132,
    0.7810048670915002,
    1.5825305843214879,
    1.3754091313969945,
    1.997253024492295,
    1.239976187334439,
    1.4558221123423296,
    1.6855505422703576,
    0.6713539625790665
]

expangles = exp(im*angles)
expcoef = r.*expangles
true_coh = mean(expcoef)./sqrt(mean(abs2(r)))
true_meanphasediff = mean(expangles)
true_ppc = 1-var(expangles)
true_pli = mean(sign(angles))
true_pli2unbiased = (length(angles)*abs2(mean(sign(angles)))-1)/(length(angles)-1)
true_wpli = abs(mean(imag(expcoef)))./mean(abs(imag(expcoef)))
pairs = Float64[]
for i1 = 1:length(angles), i2 = 1:length(angles)
    i1 != i2 || continue
    push!(pairs, imag(expcoef[i1])*imag(expcoef[i2]))
end
true_wpli2debiased = mean(pairs)./mean(abs(pairs))

oneinput = ones(Complex128, length(angles))
csinput = [oneinput expcoef]

# Single input
for (stat, output) in ((Coherency(), true_coh),
                       (Coherence(), abs(true_coh)),
                       (MeanPhaseDiff(), true_meanphasediff),
                       (PLV(), abs(true_meanphasediff)),
                       (PPC(), true_ppc),
                       (PLI(), true_pli),
                       (PLI2Unbiased(), true_pli2unbiased),
                       (WPLI(), true_wpli),
                       (WPLI2Debiased(), true_wpli2debiased))
    @test_approx_eq computestat(stat, csinput)[1, 2] output
    @test_approx_eq computestat(stat, oneinput, expcoef)[1] output
    @test_approx_eq Base.LinAlg.copytri!(computestat(stat, csinput), 'U', true) computestat(stat, csinput, csinput)
end

# Test correlation
correlation = cor(r, angles)
@test_approx_eq computestat(Correlation(), [r angles])[1, 2] correlation
@test_approx_eq computestat(Correlation(), r, angles)[1] correlation
@test_approx_eq Base.LinAlg.copytri!(computestat(Correlation(), [r angles]), 'U', true) computestat(Correlation(), [r angles], [r angles])

# Test Jammalamadaka circular correlation
angles = [
 -2.386402923969883 0.7120356435187083
 -1.2179485025557881 -0.5107153930236978
 0.28658754714190227 -0.5752672716762053
 -0.9331371629333169 2.1870089874804917
 -0.8995711841426721 5.179928197813174
 0.8887003244657483 0.5249512145071122
 -0.9390694729799546 2.1193361160731787
 2.4997562944878062 4.444408123844822
 -0.15260607000620233 0.32401343526930065
 -1.58821218340456 3.838922093624558
 0.3300644770555941 3.612382182527298
 -0.9159071492390857 1.3286113019566637
 -1.2527393658061554 1.0733970760561626
 -4.364282749017247 1.9318390452682834
 -2.102397073673845 -1.1008030808360818
]
# adiff = randn(1000000)
# angles = [randn(1000000).+adiff randn(1000000).+adiff.+pi/2]

expi = exp(im*angles)

abar, bbar = angle(sum(expi, 1))
aabar = sin(angles[:, 1] .- abar)
bbbar = sin(angles[:, 2] .- bbar)
true_ccor = sum(aabar.*bbbar)/sqrt(sum(abs2(aabar)).*sum(abs2(bbbar)))

@test_approx_eq computestat(JammalamadakaR(), expi)[1, 2] true_ccor
@test_approx_eq computestat(JammalamadakaR(), expi[:, 1], expi[:, 2])[1] true_ccor

# Test Jupp and Mardia circular correlation
am = [real(expi[:, 1]) imag(expi[:, 1])]
bm = [real(expi[:, 2]) imag(expi[:, 2])]
true_ccor = sum(abs2(cor(canoncor(am, bm))))

# println((cor(real(expi[:, 1]), real(expi[:, 2])),
#          cor(real(expi[:, 1]), imag(expi[:, 2])),
#          cor(imag(expi[:, 1]), real(expi[:, 2])),
#          cor(imag(expi[:, 1]), imag(expi[:, 2])),
#          cor(real(expi[:, 1]), imag(expi[:, 1])),
#          cor(real(expi[:, 2]), imag(expi[:, 2]))))
@test_approx_eq computestat(JuppMardiaR(), expi)[1, 2] true_ccor
@test_approx_eq computestat(JuppMardiaR(), expi[:, 1], expi[:, 2])[1] true_ccor

# Test Hurtado modulation index
# Create two signals with some phase-amplitude coupling
function dumbcfc(x, y, nbins)
    angles = angle(x)
    amps = abs(y)
    bins = linspace(-pi, pi, nbins+1)
    mean_amp = [mean(amps[bins[i] .<= angles .<= bins[i+1]]) for i = 1:nbins]
    pj = mean_amp/sum(mean_amp)
    h = -sum(pj.*log(pj))
    hmax = log(nbins)
    (hmax - h)/hmax
end
s1 = exp(im*2pi*[0.05:.1:10;]).*repmat([0.05:.1:1;], 10)
s1s = shuffle(s1)
s2 = exp(im*2pi*[1/18:1/9:11+1/18;]).*repmat([0.99; ones(9)*0.01], 10)

out = computestat(HurtadoModulationIndex(10), [s1 s2 s1s], [s1 s2])
@test_approx_eq out [dumbcfc(s1, s1, 10) dumbcfc(s1, s2, 10)
                     dumbcfc(s2, s1, 10) dumbcfc(s2, s2, 10)
                     dumbcfc(s1s, s1, 10) dumbcfc(s1s, s2, 10)]

# Test UniformScores
d = [50, 120, 192, 210, 220, 250, 262, 291, 292, 320, 321, 340,
     0, 20, 40, 60, 160, 171, 200, 221, 270, 294, 341, 350,
     10, 11, 21, 22, 31, 32, 41, 150, 151, 152, 170, 190, 293,
     30, 70, 110, 172, 180, 191, 240, 251, 260, 261, 290, 351]
groups = [fill(1, 12); fill(2, 12); fill(3, 13); fill(4, 12)]
v = computestat(UniformScores(groups), Complex128[rand(length(d)).*cis(d/180*pi) zeros(length(d))])
@test_approx_eq_eps v[1, 2] 12.81 0.01

csinput = complex(randn(25, 3), randn(25, 3))
csinput2 = complex(randn(25, 3), randn(25, 3))
nmulti = 5
corvec = randn(size(csinput, 1))

# Generate bootstrap weights
bsindices = rand(1:size(csinput, 1), size(csinput, 1), 10)
permindices = Permutation(Coherence(), size(csinput, 1), 10).indices
jcpermindices = Permutation(Coherence(), div(size(csinput, 1), nmulti), 10).indices
groups = Vector{Tuple{Int,Int}}[[(1, 2)], [(3, 2)], [(1, 2), (2, 3)]]

# Single input tests of nested stats
for stat in [Correlation(), Coherence(), Coherency(), MeanPhaseDiff(), PLV(), PPC(), PLI(), PLI2Unbiased(), WPLI(), WPLI2Debiased(), JammalamadakaR(), JuppMardiaR()]
    if isa(stat, Correlation)
        input = abs(csinput)
    else
        input = csinput
    end

    # Test GroupMean
    trueval = computestat(stat, input)
    gmstat = GroupMean(stat, size(input, 2), groups)
    gm = computestat(gmstat, input)
    @test_approx_eq gm [trueval[1, 2], trueval[2, 3], (trueval[1, 2] + trueval[2, 3])/2]

    # Test JackknifeSurrogates
    jnsurrogates = zeros(eltype(trueval), size(input, 1), size(input, 2), size(input, 2))
    for i = 1:size(input, 1)
        jnsurrogates[i, :, :] = computestat(stat, input[[1:i-1; i+1:size(input, 1)], :])
    end
    jnvar = squeeze(var(jnsurrogates, 1, corrected=false), 1)*(size(input, 1)-1)
    jn = computestat(JackknifeSurrogates(stat), input)
    @test_approx_eq jn.trueval trueval
    @test_approx_eq_eps jackknife_bias(jn) (size(input, 1)-1)*(squeeze(mean(jnsurrogates, 1), 1) - trueval) sqrt(eps())
    @test_approx_eq jackknife_var(jn) jnvar

    # Test MultiJackknifeSurrogates
    mjnsurrogates = zeros(eltype(trueval), div(size(input, 1), nmulti), size(input, 2), size(input, 2))
    for i = 1:size(mjnsurrogates, 1)
        mjnsurrogates[i, :, :] = computestat(stat, input[[1:(i-1)*nmulti; i*nmulti+1:size(input, 1)], :])
    end
    mjnvar = squeeze(var(mjnsurrogates, 1, corrected=false), 1)*(size(mjnsurrogates, 1)-1)
    jn = computestat(MultiJackknifeSurrogates(stat, nmulti), input)
    @test_approx_eq jn.trueval trueval
    @test_approx_eq_eps jackknife_bias(jn) (size(mjnsurrogates, 1)-1)*(squeeze(mean(mjnsurrogates, 1), 1) - trueval) sqrt(eps())
    @test_approx_eq jackknife_var(jn) mjnvar

    # Test Jackknife
    a = Jackknife(stat)
    jn = computestat(Jackknife(stat), input)
    @test_approx_eq jn.trueval trueval
    @test_approx_eq jn.var jnvar

    # Test JackknifeCorrelation
    if eltype(trueval) <: Real
        jns = reshape(jnsurrogates, size(input, 1), size(input, 2)*size(input, 2))
        truejc = -reshape(cor(jns, corvec, vardim=1), size(input, 2), size(input, 2))
        jc = computestat(JackknifeCorrelation(stat, corvec), input)
        @test_approx_eq jc truejc
        truejc = -reshape(cor(abs2(jns), corvec, vardim=1), size(input, 2), size(input, 2))
        jc = computestat(JackknifeCorrelation(stat, corvec, Base.Abs2Fun()), input)
        truejc[squeeze(all(jnsurrogates .== jnsurrogates[1, :, :], 1), 1)] = NaN # fixup for inaccurate cor std
        jc[squeeze(all(jnsurrogates .== jnsurrogates[1, :, :], 1), 1)] = NaN # fixup for inaccurate cor std
        @test_approx_eq jc truejc
    else
        @test_throws ArgumentError computestat(JackknifeCorrelation(stat, corvec), input)
    end

    # Test Jackknife{GroupMean}
    jn = computestat(Jackknife(gmstat), input)
    @test_approx_eq jn.trueval [trueval[1, 2], trueval[2, 3], (trueval[1, 2] + trueval[2, 3])/2]
    est2 = mean([jnsurrogates[:, 1, 2] jnsurrogates[:, 2, 3]], 2)
    @test_approx_eq jn.var [jnvar[1, 2], jnvar[2, 3], var(est2, corrected=false)*(size(input, 1)-1)]

    # Test Bootstrap
    estimates = zeros(eltype(trueval), size(input, 2), size(input, 2), size(bsindices, 2))
    for i = 1:size(bsindices, 2)
        estimates[:, :, i] = computestat(stat, input[bsindices[:, i], :])
    end
    bs = computestat(Bootstrap(stat, bsindices), input)
    @test_approx_eq bs estimates

    # Test Bootstrap{GroupMean}
    bs = computestat(Bootstrap(gmstat, bsindices), input)
    @test_approx_eq bs squeeze([estimates[1, 2, :]; estimates[2, 3, :]; (estimates[1, 2, :] + estimates[2, 3, :])/2], 2)

    # Test Permutation
    estimates = zeros(eltype(trueval), size(input, 2), size(input, 2), size(permindices, 2))
    for i = 1:size(permindices, 2)
        estimates[:, :, i] = computestat(stat, input[permindices[:, i], :], input)
    end
    perms = computestat(Permutation(stat, permindices), input)
    @test_approx_eq perms estimates

    # Test Permutation{GroupMean}
    bs = computestat(Permutation(gmstat, permindices), input)
    @test_approx_eq bs squeeze([estimates[1, 2, :]; estimates[2, 3, :]; (estimates[1, 2, :] + estimates[2, 3, :])/2], 2)

    # Test Permutation{JackknifeCorrelation}
    if eltype(trueval) <: Real
        trueperm = zeros(size(input, 2), size(input, 2), size(permindices, 2))
        for i = 1:size(permindices, 2)
            trueperm[:, :, i] = -reshape(cor(reshape(jnsurrogates, size(input, 1), size(input, 2)*size(input, 2))[permindices[:, i], :], corvec, vardim=1), size(input, 2), size(input, 2))
        end
        jcperm = computestat(Permutation(JackknifeCorrelation(stat, corvec), permindices), input)
        @test_approx_eq jcperm trueperm

        trueperm = zeros(size(input, 2), size(input, 2), size(jcpermindices, 2))
        for i = 1:size(jcpermindices, 2)
            trueperm[:, :, i] = -reshape(cor(reshape(mjnsurrogates, size(jcpermindices, 1), size(input, 2)^2)[jcpermindices[:, i], :], corvec[1:size(jcpermindices, 1)], vardim=1), size(input, 2), size(input, 2))
        end
        jcperm = computestat(Permutation(JackknifeCorrelation(stat, nmulti, corvec[1:size(jcpermindices, 1)]), jcpermindices), input)
        @test_approx_eq jcperm trueperm
    else
        @test_throws ArgumentError computestat(Permutation(JackknifeCorrelation(stat, corvec), permindices), input)
    end
end

# Two input tests of nested stats
for stat in [Correlation(), Coherence(), Coherency(), MeanPhaseDiff(), PLV(), PPC(), PLI(), PLI2Unbiased(), WPLI(), WPLI2Debiased(), JammalamadakaR(), JuppMardiaR(), HurtadoModulationIndex(5)]
    if isa(stat, Correlation)
        input = abs(csinput)
        input2 = abs(csinput2)
    else
        input = csinput
        input2 = csinput2
    end

    # Test GroupMean
    trueval = computestat(stat, input, input2)
    gmstat = GroupMean(stat, size(input, 2), groups)
    gm = computestat(gmstat, input, input2)
    @test_approx_eq gm [trueval[1, 2], trueval[2, 3], (trueval[1, 2] + trueval[2, 3])/2]

    # Test JackknifeSurrogates
    jnsurrogates = zeros(eltype(trueval), size(input, 1), size(input, 2), size(input2, 2))
    for i = 1:size(input, 1)
        jnsurrogates[i, :, :] = computestat(stat, input[[1:i-1; i+1:size(input, 1)], :], input2[[1:i-1; i+1:size(input, 1)], :])
    end
    jnvar = squeeze(var(jnsurrogates, 1, corrected=false), 1)*(size(input, 1)-1)
    jn = computestat(JackknifeSurrogates(stat), input, input2)
    @test_approx_eq jn.trueval trueval
    @test_approx_eq_eps jackknife_bias(jn) (size(input, 1)-1)*(squeeze(mean(jnsurrogates, 1), 1) - trueval) sqrt(eps())
    @test_approx_eq jackknife_var(jn) jnvar

    # Test MultiJackknifeSurrogates
    mjnsurrogates = zeros(eltype(trueval), div(size(input, 1), nmulti), size(input, 2), size(input2, 2))
    for i = 1:size(mjnsurrogates, 1)
        mjnsurrogates[i, :, :] = computestat(stat, input[[1:(i-1)*nmulti; i*nmulti+1:size(input, 1)], :], input2[[1:(i-1)*nmulti; i*nmulti+1:size(input, 1)], :])
    end
    mjnvar = squeeze(var(mjnsurrogates, 1, corrected=false), 1)*(size(mjnsurrogates, 1)-1)
    jn = computestat(MultiJackknifeSurrogates(stat, nmulti), input, input2)
    @test_approx_eq jn.trueval trueval
    @test_approx_eq_eps jackknife_bias(jn) (size(mjnsurrogates, 1)-1)*(squeeze(mean(mjnsurrogates, 1), 1) - trueval) sqrt(eps())
    @test_approx_eq jackknife_var(jn) mjnvar

    # Test Jackknife
    jn = computestat(Jackknife(stat), input, input2)
    @test_approx_eq jn.trueval trueval
    @test_approx_eq jn.var jnvar

    # Test JackknifeCorrelation
    if eltype(trueval) <: Real
        truejc = -reshape(cor(reshape(jnsurrogates, size(input, 1), size(input, 2)*size(input2, 2)), corvec, vardim=1), size(input, 2), size(input2, 2))
        jc = computestat(JackknifeCorrelation(stat, corvec), input, input2)
        @test_approx_eq jc truejc
        if !isa(stat, Correlation)
            truejc = -reshape(cor(abs2(reshape(jnsurrogates, size(input, 1), size(input, 2)*size(input2, 2))), corvec, vardim=1), size(input, 2), size(input2, 2))
            jc = computestat(JackknifeCorrelation(stat, corvec, Base.Abs2Fun()), input, input2)
            @test_approx_eq jc truejc
        end
    else
        @test_throws ArgumentError computestat(JackknifeCorrelation(stat, corvec), input, input2)
    end

    # Test Jackknife{GroupMean}
    jn = computestat(Jackknife(gmstat), input, input2)
    @test_approx_eq jn.trueval [trueval[1, 2], trueval[2, 3], (trueval[1, 2] + trueval[2, 3])/2]
    est2 = mean([jnsurrogates[:, 1, 2] jnsurrogates[:, 2, 3]], 2)
    @test_approx_eq jn.var [jnvar[1, 2], jnvar[2, 3], var(est2, corrected=false)*(size(input, 1)-1)]

    # Test Bootstrap
    estimates = zeros(eltype(trueval), size(input, 2), size(input2, 2), size(bsindices, 2))
    for i = 1:size(bsindices, 2)
        estimates[:, :, i] = computestat(stat, input[bsindices[:, i], :], input2[bsindices[:, i], :])
    end
    bs = computestat(Bootstrap(stat, bsindices), input, input2)
    @test_approx_eq bs estimates

    # Test Bootstrap{GroupMean}
    bs = computestat(Bootstrap(gmstat, bsindices), input, input2)
    @test_approx_eq bs squeeze([estimates[1, 2, :]; estimates[2, 3, :]; (estimates[1, 2, :] + estimates[2, 3, :])/2], 2)

    # Test Permutation
    estimates = zeros(eltype(trueval), size(input, 2), size(input2, 2), size(permindices, 2))
    for i = 1:size(permindices, 2)
        estimates[:, :, i] = computestat(stat, input[permindices[:, i], :], input2)
    end
    perms = computestat(Permutation(stat, permindices), input, input2)
    @test_approx_eq perms estimates

    # Test Permutation{GroupMean}
    bs = computestat(Permutation(gmstat, permindices), input, input2)
    @test_approx_eq bs squeeze([estimates[1, 2, :]; estimates[2, 3, :]; (estimates[1, 2, :] + estimates[2, 3, :])/2], 2)

    # Test Permutation{JackknifeCorrelation}
    if eltype(trueval) <: Real
        trueperm = zeros(size(input, 2), size(input, 2), size(permindices, 2))
        for i = 1:size(permindices, 2)
            trueperm[:, :, i] = -reshape(cor(reshape(jnsurrogates, size(input, 1), size(input, 2)*size(input2, 2))[permindices[:, i], :], corvec, vardim=1), size(input, 2), size(input2, 2))
        end
        jcperm = computestat(Permutation(JackknifeCorrelation(stat, corvec), permindices), input, input2)
        @test_approx_eq jcperm trueperm

        trueperm = zeros(size(input, 2), size(input, 2), size(jcpermindices, 2))
        for i = 1:size(jcpermindices, 2)
            trueperm[:, :, i] = -reshape(cor(reshape(mjnsurrogates, size(jcpermindices, 1), size(input, 2)*size(input2, 2))[jcpermindices[:, i], :], corvec[1:size(jcpermindices, 1)], vardim=1), size(input, 2), size(input2, 2))
        end
        jcperm = computestat(Permutation(JackknifeCorrelation(stat, nmulti, corvec[1:size(jcpermindices, 1)]), jcpermindices), input, input2)
        @test_approx_eq jcperm trueperm
    else
        @test_throws ArgumentError computestat(Permutation(JackknifeCorrelation(stat, corvec), permindices), input, input2)
    end
end

# Test nd
offsets = reshape(rand(10)*2pi, 1, 10)
snr = 5*reshape(rand(10), 1, 10)
csinputs = rand(100, 10, 1, 7).*exp(im*(offsets .+ randn(100, 10, 1, 7)./snr))

ngroups = 10
groups = Array(Vector{Tuple{Int, Int}}, ngroups)
for igroup = 1:ngroups
    pairs = Array(Tuple{Int, Int}, rand(1:10))
    for ipair = 1:length(pairs)
        ch1 = rand(1:size(csinputs, 2))
        ch2 = rand(1:size(csinputs, 2)-1)
        ch2 >= ch1 && (ch2 += 1)
        pairs[ipair] = (ch1, ch2)
    end
    groups[igroup] = pairs
end

nbootstraps = 10

for stat in [Correlation, Coherence, Coherency, MeanPhaseDiff, PLV, PPC, PLI, PLI2Unbiased, WPLI, WPLI2Debiased, JammalamadakaR, JuppMardiaR]
    if stat == Correlation
        inputs = abs(csinputs)
    else
        inputs = csinputs
    end

    # True statistic
    tv = mapslices(x->computestat(stat(), x), inputs, (1, 2))
    @test_approx_eq computestat(stat(), inputs) tv
    @test_approx_eq computestat_parallel(stat(), inputs) tv

    # GroupMean
    gmstat = GroupMean(stat(), size(inputs, 2), groups)
    tgm = mapslices(x->computestat(gmstat, x), inputs, (1, 2))
    @test_approx_eq computestat(gmstat, inputs) tgm
    @test_approx_eq computestat_parallel(gmstat, inputs) tgm

    # JackknifeSurrogates
    jns = computestat(JackknifeSurrogates(stat()), inputs)
    tjnvar = mapslices(x->jackknife_var(computestat(JackknifeSurrogates(stat()), x)), inputs, (1, 2))
    @test_approx_eq jns.trueval tv
    @test_approx_eq jackknife_var(jns) tjnvar
    jns = computestat_parallel(JackknifeSurrogates(stat()), inputs)
    @test_approx_eq jns.trueval tv
    @test_approx_eq jackknife_var(jns) tjnvar

    # Jackknife
    jn = computestat(Jackknife(stat()), inputs)
    @test_approx_eq jn.trueval tv
    @test_approx_eq jn.var tjnvar
    jns = computestat_parallel(Jackknife(stat()), inputs)
    @test_approx_eq jn.trueval tv
    @test_approx_eq jn.var tjnvar

    # Jackknife{GroupMean}
    jn = computestat(Jackknife(gmstat), inputs)
    @test_approx_eq jn.trueval mapslices(x->computestat(Jackknife(gmstat), x).trueval, inputs, (1, 2))
    @test_approx_eq jn.var mapslices(x->computestat(Jackknife(gmstat), x).var, inputs, (1, 2))

    # Bootstrap and Permutation
    for t in (Bootstrap, Permutation)
        bs = t(stat(), size(inputs, 1), nbootstraps)
        bstrue = zeros(eltype(tv), size(inputs, 2), size(inputs, 2), nbootstraps, size(inputs)[3:end]...)
        for i = 1:Base.trailingsize(inputs, 3)
            bstrue[:, :, :, i] = computestat(bs, inputs[:, :, i])
        end
        @test_approx_eq computestat(bs, inputs) bstrue
        @test_approx_eq computestat_parallel(bs, inputs) bstrue

        bs = t(gmstat, size(inputs, 1), nbootstraps)
        bstrue = zeros(eltype(tv), length(groups), nbootstraps, size(inputs)[3:end]...)
        for i = 1:Base.trailingsize(inputs, 3)
            bstrue[:, :, i] = computestat(bs, inputs[:, :, i])
        end
        @test_approx_eq computestat(bs, inputs) bstrue
        @test_approx_eq computestat_parallel(bs, inputs) bstrue
    end
end

# Test shift predictor
# x = 0:63
# signal = zeros(length(x), 1, 35)
# for i = 1:size(signal, 3)
#     signal[:, 1, i] = cospi(0.2*x+rand()*2)
# end
# c = multitaper([signal signal], Coherence())
# @test_approx_eq c ones(size(c))
# for lag = 1:5
#     t = multitaper([signal signal[:, :, circshift([1:size(signal, 3)], lag)]], Coherence())
#     sp = multitaper([signal signal], ShiftPredictor(Coherence(), lag))
#     @test_approx_eq t sp
# end

# # Test jackknifed shift predictor
# for stat in (Coherence, PLV), lag = 1:5
#     t = multitaper([signal signal[:, :, circshift([1:size(signal, 3)], lag)]], JackknifeSurrogates(stat()))
#     sp = multitaper([signal signal], JackknifeSurrogates(ShiftPredictor(stat(), lag)))
#     @test_approx_eq t[1] sp[1]
#     (tbias, tvar) = jackknife_bias_var(t...)
#     (spbias, spvar) = jackknife_bias_var(sp...)
#     @test_approx_eq tbias spbias
#     @test_approx_eq tvar spvar
# end


text:
export julia_expression

# to make a variable name from operators
varname(A::BaseOperator) = string(A.name,BaseOpType_expr[Int(A.t)])
varname(A::BaseOpProduct) = Symbol(varname.(A.v)...)
varname(A::Union{ExpVal,Corr}) = varname(A.ops)

# for sum indices use a valid and short symbol that cannot be expressed as a single char, i.e., no risk of collisions with "real" index names
# (the "standard" symbol '#' is not a variable name, so it becomes var"#")
indexpr(ind::QuIndex) = issumindex(ind) ? Symbol(:s̄,subscript(ind.num)) : (isintindex(ind) ? ind.num : Symbol(ind))
indexpr(A) = indexpr.(indices(A))

# the varnames optional (and ignored) argument allows to define a closure getexpr(A) = julia_expression(A,my_names)
function julia_expression(A::Param,varnames=nothing)
    # the parsing here allows "parameters" that are functions
    # e.g., with names such as "f(t)"
    x = Meta.parse(string(A.name))
    # is it a function call?
    if x isa Expr
        @assert isempty(indices(A))
    else
        x = :( $x[$(indexpr(A)...)] )
    end
    A.state == 'c' ? :( conj($x) ) : x
end
julia_expression(A::δ,varnames=nothing) = :( I[$(indexpr(A)...)] )
function julia_expression(A::Union{ExpVal,Corr},varnames=nothing)
    A = unalias(A)
    x = varname(A)
    if varnames === nothing || x ∈ varnames
        :( $x[$(indexpr(A)...)] )
    else
        # <A> = <A^†>^*
        Ap = A'
        x = varname(Ap)
        x ∈ varnames || throw(ValueError("passed allowed varnames = $(varnames), but neither $(varname(A)) nor $x are valid."))
        # since conjugation changes the order, we have to use the indices of A after conjugation
        :( conj($x[$(indexpr(Ap)...)]) )
    end
end
num_expr(s::Number) = :( $(isreal(s) ? Float64(s) : ComplexF64(s)) )
function julia_expression(A::QuTerm,varnames=nothing,s=1)
    isempty(A) && return num_expr(s)
    isempty(A.bares) || throw(ArgumentError("Cannot convert term $A with bare operators to a julia expression."))
    sexpr = isone(s) ? [] : num_expr(s)
    exprs = [sexpr; julia_expression.(A.δs); julia_expression.(A.params);
             julia_expression.(A.expvals,(varnames,)); julia_expression.(A.corrs,(varnames,))]
    length(exprs)==1 ? exprs[1] : :( *($(exprs...)) )
end
function julia_expression(A::QuExpr,varnames=nothing)
    isempty(A) && return 0
    exprs = [julia_expression(term,varnames,s) for (term,s) in sort!(Tuple.(collect(A.terms)))]
    length(exprs)==1 ? exprs[1] : :( +($(exprs...)) )
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_transformations_basic" begin
     err=1e-11
     A = SMatrix{4,4,Float64}([1.0 2.0 3.0 4.0; 5.0 6.0 7.0 8.0; 9.0 9.0 8.0 7.0; 6.0 5.0 4.0 1.0])
     invA = SMatrix{4,4,Float64}([-3.75  2.75   -1.0 0.0; 4.375  -3.875 2.0  -0.5; 0.5    0.5    -1.0 1.0; -1.375 0.875   0.0 -0.5])
     B = SMatrix{4,4,Float64}([3.0 5.0 2.0 4.0; 4.0 1.0 0.0 5.0; 6.0 3.0 2.0 0.0; 1.0 4.0 2.0 1.0])
     invB = SMatrix{4,4,Float64}([0.4 -0.2 0.2 -0.6; 2.9 -1.7 0.2 -3.1; -5.55 3.15 -0.4 6.45; -0.9 0.7 -0.2 1.1])
     C = SMatrix{4,4,Float64}([1.0 2.0 3.0 4.0; 5.0 6.0 7.0 8.0; 9.0 9.0 8.0 7.0; 0.0 0.0 0.0 1.0])
     invC = SMatrix{4,4,Float64}([-3.75 2.75 -1 0; 5.75 -4.75 2.0 1.0; -2.25 2.25 -1.0 -2.0; 0.0 0.0 0.0 1.0])
     D = SMatrix{4,4,Float64}([33.0 32.0 16.0 18.0; 89.0 84.0 40.0 58.0; 118.0 106.0 48.0 88.0; 63.0 51.0 22.0 50.0])
     invD = SMatrix{4,4,Float64}([-1.45 1.45 -1.0 0.6; -13.95 11.95 -6.5 2.6; 25.525 -22.025 12.25 -5.2; 4.825 -4.325 2.5 -1.1])

     m1 = Transformation(A, invA)
     m2 = Transformation(B, invB)
     m = Transformation(C, invC)
     exp = Transformation(D, invD)

     exp_v = Vec(14.0, 38.0, 51.0+err)
     exp_p = Point(18.0, 46.0, 58.0-2*err)
     exp_n = Normal(-8.75, 7.75+6*err, -3.0)

     # is_consistent for manual definition of Transformation
     @test Raytracing.is_consistent(m1)
     @test Raytracing.is_consistent(m2)
     @test Raytracing.is_consistent(m)
     @test Raytracing.is_consistent(exp)

     # approx for multiplications with Transformation
     @test exp ≈ (m1 * m2)
     @test exp_v ≈ (m * Vec(1.0, 2.0, 3.0))
     @test exp_p ≈ (m * Point(1.0, 2.0, 3.0))
     @test exp_n ≈ (m * Normal(3.0, 2.0, 4.0))

end

@testset "test_transformations_rotations" begin
     err = 1e-11
     @test Raytracing.is_consistent(rotation_x(0.1))
     @test Raytracing.is_consistent(rotation_y(0.1))
     @test Raytracing.is_consistent(rotation_z(0.1))
     @test (rotation_x(π/2) * Vec(0.0, 1.0, 0.0+3*err)) ≈ (Vec(0.0, 0.0, 1.0))
     @test (rotation_y(π/2) * Vec(0.0, 0.0-2*err, 1.0)) ≈ (Vec(1.0, 0.0, 0.0))
     @test (rotation_z(π/2) * Vec(1.0+err, 0.0, 0.0)) ≈ (Vec(0.0, 1.0, 0.0))
end

@testset "test_transformations_scaling" begin
     err = 1e-11
     tr1 = scaling(Vec(2.0, 5.0, 10.0+err))
     tr2 = scaling(Vec(3.0, 2.0, 4.0))
     exp = scaling(Vec(6.0, 10.0, 40.0))

     @test Raytracing.is_consistent(tr1)
     @test Raytracing.is_consistent(tr2)
     @test exp ≈ (tr1 * tr2)
end

@testset "test_transformations_translation" begin
     err=1e-11
     tr1 = translation(Vec(1.0, 2.0, 3.0))
     tr2 = translation(Vec(4.0, 6.0, 8.0))
     prd = tr1 * tr2
     exp = translation(Vec(5.0, 8.0, 11.0-7*err))

     @test Raytracing.is_consistent(tr1)
     @test Raytracing.is_consistent(tr2)
     @test Raytracing.is_consistent(prd)
     @test prd ≈ exp
end

@testset "test_transformations_inverse" begin
     A = SMatrix{4,4,Float64}([1.0 2.0 3.0 4.0 ; 5.0 6.0 7.0 8.0 ; 9.0 9.0 8.0 7.0 ; 6.0 5.0 4.0 1.0])
     invA = SMatrix{4,4,Float64}([-3.75 2.75 -1 0 ; 4.375 -3.875 2.0 -0.5 ; 0.5 0.5 -1.0 1.0 ; -1.375 0.875 0.0 -0.5])
     m1 = Transformation(A, invA)
     m2 = inverse(m1)

     @test Raytracing.is_consistent(m2)
     @test Raytracing.is_consistent(m1*m2)
     @test m1*m2 ≈ Transformation()

end


text:
abstract type Service <: Component end

include("products/reserves.jl")
include("products/transfers.jl")


text:
view_x(d::AbstractDataContainer, indices) = view(d.X, indices)
view_y(::AbstractLikelihood, y::AbstractVector, i::AbstractVector) = view(y, i)
function view_y(l::AbstractLikelihood, d::AbstractDataContainer, i::AbstractVector)
    return view_y(l, output(d), i)
end
function view_y(l::AbstractLikelihood, d::MODataContainer, i::AbstractVector)
    return view_y.(l, output(d), Ref(i))
end
function view_y(
    l::AbstractVector{<:AbstractLikelihood}, d::MODataContainer, i::AbstractVector
)
    return view_y.(l, output(d), Ref(i))
end

# Verify that the data is self-consistent and consistent with the likelihood ##
function check_data!(y::AbstractArray, likelihood::Union{Distribution,AbstractLikelihood})
    y, nLatent, likelihood = treat_labels!(y, likelihood)
    return y, nLatent, likelihood
end

# Transform Z as an OptimIP if it's not the case already
# function init_Z(Z::AbstractInducingPoints, Zoptimiser)
#     if Zoptimiser isa Bool
#         Zoptimiser = Zoptimiser ? ADAM(1e-3) : nothing
#     end
#     if Z isa OptimIP
#         return Z
#     else
#         return OptimIP(Z, Zoptimiser)
#     end
# end


text:
using IntervalSets

using Random
using ForwardDiff
using Statistics

export ParametricCurve
function ParametricCurve(parametrization::Function, D::Interval, r::Real; mesh::Int=50)
    n = mesh
    a, b = endpoints(D)
    𝒑 = parametrization
    ts = range(a,b,length=n+1)
    pts = csgUnion(Sphere(𝒑(t),r) for t in ts)
    segts = csgUnion(Cylinder(𝒑(ts[i]),𝒑(ts[i+1]),r) for i in 1:n)
    return csgUnion(pts,segts)
end

export ParametricSurface
struct ParametricSurface <: PrimitiveObject
    parametrization :: Function
    D1 :: Interval{:closed,:closed,Float64}
    D2 :: Interval{:closed,:closed,Float64}
    mesh :: Tuple{Int,Int}
    smooth :: Bool
    function ParametricSurface(parametrization::Function, D1::Interval, D2::Interval; mesh::Tuple{Int,Int}=(10,10), smooth=true)
        new(parametrization, ClosedInterval(D1), ClosedInterval(D2), mesh, smooth)
    end
end

export ColoredParametricSurface
struct ColoredParametricSurface <: PrimitiveObject
    parametrization :: Function
    D1 :: Interval{:closed,:closed,Float64}
    D2 :: Interval{:closed,:closed,Float64}
    mesh :: Tuple{Int,Int}
    smooth :: Bool
    rgb :: Function
    ft :: Function
    smooth_texture :: Bool
    function ColoredParametricSurface(parametrization::Function, D1::Interval, D2::Interval; mesh::Tuple{Int,Int}=(10,10), smooth=true, rgb=(t->RGB(1,0,0)), ft=(t->FT(0,0)), smooth_texture=false)
        new(parametrization, ClosedInterval(D1), ClosedInterval(D2), mesh, smooth, rgb, ft, smooth_texture)
    end
end

function povray_script(parametricsurface::ParametricSurface; preindent=0)
    𝒑 = parametricsurface.parametrization
    D1, D2 = parametricsurface.D1, parametricsurface.D2
    mesh = parametricsurface.mesh
    smooth = parametricsurface.smooth
    return mesh2(𝒑, D1, D2; mesh=mesh, smooth=smooth, preindent=preindent)
end

function povray_script(parametricsurface::ColoredParametricSurface; preindent=0)
    𝒑 = parametricsurface.parametrization
    D1, D2 = parametricsurface.D1, parametricsurface.D2
    mesh = parametricsurface.mesh
    smooth = parametricsurface.smooth
    rgb = parametricsurface.rgb
    ft = parametricsurface.ft
    smooth_texture = parametricsurface.smooth_texture
    return mesh2_with_texture(𝒑, D1, D2; mesh=mesh, smooth=smooth, rgb=rgb, ft=ft, smooth_texture=smooth_texture, preindent=preindent)
end

function mesh2(𝒑, D1, D2; mesh=(10,10), smooth=true, preindent=0)
    n1, n2 = mesh
    D1₋, D1₊ = endpoints(D1)
    D2₋, D2₊ = endpoints(D2)
    Δ1, Δ2 = Δ = (D1₊-D1₋)/n1, (D2₊-D2₋)/n2

    ts = [[t1, t2] for t1 in range(D1₋, D1₊, length = n1 +1), t2 in range(D2₋, D2₊, length = n2 +1)]
    tc = [mean([ts[i1,i2], ts[i1+1,i2], ts[i1,i2+1], ts[i1+1,i2+1]]) for i1 in 1:n1, i2 in 1:n2]

    𝒑′(t)=ForwardDiff.jacobian(𝒑,t)
    𝒆(t) = normalize(cross(𝒑′(t)[1:3,1],𝒑′(t)[1:3,2]))

    𝒑s = 𝒑.(ts)
    𝒑c = 𝒑.(tc)

    𝒆s = 𝒆.(ts)
    𝒆c = 𝒆.(tc)

    Ns(i1, i2) = i1 + (n1+1) * (i2-1)
    Nc(i1, i2) = i1 + n1 * (i2-1)

    F1 = [[Ns(i1,i2)-1, Ns(i1+1,i2)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]
    F2 = [[Ns(i1,i2)-1, Ns(i1,i2+1)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]
    F3 = [[Ns(i1+1,i2+1)-1, Ns(i1+1,i2)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]
    F4 = [[Ns(i1+1,i2+1)-1, Ns(i1,i2+1)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]

    np = (n1+1)*(n2+1) + n1*n2
    nf = 4*n1*n2

    script = ""
    script *= "mesh2{\n" * "  "^(preindent)
    script *= "  vertex_vectors{\n" * "  "^(preindent)
    script *= "    " * povray_script(np) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([𝒑s...]) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([𝒑c...]) * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)

    if smooth
        script *= "  normal_vectors{\n" * "  "^(preindent)
        script *= "    " * povray_script(np) * ", \n" * "  "^(preindent)
        script *= "    " * povray_script([𝒆s...]) * ", \n" * "  "^(preindent)
        script *= "    " * povray_script([𝒆c...]) * "\n" * "  "^(preindent)
        script *= "  }\n" * "  "^(preindent)
    end

    script *= "  face_indices{\n" * "  "^(preindent)
    script *= "    " * povray_script(nf) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([F1...]) * "\n" * "  "^(preindent)
    script *= "    " * povray_script([F2...]) * "\n" * "  "^(preindent)
    script *= "    " * povray_script([F3...]) * "\n" * "  "^(preindent)
    script *= "    " * povray_script([F4...]) * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)
    script *= "}"

    return script
end

function mesh2_with_texture(𝒑, D1, D2; mesh=(10,10), smooth=true, smooth_texture=true, rgb=(t->RGB(1,0,0)), ft=(t->FT(0,0)), preindent=0)
    n1, n2 = mesh
    D1₋, D1₊ = endpoints(D1)
    D2₋, D2₊ = endpoints(D2)
    Δ1, Δ2 = Δ = (D1₊-D1₋)/n1, (D2₊-D2₋)/n2

    ts = [[t1, t2] for t1 in range(D1₋, D1₊, length = n1 +1), t2 in range(D2₋, D2₊, length = n2 +1)]
    tc = [mean([ts[i1,i2], ts[i1+1,i2], ts[i1,i2+1], ts[i1+1,i2+1]]) for i1 in 1:n1, i2 in 1:n2]

    𝒑s = 𝒑.(ts)
    𝒑c = 𝒑.(tc)

    if smooth
        𝒑′(t)=ForwardDiff.jacobian(𝒑,t)
        𝒆(t) = normalize(cross(𝒑′(t)[1:3,1],𝒑′(t)[1:3,2]))

        𝒆s = 𝒆.(ts)
        𝒆c = 𝒆.(tc)
    end

    Cs = rgb.(ts)
    Cc = rgb.(tc)
    Ts = ft.(ts)
    Tc = ft.(tc)

    Ns(i1, i2) = i1 + (n1+1) * (i2-1)
    Nc(i1, i2) = i1 + n1 * (i2-1)

    F1 = [[Ns(i1,i2)-1, Ns(i1+1,i2)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]
    F2 = [[Ns(i1,i2)-1, Ns(i1,i2+1)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]
    F3 = [[Ns(i1+1,i2+1)-1, Ns(i1+1,i2)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]
    F4 = [[Ns(i1+1,i2+1)-1, Ns(i1,i2+1)-1, (n1+1)*(n2+1)+Nc(i1,i2)-1] for i1 in 1:n1, i2 in 1:n2]

    np = (n1+1)*(n2+1) + n1*n2
    nf = 4*n1*n2

    script = ""
    script *= "mesh2{\n" * "  "^(preindent)
    script *= "  vertex_vectors{\n" * "  "^(preindent)
    script *= "    " * povray_script(np) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([𝒑s...]) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([𝒑c...]) * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)

    if smooth
        script *= "  normal_vectors{\n" * "  "^(preindent)
        script *= "    " * povray_script(np) * ", \n" * "  "^(preindent)
        script *= "    " * povray_script([𝒆s...]) * ", \n" * "  "^(preindent)
        script *= "    " * povray_script([𝒆c...]) * "\n" * "  "^(preindent)
        script *= "  }\n" * "  "^(preindent)
    end

    script *= "  texture_list{\n" * "  "^(preindent)
    script *= "    " * povray_script(np) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([Cs...],[Ts...]) * ", \n" * "  "^(preindent)
    script *= "    " * povray_script([Cc...],[Tc...]) * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)

    F1_ = [povray_script(f)*","*string(f[end]) for f in F1]
    F2_ = [povray_script(f)*","*string(f[end]) for f in F2]
    F3_ = [povray_script(f)*","*string(f[end]) for f in F3]
    F4_ = [povray_script(f)*","*string(f[end]) for f in F4]

    script *= "  face_indices{\n" * "  "^(preindent)
    script *= "    " * povray_script(nf) * ", \n" * "  "^(preindent)
    script *= "    " * join(F1_,",") * "\n" * "  "^(preindent)
    script *= "    " * join(F2_,",") * "\n" * "  "^(preindent)
    script *= "    " * join(F3_,",") * "\n" * "  "^(preindent)
    script *= "    " * join(F4_,",") * "\n" * "  "^(preindent)
    script *= "  }\n" * "  "^(preindent)
    script *= "}"

    return script
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
name = "libLLVM"
version = v"10.0.0"

# Include common tools
include("../common.jl")
build_tarballs(ARGS, configure_extraction(ARGS, version, name)...; skip_audit=true)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

include("learning/tasks.jl")
include("learning/examples.jl")
include("learning/models.jl")
include("learning/traits.jl")


text:
type RingBuffer
  buffer::Vector{Float64}
  readInd::Int
  writeInd::Int
end

function enque!(buf::RingBuffer,val::Float64)
  len = length(buf.buffer)
  writeInd = buf.writeInd%len + 1
  if writeInd == buf.readInd
    return buf # silently fail
  end
  @inbounds buf.buffer[buf.writeInd] = val
  buf.writeInd = writeInd
  return buf
end

function deque!(buf::RingBuffer)
  if buf.readInd == buf.writeInd
    return NaN # fail
  end
  @inbounds val::Float64 = buf.buffer[buf.readInd]
  len = length(buf.buffer)
  readInd = buf.readInd%len + 1
  buf.readInd = readInd
  return val
end

function used(buf::RingBuffer)
  if buf.writeInd >= buf.readInd
    return buf.writeInd - buf.readInd
  else
    return length(buf.buffer) - (buf.readInd-buf.writeInd)
  end
end


text:
function isinplace_jump(p,rj)
    if p isa DiscreteProblem && p.f === DiffEqBase.DISCRETE_INPLACE_DEFAULT && rj !== nothing
        # Just a default discrete problem f, so don't use it for iip
        DiffEqBase.isinplace(rj)
    else
        DiffEqBase.isinplace(p)
    end
end

"""
$(TYPEDEF)

Defines a collection of jump processes to associate with another problem type.
- [Documentation Page](https://diffeq.sciml.ai/stable/types/jump_types/)
- [Tutorial
  Page](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/)
- [FAQ
  Page](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/#FAQ)

### Constructors

`JumpProblem`s can be constructed by first building another problem type to
which the jumps will be associated. For example, to  simulate a collection of
jump processes for which the transition rates do not explicitly depend on time
(called `ConstantRateJump`s or `MassActionJump`s), we must first construct a
`DiscreteProblem`
```julia
prob = DiscreteProblem(u0, p, tspan)
```
where `u0` is the initial condition, `p` the parameters and `tspan` the time
span. If we wanted to have the jumps coupled with a system of ODEs, or have
transition rates with explicit time dependence, we would use an `ODEProblem`
instead that defines the ODE portion of the dynamics.

Given `prob` we define the jumps via
- `JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet ;
  kwargs...)`
- `JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps...;
  kwargs...)`

Here `aggregator` specifies the underlying algorithm for calculating next jump
times and types, for example `Direct`. The collection of different
`AbstractJump` types can then be passed within a single `JumpSet` or as
subsequent sequential arguments.

### Fields

$(FIELDS)

## Keyword Arguments
- `rng`, the random number generator to use. On 1.7 and up defaults to Julia's
  builtin generator, below 1.7 uses RandomNumbers.jl's
  `Xorshifts.Xoroshiro128Star(rand(UInt64))`.
- `save_positions=(true,true)`, specifies whether to save the system's state
  (before,after) the jump occurs.
- `spatial_system`, for spatial problems the underlying spatial structure.
- `hopping_constants`, for spatial problems the spatial transition rate
  coefficients.

Please see the [tutorial
page](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/) in
the DifferentialEquations.jl [docs](https://diffeq.sciml.ai/stable/) for usage
examples and commonly asked questions.
"""
mutable struct JumpProblem{iip,P,A,C,J<:Union{Nothing,AbstractJumpAggregator},J2,J3,J4,R} <: DiffEqBase.AbstractJumpProblem{P,J}
  """The type of problem to couple the jumps to. For a pure jump process use `DiscreteProblem`, to couple to ODEs, `ODEProblem`, etc."""
  prob::P
  """The aggregator algorithm that determines the next jump times and types for `ConstantRateJump`s and `MassActionJump`s. Examples include `Direct`."""
  aggregator::A
  """The underlying state data associated with the chosen aggregator."""
  discrete_jump_aggregation::J
  """`CallBackSet` with the underlying `ConstantRate` and `VariableRate` jumps."""
  jump_callback::C
  """The `VariableRateJump`s."""
  variable_jumps::J2
  """The `RegularJump`s."""
  regular_jump::J3
  """The `MassActionJump`s."""
  massaction_jump::J4
  """The random number generator to use."""
  rng::R
end
function JumpProblem(p::P,a::A,dj::J,jc::C,vj::J2,rj::J3,mj::J4,rng::R) where {P,A,J,C,J2,J3,J4,R}
    iip = isinplace_jump(p,rj)
    JumpProblem{iip,P,A,C,J,J2,J3,J4,R}(p,a,dj,jc,vj,rj,mj,rng)
end

# for remaking
Base.@pure remaker_of(prob::T) where {T <: JumpProblem} = DiffEqBase.parameterless_type(T)
function DiffEqBase.remake(thing::JumpProblem; kwargs...)
  T = remaker_of(thing)

  errmesg = """
  JumpProblems can currently only be remade with new u0, p, tspan or prob fields. To change other fields create a new JumpProblem. Feel free to open an issue on DiffEqJump to discuss further.
  """
  !issubset(keys(kwargs),(:u0,:p,:tspan,:prob)) && error(errmesg)

  if :prob ∉ keys(kwargs)
    dprob = DiffEqBase.remake(thing.prob; kwargs...)

    # if the parameters were changed we must remake the MassActionJump too
    if (:p ∈ keys(kwargs)) && using_params(thing.massaction_jump)
      update_parameters!(thing.massaction_jump, dprob.p; kwargs...)
    end
  else
    any(k -> k in keys(kwargs), (:u0,:p,:tspan)) && error("If remaking a JumpProblem you can not pass both prob and any of u0, p, or tspan.")
    dprob = kwargs[:prob]

    # we can't know if p was changed, so we must remake the MassActionJump
    if using_params(thing.massaction_jump)
      update_parameters!(thing.massaction_jump, dprob.p; kwargs...)
    end
  end

  T(dprob, thing.aggregator, thing.discrete_jump_aggregation, thing.jump_callback,
     thing.variable_jumps, thing.regular_jump, thing.massaction_jump, thing.rng)
end

DiffEqBase.isinplace(::JumpProblem{iip}) where {iip} = iip
JumpProblem(prob::JumpProblem) = prob

JumpProblem(prob,jumps::ConstantRateJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::VariableRateJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::RegularJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::MassActionJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::AbstractJump...;kwargs...) = JumpProblem(prob,JumpSet(jumps...);kwargs...)

JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::ConstantRateJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::VariableRateJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::RegularJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::AbstractMassActionJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::AbstractJump...;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps...);kwargs...)
JumpProblem(prob,jumps::JumpSet;kwargs...) = JumpProblem(prob,NullAggregator(),jumps;kwargs...)

function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet;
                     save_positions = typeof(prob) <: DiffEqBase.AbstractDiscreteProblem ? (false,true) : (true,true),
                     rng = DEFAULT_RNG, scale_rates = true, useiszero = true, spatial_system=nothing, hopping_constants=nothing, kwargs...)

  # initialize the MassActionJump rate constants with the user parameters
  if using_params(jumps.massaction_jump)
    rates = jumps.massaction_jump.param_mapper(prob.p)
    maj = MassActionJump(rates, jumps.massaction_jump.reactant_stoch, jumps.massaction_jump.net_stoch,
                         jumps.massaction_jump.param_mapper; scale_rates=scale_rates, useiszero=useiszero,
                         nocopy=true)
  else
    maj = jumps.massaction_jump
  end

  ## Spatial jumps handling
  if spatial_system !== nothing && hopping_constants !== nothing && !is_spatial(aggregator) # check if need to flatten
    prob, maj = flatten(maj, prob, spatial_system, hopping_constants; kwargs...)
  end
  ## Constant Rate Handling
  t,end_time,u = prob.tspan[1],prob.tspan[2],prob.u0
  if (typeof(jumps.constant_jumps) <: Tuple{}) && (maj === nothing) && !is_spatial(aggregator) # check if there are no jumps
    disc = nothing
    constant_jump_callback = CallbackSet()
  else
    disc = aggregate(aggregator,u,prob.p,t,end_time,jumps.constant_jumps,maj,save_positions,rng; spatial_system = spatial_system, hopping_constants = hopping_constants, kwargs...)
    constant_jump_callback = DiscreteCallback(disc)
  end

  iip = isinplace_jump(prob, jumps.regular_jump)

  ## Variable Rate Handling
  if typeof(jumps.variable_jumps) <: Tuple{}
    new_prob = prob
    variable_jump_callback = CallbackSet()
  else
    new_prob = extend_problem(prob,jumps; rng=rng)
    variable_jump_callback = build_variable_callback(CallbackSet(),0,jumps.variable_jumps...; rng=rng)
  end
  callbacks = CallbackSet(constant_jump_callback,variable_jump_callback)

  JumpProblem{iip,typeof(new_prob),typeof(aggregator),typeof(callbacks),
              typeof(disc),typeof(jumps.variable_jumps),
              typeof(jumps.regular_jump),typeof(maj),typeof(rng)}(
                        new_prob,aggregator,disc,
                        callbacks,
                        jumps.variable_jumps,
                        jumps.regular_jump, maj, rng)
end

function extend_problem(prob::DiffEqBase.AbstractDiscreteProblem, jumps; rng=DEFAULT_RNG)
  error("VariableRateJumps require a continuous problem, like an ODE/SDE/DDE/DAE problem.")
end

function extend_problem(prob::DiffEqBase.AbstractODEProblem,jumps; rng=DEFAULT_RNG)
  function jump_f(du::ExtendedJumpArray,u::ExtendedJumpArray,p,t)
    prob.f(du.u,u.u,p,t)
    update_jumps!(du,u,p,t,length(u.u),jumps.variable_jumps...)
  end
  ttype = eltype(prob.tspan)
  u0 = ExtendedJumpArray(prob.u0,[-randexp(rng,ttype) for i in 1:length(jumps.variable_jumps)])
  remake(prob,f=ODEFunction{true}(jump_f),u0=u0)
end

function extend_problem(prob::DiffEqBase.AbstractSDEProblem,jumps; rng=DEFAULT_RNG)
  function jump_f(du,u,p,t)
    prob.f(du.u,u.u,p,t)
    update_jumps!(du,u,p,t,length(u.u),jumps.variable_jumps...)
  end

  if prob.noise_rate_prototype === nothing
    jump_g = function (du,u,p,t)
      prob.g(du.u,u.u,p,t)
    end
  else
    jump_g = function (du,u,p,t)
      prob.g(du,u.u,p,t)
    end
  end

  ttype = eltype(prob.tspan)
  u0 = ExtendedJumpArray(prob.u0,[-randexp(rng,ttype) for i in 1:length(jumps.variable_jumps)])
  remake(prob,f=SDEFunction{true}(jump_f,jump_g),g=jump_g,u0=u0)
end

function extend_problem(prob::DiffEqBase.AbstractDDEProblem,jumps; rng=DEFAULT_RNG)
  jump_f = function (du,u,h,p,t)
    prob.f(du.u,u.u,h,p,t)
    update_jumps!(du,u,p,t,length(u.u),jumps.variable_jumps...)
  end
  ttype = eltype(prob.tspan)
  u0 = ExtendedJumpArray(prob.u0,[-randexp(rng,ttype) for i in 1:length(jumps.variable_jumps)])
  ramake(prob,f=DDEFunction{true}(jump_f),u0=u0)
end

# Not sure if the DAE one is correct: Should be a residual of sorts
function extend_problem(prob::DiffEqBase.AbstractDAEProblem,jumps; rng=DEFAULT_RNG)
  jump_f = function (out,du,u,p,t)
    prob.f(out.u,du.u,u.u,t)
    update_jumps!(du,u,t,length(u.u),jumps.variable_jumps...)
  end
  ttype = eltype(prob.tspan)
  u0 = ExtendedJumpArray(prob.u0,[-randexp(rng,ttype) for i in 1:length(jumps.variable_jumps)])
  remake(prob,f=DAEFunction{true}(jump_f),u0=u0)
end

function build_variable_callback(cb,idx,jump,jumps...; rng=DEFAULT_RNG)
  idx += 1
  condition = function (u,t,integrator)
    u.jump_u[idx]
  end
  affect! = function (integrator)
    jump.affect!(integrator)
    integrator.u.jump_u[idx] = -randexp(rng,typeof(integrator.t))
  end
  new_cb = ContinuousCallback(condition,affect!;
                      idxs = jump.idxs,
                      rootfind = jump.rootfind,
                      interp_points = jump.interp_points,
                      save_positions = jump.save_positions,
                      abstol = jump.abstol,
                      reltol = jump.reltol)
  build_variable_callback(CallbackSet(cb,new_cb),idx,jumps...; rng=DEFAULT_RNG)
end

function build_variable_callback(cb,idx,jump; rng=DEFAULT_RNG)
  idx += 1
  condition = function (u,t,integrator)
    u.jump_u[idx]
  end
  affect! = function (integrator)
    jump.affect!(integrator)
    integrator.u.jump_u[idx] = -randexp(rng,typeof(integrator.t))
  end
  new_cb = ContinuousCallback(condition,affect!;
                      idxs = jump.idxs,
                      rootfind = jump.rootfind,
                      interp_points = jump.interp_points,
                      save_positions = jump.save_positions,
                      abstol = jump.abstol,
                      reltol = jump.reltol)
  CallbackSet(cb,new_cb)
end

aggregator(jp::JumpProblem{P,A,C,J,J2}) where {P,A,C,J,J2} = A

@inline function extend_tstops!(tstops,jp::JumpProblem{P,A,C,J,J2}) where {P,A,C,J,J2}
  !(typeof(jp.jump_callback.discrete_callbacks) <: Tuple{}) && push!(tstops,jp.jump_callback.discrete_callbacks[1].condition.next_jump_time)
end

@inline function update_jumps!(du,u,p,t,idx,jump)
  idx += 1
  du[idx] = jump.rate(u.u,p,t)
end

@inline function update_jumps!(du,u,p,t,idx,jump,jumps...)
  idx += 1
  du[idx] = jump.rate(u.u,p,t)
  update_jumps!(du,u,p,t,idx,jumps...)
end

### Displays
num_constant_rate_jumps(aggregator::AbstractSSAJumpAggregator) = length(aggregator.rates)

Base.summary(io::IO, prob::JumpProblem) = string(DiffEqBase.parameterless_type(prob)," with problem ",DiffEqBase.parameterless_type(prob.prob)," and aggregator ",typeof(prob.aggregator))
function Base.show(io::IO, mime::MIME"text/plain", A::JumpProblem)
  println(io,summary(A))
  println(io,"Number of constant rate jumps: ",A.discrete_jump_aggregation === nothing ? 0 : num_constant_rate_jumps(A.discrete_jump_aggregation))
  println(io,"Number of variable rate jumps: ",length(A.variable_jumps))
  if A.regular_jump !== nothing
    println(io,"Have a regular jump")
  end
  if (A.massaction_jump !== nothing) && (get_num_majumps(A.massaction_jump) > 0)
    println(io,"Have a mass action jump")
  end
end

TreeViews.hastreeview(x::JumpProblem) = true


text:
using Deldir
using Test

import DataFrames

@testset "Deldir output are expected dataframes" begin
    N = rand(5:15)
    x = rand(N)
    y = rand(N)
    del, vor, summ = deldir(x, y)

    @test DataFrames.names(del) == ["x1", "y1", "x2", "y2", "ind1", "ind2"]

    @test DataFrames.names(vor) == ["x1", "y1", "x2", "y2", "ind1", "ind2", "bp1", "bp2", "thirdv1", "thirdv2"]

    @test DataFrames.names(summ) == ["x", "y", "ntri", "del_area", "n_tside", "nbpt", "vor_area"]
    @test DataFrames.nrow(summ) == N

    @test summ[!, :x] == x
    @test summ[!, :y] == y
end

@testset "Delaunay triangle corners are indexed correctly" begin
    N = rand(5:15)
    x = rand(N)
    y = rand(N)
    del = deldir(x, y)[1]

    @test del[!, :x1] == x[del[!, :ind1]]
    @test del[!, :x2] == x[del[!, :ind2]]

    @test del[!, :y1] == y[del[!, :ind1]]
    @test del[!, :y2] == y[del[!, :ind2]]
end

@testset "Line segments of Voronoi cells are within window" begin
    N = rand(5:15)
    x = rand(N)
    y = rand(N)
    vor = deldir(x, y)[2]

    @test all(0 .<= vor[!, :x1] .<= 1)
    @test all(0 .<= vor[!, :y1] .<= 1)
    @test all(0 .<= vor[!, :x2] .<= 1)
    @test all(0 .<= vor[!, :y2] .<= 1)
end

@testset "Area of Voronoi cells sum to area of window" begin
    N = rand(5:15)
    x = rand(N)
    y = rand(N)
    A = voronoiarea(x, y)

    @test sum(A) ≈ 1 atol = 0.001

    rw = [-rand(); 1 + rand(); -rand(); 1 + rand()]
    rw_area = (rw[2] - rw[1])*(rw[4] - rw[3])

    A = voronoiarea(x, y, rw)

    @test sum(A) ≈ rw_area atol = 0.001
end

@testset "Particular outputs" begin
    @testset "Simple point pattern" begin
        x = [0.25; 0.25; 0.75; 0.75]
        y = [0.25; 0.75; 0.25; 0.75]

        del, vor, summ = deldir(x, y)

        # Expected output is from R's deldir
        expected_del = DataFrames.DataFrame(
            x1 = [0.75; 0.75; 0.75; 0.25; 0.25],
            y1 = [0.25; 0.75; 0.75; 0.75; 0.75],
            x2 = [0.25; 0.25; 0.75; 0.25; 0.75],
            y2 = [0.25; 0.25; 0.25; 0.25; 0.75],
            ind1 = [3; 4; 4; 2; 2],
            ind2 = [1; 1; 3; 1; 4]
        )

        @test del == expected_del

        expected_vor = DataFrames.DataFrame(
            x1 = [0.5; 0.5; 0.5; 0.0; 0.5],
            y1 = [0.5; 0.5; 0.5; 0.5; 0.5],
            x2 = [0.5; 0.5; 1.0; 0.5; 0.5],
            y2 = [0.0; 0.5; 0.5; 0.5; 1.0],
            ind1 = [3; 4; 4; 2; 2],
            ind2 = [1; 1; 3; 1; 4],
            bp1 = [false; false; false; true; false],
            bp2 = [true; false; true; false; true],
            thirdv1 = [4; 2; 1; -2; 1],
            thirdv2 = [-1; 3; -4; 4; -3]
        )

        # There is no ≈ for DataFrames
        @test vor[!, "x1"] ≈ [0.5; 0.5; 0.5; 0.0; 0.5]
        @test vor[!, "y1"] ≈ [0.5; 0.5; 0.5; 0.5; 0.5]
        @test vor[!, "x2"] ≈ [0.5; 0.5; 1.0; 0.5; 0.5]
        @test vor[!, "y2"] ≈ [0.0; 0.5; 0.5; 0.5; 1.0]

        @test vor[!, ["ind1", "ind2", "bp1", "bp2", "thirdv1", "thirdv2"]] ==
            expected_vor[!, ["ind1", "ind2", "bp1", "bp2", "thirdv1", "thirdv2"]]

        expected_summ = DataFrames.DataFrame(
            x = [0.25; 0.25; 0.75; 0.75],
            y = [0.25; 0.75; 0.25; 0.75],
            ntri = [2; 1; 1; 2],
            del_area = [0.08333; 0.04166; 0.04166; 0.08333],
            n_tside = [3; 2; 2; 3],
            nbpt = [2; 2; 2; 2],
            vor_area = [0.25; 0.25; 0.25; 0.25]
        )

        # There is no ≈ for DataFrames
        @test summ[!, "del_area"] ≈ [0.08333; 0.04166; 0.04166; 0.08333] atol = 0.001
        @test summ[!, "vor_area"] ≈ [0.25; 0.25; 0.25; 0.25]
        @test summ[!, ["ntri", "n_tside", "nbpt"]] == expected_summ[!, ["ntri", "n_tside", "nbpt"]]
    end

    @testset "voronoiarea and deldir gives the same Voronoi area" begin
        # Points that are *not* sorted. From R documentation
        x = [2.3, 3.0, 7.0, 1.0, 3.0, 8.0]
        y = [2.3, 3.0, 2.0, 5.0, 8.0, 9.0]

        rw = [0.0; 10; 0; 10]
        A = voronoiarea(x, y, rw)

        @test A ≈ [11.737; 10.739; 26.839; 9.503; 21.306; 19.877] atol = 0.001

        summary = deldir(x, y, rw)[3]
        @test summary[!, "vor_area"] == A
    end

    @testset "Edges for plotting" begin
        x = [0.25; 0.75; 0.5]
        y = [0.25; 0.25; 0.75]

        del, vor, _ = deldir(x, y)

        Dx, Dy = edges(del)

        @test filter(!isnan, Dx) == [0.75; 0.25; 0.5; 0.25; 0.5; 0.75]
        @test filter(!isnan, Dy) == [0.25; 0.25; 0.75; 0.25; 0.75; 0.25]

        Vx, Vy = edges(vor)

        @test filter(!isnan, Vx) == [0.5; 0.5; 0.0; 0.5; 0.5; 1.0]
        @test filter(!isnan, Vy) ≈ [0.4375; 0.0; 0.6875; 0.4375; 0.4375; 0.6875] atol = 0.001
    end
end

@testset "Errors with inappropriate input" begin
    @testset "Error when points are outside window" begin
        x = [-rand(); rand()]
        y = rand(2)

        @test_throws ErrorException deldir(x, y)
    end

    @testset "Error when number of x's and y's are not equal" begin
        x = rand(rand(2:7))
        y = rand(rand(8:12))

        @test_throws DimensionMismatch deldir(x, y)
    end
end

@testset "Fortran errors" begin
    @testset "Triangle problems" begin
        # Data extracted from deldir::deldir documentation in R
        # In Fortran code used in Deldir_jll we have error code number 12, but in the current
        # version in the R package we get a more elaborate error mesage
        x = [0.21543139749966067; 0.18676067638651864; 0.12941923416171849; 0.37808260371144037; 0.08619595005015318; 0.15808995527500894]
        y = [1.0000000000000000; 0.9981701480225297; 0.9945104441215969; 0.4421766790515620; 0.9323236302262247; 0.9963402960710632]

        @test_throws ErrorException deldir(x, y)
    end

    @testset "Error number 4" begin
        # Data from GitHub issue #17
        x = [0.4, 0.3, 0.5, 0.2406, 0.2964, 0.5498, 0.2332, 0.3, 0.5041, 0.0824, 0.0594, 0.0126, 0.4385, 0.3575, 0.7737, 0.1, 0.1997, 0.6806, 0.8219, 0.0098, 0.4568, 0.0136]
        y = [0.3856, 0.5588, 0.0, 0.0725, 0.0433, 0.0025, 0.0771, 0.2124, 0.0, 0.2251, 0.7363, 0.3885, 0.0038, 0.0207, 0.0816, 0.2124, 0.1002, 0.0338, 0.3856, 0.4017, 0.0019, 0.616]

        @test_logs (:info, "Fortran error 4. Increasing madj to 24") deldir(x, y)
    end
end

@testset "Sort points" begin
    # Data extracted from deldir::deldir documentation in R where the points *are* shuffled.
    # Also used to test errors
    x = [0.21543139749966067; 0.18676067638651864; 0.12941923416171849; 0.37808260371144037; 0.08619595005015318; 0.15808995527500894]
    y = [1.0000000000000000; 0.9981701480225297; 0.9945104441215969; 0.4421766790515620; 0.9323236302262247; 0.9963402960710632]

    x_copy = deepcopy(x)
    y_copy = deepcopy(y)

    rw = [0.0; 1.0; 0.0; 1.0]
    indices, reverse_indices = Deldir.sortperm_points!(x_copy, y_copy, rw)

    @test x == x_copy[indices]
    @test x[reverse_indices] == x_copy

    @test y == y_copy[indices]
    @test y[reverse_indices] == y_copy
end


text:
#
# This implements Agha's example 3.2.1
#
using YAActL, Printf

# t = Ref{Task}()                # this is for debugging
# lk = LinkParams(taskref=t)

mutable struct StackNode{T}
    content::T
    link::Link
end

struct Pop <: Message
    customer::Link
end

struct Push{T} <: Message
    content::T
end

struct Print <: Message end

forward!(lk::L, msg::M) where {L<:Link, M<:Message} = send!(lk, msg)

# ----- this is essentially Agha's code example written in Julia/YAActL
function stack_node(sn::StackNode, msg::Pop)
    isnothing(sn.content) || become(forward!, sn.link)
    send!(msg.customer, Response(sn.content))
end

function stack_node(sn::StackNode, msg::Push)
    P = Actor(stack_node, sn)
    become(stack_node, StackNode(msg.content, P))
end
# ----------------------------------------------------------------------

# function stack_node(sn::StackNode, msg::Print)  # for debugging
#     print(@sprintf("content: %s\n", sn.content))
#     isnothing(sn.content) || send!(sn.link, msg)
# end

mystack = Actor(stack_node, StackNode(nothing, Link()))

response = Link()

send!(mystack, Pop(response))  # new stack
receive!(response).y           # returns nothing
send!(mystack, Push(1))        # push 1 on the stack
send!(mystack, Pop(response))  # pop it
receive!(response).y           # returns 1, 1st node now forwards messages
send!(mystack, Pop(response)); # pop again
receive!(response).y           # now returns nothing

for i ∈ 1:5
    send!(mystack, Push(i))
end

# (send!(mystack, Print()); sleep(0.1))

for i ∈ 1:5
    send!(mystack, Pop(response))
    println(receive!(response).y)
end


text:
module ACEfit

# data management
include("data.jl")

include("datautils.jl")

# iteration utilities over training data
include("iterate.jl")

# managing and visualising training and test errors

# loss functions and nonlinear solvers
# should this stay in ACEflux?

# linear least squares assembly for linear models
include("llsq.jl")

# linear least squares solvers for linear models
include("linearsolvers.jl")

end


text:
using QRupdate
using Test
using LinearAlgebra

@testset "qraddcol with β = 0" begin
    m = 100
    A = randn(m, 0)
    b = randn(m)
    R = zeros(0, 0)

    for i in 1:m
        a = randn(m)
        R = qraddcol(A, R, a)
        A = [A a]
        x, r = csne(R, A, b)
        @test norm(A'r, Inf) < 1e-5
        @test norm( R'*R - A'*A ) < 1e-5
    end
end

@testset "qraddcol with β > 0" begin
    m = 100
    A = randn(m, 0)
    b = randn(m)
    R = zeros(0, 0)
    β = 0.1
    for i in 1:m
        a = randn(m)
        R = qraddcol(A, R, a, β)
        A = [A a]
        # x, r = csne(R, A, b) # TODO: doesn't yet work with β > 0
        # @test norm(A'r, Inf) < 1e-5
        @test norm( R'*R - A'*A - β^2*I ) < 1e-5
    end
end

@testset "qraddrow" begin
    m, n = 3, 3
    A = randn(m,m)
    Q, R = qr(A)
    for i in 1:100
        a = randn(m)'
        A = [A; a]
        R = qraddrow(R, a)
        @test norm( R'R - A'*A ) < 1e-5
    end

end

@testset "qrdelcol" begin
    m = 100
    A = randn(m,m)
    Q, R = qr(A)
    for i in 100:-1:1
        k = rand(1:i)
        A = A[:,1:i .!= k]
        R = qrdelcol(R, k)
        @test norm( R'*R - A'*A ) < 1e-5
    end
end


text:
using Plots
using LinearAlgebra

theta = [pi/15, -pi/9, -pi/20]
Sigma = diagm([3.0, 1.0, 0.5]) # scale x, then y, then z

Rx = [1 0 0; 0 cos(theta[1]) -sin(theta[1]);  0 sin(theta[1]) cos(theta[1])] # rotate about x-axis
Ry = [cos(theta[2]) 0 sin(theta[2]); 0 1 0; -sin(theta[2]) 0 cos(theta[2])] # rotate about y-axis
Rz = [cos(theta[3]) -sin(theta[3]) 0; sin(theta[3]) cos(theta[3]) 0; 0 0 1] # rotate about z-axis

X = Rz*Ry*Rx*Sigma

# Plot sphere
n = 100
u = range(0, stop = 2pi, length = n)
v = range(0, stop = pi, length = n)

x = cos.(u) * sin.(v)'
y = sin.(u) * sin.(v)'
z = ones(n) * cos.(v)'

p1 = plot(x, y, z, legend = false)

xR = zeros(n, n)
yR = zeros(n, n)
zR = zeros(n, n)

for i = 1:n
    for j = 1:n
        vecR = X*[x[i,j]; y[i,j]; z[i,j]]
        xR[i,j] = vecR[1]
        yR[i,j] = vecR[2]
        zR[i,j] = vecR[3]
    end
end

p2 = plot(xR, yR, zR, legend = false)


text:
cfg = NeuroEvolution.get_config("../cfg/test.yaml")

@testset "Species" begin
    @test typeof(cfg["Species"])==Dict{Any,Any}
    @test length(cfg["Species"])==0
    ind = NEATIndiv(cfg)
    find_species!(ind, cfg)
    @test length(cfg["Species"])==1
    @test cfg["species_max"] == 1
end


text:
function get_filenames(; num_max_files::Integer = 0)
  files = readdir()
  datafiles = String[]
  strengths = Float64[]
  snrs = Float64[]
  phases = Float64[]
  for f in files
    m = match(r"(\w+)_(\d+)\w+_(\d+)\w+_(-?\d+\.\d*)\w+_(-?\d+\.\d*)\.csv",f)
    if m==nothing continue end
    effect_name = m.captures[1]
    push!(strengths,parse(m.captures[2]))
	  push!(snrs,parse(m.captures[3]))
    sampling = parse(m.captures[4])
    push!(phases,parse(m.captures[5]))
    push!(datafiles,f)
    if num_max_files>0 && length(datafiles)==num_max_files break end
  end
  @assert(length(datafiles)>0)
  idx = sortperm(phases)
  return (datafiles[idx],phases[idx])
end

function read_filelist(datafiles::Vector{String})
  data = readcsv(datafiles[1],skipstart=1);
  lambda = data[:,1]
  p = size(data,1)
  n = length(datafiles)
  X = zeros(p,n);
  for i in 1:n
    X[:,i] = readcsv(datafiles[i],skipstart=1)[:,2];
  end
  return (lambda,X)
end

# Mirrors code in  Sept 8, 2016 email from Allen Davis
function make_noisy_spectrum{T<:Real}(spec::Array{T,1}, snr::T; sampling::T = 1.0)
  snr_per_pixel = float(snr) / sqrt(sampling)
  orig_mean = mean(spec)
  scaled_spec = spec*snr*snr/orig_mean
  noisy_spec = scaled_spec + randn(length(spec)).*sqrt.(scaled_spec) #  np.random.normal(scaled_spec,np.sqrt.(scaled_spec))
  return orig_mean*noisy_spec/(snr*snr)
end

function make_noisy_spectra{T<:Real}(spec::Array{T,2}, snr::T; sampling::T = 1.0)
  output = similar(spec)
  for i in 1:size(spec,2)
    output[:,i] = make_noisy_spectrum(spec[:,i], snr, sampling=sampling)
  end
  return output
end

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


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

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

# ╔═╡ 0e592c68-9167-11eb-30fd-efdb809f47f7
#Set up the working environment
begin
	using Pkg
	Pkg.activate("..")
	using Plots, PlutoUI,LaTeXStrings, DifferentialEquations
end

# ╔═╡ e81c951c-9166-11eb-1dd3-51fe9fcb57db
md"""
# Modelling Predator-Prey Dynamics - Coupled Differential Equations

So far, we've looked at a couple of simple models that described the growth or decline of a single population. Now, let's turn towards a more complicated scenario, modelling of **predator-prey** dynamics...

A famous and very simple model of pedator-prey dynamics is the so-called **Lotka-Volterra** model:

$$\begin{eqnarray}
	\frac{dN}{dt} = N (a-bP)\\
	\frac{dP}{dt} = P (cN - d)
\end{eqnarray}$$

Here, $N$ represents the prey population and $P$ represents the population of predators.

Let us take a closer look at the assumptions baked into this model:
1. In the abesence of predators, the prey grows exponentially
2. Predation has the effect that it reduces the growth rate of the prey population in a manner that depends on $N$ and $P$.

3. If no prey is present, then the predator population declines in an exponential manner.

4. The growth rate of the predator depends on prey in a directly proportional manner.
"""

# ╔═╡ 4611ebd2-3a1e-4c62-80c7-e7dfc56483e6
md"""
> **Exercise - Equilibrium/Stability:** The given parameter set gives rise oscillating (regularly fluctuating) population of predators and prey.
> 1. Try to find one parameter set that gives rise to stable populations (equilibrium) of predators and prey by systematically varying the parameters.
> 2. What would equilibrium mean in terms of the differential equations that underlie the model?
> 3. Try to find the expressions for $N$ and $P$ that correspond to their equilibrium populations.
> 4. Equilibria can be stable or unstable, depending on how they respond to tiny changes in the parameters. What happens if any of the two initial condition deviates a little bit from the equilibrium condition?
> 5. In reality, population sizes, as well as the factors that are represented by the parameters are subject to stochastic (random) variation. This randomness is not represented in the model. Discuss the potential effects of randomness on real populations that would otherwise be well described by the model.
"""

# ╔═╡ 0c6eb46c-a03e-4950-abe6-368b85a4b55e
md"""
> **Exercise - Interspecies Relationships**: Unter certain parameters, the model gives results that appear strange at first sight: Peaks in predator populations vastly outnumber those of the prey population.
>
>1. Using the concept of energy flow in ecosystems, suggest why a predator population that vastly outnumbers the prey population is generally not what we would expect to see in an ecosystem
>2. Find a parameter set that gives you a predator population that is substantially larger than the prey population.
>3. Which other type of interspecies relationship could be adequately described by this model? Would the relative population sizes make more sense given this type of interspecies relationship?
"""

# ╔═╡ f1a68086-9169-11eb-2721-652c6732d36f
md""" ## Simulation
*Change the sliders to adjust the corresponding values*:

| Parameter               | Value                  |
|-------------------------|------------------------|
| prey growth constant $a$ | $(@bind a Slider(0:0.001:0.25, default=0.2, show_value=true)) |
| prey death constant $b$ | $(@bind b Slider(0:0.001:0.25, default=0.02, show_value=true)) |
| predator growth constant $c$ | $(@bind c Slider(0:0.001:0.25, default=0.005, show_value=true)) |
| predator death constant $d$ | $(@bind d Slider(0:0.001:0.25, default=0.1, show_value=true)) |
| initial prey population size $N_0$ | $(@bind N0 Slider(1:1:100, default=50, show_value=true)) |
| initial predator population size $P_0$ | $(@bind P0 Slider(1:1:100, default=10, show_value=true)) |
| simulation time $t_{max}$ | $(@bind tmax Slider(1:5:5000, default=100, show_value=true)) |
"""

# ╔═╡ 0f573286-916e-11eb-157e-1d0139bd363c
function lotkavolterra!(du,u,p,t)
	du[1] = u[1]*(p[1] - p[2]*u[2])
	du[2] = u[2]*(p[3]*u[1] - p[4])
end

# ╔═╡ 2400191e-916e-11eb-2d01-3b3a48757ff4
p = [a, b, c, d]

# ╔═╡ 2c79a5ea-916e-11eb-1abb-7933f2a1eae1
u0 = [N0, P0]

# ╔═╡ 2c7a590c-916e-11eb-227b-97d8098a22b3
tspan = (0, tmax)

# ╔═╡ 33519b3e-916e-11eb-1cb5-7fbc17fba1d8
prob = ODEProblem(lotkavolterra!,u0,tspan,p)

# ╔═╡ 6c09dc9c-916a-11eb-1bff-fb5495955ad7
begin
	# Defining the system of ODEs in paramtetrised manner and using in-place notation 		(! after function name) to improve performance
	# u[1] = N, u[2] = P

	sol = solve(prob)
	simplot = plot(sol, labels = ["N(t)" "P(t)"], title = "Timecourse")
	phaseplot = plot(sol, vars = (1,2), xlabel = "N(t)", ylabel = "P(t)", labels= "", title = "Phaseplane")
	plot(simplot, phaseplot, layout = (1, 2),fontfamily="Computer Modern")
end

# ╔═╡ 6acb8812-fd98-492b-ae20-38efd9cf8cba
md"""
#### Notebook options
Show solutions: $(@bind show_solutions CheckBox(default=false))
"""

# ╔═╡ 3bb1ca46-0993-47bf-9ba7-ad65185bac96
begin
	if show_solutions == true
md"""
!!! tip "Solution"

	**1.** A simple equilibrium can be found by setting any of the two growth rate constants to zero. This would correspond to the extinction of the two populations.

	**2.** In terms of the equations, equilibrium corresponds to a situation where neither of the two populations grows or declines, hence both expressions ($\frac{dN}{dt}$ and $\frac{dP}{dt}$) are zero:

	$$\begin{eqnarray}
		\frac{dN}{dt} = 0 = N (a-bP)\\
		\frac{dP}{dt} = 0 = P (cN - d)
	\end{eqnarray}$$

	**3.** Using the insight from above and solving for $N$ and $P$, we find:

	$$\{N_{eq,1} = 0, P_{eq,1} = 0\}$$

	and

	$$\{N_{eq,2} = \frac{d}{c}, P_{eq,2} = \frac{a}{b}\}$$

	**4.** Tiny deviations instantly give rise to oscillations around the equilibrium point. The populations never reach equilibrium.

	**5.** Due to the tendency of the populations to oscillate around the non-zero equilibrium, we would not expect to see non-oscillation populations.
"""
	end
end

# ╔═╡ f6a109fd-232e-42b6-a2c2-88c7eb9befce
begin
	if show_solutions == true
md"""
!!! tip "Solution"
	**1.** Since only a fraction of the energy is transfered from lower to higher trophic levels, the situation described by the model does not reflect expected population sizes for typical predator-prey relationships.
	**2.** Many combinations work, e.g. $a$ = 0.2, $b$ = 0.02, $c$ = 0.2, $d$ = 0.2, $N_0$ = 10, $P_0$ = 40
	**3.** Another antagonistic (-/-) relationship that could be described by this model is parasitism. In this case the relative population sizes would make more sense.
"""
	end
end

# ╔═╡ Cell order:
# ╟─0e592c68-9167-11eb-30fd-efdb809f47f7
# ╟─e81c951c-9166-11eb-1dd3-51fe9fcb57db
# ╟─4611ebd2-3a1e-4c62-80c7-e7dfc56483e6
# ╟─3bb1ca46-0993-47bf-9ba7-ad65185bac96
# ╟─0c6eb46c-a03e-4950-abe6-368b85a4b55e
# ╟─f6a109fd-232e-42b6-a2c2-88c7eb9befce
# ╟─f1a68086-9169-11eb-2721-652c6732d36f
# ╟─6c09dc9c-916a-11eb-1bff-fb5495955ad7
# ╟─0f573286-916e-11eb-157e-1d0139bd363c
# ╟─2400191e-916e-11eb-2d01-3b3a48757ff4
# ╟─2c79a5ea-916e-11eb-1abb-7933f2a1eae1
# ╟─2c7a590c-916e-11eb-227b-97d8098a22b3
# ╟─33519b3e-916e-11eb-1cb5-7fbc17fba1d8
# ╟─6acb8812-fd98-492b-ae20-38efd9cf8cba

max_stars_repo_name: el-uhu/intro-to-modelling

max_issues_repo_name: el-uhu/intro-to-modelling

max_forks_repo_name: el-uhu/intro-to-modelling


text:
function distance_threshold(overlap_len)
    return min(6, overlap_len/10.0)
end

"""
calculate the overlap length of a pair of reads
"""
function overlap(r1::Sequence, r2::Sequence)
    len1 = length(r1)
    len2 = length(r2)
    reverse_r2 = ~r2

    overlapped = false
    overlap_len = 0
    offset = 0
    distance = 0
    # a match of less than 10 is considered as unconfident
    for offset = 0:len1-10
        # the overlap length of r1 & r2 when r2 is move right for offset
        overlap_len = min(len1-offset, len2)

        # remind that Julia is a 1-based coordination system
        distance = edit_distance(r1.seq[offset+1 : offset+overlap_len], reverse_r2.seq[1 : overlap_len])
        if distance <= distance_threshold(overlap_len)
            # now we find a good candidate
            # we verify it by moving r2 one more base to see if the distance is getting longer
            # if yes, then current is the best match, otherwise it's not
            next = offset + 1
            next_overlap_len = min(len1-next, len2)
            distance2 = edit_distance(r1.seq[next+1 : next+next_overlap_len], reverse_r2.seq[1 : next_overlap_len])
            if distance <= distance2
                overlapped = true
                break
            end
        end
    end

    if overlapped && offset == 0
        # check if distance can get smaller if offset goes negative
        # this only happens when insert DNA is shorter than sequencing read length, and some adapter/primer is sequenced but not trimmed cleanly
        # we go reversely
        for offset = 0:-1:-(len2-10)
            # the overlap length of r1 & r2 when r2 is move right for offset
            overlap_len = min(len1,  len2-abs(offset))
            distance = edit_distance(r1.seq[1:overlap_len], reverse_r2.seq[-offset + 1 : -offset + overlap_len])
            if distance <= distance_threshold(overlap_len)
                next = offset - 1
                next_overlap_len = min(len1,  len2-abs(next))
                distance2 = edit_distance(r1.seq[1:next_overlap_len], reverse_r2.seq[-next + 1 : -next + next_overlap_len])
                if distance <= distance2
                    return (offset, overlap_len, distance)
                end
            end
        end
    elseif overlapped
        return (offset, overlap_len, distance)
    end

    return (0,0,0)
end

function overlap(r1::FastqRead, r2::FastqRead)
    return overlap(r1.sequence, r2.sequence)
end

function overlap(pair::FastqPair)
    return overlap(pair.read1, pair.read2)
end

function simple_merge(r1::Sequence, r2::Sequence, overlap, offset = 0)
    reverse = ~r2
    if offset >= 0
        if overlap < length(r2)
            return Sequence(r1.seq * reverse.seq[overlap+1 : end])
        else
            return r1
        end
    else
        # offset is negative means insert DNA is shorter than read length
        # just shift the junk bases in the tail of r1
        return Sequence(r1.seq[1:overlap])
    end
end

function try_merge(r1::Sequence, r2::Sequence; distance_threshold = 2, overlap_threshold = 20)
    offset, overlap_len, distance = overlap(r1, r2)
    if overlap_len < overlap_threshold || distance > distance_threshold
        return false
    else
        return simple_merge(r1, r2, overlap_len, offset)
    end
end