text:
using PyCall
@pyimport matplotlib.pyplot as plt

@doc """
Sampler class
""" ->
type Sampler
	num_particles::Int64
	mcmc_steps::Int64
	particles::Vector{Particle}
	logl::Vector{Float64}
	tiebreakers::Vector{Float64}

	# Current iteration
	iteration::Int64

	# Current log likelihood threshold
	logx_threshold::Float64
	logl_threshold::Float64
	tb_threshold::Float64
end

@doc """
Constructor that only takes num_particles and mcmc_steps
as input
""" ->
function Sampler(num_particles::Int64, mcmc_steps::Int64)
	@assert (num_particles >= 1) & (mcmc_steps >= 1)
	return Sampler(num_particles, mcmc_steps,
								Array(Particle, (num_particles, )),
								zeros(num_particles), zeros(num_particles),
								0, 0.0, -Inf, 0.0)
end

@doc """
Generate all particles from the prior
""" ->
function initialise!(sampler::Sampler)
	for i in 1:sampler.num_particles
		sampler.particles[i] = Particle()
		from_prior!(sampler.particles[i])
		sampler.logl[i] = log_likelihood(sampler.particles[i])
		sampler.tiebreakers[i] = rand()
	end
	return nothing
end

@doc """
Find and save worst particle,
then generate replacement.
""" ->
function do_iteration!(sampler::Sampler, verbose::Bool)
	sampler.iteration += 1

    # Standard
	sampler.logx_threshold = -sampler.iteration/sampler.num_particles

    # Walter
    # sampler.logx_threshold = sampler.iteration*log(1.0 - 1.0/sampler.num_particles)

	# Find index of worst particle
	worst = find_worst_particle(sampler::Sampler)

	# Write its information to the output files
	if(sampler.iteration == 1)
		f = open("sample_info.txt", "w")
		write(f, "# num_particles, iteration, log(X), log(L)\n")
		f2 = open("sample.txt", "w")
		write(f2, "# The samples themselves. Use log(X) from sample_info.txt as un-normalised prior weights.\n")
	else
		f = open("sample_info.txt", "a")
		f2 = open("sample.txt", "a")
	end
	write(f, string(sampler.num_particles), " ", string(sampler.iteration), " ",
			string(sampler.logx_threshold, " ", string(sampler.logl[worst]), "\n"))
	close(f)
	write(f2, string(string(sampler.particles[worst]), "\n"))
	close(f2)

	# Set likelihood threshold
	sampler.logl_threshold = sampler.logl[worst]
	if(verbose)
		println("Iteration ", sampler.iteration, ", log(X) = ",
				sampler.logx_threshold,	", log(L) = ", sampler.logl_threshold)
	end

	# Clone a survivor
	if(sampler.num_particles != 1)
		which = rand(1:sampler.num_particles)
		while(which == worst)
			which = rand(1:sampler.num_particles)
		end
		sampler.particles[worst] = deepcopy(sampler.particles[which])
		sampler.logl[worst] = deepcopy(sampler.logl[which])
		sampler.tiebreakers[worst] = deepcopy(sampler.tiebreakers[which])
	end

	# Evolve
	accepted = 0::Int64
	for i in 1:sampler.mcmc_steps
		proposal = deepcopy(sampler.particles[worst])
		logH = perturb!(proposal)
		if(logH > 0.0)
			logH = 0.0
		end
		logl_proposal = log_likelihood(proposal)
		tb_proposal = sampler.tiebreakers[worst] + randh()
		tb_proposal = mod(tb_proposal, 1.0)

		if((rand() <= exp(logH)) && is_less_than(
						(sampler.logl_threshold, sampler.tb_threshold),
						(logl_proposal, tb_proposal)))
			sampler.particles[worst] = proposal
			sampler.logl[worst] = logl_proposal
			sampler.tiebreakers[worst] = tb_proposal
			accepted += 1
		end
	end
	if(verbose)
		println("Accepted ", accepted, "/", sampler.mcmc_steps, " MCMC steps.\n")
	end

	return (sampler.logx_threshold, sampler.logl_threshold)
end

@doc """
Compare based on likelihoods first. Use tiebreakers to break a tie
""" ->
function is_less_than(x::Tuple{Float64, Float64}, y::Tuple{Float64, Float64})
	if(x[1] < y[1])
		return true
	end
	if((x[1] == y[1]) && (x[2] < y[2]))
		return true
	end
	return false
end

@doc """
Find the index of the worst particle.
""" ->
function find_worst_particle(sampler::Sampler)
	# Find worst particle
	worst = 1
	for i in 2:sampler.num_particles
		if(is_less_than((sampler.logl[i], sampler.tiebreakers[i]),
						(sampler.logl[worst], sampler.tiebreakers[worst])))
			worst = i
		end
	end
	return worst
end

@doc """
Calculate the log evidence, information, and posterior weights from the output of a run
""" ->
function calculate_logZ(logX::Vector{Float64}, logL::Vector{Float64})
    # Add a zero to the beginning of logX
    logX2 = [0.0; logX]

	# Prior weights
	log_prior = Array(Float64, (length(logX), ))
    for i in 1:length(log_prior)
        log_prior[i] = logdiffexp(logX2[i], logX2[i+1])
    end

	# Unnormalised posterior weights
	log_post = log_prior + logL

	# log evidence and information
	logZ = logsumexp(log_post)
    post = exp(log_post - logZ)
	H = sum(post.*(log_post - logZ - log_prior))
	return (logZ, H, log_post)
end

@doc """
Do a Nested Sampling run.
""" ->
function do_nested_sampling(num_particles::Int64, mcmc_steps::Int64,
												depth::Float64; plot=true,
												verbose=true)

	# Number of NS iterations
	steps = Int64(max_depth*mcmc_steps)

	# Create the sampler
	sampler = Sampler(num_particles, mcmc_steps)
	initialise!(sampler)

	# Do 'steps' iterations of NS
	# Storage for results
	steps = Int64(max_depth)*num_particles
	plot_skip = num_particles

	# Store logX, logL
	keep = Array(Float64, (steps, 2))

	for i in 1:steps
		(keep[i, 1], keep[i, 2]) = do_iteration!(sampler, verbose)

		if(plot && (rem(i, plot_skip) == 0))
			(logZ, H, log_post) = calculate_logZ(keep[1:i, 1], keep[1:i, 2])

            uncertainty = NaN
            try
    			uncertainty = sqrt(H/sampler.num_particles)
            catch
            end

			plt.subplot(2, 1, 1)
            plt.cla()
			plt.plot(keep[1:i, 1], keep[1:i, 2], "ko-", markersize=1)
			plt.ylabel("\$\\ln(L)\$")
			plt.title(string("\$\\ln(Z) =\$ ", signif(logZ, 6),
						" +- ", signif(uncertainty, 3),
						", \$H = \$", signif(H, 6), " nats"))

			# Adaptive ylim (exclude bottom 5%)
			logl_sorted = sort(keep[1:i, 2])
			lower = logl_sorted[1 + Int64(floor(0.05*i))]
			plt.ylim([lower, logl_sorted[end] + 0.05*(logl_sorted[end] - lower)])

			plt.subplot(2, 1, 2)
            plt.cla()
			plt.plot(keep[1:i], exp(log_post - maximum(log_post)), "ko-", markersize=1)
			plt.xlabel("\$\\ln(X)\$")
			plt.ylabel("Relative posterior weights")

            plt.savefig("progress_plot.png", bbox_inches="tight")
		end
	end

    results = calculate_logZ(keep[:,1], keep[:,2])
    log_post = results[3] - results[1]
    post = exp(log_post)
    max_post = maximum(post)
    ESS = exp(-sum(post.*log_post))
    writedlm("weights.txt", post)

    # Do resampling do create posterior_sample.txt
    sample = readdlm("sample.txt")
    posterior_sample = Array(Float64, (Int64(ceil(ESS)), size(sample)[2]))
    for i in 1:size(posterior_sample)[1]
        which = 1
        while(true)
            which = rand(1:size(sample)[1])
            if(rand() <= post[which]/max_post)
                break
            end
        end
        posterior_sample[i, :] = sample[which, :]
    end
    writedlm("posterior_sample.txt", posterior_sample)

	if(verbose)
        (logZ, H, log_post) = results
        println("ln(Z) = ", signif(logZ, 12), " +- ", signif(sqrt(H/num_particles), 3))
        println("H = ", signif(H, 12), " nats")
        println("Effective posterior sample size = ", ESS)
		println("Done!")
	end
	if(plot)
		plt.show()
	end
    return(results)
end


text:
include("FileInput.jl")

export get_colnames
export get_f1
export get_f2
export set_colnames!
export set_f1!
export set_f2!

ext: jl

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# # Sinkhorn divergences
#
# In this tutorial we provide a minimal example for using the Sinkhorn divergence as a loss function [FSV+19] on empirical distributions.
# [FSV+19]: Feydy, Jean, et al. "Interpolating between optimal transport and MMD using Sinkhorn divergences." The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019.
#
# While entropy-regularised optimal transport $\operatorname{OT}_{\varepsilon}(\cdot, \cdot)$ is commonly used as a loss function, it suffers from a problem of *bias*: namely that $\nu \mapsto \operatorname{OT}_{\varepsilon}(\mu, \nu)$ is *not* minimised at $\nu = \mu$.
#
# A fix to this problem is proposed by Genevay et al [GPC18] and subsequently Feydy et al. [FSV+19], which introduce the *Sinkhorn divergence* between two measures $\mu$ and $\nu$, defined as
# ```math
# \operatorname{S}_{\varepsilon}(\mu, \nu) = \operatorname{OT}_{\varepsilon}(\mu, \nu) - \frac{1}{2} \operatorname{OT}_{\varepsilon}(\mu, \mu) - \frac{1}{2} \operatorname{OT}_{\varepsilon}(\nu, \nu).
# ```
# In the above, we have followed the convention taken by Feydy et al. and included the entropic regularisation in the definition of $\operatorname{OT}_\varepsilon$.
# [GPC18]: Aude Genevay, Gabriel Peyré, Marco Cuturi, Learning Generative Models with Sinkhorn Divergences, Proceedings of the Twenty-First International Conference on Artficial Intelligence and Statistics, (AISTATS) 21, 2018
# [FSV+19]: Feydy, Jean, et al. "Interpolating between optimal transport and MMD using Sinkhorn divergences." The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019.
#
# Like the Sinkhorn loss, the Sinkhorn divergence is smooth and convex in both of its arguments. However, the Sinkhorn divergence is unbiased -- i.e. $S_{\varepsilon}(\mu, \nu) = 0$ iff $\mu = \nu$.
#
# Unlike previous examples, here we demonstrate a learning problem similar to Figure 1 of Feydy et al. over *empirical measures*, i.e. measures that have the form $\mu = \frac{1}{N} \sum_{i = 1}^{N} \delta_{x_i}$ where $\delta_x$ is the Dirac delta function at $x$.
#
# We first load packages.
#
using OptimalTransport
using ReverseDiff
using Distributions
using LinearAlgebra
using Distances
using Plots
using Logging
using Optim

# As a ground truth distribution, we set $\rho$ to be a Gaussian mixture model with `k = 3` components, equally spaced around a circle, and sample an empirical distribution of size $N$, $\mu \sim \rho$.

k = 3
d = 2
θ = π * range(0, 2(1 - 1 / k); length=k)
μ = 2 * hcat(sin.(θ), cos.(θ))
ρ = MixtureModel(MvNormal[MvNormal(x, 0.25 * I) for x in eachrow(μ)])
N = 100
μ_spt = rand(ρ, N)'
scatter(μ_spt[:, 1], μ_spt[:, 2]; markeralpha=0.25, title=raw"$\mu$")

# Now, suppose we want to approximate $\mu$ with another empirical distribution $\nu$, i.e. we want to minimise $\nu \mapsto \operatorname{S}_{\varepsilon}(\mu, \nu)$ over possible empirical distributions $\nu$. In this case we have $M$ particles in $\nu$, which we initialise following a Gaussian distribution.
M = 100
ν_spt = rand(M, d);
# Assign uniform weights to the Diracs in each empirical distribution.
μ = fill(1 / N, N)
ν = fill(1 / M, M);

# Since $\mu$ is fixed, we pre-compute the cost matrix $C_{\mu}$.
C_μ = pairwise(SqEuclidean(), μ_spt');

# Define the loss function to minimise, where `x` specifies the locations of the Diracs in $\nu$.
#
# We will be using `ReverseDiff` with a precompiled tape. For this reason, we need the Sinkhorn algorithm to perform a fixed number of (e.g. 50) iterations.
# Currently, this can be achieved by setting `maxiter = 50` and `atol = rtol = 0` in calls to `sinkhorn` and `sinkhorn_divergence`.
function loss(x, ε)
    C_μν = pairwise(SqEuclidean(), μ_spt', x')
    C_ν = pairwise(SqEuclidean(), x')
    return sinkhorn_divergence(
        μ, ν, C_μν, C_μ, C_ν, ε; maxiter=50, atol=rtol = 0, regularization=true
    )
end
# Set entropy regularisation parameter
ε = 1.0;

# Use ReverseDiff with a precompiled tape and Optim.jl to minimise $\nu \mapsto \operatorname{S}_{\varepsilon}(\mu, \nu)$. Note that this is problem is *not* convex, so we find a local minimium.
const loss_tape = ReverseDiff.GradientTape(x -> loss(x, ε), ν_spt)
const compiled_loss_tape = ReverseDiff.compile(loss_tape)
opt = with_logger(SimpleLogger(stderr, Logging.Error)) do
    optimize(
        x -> loss(x, ε),
        (∇, x) -> ReverseDiff.gradient!(∇, compiled_loss_tape, x),
        ν_spt,
        GradientDescent(),
        Optim.Options(; iterations=10, g_tol=1e-6, show_trace=true),
    )
end
ν_opt = Optim.minimizer(opt)
plt1 = scatter(μ_spt[:, 1], μ_spt[:, 2]; markeralpha=0.25, title="Sinkhorn divergence")
scatter!(plt1, ν_opt[:, 1], ν_opt[:, 2]);

# For comparison, let us do the same computation again, but this time we want to minimise $\nu \mapsto \operatorname{OT}_{\varepsilon}(\mu, \nu)$.
function loss_biased(x, ε)
    C_μν = pairwise(SqEuclidean(), μ_spt', x')
    return sinkhorn2(μ, ν, C_μν, ε; maxiter=50, atol=rtol = 0, regularization=true)
end
const loss_biased_tape = ReverseDiff.GradientTape(x -> loss_biased(x, ε), ν_spt)
const compiled_loss_biased_tape = ReverseDiff.compile(loss_biased_tape)
opt_biased = with_logger(SimpleLogger(stderr, Logging.Error)) do
    optimize(
        x -> loss_biased(x, ε),
        (∇, x) -> ReverseDiff.gradient!(∇, compiled_loss_biased_tape, x),
        ν_spt,
        GradientDescent(),
        Optim.Options(; iterations=10, g_tol=1e-6, show_trace=true),
    )
end
ν_opt_biased = Optim.minimizer(opt_biased)
plt2 = scatter(μ_spt[:, 1], μ_spt[:, 2]; markeralpha=0.25, title="Sinkhorn loss")
scatter!(plt2, ν_opt_biased[:, 1], ν_opt_biased[:, 2]);

# Observe that the Sinkhorn divergence results in $\nu$ that matches $\mu$ quite well, while entropy-regularised transport is biased to producing $\nu$ that seems to concentrate around the mean of each Gaussian component.
plot(plt1, plt2)


text:
"""
    FaceContinuous(b::Int, n::Int)
    FaceContinuous(::Type{T}, b::Int, n::Int)

For `n` dimensions, use `b` bits of precision in this Hilbert curve.
If you specify a type `T`, this will be used as the type of the Hilbert encoding.
If not, the smallest unsigned integer that can hold `n*b` bits will be used
for the Hilbert index data type.

This is the Butz algorithm, as presented by Lawder. Haverkort2017
says it is face continuous. The code is in lawder.c.
The original paper had an error, and Lawder put a correction on his website.
http://www.dcs.bbk.ac.uk/~jkl/publications.html
"""
struct FaceContinuous{T} <: HilbertAlgorithm{T}
    b::Int
    n::Int
end

function FaceContinuous(b::Int, n::Int)
    T = large_enough_unsigned(b * n)
    FaceContinuous{T}(b, n)
end

FaceContinuous(::Type{T}, b::Int, n::Int) where {T} = FaceContinuous{T}(b, n)

"""
Used during testing to pick a type for the xyz coordinate.
"""
axis_type(gg::FaceContinuous) = large_enough_unsigned(gg.b)

# Mark one-based variables with a trailing o.
# Mark zero-based variables with a trailing z.
# The C code takes place in vectors, without the packed Hilbert index.
# ORDER -> b  The number of bits.
# DIM -> n  The number of dimensions.
# U_int -> A, the type for the Axis coordinate.
# C translation rules:
# * / -> ÷
# * ^ -> ⊻
# * % -> %  This is the same choice. mod() is wrong for negative numbers.
# Julia has different operator precedence among addition, subtraction, bit shifts, and and xor.
# Lawder knew C operator precedence well, so there are few parentheses. Add them liberally.
# https://en.cppreference.com/w/c/language/operator_precedence. Relevant ones:
# (++, --), (* / %), (+ -), (<< >>), (< >), (== !=), (&), (^), (= += <<= &= ^= |=)

# g_mask[x] can be replaced by (1 << DIM - 1 - x)
g_mask(::Type{A}, n, iz) where {A} = one(A) << (n - iz - 1)

function H_encode!(gg::FaceContinuous, pt::Vector{K}, h::Vector{K}) where {K}
    wordbits = 8 * sizeof(K)
    W = zero(K)
    P = zero(K)
    h .= zero(K)
    # Start from the high bit in each element and work to the low bit.
    mask = one(K) << (gg.b - 1)  # lawder puts wordbits here.
    # The H-index needs b * n bits of storage. i points to the base of the current level.
    i = gg.b * gg.n - gg.n
    # A will hold bits from all pt elements, at the current level.
    A = zero(K)
    for j = 0:(gg.n - 1)
        if (pt[j + 1] & mask) != zero(K)
            A |= g_mask(K, gg.n, j)
        end
    end
    S = tS = A
    P |= (S & g_mask(K, gg.n, 0))
    for j = 1:(gg.n - 1)
        gm = g_mask(K, gg.n, j)
        if ((S & gm) ⊻ ((P >> 1) & gm)) != 0
            P |= gm
        end
    end

    # P is the answer, but it has to be packed into a vector.
    element = i ÷ wordbits
    if (i % wordbits) > (wordbits - gg.n)
        h[element + 1] |= (P << (i % wordbits))
        h[element + 2] |= P >> (wordbits - (i % wordbits))
	else
        h[element + 1] |= (P << (i - element * wordbits))
    end

    J = gg.n
    j = 1
    while j < gg.n
        if ((P >> j) & one(K)) == (P & one(K))
            j += 1
			continue
		else
            break
        end
    end
    if j != gg.n
        J -= j
    end
    xJ = J - 1

    if P < K(3)
        T = zero(K)
    else
        if (P % K(2)) != zero(K)
            T = (P - one(K)) ⊻ ((P - one(K)) >> 1)
        else
            T = (P - K(2)) ⊻ ((P - K(2)) >> 1)
        end
    end
    tT = T
    i -= gg.n
    mask >>= 1
    while i >= 0
        # println("i=$i mask=$mask T=$T P=$P")
		# println("h[0]=$(h[1]) h[1]=$(h[2]) h[2]=$(h[3])")
        A = zero(K)
        for j = 0:(gg.n - 1)
            if pt[j + 1] & mask != zero(K)
                A |= g_mask(K, gg.n, j)
            end
        end
        W ⊻= tT
        tS = A ⊻ W
        if xJ % gg.n != 0
            temp1 = tS << (xJ % gg.n)
            temp2 = tS >> (gg.n - (xJ % gg.n))
            S = temp1 | temp2
            S &= (one(K) << gg.n) - one(K)
        else
            S = tS
        end

        P = S & g_mask(K, gg.n, 0)
        for j = 1:(gg.n - 1)
            gn = g_mask(K, gg.n, j)
            if ((S & gn) ⊻ ((P >> 1) & gn)) != 0
                P |= gn
            end
        end
        element = i ÷ wordbits
        if (i % wordbits) > (wordbits - gg.n)
            h[element + 1] |= (P << (i % wordbits))
            h[element + 2] |= (P >> (wordbits - (i % wordbits)))
        else
            h[element + 1] |= P << (i - element * wordbits)
        end
        if i > 0
            if P < K(3)
                T = 0
            else
                if P % K(2) != zero(K)
                    T = (P - one(K)) ⊻ ((P - one(K)) >> 1)
                else
                    T = (P - K(2)) ⊻ ((P - K(2)) >> 1)
                end
            end
            if xJ % gg.n != 0
                temp1 = T >> (xJ % gg.n)
                temp2 = T << (gg.n - (xJ % gg.n))
                tT = temp1 | temp2
                tT &= (one(K) << gg.n) - one(K)
            else
                tT = T
            end

            J = gg.n
            j = 1
            while j < gg.n
                if ((P >> j) & one(K)) == (P & one(K))
                    j += 1
                    continue
                else
                    break
                end
            end
            if j != gg.n
                J -= j
            end
            xJ += J - 1
        end

        i -= gg.n
        mask >>= 1
    end
    # println("h[0]=$(h[1]) h[1]=$(h[2]) h[2]=$(h[3])")
end

mutable struct FCLevel{K}
    mask::K
    i::Int
    n::Int
end

function FCLevel(fc, ::Type{K}) where {K}
    # XXX lawder uses wordbits instead of fc.b
    FCLevel{K}(one(K) << (fc.b - 1), fc.b * fc.n - fc.n, fc.n)
end

function downlevel!(l::FCLevel)
    l.mask >>= 1
    l.i -= l.n
end

function index_at_level(H::Vector{K}, l::FCLevel{K}) where {K}
    wordbits = 8 * sizeof(K)
    element = l.i ÷ wordbits
    P = H[element + 1]
    if (l.i % wordbits) > (wordbits - l.n)
        temp1 = H[element + 2]  # one-based
        P >>= l.i % wordbits
        temp1 <<= wordbits - (l.i % wordbits)
        P |= temp1
    else
        P >>= (l.i % wordbits)
    end

    # 	/* the & masks out spurious highbit values */
    if l.n < wordbits
        P &= (one(K) << l.n) - one(K)
    end
    P
end

function distribute_to_coords!(bits::K, axes::Vector{K}, l::FCLevel{K}) where K
    j = l.n - 1
    while bits > 0
        if bits & one(K) != 0
            axes[j + 1] |= l.mask
        end
        bits >>= 1
        j -= 1
    end
end

function fc_parity_match(P, n)
    J = n
    j = 1
    parity = P & one(P)
    while j < n
        if ((P >> j) & one(P)) != parity
            break
        end
        j += 1
    end
    j
end

function fc_rotation(P, n)
    j = fc_parity_match(P, n)
    if j == n
        xJ = n - 1
    else
        xJ = n - j - 1
    end
    xJ
end

function fc_flip(P::K) where {K}
    if P < K(3)
        T = zero(K)
    else
        if P % 2 != 0
            T = (P - one(K)) ⊻ ((P - one(K)) >> 1)
        else
            T = (P - K(2)) ⊻ ((P - K(2)) >> 1)
        end
    end
    T
end

function H_decode!(gg::FaceContinuous, H::Vector{K}, pt::Vector{K}) where {K}
    pt .= zero(K)
    l = FCLevel(gg, K)
    P = index_at_level(H, l)
    xJ = fc_rotation(P, gg.n)
    A = S = tS = brgc(P)
    tT = T = fc_flip(P)
    # XXX Lawder's code puts P here instead of A.
    distribute_to_coords!(A, pt, l)

    nmask = (one(K) << gg.n) - one(K)  # mask of n bits.

    W = zero(K)
    downlevel!(l)
    while l.i >= 0
        P = index_at_level(H, l)
        S = brgc(P)
        @assert S & ~nmask == 0
        tS = rotateright(S, xJ % gg.n, gg.n)

        W ⊻= tT
        A = W ⊻ tS
        distribute_to_coords!(A, pt, l)

        if l.i > 0
            T = fc_flip(P)
            @assert T & ~nmask == 0
            tT = rotateright(T, xJ % gg.n, gg.n)

            xJ += fc_rotation(P, gg.n)
        end
        downlevel!(l)
    end
end

function encode_hilbert_zero(fc::FaceContinuous{T}, X::Vector{A})::T where {A,T}
    hvec = zeros(A, fc.n)
    H_encode!(fc::FaceContinuous, X, hvec)
    # Encoding packs H into a vector of A, using all bits in the A type.
    h = zero(T)
    for i in fc.n:-1:1
        h <<= 8 * sizeof(A)
        h |= hvec[i]
    end
    h
end

function decode_hilbert_zero!(fc::FaceContinuous{T}, X::Vector{A}, h::T) where {A,T}
    # H is in a larger type T but algorithm expects it packed into a vector of A.
    hvec = zeros(A, fc.n)
    for i in 1:fc.n
        hvec[i] |= h & ~zero(A)
        h >>= 8 * sizeof(A)
    end
    H_decode!(fc::FaceContinuous, hvec, X)
end


text:
using TcpInstruments
using Test

try
    include("./fake_device.jl")
catch
    @error "External emulator NOT found"
end

function test_initialize(type)
    @info "Creating fake $type instrument at localhost:8080"
    handle = initialize(type, "127.0.0.1:8080")
    @info handle
    @test handle.connected
    @info "Successfully connected"
    return handle
end

# Start Tests
@testset "Fake Instruments" begin

@testset "Scope Interface" begin
    handle = test_initialize(AgilentDSOX4034A)

    id = info(handle)
    @test id == "1"

    lpf_on(handle)
    lpf_off(handle)

    set_impedance_1Mohm(handle)
    set_impedance_50ohm(handle)

    terminate(handle)
    @info "Successfully disconnected"
end

@testset "Keysight34465A" begin
    handle = test_initialize(KeysightDMM34465A)

    id = info(handle)
    @test id == "1"

    set_temp_unit_kelvin(handle)
    set_temp_unit_farenheit(handle)
    set_temp_unit_celsius(handle)

    terminate(handle)
    @info "Successfully disconnected"
end

@testset "SRSPS310" begin
    handle = test_initialize(SRSPS310)
    enable_output(handle)
    disable_output(handle)
    set_current_limit(handle, 1A)
    set_current_limit(handle, 2A)
    set_current_limit(handle, 5A)
    terminate(handle)
    @info "Successfully disconnected"
end

end # End Tests

try
close_emulator()
catch
end


text:
function adjoint(A)
    @static if VERSION < v"0.7-"
        A'
    else
        copy(A')
    end
end

function approximate_spectral_radius(A, tol = 0.01,
                                        maxiter = 15, restart = 5)

    symmetric = false

    # Initial guess
    v0 = rand(eltype(A), size(A,2))
    maxiter = min(size(A, 1), maxiter)
    ev = zeros(eltype(A), maxiter)
    max_index = 0
    X = zeros(eltype(A), size(A,1), maxiter)

    for i in 1:restart+1
        evect, ev, H, V, flag =
                    approximate_eigenvalues(A, tol, maxiter,
                                            symmetric, v0)
        nvecs = size(ev, 1)
        # X = hcat(V[1:end-1]...)
        copy_V!(X, V)
        # m, max_index = findmax(abs.(ev))
        m, max_index = findmaxabs(ev)
        error = H[nvecs, nvecs-1] * evect[end, max_index]
        @static if VERSION < v"0.7-"
            @views A_mul_B!(v0, X, evect[:, max_index])
        else
            @views mul!(v0, X, evect[:, max_index])
        end
        if (abs(error) / abs(ev[max_index]) < tol) || flag
            # v0 = X * evect[:, max_index]
            break
        end
    end
    ρ = abs(ev[max_index])
end

function findmaxabs(arr)
    m = abs(arr[1])
    m_i = 1
    for i = 2:size(arr, 1)
        x = abs(arr[i])
        if x > m
            m = x
            m_i  = i
        end
    end
    m, m_i
end

function copy_V!(X, V)
    n = size(V,1)
    for i = 1:n-1
        X[:,i] = V[i]
    end
end

function approximate_eigenvalues(A, tol, maxiter, symmetric, v0)

    # maxiter = min(size(A, 1), maxiter)
    v0 ./= norm(v0)
    H = zeros(eltype(A), maxiter + 1, maxiter)
    V = [v0]
    # V = Vector{Vector{eltype(A)}}(maxiter + 1)
    # V = zeros(size(A,1), maxiter)
    # V[1] = v0
    breakdown = find_breakdown(eltype(A))
    flag = false

    for j = 1:maxiter
        w = A * V[end]
        # V[j+1] = A * V[j]
        # w = V[j+1]
        # mul!(w, A, V[j])
        for (i,v) in enumerate(V)
        # for i = 1:j
            v = V[i]
            H[i,j] = dot(conj.(v), w)
            BLAS.axpy!(-H[i,j], v, w)
        end
        H[j+1,j] = norm(w)
        if H[j+1, j] < breakdown
            flag = true
            if H[j+1,j] != 0
                rmul!(w, 1/H[j+1,j])
                push!(V, w)
                break
            end
        end

        #w = w / H[j+1, j]
        rmul!(w, 1/H[j+1,j])
        push!(V, w)
    end
    @static if VERSION < v"0.7-"
        Eigs, Vects = eig(H[1:maxiter, 1:maxiter], Matrix{eltype(A)}(I, maxiter, maxiter))
    else
        Eigs, Vects = (eigen(H[1:maxiter, 1:maxiter], Matrix{eltype(A)}(I, maxiter, maxiter))...,)
    end

    Vects, Eigs, H, V, flag
end

find_breakdown(::Type{Float64}) = eps(Float64) * 10^6
find_breakdown(::Type{Float32}) = eps(Float64) * 10^3

using Base.Threads
#=function mul!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)
    A.n == size(B, 1) || throw(DimensionMismatch())
    A.m == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = A.nzval
    rv = A.rowval
    if β != 1
        β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
    end
    for k = 1:size(C, 2)
        for col = 1:A.n
            αxj = α*B[col,k]
            for j = A.colptr[col]:(A.colptr[col + 1] - 1)
                C[rv[j], k] += nzv[j]*αxj
            end
        end
    end
    C
end
spmatvec(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx} =
    (T = promote_type(TA, Tx); mul!(one(T), A, x, zero(T), similar(x, T, A.m)))=#

# export spmatvec

#=@static if haskey(ENV, "JULIA_NUM_THREADS")
    info("Initializing AMG to use $(nthreads()) threads")
    import Base: *
    function *(A::SparseMatrixCSC{T,V}, b::Vector{T}) where {T,V}
        m, n = size(A)
        ret = zeros(T, m)
        @threads for i = 1:n
            for j in nzrange(A, i)
                row = A.rowval[j]
                val = A.nzval[j]
                ret[row] += val * b[i]
            end
        end
        ret
    end
end=#


text:
"""
Module containing functionality to read and process tetrahedral meshes in gmsh or nastran format.
"""
module Meshutils
import LinearAlgebra, ProgressMeter
include("Meshutils_exports.jl")

"""
Simple struct element to store additional information on symmetry for rotational symmetric meshes.

#Fields
- `DOS::Int64`: degree of symmetry
- `naxis::Int64`: number of grid points lying on the center axis
- `nxbloch::Int64`: number of grid points lying on the Bloch plane but not on the center axis
- `nbody::Int64`: number of interior points in half cell
- `shiftbody::Int64`: difference from body point to reflected body point
- `nxsymmetry::Int64`: number of grid points on symmetry plane (without center axis)
- `nxsector::Int64`: number of grid points belonging to a unit cell (without cenetraxis, Bloch and Bloch image plane)
- `naxis_ln::Int64`: number of line segments lying on the center axis
- `nxbloch_ln::Int64`: number of line segments belonging to a unit cell (without cenetraxis, Bloch and Bloch image plane)
- `nxsector_ln::Int64`: number of line segments belonging to a unit cell (without cenetraxis, Bloch and Bloch image plane)
- `nxsector_tri::Int64`: number of surface triangles belonging to a unit cell (without Bloch and Bloch image plane)
- `nxsector_tet::Int64`: number of tetrahedra of a unit cell
- `n`: unit axis vector of the center axis
- `pnt: foot point of the center axis
- `unit::Bool`: true if mesh represents only a unit cell
"""
struct SymInfo
    DOS::Int64 #degree of symmetry
    naxis::Int64 #number of gridpoints lying on the centeraxis
    nxbloch::Int64
    nbody::Int64
    shiftbody::Int64
    nxsymmetry::Int64
    nxsector::Int64
    naxis_ln::Int64
    nxbloch_ln::Int64
    nxsector_ln::Int64
    nxsector_tri::Int64
    nxsector_tet::Int64
    n
    pnt
    unit::Bool
end

"""
Definition of the mesh type

# Fields
- `name::String`: the name of the mesh.
- `points::Array`: 3×N array containing the coordinates of the N points defining the mesh.
- `lines::List`: List of simplices defining the edges of the mesh
- `triangles::List`: List of simplices defining the surface triangles of the mesh
- `int_triangles::List`: List of simplices defining the interior triangles of the mesh
- `tetrahedra::List`: List of simplices defining the tetrahedra of the mesh
- `domains::Dict`: Dictionary defining the domains of the mesh. See comments below.
- `file::String`: path to the file containing the mesh.
- `tri2tet::Array`: Array of length `length(tetrahedra)` containing the indices of the connected tetrahedra.
- `dos`: special field meant to contain symmetry information of highly symmetric meshes.

# Notes
The meshes are supposed to be tetrahedral. All simplices (lines, triangles, and tetrahedra) are stored as lists of simplices.
Simplices are lists of integers containing the indices of the points (i.e. the column number in the `points` array) forming the simplex.
This means a line is a two-entry list, a triangle a three-entry list, and a tetrahedron a four-entry list.
For convenience certain entities of the mesh can be further defined in the `domains` dictionary. Each key defines a domain and maps to another dictionary.
This second-level dictionary contains at least two keys: `"dimension"` mapping to the dimension of the specified domain (1,2, or 3) and
`"simplices"` containing a list of integers mapping into the respective simplex lists. More keys may be added to the dictionary to
define additional and/or custom information on the domain. For instance the `compute_size!` function adds an entry with the domain size.
"""
struct Mesh
    name
    points
    lines
    triangles
    int_triangles
    tetrahedra
    domains
    file
    tri2tet
    dos
end
## constructor
"""
    mesh=Mesh(file_name::String; scale=1, inttris=false)

read a tetrahedral mesh from gmsh or nastran file into `mesh`.
The optional scaling factor `scale` may be used to scale the units of the mesh.
The optional parameter `inttris` toggles whether triangles are sortet into two
seperate lists for surface and internal triangles.
"""
function Mesh(file_name::String;scale=1,inttris=false)
    #mesh constructor
    ext=split(file_name,".")
    if ext[end]=="msh"
        #gmsh mesh file
        #check file format
        local line
        open(file_name,"r") do fid
            line=readline(fid)
            line=readline(fid)
        end
        if line[1]=='4'
            points, lines, triangles, tetrahedra, domains =read_msh4(file_name)
        else
            println("Please, convert msh-file to gmsh's .msh-format version 4!")
            #points, lines, triangles, tetrahedra, domains =read_msh2(file_name)
        end
    elseif ext[end]=="nas" || ext[end]=="bdf"
        points, lines, triangles, tetrahedra, domains =read_nastran(file_name)
        # if isfile(ext[1]*"_surfaces.txt")
        #     relabel_nas_domains!(domains,ext[1]*"_surfaces.txt")
        # end
    else
        println("mesh type not supported")
        return nothing
    end

    #sometimes elements are multiply defined, make them unique
    #TODO: This can be done more memory efficient by reading the msh-file twice
    ulines=Array{UInt32,1}[]
    utriangles=Array{UInt32,1}[]
    utetrahedra=Array{UInt32,1}[]
    uinttriangles=Array{UInt32,1}[]

    for ln in lines
        insert_smplx!(ulines,ln)
    end

    for tet in tetrahedra
        insert_smplx!(utetrahedra,tet)
    end

    if inttris
        utriangles,uinttriangles=assemble_triangles(utetrahedra)
    else
        for tri in triangles
            insert_smplx!(utriangles,tri)
        end
    end

    for dom in keys(domains)
        for (idx,smplx_idx) in enumerate(domains[dom]["simplices"])
            if domains[dom]["dimension"]==1
                new_idx=find_smplx(ulines,lines[smplx_idx])
            elseif domains[dom]["dimension"]==2
                new_idx=find_smplx(utriangles,triangles[smplx_idx])
            elseif domains[dom]["dimension"]==3
                new_idx=find_smplx(utetrahedra,tetrahedra[smplx_idx])
            end
            domains[dom]["simplices"][idx]=new_idx
        end
        unique!(domains[dom]["simplices"])
    end

    if inttris
        tri2tet=[zeros(UInt32,length(utriangles)),zeros(UInt32,length(uinttriangles),2)]
    else
        tri2tet=zeros(UInt32,length(utriangles))
        if length(tri2tet)!=0
            tri2tet[1]=0xffffffff #magic sentinel value to check whether list is linked
        end
    end
    Mesh(file_name, points*scale, ulines, utriangles, uinttriangles, utetrahedra, domains, file_name,tri2tet,1)
end
#deprecated constructor
function Mesh(file, points, lines, triangles, tetrahedra, domains, file_name, tri2tet, dos)
    return Mesh(file, points, lines, triangles, [], tetrahedra, domains, file_name, tri2tet, dos)
end

#show mesh
import Base.string
function string(mesh::Mesh)
    if isempty(mesh.int_triangles)
        txt="""mesh: $(mesh.file)
        #################
        points:     $(size(mesh.points)[2])
        lines:      $(length(mesh.lines))
        triangles:  $(length(mesh.triangles))
        tetrahedra: $(length(mesh.tetrahedra))
        #################
        domains: """
    else
        txt="""mesh: $(mesh.file)
        #################
        points:            $(size(mesh.points)[2])
        lines:             $(length(mesh.lines))
        triangles (surf):  $(length(mesh.triangles))
        triangles (int):   $(length(mesh.int_triangles))
        tetrahedra:        $(length(mesh.tetrahedra))
        #################
        domains: """
    end
    for key in sort([key for key in keys(mesh.domains)])
        txt*=string(key)*", "
    end
    return txt[1:end-2]
end
import Base.show
function show(io::IO,mesh::Mesh)
  print(io,string(mesh))
end

##
# import Base.getindex(mesh::Mesh, i::String)
#     #parse for sectors
#     splt=split(i,"_#")
#     sctr,hlf=-1,-1
#     dom=splt[1] #domain name
#     #TODO: error message of split has more than 3 entries
#     if length(splt)>1
#         splt=split(splt[2],".")
#         sctr=splt[1]
#         if length(splt)>1 #TODO: length check, see above
#             hlf=splt[2] #TODO: check that this value should be 0 or 1 (assert macro?)
#         end
#     end
#
#     #just call domain
#     domain = Dict()
#     domain["dimension"]=mesh.domains[dom]["dimension"]
#     smplx_list=[]
#     #TODO: compute volume if existent
#     if hlf>-1 #build half-cell
#         for smplx in mesh.domains[dom]["simplices"]
#             if hlf==1
#                 smplx=get_reflected_index.(smplx)
#             end
#             insert_smplx!(smplx_list, get_rotated_index.(smplx, sctr))
#         end
#
#     elseif sctr>-1 #build sector
#         for smplx in mesh.domains[dom]["simplices"]
#             rsmplx=get_reflected_index.(smplx)
#             insert_smplx!(smplx_list, get_rotated_index.(smplx, sctr))
#             insert_smplx!(smplx_list, get_rotated_index.(rsmplx, sctr))
#         end
#
#     else #build full annulus or just call domain
#         for smplx in mesh.domains[dom]["simplices"]
#             rsmplx=get_reflected_index.(smplx)
#             for sector=0:DOS-1
#                 insert_smplx!(smplx_list, get_rotated_index.(smplx, sector))
#                 insert_smplx!(smplx_list, get_rotated_index.(rsmplx, sector))
#             end
#         end
#     end
#
#
#     domain["simplices"] = smplx_idx
#     return domain
# end
#
#
# import Base.getindex(mesh::Mesh, i::Int)
#
# end

include("./Mesh/sorter.jl")
include("./Mesh/annular_meshes.jl")
include("./Mesh/vtk_write.jl")
##

"""
    points, lines, triangles, tetrahedra, domains=read_msh4(file_name)

Read gmsh's .msh-format version 4.

http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
"""
function read_msh4(file_name)
    lines=Array{UInt32,1}[]
    triangles=Array{UInt32,1}[]
    tetrahedra=Array{UInt32,1}[]
    tag2dom=Dict()
    ent2dom=Array{Dict}(undef,4)
    local points
    domains=Dict()
    open(file_name) do fid
        while !eof(fid)
            line=readline(fid)
            fldname=line[2:end]
            if fldname=="MeshFormat"
                line=readline(fid)
                #if split(line)[1] != "4.1"
                #    break
                #end

            elseif fldname=="PhysicalNames"
                line=readline(fid)
                numPhysicalNames=parse(UInt32,line)
                #Array{String}(undef,numPhysicalNames)
                for idx =1:numPhysicalNames
                    line=readline(fid)
                    dim, tag, dom=split(line)
                    dim=parse(UInt32,dim)
                    dom=String(dom[2:end-1])
                    tag2dom[tag]=dom
                    domains[dom]=Dict()
                    domains[dom]["dimension"]=dim
                    domains[dom]["simplices"]=[]
                end

            elseif fldname=="Entities"
                line=readline(fid)
                numPoints, numCurves, numSurfaces, numVolumes=parse.(UInt32,split(line))
                for idx=1:4 #inititlaize entity dictionary
                    ent2dom[idx]=Dict()
                end
                for idx=1:numPoints
                    line=readline(fid)
                    splitted=split(line)
                    entTag=splitted[1]
                    numPhysicalTags=parse(UInt32,splitted[5])
                    physicalTags=splitted[6:5+numPhysicalTags]
                    ent2dom[1][entTag]=[tag2dom[tag] for tag in physicalTags]
                end
                for idx=1:numCurves
                    line=readline(fid)
                    splitted=split(line)
                    entTag=splitted[1]
                    numPhysicalTags=parse(UInt32,splitted[8])
                    physicalTags=splitted[9:8+numPhysicalTags]
                    ent2dom[2][entTag]=[tag2dom[tag] for tag in physicalTags]
                end
                for idx=1:numSurfaces
                    line=readline(fid)
                    splitted=split(line)
                    entTag=splitted[1]
                    numPhysicalTags=parse(UInt32,splitted[8])
                    physicalTags=splitted[9:8+numPhysicalTags]
                    ent2dom[3][entTag]=[tag2dom[tag] for tag in physicalTags]
                end
                for idx=1:numVolumes
                    line=readline(fid)
                    splitted=split(line)
                    entTag=splitted[1]
                    numPhysicalTags=parse(UInt32,splitted[8])
                    physicalTags=splitted[9:8+numPhysicalTags]
                    ent2dom[4][entTag]=[tag2dom[tag] for tag in physicalTags]
                end

            elseif fldname=="Nodes"
                line=readline(fid)
                numEntityBlocks, numNodes,minNodeTag, maxNodeTag =parse.(UInt32,split(line))
                points=Array{Float64}(undef,3,numNodes)
                for idx =1:numEntityBlocks
                    line=readline(fid)
                    entityDim, entityTag, parametric, numNodesInBlock=parse.(UInt32,split(line))
                    #TODO: support parametric
                    nodeTags=Array{UInt32}(undef,numNodesInBlock)
                    for jdx =1:numNodesInBlock
                        line=readline(fid)
                        nodeTags[jdx]=parse(UInt32,line)
                    end
                    for jdx=1:numNodesInBlock
                        line=readline(fid)
                        xyz=parse.(Float64,split(line))
                        points[:,nodeTags[jdx]]=xyz[1:3]
                        #TODO: xyz[4:6] if parametric
                    end
                end
            elseif fldname=="Elements"
                line=readline(fid)
                numEntityBlocks, numElements, minElementTag, maxElementTag = parse.(UInt32,split(line))
                for idx=1:numEntityBlocks
                    line=readline(fid)
                    splitted=split(line)
                    entityDim=parse(UInt32,splitted[1])
                    entityTag=splitted[2]
                    elementType, numElementsInBlock = parse.(UInt32,splitted[3:4])
                    for jdx =1:numElementsInBlock
                        line=readline(fid)
                        nodeTags=parse.(UInt32,split(line))[2:end]
                        if elementType==1
                            append!(lines,[nodeTags]) #TODO: use simplexSorter!!!
                            for dom in ent2dom[entityDim+1][entityTag]
                                append!(domains[dom]["simplices"],length(lines))
                            end
                        elseif elementType==2
                            append!(triangles,[nodeTags]) #TODO: Use simplexSorter!!!
                            for dom in ent2dom[entityDim+1][entityTag]
                                append!(domains[dom]["simplices"],length(triangles))
                            end
                        elseif elementType==4
                            append!(tetrahedra,[nodeTags])
                            for dom in ent2dom[entityDim+1][entityTag]
                                append!(domains[dom]["simplices"],length(tetrahedra))
                            end
                        end
                    end
                end
            end
            #find closing tag and read it
            #while line[2:end]!="End"*fldname
            #    line=readline(fid)
            #end
        end#while
    end#do
    return points, lines, triangles, tetrahedra, domains
end#function

"""
    points, lines, triangles, tetrahedra, domains=read_msh2(file_name)

Read gmsh's .msh-format version 2. This is an old format and wherever possible version 4 should be used.

http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
"""
function read_msh2(file_name)
    lines=Array{UInt32,1}[]
    triangles=Array{UInt32,1}[]
    tetrahedra=Array{UInt32,1}[]
    points=[]
    domains=Dict()
    open(file_name) do f
        field=[]
        field_open=false
        i=1
        #global points, triangles, lines, tetrahedra, domains
        local field_name, domain_map
        while !eof(f)
            line=readline(f)
            if (line[1]=='$')
                if field_open==false
                    field_name=line[2:end]
                    field_open=true
                elseif field_open && (line[2:end]=="End"*field_name)
                   field_open=false
                else
                   println("Error:Field is invalid")
                   break
                end
            end
            if field_open
                if field_name=="MeshFormat"
                    MeshFormat=readline(f)
                elseif field_name=="PhysicalNames"
                    #get number of field entries
                    number_of_entries=parse(Int,readline(f))
                    domain_map=Dict()
                    for i=1:number_of_entries
                        line=split(readline(f))
                        line[3]=strip(line[3],['\"'])
                        domains[line[3]]=Dict("dimension"=>parse(UInt32,line[1])
                        ,"points"=>UInt32[])
                        domain_map[parse(UInt32,line[2])]=line[3]
                    end
                elseif field_name=="Nodes"

                    #get number of field entries
                    number_of_entries=parse(Int,readline(f))
                    points=Array{Float64}(undef,3,number_of_entries)
                        for i=1:number_of_entries
                            line=split(readline(f))
                            points[:,i]=[parse(Float64,line[2]) parse(Float64,line[3]) parse(Float64,line[4])]
                        end
                elseif field_name=="Elements"
                    #get number of field entries
                    number_of_entries=parse(Int,readline(f))
                    for i=1:number_of_entries
                        line=split(readline(f))
                        element_type=parse(Int16,line[2])
                        number_of_tags=parse(Int16,line[3])
                        tags=line[4:4+number_of_tags-1]
                        node_number_list=UInt32[]

                        for node in line[4+number_of_tags:end]
                            append!(node_number_list,parse(UInt32,node))
                        end

                        dom=domain_map[parse(UInt32,tags[1])]
                        for node in node_number_list
                            insert_smplx!(domains[dom]["points"],node)
                        end

                        if element_type==4 #4-node tetrahedron
                        insert_smplx!(tetrahedra,node_number_list)
                        elseif element_type==2 #3-node triangle
                            insert_smplx!(triangles,node_number_list)
                        elseif element_type==1 #2-node line
                            insert_smplx!(lines,node_number_list)
                        #elseif element type==15 #1-node point
                        #    insert_smplx!(      ,node_number_list)
                        else
                            println("error, element type is not defined")
                            return
                        end
                    end
                end
            end
            i+=1
        end
    end
    #triangles=Array{UInt32,1}[]#hot fix for triangle bug
    return points, lines, triangles, tetrahedra, domains
end

include("./Mesh/read_nastran.jl")

##

##

##
"""
    link_triangles_to_tetrahedra!(mesh::Mesh)

Find the tetrahedra that are connected to the triangles in mesh.triangles (and mesh.inttriangles) and store this information in mesh.tri2tet.
"""
function link_triangles_to_tetrahedra!(mesh::Mesh)
    if !isempty(mesh.int_triangles)
        #IDEA: move this to constructor
        for (idt,tet) in enumerate(mesh.tetrahedra)
            for tri in (tet[[1,2,3]],tet[[1,2,4]],tet[[1,3,4]],tet[[2,3,4]])
                idx,ins=sort_smplx(mesh.triangles,tri)
                if !ins
                    mesh.tri2tet[1][idx]=idt
                else
                    idx,ins=sort_smplx(mesh.int_triangles,tri)
                    if !ins
                        if mesh.tri2tet[2][idx,1]==0
                            mesh.tri2tet[2][idx,1]=idt
                        else
                            mesh.tri2tet[2][idx,2]=idt
                        end
                    end
                end
            end
        end

    else
        for (idt,tet) in enumerate(mesh.tetrahedra)
            for tri in (tet[[1,2,3]],tet[[1,2,4]],tet[[1,3,4]],tet[[2,3,4]])
                idx,ins=sort_smplx(mesh.triangles,tri)
                if !ins
                    mesh.tri2tet[idx]=idt
                end
            end
        end
    end
    return nothing
end
##
function assemble_triangles(tetrahedra)
    int_triangles=Array{UInt32,1}[]
    ext_triangles=Array{UInt32,1}[]
    for tet in tetrahedra#loop over all tetrahedra
        for tri_idcs in [[1,2,3],[1,2,4],[1,3,4],[2,3,4]] #loop over all 4 triangles of the tetrahedron
            tri=tet[tri_idcs]
            idx,notinside=sort_smplx(ext_triangles, tri)
            if notinside #triangle occurs for the first time
                #IDEA: there is no sanity check. If tetrahedra are ill defined a
                #triangle might show up three times. That would be fatal.
                #may implement some low-level-sanity checks ensuring data
                #integrity

                insert!(ext_triangles,idx,tri)
            else #triangle occurs for the second time, it must be an interior triangle.
                deleteat!(ext_triangles,idx)
                insert_smplx!(int_triangles,tri)
            end
        end
    end
    return ext_triangles,int_triangles
end
##
"""
    new_mesh=octosplit(mesh::Mesh)

Subdivide each tetrahedron in the mesh `mesh` into 8 tetrahedra by splitting
each edge at its center.

# Notes
The algorithm introduces `length(mesh.lines)` new vertices. This yields a finer
mesh featuring `8*size(mesh,2)` tetrahedra, `4*length(mesh.triangles)`
triangles, and `2*length(mesh.lines)` lines. From the  3 possible subdivision of
a tetrahedron the algorithm automatically chooses the one that minimizes the
edge lengths of the new tetrahedra. The point labeling of `new_mesh` is
consistent with the point labeling in `mesh`, i.e., the first
`size(mesh.points,2)` points in `new_mesh.points` are identical to the points
in `mesh.points`. Hence, `mesh` and `new_mesh` form a hierachy of meshes.
"""
function octosplit(mesh::Mesh)
    collect_lines!(mesh)
    N_points=size(mesh.points,2)
    N_lines=length(mesh.lines)
    N_tet=length(mesh.tetrahedra)
    N_tri=length(mesh.triangles)
    ##
    points=Array{Float64,2}(undef,3,N_points+N_lines)
    ## extend point list
    points[:,1:N_points]=mesh.points
    for (idx,ln) in enumerate(mesh.lines)
        points[:,N_points+idx]= sum(mesh.points[:,ln],dims=2)./2
    end
    #create new tetrahedra by splitting the old ones
    tetrahedra=Array{UInt32}[]
    for tet in mesh.tetrahedra
        A,B,C,D=tet #old vertices
        #indices of the center points of the lines also become vertices
        AB=N_points+find_smplx(mesh.lines,[A,B])
        AC=N_points+find_smplx(mesh.lines,[A,C])
        AD=N_points+find_smplx(mesh.lines,[A,D])
        BC=N_points+find_smplx(mesh.lines,[B,C])
        BD=N_points+find_smplx(mesh.lines,[B,D])
        CD=N_points+find_smplx(mesh.lines,[C,D])

        #outer tetrahedra always present...
        insert_smplx!(tetrahedra,[A, AB, AC, AD])
        insert_smplx!(tetrahedra,[B, AB, BC, BD])
        insert_smplx!(tetrahedra,[C, AC, BC, CD])
        insert_smplx!(tetrahedra,[D, AD, BD, CD])

        # split remaining octet based on minimum distance
        AB_CD=LinearAlgebra.norm(points[:,AB].-points[:,CD])
        AC_BD=LinearAlgebra.norm(points[:,AC].-points[:,BD])
        AD_BC=LinearAlgebra.norm(points[:,AD].-points[:,BC])

        if AB_CD<=AC_BD && AB_CD<=AD_BC
            insert_smplx!(tetrahedra,[AB,CD,AC,AD])
            insert_smplx!(tetrahedra,[AB,CD,AD,BD])
            insert_smplx!(tetrahedra,[AB,CD,BD,BC])
            insert_smplx!(tetrahedra,[AB,CD,BC,AC])
        elseif AC_BD<=AB_CD && AC_BD<=AD_BC
            insert_smplx!(tetrahedra,[AC,BD,AB,AD])
            insert_smplx!(tetrahedra,[AC,BD,AD,CD])
            insert_smplx!(tetrahedra,[AC,BD,CD,BC])
            insert_smplx!(tetrahedra,[AC,BD,BC,AB])
        elseif AD_BC<=AC_BD && AD_BC<=AB_CD
            insert_smplx!(tetrahedra,[AD,BC,AC,CD])
            insert_smplx!(tetrahedra,[AD,BC,CD,BD])
            insert_smplx!(tetrahedra,[AD,BC,BD,AB])
            insert_smplx!(tetrahedra,[AD,BC,AB,AC])
        end
    end

    #split the old surface triangles
    triangles=Array{UInt32}[]
    for tri in mesh.triangles
        A,B,C=tri
        AB=N_points+find_smplx(mesh.lines,[A,B])
        AC=N_points+find_smplx(mesh.lines,[A,C])
        BC=N_points+find_smplx(mesh.lines,[B,C])
        insert_smplx!(triangles,[A,AB,AC])
        insert_smplx!(triangles,[B,AB,BC])
        insert_smplx!(triangles,[C,AC,BC])
        insert_smplx!(triangles,[AB,AC,BC])
    end

    #split the old interior triangles
    #TODO:Testing
    inttriangles=Array{UInt32}[]
    # for tri in mesh.int_triangles
    #     A,B,C=tri
    #     AB=N_points+find_smplx(mesh.lines,[A,B])
    #     AC=N_points+find_smplx(mesh.lines,[A,C])
    #     BC=N_points+find_smplx(mesh.lines,[B,C])
    #     insert_smplx!(inttriangles,[A,AB,AC])
    #     insert_smplx!(inttriangles,[B,AB,BC])
    #     insert_smplx!(inttriangles,[C,AC,BC])
    #     insert_smplx!(inttriangles,[AB,AC,BC])
    # end

    ## relabel tetrahedra
    tet_labels=Array{UInt32}(undef,N_tet,8)
    for (idx,tet) in enumerate(mesh.tetrahedra)
        A,B,C,D=tet #old vertices
        #indices of the center points of the lines also become vertices
        AB=N_points+find_smplx(mesh.lines,[A,B])
        AC=N_points+find_smplx(mesh.lines,[A,C])
        AD=N_points+find_smplx(mesh.lines,[A,D])
        BC=N_points+find_smplx(mesh.lines,[B,C])
        BD=N_points+find_smplx(mesh.lines,[B,D])
        CD=N_points+find_smplx(mesh.lines,[C,D])

        #outer tetrahedra always present...
        tet_labels[idx,1]=find_smplx(tetrahedra,[A, AB, AC, AD])
        tet_labels[idx,2]=find_smplx(tetrahedra,[B, AB, BC, BD])
        tet_labels[idx,3]=find_smplx(tetrahedra,[C, AC, BC, CD])
        tet_labels[idx,4]=find_smplx(tetrahedra,[D, AD, BD, CD])

        # split remaining octet based on minimum distance
        AB_CD=LinearAlgebra.norm(points[:,AB].-points[:,CD])
        AC_BD=LinearAlgebra.norm(points[:,AC].-points[:,BD])
        AD_BC=LinearAlgebra.norm(points[:,AD].-points[:,BC])

        if AB_CD<=AC_BD && AB_CD<=AD_BC
            tet_labels[idx,5]=find_smplx(tetrahedra,[AB,CD,AC,AD])
            tet_labels[idx,6]=find_smplx(tetrahedra,[AB,CD,AD,BD])
            tet_labels[idx,7]=find_smplx(tetrahedra,[AB,CD,BD,BC])
            tet_labels[idx,8]=find_smplx(tetrahedra,[AB,CD,BC,AC])
        elseif AC_BD<=AB_CD && AC_BD<=AD_BC
            tet_labels[idx,5]=find_smplx(tetrahedra,[AC,BD,AB,AD])
            tet_labels[idx,6]=find_smplx(tetrahedra,[AC,BD,AD,CD])
            tet_labels[idx,7]=find_smplx(tetrahedra,[AC,BD,CD,BC])
            tet_labels[idx,8]=find_smplx(tetrahedra,[AC,BD,BC,AB])
        elseif AD_BC<=AC_BD && AD_BC<=AB_CD
            tet_labels[idx,5]=find_smplx(tetrahedra,[AD,BC,AC,CD])
            tet_labels[idx,6]=find_smplx(tetrahedra,[AD,BC,CD,BD])
            tet_labels[idx,7]=find_smplx(tetrahedra,[AD,BC,BD,AB])
            tet_labels[idx,8]=find_smplx(tetrahedra,[AD,BC,AB,AC])
        end
    end

    #relabel triangles
    tri_labels=Array{UInt32}(undef,N_tri,4)
    for (idx,tri) in enumerate(mesh.triangles)
        A,B,C=tri
        AB=N_points+find_smplx(mesh.lines,[A,B])
        AC=N_points+find_smplx(mesh.lines,[A,C])
        BC=N_points+find_smplx(mesh.lines,[B,C])
        tri_labels[idx,1]=find_smplx(triangles,[A,AB,AC])
        tri_labels[idx,2]=find_smplx(triangles,[B,AB,BC])
        tri_labels[idx,3]=find_smplx(triangles,[C,AC,BC])
        tri_labels[idx,4]=find_smplx(triangles,[AB,AC,BC])
    end

    #reassemble domains
    domains=Dict()
    for (dom,domain) in mesh.domains
        domains[dom]=Dict()
        dim=domain["dimension"]
        domains[dom]["dimension"]=dim
        #TODO: detect whether there optional field like "size"
        simplices=[]
        for smplx_idx in domain["simplices"]
            if dim==3
                append!(simplices,tet_labels[smplx_idx,:])
            elseif dim==2
                append!(simplices,tri_labels[smplx_idx,:])
            end
        end
        domains[dom]["simplices"]=sort(simplices)
    end

    tri2tet=zeros(UInt32,length(triangles))
    if length(tri2tet)>0
        tri2tet[1]=0xffffffff
    end #magic sentinel value to check whether list is linked
    return Mesh(mesh.file, points, [], triangles, inttriangles, tetrahedra, domains, mesh.file,tri2tet,1)
end

##
"""
    compute_size!(mesh::Mesh,dom::String)

Compute the size of the domain `dom` and store it in the field `"size"` of the domain definition.
The size will be a length, area, or volume depending on the dimension of the domain.
"""
function compute_size!(mesh::Mesh,dom::String)
    dim=mesh.domains[dom]["dimension"]
    V=0 #TODO: check whether domains are properly defined
    if dim ==3
        for tet in mesh.tetrahedra[mesh.domains[dom]["simplices"]]
            tet=mesh.points[:,tet]
            V+=abs(LinearAlgebra.det(tet[:,1:3]-tet[:,[4,4,4]]))/6
        end
    #mesh.domains[dom]["volume"]=V
    elseif dim == 2
        for tri in mesh.triangles[mesh.domains[dom]["simplices"]]
            tri=mesh.points[:,tri]
            V+=LinearAlgebra.norm(LinearAlgebra.cross(tri[:,1]-tri[:,3],tri[:,2]-tri[:,3]))/2
        end
    #mesh.domains[dom]["area"]=V
    elseif dim == 1
        for ln in in mesh.lines[mesh.domains[dom]["simplices"]]
            ln=mesh.points[:,ln]
            V+=LinearAlgebra.norm(ln[:,2]-ln[:,1])
        end

    end
    mesh.domains[dom]["size"]=V
end
##
# function get_simplices(mesh::Mesh, dom::String)
#     dim=mesh.domains[dom]["dimension"]
#     if dim==3
#         return mesh.tetrahedra[mesh.domains[dom]["simplices"]]
#     elseif dim==2
#         return mesh.triangles[mesh.domains[dom]["simplices"]]
#     end
# end
##

"""
    tet_idx=find_tetrahedron_containing_point(mesh::Mesh,point)

Find the tetrahedron containing the point `point` and return the index
of the tetrahedron as `tet_idx`. If the point lies on the interface of two or more tetrahedra, the returned `tet_idx` will be the lowest index of these tetrahedra,
i.e. the index depends on the ordering of the tetrahedra in the `mesh.tetrahedra`.
If no tetrahedron in the mesh encloses the point `tet_idx==0`.
"""
function find_tetrahedron_containing_point(mesh::Mesh,point)
    tet_idx=0
    for (idx::UInt32,tet) in enumerate(mesh.tetrahedra)
        tet=mesh.points[:,tet]
        J=tet[:,1:3]-tet[:,[4,4,4]]
        #convert to local barycentric coordinates
        ξ=J\(point-tet[:,4])
        ξ=vcat(ξ,1-sum(ξ))
        if all(0 .<= ξ .<= 1)
            tet_idx=idx
            break
        end
        #TODO impleemnt points on surfaces
    end
    return tet_idx
end

##
function keep!(mesh::Mesh,keys_to_be_kept)
    for key in keys(mesh.domains)
        if !(key in keys_to_be_kept)
            delete!(mesh.domains,key)
        end
    end
end

"""
    collect_lines!(mesh::Mesh)

Populate the list of lines in `mesh.lines`.
"""
function collect_lines!(mesh::Mesh)
    for tet in mesh.tetrahedra
        insert_smplx!(mesh.lines,[tet[1],tet[2]])
        insert_smplx!(mesh.lines,[tet[1],tet[3]])
        insert_smplx!(mesh.lines,[tet[1],tet[4]])
        insert_smplx!(mesh.lines,[tet[2],tet[3]])
        insert_smplx!(mesh.lines,[tet[2],tet[4]])
        insert_smplx!(mesh.lines,[tet[3],tet[4]])
    end
end

"""
    unify!(mesh::Mesh ,new_domain::String,domains...)

Create union of `domains` and name it `new_domain`.
`domains` must share the same dimension and `new_domain`
must be a new name otherwise errors will occur.
"""
function unify!(mesh::Mesh ,new_domain::String,domains...)
    #check new domain is non-existent
    if new_domain in keys(mesh.domains)
        println("ERROR: Domain '"*new_domain*"' already exists!\n Use a different name.")
        return
    end
    #check dimension
    dim=mesh.domains[domains[1]]["dimension"]
    for dom in domains
        if dim != mesh.domains[dom]["dimension"]
            println("Error: Domains are of different dimension!")
            return
        end
    end

    mesh.domains[new_domain]=Dict{String,Any}()
    mesh.domains[new_domain]["dimension"]=dim
    mesh.domains[new_domain]["simplices"]=UInt32[]
    for dom in domains
        #mesh.domains[new_domain]["simplices"]=
        union!(mesh.domains[new_domain]["simplices"],mesh.domains[dom]["simplices"])
    end
    sort!(mesh.domains[new_domain]["simplices"])
end

"""
    surface_points, tri_mask, tet_mask = get_surface_points(mesh::Mesh,output=true)

Get a list `surface_points` of all point indices that are on the surface of the
mesh `mesh`. The corresponding lists `tri_mask and `tet_mask` contain lists of
indices of all triangles and tetrahedra that are connected to the surface
point. The optional parameter `output` toggles whether a progressbar should be
shown or not.
"""
function get_surface_points(mesh::Mesh,output::Bool=true)
    #TODO put return values as fields into Mesh
    #step #1 find all points on the surface
    surface_points=[]
    #intialize progressbar
    if output
        p = ProgressMeter.Progress(length(mesh.triangles),desc="Find points #1... ", dt=1,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)
    end
    for tri in mesh.triangles
        for pnt in tri
            insert_smplx!(surface_points,pnt)
        end
        #update progressbar
        if output
            ProgressMeter.next!(p)
        end
    end

    #step 2 link points to simplices
    #tri_mask=Array{Int64,1}[]
    #tet_mask=Array{Int64,1}[]
    tri_mask=[[] for idx = 1:length(surface_points)]
    tet_mask=[[] for idx = 1:length(surface_points)]
    #intialize progressbar
    if output
        p = ProgressMeter.Progress(length(mesh.triangles),desc="Link to triangles... ", dt=1,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)
    end

    for (tri_idx,tri) in enumerate(mesh.triangles)
        for pnt_idx in tri
            idx=find_smplx(surface_points,pnt_idx)
            if idx>0
                append!(tri_mask[idx],tri_idx)
            end
        end
        #update progressbar
        if output
            ProgressMeter.next!(p)
        end
    end
    if output
        p = ProgressMeter.Progress(length(mesh.tetrahedra),desc="Link to tetrahedra... ", dt=1,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)
    end
    for (tet_idx,tet) in enumerate(mesh.tetrahedra)
        for pnt_idx in tet
            idx=find_smplx(surface_points,pnt_idx)
            if idx>0
                append!(tet_mask[idx],tet_idx)
            end
        end
        #update progressbar
        if output
            ProgressMeter.next!(p)
        end
    end

    if mesh.dos!=1 && mesh.dos.unit
        #blochify surface points
        for idx in surface_points
            bidx=idx -mesh.dos.naxis
            #reference surface
            if 0 < bidx <= mesh.dos.nxbloch
                append!(tri_mask[idx],tri_mask[end-mesh.dos.nxbloch+bidx])
                append!(tet_mask[idx],tet_mask[end-mesh.dos.nxbloch+bidx])
                #unique!(tri_mask)
                #unique!(tet_mask)
                unique!(tri_mask[idx])
                unique!(tet_mask[idx])
                #symmetrize
                append!(tri_mask[end-mesh.dos.nxbloch+bidx],tri_mask[idx])
                append!(tet_mask[end-mesh.dos.nxbloch+bidx],tet_mask[idx])
                unique!(tri_mask[end-mesh.dos.nxbloch+bidx])
                unique!(tet_mask[end-mesh.dos.nxbloch+bidx])
            end
        end
    end

    return surface_points, tri_mask, tet_mask
end
# function get_surface_points(mesh::Mesh,output::Bool=true)
#     #TODO put return values as fields into Mesh
#     #step #1 find all points on the surface
#     surface_points=[]
#     #intialize progressbar
#     if output
#         p = ProgressMeter.Progress(length(mesh.triangles),desc="Find points #1... ", dt=1,
#              barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
#              barlen=20)
#     end
#     for tri in mesh.triangles
#         for pnt in tri
#             insert_smplx!(surface_points,pnt)
#         end
#         #update progressbar
#         if output
#             ProgressMeter.next!(p)
#         end
#     end
#
#     #step 2 link points to simplices
#     tri_mask=Array{Int64,1}[]
#     tet_mask=Array{Int64,1}[]
#     #intialize progressbar
#     if output
#         p = ProgressMeter.Progress(length(surface_points),desc="Link to simplices... ", dt=1,
#              barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
#              barlen=20)
#     end
#     for (idx,pnt_idx) in enumerate(surface_points)
#         mask=Int64[]
#         for (tri_idx,tri) in enumerate(mesh.triangles)
#             if pnt_idx in tri
#                 append!(mask,tri_idx)
#             end
#         end
#         append!(tri_mask,[mask])
#         mask=Int64[]
#         for (tet_idx,tet) in enumerate(mesh.tetrahedra)
#             if pnt_idx in tet
#                 append!(mask,tet_idx)
#             end
#         end
#         append!(tet_mask, [mask])
#         #update progressbar
#         if output
#             ProgressMeter.next!(p)
#         end
#     end
#     return surface_points, tri_mask, tet_mask
# end

"""
    normal_vectors=get_normal_vectors(mesh::Mesh,output::Bool=true)

Compute a 3×`length(mesh.triangles)` Array containing the outward pointing
normal vectors of each of the surface triangles of the mesh `mesh`. The vectors
are not normalised but their norm is twice the area of the corresponding
triangle. The optional parameter `output` toggles whether a progressbar should be
shown or not.
"""
function get_normal_vectors(mesh::Mesh,output::Bool=true)
    #TODO put return values as fields into Mesh
    #initialize progressbar
    if output
        p = ProgressMeter.Progress(length(mesh.triangles),desc="Find points #1... ", dt=1,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)
    end
    normal_vectors=Array{Float64}(undef,3,length(mesh.triangles))
    #Identify points A,B,C, and D that form the tetrahedron connected to the
    #surface triangle. D is not part of the surface but is inside the body
    #therefore its positio can be used to find determine the outward direction.
    if length(mesh.tri2tet)!=0 && mesh.tri2tet[1]==0xffffffff #TODO: Initialize as empty instead with sentinel value
        link_triangles_to_tetrahedra!(mesh)
    end
    for (idx,tri) in enumerate(mesh.triangles)
        A,B,C=tri
        tet=mesh.tri2tet[idx]
        tet=mesh.tetrahedra[tet]
        D=0
        for idx in tet
            if !(idx in tri)
                D=idx
                break
            end
        end
        A=mesh.points[:,A]
        B=mesh.points[:,B]
        C=mesh.points[:,C]
        D=mesh.points[:,D]
        N=LinearAlgebra.cross(A.-C,B.-C) #compute normal vector (length is twice the area of teh surface triangle)
        N*=sign(LinearAlgebra.dot(N,C.-D))# make normal outward
        normal_vectors[:,idx]=N
    end
    #update progressbar
    if output
        ProgressMeter.next!(p)
    end
    return normal_vectors
end

"""
    generate_field(mesh::Mesh,func;el_type=0)

Generate field from function `func`for mesh `mesh`. The element type is either
`el_type=0` for field values associated with the mesh tetrahedra or `el_type=1`
for field values associated with the mesh vertices. The function must accept
three input arguments corresponding to the three space dimensions.
"""
function generate_field(mesh::Mesh,func;order::Symbol=:const)
    #TODO: implement el_type=2
    if order==:const
        field=zeros(Float64,length(mesh.tetrahedra))#*347.0
        for idx=1:length(mesh.tetrahedra)
            cntr=sum(mesh.points[:,mesh.tetrahedra[idx]],dims=2)/4 #centerpoint of a tetrahedron is arithmetic mean of the vertices
            field[idx]=func(cntr...)
        end
    elseif order==:lin
        field=zeros(Float64,size(mesh.points,2))
        for idx in 1:size(mesh.points,2)
            pnt=mesh.points[:,idx]
            field[idx]=func(pnt...)
        end
    else
        println("Error: Can't generate field order $order not supported!")
        return nothing
    end
    return field
end
##legacycode
# function generate_field(mesh::Mesh,func;el_type=0)
#     #TODO: implement el_type=2
#     if el_type==0
#         field=zeros(Float64,length(mesh.tetrahedra))#*347.0
#         for idx=1:length(mesh.tetrahedra)
#             cntr=sum(mesh.points[:,mesh.tetrahedra[idx]],dims=2)/4 #centerpoint of a tetrahedron is arithmetic mean of the vertices
#             field[idx]=func(cntr...)
#         end
#     else
#         field=zeros(Float64,size(mesh.points,2))
#         for idx in 1:size(mesh.points,2)
#             pnt=mesh.points[:,idx]
#             field[idx]=func(pnt...)
#         end
#     end
#     return field
# end
"""
    data,surf_keys,vol_keys=color_domains(mesh::Mesh,domains=[])

Create data fields containing an integer number corresponding to the local
domain.

# Arguments
- `mesh::Mesh`: mesh to be colored
- `domains::List`: (optional) list of domains that are to be colored. If empty  all domains will be colored.

# Returns
- `data::Dict`: Dictionary containing a key for each colored domain plus magic keys `"__all_surfaces__"` and `"__all_volumes__"` holding all colors in one field. These magic fields may not work correctly if the domains are not disjoint.
- `surf_keys::Dict`: Dictionary mapping the surface domain names to their indeces.
- `vol_keys::Dict`: Dictionary mapping the volume domain names to their indeces.

# Notes
The `data`variable is designed to flawlessly work with `vtk_write` for
visualization with paraview. Check the "Interpret Values as Categories" box in
paraview's color map editor for optimal visualization.

See also: [`vtk_write`](@ref)
"""
function color_domains(mesh::Mesh,domains=[])
    number_of_triangles=length(mesh.triangles)
    number_of_tetrahedra=length(mesh.tetrahedra)
    #IDEA: check number of domains and choose datatype accordingly (UInt8 might be possible)
    tri_color=zeros(UInt16,number_of_triangles)
    tet_color=zeros(UInt16,number_of_tetrahedra)
    data=Dict()
    surf_keys=Dict()
    vol_keys=Dict()
    if length(domains)==0
        domains=sort([key for key in keys(mesh.domains)])
    end
    surf_idx=0
    vol_idx=0
    for key in domains
        if key in keys(mesh.domains)
            dom=mesh.domains[key]
        else
            println("Warning: No domain named '$key' in mesh.")
            continue
        end
        smplcs=dom["simplices"]
        dim=dom["dimension"]
        if dim==2
            surf_idx+=1

            if any(tri_color[smplcs].!=0)
                println("domain $key is overlapping")
            end
            tri_color[smplcs].=surf_idx
            data[key]=zeros(Int64,number_of_triangles)
            data[key][smplcs].=surf_idx
            surf_keys[key]=surf_idx
            println("surface:$key-->$surf_idx")
        elseif dim==3
            vol_idx+=1
            if any(tet_color[smplcs].!=0)
                println("domain $key is overlapping")
            end
            tet_color[smplcs].=vol_idx
            data[key]=zeros(Int64,number_of_tetrahedra)
            data[key][smplcs].=vol_idx
            vol_keys[key]=vol_idx
            println("volume:$key-->$vol_idx")
        end
    end

    data["__all_surfaces__"]=tri_color
    data["__all_volumes__"]=tet_color

 return data,surf_keys,vol_keys
end

#################################################
#! legacy functions will be removed in the future
#!################################################
function build_domains!(mesh::Mesh)
    for dom in keys(mesh.domains)
        dim=mesh.domains[dom]["dimension"]
        mesh.domains[dom]["simplices"]=UInt32[]
        if dim == 3
            for (idx::UInt32,tet) in enumerate(mesh.tetrahedra)
                count=0
                for node in tet
                    count+=node in mesh.domains[dom]["points"]
                end
                if count==4
                    append!(mesh.domains[dom]["simplices"],idx)
                end
            end
        elseif dim ==2
            for (idx::UInt32,tri) in enumerate(mesh.triangles)
                count=0
                for node in tri
                    count+=node in mesh.domains[dom]["points"]
                end
                if count==3
                    append!(mesh.domains[dom]["simplices"],idx)
                end
            end
        end
    end
end

function find_tetrahedron_containing_triangle!(mesh::Mesh,dom::String)
    tetmap=Array{UInt32,1}(undef,length(mesh.domains[dom]["simplices"]))
    for (idx,tri) in enumerate(mesh.domains[dom]["simplices"])
        tri=mesh.triangles[tri]
        for (idt::UInt32,tet) in enumerate(mesh.tetrahedra)
            if all([i in tet for i=tri])
                tetmap[idx]=idt
                break
            end
        end
    end
    mesh.domains[dom]["tetmap"]=tetmap
    return nothing
end

""" read a mesh from ansys cfx file"""
function read_ansys(filename)
    lines=Array{UInt32,1}[]
    triangles=Array{UInt32,1}[]
    names=Dict()
    local points,domains,tetrahedra
    #points=[]
    open(filename) do f
        domains=Dict()
        while !eof(f)
            line=readline(f)
            splitted=split(line)
            if line==""
                continue
            end
            if splitted[1]=="(10" && splitted[2]=="(0"
                 #point declaration
                 number_of_points=parse(Int,splitted[4],base=16)
                 points=zeros(3,number_of_points)
                 pnt_idx=1
            elseif splitted[1]=="(10" && splitted[2]!="(0"
                #point field header
                pnt_field=true
                zone_id=splitted[2][2:end]
                first_index=parse(Int,splitted[3],base=16)
                last_index=parse(Int,splitted[4],base=16)
                type=splitted[5]
                #read point field
                if length(splitted)==6
                    ND=parse(UInt32,splitted[6][1:end-2],base=16)
                else
                    ND=-1
                end
                domains[zone_id]=Dict("dimension"=>ND,"points"=>UInt32[])
                for idx=first_index:last_index
                    line=readline(f)
                    splitted=split(line)
                    points[:,idx]=parse.(Float64,splitted)
                    append!(domains[zone_id]["points"],idx)
                end

            elseif splitted[1]=="(12" && splitted[2]=="(0"
                #cell declaration
                number_of_tetrahedra=parse(Int,splitted[4],base=16)
                tetrahedra=Array{Array{UInt32,1}}(undef,number_of_tetrahedra)
                for idx in 1:number_of_tetrahedra
                    tetrahedra[idx]=[]
                end

            elseif splitted[1]=="(13" && splitted[2]=="(0"
                #face declaration
                number_of_faces=parse(Int,splitted[4],base=16)
            elseif splitted[1]=="(13" && splitted[2]!="(0"
                #face field header
                zone_id=splitted[2][2:end]
                first_index=parse(Int,splitted[3],base=16)
                last_index=parse(Int,splitted[4],base=16)
                type=splitted[5]
                element_type=parse(Int,splitted[6][1:end-2],base=16)
                if type=="2"
                    ND=3
                else
                    ND=2
                end
                domains[zone_id]=Dict("dimension"=>ND,"points"=>UInt32[])

                #read face field
                for idx =first_index:last_index
                    line=readline(f)
                    splitted=split(line)
                    left=parse(UInt32,splitted[end-1],base=16)
                    right=parse(UInt32,splitted[end],base=16)
                    tri=parse.(UInt32,splitted[1:3],base=16)
                    for pnt in tri
                        #Hier können die Domain simplices abgefragt werden
                        if left!=0
                            insert_smplx!(tetrahedra[left],pnt)
                        end
                        if right!=0
                            insert_smplx!(tetrahedra[right],pnt)
                        end
                        insert_smplx!(domains[zone_id]["points"],pnt)
                    end
                    if left==0 || right==0
                        insert_smplx!(triangles,tri)
                        #Das muss Rückverfolgt werden um die Dreiecke zu identifizieren
                    end
                end
            elseif splitted[1]=="(45"
                names[splitted[2][2:end]]=splitted[end][1:end-4]
            end
        end
    end
    for (key,value) in names
        if key in keys(domains)
            domains[value]=pop!(domains,key)
        else
            #TODO: warning message
        end

    end
    #points, lines, triangles, tetrahedra, domains, names
    return Mesh(filename, points, lines, triangles, tetrahedra, domains, filename,1)
end

end#module


text:
# TODO relies on implementation details of BangBang for pure versions (_setindex)
module SetfieldImpl

using BangBang: BangBang, NoBang, implements
using DocStringExtensions

using Setfield: Setfield
using Setfield: ComposedLens, ConstIndexLens, DynamicIndexLens, IndexLens
using Setfield: Lens, PropertyLens, set

const SupportedIndexLens = Union{ConstIndexLens,DynamicIndexLens,IndexLens}

struct Lens!{T<:Lens} <: Lens
    lens::T
end

Setfield.get(obj, lens::Lens!) = get(obj, lens.lens)

function Setfield.set(obj, lens::Lens!, value)
    if issettable(obj, lens.lens)
        _set(obj, lens, value)
    else
        error("At least one of the applied lens must correspond to a mutable object")
    end
end

function _set(obj, lens::Lens!{<:ComposedLens}, value)
    inner_lens, outer_lens = Lens!(lens.lens.inner), Lens!(lens.lens.outer)
    inner_obj = get(obj, outer_lens)
    inner_val = _set(inner_obj, inner_lens, value)
    _set(obj, outer_lens, inner_val)
end

function _set(obj, ::Lens!{<:PropertyLens{name}}, value) where {name}
    if implements(setproperty!, obj)
        setproperty!(obj, name, value)
        return obj
    else
        return NoBang.setproperty(obj, name, value)
    end
end

indicesfor(lens::IndexLens, _) = lens.indices
indicesfor(::ConstIndexLens{I}, _) where {I} = I
indicesfor(lens::DynamicIndexLens, obj) = lens.f(obj)

function _set(obj, lens::Lens!{<:SupportedIndexLens}, value)
    if implements(setindex!, obj)
        setindex!(obj, value, indicesfor(lens.lens, obj)...)
        return obj
    else
        return NoBang._setindex(obj, value, indicesfor(lens.lens, obj)...)
    end
end

"""
    $(TYPEDSIGNATURES)

See also [`@set!`](@ref).
"""
prefermutation(lens::Lens) = Lens!(lens)

"""
    $(TYPEDSIGNATURES)

Like `Setfield.@set`, but always mutate when at least one of the nested objects is
mutable. Otherwise fails.

# Examples
```jldoctest
julia> using LyceumBase

julia> mutable struct Mutable
           a
           b
       end

julia> x = orig = Mutable((x=Mutable(1, 2), y=3), 4);

julia> @set! x.a.x.a = 10;

julia> @assert x.a.x.a == orig.a.x.a == 10

julia> immutable = (x=(y=(z=1,),),)

julia> # this would error: @set! y.x.a = 10;
```
"""
macro set!(ex)
    Setfield.setmacro(prefermutation, ex, overwrite = false)
end

Base.@propagate_inbounds issettable(obj, ::PropertyLens) = implements(setproperty!, obj)
Base.@propagate_inbounds issettable(obj, ::SupportedIndexLens) = implements(setindex!, obj)
Base.@propagate_inbounds function issettable(obj, lens::ComposedLens)
    issettable(obj, lens.outer) || issettable(get(obj, lens.outer), lens.inner)
end

end


text:
module Padre

using JSON

export Subject,Session,Dataset

export padre_subjects
export subj_dir,tags,labels,sessions
export sess_dir,dsets,filename

## Config

padre_dir = ENV["PADRE_ROOT"]
subj_dir(subjid::String) = "$padre_dir/Data/$subjid"
json_file(subjid::String) = "$(subj_dir(subjid))/$subjid.json"

## Objects

mutable struct Subject
    subjid::String
    json::Dict
end

struct Session
    subj::Subject
    session::String
end

mutable struct Dataset
    filename::String
    meta::Dict
    sess::Session
end

import Base.show
function show(io::IO,s::Subject)
    write(io,"Subject(\"$(s.subjid)\")")
    return nothing
end

## Global functions

"""
    padre_subjects()

Return list of all available subjects as `Subject` objects."""
function padre_subjects()
    subjs = Subject[]

    for s in readdir(joinpath(padre_dir,"Data"))
        try
            subj = Subject(s)
            push!(subjs,subj)
        catch
        end
    end
    subjs
end

## Make a Subject

Subject(subjid::String) = Subject(subjid,JSON.parsefile(json_file(subjid)))

## Access functions

sess_dict(sess::Session) = sess.subj.json["sessions"][sess.session]
function tags(sess::Session)
    sd = sess_dict(sess)
    return haskey(sd,"tags") ? sd["tags"] : []
end
function labels(sess::Session)
    sd = sess_dict(sess)
    return haskey(sd,"labels") ? sd["labels"] : []
end

function sessions(s::Subject;tag=nothing,label=nothing,include::Bool=true)
    ss = Session[]
    for session in keys(s.json["sessions"])
        if (include == false) || (haskey(s.json["sessions"][session],"include") && s.json["sessions"][session]["include"])
            push!(ss,Session(s,session))
        end
    end
    tag != nothing && (ss = filter(x->tag in tags(x),ss))
    label != nothing && (ss = filter(x->label in labels(x),ss))
    return ss
end

function dsets(sess::Session,label::String;complete::Bool=true)
    ds = Dataset[]
    sd = sess_dict(sess)["labels"]
    haskey(sd,label) || return Dataset[]
    for x in sd[label]
        if (! complete) || x["complete"]
            push!(ds,Dataset(x["filename"],x["meta"],sess))
        end
    end
    return ds
end

function dsets(subj::Subject,label::String;complete::Bool=true,tag=nothing)
    ds = Dataset[]
    for sess in sessions(subj;tag)
        sd = sess_dict(sess)["labels"]
        haskey(sd,label) || return Dataset[]
        for x in sd[label]
            if (! complete) || x["complete"]
                push!(ds,Dataset(x["filename"],x["meta"],sess))
            end
        end
    end
    return ds
end

sess_dir(sess::Session) = "$(subj_dir(sess.subj.subjid))/sessions/$(sess.session)"
filename(dset::Dataset) = "$(sess_dir(dset.sess))/$(dset.filename)"

end # module


text:
# create documentation for julia
using Revise
using Documenter, COAST

makedocs(sitename="COAST.jl")

# to host live run:
# using LiveServer; servedocs()


text:
"""
    LD(A::Diagonal)
    LD(A::HBlikDiag)
    LD(A::DenseMatrix)

Return `log(det(tril(A)))` evaluated in place.
"""
LD(d::Diagonal{T}) where {T<:Number} = sum(log, d.diag)

function LD(d::UniformBlockDiagonal{T}) where {T}
    dat = d.data
    m, n, k = size(dat)
    m == n || throw(ArgumentError("Blocks of d must be square"))
    s = log(one(T))
    @inbounds for j in 1:k, i in 1:m
        s += log(dat[i,i,j])
    end
    s
end

function LD(d::DenseMatrix{T}) where {T}
    s = log(one(T))
    for i in 1:LinearAlgebra.checksquare(d)
        s += log(d[i, i])
    end
    s
end

"""
    logdet(m::LinearMixedModel)

Return the value of `log(det(Λ'Z'ZΛ + I)) + log(det(LX*LX'))` evaluated in place.

Here LX is the diagonal term corresponding to the fixed-effects in the blocked
lower Cholesky factor.
"""
function LinearAlgebra.logdet(m::LinearMixedModel{T}) where {T}
    s = log(one(T))
    Ldat = m.L.data
    @inbounds for (i, trm) in enumerate(m.trms)
        isa(trm, AbstractFactorReTerm) && (s += LD(Ldat[Block(i, i)]))
    end
    if m.optsum.REML
        feindex = length(m.trms) - 1
        fetrm = m.trms[feindex]
        if isa(fetrm, MatrixTerm)
            lblk = Ldat[Block(feindex, feindex)]
            for i in 1:fetrm.rank
                s += log(lblk[i, i])
            end
        end
    end
    2s
end

LinearAlgebra.logdet(m::GeneralizedLinearMixedModel) = logdet(m.LMM)


text:
#=
    Figure 6
=#

using Inference
using Model

using JLD2
using Plots
using StatsPlots
using Random

include("../figure_defaults.jl")

##############################################################
## Load model and results
##############################################################

include("../../Results/MainResult_Setup.jl")
@load "Results/MainResult.jld2"

include("../../Results/DeterministicModel.jl")

##############################################################
## Figure 5
##############################################################

# Parameter limits
lims = [[0.0,3.0],[0.0,0.5],[0.0,0.2]]

# Dimension combinations for bivariate plots
dims = [[1,2],[1,3],[2,3]]

plts1 = [plot_distribution(θ -> param_dist(θ[5:end-1]).v[i],sample(C,2000),θ;lims=lims[i],xlabel=["R","λ","β"][i]) for i = 1:3]
vline!(plts1[1],[1.0],lw=2.0,c=:red,ls=:dash)
vline!(plts1[2],[θdet.λ],lw=2.0,c=:red,ls=:dash)
vline!(plts1[3],[θdet.β],lw=2.0,c=:red,ls=:dash)

plts2 = [plot_biv_distribution(param_dist(θ[5:end-1]);
            dims=dims[i],
            lims=lims[dims[i]],
            labs=["R","λ","β"][dims[i]],
            c=:RdPu_9,levels=10,box=:on
        ) for i = 1:3]
scatter!(plts2[1],[1.0],[θdet.λ],c=:white,msw=0.0)
scatter!(plts2[2],[1.0],[θdet.β],c=:white,msw=0.0)
scatter!(plts2[3],[θdet.λ],[θdet.β],c=:white,msw=0.0)

fig5 = plot(plts1[1],plts2[1],plts1[2],plts2[2:3]...,plts1[3],layout=@layout([a _ _; b c _; d e f]))

plot!(fig5,
    widen=false,
    legend=:none,
    size=(600,600)
)
[plot!(fig5,yticks=[],axis=:x,subplot=i) for i = [1,3,6]]
plot!(fig5,yticks=[],subplot=5)

##############################################################
## Save...
##############################################################

savefig(fig5,"$(@__DIR__)/fig5.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
"""
    group_element_calculate(expr::String, group::Bijection{String, T}, identity_element::T) -> key::String

Calculate the expression `expr` using the labels defined in `group`. The group `identity_element` is necessary.

Usage:
1. Let `G` be a labelled group.
2. Let `I` be the identity element from `G`.
3. Define a shorthand:
    `Gcalc(expr::String) = group_element_calculate(expr, G, I)`
4. Use `Gcalc()``.
"""
function group_element_calculate(expr::String, group::Bijection{String, T}, identity_element::T)::String where T
    # if the identifier expr is already a key, return expr
    if expr in keys(group)
        return expr
    end
    # split expr, calculate element and return identifier from keys
    ks = split(expr, "")
    el = identity_element
    for k in ks
        el *= group[k]
    end
    return group(el)
end


text:
using LinearAlgebra
using OrdinaryDiffEq
using ControlSystems
using PyPlot
using GaussianProcesses, BayesianOptimization, Distributions

const Kt = 1/(60.3e-3)
const r = 45/34*13/3        # Gear and belt ratio
const mp = 0.706
# const mr = 0.56424  # m_resin + 4*m_steel = 0.084 + 4 * 0.12946
const mr = 0.223
const m = mp + mr
const lr = 0.18217
const lp = lr/2
const l = (mp*lp + mr*lr)/m
const g = 9.81
const Ip_com = 0.00249
const I_epoxy = 1100 * 7.919e4 / 1e9 * (0.18/2)^2
const I_steel = 9.484e-4
const I_resin = 7.559e-4
const I_resin_new = 4.64e-4
# const I1 = Ip_com + mp*lp^2 + mr*lr^2
# const I2 = I_resin_new + 4*I_steel
const I2 = I_resin + I_epoxy
const I1_estimate = Ip_com + mp*lp^2 + mr*lr^2
const I1_true = I1_estimate*(1 + 0.22);
const τ_max = 0.2
const γ1 = 2. # 1.5*ceil(m*g*l)
const γ2 = 25. # 3*ceil(I1/I2) * γ1 / (γ1 - m*g*l)
const k = (1.0, 10.0)
I1 = I1_true

function compute_cost(sol::ODESolution, par::AbstractArray)
    C = [1 0 0 0; 0 1 0 0];
    Q = C' * C
    R = 100.0
    ctrl(u) = -dot(par,u)
    return sum( 1/2*dot(sol.u[i], Q*sol.u[i]) + 1/2*R*ctrl(sol.u[i])^2 for i=1:length(sol) )
end

function compute_lqr_gain(I1)
    A = [0 0 1 0; 0 0 0 1; m*g*l/I1 0 0 0; 0 0 0 0]
    B = r*[0; 0; -1/I1; 1/I2]
    C = [1 0 0 0; 0 1 0 0];
    sys = ss(A, B, C, 0)
    Q = C' * C
    # Q = 1.0*I
    R = 100.0*I
    L = lqr(sys, Q, R)
    return L
end

compute_lqr_estimate(u) = -(K_estimate * u)[1]
compute_lqr_true(u) = -(K_true * u)[1]

function compute_ida_pbc(u)
    q1, q2, q1dot, q2dot = u
    qbar = q2 + γ2*q1
    qbar_dot = q2dot + γ2*q1dot
    u_es = γ1*sin(q1) + k[1]*sin(qbar)
    u_di = k[2]*qbar_dot
    return 1/r*( u_es + u_di )
end

function compute_switching_control(u)
    q1, q2, q1dot, q2dot = u
    τ = compute_ida_pbc(u)
    if 1 - cos(q1) + 1 - cos(q2) + abs(q1) + abs(q2) < 0.1
        τ = compute_lqr_estimate(u)
        # τ = compute_lqr_true(u)
    end
    return τ
end

function iwp_eom!(du, u, p, t)
    # τ = clamp(compute_lqr_estimate(u), -τ_max, τ_max)
    # τ = clamp(compute_lqr_true(u), -τ_max, τ_max)
    # τ = clamp(compute_ida_pbc(u), -τ_max, τ_max)
    # τ = clamp(compute_switching_control(u), -τ_max, τ_max)

    τ = clamp(-dot(p,u), -τ_max, τ_max)

    du[1] = u[3]
    du[2] = u[4]
    du[3] = 1/I1 * (m*g*l*sin(u[1]) - r*τ)
    du[4] = 1/I2 * r*τ
end

function compute_control(sol::ODESolution)
    # return hcat( [clamp((-K*sol.u[i])[1], -τ_max, τ_max) for i = 1:length(sol)]... )'
    # return hcat( [clamp(compute_ida_pbc(sol.u[i]), -τ_max, τ_max) for i = 1:length(sol)]... )'
    # return hcat( [clamp(compute_lqr_estimate(sol.u[i]), -τ_max, τ_max) for i = 1:length(sol)]... )'
    return hcat( [clamp(compute_lqr_true(sol.u[i]), -τ_max, τ_max) for i = 1:length(sol)]... )'
end

K_estimate = compute_lqr_gain(I1_estimate)
K_true = compute_lqr_gain(I1_true)
u0 = [1*π/180, 0.5, 0, 0]
tspan = (0.0, 2.0)
p = K_true
prob = ODEProblem(iwp_eom!, u0, tspan, p)
sol = solve(prob, Tsit5(), saveat=0.01, dtmax=0.01)

function f(x)
    prob = ODEProblem(iwp_eom!, u0, tspan, x)
    sol = solve(prob, Tsit5(), saveat=0.01, dtmax=0.01)
    return compute_cost(sol, x)
end

# Choose as a model an elastic GP with input dimensions 2.
# The GP is called elastic, because data can be appended efficiently.
model = ElasticGPE(4,                            # 2 input dimensions
                   mean = MeanConst(0.),
                   kernel = SEArd([0., 0., 0., 0.], 5.),
                   logNoise = 0.,
                   capacity = 3000)              # the initial capacity of the GP is 3000 samples.
set_priors!(model.mean, [Normal(1, 2)])

# Optimize the hyperparameters of the GP using maximum a posteriori (MAP) estimates every 50 steps
modeloptimizer = MAPGPOptimizer(every = 50, noisebounds = [-4, 3],       # bounds of the logNoise
                                kernbounds = [[-1, -1, -1, -1, 0], [4, 4, 4, 4, 10]],  # bounds of the 3 parameters GaussianProcesses.get_param_names(model.kernel)
                                maxeval = 40)

# modeloptimizer = NoModelOptimizer()

opt = BOpt(f,
            model,
            UpperConfidenceBound(),                   # type of acquisition
            modeloptimizer,
            [-10., -10., -10., -10.], [0., 0., 0., 0.],                     # lowerbounds, upperbounds
            repetitions = 5,                          # evaluate the function for each input 5 times
            maxiterations = 1000,                      # evaluate at 100 input positions
            sense = Min,                              # minimize the function
            acquisitionoptions = (method = :LD_LBFGS, # run optimization of acquisition function with NLopts :LD_LBFGS method
                                restarts = 5,       # run the NLopt method from 5 random initial conditions each time.
                                maxtime = 0.1,      # run the NLopt method for at most 0.1 second each time
                                maxeval = 1000),    # run the NLopt methods for at most 1000 iterations (for other options see https://github.com/JuliaOpt/NLopt.jl)
            verbosity = Progress)

# x = Array{Array{Float64, 1}, 1}()
# push!(x, K_estimate[:])
# push!(x, K_true[:])
# y = -f.(x)
# append!(model, hcat(x...), y)

opt = BOpt(f,
           model,
        #    UpperConfidenceBound(),
           ThompsonSamplingSimple(),
           modeloptimizer,
           [-10., -10., -10., -10.], [0., 0., 0., 0.],
           maxiterations = 1000,
           sense = Min,
           initializer_iterations = 0
          )

result = boptimize!(opt)

################################################################################
## Plotting

# fig = figure(1)
# fig.clf()
# ax = fig.add_subplot(2,1,1)
# ax.plot(sol.t, getindex.(sol.u, 1), label=L"$q_1$")
# ax.plot(sol.t, getindex.(sol.u, 3), label=L"$\dot{q}_1$")
# ax.legend(fontsize=15)
# ax = fig.add_subplot(2,1,2)
# ax.plot(sol.t, getindex.(sol.u, 2), label=L"$q_2$")
# ax.plot(sol.t, getindex.(sol.u, 4), label=L"$\dot{q}_2$")
# ax.legend(fontsize=15)
# ax.set_xlabel("time", fontsize=15)

# fig = figure(2)
# fig.clf()
# ax = fig.add_subplot(1,1,1)
# ax.plot(sol.t, compute_control(sol), label=L"$\tau$")
# ax.set_xlabel("time", fontsize=15)
# ax.legend(fontsize=15)


text:
#########
# This file is to add support for lowering broadcast notation
#       y .= α .* Mul(A,x) .+ β .* y
# to
#       materialize!(MulAdd(α, A, x, β, y))
# which then becomes a blas call.
#########

struct ArrayMulArrayStyle{StyleA, StyleB, p, q} <: BroadcastStyle end

@inline copyto!(dest::AbstractArray, bc::Broadcasted{<:ArrayMulArrayStyle}) =
    _copyto!(MemoryLayout(dest), dest, bc)
# Use default broacasting in general
@inline _copyto!(_, dest, bc::Broadcasted) = copyto!(dest, Broadcasted{Nothing}(bc.f, bc.args, bc.axes))

const BArrayMulArray{styleA, styleB, p, q, T, V} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q}, <:Any, typeof(identity),
                <:Tuple{<:ArrayMulArray{styleA,styleB,p,q,T,V}}}
const BConstArrayMulArray{styleA, styleB, p, q, T, U, V} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q},
                    <:Any, typeof(*),
                    <:Tuple{T,<:ArrayMulArray{styleA,styleB,p,q,U,V}}}
const BArrayMulArrayPlusArray{styleA, styleB, p, q, T, U, V} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q},
                <:Any, typeof(+),
                <:Tuple{<:ArrayMulArray{styleA,styleB,p,q,T,U},<:AbstractArray{V,q}}}
const BArrayMulArrayPlusConstArray{styleA, styleB, p, q, T, U, V, W} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q},
                <:Any, typeof(+),
                <:Tuple{<:ArrayMulArray{styleA,styleB,p,q,T,U},
                Broadcasted{DefaultArrayStyle{q},<:Any,typeof(*),
                            <:Tuple{V,<:AbstractArray{W,q}}}}}
const BConstArrayMulArrayPlusArray{styleA, styleB, p, q, T, U, V, W} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q},
                <:Any, typeof(+),
                <:Tuple{Broadcasted{<:ArrayMulArrayStyle{styleA,styleB,p,q},
                                    <:Any, typeof(*),
                                    <:Tuple{T,<:ArrayMulArray{styleA,styleB,p,q,U,V}}},
                        <:AbstractArray{W,q}}}
const BConstArrayMulArrayPlusConstArray{styleA, styleB, p, q, T, U, V, W, X} =
    Broadcasted{ArrayMulArrayStyle{styleA,styleB,p,q},
                <:Any, typeof(+),
                <:Tuple{Broadcasted{<:ArrayMulArrayStyle{styleA,styleB,p,q},
                                    <:Any, typeof(*),
                                    <:Tuple{T,<:ArrayMulArray{styleA,styleB,p,q,U,V}}},
                        Broadcasted{DefaultArrayStyle{q},<:Any,typeof(*),<:Tuple{W,<:AbstractArray{X,q}}}}}

BroadcastStyle(::Type{<:ArrayMulArray{StyleA,StyleB,p,q}}) where {StyleA,StyleB,p,q} =
    ArrayMulArrayStyle{StyleA,StyleB,p,q}()
BroadcastStyle(M::ArrayMulArrayStyle, ::DefaultArrayStyle) = M
BroadcastStyle(::DefaultArrayStyle, M::ArrayMulArrayStyle) = M
similar(M::Broadcasted{<:ArrayMulArrayStyle}, ::Type{ElType}) where ElType =
    Array{ElType}(undef,size(M.args[1]))

broadcastable(M::ArrayMulArray) = M
instantiate(bc::Broadcasted{<:ArrayMulArrayStyle}) = bc

@inline function _copyto!(_, dest::AbstractArray{T}, M::ArrayMulArray) where T
    A,B = M.factors
    materialize!(MulAdd(one(T), A, B, zero(T), dest))
end

# Use copyto! for y .= Mul(A,b)
@inline function _copyto!(_, dest::AbstractArray, bc::BArrayMulArray)
    (M,) = bc.args
    copyto!(dest, M)
end

@inline function _copyto!(_, dest::AbstractArray{T}, bc::BConstArrayMulArray) where T
    α,M = bc.args
    A,B = M.factors
    materialize!(MulAdd(α, A, B, zero(T), dest))
end

@inline function _copyto!(_, dest::AbstractArray{T}, bc::BArrayMulArrayPlusArray) where T
    M,C = bc.args
    A,B = M.factors
    copyto!(dest, MulAdd(one(T), A, B, one(T), C))
end

@inline function _copyto!(_, dest::AbstractArray{T}, bc::BArrayMulArrayPlusConstArray) where T
    M,βc = bc.args
    β,C = βc.args
    A,B = M.factors
    copyto!(dest, MulAdd(one(T), A, B, β, C))
end

@inline function _copyto!(_, dest::AbstractArray{T}, bc::BConstArrayMulArrayPlusArray) where T
    αM,C = bc.args
    α,M = αM.args
    A,B = M.factors
    copyto!(dest, MulAdd(α, A, B, one(T), C))
end

@inline function _copyto!(_, dest::AbstractArray, bc::BConstArrayMulArrayPlusConstArray)
    αM,βc = bc.args
    α,M = αM.args
    A,B = M.factors
    β,C = βc.args
    copyto!(dest, MulAdd(α, A, B, β, C))
end


text:
using Markdown
using Kyulacs
using Kyulacs:LazyHelp
using Documenter
import Documenter.Utilities.MDFlatten: mdflatten

#DocMeta.setdocmeta!(Kyulacs, :DocTestSetup, :(using Kyulacs); recursive=true)

mdflatten(io::IOBuffer, h::LazyHelp, md::Markdown.MD) = nothing

makedocs(;
    modules=[Kyulacs],
    authors="Satoshi Terasaki <terasakisatoshi.math@gmail.com> and contributors",
    repo="https://github.com/AtelierArith/Kyulacs.jl/blob/{commit}{path}#{line}",
    sitename="Kyulacs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://AtelierArith.github.io/Kyulacs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
    doctest=false,
)

deploydocs(;
    repo="github.com/AtelierArith/Kyulacs.jl",
    devbranch="main",
)


text:
#!/usr/bin/julia
push!(LOAD_PATH, dirname(@__FILE__))
push!(LOAD_PATH, "/usr/share/julia/DistroHelper")

using DistroHelper
using Pkg

function main(argv)
	if length(argv) == 0
		exit(0)
	elseif argv[1] in ("slug", "depends")  # raw print
		print(getfield(DistroHelper, Symbol(argv[1]))(argv[2:end]...))
	elseif argv[1] in ("gather_stdlib_manifest", "distro_manifest_")  # requires pretty print
		result = getfield(DistroHelper, Symbol(argv[1]))(argv[2:end]...)
		Pkg.TOML.print(result)
	elseif argv[1] in ("distro_project", "distro_manifest", "distro_merge", "metadata", "download")  # no print
		getfield(DistroHelper, Symbol(argv[1]))(argv[2:end]...)
	else
		throw(Base.error("unsupported args"))
	end
end

main(ARGS)