text:
n = 3
m = 2
T = 5

nz = n * (T - 1) + m * (T - 1)
nd = n * (T - 1)

u_idx = [collect((t - 1) * (m + n) .+ (1:m)) for t = 1:T-1]
x_idx = [collect((t - 1) * (m + n) + m .+ (1:n)) for t = 1:T-1]
n_idx = [(t - 1) * n .+ (1:n) for t = 1:T-1]

A = [rand(n,n) for t = 1:T-1]
B = [rand(n,m) for t = 1:T-1]
Q = [Diagonal(rand(n)) for t = 1:T]
R = [Diagonal(rand(m)) for t = 1:T-1]

x_ref = [randn(n) for t = 1:T-1]
u_ref = [randn(m) for t = 1:T-1]

x_init = ones(n)
z0 = rand(nz + nd)

function constraints(z)
	c = zeros(eltype(z), nd)

	for t = 1:T-1
		x2 = z[x_idx[t]]
		x1 = (t == 1 ? x_init : z[x_idx[t-1]])
		u1 = z[u_idx[t]]
		c[n_idx[t]] = x2 - A[t] * x1 - B[t] * u1
	end

	return c
end

constraints(z0)

function constraints_jacobian(z)
	C = zeros(eltype(z), nd, nz)

	for t = 1:T-1
		C[n_idx[t], u_idx[t]] = -B[t]
		C[n_idx[t], x_idx[t]] = Diagonal(ones(n))
		t == 1 && continue
		C[n_idx[t], x_idx[t-1]] = -A[t]
	end

	return C
end

constraints_jacobian(z0)

function objective_gradient(z)
	j = zeros(eltype(z), nz)
	for t = 1:T-1
		u1 = z[u_idx[t]]
		x1 = z[x_idx[t]]
		j[u_idx[t]] = R[t] * (u1 - u_ref[t])
		j[x_idx[t]] = Q[t] * (x1 - x_ref[t])
	end
	return j
end

objective_gradient(z0)

function objective_hessian(z)
	J = zeros(eltype(z), nz, nz)
	for t = 1:T-1
		J[u_idx[t], u_idx[t]] = R[t]
		J[x_idx[t], x_idx[t]] = Q[t]
	end
	return J
end

objective_hessian(z0)

function hessian(z)
	J = objective_hessian(z)
	C = constraints_jacobian(z)

	return [J C'; C zeros(eltype(z), nd, nd)]
end

function gradient(z)
	y = z[nz .+ (1:nd)]
	return [objective_gradient(z) + constraints_jacobian(z)' * y; constraints(z)]
end

u = [randn(m) for t = 1:T-1]
x = [x_init]
for t = 1:T-1
	push!(x, A[t] * x[end], + B[t] * u[t])
end
z = [vcat([[u[t]; x[t]] for t = 1:T-1]...); zeros(nd)]

res_norm = norm(gradient(z), 1)

Δ = hessian(z) \ gradient(z)

α = 1.0
z_cand = z - α * Δ

res_cand_norm = norm(gradient(z_cand), 1)

z_cand[nz .+ (1:nd)]


text:
@everywhere function write(i)
  Libc.systemsleep(5)
  (myid(),i)
end

sss = Future[]
println("Write")

for i=1:nprocs() println("$i"); push!(sss, @spawn write(i)) end
println("Read - Wait ~ 5 sec")
for s in sss println(fetch(s)) end

#s = @spawnat 1 myid()
#println(fetch(s))

ext: jl

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function laplace_iter_devec(u, dx2, dy2, Niter, N)
    uout = copy(u)
    for iter = 1:Niter
        for i = 2:N-1
            for j = 2:N-1
                uout[i,j] = ( (u[i-1,j]+u[i+1,j])*dy2 + (u[i,j-1]+u[i,j+1])*dx2 ) * (1./(2*(dx2+dy2)))
            end
        end
        u, uout = uout, u
    end
    return u
end

function laplace_devec()
    N = 150
    u = zeros(N, N)
    u[1,:] = 1
    Niter = 2^10
    dx2 = dy2 = 0.1*0.1
    u = laplace_iter_devec(u, dx2, dy2, Niter, N)
end

function laplace_iter_vec(u, dx2, dy2, Niter, N)
    for i = 1:Niter
        u[2:N-1, 2:N-1] = ((u[1:N-2, 2:N-1] + u[3:N, 2:N-1])*dy2 + (u[2:N-1,1:N-2] + u[2:N-1, 3:N])*dx2) * (1./ (2*(dx2+dy2)))
    end
    return u
end

function laplace_vec()
    N = 150
    u = zeros(N,N)
    u[1,:] = 1
    Niter = 2^10
    dx2 = dy2 = 0.1*0.1
    u = laplace_iter_vec(u, dx2, dy2, Niter, N)
end

function laplace_iter_vec_sub(u, dx2, dy2, Niter, N)
    for i = 1:Niter
        u[2:N-1, 2:N-1] = ((view(u, 1:N-2, 2:N-1) + view(u,3:N, 2:N-1))*dy2 + (view(u,2:N-1,1:N-2) + view(u,2:N-1, 3:N))*dx2) * (1./ (2*(dx2+dy2)))
    end
    return u
end

function laplace_vec_sub()
    N = 150
    u = zeros(N,N)
    u[1,:] = 1
    Niter = 2^10
    dx2 = dy2 = 0.1*0.1
    u = laplace_iter_vec_sub(u, dx2, dy2, Niter, N)
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
using ConjugatePriors

data = rand(Bernoulli(), 100)

prior = Beta(1, 1)
posterior = update(prior, data)

prior = Beta(1000, 100)
posterior = update(prior, data)


text:
export Dirac

struct Dirac{X} <: AbstractMeasure
    x::X
end

function (μ::Dirac{X})(s) where {X}
    μ.x ∈ s && return 1
    return 0
end

density_def(μ::Dirac, x) = x == μ.x

logdensity_def(μ::Dirac, x) = (x == μ.x) ? 0.0 : -Inf

basemeasure(::Dirac) = CountingMeasure()

basemeasure_type(::Dirac) = CountingMeasure

Base.rand(::Random.AbstractRNG, T::Type, μ::Dirac) = μ.x

export dirac

dirac(d::AbstractMeasure) = Dirac(rand(d))


text:
# Convert gff3 to BED12 format
# Usage: julia gff3ToBed12.jl path_gff3 path_bed

# Assumption:
# - The input gff3 file has to contain rows of "exon", "start_codon", "stop_codon".
# - chromStrat and chromEnd is set to the left-most exon's start (0-based) and the right-most exon's end positions, respectively.
# - thickStart and thickEnd is set to the minimum (0-based) and maximum (1-based) coordinates of start_codon and stop_codon.
# - If a transcript_id doesn't have start_codon or stop codon, thickStart and thickEnd is set to the value of chromEnd.

# Test:
	# $ grep ENSMUST00000070533 Mouse_UCSC_mm10.Ensembl.bed
# 	# chr1    3214481 3671498 ENSMUST00000070533      0       -       3216021 3671348 0       3       2487,200,947,   0,207220,456070,
	# $ grep ENSMUST00000070533 gencode.vM9.annotation.bed
	# chr1    3214481 3671498 ENSMUST00000070533.4    0       -       3216021 3671348 0       3       2487,200,947,   0,207220,456070,

function gff3ToBed12(path_gff3, path_bed)
	transcript_id = ""
	chr = ""
	set_start = BitSet()
	set_end = BitSet()
	set_codon = BitSet()
	strand = ""

	fw = open(path_bed, "w")

	f = open(path_gff3)
	for ln in eachline(f)
		if occursin(r"^#", ln)
			continue
		end

		arr = split(chomp(ln), "\t")
		if arr[3] == "gene"
			continue
		end

		new_transcript_id = match(r"transcript_id=(.+?);", arr[9]).captures[1]
		if transcript_id != new_transcript_id
			if length(set_start) > 0
				res = calcLine(chr, set_start, set_end, transcript_id, strand, set_codon)
				write(fw, res * "\n")
			end

			transcript_id = new_transcript_id
			chr = arr[1]
			empty!(set_start)
			empty!(set_end)
			empty!(set_codon)
			strand = arr[7]
		end

		if arr[3] == "exon"
			push!(set_start, parse(Int64, arr[4]))
			push!(set_end, parse(Int64, arr[5]))
		elseif arr[3] == "stop_codon" || arr[3] == "start_codon"
			push!(set_codon, parse(Int64, arr[4]))
			push!(set_codon, parse(Int64, arr[5]))
		end

	end
	close(f)

	if length(set_start) > 0
		res = calcLine(chr, set_start, set_end, transcript_id, strand, set_codon)
		write(fw, res * "\n")
	end
	close(fw)
end

function calcBlockSizes(set_start::BitSet, set_end::BitSet)
	join([string(collect(set_end)[i] - collect(set_start)[i] + 1) * "," for i in range(1, stop=length(set_start))], "")
end

function calcBlockStarts(set_start::BitSet, chromStart)
	join([string(i - chromStart) * "," for i in set_start ], "")
end

function calcThickStart(set_codon::BitSet, chromEnd)
	if length(set_codon) == 0
		chromEnd
	else
		first(set_codon) - 1
	end
end

function calcThickEnd(set_codon::BitSet, chromEnd)
	if length(set_codon) == 0
		chromEnd
	else
		last(set_codon)
	end
end

function calcLine(
					chr::AbstractString, set_start::BitSet, set_end::BitSet,
					transcript_id::AbstractString, strand::AbstractString,
					set_codon::BitSet
					)
	chrom = chr
	if length(set_start) == 0
		print(transcript_id)
	end

	chromStart = first(set_start) - 1
	chromEnd = last(set_end)
	name = transcript_id
	score = 0
	strand = strand
	thickStart = calcThickStart(set_codon, chromEnd)
	thickEnd = calcThickEnd(set_codon, chromEnd)
	itemRgb = 0
	blockCount = length(set_start)
	blockSizes = calcBlockSizes(set_start, set_end)
	blockStarts = calcBlockStarts(set_start, chromStart + 1)

	res = join([
				chrom;
				string(chromStart);
				string(chromEnd);
				name;
				string(score);
				strand;
				string(thickStart);
				string(thickEnd);
				string(itemRgb);
				blockCount;
				blockSizes;
				blockStarts;
			], "\t")
	return res
end

# chrom - The name of the chromosome (e.g. chr3, chrY, chr2_random) or scaffold (e.g. scaffold10671).
# chromStart - The starting position of the feature in the chromosome or scaffold. The first base in a chromosome is numbered 0.
# chromEnd - The ending position of the feature in the chromosome or scaffold. The chromEnd base is not included in the display of the feature. For example, the first 100 bases of a chromosome are defined as chromStart=0, chromEnd=100, and span the bases numbered 0-99.
# name - Defines the name of the BED line. This label is displayed to the left of the BED line in the Genome Browser window when the track is open to full display mode or directly to the left of the item in pack mode.
# score - A score between 0 and 1000. If the track line useScore attribute is set to 1 for this annotation data set, the score value will determine the level of gray in which this feature is displayed (higher numbers = darker gray). This table shows the Genome Browser's translation of BED score values into shades of gray: shade score in range  	≤ 166	167-277	278-388	389-499	500-611	612-722	723-833	834-944	≥ 945
# strand - Defines the strand - either '+' or '-'.
# thickStart - The starting position at which the feature is drawn thickly (for example, the start codon in gene displays). When there is no thick part, thickStart and thickEnd are usually set to the chromStart position.
# thickEnd - The ending position at which the feature is drawn thickly (for example, the stop codon in gene displays).
# itemRgb - An RGB value of the form R,G,B (e.g. 255,0,0). If the track line itemRgb attribute is set to "On", this RBG value will determine the display color of the data contained in this BED line. NOTE: It is recommended that a simple color scheme (eight colors or less) be used with this attribute to avoid overwhelming the color resources of the Genome Browser and your Internet browser.
# blockCount - The number of blocks (exons) in the BED line.
# blockSizes - A comma-separated list of the block sizes. The number of items in this list should correspond to blockCount.
# blockStarts - A comma-separated list of block starts. All of the blockStart positions should be calculated relative to chromStart. The number of items in this list should correspond to blockCount.

path_gff3 = ARGS[1]
path_bed = ARGS[2]

gff3ToBed12(path_gff3, path_bed)

max_stars_repo_name: rikenbit/RamDAQ-prototype

max_issues_repo_name: rikenbit/RamDAQ-prototype

max_forks_repo_name: rikenbit/RamDAQ-prototype


text:
## geometry functions
"""
    pln=three_points_to_plane(A)

Compute equation for a plane from the three points defined in A. The plane is parameterized
as `a*x+b*y+c*z+d==0` with the coefficients returned as `pln=[a,b,c,d]`.

# Arguments
- A::3×3-Array : Array containing the three points defining the plane as columns.

# Notes
Code adapted from https://www.geeksforgeeks.org/program-to-find-equation-of-a-plane-passing-through-3-points/
"""
function three_points_to_plane(A)
    A,B,C=A[:,1],A[:,2],A[:,3]
    u=A-C
    v=B-C
    a=LinearAlgebra.cross(u,v)
    a/=LinearAlgebra.norm(a) #normalize
    d=-LinearAlgebra.dot(a,C)
    if d<1E-7 #TODO: find a better fix for the intersection problem
        d=0.0
    end
    a,b,c=a
    return [a,b,c,d]
end

"""
    check = is_point_in_plane(pnt, pln)

Check whether point "pnt" is in plane . Note that there is currently no handling of round off errors.
"""
function is_point_in_plane(pnt,pln)
    #TODO: implement round off handling
    return pln[1]*pnt[1]+pln[2]*pnt[2]+pln[3]*pnt[3]+pln[4]
end
"""
    p = reflect_point_at_plane(pnt,pln)

Reflect point `pnt` at plane `pln` and return as `p`.

# Notes
Code adapted from https://www.geeksforgeeks.org/mirror-of-a-point-through-a-3-d-plane/
"""
function reflect_point_at_plane(pnt,pln)
    a,d=pln[1:3],pln[4]
    k=(-LinearAlgebra.dot(a,pnt)-d)#/(LinearAlgebra.dot(a,a))
    p=a.*k.+pnt #foot of the perpendicular
    p=2*p-pnt #reflected point
    return p
end
"""
    foot = find_foot_of_perpendicular(pnt,pln)

Compute the position of the foot of the perpendicular from the point `pnt` to the plane `pln`.

# Notes
Code adapted from https://www.geeksforgeeks.org/mirror-of-a-point-through-a-3-d-plane/
"""
function find_foot_of_perpendicular(pnt,pln)
    a,d=pln[1:3],pln[4]
    k=(-LinearAlgebra.dot(a,pnt)-d)#/(LinearAlgebra.dot(a,a))
    foot=a.*k.+pnt #foot of the perpendicular
    return foot
end

"""
 make_normal_outwards!(pln,testpoint)

Ensure that the parameterization of the plane "pln" has a normal that is directed twoards `testpoint`.
If necessary reparametrize `pln`.
"""
function make_normal_outwards(pln,testpoint)
    foot=find_foot_of_perpendicular(testpoint,pln)
    scalar_product=LinearAlgebra.dot(pln[1:3],testpoint-foot)
    pln*=-sign(scalar_product)
    return pln
end

"""
    p,n = find_intersection_of_two_planes(pln1,pln2)

Find a parametrization of the axis defined by the intersection of the planes `pln1`and `pln2`.
The axis is parameterized by a point `p` and direction `n`.

# Notes
The algorithm follows John Krumm's solution for finding a common point on two intersecting planes as it is
explained in https://math.stackexchange.com/questions/475953/how-to-calculate-the-intersection-of-two-planes
To simplify the code, here the point p0 is the origin.
"""
function find_intersection_of_two_planes(pln1,pln2)
    n=LinearAlgebra.cross(pln1[1:3],pln2[1:3]) #axis vector
    n/=LinearAlgebra.norm(n)
    #TODO: check whether vector is nonzero and planes truly intersect!
    #generate points on the planes closest to origin.
    p1=find_foot_of_perpendicular([0.,0.,0.],pln1)
    p2=find_foot_of_perpendicular([0.,0.,0.],pln2)
    #Linear system from Lagrange optimization
    M=[2.0 0.0 0.0 pln1[1] pln2[1];
       2.0 0.0 0.0 pln1[2] pln2[2];
       2.0 0.0 0.0 pln1[3] pln2[3];
     pln1[1] pln1[2] pln1[3] 0 0;
     pln2[1] pln2[2] pln2[3] 0 0]

     rhs=[0;0;0; LinearAlgebra.dot(p1,pln1[1:3]); LinearAlgebra.dot(p2,pln2[1:3])]
     p=rhs\M
     return p[1:3],n
end

"""
    R = create_rotation_matrix_around_axis(n,α)

Compute the rotation matrix around the directional vector `n`rotating by an angle `α` (in rad).
"""
function create_rotation_matrix_around_axis(n,α)
    R=[n[1]^2*(1-cos(α))+cos(α) n[1]*n[2]*(1-cos(α))-n[3]*sin(α) n[1]*n[3]*(1-cos(α))+n[2]*sin(α);
       n[2]*n[1]*(1-cos(α))+n[3]*sin(α) n[2]^2*(1-cos(α))+cos(α) n[2]*n[3]*(1-cos(α))-n[1]*sin(α);
       n[3]*n[1]*(1-cos(α))-n[2]*sin(α) n[3]*n[2]*(1-cos(α))+n[1]*sin(α) n[3]^2*(1-cos(α))+cos(α)]
       return R
end

## mesh book keeping
function find_testpoint_idx(tri, tet)
    testpoint_idx=0
    for idx in tet
        if !(idx in tri)
            testpoint_idx=idx
        end
    end
    return testpoint_idx
end

"""
    ridx = get_rotated_index(idx,sector, naxis, nxsector, DOS)

Get index `ridx` of point in sector `sector` created from rotation of the point
with index `idx`. The definig mesh feature `naxis` points on the center axis,
`nxsector` points in a unit cell excluding the center axis, and has a degree of
symmetry `DOS`.
"""
function get_rotated_index(idx,sector, naxis, nxsector, DOS)
    if idx <= naxis
        idx=idx
    else
        idx=idx+nxsector*sector #sector counting starts from zero
    end
    idx=((idx-naxis-1)%(nxsector*DOS))+naxis+1
    if idx<0
        idx+=DOS*nxsector
    end

    return idx
end

"""
    ridx = get_reflected_index(idx,naxis,nxbloch,nbody,shiftbody,nxsymmetry)

Get index `ridx` of a point in the first unit cell created from reflection of
the point with index `idx`. The definig mesh feature `naxis` points on the
center axis, `nxbloch` points on the bloch plane (excluding the center axis),
`nbody` points in the interiror of the first half-cell, `nxsymmetry`points on
the symmetry plane (excluding the center axis) `nxsector` points in a unit cell
excluding the center axis. The difference of the indeces of a refernce body point
and a reflected body point is `shiftbody`. The number of points per unit cell
(excluding the center axis) is `nxsector`.
"""
function get_reflected_index(idx,naxis,nxbloch,nbody,shiftbody,nxsymmetry,nxsector)
    if idx <= naxis
        idx=idx
    elseif idx <= naxis+nxbloch
        idx=idx+nxsector
    elseif idx <= naxis+nxbloch+nbody
            idx=idx+shiftbody
    elseif idx <=  naxis+nxbloch+nbody+nxsymmetry
        idx=idx
    else
        idx=0 #error flag
    end
    return idx
end

"""
    sector = get_point_sector(pnt_idx,naxis, nxsector)

Get the sector containing point with index `pnt_idx`, in a mesh that features
`naxis` grid points on the center line and `nxsector` grid points in a unit cell
excluding the center axis.
"""
function get_point_sector(pnt_idx,naxis, nxsector)
    if pnt_idx<=naxis
        return typemax(0)
    else
        return (pnt_idx-naxis-1)÷nxsector #-1 is necessary because otherwise last point would count to next sector
    end
end

"""
    s = get_line_sector(ln)

Get index `s` of sector containing line `ln`, in a mesh that features
`naxis` grid points on the center line and `nxsector` grid points in a unit cell
excluding the center axis..
"""
function get_line_sector(ln,naxis,nxsector)
    s=typemax(0)
    for (idx,p) in enumerate(ln)
        new_s=get_point_sector(p,naxis,nxsector)
        if new_s==s #sector has not changed do nothing
        elseif s==typemax(0) #some sector is identified that is not axis
            s=new_s
        elseif abs(new_s-s)==1 #standard sector boundary
            s=min(s,new_s)
        elseif abs(new_s-s)>1# #periodic sector boundary #TODO: Does not work if DOS==2
            if new_s!=typemax(0)
                s=max(s,new_s)
            end
        else
            println("Error: Could not identify sector for line.")
            return nothing
        end
    end
    return s
end

"""
    ln_idx = get_line_idx(ln,naxis_ln,nxsector_ln)

Get index `ln_idx` of line `ln`  in a mesh that features `naxis` grid points and
`naxis_ln' lines on the center line, `nxsector` grid points and `nxsector_ln'
lines in a unit cell excluding the center axis.
"""
function get_line_idx(lines,ln,naxis,nxsector, naxis_ln,nxsector_ln,DOS)
    s=get_line_sector(ln,naxis, nxsector)
    if s==typemax(0)
        return find_smplx(lines[1:naxis_ln],ln)
    else
        return find_smplx(lines[naxis_ln+1:naxis_ln+nxsector_ln],get_rotated_index.(ln,-s,naxis,nxsector, DOS))+s*nxsector_ln+naxis_ln
    end
end

##
"""
    full_mesh=extend_mesh(mesh::Mesh, doms; sym_name="Symmetry", blch_name="Bloch", unit=false)

Create full mesh or unit cell from half cell represented in `mesh`.

# Arguments
- `mesh::Mesh`: mesh representing the half cell. The mesh must span a sector of `2π/2N`, where `N` is the (integer) degree of symmetry of the full mesh.
- `doms::List`: list of 2-tuples containing the domain names and their `copy_degree` (see notes below).
- `sym_name::String="Symmetry"`: name of the domain of `mesh` that forms the symmetry plane of the half-mesh.
- `blch_name::String="Bloch"`: name of the domain of `mesh` that forms the remaining azimuthal plane of the half-mesh.
- `unit::Bool=false`: toggle whether extend the mesh to unit cell only.

# Returns
- `full_mesh::Mesh`: representation of the full mesh.

# Notes
The routine copies only domains that are specified  in `doms`. These domains are extended according to the specified `copy_degree`.
The following are available:

- `:full`: extent the domain and save it under the same name in `full_mesh`.
- `:unit`: extent the domain and save the individual unit cells labeled from `0` to `N-1` in `full_mesh`.
- `:half`: extent the domain and save the individual unit cells as half cells labeled from `0` to `N-1` in `full_mesh` where one half-cell  contains `_img` in its name.

Multiple styles can be mixed in one domain specification. An example for `doms` would be
    doms=[("Interior", :full), ("Outlet", :unit), ("Flame", :half)]
"""
function extend_mesh(mesh::Mesh, doms; sym_name="Symmetry", blch_name="Bloch", unit=false)
    collect_lines!(mesh) #TODO: make it independent of line collection
    npoints=size(mesh.points,2)

    bloch=[]
    for smplx_idx in mesh.domains[blch_name]["simplices"]
        for pnt in mesh.triangles[smplx_idx]
            insert_smplx!(bloch,pnt)
        end
    end

    symmetry=[]
    for smplx_idx in mesh.domains[sym_name]["simplices"]
        for pnt in mesh.triangles[smplx_idx]
            insert_smplx!(symmetry,pnt)
        end
    end

    axis=intersect(bloch,symmetry)

    #check whether two methods yield same result
    #caxis==axis
    #resort points
    new_order=[]

    for idx in axis
        append!(new_order,idx)
    end

    for idx in bloch
        if !(idx in new_order)
            append!(new_order,idx)
        end
    end

    for idx in 1:npoints
        if !(idx in new_order) && !(idx in symmetry)
            append!(new_order,idx)
        end
    end

    for idx in symmetry
        if !(idx in new_order)
            append!(new_order,idx)
        end
    end

    trace_order=zeros(Int64,npoints) #TODO: loop is not so nice but does the job...
    for idx = 1: length(new_order)
        trace_order[idx]=findfirst(x->x==idx,new_order)
    end
    #next mirror at symmetry

    # compute symmetry plane
    link_triangles_to_tetrahedra!(mesh)
    tetmap=mesh.tri2tet

    smplx_idx=mesh.domains[sym_name]["simplices"][1]
    smplx=mesh.triangles[smplx_idx]
    tri=mesh.points[:,smplx]
    pln=three_points_to_plane(tri)
    pnt=find_foot_of_perpendicular([0.,0.,0.],pln)
    testpoint_idx=find_testpoint_idx(smplx,mesh.tetrahedra[tetmap[smplx_idx]])
    pln=make_normal_outwards(pln,mesh.points[:,testpoint_idx])

    # compute bloch plane
    smplx_idx=mesh.domains[blch_name]["simplices"][1]
    smplx=mesh.triangles[smplx_idx]
    tri2=mesh.points[:,smplx]
    bpln=three_points_to_plane(tri2)
    testpoint_idx=find_testpoint_idx(smplx,mesh.tetrahedra[tetmap[smplx_idx]])
    bpln=make_normal_outwards(bpln,mesh.points[:,testpoint_idx])
    #little tests
    #check=is_point_in_plane(tri2[:,2],pln)
    #tri_refl=reflect_point_at_plane(tri2[:,2],pln)

    # do the reflection
    nbloch=length(bloch)
    naxis=length(axis)
    nsymmetry=length(symmetry)
    nxsymmetry=nsymmetry-naxis
    nxbloch=nbloch-naxis
    nbody=npoints-nbloch-nxsymmetry
    shiftbody=npoints-nbloch
    nxsector=nxbloch+nbody+nxsymmetry+nbody
    nsector=nxsector+naxis
    #initialize point array
    points=zeros(3,2*npoints-nsymmetry)
    points[:,1:npoints]=mesh.points[:,new_order]
    # reflect body points
    for idx=nbloch+1:npoints-nxsymmetry
        pnt=points[:,idx]
        points[:,idx+shiftbody]=reflect_point_at_plane(pnt,pln)
    end
    #reflect exclusive bloch points
    for idx=naxis+1:nbloch
        pnt=points[:,idx]
        points[:,idx+nxsector]=reflect_point_at_plane(pnt,pln)
    end

    #get degree of symmetry
    phi=LinearAlgebra.dot(pln[1:3],-bpln[1:3])#/(LinearAlgebra.norm(pln[1:3])*LinearAlgebra.norm(bpln[1:3]))
    phi=acos(phi)##TODO: get direction
    DOS=round(Int,pi/phi)
    #create/virtualize full mesh

    p,n=find_intersection_of_two_planes(pln,bpln)
    phi=2pi/DOS
    R=create_rotation_matrix_around_axis(n,phi)
    #
    #rot_points=R*(points.-p).+p #Test
    #Full mesh
    nfpoints=naxis+(nxbloch+nbody+nxsymmetry+nbody)*DOS

    if unit
        fpoints=points
        DOS_lim=1
        file_txt="unit from $(mesh.file)"
    else
        DOS_lim=DOS
        fpoints=zeros(3,nfpoints)

        fpoints[:,1:nsector]=points[:,1:nsector]
        for idx in 1:DOS-1
            R=create_rotation_matrix_around_axis(n,idx*phi)
            fpoints[:,naxis+nxsector*idx+1:naxis+nxsector*(idx+1)]=R*(points[:,naxis+1:naxis+nxsector].-p).+p
        end
        file_txt="extended from $(mesh.file)"
    end
    #domains=Dict()
    #triangles=Array{UInt32,1}[]

    reflected_index(idx) = get_reflected_index(idx,naxis,nxbloch,nbody,shiftbody,nxsymmetry,nxsector)
    # function get_reflected_index(idx)
    #     if idx <= naxis
    #         idx=idx
    #     elseif idx <= nbloch
    #         idx=idx+nxsector
    #     elseif idx <= nbloch+nbody
    #             idx=idx+shiftbody
    #     elseif idx <=  nsector-nbody #naxis+nxbloch+nbody+nxsymmetry
    #         idx=idx
    #     else
    #         idx=0 #error flag
    #     end
    #     return idx
    # end

    rotated_index(idx,sector)=get_rotated_index(idx,sector, naxis, nxsector, DOS)
    # nasector=nxsector*DOS
    # function get_rotated_index(idx,sector)
    #     if idx <= naxis
    #         idx=idx
    #     else
    #         idx=idx+nxsector*sector
    #     end
    #     #modulo for periodicity
    #     idx=((idx-naxis-1)%nasector)+1+naxis
    #     if idx<0
    #         idx+=nasector
    #     end
    #     return idx
    # end

    tetrahedra=Array{UInt32,1}[] #TODO: sectorwise arrangement of tetrahedra
    for (smplx_idx,tet) in enumerate(mesh.tetrahedra)
        #println(idx," ", length(tetrahedra))
        tet=trace_order[tet]
        rtet=reflected_index.(tet)

        for sector=0:DOS_lim-1
            insert_smplx!(tetrahedra,rotated_index.(tet, sector))
            insert_smplx!(tetrahedra,rotated_index.(rtet, sector))
        end
    end

    triangles=Array{UInt32,1}[] #TODO: sectorwise arrangement of triangles
    for (smplx_idx,smplx) in enumerate(mesh.triangles)
        if !(smplx_idx in mesh.domains[sym_name]["simplices"]) && (!(smplx_idx in mesh.domains[blch_name]["simplices"]) || unit)
            smplx=trace_order[smplx]
            rsmplx=reflected_index.(smplx)
            for sector=0:DOS_lim-1
                insert_smplx!(triangles,rotated_index.(smplx, sector))
                insert_smplx!(triangles,rotated_index.(rsmplx, sector))
            end
        end
    end

    ## sectorwise line handling
    lines=Array{UInt32,1}[]
    for (smplx_idx, smplx) in enumerate(mesh.lines)
        smplx=trace_order[smplx]
        insert_smplx!(lines,smplx)
        if !all(smplx.<=nbloch)
            rsmplx=reflected_index.(smplx)
            insert_smplx!(lines,rsmplx)
        end
    end

    #count second order dof on axis and bloch boundary
    naxis_ln=0
    nbloch_ln=0
    for (idx,ln) in enumerate(lines)
            if naxis_ln==0 &&  !all(ln.<=naxis)
                naxis_ln=idx-1
            end
            if nbloch_ln==0 &&  !all(ln.<=naxis+nxbloch)
                nbloch_ln=idx-1
            end
    end
    nxbloch_ln=nbloch_ln-naxis_ln
    nsector_ln=length(lines)#-naxis_ln
    nxsector_ln=nsector_ln-naxis_ln

    # extend to full mesh
    if unit #just create bloch image boundary
        for ln in lines[naxis_ln+1:nbloch_ln]
            append!(lines,[rotated_index.(ln,1)])
        end
    else #create full domain
        for sector=1:DOS-1
            for lns in lines[naxis_ln+1:nsector_ln]
                append!(lines,[rotated_index.(lns,sector)])
            end
        end
    end
    ##build domains
    domains=Dict()
    #build unit cell

    #build full circumference #TODO: virtualize domains to save memeory and building time
    for (dom, copy_degree) in doms
        dim=mesh.domains[dom]["dimension"]
        if dim == 3
            list_of_simplices=mesh.tetrahedra
            full_list_of_simplices=tetrahedra
        elseif dim == 2
            list_of_simplices=mesh.triangles
            full_list_of_simplices=triangles
        elseif dim == 1
            list_of_simplices=mesh.lines
            full_list_of_simplices=nothing#lines
            println("Error: Dimension 1 not implemented for extensible domain")
            #TODO: Implement
        end

        if copy_degree==:full
            domains[dom]=Dict()
            domains[dom]["dimension"]=dim
            domains[dom]["simplices"]=UInt64[]
        elseif copy_degree==:unit
            for sector=0:DOS_lim-1
                domains[dom*"#$sector"]=Dict()
                domains[dom*"#$sector"]["dimension"]=dim
                domains[dom*"#$sector"]["simplices"]=UInt64[]
            end
        elseif copy_degree==:half
            for sector=0:DOS_lim-1
                domains[dom*"#$sector.0"]=Dict()
                domains[dom*"#$sector.0"]["dimension"]=dim
                domains[dom*"#$sector.0"]["simplices"]=UInt64[]
                domains[dom*"#$sector.1"]=Dict()
                domains[dom*"#$sector.1"]["dimension"]=dim
                domains[dom*"#$sector.1"]["simplices"]=UInt64[]
            end
        else
            println("Error: copy_degree '$copy_degree' not supported! Use either
            'full', 'unit', or 'half'.")
            return
        end
        for smplx_idx in mesh.domains[dom]["simplices"]
            smplx=list_of_simplices[smplx_idx]
            smplx=trace_order[smplx]
            rsmplx=reflected_index.(smplx)
            for sector=0:DOS_lim-1
                idx=find_smplx(full_list_of_simplices,rotated_index.(smplx, sector))
                ridx=find_smplx(full_list_of_simplices,rotated_index.(rsmplx, sector))
                if copy_degree==:full
                    append!(domains[dom]["simplices"],idx)
                    append!(domains[dom]["simplices"],ridx)
                elseif copy_degree==:unit
                    append!(domains[dom*"#$sector"]["simplices"],idx)
                    append!(domains[dom*"#$sector"]["simplices"],ridx)
                elseif copy_degree==:half
                    append!(domains[dom*"#$sector.0"]["simplices"],idx)
                    append!(domains[dom*"#$sector.1"]["simplices"],ridx)
                end
            end
        end
    end

    #create new mesh object
    #symmetry_info=(DOS,naxis,nxbloch,nbody,nxsymmetry,n,naxis_ln,nxbloch_ln,nxsector_ln)
    symmetry_info=SymInfo(DOS,naxis,nxbloch,nbody,shiftbody,nxsymmetry,nxsector,naxis_ln,nxbloch_ln,nxsector_ln,0,0,n,p,unit)
    tri2tet=zeros(UInt32,length(triangles))
    tri2tet[1]=0xffffffff #magic sentinel value to check whether list is linked
    full_mesh=Mesh(mesh.file,fpoints,lines,triangles, tetrahedra, domains, file_txt, tri2tet, symmetry_info)
    return full_mesh
end
## wrapper for finder routines to Mesh type
"""
    sidx=get_point_sector(mesh::Mesh,pnt_idx)

Get sector idx `sidx`of point with index `pnt_idx`in mesh `mesh`.
"""
function get_point_sector(mesh::Mesh,pnt_idx)
    if mesh.dos==1
        return 0
    else
        return get_point_sector(pnt_idx, mesh.dos.naxis, mesh.dos.nxsector)
    end
end

"""
    sidx = get_line_sector(mesh::Mesh, ln)

Get sector idx `sidx`of line `ln` in mesh `mesh`.
"""
function get_line_sector(mesh::Mesh, ln)
    if mesh.dos==1
        return 0
    else
        return get_line_sector(ln,mesh.dos.naxis,mesh.dos.nxsector)
    end
end
# line finder for mesh TODO: smplx_finder for mesh
"""
    ln_idx = get_line_idx(mesh::Mesh, ln)

Get index `ln_idx` of line `ln` in mesh `mesh`.
"""
function get_line_idx(mesh::Mesh, ln)
    if mesh.dos==1
        ln_idx=find_smplx(mesh.lines,ln)
    else
        ln_idx=get_line_idx(mesh.lines,ln,mesh.dos.naxis,mesh.dos.nxsector,mesh.dos.naxis_ln,mesh.dos.nxsector_ln,mesh.dos.DOS)
    end
    return ln_idx
end


text:
function get_friends_ids(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friends/ids.json", options))
end

function get_friends_list(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/friends/list.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
###
# helpers
###

# make sure a passed distance matrix is a square
function check_square(m, msg)
    n = size(m, 1)
    if n != size(m, 2)
        error(msg)
    end
    return n
end

"""
    legal_circuit(circuit::AbstractArray{<:Integer})

Check that an array of integers is a valid circuit. A valid circuit over `n` locations has
length `n+1`. The first `n` entries are a permutation of `1, ..., n`, and the `(n+1)`-st entry
is equal to the first entry.
"""
function legal_circuit(circuit)
    n = length(circuit) - 1
    return circuit[1] == circuit[end] && sort(circuit[1:(end-1)]) == 1:n
end

"""
    repetitive_heuristic(distmat::Matrix, heuristic::Function, repetitive_kw::Symbol; keywords...)

For each `i` in `1:n`, run `heuristic` with keyword `repetitive_kw` set to `i`. Return the tuple
`(bestpath, bestcost)`. By default, `repetitive_kw` is `:firstcity`, which is appropriate for the
`farthest_insertion`, `cheapest_insertion`, and `nearest_neighbor` heuristics.

Any keywords passed to `repetitive_heuristic` are passed along to each call of `heuristic`. For
example, `repetitive_heuristic(dm, nearest_neighbor; do2opt = true)` will perform 2-opt for
each of the `n` nearest neighbor paths.

!!! note
    The repetitive heuristic calls are parallelized with threads. For optimum speed make sure
    Julia is running with multiple threads.
"""
function repetitive_heuristic(dm::AbstractMatrix{T},
                              heuristic::Function,
                              repetitive_kw = :firstcity;
                              kwargs...) where {T<:Real}
    # call the heuristic with varying starting cities
    n = size(dm, 1)
    results_list = Vector{Tuple{Vector{Int}, T}}(undef, n)
    Threads.@threads for i in 1:n
        results_list[i] = heuristic(dm; kwargs..., repetitive_kw => i)
    end

    bestind, bestcost = 1, results_list[1][2]
    for i in 2:n
        if results_list[i][2] < bestcost
            bestind, bestcost = i, results_list[i][2]
        end
    end
    return results_list[bestind]
end

# helper for readable one-line path costs
# optionally specify the bounds for the subpath we want the cost of
# defaults to the whole path
# but when calculating reversed path costs can help to have subpath costs
function pathcost(distmat::AbstractMatrix{T}, path::AbstractArray{S},
                  lb::Int = 1, ub::Int = length(path)) where {T<:Real, S<:Integer}
    cost = zero(T)
    for i in lb:(ub - 1)
        @inbounds cost += distmat[path[i], path[i+1]]
    end
    return cost
end

"Compute the change in cost from reversing the subset of the path from indices
`revLow` to `revHigh`, inclusive."
function pathcost_rev_delta(distmat::AbstractMatrix{T}, path::AbstractArray{S},
                            revLow::Int, revHigh::Int) where {T<:Real, S<:Integer}
    cost_delta = zero(T)
    # if there's an initial unreversed section
    if revLow > 1
        # new step onto the reversed section
        @inbounds cost_delta += distmat[path[revLow - 1], path[revHigh]]
        # no longer pay the cost of old step onto the reversed section
        @inbounds cost_delta -= distmat[path[revLow - 1], path[revLow]]
    end
    # new cost of the reversed section
    for i in revHigh:-1:(revLow + 1)
        @inbounds cost_delta += distmat[path[i], path[i-1]]
    end
    # no longer pay the forward cost of the reversed section
    for i in revLow:(revHigh - 1)
        @inbounds cost_delta -= distmat[path[i], path[i+1]]
       end
    # if there's an unreversed section after the reversed bit
    if revHigh < length(path)
        # new step out of the reversed section
        @inbounds cost_delta += distmat[path[revLow], path[revHigh + 1]]
        # no longer pay the old cost of stepping out of the reversed section
        @inbounds cost_delta -= distmat[path[revHigh], path[revHigh + 1]]
    end
    return cost_delta
end
"Specialized for symmetric matrices: compute the change in cost from
reversing the subset of the path from indices `revLow` to `revHigh`,
inclusive."
function pathcost_rev_delta(distmat::Symmetric, path::AbstractArray{S},
                      revLow::Int, revHigh::Int) where {S<:Integer}
    cost_delta = zero(eltype(distmat))
    # if there's an initial unreversed section
    if revLow > 1
        # new step onto the reversed section
        @inbounds cost_delta += distmat[path[revLow - 1], path[revHigh]]
        # no longer pay the cost of old step onto the reversed section
        @inbounds cost_delta -= distmat[path[revLow - 1], path[revLow]]
    end
    # The actual cost of walking along the reversed section doesn't change
    # because the distance matrix is symmetric.
    # if there's an unreversed section after the reversed bit
    if revHigh < length(path)
        # new step out of the reversed section
        @inbounds cost_delta += distmat[path[revLow], path[revHigh + 1]]
        # no longer pay the old cost of stepping out of the reversed section
        @inbounds cost_delta -= distmat[path[revHigh], path[revHigh + 1]]
    end
    return cost_delta
end

"Due to floating point imprecision, various path refinement heuristics may get
stuck in infinite loops as both doing and un-doing a particular change apparently
improves the path cost. For example, floating point error may suggest that both
reversing and un-reversing a path for a symmetric TSP instance improves the cost
slightly. To avoid such false-improvement infinite loops, limit refinement heuristics
to improvements with some minimum magnitude defined by the element type of the
distance matrix."
improvement_threshold(T::Type{<:Integer}) = one(T)
improvement_threshold(T::Type{<:AbstractFloat}) = sqrt(eps(one(T)))
improvement_threshold(T::Type{<:Real}) = sqrt(eps(1.0))

"Cost of inserting city `k` after index `after` in path `path` with costs `distmat`."
function inscost(
    k::Int,
    after::Int,
    path::AbstractArray{<:Integer},
    distmat::AbstractMatrix{<:Real}
)
    return distmat[path[after], k] +
           distmat[k, path[after + 1]] -
           distmat[path[after], path[after + 1]]
end


text:
function printMatrix( A::Array{ComplexF64,2} )
    Nrows = size(A)[1]
    Ncols = size(A)[2]
    for ir = 1:Nrows
        for ic = 1:Ncols
            @printf("(%6.3f,%6.3f) ", real(A[ir,ic]), imag(A[ir,ic]))
        end
        @printf("\n")
    end
end

function printMatrix( A::Array{Float64,2} )
    Nrows = size(A)[1]
    Ncols = size(A)[2]
    for ir = 1:Nrows
        for ic = 1:Ncols
            @printf("%6.3f ", A[ir,ic])
        end
        @printf("\n")
    end
end


text:
function generate_scenario(scenario, ego_v, hit_point)

    ego_y = 0.0
    ped_x = 100.0
    ped_y = 0.0

    if scenario == "CPCN"
        scenario_id = 1
        ped_v = 5/3.6
        obstacles = [ConvexPolygon([VecE2(ped_x-VehicleDef().length-1, ped_y-2),
                VecE2(ped_x-VehicleDef().length-1, ped_y-2-VehicleDef().width),
                VecE2(ped_x-1, ped_y-2-VehicleDef().width),
                VecE2(ped_x-1, ped_y-2)],4)]

    elseif scenario == "CPAN25"
        scenario_id = 2
        ped_v = 5 /3.6
        obstacles = []

    elseif scenario == "CPAN75"
        scenario_id = 3
        ped_v = 5 /3.6
        hit_point = hit_point + 50
        obstacles = []

    elseif scenario == "CPFA"
        scenario_id = 4
        ped_v = 8 / 3.6
        obstacles = []

    elseif scenario == "FP"
        scenario_id = 9
        ped_v = 5 /3.6
        obstacles = []
    else
        scenario_id = -1
        ped_v = 5 / 3.6
        obstacles = []
    end

    # fix values
    ped_right_side = true
    ped_theta = 1.5707
    ped_t_collision = 3.0

    ped_y_end = VehicleDef().width * hit_point / 100 - VehicleDef().width/2
    if (ped_right_side)
        ped_y_start = ped_y_end - ped_t_collision * ped_v
    else
        ped_y_start = ped_y_end + ped_t_collision * ped_v
    end
    ped_y = ped_y_start
    ego_x = ped_x - ego_v * ped_t_collision - VehicleDef().length/2;
    #println("Scenario: ", scenario, " v_ego=", ego_v*3.6, "km/h v_ped=", ped_v*3.6, "km/h HP=", hit_point, " ped_y: ", ped_y)
    return (ego_x, ego_y, ego_v, ped_x, ped_y, ped_v, ped_theta, obstacles, scenario_id)
end

function evaluate_scenario(ego_x, ego_y, ego_v, ped_x, ped_y, ped_v, ped_theta, obstacles)

    timestep = 0.05

    params = CrosswalkParams()
    if (length(obstacles) > 0)
        params.obstacles = obstacles
        params.obstacles_visible = true
    else
        params.obstacles_visible = false
    end
    params.roadway_length = 400.0

    env = CrosswalkEnv(params);

    ego_id = EGO_ID
    ped_id = 2
    ped2_id = 3
    ped3_id = 4

    # Car definition
    ego_initial_state = VehicleState(VecSE2(ego_x, ego_y, 0.), env.roadway.segments[1].lanes[1], env.roadway, ego_v)
    ego = Vehicle(ego_initial_state, VehicleDef(), ego_id)

    # Pedestrian definition using our new Vehicle definition
    ped_initial_state = VehicleState(VecSE2(ped_x, ped_y, ped_theta), env.crosswalk, env.roadway, ped_v)
    ped = Vehicle(ped_initial_state, AutomotivePOMDPs.PEDESTRIAN_DEF, ped_id)

    ped2 = Vehicle(VehicleState(VecSE2(90., 3., -1.57), env.crosswalk, env.roadway, 0.), AutomotivePOMDPs.PEDESTRIAN_DEF, ped2_id)
    ped3 = Vehicle(VehicleState(VecSE2(103., 10., -1.57), env.crosswalk, env.roadway, 1.), AutomotivePOMDPs.PEDESTRIAN_DEF, ped3_id)

    scene = Scene()
    push!(scene, ego)
    push!(scene, ped)
    push!(scene, ped2)
    push!(scene, ped3)

    pos_noise = 0.0
    vel_noise = 0.0
    false_positive_rate = 0.0
    false_negative_rate = 0.0
    rng = MersenneTwister(1);
    sensor = AutomotiveSensors.GaussianSensor(AutomotiveSensors.LinearNoise(10, pos_noise, 0.0),
                     AutomotiveSensors.LinearNoise(10, vel_noise, 0.0), false_positive_rate, false_negative_rate, rng)

    # define a model for each entities present in the scene
    models = Dict{Int, DriverModel}()

    models[ego_id] = EmergencyBrakingSystem.EmergencySystem(a=LatLonAccel(0.0, 0.0),
        env=env,
        sensor=sensor,
        obstacles=env.obstacles,
        SAFETY_DISTANCE_LON=1.0,
        AX_MAX=-10.0,
        THRESHOLD_COLLISION_RATE = 0.5,
        THRESHOLD_TIME_TO_REACT = 0.99,
        timestep=timestep)

    models[ped_id] = ConstantPedestrian(v_desired=ped_v, dawdling_amp=0.0) # dumb model
    models[ped2_id] = ConstantPedestrian(v_desired=0.0, dawdling_amp=0.05) # dumb model
    models[ped3_id] = ConstantPedestrian(v_desired=1.0, dawdling_amp=0.05) # dumb model

    nticks = 100
    sensor_observations = [Vehicle[]]
    ego_vehicle = Vehicle[]
    ego_a = Float64[]
    collision = Bool[]

    collision_rate = Float64[]  # not used
    ttc = Float64[]             # not used
    risk = Float64[]
    emergency_brake_request = Bool[]      # not used
    prediction_obstacle = Vector{Array{Float64}}()  # not used

    obs_callback = (EmergencyBrakingSystem.ObservationCallback(sensor_observations, ego_vehicle, ego_a, collision, collision_rate, ttc, risk, emergency_brake_request, prediction_obstacle),)

    scenes = simulate(scene, env.roadway, models, nticks, timestep, callbacks=obs_callback)

    return (scenes, timestep, env, sensor, sensor_observations, ego_vehicle, ego_a, collision, collision_rate, ttc, risk, emergency_brake_request, prediction_obstacle)

end

function evaluateScenarioMetric(ego_vehicle, emergency_brake_request, ego_a, collision_vector, ped_x)

    emergency_brake_intervention = false
    collision = false
    dv_collision = 0.

    v = []
    a = []
    a_jerk = 0.
    a_last = 0.
    for i=1:length(ego_vehicle)
        # do not consider acceleration part, only until the pedestrian is reached
        if (ego_vehicle[i].state.v > 0 && ego_vehicle[i].state.posG.x < ped_x - 3 )
        #if (ego_vehicle[i].state.v > 0)
            push!(v,ego_vehicle[i].state.v)
        end

        if (ego_a[i] != 0. && ego_a[i] < 0.)
            push!(a,ego_a[i])
        end

        if ( a_last != ego_a[i] )
            a_jerk = a_jerk + abs(ego_a[i] - a_last)
        end
        a_last = ego_a[i]

        if length(emergency_brake_request) > 0 && emergency_brake_request[i] == 1
            emergency_brake_intervention = true
        end

    end
    v_mean = mean(v)
    if ( length(a) > 0)
        a_mean = mean(a)
        a_min = minimum(a)
    else
        a_mean = 0.
        a_min = 0.
    end

    if (collision_vector[end] == true && ego_vehicle[end].state.posG.x < ped_x - 2.)   # ignore collisions when the pedestrian runs into the side of the car (ego vehicle too slow)
        collision = true
        dv_collision = ego_vehicle[end].state.v
    end

    return (collision, emergency_brake_intervention, dv_collision, v_mean, a_mean, a_jerk, a_min)
end

function evaluateScenariosMetric(results)

    sum_collisions = 0
    sum_eb = 0

    dv = []
    v_mean = []
    a_mean = []
    a_jerk = []
    a_min = []

    for i=1:length(results)
        #println(results[i])
        if ( results[i][4] == 1.0 )
            sum_collisions += 1
        end
        if ( results[i][5] == 1.0 )
            sum_eb += 1
        end

        push!(dv,results[i][6])
        push!(v_mean,results[i][7])
        push!(a_mean,results[i][8])
        push!(a_jerk,results[i][9])
        push!(a_min,results[i][10])

    end
    dv = mean(dv)
    v_mean = mean(v_mean)
    a_mean = mean(a_mean)
    a_jerk = mean(a_jerk)
    a_min = mean(a_min)

    return (sum_collisions, sum_eb, dv, v_mean, a_mean, a_jerk, a_min)
end


text:
"""
    HorizontalAverage

Targeted at AtmosLES configurations, produces a vector of length `z`
of averages for the horizontal planes in the domain, where `z` is the
interpolated grid height or the DG grid`s `x3id` over the nodes and
elements (duplicates are not removed).
"""
abstract type HorizontalAverage <: DiagnosticVar end
function dv_HorizontalAverage end

# TODO: replace these with a `dv_array_dims` that takes `nvars`
# and returns the dims for the array? Or create the array? Use
# `Array`? `similar`?

function dv_dg_points_length(
    ::ClimateMachineConfigType,
    ::Type{HorizontalAverage},
)
    :(Nqk)
end
function dv_dg_points_index(
    ::ClimateMachineConfigType,
    ::Type{HorizontalAverage},
)
    :(k)
end

function dv_dg_elems_length(
    ::ClimateMachineConfigType,
    ::Type{HorizontalAverage},
)
    :(nvertelem)
end
function dv_dg_elems_index(
    ::ClimateMachineConfigType,
    ::Type{HorizontalAverage},
)
    :(ev)
end

function dv_dg_dimnames(::ClimateMachineConfigType, ::Type{HorizontalAverage})
    ("z",)
end
function dv_dg_dimranges(::ClimateMachineConfigType, ::Type{HorizontalAverage})
    z = quote
        ijk_range = 1:Nqh:npoints
        e_range = 1:nvertelem
        reshape(grid.vgeo[ijk_range, grid.x3id, e_range], :)
    end
    (z,)
end

function dv_i_dimnames(::AtmosLESConfigType, ::Type{HorizontalAverage})
    ("z",)
end
function dv_i_dimnames(::AtmosGCMConfigType, ::Type{HorizontalAverage})
    ("level",)
end

function dv_op(::ClimateMachineConfigType, ::Type{HorizontalAverage}, lhs, rhs)
    :($lhs += MH * $rhs)
end
function dv_reduce(
    ::ClimateMachineConfigType,
    ::Type{HorizontalAverage},
    array_name,
)
    quote
        MPI.Reduce!($array_name, +, 0, mpicomm)
        if mpirank == 0
            for v in 1:size($array_name, 2)
                $(array_name)[:, v, :] ./= DiagnosticsMachine.Collected.ΣMH_z
            end
        end
    end
end

macro horizontal_average(impl, config_type, name, scale = nothing)
    iex = quote
        $(generate_dv_interface(:HorizontalAverage, config_type, name))
        $(generate_dv_function(:HorizontalAverage, config_type, name, impl))
        $(generate_dv_scale(:HorizontalAverage, config_type, name, scale))
    end
    esc(MacroTools.prewalk(unblock, iex))
end

"""
    @horizontal_average(
        impl,
        config_type,
        name,
        units,
        long_name,
        standard_name,
        scale = nothing,
    )

Define `name`, a horizontal average diagnostic variable for `config_type`
with the specified attributes and the given implementation. In order to
produce a density-weighted average, `scale` must be specified as the name
of the horizontal average diagnostic variable for density.

# Example

```julia
@horizontal_average(
    AtmosLESConfigType,
    w_ht_sgs,
    "kg kg^-1 m s^-1",
    "vertical sgs flux of total specific enthalpy",
    "",
    rho,
) do (atmos::AtmosModel, states::States, curr_time, cache)
    ts = get!(cache, :ts) do
        recover_thermo_state(atmos, states.prognostic, states.auxiliary)
    end
    D_t = get!(cache, :D_t) do
        _, D_t, _ = turbulence_tensors(
            atmos,
            states.prognostic,
            states.gradient_flux,
            states.auxiliary,
            curr_time,
        )
        D_t
    end
    d_h_tot = -D_t .* states.gradient_flux.∇h_tot
    d_h_tot[end] * states.prognostic.ρ
end
```
"""
macro horizontal_average(
    impl,
    config_type,
    name,
    units,
    long_name,
    standard_name,
    scale = nothing,
)
    iex = quote
        $(generate_dv_interface(
            :HorizontalAverage,
            config_type,
            name,
            units,
            long_name,
            standard_name,
        ))
        $(generate_dv_function(:HorizontalAverage, config_type, name, impl))
        $(generate_dv_scale(:HorizontalAverage, config_type, name, scale))
    end
    esc(MacroTools.prewalk(unblock, iex))
end


text:
using Turing, ContinuousBenchmarks

data = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]

@model constrained_test(obs) = begin
    p ~ Beta(2,2)
    for i = 1:length(obs)
        obs[i] ~ Bernoulli(p)
    end
    p
end

bench_res = @tbenchmark(HMC(1000, 1.5, 3), constrained_test, data)

# bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"]
logd = build_logd("Dummy-Benchmark", bench_res...)

print_log(logd)
send_log(logd)


text:
# Run `asciinema rec example_slow.json`, make the recording, then run the following to speed it up and make an svg.
# Uses https://github.com/marionebl/svg-term-cli and https://asciinema.org/

using Printf

function speed_up_asciicinema(io, input_file_name; speedup=10.0)
    f = str -> begin
        num = parse(Float64,str)
        "[$(@sprintf("%2.6f", num / speedup)), \"o\","
    end
    for line in readlines(input_file_name)
        out =  findfirst(r"\[.*, \"o\",", line)
        if !isnothing(out) && out.start == 1
            m = match(r"\[(.*), \"o\",", line)
            new_line = replace(line, r"\[(.*), \"o\"," => f(m.captures[]), count=1 )
        else
            new_line = line
        end
        println(io, new_line)
    end
end

open("example.json", "w") do f
    speed_up_asciicinema(f, "example_slow.json")
end

run(`cat example.json | svg-term --out=example.svg`)


text:
using Test: @test, @testset
import Zygote
import WeightNormalizedRBMs as wRBMs

@testset "gu2w, w2gu" begin
    g = abs.(randn(1,1,4,2))
    u = randn(5,3,4,2)
    nt = wRBMs.gu2w(g, u)
    @test nt.un ≈ sqrt.(sum!(similar(g), abs2.(u)))
    @test nt.w ≈ g .* u ./ nt.un
    @test sqrt.(sum!(similar(g), abs2.(nt.w))) ≈ g

    nt1 = wRBMs.w2gu(nt.w, nt.un)
    @test nt1.g ≈ g
    @test nt1.u ≈ u
end

@testset "∂wnorm" begin
    foo(w) = sum(sin.(w))
    function foo(g, u)
        un = sqrt.(sum(abs2.(u); dims=findall(isone, size(g))))
        w = g .* u ./ un
        return foo(w)
    end
    g = abs.(randn(1,1,4,2))
    u = randn(5,3,4,2)
    w = wRBMs.gu2w(g, u).w
    ∂w, = Zygote.gradient(foo, w)
    ∂g, ∂u = Zygote.gradient(foo, g, u)
    ∂ = wRBMs.∂wnorm(∂w, g, u)
    @test ∂.g ≈ ∂g
    @test ∂.u ≈ ∂u
end


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "6-5-Y scenario" begin

@time for Y in 7:9
    filepath = "./data/quick_adjacency_decomposition/"
    filename = "6-5-$Y"

    println("$filename scenario")

    PM = LocalSignaling(6,Y,5)
    BG_seed = BellGame(cat(
        cat(
            ones(Int64,Y-6),
            zeros(Int64,Y-6,5),
            dims=2
        ),
        [1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1],
        dims=1
    ), 5)

    vertices = LocalPolytope.vertices(PM)
    canonical_facets = LocalPolytope.adjacency_decomposition(vertices, BG_seed, PM, max_vertices = 100)
    games = collect(keys(canonical_facets))

    pretty_print_txt(games, filepath*filename)
    write_ieq(filepath*"ieq/"*filename, convert(IEQ, games))
end

end


text:
function fit_SAIM(file_path::String, file_name::String, optic::SAIMOptics, angles::AbstractArray,
	init_params::AbstractArray, lower_bounds::AbstractArray, upper_bounds::AbstractArray; glb::Bool=false,
	step::Float64=40.0, show::Bool=false, color::Int64=1)

    angles_rad = angles*2*pi/360

    if color == 1
        if glb	#Run a global grid search when true
            fit_1c_global(file_path, file_name, optic, angles_rad, init_params, lower_bounds,
            upper_bounds, step, show)
        else	#Conduct a local search
            fit_1c_local(file_path, file_name, optic, angles_rad, init_params, lower_bounds,
            upper_bounds, show)
        end
    elseif color == 2
        if glb	#Run a global grid search when true
            fit_2c_global(file_path, file_name, optic, angles_rad, init_params, lower_bounds,
            upper_bounds, step, show)
        else	#Conduct a local search
            fit_2c_local(file_path, file_name, optic, angles_rad, init_params, lower_bounds,
            upper_bounds, show)
        end
    elseif color == 3
        if glb	#Run a global grid search when true
            fit_3c_global(file_path, file_name, optic, angles_rad, init_params,
            lower_bounds, upper_bounds, step, show)
        else	#Conduct a local search
            fit_3c_local(file_path, file_name, optic, angles_rad, init_params,
            lower_bounds, upper_bounds, show)
        end
    else
        println("Only 1, 2, and 3 colors supported")
    end
end


text:
# ******************************************************************************
# TXT
# ******************************************************************************
struct PlotVar{S, T}
    name :: S
    siunit :: S
    unit :: T
end

struct OutputTXT{S}
    fname :: S
end

# ------------------------------------------------------------------------------
# TXT T
# ------------------------------------------------------------------------------
function OutputTXT(fname::String, grid::GridT, field::FieldT, neu, zu)
    (; tu) = grid
    (; Iu) = field
    plotvars = (
        PlotVar("z", "m", zu),
        PlotVar("Imax", "W/m^2", Iu),
        PlotVar("nemax", "1/m^3", neu),
        PlotVar("tau", "s", tu),
        PlotVar("F", "J/m^2", tu * Iu),
    )
    write_txt_header(fname, plotvars)
    return OutputTXT(fname)
end

function write_txt(out::OutputTXT, analyzer::FieldAnalyzerT)
    (; z, Imax, nemax, tau, F) = analyzer
    fp = open(out.fname, "a")
    write(fp, "  ")
    @printf(fp, "%18.12e ", z)
    @printf(fp, "%18.12e ", Imax)
    @printf(fp, "%18.12e ", nemax)
    @printf(fp, "%18.12e ", tau)
    @printf(fp, "%18.12e ", F)
    write(fp, "\n")
    close(fp)
    return nothing
end

# ------------------------------------------------------------------------------
# TXT RT
# ------------------------------------------------------------------------------
function OutputTXT(fname::String, grid::GridRT, field::FieldRT, neu, zu)
    (; ru, tu) = grid
    (; Iu) = field
    plotvars = (
        PlotVar("z", "m", zu),
        PlotVar("Imax", "W/m^2", Iu),
        PlotVar("Fmax", "J/m^2", tu * Iu),
        PlotVar("nemax", "1/m^3", neu),
        PlotVar("rad", "m", ru),
        PlotVar("tau", "s", tu),
        PlotVar("W", "J", tu * ru^2 * Iu),
    )
    write_txt_header(fname, plotvars)
    return OutputTXT(fname)
end

function write_txt(out::OutputTXT, analyzer::FieldAnalyzerRT)
    (; z, Imax, Fmax, nemax, rad, tau, W) = analyzer
    fp = open(out.fname, "a")
    write(fp, "  ")
    @printf(fp, "%18.12e ", z)
    @printf(fp, "%18.12e ", Imax)
    @printf(fp, "%18.12e ", Fmax)
    @printf(fp, "%18.12e ", nemax)
    @printf(fp, "%18.12e ", rad)
    @printf(fp, "%18.12e ", tau)
    @printf(fp, "%18.12e ", W)
    write(fp, "\n")
    close(fp)
    return nothing
end

# ------------------------------------------------------------------------------
# TXT tools
# ------------------------------------------------------------------------------
function write_txt_header(fname, plotvars)
    fp = open(fname, "w")

    # write names:
    write(fp, "#")
    for pvar in plotvars
        @printf(fp, " %-18s", pvar.name)
    end
    write(fp, "\n")

    # write SI units:
    write(fp, "#")
    for pvar in plotvars
        @printf(fp, " %-18s", pvar.siunit)
    end
    write(fp, "\n")

    # write dimensionless units:
    write(fp, "#")
    for pvar in plotvars
        @printf(fp, " %-18s", pvar.unit)
    end
    write(fp, "\n")

    close(fp)
    return nothing
end

# ******************************************************************************
# BIN
# ******************************************************************************
mutable struct OutputBIN{T}
    fname :: String
    iout :: Int
    zout :: T
    dzout :: T
    zdata :: Bool
    izdata :: Int
end

# ------------------------------------------------------------------------------
# BIN T
# ------------------------------------------------------------------------------
function OutputBIN(fname::String, grid::GridT{T}, zu, z, dzout) where T
    (; tu, t) = grid

    fp = HDF5.h5open(fname, "w")
    group = HDF5.create_group(fp, "units")
    group["t"] = tu
    group["z"] = zu
    group = HDF5.create_group(fp, "grid")
    group["t"] = collect(t)
    group = HDF5.create_group(fp, "field")
    HDF5.close(fp)

    iout = 1
    zout = z

    zdata = false
    izdata = 1
    return OutputBIN{T}(fname, iout, zout, dzout, zdata, izdata)
end

function write_field(field::FieldT, group, dset)
    group[dset] = collect(field.E)
    return nothing
end

# ------------------------------------------------------------------------------
# BIN RT
# ------------------------------------------------------------------------------
function OutputBIN(fname::String, grid::GridRT{T}, zu, z, dzout) where T
    (; Nr, Nt, ru, tu, r, t) = grid
    Nwr = iseven(Nt) ? div(Nt, 2) : div(Nt+1, 2)

    fp = HDF5.h5open(fname, "w")
    group = HDF5.create_group(fp, "units")
    group["r"] = ru
    group["t"] = tu
    group["z"] = zu
    group = HDF5.create_group(fp, "grid")
    group["r"] = r
    group["t"] = collect(t)
    group = HDF5.create_group(fp, "field")
    group = HDF5.create_group(fp, "zdata")
    create_dataset(group, "z", T, ((1,), (-1,)), chunk=(100,))
    create_dataset(group, "Fr", T, ((1,Nr), (-1,Nr)), chunk=(100,Nr))
    create_dataset(group, "Si", T, ((1,Nwr+1), (-1,Nwr+1)), chunk=(100,Nwr+1))
    create_dataset(group, "ne", T, ((1,Nr), (-1,Nr)), chunk=(100,Nr))
    HDF5.close(fp)

    iout = 1
    zout = z

    zdata = true
    izdata = 1
    return OutputBIN{T}(fname, iout, zout, dzout, zdata, izdata)
end

function write_field(field::FieldRT, group, dset)
    group[dset, shuffle=true, deflate=9] = real.(collect(field.E))
    return nothing
end

function write_zdata(analyzer::FieldAnalyzerRT, group, iz)
    data = group["z"]
    HDF5.set_extent_dims(data, (iz,))
    data[iz] = analyzer.z

    data = group["Fr"]
    HDF5.set_extent_dims(data, (iz, length(analyzer.Fr)))
    data[iz, :] = collect(analyzer.Fr)

    data = group["Si"]
    Nt = length(analyzer.Si)
    Nwr = iseven(Nt) ? div(Nt, 2) : div(Nt+1, 2)
    HDF5.set_extent_dims(data, (iz, Nwr+1))
    data[iz, :] = collect(analyzer.Si)[Nwr:end]

    data = group["ne"]
    HDF5.set_extent_dims(data, (iz, length(analyzer.ne)))
    data[iz, :] = collect(analyzer.ne)
    return nothing
end

# ------------------------------------------------------------------------------
# BIN tools
# ------------------------------------------------------------------------------
function write_bin(out::OutputBIN, field::Field, analyzer::FieldAnalyzer, z)
    if z >= out.zout
        dset = @sprintf("%03d", out.iout)
        @printf("Writing dataset %s...\n", dset)

        fp = HDF5.h5open(out.fname, "r+")
        group = fp["field"]
        write_field(field, group, dset)
        HDF5.attributes(group[dset])["z"] = z
        HDF5.close(fp)

        out.iout = out.iout + 1
        out.zout = out.zout + out.dzout
    end

    if out.zdata
        fp = HDF5.h5open(out.fname, "r+")
        group = fp["zdata"]
        write_zdata(analyzer, group, out.izdata)
        HDF5.close(fp)

        out.izdata = out.izdata + 1
    end
    return nothing
end

# ******************************************************************************
# TXT + BIN
# ******************************************************************************
struct Output{TT<:OutputTXT, TB<:OutputBIN}
    txt :: TT
    bin :: TB
end

function Output(
    grid::Grid, field::Field, neu, zu, z, dzout;
    prefix::String="",
)
    if dirname(prefix) != ""
        mkpath(dirname(prefix))
    end

    fname = joinpath(dirname(prefix), basename(prefix) * "out.dat")
    txt = OutputTXT(fname, grid, field, neu, zu)

    fname = joinpath(dirname(prefix), basename(prefix) * "out.h5")
    bin = OutputBIN(fname, grid, zu, z, dzout)

    return Output(txt, bin)
end

function write_output(out::Output, field::Field, analyzer::FieldAnalyzer, z)
    write_txt(out.txt, analyzer)
    write_bin(out.bin, field, analyzer, z)
    return nothing
end


text:
################################################################################
### Linalg additions
################################################################################

# Notes:
# * det(A B) = det(A) det(B)
# * det(A^-1) = 1/det(A)
# * |det(U)| = 1 (complex norm, U unitary)
# * det(T) = 1 (T unit-triangular like our T's)
# * our UDT decomposition always makes D real positive
# * all weights (det(G)) should be real positive
# * local updates already check ratios (det(G)/det(G')) so it's probably OK to
#   ignore phases here!?
# * With the introduction of fields we now pass it through slice matrix
#   functions because we dispatch on it down the line. (interaction_matrix_exp!)
#   These functions now assume the conf that matters to be field.conf, so we
#   need to update that instead of a temp_conf

# TODO WARNING
# Global updates assume that U and T only ever contribute a complex phase. I did
# not verify this thoroughly. Furthermore it assumes that the probability
# always simplifies to ||p||.

@bm function calculate_inv_greens_udt(Ul, Dl, Tl, Ur, Dr, Tr, G, pivot, temp)
    # G = [I + Ul Dl Tl Tr^† Dr Ur^†]^-1
    vmul!(G, Tl, adjoint(Tr))
    vmul!(Tr, G, Diagonal(Dr))
    vmul!(G, Diagonal(Dl), Tr)
    #   = [I + Ul G Ur^†]^-1
    udt_AVX_pivot!(Tr, Dr, G, pivot, temp, Val(false)) # Dl available
    #   = [I + Ul Tr Dr G Ur^†]^-1  w/ Tr unitary, G triangular
    #   = Ur G^-1 [(Ul Tr)^† Ur G^-1 + Dr]^-1 (Ul Tr)^†
    vmul!(Tl, Ul, Tr)
    rdivp!(Ur, G, Ul, pivot) # requires unpivoted udt decompostion (Val(false))
    #   = Ur [Tl^† Ur + Dr]^-1 Tl^†  w/ Tl unitary, Ur not
    vmul!(Tr, adjoint(Tl), Ur)
    rvadd!(Tr, Diagonal(Dr))
    #   = Ur Tr^-1 Tl^†
    udt_AVX_pivot!(Ul, Dr, Tr, pivot, temp, Val(false)) # Dl available
    #   = Ur Tr^-1 Dr^-1 Ul^† Tl^†
    #   = (old_Ur / G) Tr^-1 Dr^-1 Ul^† Tl^†
    # with old_Ur, Ul, Tl unitary and G, Tr triangular
    # det(G) = phase1 / 1 / 1 / det(Dr) / phase2 / phase3
    # where we ignore phases because they should be 1 and we already check this
    # in local updates.
    return
end

# after the above without modifying Ur, Tr, Tl, Ul, Dr
@bm function finish_calculate_greens(Ul, Dl, Tl, Ur, Dr, Tr, G, pivot, temp)
    # G = Ur Tr^-1 Dr^-1 Ul^† Tl^†
    rdivp!(Ur, Tr, G, pivot) # requires unpivoted udt decompostion (false)
    vmul!(Tr, Tl, Ul)
    #   = Ur Dr^-1 Tr^†
    vinv!(Dl, Dr)
    vmul!(Ul, Ur, Diagonal(Dl))
    vmul!(G, Ul, adjoint(Tr))
    return G
end

# This calculates the UDT stack stuff from scratch, but doesn't calculate greens
# fully. We use that det(UDT) = prod(D), i.e. that det(U) = 1 by definition and
# det(T) = 1 because T is unit-triangular by construction
@bm function inv_det(
        mc::DQMC, slice::Int, field::AbstractField,
        safe_mult::Int = mc.parameters.safe_mult
    )
    copyto!(mc.stack.curr_U, I)
    copyto!(mc.stack.Ur, I)
    mc.stack.Dr .= one(eltype(mc.stack.Dr))
    copyto!(mc.stack.Tr, I)

    # Calculate Ur,Dr,Tr=B(slice)' ... B(M)'
    if slice+1 <= mc.parameters.slices
        start = slice+1
        stop = mc.parameters.slices
        for k in reverse(start:stop)
            if mod(k,safe_mult) == 0
                multiply_daggered_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, field)
                vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dr))
                udt_AVX_pivot!(mc.stack.curr_U, mc.stack.Dr, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
                copyto!(mc.stack.tmp2, mc.stack.Tr)
                vmul!(mc.stack.Tr, mc.stack.tmp1, mc.stack.tmp2)
            else
                multiply_daggered_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, field)
            end
        end
        vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dr))
        udt_AVX_pivot!(mc.stack.Ur, mc.stack.Dr, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
        copyto!(mc.stack.tmp2, mc.stack.Tr)
        vmul!(mc.stack.Tr, mc.stack.tmp1, mc.stack.tmp2)
    end

    copyto!(mc.stack.curr_U, I)
    copyto!(mc.stack.Ul, I)
    mc.stack.Dl .= one(eltype(mc.stack.Dl))
    copyto!(mc.stack.Tl, I)

    # Calculate Ul,Dl,Tl=B(slice-1) ... B(1)
    if slice >= 1
        start = 1
        stop = slice
        for k in start:stop
            if mod(k,safe_mult) == 0
                multiply_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, field)
                vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dl))
                udt_AVX_pivot!(mc.stack.curr_U, mc.stack.Dl, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
                copyto!(mc.stack.tmp2, mc.stack.Tl)
                vmul!(mc.stack.Tl, mc.stack.tmp1, mc.stack.tmp2)
            else
                multiply_slice_matrix_left!(mc, mc.model, k, mc.stack.curr_U, field)
            end
        end
        vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(mc.stack.Dl))
        udt_AVX_pivot!(mc.stack.Ul, mc.stack.Dl, mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv)
        copyto!(mc.stack.tmp2, mc.stack.Tl)
        vmul!(mc.stack.Tl, mc.stack.tmp1, mc.stack.tmp2)
    end

    calculate_inv_greens_udt(
        mc.stack.Ul, mc.stack.Dl, mc.stack.Tl, mc.stack.Ur, mc.stack.Dr, mc.stack.Tr,
        mc.stack.greens_temp, mc.stack.pivot, mc.stack.tempv
    )
    return mc.stack.Dr
end

################################################################################
### Global Metropolis update
################################################################################

@bm function propose_global_from_conf(mc::DQMC, m::Model, field::AbstractField)
    # I don't think we need this...
    @assert mc.stack.current_slice == 1
    @assert mc.stack.direction == 1

    # This should be just after calculating greens, so mc.stack.Dl is from the UDT
    # decomposed G
    # We need an independent temp vector here as inv_det changes Dl, Dr and tempv
    mc.stack.tempvf .= mc.stack.Dl

    # This is essentially reverse_build_stack + partial calculate_greens
    # after this: G = Ur Tr^-1 Dr^-1 Ul^† Tl^†
    # where Ur, Ul, Tl only contribute complex phases and Tr contributes 1
    # Since we are already checking for sign problems in local updates we ignore
    # the phases here and set det(G) = 1 / det(Dr)
    inv_det(mc, current_slice(mc)-1, field)

    # This loop helps with stability - it multiplies large and small numbers
    # whihc avoid reaching extremely large or small (typemin/max) floats
    detratio = 1.0
    for i in eachindex(mc.stack.tempvf)
        detratio *= mc.stack.tempvf[i] * mc.stack.Dr[i]
    end
    ΔE_Boson = energy_boson(field, conf(field)) - energy_boson(field, temp_conf(field))

    if max(nflavors(m), nflavors(field)) == 1
        detratio = detratio * detratio
    end

    return detratio, ΔE_Boson, nothing
end

@bm function accept_global!(mc::DQMC, m::Model, field, passthrough)
    # for checking
    # new_G = finish_calculate_greens(
    #     mc.stack.Ul, mc.stack.Dl, mc.stack.Tl, mc.stack.Ur, mc.stack.Dr, mc.stack.Tr,
    #     mc.stack.greens_temp, mc.stack.pivot, mc.stack.tempv
    # )

    # Need a full stack rebuild
    reverse_build_stack(mc, mc.stack)
    # This calculates greens
    propagate(mc)

    # @info mc.stack.current_slice, mc.stack.direction
    # which should match new_G
    # display(new_G .- mc.stack.greens)
    # @assert new_G ≈ mc.stack.greens
    nothing
end

# This does a MC update with the given temp_conf as the proposed new_conf
function global_update(mc::DQMC, model::Model, field::AbstractField)
    detratio, ΔE_boson, passthrough = propose_global_from_conf(mc, model, field)

    p = abs(exp(- ΔE_boson) * detratio)

    # Gibbs/Heat bath
    # p = p / (1.0 + p)
    # Metropolis
    if p > 1 || rand() < p
        accept_global!(mc, model, field, passthrough)
        return 1
    else
        copyto!(field.conf, field.temp_conf)
    end

    return 0
end

################################################################################
### Update Interface
################################################################################

abstract type AbstractGlobalUpdate <: AbstractUpdate end
requires_temp_conf(update::AbstractGlobalUpdate) = true
@inline function update(u::AbstractGlobalUpdate, mc, model, field)
    propose_conf!(u, mc, model, field)
    return global_update(mc, model, field)
end

"""
    GlobalFlip([mc, model])

A global update that flips the configuration (±1 -> ∓1).
"""
struct GlobalFlip <: AbstractGlobalUpdate end
GlobalFlip(mc, model) = GlobalFlip()
name(::GlobalFlip) = "GlobalFlip"

@bm function propose_conf!(::GlobalFlip, mc, model, field::AbstractHirschField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    c .*= -1
    return nothing
end
@bm function propose_conf!(::GlobalFlip, mc, model, field::AbstractGHQField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    c .= Int8(5) .- c
    return nothing
end

"""
    GlobalShuffle([mc, model])

A global update that shuffles the current configuration. Note that this is not
local to a time slice.
"""
struct GlobalShuffle <: AbstractGlobalUpdate end
GlobalShuffle(mc, model) = GlobalShuffle()
name(::GlobalShuffle) = "GlobalShuffle"

@bm function propose_conf!(::GlobalShuffle, mc, model, field::AbstractField)
    copyto!(temp_conf(field), conf(field))
    shuffle!(conf(field))
    return nothing
end

"""
    SpatialShuffle([mc, model])

A global update that randomly swaps spatial indices of a configuration without
changing the temporal indices. I.e. it may set
`new_conf[i, :] = current_conf[j, :]` where the second index is a time slice
index.
"""
struct SpatialShuffle <: AbstractGlobalUpdate
    indices::Vector{Int}
    SpatialShuffle() = new(Vector{Int}(undef, 0))
end
SpatialShuffle(mc, model) = SpatialShuffle()
function init!(mc, u::SpatialShuffle)
    resize!(u.indices, length(lattice(mc)))
    u.indices .= 1:length(lattice(mc))
    nothing
end
name(::SpatialShuffle) = "SpatialShuffle"

@bm function propose_conf!(u::SpatialShuffle, mc, model, field::AbstractHirschField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    shuffle!(u.indices)
    for slice in 1:nslices(mc), (i, j) in enumerate(u.indices)
        c[i, slice] = tc[j, slice]
    end
    return nothing
end

"""
    TemporalShuffle([mc, model])

A global update that randomly swaps time slice indices of a configuration
without changing the spatial indices. I.e. it may set
`new_conf[:, k] = current_conf[:, l]` where the second index is a time slice
index.
"""
struct TemporalShuffle <: AbstractGlobalUpdate
    indices::Vector{Int}
    TemporalShuffle() = new(Vector{Int}(undef, 0))
end
TemporalShuffle(mc, model) = TemporalShuffle()
function init!(mc, u::TemporalShuffle)
    resize!(u.indices, nslices(mc))
    u.indices .= 1:nslices(mc)
    nothing
end
name(::TemporalShuffle) = "TemporalShuffle"

@bm function propose_conf!(u::TemporalShuffle, mc, model, field::AbstractHirschField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    shuffle!(u.indices)
    for (k, l) in enumerate(u.indices), i in 1:length(lattice(mc))
        c[i, k] = tc[i, l]
    end
    return nothing
end

"""
    Denoise([mc, model])

This global update attempts to remove noise, i.e. it attempts to build spatial
domains. This is done by setting each site to dominant surrounding value.
"""
struct Denoise <: AbstractGlobalUpdate end
Denoise(mc, model) = Denoise()
name(::Denoise) = "Denoise"

@bm function propose_conf!(::Denoise, mc, model, field::AbstractHirschField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    for slice in 1:nslices(mc), i in 1:length(lattice(mc))
        average = tc[i, slice]
        for j in neighbors(lattice(model), i)
            average += 2 * tc[j, slice]
        end
        c[i, slice] = sign(average)
    end
    return nothing
end

"""
    DenoiseFlip([mc, model])

This update is similar to `Denoise` but sets each site to the opposite value to
its surrounding.
"""
struct DenoiseFlip <: AbstractGlobalUpdate end
DenoiseFlip(mc, model) = DenoiseFlip()
name(::DenoiseFlip) = "DenoiseFlip"

@bm function propose_conf!(::DenoiseFlip, mc, model, field::AbstractHirschField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    for slice in 1:nslices(mc), i in 1:length(lattice(mc))
        average = tc[i, slice]
        for j in neighbors(lattice(model), i)
            average += 2 * tc[j, slice]
        end
        c[i, slice] = -sign(average)
    end
    return nothing
end

"""
    StaggeredDenoise([mc, model])

This update is similar to `Denoise` but adds a multiplier based on lattice site
index. Even sites get multiplied by `+1`, odd sites by `-1`.
"""
struct StaggeredDenoise <: AbstractGlobalUpdate end
StaggeredDenoise(mc, model) = StaggeredDenoise()
name(::StaggeredDenoise) = "StaggeredDenoise"

@bm function propose_conf!(::StaggeredDenoise, mc, model, field::AbstractHirschField)
    c = conf(field); tc = temp_conf(field)
    copyto!(tc, c)
    for slice in 1:nslices(mc), i in 1:length(lattice(mc))
        average = tc[i, slice]
        for j in neighbors(lattice(model), i)
            average += 2 * tc[j, slice]
        end
        c[i, slice] = (1 - 2 * (i % 2)) * sign(average)
    end
    return nothing
end


text:
using ColorTypes, FixedPointNumbers
using Test

@testset "single color" begin
    iob = IOBuffer()
    cf  = RGB{Float32}(0.32218,0.14983,0.87819)
    c   = convert(RGB{N0f8}, cf)
    c16 = RGB{N0f16}(0.32218,0.14983,0.87819)
    ca  = RGBA{N0f8}(0.32218,0.14983,0.87819,0.99241)
    ac  = alphacolor(ca)

    show(iob, c)
    @test String(take!(iob)) == "RGB{N0f8}(0.322,0.149,0.878)"
    show(iob, c16)
    @test String(take!(iob)) == "RGB{N0f16}(0.32218,0.14983,0.87819)"
    show(iob, cf)
    @test String(take!(iob)) == "RGB{Float32}(0.32218f0,0.14983f0,0.87819f0)"
    show(IOContext(iob, :compact => true), cf)
    @test String(take!(iob)) == "RGB{Float32}(0.32218,0.14983,0.87819)"
    show(iob, ca)
    @test String(take!(iob)) == "RGBA{N0f8}(0.322,0.149,0.878,0.992)"
    show(IOContext(iob, :compact => true), ca)
    @test String(take!(iob)) == "RGBA{N0f8}(0.322,0.149,0.878,0.992)"
    show(iob, ac)
    @test String(take!(iob)) == "ARGB{N0f8}(0.322,0.149,0.878,0.992)"
    show(IOContext(iob, :compact => true), ac)
    @test String(take!(iob)) == "ARGB{N0f8}(0.322,0.149,0.878,0.992)"

    show(iob, RGB24(0.4,0.2,0.8))
    @test String(take!(iob)) == "RGB24(0.4N0f8,0.2N0f8,0.8N0f8)"
    show(iob, ARGB32(0.4,0.2,0.8,1.0))
    @test String(take!(iob)) == "ARGB32(0.4N0f8,0.2N0f8,0.8N0f8,1.0N0f8)"

    show(iob, Gray(0.8))
    @test String(take!(iob)) == "Gray{Float64}(0.8)"
    show(iob, GrayA(0.8))
    @test String(take!(iob)) == "GrayA{Float64}(0.8,1.0)"
    show(iob, AGray(0.8))
    @test String(take!(iob)) == "AGray{Float64}(0.8,1.0)"
    show(iob, Gray24(0.4))
    @test String(take!(iob)) == "Gray24(0.4N0f8)"
    show(iob, AGray32(0.8))
    @test String(take!(iob)) == "AGray32{N0f8}(0.8,1.0)"
end

@testset "collection of colors" begin
    iob = IOBuffer()
    summary(iob, Gray{N0f8}[0.2, 0.4, 0.6])
    @test String(take!(iob)) == "3-element Array{Gray{N0f8},1} with eltype Gray{Normed{UInt8,8}}"
    mat = RGB{Float64}[RGB(1,0,0) RGB(0,1,0); RGB(0,0,1) RGB(1,1,1)]
    summary(iob, mat)
    @test String(take!(iob)) == "2×2 Array{RGB{Float64},2} with eltype RGB{Float64}"
    summary(iob, view(mat,:,:))
    @test String(take!(iob)) == "2×2 view(::Array{RGB{Float64},2}, :, :) with eltype RGB{Float64}"
end


text:
module SWI

using MRI, Statistics

include("tissue.jl")
include("utility.jl")
include("functions.jl")

export calculateSWI, createMIP, saveconfiguration, Data, Options

"""
    SWI.dir(path...)
Construct a path relative to SWI root.
# Example
```julia
julia> SWI.dir("test","testData","small","Mag.nii")
"/home/korbinian90/.julia/dev/SWI/test/testData/small/Mag.nii"
```
"""
dir(path...) = joinpath(dirname(@__DIR__),path...)

end


text:
@resumable function customerLCFS(env::Environment, server::Resource, id::Integer, t_a::Float64, d_s::Distribution, starttime::Vector{Float64}, departtime::Vector{Float64}, NA::Vector{Int}, ND::Vector{Int}, NC::Vector{Int})
    @yield timeout(env, t_a) # customer arrives
    NA[id] = NC[1]
    NC[1] += 1
    @yield request(server, priority = -id) # customer starts service
    starttime[id] = now(env)
    @yield timeout(env, rand(d_s)) # server is busy
    @yield release(server) # customer exits service
    departtime[id] = now(env)
    NC[1] -= 1
    ND[id] = NC[1]
end

function runsimLCFS(s::Sim, q::Queue)
    ncust = s.ncust
    seed = s.seed
    c = q.c
    adist = q.adist
    sdist = q.sdist
    srand(seed) # set rng seed
    # output of queue simulation
    # arrive order, arrive time, start service time, depart order, depart time
    arriveorder = collect(1:ncust)
    arrivetime = cumsum(rand(adist, ncust))
    starttime = zeros(ncust)
    departtime = zeros(ncust)
    NA = zeros(Int, ncust)
    ND = zeros(Int, ncust)
    NC = zeros(Int, 1)
    # setup and run simulation
    queuesim = Simulation() # initialize simulation environment
    server = Resource(queuesim, c) # initialize servers
    for i = 1:ncust # initialize customers
        @process customerLCFS(queuesim, server, i, arrivetime[i], sdist, starttime, departtime, NA, ND, NC)
    end
    run(queuesim) # run simulation
    depart = hcat(sortrows(hcat(departtime, arriveorder)), collect(1:ncust))
    departorder = convert(Vector{Int64}, sortrows(depart, by = x -> (x[2],x[1],x[3]))[:,3])
    DataFrame(aorder = arriveorder, dorder = departorder, atime = arrivetime, stime = starttime, dtime = departtime, na = NA, nd = ND)
end

function lcfs(df::DataFrame, dtype::Symbol)
    arrivals = DataFrame(aorder = df[:aorder], time = df[:atime])
    departs = DataFrame(daorder = df[:aorder], time = df[dtype])
    dfl = join(arrivals, departs, on = :time, kind = :outer)
    sort!(dfl, :time)
    len = size(dfl,1)
    lcfsorder = zeros(Int, len)
    for i in 1:len
        a = dfl[:aorder][i]
        lcfsorder[i] = ismissing(a) ? dfl[:daorder][i] : a
    end
    lcfsorder, dfl
end

# LCFS Disorder Function
# Create noisy data for robust methods
# Input: queue sim output
# Output: true and noisy arrival and departure data
function disorderLCFS(simout::DataFrame)
  len = size(simout, 1)
  D_meas = zeros(len) # measure error in depart times (original order)
  lcfsorder, dfl = lcfs(simout, :dtime)
  len2 = size(dfl,1)
  ismissing(dfl[:daorder][len2]) ? nothing : D_meas[dfl[:daorder][len2]] = dfl[:time][len2]
  for i = 2:(len2 - 1)
      ind = dfl[:daorder][i]
      if !ismissing(ind)
          a = mean(dfl[:time][(i-1):i])
          b = mean(dfl[:time][i:(i+1)])
          D_meas[ind] = a + rand() * (b - a)
      end
  end
  simout[:dmeas] = D_meas
  simout
end

function firstdepart(df::DataFrame, n::Int, t::Float64)
    dfdepart = sort(df, :dtime)
    indt = findfirst(dfdepart[:dtime] .> t, true)
    indt == 0 ? indt = size(dfdepart[:dtime],1) : nothing
    indn = findfirst(dfdepart[:nd][indt:end], n)
    indn == 0 ? indn = size(dfdepart[:nd][indt:end],1) : nothing
    dfdepart[:dtime][indt + indn - 1]
end

function c_order_LCFS(lcfsorder::Vector{Int64})
    nservers = 0
    N = size(lcfsorder, 1)
    est = zeros(Int, N)
    system = Set{Int}()
    service = Set{Int}()
    for i = 1:N
        if lcfsorder[i] in system
            setdiff!(system, lcfsorder[i])
            lcfsorder[i] in service && setdiff!(service, lcfsorder[i])
            if !isempty(system)
                notinservice = setdiff(system, service)
                !isempty(notinservice) && push!(service, maximum(notinservice))
            end
        else
            push!(system, lcfsorder[i])
        end
        est[i] = length(service)
    end
    maximum(est)
end

# Summary Estimator
# Estimate the number of servers using the Variance and Uninformed methods
# Input: arrival times, departure times, and the max number of servers
# Output: estimated number of servers for each method
function c_var_unf_LCFS(df::DataFrame, A::Vector{Float64}, D::Vector{Float64}, cmax::Int64)
  # Initilaize data structures and functions
  N = size(D, 1)
  B̂ = zeros(N)
  temp = zeros(N)
  VS = zeros(cmax, 2) # var, unf
  max_undef = 0
  NA = df[:na]
  ND = df[:nd]
  ## Grid Search
  for c = 1:cmax
    for i = 1:N
        if NA[i] < c
            B̂[i] = A[i]
        else
            B̂[i] = firstdepart(df, NA[i], A[i])
        end
    end
    Ŝ = D - B̂
    if any(Ŝ .< 0)
      VS[c, 1:2] = -1
      max_undef = c
    else
      VS[c, 1] = var(Ŝ) * (N-1) # var
    end
  end
  (r1, r2) = (99, 99) # initialize result to 99 for debugging
  # var method
  try
    r1 = findin(VS[:, 1], minimum(filter(x -> x >= 0, VS[:, 1])))[1]
  catch
    r1 = 1
  end
  # uninformed method
  r2 = max_undef + 1
  (r1, r2, VS)
end

# Convergence Estimator
# Calculate aprpoximate convergence
# Input: parametrs, queueing output data, number of servers
# Output: number of observations to approximately converge
function convergenceLCFS(param_inf::Paraminf, indrun::Int)
    settings = param_inf.settings
    aname = settings[indrun, 2] # arrival distribution
    sname = settings[indrun, 3] # service distribution
    c_true = settings[indrun, 5] # actual number of servers
    rho = settings[indrun, 6] # traffic intensity
    max_servers = param_inf.max_servers
    lim = param_inf.max_servers + 1
    n_methods = param_inf.n_methods
    obs_max = param_inf.obs_max
    time_limit = param_inf.time_limit
    window = param_inf.window
    step = param_inf.step
    ϵ = param_inf.ϵ
    bool_true = falses(n_methods)
    bool_meas = falses(n_methods)
    conv_true = zeros(n_methods)
    conv_meas = zeros(n_methods)
    array_size = floor(Int64, (obs_max - step) / step)
    ests_true = zeros(Int64, array_size, n_methods) # c estimates by obs by run
    ests_meas = zeros(Int64, array_size, n_methods) # c estimates by obs by run
    (adist, sdist) = builddist(c_true, aname, sname, rho)
    q = Queue(c_true, adist, sdist)
    seed = param_inf.seed
    for (ind, obs) in enumerate(lim:step:(obs_max - step))
        s = Sim(obs, time_limit, seed)
        q_out = runsimLCFS(s,q)
        output = disorderLCFS(q_out) # queue sim results with noise
        # true departure
        lcfsorder, dfl = lcfs(output, :dtime)
        rol = c_order_LCFS(lcfsorder)
        A1, D1 = output[:atime], output[:dtime]
        (rvl, rul, VSl) = c_var_unf_LCFS(output, A1, D1, max_servers)
        # time error
        lcfsorder, dfl = lcfs(output, :dmeas)
        rol_meas = c_order_LCFS(lcfsorder)
        A1, D1 = output[:atime], output[:dmeas]
        (rvl_meas, rul_meas, VSl_meas) = c_var_unf_LCFS(output, A1, D1, max_servers)
        # record results
        ests_true[ind, :] = [rol rvl rul]
        ests_meas[ind, :] = [rol_meas rvl_meas rul_meas]
        if ind >= window
        for j = 1:n_methods
          if !bool_true[j] && mean(abs.(ests_true[(ind - window + 1):(ind), j] - c_true)) < ϵ
            conv_true[j] = (ind - window + 1) * step
            bool_true[j] = true
          end
          if !bool_meas[j] && mean(abs.(ests_meas[(ind - window + 1):(ind), j]- c_true)) < ϵ
            conv_meas[j] = (ind - window + 1) * step
            bool_meas[j] = true
          end
        end
        end
    end
    (conv_true, conv_meas, ests_true, ests_meas)
end

function esterrorLCFS(param_inf::Paraminf, indrun::Int)
    settings = param_inf.settings
    aname = settings[indrun, 2] # arrival distribution
    sname = settings[indrun, 3] # service distribution
    obs_limit = settings[indrun, 4] # number of runs
    c_true = settings[indrun, 5] # actual number of servers
    rho = settings[indrun, 6] # traffic intensity
    window_detail = param_inf.window_detail
    n_methods = param_inf.n_methods
    max_servers = param_inf.max_servers
    time_limit = param_inf.time_limit
    ind_l = Int64(obs_limit)
    ind_u = Int64(obs_limit  + window_detail - 1)
    detail_true = zeros(Int64, window_detail, n_methods) - 1 # precise c estimates
    detail_meas = zeros(Int64, window_detail, n_methods) - 1 # precise c estimates
    error_true = zeros(n_methods)
    error_meas = zeros(n_methods)
    (adist, sdist) = builddist(c_true, aname, sname, rho)
    q = Queue(c_true, adist, sdist)
    seed = param_inf.seed
    for (ind, obs) in enumerate(ind_l:ind_u)
        s = Sim(obs, time_limit, seed)
        q_out = runsimLCFS(s,q)
        output = disorderLCFS(q_out) # queue sim results with noise
        # true departure
        lcfsorder, dfl = lcfs(output, :dtime)
        rol = c_order_LCFS(lcfsorder)
        A1, D1 = output[:atime], output[:dtime]
        (rvl, rul, VSl) = c_var_unf_LCFS(output, A1, D1, max_servers)
        # time error
        lcfsorder, dfl = lcfs(output, :dmeas)
        rol_meas = c_order_LCFS(lcfsorder)
        A1, D1 = output[:atime], output[:dmeas]
        (rvl_meas, rul_meas, VSl_meas) = c_var_unf_LCFS(output, A1, D1, max_servers)
        # record results
        detail_true[ind, :] = [rol rvl rul]
        detail_meas[ind, :] = [rol_meas rvl_meas rul_meas]
    end
    for j = 1:n_methods
        error_true[j] = mean(abs.(detail_true[:, j] - c_true))
        error_meas[j] = mean(abs.(detail_meas[:, j] - c_true))
    end
    (error_true, error_meas, detail_true, detail_meas)
end

function inferLCFS(param_inf::Paraminf)
    settings = param_inf.settings
    n_runs = param_inf.n_runs
    n_methods = param_inf.n_methods
    max_servers = param_inf.max_servers
    obs_max = param_inf.obs_max
    time_limit = param_inf.time_limit
    ϵ = param_inf.ϵ
    window = param_inf.window
    step = param_inf.step
    seed = param_inf.seed
    srand(seed)
    raw_true = Array{Array{Int64}}(n_runs, n_methods) # c estimates by obs
    raw_meas = Array{Array{Int64}}(n_runs, n_methods) # c estimates by obs
    conv_true = zeros(n_runs, n_methods) + obs_max - window * step # num obs to conv
    conv_meas = zeros(n_runs, n_methods) + obs_max - window * step # num obs to conv
    error_true = zeros(n_runs, n_methods) # c estimate with true departure times
    error_meas = zeros(n_runs, n_methods) # c estimate with noisy departure times
    # caculate each metric for each DOE run
    for i = 1:n_runs
        print("Run $i: ") # status update
        array_size = floor(Int64, (obs_max - 20) / step)
        ests_true = zeros(Int64, array_size, n_methods) # c estimates by obs by run
        ests_meas = zeros(Int64, array_size, n_methods) # c estimates by obs by run
        ests_err = zeros(Int64, array_size, n_methods) # c estimates by obs by run
        print("conv, ")
        (conv_true[i,:], conv_meas[i,:], ests_single, ests_meas_single) = convergenceLCFS(param_inf, i)
        println("err")
        (error_true[i,:], error_meas[i,:], detail_single, detail_meas_single) = esterrorLCFS(param_inf, i)
        for j = 1:n_methods
            raw_true[i,j] = ests_single[:,j]
            raw_meas[i,j] = ests_meas_single[:,j]
        end
  end
  resp_error = (error_true, error_meas)
  resp_conv = (conv_true, conv_meas)
  resp_raw = (raw_true, raw_meas)
  (resp_error, resp_conv, resp_raw)
end


text:
# Autogenerated wrapper script for SOLAR_jll for armv7l-linux-gnueabihf-cxx03
export solar

JLLWrappers.@generate_wrapper_header("SOLAR")
JLLWrappers.@declare_executable_product(solar)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        solar,
        "bin/solar",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Compare Jain-Neal sampler to Reversible Jump on simulated data of increasing dimension.
module SimCompareJNtoRJ

using Statistics
using Random
using BayesianMixtures
B = BayesianMixtures

# Settings
ds = [4,8,12]  # increasing dimensions
reset_random = true  # reset the random number generator before each run
save_figures = true  # save the figures to file
n_total = 1000000  # number of MCMC sweeps to run each algorithm

fignum = 0
for (i_d,d) in enumerate(ds)
    B.open_figure(9+i_d)
    for (label,mode,color) in [("Collapsed Jain-Neal","MVNaaC","g"),
                               ("Uncollapsed Jain-Neal","MVNaaN","r"),
                               ("RJMCMC","MVNaaRJ","b")]

        if reset_random; Random.seed!(0); end

        # Generate data
        n = 100
        xs = [randn(d) for i = 1:n]
        zs = [Int(ceil(rand()*3)) for i = 1:n]
        shift = (3/sqrt(d))*[-1,0,1]
        x = convert(Array{Array{Float64,1},1},[xs[i].+shift[zs[i]] for i = 1:n])
        mu = mean(x)  # sample mean
        v = mean([xi.*xi for xi in x]) - mu.*mu  # sample variance
        x = [((xi-mu)./sqrt.(v))::Array{Float64,1} for xi in x] # normalized to zero mean, unit variance

        # Run sampler
        println(label)
        options = B.options(mode,"MFM",x,n_total; n_keep=1000,t_max=20)
        result = B.run_sampler(options)

        # Traceplot of t
		global fignum
        B.open_figure(fignum+=1)
        t_show = 20  # number of seconds to show
        B.traceplot_timewise(result,t_show)
        B.PyPlot.title("$label, d = $d")
        if save_figures; B.PyPlot.savefig("simCompareJNtoRJ-traceplot-$mode-d=$d.png",dpi=200); end

        # Plot autocorrelation
        B.open_figure(9+i_d; clear_figure=false)
        t_show = 0.025*d  # number of seconds to show
        B.plot_autocorrelation(result,t_show; color=color,label=label)
        B.PyPlot.title("Autocorrelation, d = $d")
    end
    if i_d==1; B.PyPlot.legend(loc="upper right",fontsize=12); end
    if save_figures; B.PyPlot.savefig("simCompareJNtoRJ-autocorrelation-d=$d.png",dpi=200); end
end

end # module SimCompareJNtoRJ