text:
mutable struct DebuggerState
    stack::Vector{JuliaStackFrame}
    level::Int
    repl
    terminal
    main_mode
    julia_prompt::Ref{LineEdit.Prompt}
    standard_keymap
    overall_result
end
DebuggerState(stack, repl, terminal) = DebuggerState(stack, 1, repl, terminal, nothing, Ref{LineEdit.Prompt}(), nothing, nothing)
DebuggerState(stack, repl) = DebuggerState(stack, repl, nothing)

function sparam_syms(meth::Method)
    s = Symbol[]
    sig = meth.sig
    while sig isa UnionAll
        push!(s, Symbol(sig.var.name))
        sig = sig.body
    end
    return s
end

function print_backtrace(state::DebuggerState)
    for (num, frame) in enumerate(state.stack)
        print_frame(Base.pipe_writer(state.terminal), num, frame)
    end
end

print_backtrace(state::DebuggerState, _::Nothing) = nothing

function execute_command(state::DebuggerState, frame, ::Val{:bt}, cmd)
    print_backtrace(state)
    println()
    return false
end

function execute_command(state::DebuggerState, frame, _, cmd)
    println("Unknown command `$cmd`. Executing `?` to obtain help.")
    execute_command(state, frame, Val{Symbol("?")}(), "?")
end

function execute_command(state::DebuggerState, _::JuliaStackFrame, ::Union{Val{:f},Val{:fr}}, cmd)
    subcmds = split(cmd,' ')[2:end]
    if isempty(subcmds) || subcmds[1] == "v"
        print_frame(Base.pipe_writer(state.terminal), state.level, state.stack[state.level])
        return false
    else
        new_level = parse(Int, subcmds[1])
        new_stack_idx = length(state.stack)-(new_level-1)
        if new_stack_idx > length(state.stack) || new_stack_idx < 1
            printstyled(stderr, "Not a valid frame index\n"; color=:red)
            return false
        end
        state.level = new_level
    end
    return true
end

"""
    Start debugging the specified code in the the specified environment.
    The second argument should default to the global environment if such
    an environment exists for the language in question.
"""
function debug(meth::Method, args...)
    stack = [enter_call(meth, args...)]
    RunDebugger(stack)
end

struct FileLocInfo
    filepath::String
    line::Int
    # 0 if unknown
    column::Int
    # The line at which the current context starts, 0 if unknown
    defline::Int
end

struct BufferLocInfo
    data::String
    line::Int
    # 0 if unknown
    column::Int
    defline::Int
end

function loc_for_fname(file::Symbol, line::Integer, defline::Integer)
    if startswith(string(file),"REPL[")
        hist_idx = parse(Int,string(file)[6:end-1])
        isdefined(Base, :active_repl) || return nothing, ""
        hp = Base.active_repl.interface.modes[1].hist
        return BufferLocInfo(hp.history[hp.start_idx+hist_idx], line, 0, defline)
    else
        for path in SEARCH_PATH
            fullpath = joinpath(path,string(file))
            if isfile(fullpath)
                return FileLocInfo(fullpath, line, 0, defline)
            end
        end
    end
    return nothing
end

function locinfo(frame::JuliaStackFrame)
    if frame.code.scope isa Method
        meth = frame.code.scope
        loc_for_fname(meth.file, linenumber(frame), meth.line)
    else
        println("not yet implemented")
    end
end

function locdesc(frame::JuliaStackFrame)
    sprint() do io
        if frame.code.scope isa Method
            meth = frame.code.scope
            argnames = frame.code.code.slotnames[2:meth.nargs]
            spectypes = Any[Any for i=1:length(argnames)]
            print(io, meth.name,'(')
            first = true
            for (argname, argT) in zip(argnames, spectypes)
                first || print(io, ", ")
                first = false
                print(io, argname)
                !(argT === Any) && print(io, "::", argT)
            end
            print(io, ") at ",
                frame.code.fullpath ? meth.file :
                basename(String(meth.file)),
                ":",meth.line)
        else
            println("not yet implemented")
        end
    end
end

"""
Determine the offsets in the source code to print, based on the offset of the
currently highlighted part of the code, and the start and stop line of the
entire function.
"""
function compute_source_offsets(code::String, offset::Integer, startline::Integer, stopline::Integer; file::SourceFile = SourceFile(code))
    offsetline = compute_line(file, offset)
    if offsetline - 3 > length(file.offsets) || startline > length(file.offsets)
        return -1, -1
    end
    startoffset = max(file.offsets[max(offsetline-3,1)], file.offsets[startline])
    stopoffset = lastindex(code)-1
    if offsetline + 3 < lastindex(file.offsets)
        stopoffset = min(stopoffset, file.offsets[offsetline + 3]-1)
    end
    if stopline + 1 < lastindex(file.offsets)
        stopoffset = min(stopoffset, file.offsets[stopline + 1]-1)
    end
    startoffset, stopoffset
end

function print_sourcecode(io::IO, code::String, line::Integer, defline::Integer; file::SourceFile = SourceFile(code))
    startoffset, stopoffset = compute_source_offsets(code, file.offsets[line], defline, line+3; file=file)

    if startoffset == -1
        printstyled(io, "Line out of file range (bad debug info?)", color=:bold)
        return
    end

    # Compute necessary data for line numbering
    startline = compute_line(file, startoffset)
    stopline = compute_line(file, stopoffset)
    current_line = line
    stoplinelength = length(string(stopline))

    code = split(code[(startoffset+1):(stopoffset+1)],'\n')
    lineno = startline

    if !isempty(code) && isempty(code[end])
        pop!(code)
    end

    for textline in code
        printstyled(io,
            string(lineno, " "^(stoplinelength-length(lineno)+1));
            color = lineno == current_line ? :yellow : :bold)
        println(io, textline)
        lineno += 1
    end
    println(io)
end

function maybe_quote(x)
    (isa(x, Expr) || isa(x, Symbol)) ? QuoteNode(x) : x
end

function print_next_state(io::IO, state::DebuggerState, frame::JuliaStackFrame)
    print(io, "About to run: ")
    expr = pc_expr(frame, frame.pc[])
    isa(expr, Expr) && (expr = copy(expr))
    if isexpr(expr, :(=))
        expr = expr.args[2]
    end
    if isexpr(expr, :call) || isexpr(expr, :return)
        expr.args = map(var->maybe_quote(@lookup(frame, var)), expr.args)
    end
    if isa(expr, Expr)
        for (i, arg) in enumerate(expr.args)
            try
                nbytes = length(repr(arg))
                if nbytes > max(40, div(200, length(expr.args)))
                    expr.args[i] = Suppressed("$nbytes bytes of output")
                end
            catch
                expr.args[i] = Suppressed("printing error")
            end
        end
    end
    print(io, expr)
    println(io)
end

print_status(io::IO, state::DebuggerState) = print_status(io, state, state.stack[state.level])
function print_status(io::IO, state::DebuggerState, frame::JuliaStackFrame)
    # Buffer to avoid flickering
    outbuf = IOContext(IOBuffer(), io)
    printstyled(outbuf, "In ", locdesc(frame), "\n"; color=:bold)
    loc = locinfo(frame)
    if loc !== nothing
        data = if isa(loc, BufferLocInfo)
                loc.data
            else
                VERSION < v"0.7" ? read(loc.filepath, String) :
                read(loc.filepath, String)
            end
        print_sourcecode(outbuf, data,
            loc.line, loc.defline)
    else
        buf = IOBuffer()
        # TODO: look at the = 0
        active_line = 0
        code = split(String(take!(buf)),'\n')
        @assert active_line <= length(code)
        for (lineno, line) in enumerate(code)
            if lineno == active_line
                printstyled(outbuf, "=> ", bold = true, color=:yellow); println(outbuf, line)
            else
                printstyled(outbuf, "?  ", bold = true); println(outbuf, line)
            end
        end
    end
    print_next_state(outbuf, state, frame)
    print(io, String(take!(outbuf.io)))
end

function eval_code(state::DebuggerState, frame::JuliaStackFrame, command::AbstractString)
    expr = Base.parse_input_line(command)
    if isexpr(expr, :toplevel)
        expr = expr.args[end]
    end
    local_vars = Any[]
    local_vals = Any[]
    for i = 1:length(frame.locals)
        if !isa(frame.locals[i], Nothing)
            push!(local_vars, frame.code.code.slotnames[i])
            push!(local_vals, QuoteNode(something(frame.locals[i])))
        end
    end
    ismeth = frame.code.scope isa Method
    ismeth && (syms = sparam_syms(frame.code.scope))
    for i = 1:length(frame.sparams)
        ismeth && push!(local_vars, syms[i])
        push!(local_vals, QuoteNode(frame.sparams[i]))
    end
    res = gensym()
    eval_expr = Expr(:let,
        Expr(:block, map(x->Expr(:(=), x...), zip(local_vars, local_vals))...),
        Expr(:block,
            Expr(:(=), res, expr),
            Expr(:tuple, res, Expr(:tuple, local_vars...))
        ))
    eval_res, res = Core.eval(moduleof(frame), eval_expr)
    j = 1
    for i = 1:length(frame.locals)
        if !isa(frame.locals[i], Nothing)
            frame.locals[i] = Some{Any}(res[j])
            j += 1
        end
    end
    for i = 1:length(frame.sparams)
        frame.sparams[i] = res[j]
        j += 1
    end
    eval_res
end

@static if VERSION >= v"1.2.0-DEV.253"
    function eval_code(state::DebuggerState, code::AbstractString)
        try
            return eval_code(state, state.stack[1], code), false
        catch
            return Base.catch_stack(), true
        end
    end
else
    function eval_code(state::DebuggerState, code::AbstractString)
        try
            return true, eval_code(state, state.stack[1], code)
        catch err
            return false, (err, catch_backtrace())
        end
    end
end


text:
using WellIO
using Test

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

ext: jl


text:
@marginalrule AR(:y_x) (m_y::NormalDistributionsFamily, m_x::NormalDistributionsFamily, q_θ::NormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta) = begin
    return ar_y_x_marginal(getstype(meta), m_y, m_x, q_θ, q_γ, meta)
end

function ar_y_x_marginal(::ARsafe, m_y::NormalDistributionsFamily, m_x::NormalDistributionsFamily, q_θ::NormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta)
    mθ, Vθ = mean_cov(q_θ)
    mγ = mean(q_γ)

    mA = as_companion_matrix(mθ)
    mW = ar_precision(getvform(meta), getorder(meta), mγ)

    b_my, b_Vy = mean_cov(m_y)
    f_mx, f_Vx = mean_cov(m_x)

    inv_b_Vy = cholinv(b_Vy)
    inv_f_Vx = cholinv(f_Vx)

    D = inv_f_Vx + mγ * Vθ

    W_11 = add_precision(inv_b_Vy, mW)

    # Equvalent to -(mW * mA)
    W_12 = negate_inplace!(mW * mA)

    # Equivalent to (-mA' * mW)
    W_21 = negate_inplace!(mA' * mW)

    W_22 = D + mA' * mW * mA

    W = [ W_11 W_12; W_21 W_22 ]
    ξ = [ inv_b_Vy * b_my; inv_f_Vx * f_mx ]

    return MvNormalWeightedMeanPrecision(ξ, W)
end

function ar_y_x_marginal(::ARunsafe, m_y::NormalDistributionsFamily, m_x::NormalDistributionsFamily, q_θ::NormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta)
    mθ, Vθ = mean(q_θ), cov(q_θ)

    mA = as_companion_matrix(mθ)
    invmA = inv(mA)
    mγ = mean(q_γ)
    mV = ar_transition(getvform(meta), getorder(meta), mγ)

    b_my, b_Vy = mean(m_y), cov(m_y)
    f_mx, f_Vx = mean(m_x), cov(m_x)

    E = mV - mV*inv(b_Vy + mV)*mV
    F = mV + mV*invmA'*(inv(f_Vx) + mγ*Vθ)*invmA'*mV
    ABDC = E - E*inv(F + E)*E
    BD = -invmA' + invmA'*inv(invmA*mV*invmA' + inv((inv(f_Vx) + mγ*Vθ)))*invmA*mV*invmA'
    DC =  -invmA + invmA*mV*invmA'*inv(invmA*mV*invmA' + inv((inv(f_Vx) + mγ*Vθ)))*invmA
    D = invmA*mV*invmA' - invmA*mV*invmA'*inv(invmA*mV*invmA' + inv((inv(f_Vx) + mγ*Vθ)))*invmA*mV*invmA'
    invW = [ABDC -ABDC*BD; -DC*ABDC D+DC*ABDC*BD]

    m = invW*[inv(b_Vy)*b_my; inv(f_Vx)*f_mx]

    return MvNormalMeanCovariance(m, invW)
end


text:
module OTSM

export init_eye, init_lww1, init_sb, init_tb
export generate_procrustes_data
export otsm_pba, otsm_sdp, portwine_data, test_optimality

include("tracesum.jl")
include("initialization.jl")
include("portwine.jl")
include("generate_data.jl")

end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function calc_acoustic_sources!(acoustic_source_terms, u_euler, u_acoustics, vorticity_mean,
                                coupled_element_ids, mesh,
                                equations::AbstractCompressibleEulerEquations{2}, dg::DGSEM, cache)

  acoustic_source_terms .= zero(eltype(acoustic_source_terms))

  @threaded for k in 1:length(coupled_element_ids)
    element = coupled_element_ids[k]

    for j in eachnode(dg), i in eachnode(dg)
      vorticity = calc_vorticity_node(u_euler, mesh, equations, dg, cache, i, j, element)

      prim_euler = cons2prim(get_node_vars(u_euler, equations, dg, i, j, element), equations)
      v1 = prim_euler[2]
      v2 = prim_euler[3]
      v1_mean = u_acoustics[4, i, j, element]
      v2_mean = u_acoustics[5, i, j, element]

      vorticity_prime = vorticity - vorticity_mean[i, j, element]
      v1_prime = v1 - v1_mean
      v2_prime = v2 - v2_mean

      acoustic_source_terms[1, i, j, k] -= -vorticity_prime * v2_mean -
                                            vorticity_mean[i, j, element] * v2_prime
      acoustic_source_terms[2, i, j, k] -=  vorticity_prime * v1_mean +
                                            vorticity_mean[i, j, element] * v1_prime
    end
  end

  return nothing
end

end # @muladd


text:
using NixSourcerer.JuNix

opts = JuNix.Options(;
    nworkers=8, arch=Set(["x86_64"]), os=Set(["linux"]), libc=Set(["glibc"]), force_overwrite=true, check_store=true
)
x = JuNix.main(joinpath(@__DIR__, ".."), opts)
nothing


text:
function guess_struct(s::String)
	a=decode_get_parts(s)
	name=popfirst!(a)
	l=length(a)
	b=Vector{Any}(undef,l)
	for i in 1:l
		b[i]=decode(a[i])
	end
	f=eval(Symbol(name))
	return f(b...)
end
function work_command(s::String;commander=ply)
	v::Vector{String}=split_cmd(s)
	name=popfirst!(v)
	if haskey(command_alias,name)
		name=@inbounds command_alias[name]
	end
	if haskey(commands,name)
		pa::Tuple=@inbounds commands[name]
		types=pa[1]
		lt=length(types)
		lv=length(v)
		ld=pa[2]
		if lv>lt
			info_error(:cmd,:too_much_parameters;extra="$lv>$lt")
			return
		elseif lv<lt-ld
			info_error(:cmd,:too_little_parameters;extra="$lv<$lt-$ld")
			return
		end
		for i in 1:lv
			ty=types[i]
			if ty==:s # 字符串
			elseif ty==:i # 整数
				v[i]=parse(Int,v[i])
			elseif ty==:px # x坐标
				if v[i][1]=='~'
					v[i]=parse(Int,v[i][2:end])+commander.x
				else
					v[i]=parse(Int,v[i])
				end
			elseif ty==:py # y坐标
				if v[i][1]=='~'
					v[i]=parse(Int,v[i][2:end])+commander.y
				else
					v[i]=parse(Int,v[i])
				end
			elseif ty==:blk
				blk=guess_struct(v[i])
				if !isa(blk,Block)
					info_error(:cmd,:request_block;extra="#$i")
					return
				end
				v[i]=blk
			elseif ty==:itm
				sth=guess_struct(v[i])
				if isa(sth,Block)
					v[i]=IFB(sth)
				elseif !isa(sth,Item)
					info_error(:cmd,:request_item;extra="#$i")
					return
				end
			end
		end
		pa[3](v...)
	else
		info_error(:cmd,:unfound_command)
	end
end
global commands=Dict{Symbol,Tuple{Tuple,UInt8,Function}}(
	:clear=>((),0x0,info_clear),
	:eval=>((:s),0x0,(s::String)->eval(Meta.parse(s))),
	:fill=>((:blk,:px,:py,:i,:i),0x0,(b::Block,lx::Int,ly::Int,nx::Int,ny::Int;restricted::Bool=true)->begin
		if restricted&&nx>64&&ny>64
			info_help(:cmd,:restriction;extra=">64")
			return
		end
		for i in 1:nx
			for j in 1:ny
				setblk(i+lx,j+ly,deepcopy(b))
			end
		end
	end),
	:give=>((:itm,:i),0x1,(i::Item,num::Int=1)->give!(ply.ib,i,num))
	:say=>((:s),0x0,(s::String)->info_log("$(gsetting[:username]): $s")),
	:seed=>((),0x0,()->info_log(lsetting[:seed]))
)
const command_alias=Dict{Symbol,Symbol}(
	:clear_info=>:clear,
	:fill_override=>:fill,
	:show_seed=>:seed,
)


text:
using Distributed, StatsBase, Distributions, PyPlot
using TransitionalMCMC

lb  = -15
ub  = 15

fT(x) = logpdf(Uniform(lb,ub), x[1])
sample_fT(Nsamples) = rand(Uniform(lb,ub), Nsamples, 1)

log_fD_T(x) = log(pdf(Normal(0, 1), x[1]) + pdf(Normal(5, 0.2), x[1]))

Nsamples = 2000
samps, acc = tmcmc(log_fD_T, fT, sample_fT, Nsamples)

plt.hist(samps,50)


text:
@testset "HSAArray" begin
    agent = get_default_agent()

    arr = HSAArray(agent, UInt32, (8, ))
    fill!(arr, UInt32(1))
    _arr = Array(arr)
    @test length(arr) == length(_arr)
    for idx in 1:length(arr)
        _arr[idx] = arr[idx]
    end
    @test arr == _arr
    fill!(_arr, UInt32(2))
    for idx in 1:length(arr)
        arr[idx] = _arr[idx]
    end
    @test arr == _arr

    arr1 = HSAArray(agent, UInt32, (8,))
    arr2 = HSAArray(agent, UInt32, (8,))
    fill!(arr1, 1)
    fill!(arr2, 2)
    arr3 = similar(arr1)
    for i in 1:length(arr1)
        arr3[i] = arr1[i] + arr2[i]
    end
    @test all(arr3 .== 3)
    @test all(isa.((arr1, arr2, arr3), Ref(HSAArray)))
    @test size(arr3) == size(arr1)
end


text:
# this is the old upper triangular code

################################################################################
#
#  One echelonization step w/o transformations
#
################################################################################

# Perform one step in the echelonization, that is, find a pivot and clear
# elements below.
# sr = starting row
function _one_step(A::SMat{T}, sr = 1) where T
  nr, trafo = _one_step_with_transform(A, sr)
  return nr
end

function _one_step_with_transform(A::SMat{T}, sr = 1) where T
  trafos = SparseTrafoElem[]
  i = sr
  @assert i > 0
  all_r = Int[]
  min = A.c

  while i <= length(A.rows)
    # if we encounter zero row, push it to the end
    if iszero(A.rows[i])
      deleteat!(A.rows, i)
      push!(A.rows, sparse_row(base_ring(A)))
      push!(trafos, sparse_trafo_move_row(T, i, nrows(A)))
      i += 1
      continue
    end
    @assert length(A.rows[i].pos)>0
    @assert A.rows[i].pos[1] >= sr
    if A.rows[i].pos[1] < min
      min = A.rows[i].pos[1]
      all_r = Int[i]
    elseif A.rows[i].pos[1] == min
      push!(all_r, i)
    end
    i += 1
  end

  if length(all_r) == 0
    return A.r+1, trafos
  end

  #sorting
  sort!(all_r, lt = (a,b) -> (length(A.rows[a].pos) < length(A.rows[b].pos)))
  # we need to have rows[all_r[1]] == sr, so that the new pivot will be in
  # row sr

  if !(sr in all_r)
    push!(trafos, sparse_trafo_swap(T, sr, all_r[1]))
    q = A.rows[sr]
    A.rows[sr] = A.rows[all_r[1]]
    A.rows[all_r[1]] = q
    all_r[1] = sr
  else
    p = something(findfirst(isequal(sr), all_r), 0)
    if p > 1
      push!(trafos, sparse_trafo_swap(T, sr, all_r[1]))

      q = A.rows[sr]
      A.rows[sr] = A.rows[all_r[1]]
      A.rows[all_r[1]] = q

      all_r[p] = all_r[1]
      all_r[1] = sr
    end
  end

  if length(all_r) == 1
    return sr + 1, trafos
  end

  for j = 2:length(all_r)
    x = A.rows[sr].values[1]
    y = A.rows[all_r[j]].values[1]
    t = echelonize_basecase!(x, y, sr, all_r[j], A)
    push!(trafos, t)
  end

  sort!(all_r)

  for j = length(all_r):-1:2
    if length(A.rows[all_r[j]].pos) == 0
      deleteat!(A.rows, all_r[j])
      push!(A.rows, sparse_row(base_ring(A)))
      push!(trafos, sparse_trafo_move_row(T, all_r[j]), nrows(A))
o   end
  end
  return sr + 1, trafos
end

################################################################################
#
#  Echelonization via dense matrix
#
################################################################################

# The echelonization_via_dense must not remove zero rows automatically
function echelonize_via_dense(h::SMat{nmod})
  R = parent(h.rows[1].values[1])
  # turn into dense structure
  hdense = nmod_mat(h)
  # echelonize
  # I need to the following trick since we do not have the transformation
  rref!(hdense)
  # put back into sparse matrix over R (!= base_ring h)
  return sparse_matrix(hdense)
end

function echelonize_via_dense_with_transform(h::SMat{fmpz})
  hdense = matrix(h)
  # echelonize
  hdense, U = hnf_with_transform(hdense)
  # put back into sparse matrix
  h = sparse_matrix(hdense)
  return h, U
end

function echelonize_via_dense(h::SMat{fmpz})
  hdense = matrix(h)
  # echelonize
  hdense = hnf(hdense)
  # put back into sparse matrix
  h = sparse_matrix(hdense)
  return h
end

################################################################################
#
#  Base case of echelonization
#
################################################################################

# Given two rows i, j with entries x, y (x != 0) in the same column, this
# function must produce a zero at entry y.

function echelonize_basecase!(x::T, y::T, i::Int, j::Int, A::Hecke.SMat{T}) where T
  c = -divexact(y, x)
  add_scaled_row!(A, i, j, c)
  return sparse_trafo_addscaled(i, j, c)
end

function echelonize_basecase!(x::fmpz, y::fmpz, i::Int, j::Int,
                              A::Hecke.SMat{fmpz})
  if y % x == 0
    c = -div(y, x)
    add_scaled_row!(A, i, j, c)
    return sparse_trafo_add_scaled(i, j, c)
  else
    g,a,b = gcdx(x, y)
    c = -div(y, g)
    d = div(x, g)
    Hecke.transform_row!(A, i, j, a, b, c, d)
    return sparse_trafo_para_add_scaled(i, j, a, b, c, d)
  end
end

################################################################################
#
#  Echelonization for sparse matrices
#
################################################################################

# there is some density limit at level i, at which dense echelonization
# will be invoked
#

function upper_triangular_with_transform(M::SMat, density_limit = 0.5)
  B = deepcopy(M)
  T = upper_triangular_with_transform!(B, density_limit)
  return B, T
end

function upper_triangular_with_transform!(M::SMat, density_limit::Float64 = 0.5)
  f = (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit)
  return _upper_triangular_with_transform!(M, f)
end

function upper_triangular_with_transform!(M::SMat{fmpz},
                                      density_limit::Float64 = 0.5,
                                      size_limit::Int = 200)
  f =  (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit || nbits(maximum(abs, A)) > size_limit)
  return _upper_triangular_with_transform!(M, f)
end

# is_dense_enough is a function (::SMat{T}, i::Int) -> Bool
# At each level i, is_dense_enough(A, i) is called.
# If it evaluates to true, then dense echelonization will be called.
function _upper_triangular_with_transform!(A::SMat{T}, is_dense_enough::Function) where T
  trafo = SparseTrafoElem[]

  for i = 1:min(nrows(A), ncols(A))
    x, t = _one_step_with_transform(A, i)
    append!(trafo, t)

    if x > A.r
      return trafo
    end

    if is_dense_enough(A, i)
      h = sub(A, i:A.r, i:A.c)

      h, U = echelonize_via_dense_with_transform(h)

      # h will have zero rows

      push!(trafo, sparse_trafo_partial_dense(i, i:A.r, i:A.c, U))

      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r

      # append all rows of h to A
      # including zero rows

      for k in 1:length(h.rows)
        j = h.rows[k]
        if length(j.pos) == 0
          push!(A.rows, SRow{T}())
          A.r += 1
          continue
        end
        rw = SRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i - 1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      k = A.r

      # Now delete trailing zero rows.
      while length(A.rows[k].pos) == 0
        deleteat!(A.rows, k)
        push!(A.rows, sparse_row(base_ring(A)))
        push!(trafo, sparse_trafo_move_row(T, k, nrows(A)))
        k -= 1
      end

      # I guess I can do something more clever then the following:
      A.nnz = 0
      for r in A.rows
        A.nnz += length(r.pos)
      end

      return trafo
    end
  end
  return trafo
end

# now the same without transformation

function upper_triangular!(M::SMat, density_limit::Float64 = 0.5)
  f = (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit)
  return _upper_triangular!(M, f)
end

function upper_triangular!(M::SMat{fmpz}, density_limit::Float64 = 0.5, size_limit::Int = 200)
  f =  (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit || nbits(maximum(abs, A)) > size_limit)
  return _upper_triangular!(M, f)
end

function _upper_triangular!(A::SMat{T}, is_dense_enough) where T
  for i = 1:min(nrows(A), ncols(A))
    x = _one_step(A, i)

    if x>A.r
      return nothing
    end

    if is_dense_enough(A, i)

      h = sub(A, i:A.r, i:A.c)
      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r

      # make dense matrix
      h = echelonize_via_dense(h)

      for j in h.rows
        rw = SRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i - 1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      # I guess I can do something more clever then the following:
      A.nnz = 0
      for r in A.rows
        A.nnz += length(r.pos)
      end

      return nothing
    end
  end
  return nothing
end

################################################################################
#
#   Smith normal form
#
################################################################################

function _snf_upper_triangular_with_transform(A::SMat{fmpz})
  # We assume that A is quadratic and upper triangular
  essential_index = 1

  for i in 1:nrows(A)
    @assert A.rows[i].pos[1] == i
    if A.rows[i].values[1] != 1
      essential_index = i
      break
    end
  end

  # Let us pick up the transformations

  trafos_right = []
  trafos_left = []

  for i in 1:essential_index-1
    for (j, a) in zip(A.rows[i].pos, A.rows[i].values)
      if j == i
        continue
      end
      push!(trafos_right, sparse_trafo_add_scaled(j, i, -a))
    end
  end

  essential_part = fmpz_mat(sub(A, essential_index:nrows(A), essential_index:ncols(A)))

  snfofess, ltr, rtr = snf_with_transform(essential_part, true, true)

  push!(trafos_left, sparse_trafo_partial_dense(essential_index, essential_index:nrows(A), essential_index:ncols(A), ltr))
  push!(trafos_right, sparse_trafo_partial_dense(essential_index, essential_index:nrows(A), essential_index:ncols(A), rtr))

  return snfofess, trafos_left, trafos_right
end

@doc Markdown.doc"""
    elementary_divisors(A::SMat{fmpz}) -> Array{fmpz, 1}

The elementary divisors of $A$, ie. the diagonal elements of the Smith normal
form of $A$. $A$ needs to be of full rank - currently.
"""
function elementary_divisors(A::SMat{fmpz})
  A = copy(A)
  upper_triangular(A)

  essential_index = 1

  for i in 1:nrows(A)
    @assert A.rows[i].pos[1] == i
    if A.rows[i].values[1] != 1
      essential_index = i
      break
    end
  end

  essential_part = fmpz_mat(sub(A, essential_index:nrows(A), essential_index:ncols(A)))

  s = snf(essential_part)

  return vcat(fmpz[1 for i=1:essential_index-1], fmpz[s[i,i] for i=1:nrows(s)])
end


text:
using CUDAnative
# using Distributions
# using Distributions: zval, log2π

Base.@irrational log2π  1.8378770664093454836 log(big(2.)*π)
Base.@irrational sqrt2π 2.5066282746310005024 sqrt(big(π)*2.0)
# Base.@irrational invsqrt2π 0.3989422804014326779 inv(big(sqrt2π))

zval(μ::Real, σ::Real, x::Number) = (x - μ) / σ

# # pdf
# cunormpdf(z::Number) = CUDAnative.exp(-CUDAnative.abs2(z)/2) * invsqrt2π
# function cunormpdf(μ::Real, σ::Real, x::Number)
#     if iszero(σ)
#         if x == μ
#             z = zval(μ, one(σ), x)
#         else
#             z = zval(μ, σ, x)
#             σ = one(σ)
#         end
#     else
#         z = zval(μ, σ, x)
#     end
#     cunormpdf(z) / σ
# end

# logpdf
cunormlogpdf(x::Real) = -(x^2 + log2π)/2

function cunormlogpdf(p::Tuple{Real,Real}, x::Real)
    μ, σ = p
    if iszero(σ)
        if x == μ
            z = zval(μ, one(σ), x)
        else
            z = zval(μ, σ, x)
            σ = one(σ)
        end
    else
        z = zval(μ, σ, x)
    end
    cunormlogpdf(z) - CUDAnative.log(σ)
end


text:
__precompile__(true)

module Gnuplot

using StructC14N, ColorTypes, Printf, StatsBase, ReusePatterns, DataFrames

import Base.reset
import Base.println
import Base.iterate
import Base.convert

export gnuplot, quit, quitall, setverbose,
    @gp, @gsp, exec, save, contourlines, hist

#_____________________________________________________________________
#                          TYPE DEFINITIONS
#_____________________________________________________________________

# --------------------------------------------------------------------
mutable struct DataSource
    name::String
    lines::Vector{String}
end

# --------------------------------------------------------------------
mutable struct SinglePlot
    cmds::Vector{String}
    elems::Vector{String}
    flag3d::Bool
    SinglePlot() = new(Vector{String}(), Vector{String}(), false)
end

# --------------------------------------------------------------------
@quasiabstract mutable struct DrySession
    sid::Symbol                # session ID
    datas::Vector{DataSource}  # data sources
    plots::Vector{SinglePlot}  # commands and plot commands (one entry for eahelemec plot of the multiplot)
    curmid::Int                # current multiplot ID
end

# --------------------------------------------------------------------
@quasiabstract mutable struct Session <: DrySession
    pin::Base.Pipe;
    pout::Base.Pipe;
    perr::Base.Pipe;
    proc::Base.Process;
    channel::Channel{String};
end

# --------------------------------------------------------------------
mutable struct State
    sessions::Dict{Symbol, DrySession};
    dry::Bool
    cmd::String
    default::Symbol;           # default session name
    initcmd::Vector{String}    # Commands to initialize the gnuplot session (e.g., to set default terminal)
    verbose::Bool;             # verbosity level (true/false)
    silent::Bool;              # Silent flag
    printlines::Int;           # How many data lines are printed in log
    State() = new(Dict{Symbol, DrySession}(), true, "gnuplot", :default, Vector{String}(), false, false, 4)
end
const state = State()

#_____________________________________________________________________
#                 "PRIVATE" (NON-EXPORTED) FUNCTIONS
#_____________________________________________________________________

# --------------------------------------------------------------------
"""
  # CheckGnuplotVersion

  Check whether gnuplot is runnable with the command given in `cmd`.
  Also check that gnuplot version is >= 4.7 (required to use data
  blocks).
"""
function CheckGnuplotVersion(cmd::AbstractString)
    icmd = `$(cmd) --version`

    proc = open(`$icmd`, read=true)
    s = String(read(proc))
    if !success(proc)
        error("An error occurred while running: " * string(icmd))
    end

    s = split(s, " ")
    ver = ""
    for token in s
        try
            ver = VersionNumber("$token")
            break
        catch
        end
    end

    if ver < v"4.7"
        error("gnuplot ver. >= 4.7 is required, but " * string(ver) * " was found.")
    end
    #@info "  Gnuplot version: " * string(ver)
    return ver
end

# --------------------------------------------------------------------
function __init__()
    global state
    try
        ver = CheckGnuplotVersion(state.cmd)
        state.dry = false
    catch err

    end
end

# --------------------------------------------------------------------
function parseKeywords(; kwargs...)
    template = (xrange=NTuple{2, Real},
                yrange=NTuple{2, Real},
                zrange=NTuple{2, Real},
                cbrange=NTuple{2, Real},
                title=AbstractString,
                xlabel=AbstractString,
                ylabel=AbstractString,
                zlabel=AbstractString,
                xlog=Bool,
                ylog=Bool,
                zlog=Bool)

    kw = canonicalize(template; kwargs...)
    out = Vector{String}()
    ismissing(kw.xrange ) || (push!(out, "set xrange  [" * join(kw.xrange , ":") * "]"))
    ismissing(kw.yrange ) || (push!(out, "set yrange  [" * join(kw.yrange , ":") * "]"))
    ismissing(kw.zrange ) || (push!(out, "set zrange  [" * join(kw.zrange , ":") * "]"))
    ismissing(kw.cbrange) || (push!(out, "set cbrange [" * join(kw.cbrange, ":") * "]"))
    ismissing(kw.title  ) || (push!(out, "set title  '" * kw.title  * "'"))
    ismissing(kw.xlabel ) || (push!(out, "set xlabel '" * kw.xlabel * "'"))
    ismissing(kw.ylabel ) || (push!(out, "set ylabel '" * kw.ylabel * "'"))
    ismissing(kw.zlabel ) || (push!(out, "set zlabel '" * kw.zlabel * "'"))
    ismissing(kw.xlog   ) || (push!(out, (kw.xlog  ?  ""  :  "un") * "set logscale x"))
    ismissing(kw.ylog   ) || (push!(out, (kw.ylog  ?  ""  :  "un") * "set logscale y"))
    ismissing(kw.zlog   ) || (push!(out, (kw.zlog  ?  ""  :  "un") * "set logscale z"))
    return out
end

# --------------------------------------------------------------------
function DrySession(sid::Symbol)
    global state
    (sid in keys(state.sessions))  &&
        error("Gnuplot session $sid is already active")
    out = DrySession(sid, Vector{DataSource}(), [SinglePlot()], 1)
    return out
end

# --------------------------------------------------------------------
function getsession(sid::Symbol)
    global state
    if !(sid in keys(state.sessions))
        #@info "Creating session $sid..."
        if state.dry
            gnuplot(sid)
        else
            gnuplot(sid, state.cmd)
        end
    end
    return state.sessions[sid]
end
function getsession()
    global state
    return getsession(state.default)
end

# --------------------------------------------------------------------
"""
  # println

  Send a string to gnuplot's STDIN.

  The commands sent through `send` are not stored in the current
  session (use `newcmd` to save commands in the current session).

  ## Arguments:
  - `gp`: a `Session` object;
  - `str::String`: command to be sent;
  - `capture=false`: set to `true` to capture and return the output.
"""
println(gp::DrySession, str::AbstractString) = nothing
function println(gp::Session, str::AbstractString)
    global state
    if !state.silent  &&  state.verbose
        printstyled(color=:yellow, "GNUPLOT ($(gp.sid)) $str\n")
    end
    w = write(gp.pin, strip(str) * "\n")
    w <= 0  &&  error("Writing on gnuplot STDIN pipe returned $w")
    #flush(gp.pin)
    return w
end

println(gp::DrySession, d::DataSource) = nothing

function println(gp::Session, d::DataSource)
    if typeof(gp) == concretetype(Session)
        if !state.silent  &&  state.verbose
            for ii in 1:length(d.lines)
                v = d.lines[ii]
                printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) $v\n")
                if ii == state.printlines
                    printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) ...\n")
                    if ii < length(d.lines)
                        v = d.lines[end]
                        printstyled(color=:light_black, "GNUPLOT ($(gp.sid)) $v\n")
                    end
                    break
                end
            end
        end
        w = write(gp.pin, "\n")
        w = write(gp.pin, join(d.lines, "\n"))
        w = write(gp.pin, "\n")
        w = write(gp.pin, "\n")
        flush(gp.pin)
    end
    return nothing
end

# --------------------------------------------------------------------
writeread(gp::DrySession, str::AbstractString) = [""]
function writeread(gp::Session, str::AbstractString)
    global state
    silent = state.silent
    state.silent = true
    write(gp.pin, "print 'GNUPLOT_CAPTURE_BEGIN'\n")
    println(gp, strip(str))
    write(gp.pin, "print 'GNUPLOT_CAPTURE_END'\n")
    flush(gp.pin)
    out = Vector{String}()
    while true
        l = take!(gp.channel)
        l == "GNUPLOT_CAPTURE_END"  &&  break
        push!(out, l)
    end
    state.silent = silent
    return out
end

# --------------------------------------------------------------------
setWindowTitle(session::DrySession) = nothing
function setWindowTitle(session::Session)
    term = writeread(session, "print GPVAL_TERM")[1]
    if term in ("aqua", "x11", "qt", "wxt")
        opts = writeread(session, "print GPVAL_TERMOPTIONS")[1]
        if findfirst("title", opts) == nothing
            writeread(session, "set term $term $opts title 'Gnuplot.jl: $(session.sid)'")
        end
    end
end

# --------------------------------------------------------------------
function reset(gp::DrySession)
    gp.datas = Vector{DataSource}()
    gp.plots = [SinglePlot()]
    gp.curmid = 1
    println(gp, "reset session")
    return nothing
end

# --------------------------------------------------------------------
function setmulti(gp::DrySession, mid::Int)
    @assert mid >= 0 "Multiplot ID must be a >= 0"
    for i in length(gp.plots)+1:mid
        push!(gp.plots, SinglePlot())
    end
    (mid > 0)  &&  (gp.curmid = mid)
end

# --------------------------------------------------------------------
function newdatasource(gp::DrySession, v::Vector{T}; name="") where T <: AbstractString
    (name == "")  &&  (name = string("data", length(gp.datas)))
    name = "\$$name"

    accum = Vector{String}()
    push!(accum, "$name << EOD")
    append!(accum, v)
    push!(accum, "EOD")
    d = DataSource(name, accum)
    push!(gp.datas, d)

    println(gp, d) # Send directly to gnuplot process
    return name
end

# --------------------------------------------------------------------
function newdatasource(gp::DrySession, args...; name="")
    (name == "")  &&  (name = string("data", length(gp.datas)))
    name = "\$$name"

    # Check dimensions
    maxDim = 0
    for iarg in 1:length(args)
        d = args[iarg]

        ok = false
        if typeof(d) <: Number
            ok = true
        elseif typeof(d) <: AbstractArray
            if typeof(d[1]) <: Number
                ok = true
            end
            if typeof(d[1]) <: ColorTypes.RGB
                ok = true
            end
        end
        if ndims(d) > maxDim
            maxDim = ndims(d)
        end

        @assert ok "Invalid argument at position $iarg"
        @assert maxDim <= 3 "Array dimensions must be <= 3"
    end

    dimX = 0
    dimY = 0
    dimZ = 0
    count1D = 0
    for iarg in 1:length(args)
        d = args[iarg]
        if ndims(d) == 0
            @assert maxDim == 0 "Input data are ambiguous: use use all scalar floats or arrays of floats"
        elseif ndims(d) == 1
            count1D += 1
            if maxDim == 1
                (iarg == 1)  &&  (dimX = length(d))
                @assert dimX == length(d) "Array size are incompatible"
            else
                (iarg == 1)  &&  (dimX = length(d))
                (maxDim == 2)  &&  (iarg == 2)  &&  (dimY = length(d))
                (maxDim == 3)  &&  (iarg == 3)  &&  (dimZ = length(d))
                @assert iarg <= maxDim "2D and 3D data must be given at the end of argument list"
            end
        elseif ndims(d) == 2
            if iarg == 1
                dimX = (size(d))[1]
                dimY = (size(d))[2]
            end
            @assert dimX == (size(d))[1] "Array size are incompatible"
            @assert dimY == (size(d))[2] "Array size are incompatible"
            @assert dimZ == 0 "Mixing 2D and 3D data is not allowed"
        elseif ndims(d) == 3
            if iarg == 1
                dimX = (size(d))[1]
                dimY = (size(d))[2]
                dimZ = (size(d))[3]
            end
            @assert dimX == (size(d))[1] "Array size are incompatible"
            @assert dimY == (size(d))[2] "Array size are incompatible"
            @assert dimZ == (size(d))[3] "Array size are incompatible"
        end
    end
    if (dimZ > 0)  &&  (count1D != 0)  &&  (count1D != 3)
        error("Either zero or three 1D arrays must be given before 3D data")
    elseif (dimY > 0)  &&  (count1D != 0)  &&  (count1D != 2)
        error("Either zero or two 1D arrays must be given before 2D data")
    end

    # Prepare data
    accum = Vector{String}()
    v = "$name << EOD"
    push!(accum, v)

    if dimZ > 0 # 3D
        for ix in 1:dimX
            for iy in 1:dimY
                for iz in 1:dimZ
                    if count1D == 0
                        v = string(ix) * " " * string(iy) * " " * string(iz)
                    else
                        v = string(args[1][ix]) * " " * string(args[2][iy]) * " " * string(args[3][iz])
                    end
                    for iarg in count1D+1:length(args)
                        d = args[iarg]
                        v *= " " * string(d[ix,iy,iz])
                    end
                    push!(accum, v)
                end
            end
            push!(accum, "")
        end
    elseif dimY > 0  # 2D
        for ix in 1:dimX
            for iy in 1:dimY
                if count1D == 0
                    v = string(ix) * " " * string(iy)
                else
                    v = string(args[1][ix]) * " " * string(args[2][iy])
                end
                for iarg in count1D+1:length(args)
                    d = args[iarg]
                    if typeof(d[ix,iy]) <: ColorTypes.RGB
                        tmp = d[ix,iy]
                        v *= " " * string(float(tmp.r)*255) * " " * string(float(tmp.g)*255) * " " * string(float(tmp.b)*255)
                    else
                        v *= " " * string(d[ix,iy])
                    end
                end
                push!(accum, v)
            end
            push!(accum, "")
        end
    elseif dimX > 0  # 1D
        for ix in 1:dimX
            v = ""
            for iarg in 1:length(args)
                d = args[iarg]
                v *= " " * string(d[ix])
            end
            push!(accum, v)
        end
    else # scalars
        v = ""
        for iarg in 1:length(args)
            d = args[iarg]
            v *= " " * string(d)
        end
        push!(accum, v)
    end

    push!(accum, "EOD")
    d = DataSource(name, accum)
    push!(gp.datas, d)

    println(gp, d) # Send directly to gnuplot process
    return name
end

# --------------------------------------------------------------------
"""
  # newcmd

  Send a command to gnuplot process and store it in the current session.
"""
function newcmd(gp::DrySession, v::String; mid::Int=0)
    setmulti(gp, mid)
    (v != "")  &&  (push!(gp.plots[gp.curmid].cmds, v))
    (length(gp.plots) == 1)  &&  (println(gp, v))
    return nothing
end

function newcmd(gp::DrySession; mid::Int=0, args...)
    for v in parseKeywords(;args...)
        newcmd(gp, v, mid=mid)
    end
    return nothing
end

# --------------------------------------------------------------------
function newplotelem(gp::DrySession, name, opt=""; mid=0)
    setmulti(gp, mid)
    push!(gp.plots[gp.curmid].elems, "$name $opt")
end

# --------------------------------------------------------------------
function quitsession(gp::DrySession)
    global state
    delete!(state.sessions, gp.sid)
    return 0
end

function quitsession(gp::Session)
    close(gp.pin)
    close(gp.pout)
    close(gp.perr)
    wait( gp.proc)
    exitCode = gp.proc.exitcode
    invoke(quitsession, Tuple{DrySession}, gp)
    return exitCode
end

# --------------------------------------------------------------------
iterate(gp::DrySession) = ("#ID: $(gp.sid)", (true, 1, 1))
function iterate(gp::DrySession, state)
    (onData, mid, ii) = state
    if onData
        if mid <= length(gp.datas)
            if ii <= length(gp.datas[mid].lines)
                return (gp.datas[mid].lines[ii], (true, mid, ii+1))
            end
            return iterate(gp, (true, mid+1, 1))
        end
        return ("", (false, 1, 1))
    end

    if mid <= length(gp.plots)
        if ii <= length(gp.plots[mid].cmds)
            return (gp.plots[mid].cmds[ii], (false, mid, ii+1))
        end
        if length(gp.plots[mid].elems) == 0
            return ("", (false, mid+1, 1))
        end
        s = (gp.plots[mid].flag3d  ?  "splot "  :  "plot ") * " \\\n  " *
            join(gp.plots[mid].elems, ", \\\n  ")
        return (s, (false, mid+1, 1))
    end

    if mid == length(gp.plots)+1
        s = ""
        (length(gp.plots) > 1)  &&  (s *= "unset multiplot;")
        return (s, (false, mid+1, 1))
    end

    return nothing
end

# --------------------------------------------------------------------
# function convert(::Type{Vector{String}}, gp::DrySession)
#     out = Vector{String}()
#     for l in gp
#         push!(out, l)
#     end
#     return out
# end

# --------------------------------------------------------------------
internal_save(gp::DrySession; kw...) = internal_save(gp, gp; kw...)
function internal_save(gp::DrySession, stream; term::AbstractString="", output::AbstractString="")
    println(stream, "reset")
    if term != ""
        former_term = writeread(gp, "print GPVAL_TERM")[1]
        former_opts = writeread(gp, "print GPVAL_TERMOPTIONS")[1]
        println(stream, "set term $term")
    end
    (output != "")  &&  println(stream, "set output '$output'")
    i = (!(typeof(stream) <: DrySession), 1, 1) # Skip data sources ?
    while (next = iterate(gp, i)) != nothing
        (s, i) = next
        println(stream, s)
    end
    (output != "")  &&  println(stream, "set output")
    if term != ""
        println(stream, "set term $former_term $former_opts")
    end
    output
end

# --------------------------------------------------------------------
function driver(args...; flag3d=false)
    if length(args) == 0
        gp = getsession()
        internal_save(gp)
        return nothing
    end

    data = Vector{Any}()
    dataname = ""
    dataplot = nothing

    function dataCompleted()
        if length(data) > 0
            AllArraysAreNotEmpty = true
            for i in 1:length(data)
                if (typeof(data[i]) <: AbstractArray)  &&  (length(data[i]) == 0)
                    @warn "Input array is empty"
                    AllArraysAreNotEmpty = false
                    break
                end
            end
            if AllArraysAreNotEmpty
                last = newdatasource(gp, data...; name=dataname)
                (dataplot != nothing)  &&  (newplotelem(gp, last, dataplot))
            end
        end
        data = Vector{Any}()
        dataname = ""
        dataplot = nothing
    end
    function isPlotCmd(s::String)
        (length(s) >= 2)  &&  (s[1:2] ==  "p "    )  &&  (return (true, false, strip(s[2:end])))
        (length(s) >= 3)  &&  (s[1:3] ==  "pl "   )  &&  (return (true, false, strip(s[3:end])))
        (length(s) >= 4)  &&  (s[1:4] ==  "plo "  )  &&  (return (true, false, strip(s[4:end])))
        (length(s) >= 5)  &&  (s[1:5] ==  "plot " )  &&  (return (true, false, strip(s[5:end])))
        (length(s) >= 2)  &&  (s[1:2] ==  "s "    )  &&  (return (true, true , strip(s[2:end])))
        (length(s) >= 3)  &&  (s[1:3] ==  "sp "   )  &&  (return (true, true , strip(s[3:end])))
        (length(s) >= 4)  &&  (s[1:4] ==  "spl "  )  &&  (return (true, true , strip(s[4:end])))
        (length(s) >= 5)  &&  (s[1:5] ==  "splo " )  &&  (return (true, true , strip(s[5:end])))
        (length(s) >= 6)  &&  (s[1:6] ==  "splot ")  &&  (return (true, true , strip(s[6:end])))
        return (false, false, "")
    end

    gp = nothing
    doDump  = true
    doReset = true

    for loop in 1:2
        if loop == 2
            (gp == nothing)  &&  (gp = getsession())
            doReset  &&  reset(gp)
            gp.plots[gp.curmid].flag3d = flag3d
        end

        for iarg in 1:length(args)
            arg = args[iarg]

            if typeof(arg) == Symbol
                if arg == :-
                    (loop == 1)  &&  (iarg < length(args)) &&  (doReset = false)
                    (loop == 1)  &&  (iarg >  1)           &&  (doDump  = false)
                else
                    (loop == 1)  &&  (gp = getsession(arg))
                end
            elseif isa(arg, Tuple)  &&  length(arg) == 2  &&  isa(arg[1], Symbol)
                if arg[1] == :term
                    if loop == 1
                        if typeof(arg[2]) == String
                            term = (deepcopy(arg[2]), "")
                        elseif length(arg[2]) == 2
                            term = deepcopy(arg[2])
                        else
                            error("The term tuple must contain at most two strings")
                        end
                    end
                else
                    (loop == 2)  &&  newcmd(gp; [arg]...) # A cmd keyword
                end
            elseif isa(arg, Int)
                (loop == 2)  &&  (@assert arg > 0)
                (loop == 2)  &&  (dataplot = ""; dataCompleted())
                (loop == 2)  &&  setmulti(gp, arg)
            elseif isa(arg, String)
                # Either a dataname, a plot or a command
                if loop == 2
                    if arg[1] == '$'
                        dataname = arg[2:end]
                        dataCompleted()
                    elseif length(data) > 0
                        dataplot = arg
                        dataCompleted()
                    else
                        (isPlot, flag3d, cmd) = isPlotCmd(arg)
                        if isPlot
                            gp.plots[gp.curmid].flag3d = flag3d
                            newplotelem(gp, cmd)
                        else
                            newcmd(gp, arg)
                        end
                    end
                end
            else
                (loop == 2)  &&  push!(data, arg) # a data set
            end
        end
    end

    dataplot = ""
    dataCompleted()
    (doDump)  &&  (internal_save(gp))

    return nothing
end

#_____________________________________________________________________
#                         EXPORTED FUNCTIONS
#_____________________________________________________________________

# --------------------------------------------------------------------
"""
  `gnuplot(sid::Symbol)`

  `gnuplot(sid::Symbol, cmd::AbstractString)`

  Initialize a new session and (optionally) the associated Gnuplot process

  ## Arguments:
  - `sid`: the session name (a Julia symbol);
  - `cmd`: a string specifying the complete file path to a gnuplot executable.  If not given a *dry* session will be created, i.e. a session without underlying gnuplot process.
"""
function gnuplot(sid::Symbol)
    out = DrySession(sid)
    state.sessions[sid] = out
    return out
end

function gnuplot(sid::Symbol, cmd::AbstractString)
    function readTask(sid, stream, channel)
        global state
        saveOutput = false

        while isopen(stream)
            line = readline(stream)

            if (length(line) >= 1)  &&  (line[1] == Char(0x1b)) # Escape (xterm -ti vt340)
                buf = Vector{UInt8}()
                append!(buf, convert(Vector{UInt8}, [line...]))
                push!(buf, 0x0a)
                c = 0x00
                while c != 0x1b
                    c = read(stream, 1)[1]
                    push!(buf, c)
                end
                c = read(stream, 1)[1]
                push!(buf, c)
                write(stdout, buf)
                continue
            end
            if line == "GNUPLOT_CAPTURE_BEGIN"
                saveOutput = true
            else
                if (line != "")  &&  (line != "GNUPLOT_CAPTURE_END")  &&  (!state.silent)
                    printstyled(color=:cyan, "GNUPLOT ($sid) -> $line\n")
                end
                (saveOutput)  &&  (put!(channel, line))
                (line == "GNUPLOT_CAPTURE_END")  &&  (saveOutput = false)
            end
        end
        delete!(state.sessions, sid)
        return nothing
    end

    global state

    CheckGnuplotVersion(cmd)
    session = DrySession(sid)

    pin  = Base.Pipe()
    pout = Base.Pipe()
    perr = Base.Pipe()
    proc = run(pipeline(`$cmd`, stdin=pin, stdout=pout, stderr=perr), wait=false)
    chan = Channel{String}(32)

    # Close unused sides of the pipes
    Base.close(pout.in)
    Base.close(perr.in)
    Base.close(pin.out)
    Base.start_reading(pout.out)
    Base.start_reading(perr.out)

    # Start reading tasks
    @async readTask(sid, pout, chan)
    @async readTask(sid, perr, chan)

    out = Session(getfield.(Ref(session), fieldnames(concretetype(DrySession)))..., pin, pout, perr, proc, chan)
    state.sessions[sid] = out

    setWindowTitle(out)
    for l in state.initcmd
        writeread(out, l)
    end

    return out
end

function gnuplot()
    global state
    return gnuplot(state.default)
end
function gnuplot(cmd::AbstractString)
    global state
    return gnuplot(state.default, cmd)
end

# --------------------------------------------------------------------
"""
  `quit()`

  Quit the session and the associated gnuplot process (if any).
"""
function quit(sid::Symbol)
    global state
    if !(sid in keys(state.sessions))
        error("Gnuplot session $sid do not exists")
    end
    return quitsession(state.sessions[sid])
end

"""
  `quitall()`

  Quit all the sessions and the associated gnuplot processes.
"""
function quitall()
    global state
    for sid in keys(state.sessions)
        quit(sid)
    end
    return nothing
end

# --------------------------------------------------------------------
"""
`@gp args...`

The `@gp` macro, and its companion `@gsp` (for `splot` operations) allows to exploit all of the **Gnuplot** package functionalities using an extremely efficient and concise syntax.  Both macros accept the same syntax, as described below.

The macros accepts any number of arguments, with the following meaning:
- a symbol: the name of the session to use;
- a string: a command (e.g. "set key left") or plot specification (e.g. "with lines");
- a string starting with a `\$` sign: a data set name;
- an `Int` > 0: the plot destination in a multiplot session;
- a keyword/value pair: a keyword value (see below);
- any other type: a dataset to be passed to Gnuplot.  Each dataset must be terminated by either:
  - a string starting with a `\$` sign (i.e. the data set name);
  - or a string with the plot specifications (e.g. "with lines");
- the `:-` symbol, used as first argument, avoids resetting the Gnuplot session.  Used as last argument avoids immediate execution  of the plot/splot command.  This symbol can be used to split a  single call into multiple ones.

All entries are optional, and there is no mandatory order.  The plot specification can either be:
 - a complete plot/splot command (e.g., "plot sin(x)", both "plot" and "splot" can be abbreviated to "p" and "s" respectively);
 - or a partial specification starting with the "with" clause (if it follows a data set).

The list of accepted keyword is as follows:
- `title::String`: plot title;
- `xlabel::String`: X axis label;
- `ylabel::String`: Y axis label;
- `zlabel::String`: Z axis label;
- `xlog::Bool`: logarithmic scale for X axis;
- `ylog::Bool`: logarithmic scale for Y axis;
- `zlog::Bool`: logarithmic scale for Z axis;
- `xrange::NTuple{2, Number}`: X axis range;
- `yrange::NTuple{2, Number}`: Y axis range;
- `zrange::NTuple{2, Number}`: Z axis range;
- `cbrange::NTuple{2, Number}`: Color box axis range;

The symbol for the above-mentioned keywords may also be used in a shortened form, as long as there is no ambiguity with other keywords.  E.g. you can use: `xr=(1,10)` in place of `xrange=(1,10)`.

# Examples:

## Simple examples with no data:
```
@gp "plot sin(x)"
@gp "plot sin(x)" "pl cos(x)"
@gp "plo sin(x)" "s cos(x)"

# Split a `@gp` call in two
@gp "plot sin(x)" :-
@gp :- "plot cos(x)"

# Insert a 3 second pause between one plot and the next
@gp "plot sin(x)" 2 xr=(-2pi,2pi) "pause 3" "plot cos(4*x)"
```

### Simple examples with data:
```
@gp "set key left" tit="My title" xr=(1,12) 1:10 "with lines tit 'Data'"

x = collect(1.:10)
@gp x
@gp x x
@gp x -x
@gp x x.^2
@gp x x.^2 "w l"

lw = 3
@gp x x.^2 "w l lw \$lw"
```

### A more complex example
```
@gp("set grid", "set key left", xlog=true, ylog=true,
    title="My title", xlab="X label", ylab="Y label",
    x, x.^0.5, "w l tit 'Pow 0.5' dt 2 lw 2 lc rgb 'red'",
    x, x     , "w l tit 'Pow 1'   dt 1 lw 3 lc rgb 'blue'",
    x, x.^2  , "w l tit 'Pow 2'   dt 3 lw 2 lc rgb 'purple'")
```

### Multiplot example:
```
@gp(xr=(-2pi,2pi), "unset key",
    "set multi layout 2,2 title 'Multiplot title'",
    1, "p sin(x)"  ,
    2, "p sin(2*x)",
    3, "p sin(3*x)",
    4, "p sin(4*x)")
```
or equivalently
```
@gp xr=(-2pi,2pi) "unset key" "set multi layout 2,2 title 'Multiplot title'" :-
for i in 1:4
  @gp :- i "p sin(\$i*x)" :-
end
@gp
```

### Multiple gnuplot sessions
```
@gp :GP1 "plot sin(x)"
@gp :GP2 "plot sin(x)"

quitall()
```

### Further examples
```
x = range(-2pi, stop=2pi, length=100);
y = 1.5 * sin.(0.3 .+ 0.7x) ;
noise = randn(length(x))./2;
e = 0.5 * fill(1, size(x));

name = "\\\$MyDataSet1"
@gp x y name "plot \$name w l" "pl \$name u 1:(2*\\\$2) w l"

@gsp randn(Float64, 30, 50)
@gp randn(Float64, 30, 50) "w image"
@gsp x y y

@gp("set key horizontal", "set grid",
    xrange=(-7,7), ylabel="Y label",
    x, y, "w l t 'Real model' dt 2 lw 2 lc rgb 'red'",
    x, y+noise, e, "w errorbars t 'Data'")

@gp "f(x) = a * sin(b + c*x); a = 1; b = 1; c = 1;"   :-
@gp :- x y+noise e name                               :-
@gp :- "fit f(x) \$name u 1:2:3 via a, b, c;"         :-
@gp :- "set multiplot layout 2,1"                     :-
@gp :- "plot \$name w points" ylab="Data and model"   :-
@gp :- "plot \$name u 1:(f(\\\$1)) w lines"           :-
@gp :- 2 xlab="X label" ylab="Residuals"              :-
@gp :- "plot \$name u 1:((f(\\\$1)-\\\$2) / \\\$3):(1) w errorbars notit"

# Retrieve values for a, b and c
a = parse(Float64, exec("print a"))
b = parse(Float64, exec("print b"))
c = parse(Float64, exec("print c"))

# Save to a PDF file
save(term="pdf", output="gnuplot.pdf")
```

### Display an image
```
using TestImages
img = testimage("lena");
@gp img "w image"
@gp "set size square" img "w rgbimage" # Color image with correct proportions
@gp "set size square" img "u 2:(-\\\$1):3:4:5 with rgbimage" # Correct orientation
```
"""
macro gp(args...)
    out = Expr(:call)
    push!(out.args, :(Gnuplot.driver))
    for iarg in 1:length(args)
        arg = args[iarg]
        if (isa(arg, Expr)  &&  (arg.head == :(=)))
            sym = string(arg.args[1])
            val = arg.args[2]
            push!(out.args, :((Symbol($sym),$val)))
        else
            push!(out.args, arg)
        end
    end
    return esc(out)
end

"""
  # @gsp

  See documentation for `@gp`.
"""
macro gsp(args...)
    out = Expr(:macrocall, Symbol("@gp"), LineNumberNode(1, "Gnuplot.jl"))
    push!(out.args, args...)
    push!(out.args, Expr(:kw, :flag3d, true))
    return esc(out)
end

# # --------------------------------------------------------------------
# """
#   # repl
#
#   Read/evaluate/print/loop
# """
# function repl(sid::Symbol)
#     verb = state.verbose
#     state.verbose = 0
#     gp = getsession(sid)
#     while true
#         line = readline(stdin)
#         (line == "")  &&  break
#         answer = send(gp, line, true)
#         for line in answer
#             println(line)
#         end
#     end
#     state.verbose = verb
#     return nothing
# end
# function repl()
#     global state
#     return repl(state.default)
# end

# --------------------------------------------------------------------
"""
`exec(sid::Symbol, s::Vector{String})`

Directly execute commands on the underlying Gnuplot process, and return the result(s).

## Examples:
```julia
exec("print GPVAL_TERM")
exec("plot sin(x)")
```
"""
function exec(sid::Symbol, s::Vector{String})
    global state
    gp = getsession(sid)
    answer = Vector{String}()
    for v in s
        push!(answer, writeread(gp, v)...)
    end
    return join(answer, "\n")
end
function exec(s::String)
    global state
    exec(state.default, [s])
end
exec(sid::Symbol, s::String) = exec(sid, [s])

# --------------------------------------------------------------------
"""
`setverbose(b::Bool)`

Set verbose flag to `true` or `false` (default: `false`).
"""
function setverbose(b::Bool)
    global state
    state.verbose = b
end

function initcmd()
    global state
    return state.initcmd
end

# --------------------------------------------------------------------
"""
`save(...)`

Save the data and commands in the current session to either:
- the gnuplot process (i.e. produce a plot): `save(term="", output="")`;
- an IO stream: `save(stream::IO; term="", output="")`;
- a file: `save(file::AbstractStrings; term="", output="")`.

To save the data and command from a specific session pass the ID as first argument, i.e.:
- `save(sid::Symbol, term="", output="")`;
- `save(sid::Symbol, stream::IO; term="", output="")`;
- `save(sid::Symbol, file::AbstractStrings; term="", output="")`.

In all cases the `term` keyword allows to specify a gnuplot terminal, and the `output` keyword allows to specify an output file.
"""
save(                                 ; kw...) =                            internal_save(getsession()           ; kw...)
save(sid::Symbol                      ; kw...) =                            internal_save(getsession(sid)        ; kw...)
save(             stream::IO          ; kw...) =                            internal_save(getsession()   , stream; kw...)
save(sid::Symbol, stream::IO          ; kw...) =                            internal_save(getsession(sid), stream; kw...)
save(             file::AbstractString; kw...) = open(file, "w") do stream; internal_save(getsession()   , stream; kw...); end
save(sid::Symbol, file::AbstractString; kw...) = open(file, "w") do stream; internal_save(getsession(sid), stream; kw...); end

# --------------------------------------------------------------------
#=
Example:
v = randn(1000)
h = hist(v, bs=0.2)
@gp h.loc h.counts "w histep" h.loc h.counts "w l"
=#
function hist(v::Vector{T}; range=[NaN,NaN], bs=NaN, nbins=0, pad=true) where T <: Number
    i = findall(isfinite.(v))
    isnan(range[1])  &&  (range[1] = minimum(v[i]))
    isnan(range[2])  &&  (range[2] = maximum(v[i]))
    i = findall(isfinite.(v)  .&  (v.>= range[1])  .&  (v.<= range[2]))
    (nbins > 0)  &&  (bs = (range[2] - range[1]) / nbins)
    if isfinite(bs)
        rr = range[1]:bs:range[2]
        if maximum(rr) < range[2]
            rr = range[1]:bs:(range[2]+bs)
        end
        hh = fit(Histogram, v[i], rr, closed=:left)
        if sum(hh.weights) < length(i)
            j = findall(v[i] .== range[2])
            @assert length(j) == (length(i) - sum(hh.weights))
            hh.weights[end] += length(j)
        end
    else
        hh = fit(Histogram, v[i], closed=:left)
    end
    @assert sum(hh.weights) == length(i)
    x = collect(hh.edges[1])
    x = (x[1:end-1] .+ x[2:end]) ./ 2
    h = hh.weights
    binsize = x[2] - x[1]
    if pad
        x = [x[1]-binsize, x..., x[end]+binsize]
        h = [0, h..., 0]
    end
    return (loc=x, counts=h, binsize=binsize)
end

function hist(v1::Vector{T1}, v2::Vector{T2};
              range1=[NaN,NaN], bs1=NaN, nbins1=0,
              range2=[NaN,NaN], bs2=NaN, nbins2=0) where {T1 <: Number, T2 <: Number}
    i = findall(isfinite.(v2))
    isnan(range1[1])  &&  (range1[1] = minimum(v1[i]))
    isnan(range1[2])  &&  (range1[2] = maximum(v1[i]))
    i = findall(isfinite.(v2))
    isnan(range2[1])  &&  (range2[1] = minimum(v2[i]))
    isnan(range2[2])  &&  (range2[2] = maximum(v2[i]))

    i1 = findall(isfinite.(v1)  .&  (v1.>= range1[1])  .&  (v1.<= range1[2]))
    i2 = findall(isfinite.(v2)  .&  (v2.>= range2[1])  .&  (v2.<= range2[2]))
    (nbins1 > 0)  &&  (bs1 = (range1[2] - range1[1]) / nbins1)
    (nbins2 > 0)  &&  (bs2 = (range2[2] - range2[1]) / nbins2)
    if isfinite(bs1) &&  isfinite(bs2)
        hh = fit(Histogram, (v1[i1], v2[i2]), (range1[1]:bs1:range1[2], range2[1]:bs2:range2[2]), closed=:left)
    else
        hh = fit(Histogram, (v1[i1], v2[i2]), closed=:left)
    end
    x1 = collect(hh.edges[1])
    x1 = (x1[1:end-1] .+ x1[2:end]) ./ 2
    x2 = collect(hh.edges[2])
    x2 = (x2[1:end-1] .+ x2[2:end]) ./ 2

    binsize1 = x1[2] - x1[1]
    binsize2 = x2[2] - x2[1]
    return (loc1=x1, loc2=x2, counts=hh.weights, binsize1=binsize1)
end

function contourlines(args...; cntrparam="level auto 10")
    tmpfile = Base.Filesystem.tempname()
    sid = Symbol("j", Base.Libc.getpid())
    if !haskey(Gnuplot.state.sessions, sid)
        gp = gnuplot(sid, Gnuplot.state.cmd)
    end

    Gnuplot.exec(sid, "set term unknown")
    @gsp    sid "set contour base" "unset surface" :-
    @gsp :- sid "set cntrparam $cntrparam" :-
    @gsp :- sid "set table '$tmpfile'" :-
    @gsp :- sid args...
    Gnuplot.exec(sid, "unset table")
    Gnuplot.exec(sid, "reset")

    out = DataFrame()
    curlevel = NaN
    curx = Vector{Float64}()
    cury = Vector{Float64}()
    curid = 1
    elength(x, y) = sqrt.((x[2:end] .- x[1:end-1]).^2 .+
                          (y[2:end] .- y[1:end-1]).^2)
    function dump()
        ((length(curx) < 2)  ||  isnan(curlevel))  &&  return nothing
        tmp = DataFrame([Int, Float64, Vector{Float64}, Vector{Float64}, Vector{Float64}],
                        [:id, :level , :len           , :x             , :y])
        push!(tmp, (curid, curlevel, [0.; elength(curx, cury)], [curx...], [cury...]))
        append!(out, tmp)
        curid += 1
        # d = cumsum(elength(curx, cury))
        # i0 = findall(d .<= width); sort!(i0)
        # i1 = findall(d .>  width)
        # if (length(i0) > 0)  &&  (length(i1) > 0)
        #     rot1 = atan(cury[i0[end]]-cury[i0[1]], curx[i0[end]]-curx[i0[1]]) * 180 / pi
        #     rot = round(mod(rot1, 360))
        #     x = mean(curx[i0])
        #     y = mean(cury[i0])
        #     push!(outl, "set label " * string(length(outl)+1) * " '$curlevel' at $x, $y center front rotate by $rot")
        #     curx = curx[i1]
        #     cury = cury[i1]
        # end
        empty!(curx)
        empty!(cury)
    end

    for l in readlines(tmpfile)
        if length(strip(l)) == 0
            dump()
            continue
        end
        if !isnothing(findfirst("# Contour ", l))
            dump()
            curlevel = Meta.parse(strip(split(l, ':')[2]))
            continue
        end
        (l[1] == '#')  &&  continue

        n = Meta.parse.(split(l))
        @assert length(n) == 3
        push!(curx, n[1])
        push!(cury, n[2])
    end
    rm(tmpfile)

    if nrow(out) > 0
        levels = unique(out.level)
        sort!(levels)
        out[!, :levelcount] .= 0
        for i in 1:length(levels)
            j = findall(out.level .== levels[i])
            out[j, :levelcount] .= i
        end
    end
    return out
end

function boxxyerror(x, y; xmin=NaN, ymin=NaN, xmax=NaN, ymax=NaN, cartesian=false)
    @assert length(x) == length(y)
    @assert issorted(x)
    @assert issorted(y)
    xlow  = Vector{Float64}(undef, length(x))
    xhigh = Vector{Float64}(undef, length(x))
    ylow  = Vector{Float64}(undef, length(x))
    yhigh = Vector{Float64}(undef, length(x))
    for i in 2:length(x)-1
        xlow[i]  = (x[i-1] + x[i]) / 2
        ylow[i]  = (y[i-1] + y[i]) / 2
        xhigh[i] = (x[i+1] + x[i]) / 2
        yhigh[i] = (y[i+1] + y[i]) / 2
    end
    xlow[1]    = (isfinite(xmin)  ?  xmin  :  (x[1] - (x[2]-x[1])/2))
    ylow[1]    = (isfinite(ymin)  ?  ymin  :  (y[1] - (y[2]-y[1])/2))
    xlow[end]  = (x[end] - (x[end]-x[end-1])/2)
    ylow[end]  = (y[end] - (y[end]-y[end-1])/2)
    xhigh[1]   = (x[1] + (x[2]-x[1])/2)
    yhigh[1]   = (y[1] + (y[2]-y[1])/2)
    xhigh[end] = (isfinite(xmax)  ?  xmax  :  (x[end] + (x[end]-x[end-1])/2))
    yhigh[end] = (isfinite(ymax)  ?  ymax  :  (y[end] + (y[end]-y[end-1])/2))
    if !cartesian
        return (x, y, xlow, xhigh, ylow, yhigh)
    end
    n = length(x)
    i = repeat(1:n, outer=n)
    j = repeat(1:n, inner=n)
    return (x[i], y[j], xlow[i], xhigh[i], ylow[j], yhigh[j])
end

function histo2segments(in_x, counts)
    @assert length(in_x) == length(counts)
    x = Vector{Float64}()
    y = Vector{Float64}()
    push!(x, in_x[1])
    push!(y, counts[1])
    for i in 2:length(in_x)
        xx = (in_x[i-1] + in_x[i]) / 2.
        push!(x, xx)
        push!(y, counts[i-1])
        push!(x, xx)
        push!(y, counts[i])
    end
    push!(x, in_x[end])
    push!(y, counts[end])
    return (x, y)
end

end #module


text:
"""
Vibration analysis of simply-supported thick  (solid) plate

Free-vibration problem is solved for a homogeneous square plate,
simply-supported on the circumference.
This is the NAFEMS Benchmark, Test No. FV52.

The plate is discretized with shell elements.
The in-plane displacements are prevented.

The mass matrix is computed without consideration for the rotations. Hence, the
frequencies are underestimated.

The nonzero benchmark frequencies are (in hertz):
45.897  109.44  109.44  167.89
Note: these are only for the out of plane
mode shapes. The in-plane modes are not included.
"""
module FV52_vibration_examples

using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

function _execute_dsg_model(formul, n, visualize)
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 1.0*phun("m");
    L = 10.0*phun("m");

    # Report
    @info "Mesh: $n elements per side"

    tolerance = L/n/1000
    fens, fes = T3block(L,L,n,n);
    fens.xyz[:, 1] .-= L/2
    fens.xyz[:, 2] .-= L/2
    fens.xyz = xyz3(fens)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    associate = formul.associategeometry!
    stiffness = formul.stiffness
    mass = formul.mass

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # simple support
    l1 = selectnode(fens; box = Float64[L/2 L/2 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [3,4,6]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-L/2 -L/2 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [3,4,6]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [3,5,6]
        setebc!(dchi, l1, true, i)
    end
    l1 = selectnode(fens; box = Float64[-Inf Inf -L/2 -L/2 -Inf Inf], inflate = tolerance)
    for i in [3,5,6]
        setebc!(dchi, l1, true, i)
    end
    # in-plane, rotations
    l1 = collect(1:count(fens))
    for i in [1, 2, 6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, dchi);

    # Solve
    OmegaShift = 0.1*2*pi
    neigvs = 4
    d, v, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    d[:] = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    @show fs

    # Visualization
    U = v[:, 4]
    scattersysvec!(dchi, (L/4)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
        #plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
end

function test_convergence()
    formul = FEMMShellT3FFModule
    @info "FV12 free vibration, formulation=$(formul)"
    @info "(out of plane only): 45.897  109.44  109.44  167.89  "
    for n in [2, 4, 8, 16, 32, 64]
        _execute_dsg_model(formul, n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
using Scallop.Walkforward

@assert 12 == xwalkforward(6)


text:
using SafeTestsets

@safetestset "Earthmover distance" begin include("earthmover.jl") end


text:
function TemperatureSensor(; name)
    @named hp = HeatPort()
    @variables T(t)

    eqs = [
        T ~ hp.T
        hp.Q_flow ~ 0
    ]
    ODESystem(eqs, t, [T], [], systems=[hp], name=name)
end

function RelativeTemperatureSensor(; name)
    @named hp1 = HeatPort()
    @named hp2 = HeatPort()
    @variables T(t)

    eqs = [
        T ~ hp1.T - hp2.T
        hp1.Q_flow ~ 0
        hp2.Q_flow ~ 0
    ]
    ODESystem(eqs, t, [T], [], systems=[hp1, hp2], name=name)
end

function HeatFlowSensor(; name)
    @named hp1 = HeatPort()
    @named hp2 = HeatPort()
    @variables Q_flow(t)

    eqs = [
        hp1.T ~ hp2.T
        hp1.Q_flow + hp2.Q_flow ~ 0
        Q_flow ~ hp1.Q_flow
    ]
    ODESystem(eqs, t, [Q_flow], [], systems=[hp1, hp2], name=name)
end


text:
using StaticArrays, LinearAlgebra

function bands(el::ElecHam, kpt::SVector)
   nband = el.basic_data[:num_wann]
   Hk = zeros(ComplexF64, nband, nband)
   _band!(Hk, el, kpt)
end

function bands(el::ElecHam, kpts::AVec{<:SVector})
   nband = el.basic_data[:num_wann]
   Hk = zeros(ComplexF64, nband, nband)

   eigs = zeros(Float64, nband, length(kpts))

   for (i, k) in enumerate(kpts)
      eigs[:,i] = _band!(Hk, el, k)
   end
   return eigs .* Rydberg2eV
end

function bands(el::ElecHam, kpts::AbstractArray{T,2}) where T
   size(kpts, 1) != 3 && throw(error("1-dimension of kpts should be 3"))
   bands(el, vec(reinterpret(SVector{3,T}, kpts)))
end

function _band!(Hk::AbstractMatrix, el::ElecHam, kpt::SVector)
   hamiltonian!(Hk, el, kpt)
   eigvals( Hermitian(Hk) ) .* Rydberg2eV
end

exp_ikr(rvecs::AbstractVector{<:SVector}, kpt::SVector) = cis.(Ref(kpt') .* rvecs)

_ham_elem(hopping, rvec_idx, e_ikr) = sum( e_ikr[val] * hopping[n] for (n, val) in enumerate(rvec_idx) )

function hamiltonian!(H::AbstractMatrix, el::ElecHam, kpt::SVector)
   nband = LinearAlgebra.checksquare(H)
   nband == el.basic_data[:num_wann] || throw(DimensionMismatch("H is not square: dimensions are $(size(H))"))

   #compute exp(i 2\pi*K*R)
   e_ikr = exp_ikr(el.ws_rvecs, 2π*kpt)

   # compute upper triangle of ham
   k = 0
   @inbounds for j = 1:nband, i = 1:j
      k += 1
      H[i,j] = _ham_elem(el.hopping[k], el.rvec_idx[k], e_ikr)
   end
end

#function bands(ph::LatticeIFC, qpt::SVector)
#   nsize = 3*ph.basic_data[:nat]
#   dynmat = zeros(ComplexF64, nsize, nsize)
#   _dynamic_matrix!(dynmat, ph, qpt)
#end

#function _dynamic_matrix!(dynmat::AbstractMatrix, ph::LatticeIFC, qpt::SVector) where T
#   n, nsize = size(dynmat)
#   (n == nsize && nsize == 3*ph.basic_data[:nat]) ||
#      throw(DimensionMismatch("dynmat dimensions are $(size(dynmat))"))
#
#   #compute exp(i 2\pi*K*R)
#   e_ikr = exp_ikr(ph.ws_rvecs, 2π*qpt)
#
#   k = 0
#   R = CartesianIndices((-2:0, -2:0))
#   mass = ph.basic_data[:mass]
#   @inbounds for j = 1:ph.basic_data[:nat], i = 1:j
#      k += 1
#      mfactor = 1.0 / sqrt( mass[i]*mass[j] )
#      dynmat[R .+ CartesianIndex(3*i,3*j)] .=
#         _dyn_mat(ph.forceconst[k], ph.rvec_idx[k], e_ikr) .* mfactor
#   end
#end
#
#
#function _dyn_mat(ifc, rvec_idx, e_ikr)
#   dmat = zeros(eltype(e_ikr), 3, 3)
#   for (n, val) in enumerate(rvec_idx)
#      dmat .+= e_ikr[val] .* view(ifc, :, :, n)
#   end
#   return dmat
#end
#


text:
export install_adept, install_blas

function install_blas(blas_binary)
    if Sys.iswindows()
        if isfile(joinpath(ADCME.LIBDIR, "openblas.lib"))
            return
        end
        if blas_binary
            @info "Downloading prebuilt blas from Github. If you encounter any problem with openblas when using adept, run `install_adept(blas_binary=false)` to compile from source"
            download("https://github.com/kailaix/tensorflow-1.15-include/releases/download/v0.1.0/openblas.lib", joinpath(ADCME.LIBDIR, "openblas.lib"))
            return
        end
        @info "You are building openblas from source on Windows, and this process may take a long time.
Alternatively, you can place your precompiled binary to $(joinpath(ADCME.LIBDIR, "openblas.lib"))"
        PWD = pwd()
        download("https://github.com/xianyi/OpenBLAS/archive/v0.3.9.zip", joinpath(ADCME.LIBDIR, "OpenBlas.zip"))
        cd(ADCME.LIBDIR)
        run(`cmd /c unzip OpenBLAS.zip`)
        rm("OpenBlas", force=true, recursive=true)
        run(`cmd /c ren OpenBlas-0.3.9 OpenBlas`)
        rm("OpenBlas.zip")
        cd("OpenBlas")
        mkdir("build")
        cd("build")
        ADCME.cmake(CMAKE_ARGS="-DCMAKE_Fortran_COMPILER=flang -DBUILD_WITHOUT_LAPACK=no -DNOFORTRAN=0 -DDYNAMIC_ARCH=ON")
        ADCME.make()
        cd("../build/lib/Release")
        mv("openblas.lib", joinpath(ADCME.LIBDIR, "openblas.lib"))
        cd(PWD)
    else
        required_file = Sys.isapple() ? ".dylib" : ".so"
        required_file = joinpath(ADCME.LIBDIR, "libopenblas")*required_file
        if !isfile(required_file)
            files = readdir(ADCME.LIBDIR)
            files = filter(x->!isnothing(x), match.(r"(libopenblas\S*.dylib)", files))[1]
            target = joinpath(ADCME.LIBDIR, files[1])
            symlink(target, required_file)
            @info "Symlink $(required_file) --> $(files[1])"
        end
    end
end

"""
    install_adept(force::Bool=false)

Install adept-2 library: https://github.com/rjhogan/Adept-2
"""
function install_adept(force::Bool=false; blas_binary::Bool = true)
    PWD = pwd()
    cd(ADCME.LIBDIR)
    if force
        @info "Removing Adept-2 by force..."
        rm("Adept-2", force=true, recursive=true)
    end
    if !isdir("Adept-2")
        LibGit2.clone("https://github.com/ADCMEMarket/Adept-2", "Adept-2")
    end
    cd("Adept-2/adept")
    install_blas(blas_binary)
    try
        if (!isfile("$(LIBDIR)/libadept.so") && !isfile("$(LIBDIR)/libadept.dylib") && !isfile("$(LIBDIR)/adept.lib")) || force
            @info """Copy "$(@__DIR__)/../deps/AdeptCMakeLists.txt" to "$(joinpath(pwd(), "CMakeLists.txt"))" ... """
            cp("$(@__DIR__)/../deps/AdeptCMakeLists.txt", "./CMakeLists.txt", force=true)
            @info """Remove $(joinpath(pwd(), "build")) ... """
            rm("build", force=true, recursive=true)
            @info "Make $(joinpath(pwd(), "build")) ... "
            mkdir("build")
            @info "Change directory into $(joinpath(pwd(), "build")) ... "
            cd("build")
            @info "Cmake ... "
            ADCME.cmake()
            @info "Make ... "
            ADCME.make()
        end
        printstyled("""
∘ Add the following lines to CMakeLists.txt

include_directories(\${LIBDIR}/Adept-2/include)
find_library(ADEPT_LIB_FILE adept HINTS \${LIBDIR})
find_library(LIBOPENBLAS openblas HINTS \${LIBDIR})
message("ADEPT_LIB_FILE=\${ADEPT_LIB_FILE}")
message("LIBOPENBLAS=\${LIBOPENBLAS}")

∘ Add `\${ADEPT_LIB_FILE}` and `\${LIBOPENBLAS}` to `target_link_libraries`
""", color=:green)
    catch
        printstyled("Compliation failed\n", color=:red)
    finally
        cd(PWD)
    end
end

CONDA_ROOT = abspath(joinpath(LIBDIR, ".."))
function maybe_download(URL)
    filename = splitdir(URL)[2]
    file = joinpath(PREFIXDIR, filename)
    if ispath(file)
        return file
    end
    download(URL, file)
    return file
end

function maybe_tar_gz(file)
    filename = file[1:end-7]
    if isdir(filename)
        return
    else
        run(`tar zxvf $file`)
    end
end


text:
using HDF5,JLD
using GZip
using Lumberjack

import HDF5.dump

function dump(io::IO, x::Union{HDF5File, HDF5Group}, n::Int, indent, display_num=10)
    println(typeof(x), " len ", length(x))
    if n > 0
        i = 1
        for k in names(x)
            print(io, indent, "  ", k, ": ")
            v = o_open(x, k)
            dump(io, v, n - 1, string(indent, "  "))
            close(v)
            if i > display_num
                println(io, indent, "  ...")
                break
            end
            i += 1
        end
    end
end

function store_info_dict(genomic_store_path)
  d=Dict{String, Array{Symbol} }()
  fid=h5open(genomic_store_path,"r")
  for seq in names(fid)
     track_array=names( fid[seq] )
     d[seq]=track_array
  end
  return d
end

function tracks(genomic_store_path)
  # query tracks
end

function store_info_str(genomic_store_path)
  buf=IOBuffer()
  fid=h5open(genomic_store_path,"r")
  dump(buf,fid,10,"\n",20)
  close(fid)
  return takebuf_string(buf)
end

function get_chr_sizes_dict(chrom_sizes_path)
    chrom_sizes=readdlm(chrom_sizes_path)
    chrom_sizes_dict = Dict{ASCIIString,Int64}()
    size(chrom_sizes,1)
    for i=1:size(chrom_sizes,1)
        chrom_sizes_dict[ chrom_sizes[i,1] ] = chrom_sizes[i,2]
    end
    return chrom_sizes_dict
end

function read_gzip_file(filename)
   gzio=gzopen(filename,"r")
   line_num =0
   lines=ASCIIString[]
   for line in eachline(gzio)
        line=chomp(line)
        push!(lines,line)
        line_num +=1
        if line_num % 1000000 == 0
            Lumberjack.info("full on speed like road runner $line_num already!")
        end
    end
    #pop!(lines) # clears out the final row which due to a bug in GZlib gives and empty string
    return lines
end

function memory_read_file(filename)
    io = open(filename)
    Lumberjack.info("reading all")
    file=readall(io)
    Lumberjack.info("finished reading all")
    Lumberjack.info("split line")
    lines=split(file,'\n')
    Lumberjack.info("done split")
    close(io)
    return lines
end

# This file has first three columns chr start strand or something else
# but we ignore them.
function memory_read_and_integer_parse_point_file(filename;strand_filter_char=nothing)
    io = open(filename)
    Lumberjack.info("reading all")
    file=readall(io)
    Lumberjack.info("finished reading all")
    Lumberjack.info("split line")
    lines=split(file,'\n')
    Lumberjack.info("done split")

    num_lines=length(lines)-1
    seq_ids = ASCIIString[] #fill("",num_lines)
    starts=   Int64[] #fill( 0, num_lines)
    num_kept = 0
    Lumberjack.info("parse and assign to array")
    #parseint takes for ever, so annoying about this.
    for idx=1:num_lines
        (seq_id,start,strand)= split(strip(lines[idx]),'\t')
        if strand == strand_filter_char
            continue
        end
        push!(seq_ids,seq_id)
        push!(starts,int64(start))
        if idx % 1000000 == 0
            println(idx)
        end
        num_kept+=1
    end
    Lumberjack.info("finished parsing integer and assigning to array ( number kept after filter if applied $num_kept). first seq_id is $(seq_ids[1]). Last seq_id $(seq_ids[end])")
    return (seq_ids,starts)
end

# -
# - This file has first three columns chr start stop
# -  it could have others but we ignore them.
function memory_read_and_integer_parse_interval_file(filename; only_start=false,strand_filter_char=nothing)
    io = open(filename)
    Lumberjack.info("reading all")
    file=readall(io)
    Lumberjack.info("finished reading all")
    Lumberjack.info("split line")
    lines=split(file,'\n')
    Lumberjack.info("done split")
    # now we know what size thigns should be
    # get size and two columns and then
    # lines -1 as it splits the final row
    num_lines=length(lines)-1
    seq_ids = fill("",num_lines)
    starts=fill( 0, num_lines)
    if !only_start
       stops=fill( 0,  num_lines )
    end
    Lumberjack.info("parse and assign to array")
    #parseint takes for ever, so annoying about this.
    for idx=1:num_lines
        if only_start
            (seq_id,start)= split(lines[idx],'\t')
        else
            (seq_id,start,stop)= split(lines[idx],'\t')
        end
        seq_ids[idx]=seq_id
        starts[idx]=int64(start)
        if !only_start
            stops[idx]=int64(stop)
        end

        if idx % 1000000 == 0
            println(idx)
        end
    end
    Lumberjack.info("finished parsing integer and assigning to array")
    if only_start
        return (seq_ids,starts)
    end
    return  (seq_ids,starts,stops)
end

# the integer_parse part of this method name needs to be renamed
function memory_read_and_integer_parse_methpipe_cg_bed(filename;gzip=false)
    if gzip
        lines=read_gzip_file(filename)
    else
        lines = memory_read_file(filename)
    end
    num_lines=length(lines)-1
    seq_ids = fill("",num_lines)
    starts=fill( 0, num_lines)
    stops=fill( 0,  num_lines )
    scores=fill(float32(0.0), num_lines )
    Lumberjack.info("parse and assign to array")
    #parseint takes for ever, and that is what is so annoying about this.
    for idx=1:num_lines
        (seq_id,start,strand,context,score,coverage)= split(lines[idx],'\t')
        seq_ids[idx] = seq_id
        starts[idx]  = int64(start)
        stops[idx]   = int64(start) +1
        scores[idx]  = float32(score)
        if idx % 1000000 == 0
            println(idx)
        end
    end
    Lumberjack.info("finished parsing integer and assigning to array")
    return (seq_ids,starts,stops,scores)
end

function memory_read_and_integer_parse_methpipe_cg_bed_coverage(filename;gzip=false)
    if gzip
        lines=read_gzip_file(filename)
    else
        lines = memory_read_file(filename)
    end
    num_lines=length(lines)-1
    seq_ids = fill("",num_lines)
    starts=fill( 0, num_lines)
    stops=fill( 0,  num_lines )
    coverages=fill(int32(0), num_lines )
    Lumberjack.info("parse and assign to array")
    #parseint takes for ever, and that is what is so annoying about this.
    for idx=1:num_lines
        (seq_id,start,strand,context,score,coverage)= split(lines[idx],'\t')
        try
            seq_ids[idx] = seq_id
            starts[idx]  = int64(start)
            stops[idx]   = int64(start) +1
            coverages[idx]  = int32(coverage)
            if idx % 1000000 == 0
                println(idx)
            end
       catch y
          println("processing this line($idx):\n")
          println(lines[idx])
          println("parsed to:")
          println("$seq_id $start $strand $context $score $coverage")
          error(y)
       end
    end
    Lumberjack.info("finished parsing integer and assigning to array")
    return (seq_ids,starts,stops,coverages)
end

function memory_read_and_integer_parse_bedgraph(filename;gzip=false)
    if gzip
        lines=read_gzip_file(filename)
    else
        lines = memory_read_file(filename)
    end
    #io = open(filename)
    #Lumberjack.info("reading all")
    #file=readall(io)
    #Lumberjack.info("finished reading all")
    #Lumberjack.info("split line")
    #lines=split(file,'\n')
    #Lumberjack.info("done split")
    # now we know what size thigns should be
    # get size and two columns and then
    # lines -1 as it splits the final row
    num_lines=length(lines)-1
    seq_ids = fill("",num_lines)
    starts=fill( 0, num_lines)
    stops=fill( 0,  num_lines )
    scores=fill(float32(0.0), num_lines )
    Lumberjack.info("parse and assign to array")
    #parseint takes for ever, and that is what is so annoying about this.
    for idx=1:num_lines
        (seq_id,start,stop,score)= split(lines[idx],'\t')
        seq_ids[idx] = seq_id
        starts[idx]  = int64(start)
        stops[idx]   = int64(stop)
        scores[idx]  = float32(score)
        if idx % 1000000 == 0
            println(idx)
        end
    end
    Lumberjack.info("finished parsing integer and assigning to array")
    return (seq_ids,starts,stops,scores)
end

# writes a track - that must have a sequence id
function write_track(file_handle,seq_id,track_id,values)
    Lumberjack.info("Writing $seq_id")
    write(file_handle, bytestring("$seq_id/$track_id"), values)
    Lumberjack.info("Finished writing $seq_id")
end

#- We have different tracks for chromosomes so the files do not have
#- to be consistently sorted by chromosome

function save_bedgraph_track(genomic_store_path,input_file,track_id,chr_sizes_path;
                             start_coord_shift=0,
                             OUT_OF_RANGE_VAL=0.0,
                             gzip=false,
                             val_type="float32",
                             bedtype="bedgraph")

    return save_bed_track(genomic_store_path,input_file,track_id,chr_sizes_path,
                             start_coord_shift=start_coord_shift,
                             OUT_OF_RANGE_VAL=OUT_OF_RANGE_VAL,
                             gzip=gzip,
                             val_type=val_type,
                             bedtype="bedgraph")
end

function save_bed_track(genomic_store_path,input_file,track_id,chr_sizes_path;
                             start_coord_shift=0,
                             OUT_OF_RANGE_VAL=0.0,
                             gzip=false,
                             val_type="float32",
                             bedtype="bedgraph")
    chr_sizes_dict=get_chr_sizes_dict(chr_sizes_path)
    #--- sort in temporary dir
    #--- load
    if bedtype == "bedgraph"
        (seq_ids,starts,stops,scores)=memory_read_and_integer_parse_bedgraph(input_file,gzip=gzip)
    elseif bedtype == "methpipe_cpg_bed"
        (seq_ids,starts,stops,scores)=memory_read_and_integer_parse_methpipe_cg_bed(input_file,gzip=gzip)
    elseif bedtype == "methpipe_cpg_bed_coverage"
        (seq_ids,starts,stops,scores)=memory_read_and_integer_parse_methpipe_cg_bed_coverage(input_file,gzip=gzip)
    else
         Lumberjack.error("Invalid bedtype specified $bedtype")
    end
    bedgraph_file_length = length(seq_ids)
    Lumberjack.info("Read bedgraph file, total length: $bedgraph_file_length")
    if !isfile(genomic_store_path)
        fid=jldopen(genomic_store_path,"w",compress=true)
    else
        fid=jldopen(genomic_store_path,"r+",compress=true)
    end
    seq_id="FIRSTROW"
    if seq_ids[1] == "FIRSTROW"
       seq_id="FIRSTROW1"
    end
    seq_vals=nothing
    #seen_seq_id = Dict()
    #curr_seq_id=""
    for i=1:bedgraph_file_length
        if seq_id != seq_ids[i]
            if i != 1
                write_track(fid,seq_id,track_id,seq_vals)
            end
            seq_id   = seq_ids[i]
            Lumberjack.info("Next sequence name: $seq_id")

            # -- check the id is valid --

            if !haskey(chr_sizes_dict,seq_id)
                Lumberjack.error("Invalid ID: ( $seq_id ), this file only saved up to line $(i-1)")
            end

            # -- TO DO check we have sorted data by chr
            #if curr_seq_id != seq_id
            #    seen_seq_id[curr_seq_id]="DONE"
            #    curr_seq_id = seq_id
            #    seen_seq_id[curr_seq_id]="DOING"
            #end

            seq_len  = chr_sizes_dict[seq_id]
            seq_vals = fill(OUT_OF_RANGE_VAL,seq_len)
            Lumberjack.info("Next track to be $seq_len in length")
        end

        #--- will slow things down a little but check that the id is sorted, bare in mind a chr change.
        if val_type == "float32"
            try
                seq_vals[ starts[i] + start_coord_shift ] = float32( scores[i] )
            catch y
                println("processing this line:\n")
                curr_row = join([ seq_ids[i], starts[i],stops[i], scores[i] ],"\t" )
                println(curr_row)
                println("shift: $start_coord_shift")
                error(y)
            end
        end
        if val_type == "int32"
            seq_vals[ starts[i] + start_coord_shift ] = int32( scores[i] )
        end
        if val_type == "int16"
            seq_vals[ starts[i] + start_coord_shift ] = int16( scores[i] )
        end
        if val_type == "int8"
            seq_vals[ starts[i] + start_coord_shift ] = int8( scores[i] )
        end
    end
    # write the final track
    write_track(fid,seq_id,track_id,seq_vals)
    Lumberjack.info("Written track $seq_id")
    close(fid)
end

function get_track_interval( genomic_store_path, id, seq_id, start_pos, stop_pos )
    fid=h5open(genomic_store_path,"r")
    if has(fid,"$seq_id/$id")
        if start_pos == 0 && stop_pos == 0
            track_slice=h5read(genomic_store_path,bytestring("$seq_id/$id"))
        else
            track_slice=h5read(genomic_store_path,bytestring("$seq_id/$id"),(start_pos:stop_pos,))
        end
        close(fid)
        return track_slice
    end
    close(fid)
    return []
end

# this assumes each interval needs to be recorded
# as existing of not existing - use 1 or 0
# this function just records the start position at the moment.
function save_point_start_track(genomic_store_path,point_file,track_id,chr_sizes_path;start_coord_shift=0,strand_filter_char=nothing)
    chr_sizes_dict=get_chr_sizes_dict(chr_sizes_path)
    #--- sort in temporary dir
    #--- load
    (seq_ids,starts)=memory_read_and_integer_parse_point_file(point_file, strand_filter_char=strand_filter_char)
    bedgraph_file_length = length(seq_ids)
    Lumberjack.info("Read bedgraph file, total length: $bedgraph_file_length")

    if !isfile(genomic_store_path)
        fid=jldopen(genomic_store_path,"w",compress=true)
    else
        fid=jldopen(genomic_store_path,"r+",compress=true)
    end
    seq_id="FIRSTROW"
    if seq_ids[1] == "FIRSTROW"
       seq_id="FIRSTROW1"
    end
    seq_vals=nothing
    for i=1:bedgraph_file_length
        if seq_id != seq_ids[i]
            if i != 1
                write_track(fid,seq_id,track_id,seq_vals)
            end
            seq_id   = seq_ids[i]
            Lumberjack.info("Next sequence name: $seq_id")

            # -- check the id is valid --
            if !haskey(chr_sizes_dict,seq_id)
                Lumberjack.error("Invalid ID: ( $seq_id ), this file only saved up to line $(i-1)")
            end

            # -- check we have sorted data by chr
            seq_len  = chr_sizes_dict[seq_id]
            seq_vals = fill(int8(0),seq_len)
            Lumberjack.info("Next track to be $seq_len in length")
        end
        #---TODO will slow things down a little but check that the id is sorted baring in mind a chr change.
        seq_vals[ starts[i] + start_coord_shift ] = 1
    end
    # write the final track
    write_track(fid,seq_id,track_id,seq_vals)
    close(fid)
end

# this assumes each interval needs to be recorded
# as existing or not existing - use 1 or 0
# this function just records the start position at the moment.
# interval_label_type may be either "binary" uses 1 or 0 to indicate an interval covers a basepair
# or multiclass_int8 which allows an 256 possible integer values (-128 .. 127)
function save_interval_track(genomic_store_path,interval_file,track_id,chr_sizes_path;
                             start_coord_shift=0,strand_filter_char=nothing, interval_label_type="binary")

    chr_sizes_dict=get_chr_sizes_dict(chr_sizes_path)
    #--- sort in temporary dir
    #--- load
    if interval_label_type == "binary"
        (seq_ids,starts,stops)=memory_read_and_integer_parse_interval_file(interval_file, strand_filter_char=strand_filter_char)
    elseif interval_label_type == "multiclass_int8"
        (seq_ids,starts,stops,scores)=memory_read_and_integer_parse_bedgraph(interval_file)
    else
        Lumberjack.error("Interval Label Type not supported")
    end

    bedgraph_file_length = length(seq_ids)
    Lumberjack.info("Read bedgraph file, total length: $bedgraph_file_length")

    if !isfile(genomic_store_path)
        fid=jldopen(genomic_store_path,"w",compress=true)
    else
        fid=jldopen(genomic_store_path,"r+",compress=true)
    end
    seq_id="FIRSTROW"
    if seq_ids[1] == "FIRSTROW"
       seq_id="FIRSTROW1"
    end
    seq_vals=nothing
    for i=1:bedgraph_file_length
        if seq_id != seq_ids[i]
            if i != 1
                write_track(fid,seq_id,track_id,seq_vals)
            end
            seq_id   = seq_ids[i]
            Lumberjack.info("Next sequence name: $seq_id")

            # -- check the id is valid --
            if !haskey(chr_sizes_dict,seq_id)
                Lumberjack.error("Invalid ID: ( $seq_id ), this file only saved up to line $(i-1)")
            end

            # -- check we have sorted data by chr
            seq_len  = chr_sizes_dict[seq_id]
            seq_vals = fill(int8(0),seq_len)
            Lumberjack.info("Next track to be $seq_len in length")
        end
        #---TODO will slow things down a little but check that the id is sorted baring in mind a chr change.
        starting_pos=starts[i] + start_coord_shift
        stopping_pos=stops[i]  + start_coord_shift
        if interval_label_type == "binary"
            score=1
        elseif interval_label_type == "multiclass_int8"
            score = scores[i]
        end
        for pos=starting_pos:stopping_pos
            seq_vals[ pos ] = score
        end
    end
    # write the final track
    write_track(fid,seq_id,track_id,seq_vals)
    close(fid)
end