text:
##########################################################################################
#
#    function 'derive' returning the expr of gradient
#
#
##########################################################################################
# TODO : add operators : hcat, vcat, ? : , map, mapreduce, if else

## macro to simplify derivation rules creation
macro dfunc(func::Expr, dv::Symbol, diff::Expr)
	argsn = map(e-> isa(e, Symbol) ? e : e.args[1], func.args[2:end])
	index = find(dv .== argsn)[1]

	# change var names in signature and diff expr to x1, x2, x3, ..
	smap = { argsn[i] => symbol("x$i") for i in 1:length(argsn) }
	args2 = substSymbols(func.args[2:end], smap)

	# diff function name
	fn = symbol("d_$(func.args[1])_x$index")

	fullf = Expr(:(=), Expr(:call, fn, args2...), Expr(:quote, substSymbols(diff, smap)) )
	eval(fullf)
end

## common operators
# TODO : check if ds can be removed for clarity

@dfunc +(x::Real, y)     x     sum(ds)
@dfunc +(x::Array, y)    x     +ds
@dfunc +(x, y::Real)     y     sum(ds)
@dfunc +(x, y::Array)    y     +ds

@dfunc -x                x     -ds
@dfunc -(x::Real, y)     x     sum(ds)
@dfunc -(x::Array, y)    x     +ds
@dfunc -(x, y::Real)     y     -sum(ds)
@dfunc -(x, y::Array)    y     -ds

@dfunc sum(x)       x     +ds
@dfunc dot(x, y)    x     y .* ds
@dfunc dot(x, y)    y     x .* ds

@dfunc log(x)       x     ds ./ x
@dfunc exp(x)       x     exp(x) .* ds

@dfunc sin(x)       x     cos(x) .* ds
@dfunc cos(x)       x     -sin(x) .* ds

@dfunc abs(x)       x     sign(x) .* ds

@dfunc *(x::Real, y)     x     sum(ds .* y)
@dfunc *(x::Array, y)    x     ds * transpose(y)
@dfunc *(x, y::Real)     y     sum(ds .* x)
@dfunc *(x, y::Array)    y     transpose(x) * ds

@dfunc .*(x::Real, y)    x     sum(ds .* y)
@dfunc .*(x::Array, y)   x     ds .* y
@dfunc .*(x, y::Real)    y     sum(ds .* x)
@dfunc .*(x, y::Array)   y     ds .* x

@dfunc ^(x::Real, y::Real)  x     y * x ^ (y-1) * ds # Both args reals
@dfunc ^(x::Real, y::Real)  y     log(x) * x ^ y * ds # Both args reals

@dfunc .^(x::Real, y)    x     sum(y .* x .^ (y-1) .* ds)
@dfunc .^(x::Array, y)   x     y .* x .^ (y-1) .* ds
@dfunc .^(x, y::Real)    y     sum(log(x) .* x .^ y .* ds)
@dfunc .^(x, y::Array)   y     log(x) .* x .^ y .* ds

@dfunc /(x::Real, y)          x     sum(ds ./ y)
@dfunc /(x::Array, y::Real)   x     ds ./ y
@dfunc /(x, y::Real)          y     sum(- x ./ (y .* y) .* ds)
@dfunc /(x::Real, y::Array)   y     (- x ./ (y .* y)) .* ds

@dfunc ./(x::Real, y)        x     sum(ds ./ y)
@dfunc ./(x::Array, y)       x     ds ./ y
@dfunc ./(x, y::Real)        y     sum(- x ./ (y .* y) .* ds)
@dfunc ./(x, y::Array)       y     (- x ./ (y .* y)) .* ds

@dfunc max(x::Real, y)    x     sum((x .> y) .* ds)
@dfunc max(x::Array, y)   x     (x .> y) .* ds
@dfunc max(x, y::Real)    y     sum((x .< y) .* ds)
@dfunc max(x, y::Array)   y     (x .< y) .* ds

@dfunc min(x::Real, y)    x     sum((x .< y) .* ds)
@dfunc min(x::Array, y)   x     (x .< y) .* ds
@dfunc min(x, y::Real)    y     sum((x .> y) .* ds)
@dfunc min(x, y::Array)   y     (x .> y) .* ds

@dfunc transpose(x::Real)   x   +ds
@dfunc transpose(x::Array)  x   transpose(ds)

## Normal distribution
@dfunc logpdfNormal(mu::Real, sigma, x)    mu     sum((x - mu) ./ (sigma .* sigma)) * ds
@dfunc logpdfNormal(mu::Array, sigma, x)   mu     (x - mu) ./ (sigma .* sigma) * ds
@dfunc logpdfNormal(mu, sigma::Real, x)    sigma  sum(((x - mu).*(x - mu) ./ (sigma.*sigma) - 1.) ./ sigma) * ds
@dfunc logpdfNormal(mu, sigma::Array, x)   sigma  ((x - mu).*(x - mu) ./ (sigma.*sigma) - 1.) ./ sigma * ds
@dfunc logpdfNormal(mu, sigma, x::Real)    x      sum((mu - x) ./ (sigma .* sigma)) * ds
@dfunc logpdfNormal(mu, sigma, x::Array)   x      (mu - x) ./ (sigma .* sigma) * ds

## Uniform distribution
@dfunc logpdfUniform(a::Real, b, x)      a   sum((a .<= x .<= b) ./ (b - a)) * ds
@dfunc logpdfUniform(a::Array, b, x)     a   ((a .<= x .<= b) ./ (b - a)) * ds
@dfunc logpdfUniform(a, b::Real, x)      b   sum((a .<= x .<= b) ./ (a - b)) * ds
@dfunc logpdfUniform(a, b::Array, x)     b   ((a .<= x .<= b) ./ (a - b)) * ds
@dfunc logpdfUniform(a, b, x)            x   zero(x)

## Weibull distribution
@dfunc logpdfWeibull(sh::Real, sc, x)    sh  (r = x./sc ; sum(((1. - r.^sh) .* log(r) + 1./sh)) * ds)
@dfunc logpdfWeibull(sh::Array, sc, x)   sh  (r = x./sc ; ((1. - r.^sh) .* log(r) + 1./sh) * ds)
@dfunc logpdfWeibull(sh, sc::Real, x)    sc  sum(((x./sc).^sh - 1.) .* sh./sc) * ds
@dfunc logpdfWeibull(sh, sc::Array, x)   sc  ((x./sc).^sh - 1.) .* sh./sc * ds
@dfunc logpdfWeibull(sh, sc, x::Real)    x   sum(((1. - (x./sc).^sh) .* sh - 1.) ./ x) * ds
@dfunc logpdfWeibull(sh, sc, x::Array)   x   ((1. - (x./sc).^sh) .* sh - 1.) ./ x * ds

## Beta distribution
@dfunc logpdfBeta(a, b, x::Real)      x     sum((a-1) ./ x - (b-1) ./ (1-x)) * ds
@dfunc logpdfBeta(a, b, x::Array)     x     ((a-1) ./ x - (b-1) ./ (1-x)) * ds
@dfunc logpdfBeta(a::Real, b, x)      a     sum(digamma(a+b) - digamma(a) + log(x)) * ds
@dfunc logpdfBeta(a::Array, b, x)     a     (digamma(a+b) - digamma(a) + log(x)) * ds
@dfunc logpdfBeta(a, b::Real, x)      b     sum(digamma(a+b) - digamma(b) + log(1-x)) * ds
@dfunc logpdfBeta(a, b::Array, x)     b     (digamma(a+b) - digamma(b) + log(1-x)) * ds

## TDist distribution
@dfunc logpdfTDist(df, x::Real)     x     sum(-(df+1).*x ./ (df+x.*x)) .* ds
@dfunc logpdfTDist(df, x::Array)    x     (-(df+1).*x ./ (df+x.*x)) .* ds
@dfunc logpdfTDist(df::Real, x)     df    (tmp2 = (x.*x + df) ; sum( (x.*x-1)./tmp2 + log(df./tmp2) + digamma((df+1)/2) - digamma(df/2) ) / 2 .* ds )
@dfunc logpdfTDist(df::Array, x)    df    (tmp2 = (x.*x + df) ; ( (x.*x-1)./tmp2 + log(df./tmp2) + digamma((df+1)/2) - digamma(df/2) ) / 2 .* ds )

## Exponential distribution
@dfunc logpdfExponential(sc, x::Real)    x   sum(-1/sc) .* ds
@dfunc logpdfExponential(sc, x::Array)   x   (- ds ./ sc)
@dfunc logpdfExponential(sc::Real, x)    sc  sum((x-sc)./(sc.*sc)) .* ds
@dfunc logpdfExponential(sc::Array, x)   sc  (x-sc) ./ (sc.*sc) .* ds

## Gamma distribution
@dfunc logpdfGamma(sh, sc, x::Real)    x   sum(-( sc + x - sh.*sc)./(sc.*x)) .* ds
@dfunc logpdfGamma(sh, sc, x::Array)   x   (-( sc + x - sh.*sc)./(sc.*x)) .* ds
@dfunc logpdfGamma(sh::Real, sc, x)    sh  sum(log(x) - log(sc) - digamma(sh)) .* ds
@dfunc logpdfGamma(sh::Array, sc, x)   sh  (log(x) - log(sc) - digamma(sh)) .* ds
@dfunc logpdfGamma(sh, sc::Real, x)    sc  sum((x - sc.*sh) ./ (sc.*sc)) .* ds
@dfunc logpdfGamma(sh, sc::Array, x)   sc  ((x - sc.*sh) ./ (sc.*sc)) .* ds

## Cauchy distribution
@dfunc logpdfCauchy(mu, sc, x::Real)    x   sum(2(mu-x) ./ (sc.*sc + (x-mu).*(x-mu))) .* ds
@dfunc logpdfCauchy(mu, sc, x::Array)   x   (2(mu-x) ./ (sc.*sc + (x-mu).*(x-mu))) .* ds
@dfunc logpdfCauchy(mu::Real, sc, x)    mu  sum(2(x-mu) ./ (sc.*sc + (x-mu).*(x-mu))) .* ds
@dfunc logpdfCauchy(mu::Array, sc, x)   mu  (2(x-mu) ./ (sc.*sc + (x-mu).*(x-mu))) .* ds
@dfunc logpdfCauchy(mu, sc::Real, x)    sc  sum(((x-mu).*(x-mu) - sc.*sc) ./ (sc.*(sc.*sc + (x-mu).*(x-mu)))) .* ds
@dfunc logpdfCauchy(mu, sc::Array, x)   sc  (((x-mu).*(x-mu) - sc.*sc) ./ (sc.*(sc.*sc + (x-mu).*(x-mu)))) .* ds

## Log-normal distribution
@dfunc logpdfLogNormal(lmu, lsc, x::Real)   x    ( tmp2=lsc.*lsc ; sum( (lmu - tmp2 - log(x)) ./ (tmp2.*x) ) .* ds )
@dfunc logpdfLogNormal(lmu, lsc, x::Array)  x    ( tmp2=lsc.*lsc ; ( (lmu - tmp2 - log(x)) ./ (tmp2.*x) ) .* ds )
@dfunc logpdfLogNormal(lmu::Real, lsc, x)   lmu  sum((log(x) - lmu) ./ (lsc .* lsc)) .* ds
@dfunc logpdfLogNormal(lmu::Array, lsc, x)  lmu  ((log(x) - lmu) ./ (lsc .* lsc)) .* ds
@dfunc logpdfLogNormal(lmu, lsc::Real, x)   lsc  ( tmp2=lsc.*lsc ; sum( (lmu.*lmu - tmp2 - log(x).*(2lmu-log(x))) ./ (lsc.*tmp2) ) .* ds )
@dfunc logpdfLogNormal(lmu, lsc::Array, x)  lsc  ( tmp2=lsc.*lsc ; ( (lmu.*lmu - tmp2 - log(x).*(2lmu-log(x))) ./ (lsc.*tmp2) ) .* ds )

# TODO : find a way to implement multi variate distribs that goes along well with vectorization (Dirichlet, Categorical)
# TODO : other continuous distribs ? : Pareto, Rayleigh, Logistic, Levy, Laplace, Dirichlet, FDist
# TODO : other discrete distribs ? : NegativeBinomial, DiscreteUniform, HyperGeometric, Geometric, Categorical

## Bernoulli distribution (Note : no derivation on x parameter as it is an Int)
@dfunc logpdfBernoulli(p::Real, x)     p     sum(1. ./ (p - (1. - x))) * ds
@dfunc logpdfBernoulli(p::Array, x)    p     (1. ./ (p - (1. - x))) * ds

## Binomial distribution (Note : no derivation on x and n parameters as they are integers)
@dfunc logpdfBinomial(n, p::Real, x)   p    sum(x ./ p - (n-x) ./ (1 - p)) * ds
@dfunc logpdfBinomial(n, p::Array, x)  p    (x ./ p - (n-x) ./ (1 - p)) * ds

## Poisson distribution (Note : no derivation on x parameter as it is an Int)
@dfunc logpdfPoisson(lambda::Real, x)   lambda   sum(x ./ lambda - 1) * ds
@dfunc logpdfPoisson(lambda::Array, x)  lambda   (x ./ lambda - 1) * ds

## returns sample value for the given Symbol or Expr (for refs)
hint(v::Symbol) = vhint[v]
hint(v) = v  # should be a value if not a Symbol or an Expression
function hint(v::Expr)
	@assert v.head == :ref "[hint] unexpected variable $v"
	v.args[1] = :( vhint[$(Expr(:quote, v.args[1]))] )
	eval(v)
end

## Returns gradient expression of opex
function derive(opex::Expr, index::Int, dsym::Union(Expr,Symbol))  # opex=:(z^x);index=2;dsym=:y
	vs = opex.args[1+index]
	ds = symbol("$DERIV_PREFIX$dsym")
	args = opex.args[2:end]

	val = map(hint, args)  # get sample values of args to find correct gradient statement

	fn = symbol("d_$(opex.args[1])_x$index")

	try
		dexp = eval(Expr(:call, fn, val...))

		smap = { symbol("x$i") => args[i] for i in 1:length(args)}
		smap[:ds] = ds
		dexp = substSymbols(dexp, smap)

		# unfold for easier optimization later
	    m = ParsingStruct()
	    m.source = :(dummy = $dexp )
		unfold!(m)
		m.exprs[end] = m.exprs[end].args[2] # remove last assignment

		m.exprs[end] = :( $(symbol("$DERIV_PREFIX$vs")) = $(symbol("$DERIV_PREFIX$vs")) + $(m.exprs[end]) )
		return m.exprs
	catch e
		error("[derive] Doesn't know how to derive $opex by argument $vs")
	end

end


text:
"""
Все клеткт поля должны быть промаркерованы
"""

using HorizonSideRobots
originRobot = Robot(animate = false, "23example.sit")
include("../AbstractType.jl")
include("../functions.jl")

function mod_moveLikeSnake!(robot::Robot, moveSides::NTuple{2,HorizonSide}, borderSide::HorizonSide)::Nothing
    flag = true
    while flag
        for side in moveSides
            while try_move!(robot, side)
                putmarker!(robot)
            end
            if !(isborder(robot, borderSide))
                move!(robot, borderSide)
                putmarker!(robot)
            else
                flag = false
            end
        end
    end
end

function master!(robot)::Nothing
    backPath = moving_in_angle!(robot, West, Sud)

    putmarker!(robot)

    mod_moveLikeSnake!(robot, (Ost, West), Nord)

    moving_in_angle!(robot, West, Sud)

    moving_back_to_start!(robot, backPath)
end

master!(originRobot)
show!(originRobot)

ext: jl

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
for N in [Float64, Float32]
    # dimension (choose a multiple of 3)
    n = 3*2

    # sparse matrix
    m = spzeros(N, n, n)
    m[1, 1] = N(1//25)
    m[1, 2] = N(1//5)
    m[3, 2] = N(-8//5)
    m[4, 2] = N(-3//2)
    m[5, 2] = N(-6//5)
    m[6, 3] = N(-3//10)
    m[2, 4] = N(1//2)
    m[5, 4] = N(-3//4)
    m[6, 4] = N(-1//2)
    m[3, 5] = N(-3//4)
    m[2, 6] = N(4/5)
    m[3, 6] = N(7//10)
    m[4, 6] = N(-1//2)
    m[5, 6] = N(3//10)

    # a set
    b = BallInf(ones(N, n), N(0.1))

    # linear map
    lm = m * b

    # the (lazy) exponential of a sparse matrix
    me = SparseMatrixExp(m)

    # constructor from a dense matrix
    @test_throws ErrorException SparseMatrixExp(Matrix{N}(I, 2, 2))

    # constructor from a non-square matrix
    # we would use sprandn(N, 2, 3, 0.1) but it is available from Julia v1.1
    mr = randn(N, 2, 3)
    @test_throws AssertionError SparseMatrixExp(sparse(mr))

    # size
    @test size(me, 1) == n
    @test size(me) == (n, n)

    # columns & rows
    me2 = SparseMatrixExp(sparse(N(1) * I, n, n))
    v = zeros(N, n)
    v[1] = N(ℯ)
    @test get_columns(me2, [1, 2])[:, 1] == get_column(me2, 1) == v
    @test get_rows(me2, [1, 2])[1, :] == transpose(get_row(me2, 1)) == v

    # the exponential map of a convex set, same as ExponentialMap(me, b)
    emap = me * b
    @test emap isa ExponentialMap{N}

    # absorbing elements
    for neutral in [ZeroSet{N}(6), EmptySet{N}(6)]
        @test me * neutral == neutral
    end

    # dimension
    @test dim(emap) == n

    # the support vector of an exponential map
    d = N[2, -1, -1, 0, 1, -1]
    svec = σ(d, emap)
    # check that it works with sparse vectors
    σ(sparsevec(d), emap)

    # check consistency with respect to explicit computation of the matrix exponential
    svec_explicit = σ(d, exp(Matrix(m)) * b)
    @test svec ≈ svec_explicit

    # concretize
    X = concretize(emap)
    @test σ(d, X) ≈ svec

    # boundedness
    @test isbounded(emap) && isboundedtype(typeof(emap))
    emap2 = me * HalfSpace(ones(N, n), N(1))
    @test !isbounded(emap2) && !isboundedtype(typeof(emap2))

    # is_polyhedral
    @test is_polyhedral(emap)
    if N isa AbstractFloat
        emap2 = me * Ball2(zeros(N, n), N(1))
        @test !is_polyhedral(emap2)
    end

    # isempty
    @test !isempty(emap)

    # membership
    x = ones(N, n)
    @test x ∈ b
    @test me2 * x ∈ me2 * b
    @test -x ∉ b
    @test me2 * -x ∉ me2 * b

    # construct an exponential map where we only pass the matrix
    #ExponentialMap(m, b) # ExponentialMap

    # the exponential map of an exponential map falls back to a Cartesian product
    cpem = emap * emap # CartesianProduct of ExponentialMap

    # the exponential map applied to an exponential map
    ExponentialMap(SparseMatrixExp(2*m), emap)

    # for projection tests, let's assume that n is divisible by three
    @assert mod(n, 3) == 0
    nb = div(n, 3)
    # the projection of exp(A) on the (m, m)-dimensional right-most upper block
    R = [spzeros(N, nb, nb); spzeros(N, nb, nb); sparse(N(1) * I, nb, nb)]
    L = [sparse(N(1) * I, nb, nb) spzeros(N, nb, nb) spzeros(N, nb, nb)]
    proj = ProjectionSparseMatrixExp(L, me, R)

    # build an exponential projection map : it is the application of the projection
    # of exp(A) over a given set
    b = BallInf(ones(N, nb), N(1))
    projmap = proj * b

    # query the ambient dimension of projmap (hint: it is the output dimension)
    @test dim(projmap) == nb

    # boundedness
    @test isbounded(projmap)
    @test isbounded(ProjectionSparseMatrixExp(spzeros(N, nb, nb), me, R) * HalfSpace(ones(N, nb), N(1)))
    @test isbounded(ProjectionSparseMatrixExp(L, me, spzeros(N, nb, nb)) * HalfSpace(ones(N, nb), N(1)))
    @test !isbounded(proj * HalfSpace(ones(N, nb), N(1)))

    # isempty
    @test !isempty(projmap)

    #compute the support vector of the projection of an exponential map
    d = N[3//20, 11//50]
    svec = σ(d, projmap)
    # check that it works with sparse vectors
    svec == σ(sparsevec(d), projmap)

    # check consistency with respect to explicit computation of the matrix exponential
    P = L * exp(Matrix(m)) * R
    svec_explicit = σ(d, P*b)
    @test svec ≈ svec_explicit

    # vertices_list
    b = BallInf(N[0, 0], N(1))
    M = SparseMatrixExp(spzeros(N, 2, 2))
    vlist = vertices_list(ExponentialMap(M, b))
    @test ispermutation(vlist, [N[1, 1], N[-1, 1], N[1, -1], N[-1, -1]])
end


text:
# delegate methods for SubArrays to support view

DataAPI.levels(sa::SubArray{T,N,P}) where {T,N,P<:CategoricalArray} = levels(parent(sa))
isordered(sa::SubArray{T,N,P}) where {T,N,P<:CategoricalArray} = isordered(parent(sa))
# This method cannot support allowmissing=true since that would modify the parent
levels!(sa::SubArray{T,N,P}, newlevels::Vector) where {T,N,P<:CategoricalArray} =
    levels!(parent(sa), newlevels)

function unique(sa::SubArray{T,N,P}) where {T,N,P<:CategoricalArray}
    A = parent(sa)
    refs = view(A.refs, sa.indices...)
    S = eltype(P) >: Missing ? Union{eltype(levels(A.pool)), Missing} : eltype(levels(A.pool))
    _unique(S, refs, A.pool)
end

refs(A::SubArray{<:Any, <:Any, <:CategoricalArray}) =
    view(parent(A).refs, parentindices(A)...)

pool(A::SubArray{<:Any, <:Any, <:CategoricalArray}) = A.parent.pool

function Base.fill!(A::SubArray{<:Any, <:Any, <:CategoricalArray}, v::Any)
    r = get!(A.parent.pool, convert(leveltype(A.parent), v))
    fill!(refs(A), r)
    A
end

Base.fill!(A::SubArray{<:Any, <:Any, <:CategoricalArray{>:Missing}}, ::Missing) =
    (fill!(refs(A), 0); A)

Base.Broadcast.broadcasted(::typeof(ismissing),
                           A::SubArray{<:Any, <:Any, <:CategoricalArray{T}}) where {T} =
    T >: Missing ? Base.Broadcast.broadcasted(==, refs(A), 0) :
                   Base.Broadcast.broadcasted(_ -> false, refs(A))

Base.Broadcast.broadcasted(::typeof(!ismissing),
                           A::SubArray{<:Any, <:Any, <:CategoricalArray{T}}) where {T} =
    T >: Missing ? Base.Broadcast.broadcasted(>, refs(A), 0) :
                   Base.Broadcast.broadcasted(_ -> true, refs(A))


text:
# This file contains the code to generate the PWS data used for
# the paper.

# WARNING: This script takes a long time (hours to days) to complete
# This file should be executed on multiple CPU's (e.g.
# on a computing cluster). To use all cores of the current
# computer, start julia with `julia -p auto generate_data.jl`.

using Logging

@info "Loading Packages"
using HDF5
@everywhere using PathWeightSampling
using DrWatson
using Distributions
using ProgressMeter

duration = 10
parameters = (;
    mean_s = 50,
    corr_time_s = 1.0,
    corr_time_x = 0.1,
    dtimes = collect(0.0:0.1:duration)
)

mean_s = parameters.mean_s
κ = parameters.mean_s / parameters.corr_time_s
λ = 1 / parameters.corr_time_s
ρ = 1 / parameters.corr_time_x
μ = 1 / parameters.corr_time_x

# Use a Gaussian approximation for the initial condition
@info "Compute Initial Condition"
covariance = [mean_s (ρ*mean_s / (λ + μ)); (ρ*mean_s / (λ + μ)) (mean_s * (1 + ρ/(λ + μ)))]
initial_dist = MvNormal([mean_s, mean_s], covariance)
s_axis = 25:75
x_axis = 25:75
z = pdf.(Ref(initial_dist), [[i,j] for i=s_axis, j=x_axis])
initial_condition = PathWeightSampling.empirical_dist(z, s_axis, x_axis)

function run_simulation(systemfn, algorithm, num_responses)
    @everywhere begin
        system = $systemfn()
        global compiled_system = PathWeightSampling.compile(system)
    end

    result = @showprogress pmap(1:num_responses) do batch
        PathWeightSampling._mi_inner(Main.compiled_system, algorithm, 1, false)
    end
    vcat(result...)
end

# This function runs a simulation and saves it in the data directory
function save_simulation(alg::Symbol, N::Integer, M::Integer, parameters)
    systemfn = () -> PathWeightSampling.gene_expression_system(; u0=initial_condition, parameters...)

    algorithm = if alg == :directmc
        DirectMCEstimate(M)
    elseif alg == :smc
        SMCEstimate(M)
    elseif alg == :ti
        TIEstimate(0, 8, M÷8)
    else
        error("No algorithm $alg found.")
    end

    mi = run_simulation(systemfn, algorithm, N)
    result = Dict("Samples" => mi, "Parameters" => Dict(pairs(parameters)))

    Alg = PathWeightSampling.name(algorithm)
    Duration = duration
    filename = savename((@dict Alg Duration M), "hdf5")
    local_path = datadir("gene_expression", filename)
    mkpath(dirname(local_path))
    h5open(local_path, "w") do file
        PathWeightSampling.write_hdf5!(file, result)
    end
end

# Reduce these numbers to make the computations complete
# faster
N = 10_000
M_vals = [128, 1024, 4096]

@info "Starting DPWS simulations"
for M in M_vals
    save_simulation(:directmc, N, M, parameters)
end

@info "Starting TI-PWS simulations"
for T in [1:10], M in M_vals
    parameters = (;
        mean_s = 50,
        corr_time_s = 1.0,
        corr_time_x = 0.1,
        dtimes = collect(0.0:0.1:T)
    )
    save_simulation(:ti, N, M, parameters)
end

@info "Starting RR-PWS simulations"
for M in [128]
    save_simulation(:smc, N, M, parameters)
end


text:
module Config

using QuantumESPRESSO.Inputs.PWscf: PWInput
using QuantumESPRESSO.Commands: QuantumESPRESSOConfig, PwxConfig
using AbInitioSoftwareBase.Commands: CommandConfig, MpiexecConfig

import Configurations: convert_to_option
import Express.EquationOfStateWorkflow.Config: RuntimeConfig, ExpandConfig

function (::ExpandConfig)(template::AbstractString)
    str = read(expanduser(template), String)
    return parse(PWInput, str)
end

convert_to_option(::Type{RuntimeConfig}, ::Type{CommandConfig}, dict) =
    QuantumESPRESSOConfig(;
        mpi = get(dict, "mpi", MpiexecConfig()),
        pw = get(dict, "pw", PwxConfig()),
    )

end


text:
### ScRNAseq.jl
###
### A julia package for scRNA-seq exploratory data analysis
##
### This file is NOT YET a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

module ScRNAseq

export
    FileIO,
    QualityControl,
    Transformation,
    Normalization,
    Embedding,
    GeneRegulatoryNetwork

    include("helper.jl")
    include("fileio/fileio.jl")
    include("qualitycontrol/qualitycontrol.jl")
    include("transformation/transformation.jl")
    include("normalization/normalization.jl")
    include("embedding/embedding.jl")
    include("differentialexpression/differentialexpression.jl")
    include("generegulatorynetwork/generegulatorynetwork.jl")

using .FileIO
using .QualityControl
using .Transformation
using .Normalization
using .Embedding
using .DifferentialExpression
using .GeneRegulatoryNetwork

end  # module ScRNAseq


text:
module StructOpt

using Printf, DataStructures

using LinearAlgebra, DelimitedFiles, Random, Distributions
using JuMP, Ipopt

## For plotting
using PGFPlotsX, LaTeXStrings, Colors, Contour

export f, g, F, ∇f, prox_αg, get_gradlips, get_μ_cvx, get_optimsol, get_support

export solve_proxgrad, extra_ISTA, extra_FISTA, extra_CondInertia, extra_CondPredInertia, extra_altin_cond, extra_ADAPTsubspace, extra_FISTArestart, extra_FISTArestart_cond, extra_MFISTA, extra_AltIntertia

export LassoPb, LogisticPb
export regularizer_l1, regularizer_l12, regularizer_TV, regularizer_lnuclear, regularizer_distball

export get_randomlasso, get_IM_testcase, get_randomlassopb_TV, build_2dlasso_normdsol, get_ionopb, get_logit

export make_binbounds, make_bins, add_identif_point!, build_save_suboptimality, build_save_steplength, add_contour!, add_traj!, add_point!, pgf_build_iteratesfig, add_manifold!, get_alg_color, get_traj_params, get_legend, runexpnums

include("TV_solver.jl")

###############################################################################
### Problem abstract structure

abstract type AbstractRegularizer end
abstract type AbstractProblem{T} end

# function
f(pb::AbstractProblem, x) = error("NotImplementedError: f for type $(typeof(pb)).")
g(pb::AbstractProblem{T}, x) where T = error("NotImplementedError: g for problem type $T.")

F(pb::AbstractProblem{T}, x) where {T <: AbstractRegularizer} = f(pb, x) + g(pb, x)

# 1st order
∇f(pb::AbstractProblem, x) = error("NotImplementedError: ∇f for type $(typeof(pb)).")
prox_αg(pb::AbstractProblem{T}, x, α) where {T<:AbstractRegularizer} = error("NotImplementedError: prox_αg for problem type $T.")

# conditioning
get_gradlips(pb::AbstractProblem) = error("NotImplementedError: get_gradlips for type $(typeof(pb)).")
get_μ_cvx(pb::AbstractProblem) = -Inf

###############################################################################
## LASSO Problem parametrized by regularizer reg
#    min_x  0.5 * ||Ax-y||² + λ reg(x)

mutable struct LassoPb{T} <: AbstractProblem{T}
    A::Matrix{Float64}
    y::Vector{Float64}
    λ::Float64
    n::Int
    x0::Vector{Float64}
end

get_gradlips(pb::LassoPb) = opnorm(pb.A)^2
get_μ_cvx(pb::LassoPb) = (svdvals(pb.A)[end])^2

f(pb::LassoPb, x) = 0.5 * norm(pb.A*x-pb.y)^2
∇f(pb::LassoPb, x) = transpose(pb.A)*(pb.A*x - pb.y)

function get_optimsol(pb::LassoPb{T}) where {T <: AbstractRegularizer}
    model = Model(with_optimizer(Ipopt.Optimizer, tol=1e-13, print_level=0))
    n = pb.n

    t = @variable model t
    x = @variable model x[1:n]

    p = build_reg!(model, T, x, n)

    @objective model Min 0.5 * transpose(pb.A*x-pb.y) * (pb.A*x-pb.y) + pb.λ * t
    @constraint model t >= sum(p)

    optimize!(model)

    objopt = objective_value(model)
    xopt = value.(x)

    return objopt, xopt
end

###############################################################################
## Logistic Problem
#    min_x  1/m * ∑_{i=1}^m log(1 + exp(-y_i dot(A_i, x))) + λ reg(x)
#
#    x : R^n
#    y : R^m (observations)
#    A : mxn matrix of samples

mutable struct LogisticPb{T} <: AbstractProblem{T}
    A::Matrix{Float64}
    y::Vector{Float64}
    λ::Float64
    n::Int
    x0::Vector{Float64}
end

function f(pb::LogisticPb, x)
    m = size(pb.A, 1)

    return (1/m) * sum(log(1 + exp(-pb.y[i] * dot(pb.A[i, :], x))) for i in 1:m)
end

σ(x) = 1/(1+exp(-x))
function ∇f(pb::LogisticPb, x)
    m = size(pb.A, 1)

    vec = zeros(m)
    for i in 1:m
        vec[i] = (1/m) * σ(-pb.y[i] * dot(pb.A[i, :], x)) * (-pb.y[i])
    end
    return transpose(pb.A) * vec
end

function get_gradlips(pb::LogisticPb)
    m = size(pb.A, 1)
    return opnorm(pb.A)^2 / m
end

function get_optimsol(pb::LogisticPb)
    model = Model(with_optimizer(Ipopt.Optimizer, tol=1e-13, print_level=0))
    n = pb.n
    m = size(pb.A, 1)

    t = @variable model t
    p = @variable model p[1:n] >= 0
    x = @variable model x[1:n]

    ## Linear part of objective...
    lins = @variable model lins[1:m]

    @NLobjective model Min (1/m) * sum(log(1 + exp(-pb.y[i] * lins[i])) for i in 1:m) + pb.λ * t

    my_expr = pb.A * x
    @constraint(model, [i=1:m], lins[i] == dot(pb.A[i, :], x))

    ## l1 norm modeling
    @constraint model t >= sum(p)
    @constraint(model, [i = 1:n], 0 <= x[i]+p[i])
    @constraint(model, [i = 1:n], 0 <= -x[i]+p[i])

    optimize!(model)

    objopt = objective_value(model)
    xopt = value.(x)

    return F(pb, xopt), xopt
end

include("alg_proxgrad.jl")
include("algs_extrapolation.jl")
include("regularizers.jl")
include("build_problems.jl")

include("runexpnums.jl")
include("plot_utils.jl")

end # module


text:
function apsf(::Type{MethodSincR}, sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)
    check_amp_sampling(sz, pp, sampling)

    if isnothing(sampling)
        sampling = get_required_amp_sampling(sz, pp) # this is knowingly too small along kz but fixed down below, if needed.
    end
    Ewald_sampling = get_Ewald_sampling(sz, pp)
    undersampling_factor = sampling[3] ./ Ewald_sampling[3]
    # The big_sz is used to avoid interference effects by the application of the multiplication with the pupil maks in Fourier space
    # The size if 4/3 is such that the diagonals are covered but a wrap-around into the original size is avoided.
    # big_sz = ((sz[1:2].*2)...,sz[3]+2)  # The 2 extra slices in z seem to be necessary to avoid a problem at the fist position
    big_sz = ((sz[1:2].*2)...,sz[3]+4)  # The 2 extra slices in z seem to be necessary to avoid a problem at the fist positions

    if undersampling_factor > 1.0 # the current sampling is insufficient for the SincR method. We need to properly sample the Ewald sphere.
        # the size sz[3] will be cut out in Fourier-spage around the cleaned and processed part of the McCutchen pupil, which should yield the sampling as currently specified.
        nowrap_sz = (big_sz[1:2]..., ceil(Int, undersampling_factor .* big_sz[3]))
        big_sampling = (sampling[1:2]..., sampling[3] * big_sz[3] / nowrap_sz[3]) # calculate the sampling needed to cover the Ewald sphere.
    else
        nowrap_sz = big_sz
        big_sampling = sampling
    end

    pupil = pupil_xyz(nowrap_sz, pp, big_sampling) # field_xyz(big_sz,pp, sampling) .* aplanatic_factor(big_sz,pp,sampling) .* ft(jinc_r_2d(big_sz[1:2],pp, sampling=sampling) .* my_disc(big_sz[1:2],pp)) #
    sinc_r_big = sinc_r(nowrap_sz,pp, sampling=big_sampling) .* my_disc(nowrap_sz[1:2],pp)  # maybe this should rather already be apodized by angle?
    if pp.FFTPlan != nothing
        P3d = plan_fft(sinc_r_big, flags=pp.FFTPlan)
    end

    if big_sampling[3] != sampling[3] # we need to extract (to reduce the size to the user-defined) and circshift (to get the correct phases)
        kzc, rel_kz = get_McCutchen_kz_center(nowrap_sz,pp,big_sampling)
        shell = select_region(theta_z(nowrap_sz) .* ft(sinc_r_big, (1,2,3)), new_size=(nowrap_sz[1:2]...,big_sz[3]), center = kzc) # centers the Pupil along kz
        if !center_kz
            shell = circshift(shell, (0,0,rel_kz)) #
        end
    else
        shell =  theta_z(nowrap_sz) .* ft(sinc_r_big, (1,2,3))
        if center_kz
            kzc, rel_kz = get_McCutchen_kz_center(nowrap_sz,pp,big_sampling)
            shell = circshift(shell, (0,0,-rel_kz)) # needed for subsequent upsampling along kz
        end
    end
    # shell, sampling =  limit_kz(theta_z(nowrap_sz) .* ft(sinc_r_big, (1,2,3)), pp, sampling) # remove negative frequencies and limit to useful range

    # check_amp_sampling_sincr(nowrap_sz, pp, sampling)
    # print("Amplitude sampling is $sampling \n")
    pupils = (cos.(pupil_θ(nowrap_sz,pp,sampling)) .* pupil) .* iftz(shell)
    if pp.FFTPlan != nothing
        Pm2d = plan_fft(pupils,(1,2), flags=pp.FFTPlan)
    end

    res=ift2d(pupils) # This should really be a zoomed iFFT
    return select_region(res, new_size=sz[1:3]) # extract the central bit, which avoids the wrap-around effects
    # , sampling
end

# This function is needed for the propagate method
"""
    apply_propagators(pupil, z_planes, pp::PSFParams; sampling=nothing)

propagates a given `pupil` by a number of `z_planes` (almost) symmetrically in both directions.
The result is a stack of propagated pupils. The slice-to-slice propagator is obtained via the
`get_propagator` method.

See also:
+ get_propagator

"""
function apply_propagators(pupil, z_planes, pp::PSFParams; sampling=nothing)
    sz = (size(pupil)[1:2]...,z_planes,size(pupil)[4])
    # calculate the phase derivatives
    pupils = Array{Complex{pp.dtype}}(undef, sz)
    start_z = -z_planes÷2-1
    prop_kz , scalar, xy_scale = get_propagator(sz, pp, sampling)
    # dx_prop_phase, dy_prop_phase = get_propagator_gradient(prop_phase, scalar, xy_scale)
    for z in 1:sz[3]
        z_pos = z+start_z # reaches zero at center
        # z_pos_in_λ =  z_pos * sampling[3] / (pp.λ / pp.n)
        # pupils[:,:,z:z,:] .= pupil .* (prop.^z_pos_in_λ)
        # pupils[:,:,z:z,:] .= pupil .* cis.(z_pos .* prop_phase)

        # the integral of exp(i z prop(x)) over x is
        # exp(i z prop(x))   (-2/(z^2) + i 2 prop(x)/z)
        pupils[:,:,z:z,:] .= pupil .* cis.(z_pos .* prop_kz) # .*  # cis means: exp.(i * z_pos .* prop_phase)
             # (sinc.(z_pos .* dx_prop_phase) .* sinc.(z_pos .* dy_prop_phase))  # This dampling accounts for the aliasing problem
    end
    return pupils
end

function apsf(::Type{MethodPropagate}, sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)
    if isnothing(sampling)
        sampling = get_Ewald_sampling(sz, pp)
        print("Sampling is $sampling \n")
    end
    check_amp_sampling(sz, pp, sampling)
    pupil = pupil_xyz(sz, pp, sampling) # field_xyz(big_sz,pp, sampling) .* aplanatic_factor(big_sz,pp,sampling) .* ft(jinc_r_2d(big_sz[1:2],pp, sampling=sampling) .* my_disc(big_sz[1:2],pp)) #
    pupils = apply_propagators(pupil, sz[3], pp, sampling=sampling)
    # return pupils
    # pupils .*= window_hanning((1,1,size(pupils,3)),border_in=0.8,border_out=1.0,dims=(3,))
    if center_kz
        _, rel_kz = get_McCutchen_kz_center(sz,pp,sampling)
        pupils .*= cispi.((-2*rel_kz/sz[3]) .* zz((1,1,sz[3]))) # centers the McCutchen pupil to be able to correctly resample it along kz
    end
    if pp.FFTPlan != nothing
        Pm2d = plan_fft(pupils,(1,2), flags=pp.FFTPlan)
    end
    res=ift2d(pupils) # This should really be a zoomed iFFT
    return res # extract the central bit, which avoids the wrap-around effects
end

# This function is needed for the propagate method, but iterates between real and Fourier space always smoothly deleting "out-of-bound" waves.
# The final result is already in real space.
function apply_propagator_iteratively(sz, pp::PSFParams; sampling=nothing, center_kz=false)
    z_planes = sz[3]
    # calculate the phase derivatives
    start_z = z_planes÷2+1
    max_pix_travel = (tan(asin(pp.NA / pp.n)) * sampling[3]) ./ sampling[1:2] # how much does the maximal anlge travel geometrically
    PMLsz = (ceil.(Int, max_pix_travel .* 8)...,0) # To get the number of perfectly matched layer (PML) pixels to append
    psz = sz .+ 2 .*PMLsz # total size for propagation pupil
    pupil = pupil_xyz(psz[1:2], pp, sampling) # field_xyz(big_sz,pp, sampling) .* aplanatic_factor(big_sz,pp,sampling) .* ft(jinc_r_2d(big_sz[1:2],pp, sampling=sampling) .* my_disc(big_sz[1:2],pp)) #
    sz = (sz[1:3]...,size(pupil)[4])
    slices = Array{Complex{pp.dtype}}(undef, sz)

    prop_kz, scalar, xy_scale = get_propagator(psz, pp, sampling)
    start_pupil = pupil # copy(pupil)
    prop_pupil = cis.(prop_kz)
    if center_kz
        _, rel_kz = get_McCutchen_kz_center(psz[1:3],pp,sampling)
        prop_pupil .*= cispi(-2*rel_kz/psz[3])
    end
    border_in = 1.0 .- PMLsz[1:2] ./ psz[1:2]
    real_window = exp.(1.75 .* (window_linear(size(start_pupil)[1:2],border_in=border_in,border_out=1) .-1))  # This is maybe not the best PML?
    # real_window = window_gaussian(size(start_pupil)[1:2], border_in=border_in,border_out=border_in.*0.5 .+ 0.5)  # This is maybe not the best PML?
    # real_window = window_hanning(size(start_pupil)[1:2],border_in=border_in,border_out=1)
    prop_pupil = conj(prop_pupil) # from now the advancement is in the opposite direction
    if pp.FFTPlan !== nothing
        P2d = plan_fft(pupil,(1,2), flags=pp.FFTPlan)
    end
    if pp.FFTPlan !== nothing
        Pi2d = plan_ifft(pupil,(1,2), flags=pp.FFTPlan)
    end
    pupil = start_pupil
    for z in start_z:-1:2 # from the middle to the start
        slice = ift2d(pupil) # should be ifft2d for speed reasons
        dst = @view slices[:,:,z:z,:]
        select_region!(slice, dst, new_size=sz[1:2])
        pupil = ft2d(slice .* real_window)
        pupil .*= prop_pupil
    end
    dst = @view slices[:,:,1:1,:]
    select_region!(ift2d(pupil), dst, new_size=sz[1:2])

    if has_z_symmetry(pp) # to save some speed
        dz = sz[3] - (start_z+1)
        slices[:,:,start_z+1:start_z+1+dz,:] .= conj(slices[:,:,start_z-1:-1:start_z-1-dz,:]);
    else
        prop_pupil = conj(prop_pupil) # from now the advancement is in the opposite direction
        pupil = start_pupil .* prop_pupil
        for z in start_z+1:sz[3]-1  # from the middle forward
            slice = ift2d(pupil)
            dst = @view slices[:,:,z:z,:]
            select_region!(slice, dst, new_size=sz[1:2])
            pupil = ft2d(slice .* real_window)
            pupil .*= prop_pupil
        end
        dst = @view slices[:,:,sz[3]:sz[3],:]
        select_region!(ift2d(pupil), dst, new_size=sz[1:2])
    end
    return slices
end

function apsf(::Type{MethodPropagateIterative}, sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)
    if isnothing(sampling)
        sampling = get_Ewald_sampling(sz, pp)
        print("Sampling is $sampling \n")
    end
    check_amp_sampling(sz, pp, sampling)
    slices = apply_propagator_iteratively(sz, pp, sampling=sampling, center_kz=center_kz)
    return slices # extract the central bit, which avoids the wrap-around effects
end

# just a different way of writing the propagation down
function apply_kz_to_pupil(pupil, z_planes, pp; sampling)
    sz = (size(pupil)[1:2]...,z_planes,size(pupil)[4])
    if isnothing(sampling)
        sampling = get_Ewald_sampling(sz, pp)
    end
    k_max_rel = sampling[1:2] ./ (pp.λ / pp.n)
    scalar = (2π*sampling[3] / (pp.λ / pp.n))
    xy_scale = 1 ./ (k_max_rel .* sz[1:2])
    sqrt_term = phase_kz(pp.dtype, sz[1:2], scale = xy_scale)
    prop_phase = scalar .* sqrt_term
    pupils = pupil .* cis.(zz((1,1,z_planes)) .* prop_phase) # if you fourier-transform this along kz you should get the shifted pupil
    return pupils
end

# Basic idea of MethodShell: Displace each kxy by the appropriate kz
function apsf(::Type{MethodShell}, sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)
    if isnothing(sampling)
        sampling = get_Ewald_sampling(sz, pp)
        print("Sampling is $sampling \n")
    end
    check_amp_sampling(sz, pp, sampling)
    pupil = pupil_xyz(sz, pp, sampling) # field_xyz(big_sz,pp, sampling) .* aplanatic_factor(big_sz,pp,sampling) .* ft(jinc_r_2d(big_sz[1:2],pp, sampling=sampling) .* my_disc(big_sz[1:2],pp)) #

    pupils = apply_kz_to_pupil(pupil, sz[3], pp, sampling=sampling)
    res=ift2d(pupils) # This should really be a zoomed iFFT
    return res # extract the central bit, which avoids the wrap-around effects
end

function simpson!(f,α,N=80)
    h = α / N # N+1 points are really used...
    I012 = h/6*f(0) # middle values count 4 times
    for n in 0:N-1 # integrate using Simpsons rule
        theta = h*(n+0.5) # middle value(s)
        I012 += 4*h/6*f(theta) # middle values count 4 times
        if n<(N-1)
            theta = h*(n+1)
            I012 += 2*h/6*f(theta) # intermediate values (count twice in the sum)
        end
    end
    I012 += h/6*f(α) # last value
    return I012
end

# Calculates I0, I1 and I2 according to the Richards and Wolf paper
# The calculation is done on an RZ plane and then interpolated to a 3D volume.
# This is probably still buggy! At least the Pupil aplanatic factor looks wrong!
function apsf(::Type{MethodRichardsWolf}, sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)
    sampling_r = min(sampling[1],sampling[2])/2.0
    diagonal = sqrt(sum(abs2.(sz[1:2] .* sampling[1:2]))) / 2.0;
    sr = ceil(Int64, diagonal/sampling_r) + 1
    r = xx((sr,), scale=sampling_r, offset=CtrCorner) # a generator for radius
    k = 2pi / (pp.λ / pp.n);
    kz = yy((1,sz[3]), scale=k*sampling[3])
    function integrate!(I, w, theta)  # w includes the aplanatic factor
        (sint,cost) = sincos(theta);
        krsint =  (k*sint).*r
        ciskzcost =  cis.(kz.*cost)
        I[:,:,1] .+= (w* sint*(1+cost)).*besselj0.(krsint) .* ciskzcost
        I[:,:,2] .+= (w* sint^2) .*besselj1.(krsint) .* ciskzcost
        I[:,:,3] .+= (w* sint.*(1-cost)).*besselj.(2,krsint) .* ciskzcost
    end
    α = asin(pp.NA / pp.n); # maximal aperture angle
    # Now perform the integration according to Simpsons rule
    N = 50
    h = α / N # since these are really N+1 points
    I012 = zeros(Complex{pp.dtype}, (sr,sz[3],3)) # I0, I1 and I2
    integrate!(I012, h/6*pp.aplanatic(0), 0) # first value
    for n in 0:N-1 # integrate using Simpsons rule
        theta = h*(n+0.5) # middle value(s)
        integrate!(I012, 4*h/6*pp.aplanatic(theta), theta) # middle values count 4 times
        if n<(N-1)
            theta = h*(n+1)
            integrate!(I012, 2*h/6*pp.aplanatic(theta), theta) # intermediate values (count twice in the sum)
        end
    end
    integrate!(I012, h/6*pp.aplanatic(α), α) # last value

    # Interpolate the RZ-results onto the 3D grid
    phi = phiphi((sz[1],sz[2]), scale=(sampling[1],sampling[2]))
    sinphi = sin.(phi)
    cosphi = cos.(phi)
    cos2phi = cos.(2 .*phi)
    sin2phi = sin.(2 .*phi)
    rpos = rr((sz[1],sz[2]), scale=(sampling[1],sampling[2]))
    r_idx = 1 .+ floor.(Int64, rpos ./ sampling_r) # index position
    w = 1.0 .- (rpos/ sampling_r .+ 1 .- r_idx) # index position
    E = zeros(Complex{pp.dtype}, (sz...,3)) # I0, I1 and I2
    for z = 1:sz[3]
        I0 = @view I012[:,z,1] # create views to be able to write by 1D indexing into the appropriate slices.
        I1 = @view I012[:,z,2]
        I2 = @view I012[:,z,3]

        if pp.polarization == pol_x
            E[:,:,z,1] .= -im.*((w.*I0[r_idx].+(1 .-w).*I0[r_idx.+1]).+cos2phi.*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1]))
            E[:,:,z,2] .= -im.*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1]).*sin2phi
            E[:,:,z,3] .= -2 .*(w.*I1[r_idx].+(1 .-w).*I1[r_idx.+1]).*cosphi
        elseif pp.polarization == pol_y
            E[:,:,z,1] .= im.*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1]).*sin2phi # sin2phi with pi/2 phase shift becomes -sin2phi
            E[:,:,z,2] .= -im.*((w.*I0[r_idx].+(1 .-w).*I0[r_idx.+1]).-cos2phi.*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1])) # cos2phi becomes -cos2phi
            E[:,:,z,3] .= 2 .*(w.*I1[r_idx].+(1 .-w).*I1[r_idx.+1]).*sinphi # cosphi becomes -sinphi
        elseif pp.polarization == pol_circ # pol_x + im* pol_y
            E[:,:,z,1] .= -im/sqrt(2).*((w.*I0[r_idx].+(1 .-w).*I0[r_idx.+1]).+cos2phi.*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1])) .-
            1/sqrt(2)*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1]).*sin2phi
            E[:,:,z,2] .= -im/sqrt(2).*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1]).*sin2phi .+
            1/sqrt(2)*((w.*I0[r_idx].+(1 .-w).*I0[r_idx.+1]).-cos2phi.*(w.*I2[r_idx].+(1 .-w).*I2[r_idx.+1]))
            E[:,:,z,3] .= 1/sqrt(2)*(w.*I1[r_idx].+(1 .-w).*I1[r_idx.+1]).*(-2 .*cosphi .+ 2im * sinphi)
        else
            error("unsupported polarization for Richards-Wolf method")
        end
    end
    return E
end

"""
    apsf(sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)

dispatches to various amplitude point spread function calculation routines.
Note that the `method` entry in `pp` defines which calculation method to be used.
Alternatively a different method can be chosen like this: `apsf(::Type{MethodShell}, sz, ..)`.

Arguments:
+ sz: NTuple of size to generate
+ pp: PSF parameter structure, also contains the dtype. See PSFParam() for details
+ sampling: NTuple for pixel pitch information
+ center_kz: if true, the McCutchen pupil will be centered along the kz direction. This is important to be able to apply a consecutive resampling without errors.
             However, the phase values are then not correct, which does not matter for intensity PSFs.

See also:
+ psf():    calculates the intensity point spread function (psf) by taking (sum along field components of) the absolute square of the corresponding apsf.

Example:
```jdoctest
julia> pp = PSFParams(500.0,1.4,1.52)
julia> p = apsf((128,128,128),pp; sampling=(50,50,100));
```
"""
function apsf(sz::NTuple, pp::PSFParams; sampling=nothing, center_kz=false)
    apsf(pp.method, sz, pp, sampling=sampling, center_kz=center_kz)
end


text:
using HyperGraphs

# adding and deleting vertices -- unoriented hypergraph
x = HyperGraph([HyperEdge([1, 2, 3]), HyperEdge([2, 3, 4])])
@test add_vertex!(x, 5, check=false) && has_vertex(x, 5)
@test_throws ErrorException add_vertex!(x, 5) # vertex already in hypergraph
@test_throws ErrorException add_vertex!(x, 5.0) # type error
@test add_vertices!(x, [6, 7]) && has_vertices(x, [6, 7])
@test del_vertex!(x , 1, weak=true) && has_hyperedge(x, HyperEdge([2, 3])) # weak deletion of vertex retains hyperedges
@test del_vertex!(x, 4) && !has_hyperedge(x, HyperEdge([2, 3, 4])) # strong deletion of vertex removes hyperedges
@test del_vertices!(x, [2, 3], weak=true) && !any(has_vertex(x, [2, 3])) && has_empty_hyperedges(x) # removing all vertices of a hyperedge weakly results in an empty hyperedge
@test del_vertices!(x, vertices(x)) && !any(has_vertex(x, [5, 6, 7])) && has_hyperedge(x, HyperEdge(Int[])) # only an empty hyperedge left
@test del_empty_hyperedges!(x) && !has_empty_hyperedges(x) && isempty(x)
@test add_vertices!(x, collect(1:4))
@test del_vertices!(x, [1]) && add_vertices!(x, [1]) # testing vectors with length 1
@test add_hyperedge!(x, HyperEdge(vertices(x)), check=false) && has_hyperedge(x, HyperEdge(vertices(x))) # check is true by default but not needed here
@test_throws ErrorException add_hyperedge!(x, HyperEdge([10])) # can't add hyperedge if vertex not already in hypergraph
@test add_hyperedges!(x, [HyperEdge([1 , 2, 3]), HyperEdge([1, 2, 3])]) && has_parallel_hyperedges(x) # add parallel hyperedge
@test del_hyperedge!(x, hyperedges(x)[end]) && !has_parallel_hyperedges(x) # removed parallel hyperedges
@test del_hyperedges!(x, hyperedges(x)) && isempty(x)
@test del_vertices!(x, vertices(x)) && nv(x) == 0

# removing parallel hyperedges in hypergraph with only parallel hyperedges -- unoriented hypergraph
x = HyperGraph([HyperEdge([1, 2, 3]), HyperEdge([1, 2, 3])])
@test del_hyperedges!(x, hyperedges(x)) && isempty(x)

# dealing with loops -- unoriented hypergraph
x = HyperGraph{Int}([1, 2, 3], HyperEdge[]) # initialise an empty hypergraph with some vertices
@test add_hyperedge!(x, HyperEdge([1, 1])) && has_loops(x) && num_loops(x) == 1
@test add_hyperedges!(x, [HyperEdge([2, 2, 2]), HyperEdge([3, 3])]) && num_loops(x) == 3
@test del_hyperedge!(x, HyperEdge([1, 1])) && num_loops(x) == 2
@test del_hyperedges!(x, hyperedges(x)) && !has_loops(x) && isempty(x)

# adding and deleting vertices -- oriented hypergraph
chx = ChemicalHyperGraph([ChemicalHyperEdge([1], [2, 3]), ChemicalHyperEdge([2], [3, 4])])
@test add_vertex!(chx, 5, check=false) && has_vertex(chx, 5)
@test_throws ErrorException add_vertex!(chx, 5) # vertex already in hypergraph
@test_throws ErrorException add_vertex!(chx, 5.0) # type error (might be bad for later developments)
@test add_vertices!(chx, [6, 7]) && has_vertices(chx, [6, 7])
@test del_vertex!(chx, 1, weak=true) && has_hyperedge(chx, ChemicalHyperEdge(Int[], [2, 3])) # weak deletion of vertex retains hyperedges
@test del_vertex!(chx, 4) && !has_hyperedge(chx, ChemicalHyperEdge([2], [3, 4])) # strong deletion of vertex removes hyperedges
@test del_vertices!(chx, [2, 3], weak=true) && !any(has_vertex(chx, [2, 3])) && has_empty_hyperedges(chx) # removing all vertices of a hyperedge weakly results in an empty hyperedge
@test del_vertices!(chx, vertices(chx)) && !any(has_vertex(chx, [5, 6, 7])) && has_hyperedge(chx, ChemicalHyperEdge{Int}()) # only an empty hyperedge left
@test del_empty_hyperedges!(chx) && !has_empty_hyperedges(chx) && isempty(chx)
@test add_vertices!(chx, collect(1:4))
@test del_vertices!(chx, [1]) && add_vertices!(chx, [1]) # testing vectors with length 1
@test add_hyperedge!(chx, ChemicalHyperEdge(vertices(chx), vertices(chx)), check=false) && has_hyperedge(chx, ChemicalHyperEdge(vertices(chx), vertices(chx))) # check is true by default but not needed here
@test_throws ErrorException add_hyperedge!(chx, ChemicalHyperEdge([10], [11])) # can't add hyperedge if vertex not already in hypergraph
@test add_hyperedges!(chx, [ChemicalHyperEdge([1], [2, 3]), ChemicalHyperEdge([1], [2, 3])]) && has_parallel_hyperedges(chx) # add parallel hyperedge
@test del_hyperedge!(chx, hyperedges(chx)[end]) && !has_parallel_hyperedges(chx) # removed parallel hyperedges
@test del_hyperedges!(chx, hyperedges(chx)) && isempty(chx)
@test del_vertices!(chx, vertices(chx)) && nv(chx) == 0

# removing parallel hyperedges in hypergraph with only parallel hyperedges -- oriented hypergraph
chx = ChemicalHyperGraph([ChemicalHyperEdge([1], [2, 3]), ChemicalHyperEdge([1], [2, 3])])
@test del_hyperedges!(chx, hyperedges(chx)) && isempty(chx)

# dealing with loops - oriented hypergraph
chx = ChemicalHyperGraph{Int}([1, 2, 3], ChemicalHyperEdge[]) # initialise an empty hypergraph with some vertices
@test add_hyperedge!(chx, ChemicalHyperEdge([1], [1])) && has_loops(chx) && num_loops(chx) == 1
@test add_hyperedges!(chx, [ChemicalHyperEdge([2], [2, 2]), ChemicalHyperEdge([3], [3])]) && num_loops(chx) == 3
@test del_hyperedge!(chx, ChemicalHyperEdge([1], [1])) && num_loops(chx) == 2
@test del_hyperedges!(chx, hyperedges(chx)) && !has_loops(chx) && isempty(chx)

# modifying hypergraphs of type Any
x = HyperGraph{Any}(); v = 1; e = HyperEdge(Any[v, v])
@test add_vertex!(x, v)
@test add_hyperedge!(x, e)
@test vertices(x) == [v] && nv(x) == 1
@test hyperedges(x) == [e] && nhe(x) == 1
chx = ChemicalHyperGraph{Any}(); v = 1; che = ChemicalHyperEdge(Any[v], Any[v])
@test add_vertex!(chx, v)
@test add_hyperedge!(chx, che)
@test vertices(chx) == [v] && nv(chx) == 1
@test hyperedges(chx) == [che] && nhe(chx) == 1

e = HyperEdge([1, 2, 3])
replace_vertex!(e, (2 => 1))
@test vertices(e) == [1, 1, 3]
che = ChemicalHyperEdge([1, 2, 3], [2, 3, 4])
replace_vertex!(che, (2 => 1))
@test vertices(che) == [1, 1, 3, 1, 3, 4]
x = HyperGraph(e)
replace_vertex!(x, (2 => 1))
@test vertices(x) == [1, 3]
@test vertices(hyperedges(x)) == [[1, 1, 3]]
chx = ChemicalHyperGraph(che)
replace_vertex!(chx, (2 => 1))
@test vertices(x) == [1, 3]
@test vertices(hyperedges(chx)) == [[1, 1, 3, 1, 3, 4]]

# replace hyperedge, replace hyperedges -- unoriented hypergraph
x = HyperGraph([HyperEdge([1]), HyperEdge([2])])
new_e = HyperEdge([1, 2])
@test replace_hyperedge!(x, hyperedges(x)[1], new_e) && has_hyperedge(x, new_e)
new_es = [HyperEdge([1, 1]), HyperEdge([2, 2])]
@test replace_hyperedges!(x, hyperedges(x), new_es) && has_hyperedges(x, new_es)

# replace hyperedge, replace hyperedges -- oriented hypergraph
chx = ChemicalHyperGraph([ChemicalHyperEdge([1], [2]), ChemicalHyperEdge([2], [1])])
new_che = ChemicalHyperEdge([1], [1])
@test replace_hyperedge!(chx, hyperedges(chx)[1], new_che) && has_hyperedge(chx, new_che)
new_ches = [ChemicalHyperEdge([1], [1]), ChemicalHyperEdge([2], [2])]
@test replace_hyperedges!(chx, hyperedges(chx), new_ches) && has_hyperedges(chx, new_ches)

# weight updating
chx = ChemicalHyperGraph([ChemicalHyperEdge(src, tgt) for (src, tgt) in zip([[1], [4]], [[1, 2, 3], [2, 3, 4]])])
che = hyperedges(chx)[1]
@test update_weight!(hyperedges(chx)[1], 2) == 2 && weight(che) == 2
@test update_weight!(chx, che, 3) == 3 && weight(che) == 3
@test set_all_weights!(chx, 4) == repeat([4], nhe(chx))
@test all(weight.(hyperedges(chx)) .== 4)

# hypergraphs merging
e1 = HyperEdge([1, 1])
e2 = HyperEdge([1, 2])
@test merge(HyperGraph(e1), HyperGraph(e2)) == HyperGraph([e1, e2])
che1 = ChemicalHyperEdge([1], [1])
che2 = ChemicalHyperEdge([1], [2])
@test merge(ChemicalHyperGraph(che1), ChemicalHyperGraph(che2)) == ChemicalHyperGraph([che1, che2])

# vertex-induced subhypergraph
x = HyperGraph([HyperEdge([1, 2, 3]), HyperEdge([2, 3, 4])])
@test nv(subhypergraph(x, 1)) == 3
@test subhypergraph(x, 2) == subhypergraph(x, [2]) == x
@test subhypergraph(x, 4) == subhypergraph(x, [4])
@test subhypergraph(x, vertices(x)) == x
chx = ChemicalHyperGraph([ChemicalHyperEdge(src, tgt) for (src, tgt) in zip([[1], [4]], [[1, 2, 3], [2, 3, 4]])])
@test nv(subhypergraph(chx, 1)) == 3
@test subhypergraph(chx, 2) == subhypergraph(chx, [2]) == chx
@test subhypergraph(chx, 4) == subhypergraph(chx, [4])
@test subhypergraph(chx, vertices(chx)) == chx

# hyperedge-induced subhypergraph
@test nhe(subhypergraph(x, hyperedges(x)[1])) == 1
@test nv(subhypergraph(x, HyperEdge([1, 2]), relaxed=true)) == 2
@test_throws ErrorException subhypergraph(x, HyperEdge([2, 3]), relaxed=true) # subhyperedge is not unique
@test subhypergraph(x, hyperedges(x)) == x
@test vertices(subhypergraph(x, [HyperEdge([1, 2]), HyperEdge([3, 4])], relaxed=true)) == vertices(x)
@test nhe(subhypergraph(chx, hyperedges(chx)[1])) == 1
@test nv(subhypergraph(chx, ChemicalHyperEdge([1], [2]), relaxed=true)) == 2
@test_throws ErrorException subhypergraph(chx, ChemicalHyperEdge(Int[], [2, 3]), relaxed=true) # subhyperedge is not unique
@test subhypergraph(chx, hyperedges(chx)) == chx
@test vertices(subhypergraph(chx, [ChemicalHyperEdge(Int[], [1, 2]), ChemicalHyperEdge(Int[], [3, 4])], relaxed=true)) == vertices(chx)

# hyperedge-induced subhypergraph with parallel hyperedges
e = HyperEdge([1, 2, 3])
parallel_x = HyperGraph([e, e])
@test subhypergraph(parallel_x, e) == parallel_x
che = ChemicalHyperEdge([1, 2], [3])
parallel_chx = ChemicalHyperGraph([che, che])
@test subhypergraph(parallel_chx, che) == parallel_chx

# edge switching
che = ChemicalHyperEdge([1, 2], [3], 5)
switched_che = switch(che)
@test src(che) == tgt(switched_che) && tgt(che) == src(switched_che)
@test weight(che) == weight(switched_che)


text:
# Conversion between model types

"""
    LinearLayeredModel(m::PREMPolyModel, spacing=20)

Convert the model `m` into a `LinearLayeredModel`, with radial
knots spaced a maximum of `spacing` km apart.
"""
function LinearLayeredModel(m::PREMPolyModel, spacing=20)
    R = surface_radius(m)
    # Radii for evaluation of model points at knots
    r_eval = Float64[]
    # Radii of knots, actually saved in new model
    radii = Float64[]
    for (ilayer, rtop) in enumerate(m.r)
        rtop_minus = prevfloat(rtop)
        # Bottom of layer
        rbot = ilayer == 1 ? 0.0 : m.r[ilayer - 1]
        rbot_plus = ilayer == 1 ? rbot : nextfloat(rbot)
        push!(radii, rbot)
        push!(r_eval, rbot_plus)
        r = rbot + spacing
        # Add knots
        while r < rtop
            push!(r_eval, r)
            push!(radii, r)
            r += spacing
        end
        # Top of layer
        push!(radii, rtop)
        push!(r_eval, rtop == R ? rtop : rtop_minus)
    end
    Vp, Vs = vp.(m, r_eval), vs.(m, r_eval)
    Rho = hasdensity(m) ? density.(m, r_eval) : []
    Vpv, Vph, Vsv, Vsh, Eta = isanisotropic(m) ?
        (vpv.(m, r_eval), vph.(m, r_eval), vsv.(m, r_eval), vsh.(m, r_eval), eta.(m, r_eval)) :
        ([], [], [], [], [])
    qκ, qμ = hasattenuation(m) ? (Qκ.(m, r_eval), Qμ.(m, r_eval)) : ([], [])
    LinearLayeredModel(R, length(radii), radii, Vp, Vs, Rho,
        isanisotropic(m), Vph, Vpv, Vsh, Vsv, Eta,
        hasattenuation(m), qμ, qκ)
end

"Set of fields we copy to the LinearLayeredModel struct"
const LinearLayeredModel_CONVERSION_FIELDS =
    (x for x in fieldnames(LinearLayeredModel)
     if x ∉ (:a, :n, :r, :aniso, :attenuation))

"""
    LinearLayeredModel(m::SteppedLayeredModel)

Convert the model `m` into a `LinearLayeredModel`, with radial
knots spaced at the bottom and top of each layer in `m`.
"""
function LinearLayeredModel(m::SteppedLayeredModel)
    R = surface_radius(m)
    d = Dict{Symbol,Vector{Float64}}()
    radii = Iterators.flatten((i == 1 ? 0.0 : m.r[i-1], m.r[i]) for i in eachindex(m.r))
    for field in LinearLayeredModel_CONVERSION_FIELDS
        val = getfield(m, field)
        isempty(val) && continue
        d[field] = repeat(val, inner=2)
    end
    d[:r] = collect(radii)
    LinearLayeredModel(; d...)
end

"Set of fields we copy to the PREMPolyModel struct"
const PREMPolyModel_CONVERSION_FIELDS = (x for x in fieldnames(PREMPolyModel)
                                         if x ∉ (:r, :a, :n, :fref, :aniso, :attenuation))

"""
    PREMPolyModel(m::SteppedLayeredModel, order=0; fref=1.0)

Convert the model `m` into a `PREMPolyModel`, with maximum polynomial
`order`.  The reference frequency `fref` in Hz can be given.
"""
function PREMPolyModel(m::SteppedLayeredModel, order=0; fref=1.0)
    order < 0 &&
        throw(ArgumentError("cannot construct polynomials with negative order"))
    nlayers = length(m.r)
    d = Dict{Symbol,Any}()
    for field in PREMPolyModel_CONVERSION_FIELDS
        val = getfield(m, field)
        isempty(val) && continue
        if val isa AbstractArray && order > 0
            arr = fill(0.0, order + 1, nlayers)
            arr[1,:] .= val
            d[field] = arr
        else
            d[field] = val
        end
    end
    d[:fref] = fref
    d[:r] = m.r
    PREMPolyModel(; d...)
end

"""
    PREMPolyModel(m::LinearLayeredModel, order=1; fref=1.0)

Convert the model `m` into a `PREMPolyModel`, with maximum polynomial
`order`.  The reference frequency `fref` in Hz can be given.
"""
function PREMPolyModel(m::LinearLayeredModel, order=1; fref=1.0)
    order < 1 &&
        throw(ArgumentError("minimum polynomial order is 1 when " *
                            "converting `LinearLayeredModel`s"))
    d = Dict{Symbol,Any}()
    R = surface_radius(m)
    # Indices of nodes on top of finite layers (i.e., not discontinuity layers)
    layer_inds = [i for i in 2:length(m.r) if m.r[i] > m.r[i-1]]
    nlayers = length(layer_inds)
    for field in PREMPolyModel_CONVERSION_FIELDS
        val = getfield(m, field)
        isempty(val) && continue
        if val isa AbstractArray
            d[field] = fill(0.0, order + 1, nlayers)
            for (i, ind) in enumerate(layer_inds)
                d[field][1:2,i] .= _fit_prem_polynomial(R,
                    m.r[ind-1], m.r[ind], val[ind-1], val[ind])
            end
        else
            d[field] = val
        end
    end
    d[:fref] = fref
    d[:r] = m.r[layer_inds]
    PREMPolyModel(; d...)
end

"""
    _fit_prem_polynomial(a, r1, r2, v1, v2) -> (a0, a1)

Return PREM polynomial coefficients (where `a0` is the constant and `a1`
the linear coefficient) for the line connecting two values `v1` and
`v2` at radii `r1` and `r2`.  The radii are unnormalised compared to
the planetary radius `a`.
"""
function _fit_prem_polynomial(a, r1, r2, v1, v2)
    x1, x2 = r1/a, r2/a
    a1 = (v2 - v1)/(x2 - x1)
    a0 = v1 - a1*x1
    a0, a1
end

"Set of fields we copy into the SteppedLayeredModel struct"
const SteppedLayeredModel_CONVERSION_FIELDS =
    (x for x in fieldnames(SteppedLayeredModel)
     if x ∉ (:r, :a, :n, :fref, :aniso, :attenuation))

"""
    SteppedLayeredModel(m::Union{LinearLayeredModel, PREMPolyModel}, spacing=10)

Convert the model `m` into a `SteppedLayeredModel`, with layer tops
spaced a maximum of `spacing` km apart.  Where layer tops of `m` lie
between the regular set of `spacing` layer tops, extra layers are
inserted.  Thus layer tops in the new model always coincide with
discontinuities in `m`.

Layers take the value of the properties at the midpoint of the
`spacing` km interval in the original `PREMPolyModel`.
"""
function SteppedLayeredModel(m::Union{LinearLayeredModel, PREMPolyModel}, spacing=10)
    spacing <= 0 && throw(ArgumentError("layer spacing must be positive"))
    d = Dict{Symbol,Any}()
    R = surface_radius(m)
    # Regular set of radii spaced `spacing` km apart
    spaced_radii = reverse(range(R, stop=spacing, step=-spacing))
    # New model layer top radii with m's layers interspersed.
    radii = unique!(sort!(vcat(m.r, spaced_radii)))
    # Radii from which to take layer properties in the new model
    eval_radii = [((i == 1 ? 0.0 : radii[i-1]) + radii[i])/2 for i in eachindex(radii)]
    for field in SteppedLayeredModel_CONVERSION_FIELDS
        val = getfield(m, field)
        isempty(val) && continue
        d[field] = evaluate.(m, field, eval_radii)
    end
    d[:r] = radii
    SteppedLayeredModel(; d...)
end


text:
using XLStats
using Test

data_dir="./data"
xml_data="$data_dir/hitsDetail.xls"
links = read_all(xml_file=xml_data,
                 topolink_log="$data_dir/salbiii_topolink.log",
                 topolink_input="$data_dir/topolink.inp",
                 xic_file_name="$data_dir/salbiii_xic.dat",
                 domain=2:134)

@testset "XLStats.jl" begin

  @test name(links[1]) == "D3-E53"

  @test indexes(links[13]) == (6, 17)

  @test resnames(links[45]) == ("GLU","GLU")

  @test length(links) == 113

  @test sum(nscans(links)) == 761

  @test count(consistency(links)) == 26

  @test sum(skipmissing(maxscore1(links))) ≈ 297.84999999999997

  @test sum(skipmissing(maxscore2(links))) ≈ 86.866

  @test sum(skipmissing(avgscore1(links))) ≈ 120.55580639396922

  @test sum(skipmissing(avgscore2(links))) ≈ 30.502086812287494

  @test sum(deuc(links)) ≈ 2208.976

  @test sum(dtop(links)) ≈ 3008.168

  @test sum(dmax(links)) ≈ 1680.1999999999998

  @test count(hasxic.(links)) == 20

  @test sum(skipmissing(maxxic(links))) ≈ 6.2237511644e8

  @test sum(avgxic(links)) ≈ 2.686064260940843e8

  @test sum(nspecies.(links)) == 159

  @test sum(count_nspecies.(links,mtol=20.0)) == 154

  @test count(consistency(links,tol=2.0)) == 29

  @test ismatch(resnames(links[2]),("ASP","GLU"))
  @test ismatch(resnames(links[2]),("GLU","ASP"))

  @test ismatch(indexes(links[2]),(3,111))
  @test ismatch(indexes(links[2]),(111,3))

  @test ismatch("D143-E108","D143-E108")
  @test ismatch("D108-E143","D143-E108")
  @test ismatch("D108-D143","D143-E108") == false
  @test ismatch("D108-D143","D100-E108") == false

  @test indomain(links[2],1:131)
  @test indomain(links[112],1:120) == false

  @test point_biserial(consistency(links),deuc(links)) ≈ -0.5950247343012209

end


text:
using CpuId
using Libdl

## this lets us load CpuId only once

if cpufeature(:AVX2)
    libsuffix = :avx2
    println("AVX2 support detected, vml_avx2 selected")
else
    libsuffix = :avx
    println("AVX2 support missing, vml_avx selected")
end

if Sys.iswindows()
    rtlib = :mkl_rt
    corelib = :mkl_core
    lib = Symbol(:mkl_vml_, libsuffix)
else
    rtlib = :libmkl_rt
    corelib = :libmkl_core
    lib = Symbol(:libmkl_vml_, libsuffix)
end

depsjl_path = joinpath(@__DIR__, "deps.jl")
open(depsjl_path, "w") do depsjl_file
    println(depsjl_file, strip("""
        ## This file was autogenerated by build.jl.
        ## Do not edit.
        import Libdl

        const lib = :$lib
        const rtlib = :$rtlib
        const corelib = :$corelib
        """))
end