text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V,(VV,EV,FV) = Lar.simplex(2, true)
triangle = (V,EV,FV)
model = Lar.Struct([ triangle, Lar.t(.15,.15), Lar.s(.5,.5), triangle ])
V,EV,FV = Lar.struct2lar(model)
VV = [[k] for k=1:size(V,2)]
GL.VIEW(GL.numbering(.4)((V,[VV, EV, FV]),GL.COLORS[1],0.2));

triangles = Lar.triangulate2d(V, EV)
GL.VIEW([ GL.GLGrid(V,triangles) ])

# generate edges
ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
# remove duplicated edges from triangulation
nodupev = collect(Set(ev))
GL.VIEW(GL.numbering(.4)((V,[VV, nodupev])));


text:
include("test_common.jl")
using Random

@testset "external" begin
    new_context()

    x = Ref(123)
    o = @new node"Object"()
    o.x = NodeExternal(x)
    @test o.x ≡ x
    ox = get(o, "x", result=RESULT_NODE)
    @test sprint(show, ox) == "NodeExternal: Base.RefValue{Int64}(123)"
    @test node_value(ox) ≡ ox

    delete_context()
end

ext: jl


text:
""" `make_manifest(data_path::String, r::Regex; [opts] )`
Returns a dataframe containing a list of files in datapath with file names that match r.

# Optional arguements
- verbose = true
"""
function make_manifest(data_path::String, r::Regex; verbose::Bool = true)
   if !isdir(data_path)
      @error "Can't access data directory ", data_path, "."
   end
   if verbose  println("# Creating manifest of files to process.")    end

   dir_filelist = readdir(data_path,join=true)
   idx_spectra = map(fn->occursin(r, last(split(fn,'/')) ),dir_filelist)
   spectra_filelist = dir_filelist[idx_spectra]

   df_files = DataFrame(:Filename => dir_filelist[idx_spectra])

   if size(df_files,1) < 1
      @error("Did not find any files in ", data_path ,".")
   end

   if match(r"(\d{8})[T_](\d{6})", last(split(df_files.Filename[1],'/')) ) != nothing
      results = map(fn->match(r"(\d{8})[T_](\d{6})", last(split(fn,'/')) ), df_files.Filename)
      @assert all(map(r->length(r.captures),results).==2)
      df_files.date_str = map(r->String(r.captures[1]),results)
      df_files.time_str = map(r->String(r.captures[2]),results)
   end

   return df_files
end

""" `make_manifest(data_path::String, target_subdir::String, Inst::Module; [opts] )`
Returns a dataframe containing a list of files to be read and some metadata (e.g., observation times)

# Optional arguements
- verbose = true
"""
function make_manifest(data_path::String, Inst::Module; verbose::Bool = true)
   #=
   if Inst == EXPRES
      #@assert isdefined(Main,:expres_data_path)
      data_path = joinpath(expres_data_path,target_subdir)
   elseif Inst == NEID
      #  @assert isdefined(Main,:solar_data_path)
      data_path = joinpath(solar_data_path,target_subdir)
   else
      @error "Specified instrument didn't match a module name, ", string(Inst), "."
   end
   =#
   if !isdir(data_path)
      @error "Can't access data directory ", data_path, "."
   end
   if verbose  println("# Creating manifest of files to process.")    end
   df_files = Inst.make_manifest(data_path)
   if size(df_files,1) < 1
      @error("Did not find any files in ", data_path ,".")
   end
   #=
   df_files = df_files |> @filter( _.target == target_fits ) |> DataFrame
   if size(df_files,1) < 1
      @error("Did not find any files with target field matching ", target_fits,".")
   end
   if size(df_files,1) > max_spectra_to_use
      df_files = df_files |> @take(max_spectra_to_use) |> DataFrame
   end
   =#
   return df_files
end

""" `make_manifest(data_path::String, target_subdir::String, Inst::Module; [opts] )`
Returns a dataframe containing a list of files to be read and some metadata (e.g., observation times)

# Optional arguements
- verbose = true
"""
function make_manifest(data_path::String, target_subdir::String, Inst::Module; verbose::Bool = true)
   return make_manifest(joinpath(data_path,target_subdir), Inst, verbose=verbose)
end

"""Read manifest containing `filename`, `bjd`, `target`, and optionally additional metadata from CSV file. """
function read_manifest(fn::String)
    df = CSV.read(fn,DataFrame,threaded=false)
    @assert hasproperty(df,:filename)
    @assert hasproperty(df,:bjd)
    @assert hasproperty(df,:target)
    @assert size(df,1) >= 1
    return df
end

"""Write manifest containing `filename`, `bjd`, `target`, and optionally additional metadata from CSV file. """
function write_manifest(fn::String, df)
    @assert hasproperty(df,:filename)
    @assert hasproperty(df,:bjd)
    @assert hasproperty(df,:target)
    @assert size(df,1) >= 1
    CSV.write(fn,df)
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using Stingray
using Documenter

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

makedocs(;
    modules=[Stingray],
    authors="Aman Pandey",
    repo="https://github.com/matteobachetti/Stingray.jl/blob/{commit}{path}#{line}",
    sitename="Stingray.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/matteobachetti/Stingray.jl",
        edit_link="master",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="https://github.com/matteobachetti/Stingray.jl",
    devbranch="main",
)


text:
"""
BitSet-based container that could be used to e.g. store the ids of completed
jobs by `MultithreadEvaluator`.
Assumes there is `max_seq_el` element, so that all *1:max_seq_el* elements
are in the set.
"""
mutable struct SlidingBitset
    max_seq_el::Int # the set contains all elements from 1 to max_seq_el
    max_el::Int     # max element contained in the set
    els::BitSet     # elements greater than max_seq_el that are in the set

    SlidingBitset(; max_seq_el::Integer=0) = new(max_seq_el, max_seq_el, BitSet())
end

# how much the BitSet offset could differ from max_seq_el/64
# might need to be increased in massively parallel environments
# (more than SLIDE_OFFSET*64 threads)
const SLIDE_OFFSET = 8

function Base.push!(set::SlidingBitset, el::Integer)
    if el > set.max_el
        set.max_el = el
    end
    if el == set.max_seq_el+1
        # update the next sequential el
        set.max_seq_el = el
        # see if max_seq_el could be further advanced using els
        while (set.max_el > set.max_seq_el) && in(set.els, set.max_seq_el+1)
            set.max_seq_el += 1
        end
        if !isempty(set.els) && Base._div64(set.max_seq_el) > set.els.offset + SLIDE_OFFSET
            # slide els bitset
            len = length(set.els.bits)
            if len > SLIDE_OFFSET
                copyto!(set.els.bits, 1, set.els.bits, SLIDE_OFFSET+1, len-SLIDE_OFFSET)
                fill!(view(set.els.bits, max(1, len-SLIDE_OFFSET+1):len), 0)
            else
                fill!(set.els.bits, 0)
            end
            set.els.offset += SLIDE_OFFSET
        end
    else
        push!(set.els, el)
    end
    return set
end

Base.in(set::SlidingBitset, el::Integer) =
    el <= set.max_seq_el || in(set.els, el)


text:
"""
    @enum PrintType

`PrintType` to tell lower level printing some useful context.
Currently only supports `Unknown` and `StructField`.
"""
@enum PrintType begin
    Unknown
    StructField
end

mutable struct PrintState
    type::PrintType
    typeinfo
    noindent_in_first_line::Bool
    level::Int
    # the offset that should be applied
    # to whoever cares, e.g for the field
    # values
    offset::Int
end

PrintState() = PrintState(Unknown, Any, false, 0, 0)

"""
    GarishIO{IO_t <: IO} <: IO

`GarishIO` contains the pretty printing preference and states.

# Members

- `bland_io::IO_t`: the original io.
- `indent::Int`: indentation size.
- `compact::Bool`: whether the printing should be compact.
- `displaysize::Tuple{Int, Int}`: the terminal displaysize.
- `show_indent`: print the indentation hint or not.
- `color`: color preference, either `ColorPreference` or `nothing` for no color.
- `state`: the state of the printer, see [`PrintState`](@ref).
"""
struct GarishIO{IO_t <: IO} <: Base.AbstractPipe
    # the bland io we want look nice
    bland_io::IO_t
    indent::Int
    compact::Bool
    limit::Bool
    displaysize::Tuple{Int, Int}
    show_indent::Bool
    # option type
    include_defaults::Bool
    # use nothing for no color print
    color::Union{Nothing, ColorPreference}
    state::PrintState
end

"""
    GarishIO(io::IO; kw...)

See [`pprint`](@ref) for available keywords.
"""
function GarishIO(io::IO;
        indent::Int=2,
        compact::Bool=get(io, :compact, false),
        limit::Bool=get(io, :limit, false),
        displaysize::Tuple{Int, Int}=displaysize(io),
        color::Bool=get(io, :color, true),
        # indent is similar to color
        show_indent::Bool=get(io, :color, true),
        include_defaults::Bool=get(io, :include_defaults, false),
        kw...
    )

    if color
        color_prefs = ColorPreference(;kw...)
    else
        color_prefs = nothing
    end
    return GarishIO(
        io, indent,
        compact, limit,
        displaysize,
        show_indent,
        include_defaults,
        color_prefs,
        PrintState()
    )
end

function GarishIO(io::GarishIO;
    indent::Int=io.indent,
    compact::Bool=io.compact,
    limit::Bool=io.limit,
    displaysize::Tuple{Int, Int}=io.displaysize,
    show_indent::Bool=io.show_indent,
    include_defaults::Bool=get(io, :include_defaults, false),
    color=io.color,
    state=io.state,
    kw...)

    if isempty(kw)
        color_prefs = color
    else
        color_prefs = ColorPreference(;kw...)
    end

    return GarishIO(
        io, indent, compact, limit,
        displaysize, show_indent,
        include_defaults,
        color_prefs,
        state
    )
end

"""
    GarishIO(io::IO, garish_io::GarishIO; kw...)

Create a new similar `GarishIO` with new bland IO object `io`
based on an existing garish io preference. The preference can
be overloaded by `kw`. See [`pprint`](@ref) for the available
keyword arguments.
"""
function GarishIO(io::IO, garish_io::GarishIO;
        indent::Int=garish_io.indent,
        compact::Bool=garish_io.compact,
        limit::Bool=garish_io.limit,
        displaysize::Tuple{Int, Int}=garish_io.displaysize,
        show_indent::Bool=garish_io.show_indent,
        include_defaults::Bool=get(io, :include_defaults, false),
        color_prefs=garish_io.color,
        state=garish_io.state
    )

    if haskey(io, :color) && io[:color] == false
        color = nothing
    else
        color = garish_io.color
    end

    return GarishIO(
        io, indent,
        compact, limit,
        displaysize,
        show_indent,
        include_defaults,
        color_prefs,
        state,
    )
end

Base.pipe_reader(io::GarishIO) = io.bland_io
Base.pipe_writer(io::GarishIO) = io.bland_io
Base.lock(io::GarishIO) = lock(io.bland_io)
Base.unlock(io::GarishIO) = unlock(io.bland_io)

function Base.IOContext(io::GarishIO, KVs::Pair...)
    IOContext(
        io.bland_io,
        :compact=>io.compact,
        :limit=>io.limit,
        :color=>!isnothing(io.color),
        :displaysize=>io.displaysize,
        :typeinfo=>io.state.typeinfo,

        :pprint_indent=>io.indent,
        :color_preference=>io.color,
        :show_indent=>io.show_indent,
        :include_defaults=>io.include_defaults,
        :pprint_type=>io.state.type,
        :pprint_level=>io.state.level,
        :pprint_offset=>io.state.offset,
        :noindent_in_first_line=>io.state.noindent_in_first_line,

        KVs...
    )
end

Base.displaysize(io::GarishIO) = io.displaysize
Base.in(key_value::Pair, io::GarishIO) = in(key_value, IOContext(io).dict, ===)
Base.haskey(io::GarishIO, key) = haskey(IOContext(io).dict, key)
Base.getindex(io::GarishIO, key) = getindex(IOContext(io).dict, key)
Base.get(io::GarishIO, key, default) = get(IOContext(io).dict, key, default)
Base.keys(io::GarishIO) = keys(IOContext(io).dict)

"""
    within_nextlevel(f, io::GarishIO)

Run `f()` within the next level of indentation where `f` is a function
that print into `io`.
"""
function within_nextlevel(f, io::GarishIO)
    io.state.level += 1
    ret = f()
    io.state.level -= 1

    # upperlevel_type = io.state.type
    #     upperlevel_noindent_in_first_line = io.state.noindent_in_first_line
    #     io.state.type = StructField
    #     io.state.noindent_in_first_line = true
    #     pprint(io, s.second)
    #     io.state.noindent_in_first_line = upperlevel_noindent_in_first_line
    #     io.state.type = upperlevel_type
    return ret
end

"""
    print_token(io::GarishIO, type::Symbol, xs...)

Print `xs` to a `GarishIO` as given token type. The token type
should match the field name of `ColorPreference`.
"""
function print_token(io::GarishIO, type::Symbol, xs...)
    print_token(print, io, type, xs...)
end

"""
    print_token(f, io::GarishIO, type::Symbol, xs...)

Print `xs` to a `GarishIO` as given token type using `f(io, xs...)`
"""
function print_token(f, io::GarishIO, type::Symbol, xs...)
    isnothing(io.color) && return f(io, xs...)
    Base.with_output_color(f, getfield(io.color, type), io, xs...)
end

function max_indent_reached(io::GarishIO, offset::Int)
    io.indent * io.state.level + io.state.offset + offset > io.displaysize[2]
end


text:
# SPDX-License-Identifier: MIT

using CUDA
using Documenter
using ExpFamilyDistributions
using LinearAlgebra
using SpecialFunctions: loggamma, digamma
using Test

CUDA.allowscalar(false)

const types = [Float32, Float64]
const arraytypes = [Array]

if CUDA.functional()
    @info "gpu available testing on CuArray"
    CUDA.versioninfo()
    push!(arraytypes, CuArray)
end

for T in types

    @testset "PD matrix operations ($T)" begin
        M = T[1 0.5; 0.5 1]
        @test all(inv(M) .≈ pdmat_inverse(M))
        @test all(logdet(M) .≈ pdmat_logdet(M))
    end

    for AT in arraytypes

        @testset "DefaultParameter ($AT | $T)" begin
            ξ = AT(T[1.0, 2.0, 3.0])
            p = DefaultParameter(ξ)

            @test all(naturalform(p) .≈ ξ)
            @test all(realform(p) .≈ ξ)
            @test all(Array(jacobian(p)) .≈ Matrix{T}(I, length(ξ), length(ξ)))
            p̄ = reallocate(p, Tuple)
            @test typeof(realform(p̄)) == Tuple{T, T, T}

            dict = todict(p)
            np = fromdict(typeof(p), dict)
            @test all(realform(p) .≈ realform(np))
            @test typeof(realform(p)) == typeof(realform(np))
        end

        @testset "Normal ($AT | $T)" begin
            n = Normal(AT(T[1, 2]), AT(T[2 0; 0 3]))
            @test typeof(n).mutable

            η = AT(T[1/2, 2/3, -1/4, -1/6, 0])
            @test all(naturalform(n.param) .≈ η)

            x = AT(T[2, 3])
            Tx = AT(T[2, 3, 4, 9, 6])
            @test all(stats(n, x) .≈ Tx)

            A = .5 * pdmat_logdet(n.Σ) + .5 * dot(n.μ, pdmat_inverse(n.Σ) * n.μ)
            @test typeof(lognorm(n)) == T
            @test lognorm(n) ≈ T(A)

            @test basemeasure(n, x) ≈ T(-log(2π))

            ETx = AT(T[1, 2, 3, 7, 2])
            @test all(gradlognorm(n) .≈ ETx)

            s1, s2, s3 = splitgrad(n, gradlognorm(n))
            @test size(s1) == (2,)
            @test size(s2) == (2,)
            @test size(s3) == (1,)

            n2 = Normal(AT(zeros(T, 2)), AT(Matrix{T}(I, 2, 2)))
            @test isapprox(kldiv(n, n), 0; atol = (T == Float32 ? 1e-6 : 1e-16))
            @test kldiv(n, n2) >= 0
            @test kldiv(n2, n) >= 0

            n2 = Normal(stdparam(n, naturalform(n.param))...)
            @test isapprox(kldiv(n, n2),  0; atol = (T == Float32 ? 1e-6 : 1e-16))

            @test size(sample(n2, 1)) == (2, 1)
            @test size(sample(n2, 10)) == (2, 10)
            @test eltype(sample(n2, 2)) == T
        end

        @testset "NormalDiag ($AT | $T)" begin
            n = NormalDiag(AT(T[1, 2]), AT(T[2, 2]))
            @test typeof(n).mutable

            η = AT(T[.5, 1, -.25, -.25])
            @test all(naturalform(n.param) .≈ η)

            x = AT(T[2, 3])
            Tx = AT(T[2, 3, 4, 9])
            @test all(stats(n, x) .≈ Tx)

            A = .5 * pdmat_logdet(n.Σ) + .5 * dot(n.μ, pdmat_inverse(n.Σ) * n.μ)
            @test typeof(lognorm(n)) == T
            @test lognorm(n) ≈ T(A)

            @test basemeasure(n, x) ≈ T(-log(2π))

            ETx = AT(T[1, 2, 3, 6])
            @test all(gradlognorm(n) .≈ ETx)

            s1, s2 = splitgrad(n, gradlognorm(n))
            @test size(s1) == (2,)
            @test size(s2) == (2,)

            n2 = NormalDiag(AT(zeros(T, 2)), AT(ones(T, 2)))
            @test kldiv(n, n) ≈  0
            @test kldiv(n, n2) >= 0
            @test kldiv(n2, n) >= 0

            n2 = NormalDiag(stdparam(n, naturalform(n.param))...)
            @test kldiv(n, n2) ≈  0

            @test size(sample(n, 10)) == (2, 10)
            @test eltype(sample(n, 2)) == T
       end

        @testset "Gamma ($AT | $T)" begin
            g = Gamma(AT(T[1, 2]), AT(T[2, 3]))
            @test typeof(g).mutable

            η = AT(T[-2, -3, 1, 2])
            @test all(naturalform(g.param) .≈ η)

            x = AT(T[1, 2])
            Tx = vcat(x, log.(x))
            @test all(stats(g, x) .≈ Tx)

            A = loggamma(1) + loggamma(2) - 1*log(2) - 2*log(3)
            @test lognorm(g) ≈ T(A)

            @test basemeasure(g, x) ≈ sum(-log.(x))

            ETx = AT(T[1/2, 2/3, digamma(1) - log(2), digamma(2) - log(3)])
            @test all(gradlognorm(g) .≈ ETx)

            s1, s2 = splitgrad(g, gradlognorm(g))
            @test length(s1) == 2
            @test length(s2) == 2

            g2 = Gamma(AT(T[1, 1]), AT(T[1, 1]))
            @test kldiv(g, g) ≈  0
            @test kldiv(g, g2) >= 0
            @test kldiv(g2, g) >= 0

            g2 = Gamma(stdparam(g, naturalform(g.param))...)
            @test kldiv(g, g2) ≈  0

            @test size(sample(g, 10)) == (2,10)
        end

        @testset "Dirichlet ($AT | $T)" begin
            d = Dirichlet(AT(T[1, 2, 3]))
            @test typeof(d).mutable

            η = AT(T[1, 2, 3])
            @test all(naturalform(d.param) .≈ η)

            x = AT(T[1, 2, 3])
            Tx = AT(T[log(1), log(2), log(3)])
            @test all(stats(d, x) .≈ Tx)

            A = sum(loggamma_dot(d.α)) - loggamma(sum(d.α))
            @test lognorm(d) ≈ T(A)

            @test basemeasure(d, x) ≈ -log.(x)

            ETx = digamma.(d.α) .- digamma(sum(d.α))
            @test all(gradlognorm(d) .≈ ETx)

            s = splitgrad(d, gradlognorm(d))
            @test size(s) == (3,)

            d2 = Dirichlet(AT(ones(T, 3)))
            @test kldiv(d, d) ≈  0
            @test kldiv(d, d2) >= 0
            @test kldiv(d2, d) >= 0

            d2 = Dirichlet(stdparam(d, naturalform(d.param))...)
            @test kldiv(d, d2) ≈  0

            @test size(sample(d, 10)) == (3, 10)
        end

        @testset "Wishart($AT | $T)" begin
            X = AT(AT(T[1 0.5; 0.5 2]))

            D = T(2)
            W = AT(Matrix{T}(I, 2, 2))
            v = T(2)
            w = Wishart(W, v)
            @test typeof(w).mutable

            invW = pdmat_inverse(W)
            η = vcat(-.5 * diag(invW), vec_tril(invW), v/2)
            @test eltype(naturalform(w.param)) == T
            @test all(naturalform(w.param) .≈ η)

            TX = vcat(diag(X), vec_tril(X), pdmat_logdet(X))
            @test all(stats(w, X) .≈ TX)

            A = .5*( v*pdmat_logdet(W) + v*D*log(2) ) + sum([loggamma((v+1-i)/2) for i in 1:D])
            @test lognorm(w) ≈ T(A)

            B = -.5*( (D-1)*pdmat_logdet(X) + .5*D*(D-1)*log(π) )
            @test basemeasure(w, X) ≈ B

            vW = v*W
            ETX = vcat(diag(vW), vec_tril(vW), sum([digamma(T(0.5)*(v+1-i)) for i in 1:D]) + D*T(log(2)) + pdmat_logdet(W))
            @test eltype(gradlognorm(w)) == T
            @test all(gradlognorm(w) .≈ ETX)

            s = splitgrad(w, gradlognorm(w))
            @test length(s) == 3
            @test all(s[1] .≈ diag(vW))
            @test all(s[2] .≈ vec_tril(vW))
            @test s[3] .≈ sum([digamma((v+1-i)/2) for i in 1:D]) + D*log(2) + pdmat_logdet(W)

            w2 = Wishart(AT(Matrix{T}(I, 2, 2)), 2)
            @test kldiv(w, w) ≈  0
            @test kldiv(w, w2) >= 0
            @test kldiv(w2, w) >= 0

            w2 = Wishart(stdparam(w, naturalform(w.param))...)
            @test kldiv(w, w2) ≈  0

            @test size(sample(w, 10)) == (2, 2, 10)
        end
    end
end

# Finally, test the documentation
DocMeta.setdocmeta!(ExpFamilyDistributions, :DocTestSetup,
                    :(using ExpFamilyDistributions), recursive = true)
doctest(ExpFamilyDistributions)

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions


text:
# Second version of the new_model.jl
# Parameter iterators
postagv(m)=m[1];
bufferm(m)=m[2];
stackm(m)=m[3];
actm(m)=m[4];
actembed(m)=m[5];
parserv(m)=m[6];

# To initialize parser model from scratch
function initmodel1(o, s; ftype=Float32)
    model = Any[]

    # postag-level initialization, KnetArray initialization
    dpostag = o[:posembed] # 1
    for (k, n, d) in ((:postag, 17, dpostag),)
        push!(model, [ initr(d, GPUFEATS=true) for i=1:n ])
    end

    # buffer-lstm initialization # 2
    r_b, wr_b = initbuf(o[:bufembed], o[:bufhidden]) # you may change lstm - gru with kwargs
    push!(model, [r_b, wr_b])

    # stack-lstm initialization # 3
    r_s, wr_s = initbuf(o[:stembed], o[:sthidden])
    push!(model, [r_s, wr_s])

    # action-lstm initialization # 4
    r_a, wr_a = initbuf(o[:actembed], o[:acthidden])
    push!(model, [r_a, wr_a])

    # action embeddings
    p = o[:arctype](s)
    push!(model, initx(o[:actembed], p.nmove)) # 5

    # decision module initialization # 6
    mlpdims = (r_s.hiddenSize + r_b.hiddenSize + r_a.hiddenSize, o[:hidden]..., p.nmove)
    info("mlpdims: $mlpdims")
    decider = Any[]
    for i=2:length(mlpdims)
        push!(decider, initx(mlpdims[i], mlpdims[i-1]))
        push!(decider, initx(mlpdims[i], 1))
    end
    push!(model, decider)
    optims = optimizers(model, o[:optimization])
    return model, optims
end

function mlp(w,x; pdrop=(0,0))
    x = dropout(x,pdrop[1])
    for i=1:2:length(w)-2
        x = relu.(w[i]*x .+ w[i+1])
        x = dropout(x,pdrop[2])
    end
    return w[end-1]*x .+ w[end]
end

function scan_buffer(model, parser, hidden)
    bufmodel, pvecs = bufferm(model), postagv(model)
    rbuf, wbuf = bufmodel[1], bufmodel[2]
    sentence = parser.sentence
    seqe = length(sentence)
    seqs = parser.wptr
    range = seqe:-1:seqs
    instate = zeros(Float32, hidden, 1) # TODO: not to hardcode hidden
    ybuf = hout = cout = (gpu() >=0 ? KnetArray(instate) : Array(instate))
    for i in range
        #word = sentence.word[i]; print(" $word"); # dbg line
        x = (gpu() >=0 ? KnetArray(sentence.cavec[i]) : Array(sentence.cavec[i]))
        input = cat1d(x, pvecs[sentence.postag[i]])
        ybuf, hout, cout, _ = rnnforw(rbuf, wbuf, input, hout, cout, hy=true, cy=true)
    end
    return ybuf
end

function scan_bufbatch(model, parsers, hidden)
    # Todo: we may cache some parsers in future
    yalls = Any[]
    for p in parsers
        ybuf = scan_buffer(model, p, hidden)
        push!(yalls, ybuf)
    end
    yalls = cat1d(yalls...)
    ncols = length(parsers)
    nrows = div(length(yalls), ncols)
    return reshape(yalls, nrows, ncols)
end

function scan_stack(model, parser, hstack)
    stmodel, pvecs = stackm(model), postagv(model)
    rst, wst = stmodel[1], stmodel[2]
    sentence = parser.sentence
    seqe = parser.sptr
    instate = zeros(Float32, hstack, 1)
    yst = hout = cout = (gpu() >= 0 ? KnetArray(instate) : Array(instate))
    range = 1:seqe # words in stack p.stack[1:p.sptr]
    for i in range
        indx  = parser.stack[i]
        #word = sentence.word[indx]; print(" $word") # dbg line
        x = (gpu() >= 0 ? KnetArray(sentence.cavec[indx]) : Array(sentence.cavec[indx]))
        input = cat1d(x, pvecs[sentence.postag[indx]])
        yst, hout, cout, _ = rnnforw(rst, wst, input, hout, cout, hy=true, cy=true)
    end
    return yst
end

function scan_stackbatch(model, parsers, hstack)
    yalls = Any[]
    for p in parsers
        yst = scan_stack(model, p, hstack)
        push!(yalls, yst)
    end
    yalls = cat1d(yalls...)
    ncols = length(parsers)
    nrows = div(length(yalls), ncols)
    return reshape(yalls, nrows, ncols)
end

# compute batched actions
function scan_action(model, actions, state)
    acemb, actmodel = actembed(model), actm(model)
    yout0, hout0, cout0 = state[1], state[2], state[3]
    input = acemb[:, actions]
    yout, hout, cout = rnnforw(actmodel[1], actmodel[2], input, hout0, cout0, hy=true, cy=true)
    return (yout, hout, cout)
end

function oracleloss(allmodel, sentences, arctype, hiddens; losses=nothing, pdrop=(0.0, 0.0))
    hidstack, hidact, hidbuf = hiddens
    parsers = map(arctype, sentences)
    mcosts  = Array{Cost}(parsers[1].nmove)
    parserdone = falses(length(parsers))
    mlpmodel = parserv(allmodel)
    B = length(parsers) # bathcsize
    t0 = zeros(Float32, hidact, B, 1)
    h0 = c0 = (gpu()>=0 ? KnetArray(t0) : Array(t0))
    y0 = reshape(h0, hidact, B)
    totalloss = 0.0

    # action lstm-init
    nstate = (y0, h0, c0);
    inacts = map(t->1, 1:B)

    while !all(parserdone)
        # batched score
        yact, hact, cact = scan_action(allmodel, inacts, nstate)
        nstate = (yact, hact, cact);
        ybufs    = scan_bufbatch(allmodel, parsers, hidbuf)
        ysts     = scan_stackbatch(allmodel, parsers, hidstack)
        encoded  = vcat(yact, ybufs, ysts)
        scores   = mlp(mlpmodel, encoded, pdrop=pdrop)
        logprobs = logp(scores, 1)

        inacts = Int[] # to get the next step's actions
        # iterative loss-val
        for (i, p) in enumerate(parsers)
            if parserdone[i]
                push!(inacts, 1); # not changing batchsize
                continue
            end
            movecosts(p, p.sentence.head, p.sentence.deprel, mcosts)
            goldmove = indmin(mcosts)

            if mcosts[goldmove] == typemax(Cost)
                parserdone[i] = true
                p.sentence.parse = p
                push!(inacts, 1); # not changing batchsize
            else
                totalloss -= logprobs[goldmove, i]
                move!(p, goldmove);
                push!(inacts, goldmove);
                if losses != nothing
                    loss1 = -getval(logprobs)[goldmove,i]
                    losses[1] += loss1
                    losses[2] += 1
                    if losses[2] < 1000
                        losses[3] = losses[1]/losses[2]
                    else
                        losses[3] = 0.999 * losses[3] + 0.001 * loss1
                    end
                end
            end
        end
    end
    return totalloss / length(sentences)
end

oraclegrad = grad(oracleloss)

# Implement test related

# endofparse to use in test
endofparse(p)=(p.sptr == 1 && p.wptr > p.nword)

function oracletest(allmodel, corpus, arctype, hiddens, batchsize; pdrop=(0.0, 0.0))
    sentbatches = minibatch(corpus, batchsize)
    hidstack, hidact, hidbuf = hiddens
    mlpmodel = parserv(allmodel)

    for sentences in sentbatches
        parsers = map(arctype, sentences)
        B = length(parsers) # bathcsize
        mcosts  = Array{Cost}(parsers[1].nmove)
        parserdone = map(endofparse, parsers)
        #parserdone = falses(length(parsers))
        t0 = zeros(Float32, hidact, B, 1)
        h0 = c0 = (gpu()>=0 ? KnetArray(t0) : Array(t0))
        y0 = reshape(h0, hidact, B)
        nstate = (y0, h0, c0);
        inacts = map(t->1, 1:B)
        yact, hact, cact = scan_action(allmodel, inacts, nstate)
        iters = 1#dbg
        while !all(parserdone)
            yact, hact, cact = scan_action(allmodel, inacts, nstate)
            nstate = (yact, hact, cact);
            ybufs    = scan_bufbatch(allmodel, parsers, hidbuf)
            ysts     = scan_stackbatch(allmodel, parsers, hidstack)
            encoded  = vcat(yact, ybufs, ysts)
            scores   = mlp(mlpmodel, encoded, pdrop=pdrop)
            logprobs = Array(logp(scores, 1))

            inacts = Int[] # to get the next step's actions
            for (i, p) in enumerate(parsers)
                if parserdone[i]
                    if p.sentence.parse == nothing # to avoid from the first time
                        p.sentence.parse = p
                    end
                    push!(inacts, 1) # not to break batchsize
                    continue
                end
                isorted = sortperm(logprobs[:,i], rev=true) # ith col ith inst.
                for m in isorted # best&val action
                    if moveok(p, m)
                        move!(p, m)
                        push!(inacts, m) # not to break batchsize
                        break
                    end
                end
                if endofparse(p)
                    parserdone[i] = true
                    p.sentence.parse = p
                end
            end
        end
    end
    return las(corpus)
end

# labeled-attachment score calculator
function las(corpus)
    nword = ncorr = 0
    for s in corpus
        p = s.parse
        nword += length(s)
        ncorr += sum((s.head .== p.head) .& (s.deprel .== p.deprel))
    end
    ncorr / nword
end

function empty_parses!(corpus)
    for s in corpus
        s.parse = nothing
    end
end


text:
# common utilities

## convenient type alias
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}
#

const RealArray{T<:Real,N} = AbstractArray{T,N}
const RealVector{T<:Real} = AbstractArray{T,1}
const RealMatrix{T<:Real} = AbstractArray{T,2}

const IntegerArray{T<:Integer,N} = AbstractArray{T,N}
const IntegerVector{T<:Integer} = AbstractArray{T,1}
const IntegerMatrix{T<:Integer} = AbstractArray{T,2}

const RealFP = Union{Float32, Float64}

## conversion from real to fp types

fptype(::Type{T}) where {T<:Union{Float32,Bool,Int8,UInt8,Int16,UInt16}} = Float32
fptype(::Type{T}) where {T<:Union{Float64,Int32,UInt32,Int64,UInt64,Int128,UInt128}} = Float64
fptype(::Type{Complex{Float32}}) = Complex{Float32}
fptype(::Type{Complex{Float64}}) = Complex{Float64}

# A convenient typealias for deprecating default corrected Bool
const DepBool = Union{Bool, Nothing}

function depcheck(fname::Symbol, b::DepBool)
    if b == nothing
        msg = "$fname will default to corrected=true in the future. Use corrected=false for previous behaviour."
        Base.depwarn(msg, fname)
        false
    else
        b
    end
end


text:
update{N, T, T2}(b::HyperRectangle{N, T}, v::Vec{N, T2}) =
    update(b, Vec{N, T}(v))
function update{N, T}(b::HyperRectangle{N, T}, v::Vec{N, T})
    m = min(minimum(b), v)
    mm = max(maximum(b), v)-m
    HyperRectangle{N, T}(m, mm)
end

# Min maximum distance functions between hrectangle and point for a given dimension
@inline function min_dist_dim{N, T}(rect::HyperRectangle{N, T}, p::Vec{N, T}, dim::Int)
    max(zero(T), max(minimum(rect)[dim] - p[dim], p[dim] - maximum(rect)[dim]))
end

@inline function max_dist_dim{N, T}(rect::HyperRectangle{N, T}, p::Vec{N, T}, dim::Int)
    max(maximum(rect)[dim] - p[dim], p[dim] - minimum(rect)[dim])
end

@inline function min_dist_dim{N, T}(rect1::HyperRectangle{N, T},
                                    rect2::HyperRectangle{N, T},
                                    dim::Int)
    max(zero(T), max(minimum(rect1)[dim] - maximum(rect2)[dim],
                     minimum(rect2)[dim] - maximum(rect1)[dim]))
end

@inline function max_dist_dim{N, T}(rect1::HyperRectangle{N, T},
                                    rect2::HyperRectangle{N, T},
                                    dim::Int)
    max(maximum(rect1)[dim] - minimum(rect2)[dim],
        maximum(rect2)[dim] - minimum(rect1)[dim])
end

# Total minimum maximum distance functions
@inline function min_euclideansq{N, T}(rect::HyperRectangle{N, T},
                                       p::Union{Vec{N, T},
                                       HyperRectangle{N, T}})
    minimum_dist = T(0.0)
    for dim in 1:length(p)
        d = min_dist_dim(rect, p, dim)
        minimum_dist += d*d
    end
    return minimum_dist
end

# could add an @symmetric macro, which defines min_euclidean(pt, s) from
# min_euclidean(s, pt) automatically etc.
@inline min_euclidean(pt1::Vec, pt2::Vec) = norm(pt1-pt2)
min_euclidean(pt::Vec, s::Simplex) = √(sqdist(pt, s))
min_euclidean(s::Simplex, pt::Vec) = min_euclidean(pt, s)
min_euclidean(r1::HyperRectangle, r2::Vec) = sqrt(min_euclideansq(r1, r2))
min_euclidean(r1::Vec, r2::HyperRectangle) = sqrt(min_euclideansq(r1, r2))
min_euclidean(r1::HyperRectangle, r2::HyperRectangle) = sqrt(min_euclideansq(r1, r2))
min_euclidean(c1::AbstractConvexHull, pt::Vec) = gjk(c1,pt)
min_euclidean(pt::Vec, c2::AbstractConvexHull) = gjk(pt,c2)
min_euclidean(c1::AbstractConvexHull, c2::AbstractConvexHull) = gjk(c1,c2)

@inline function max_euclideansq{N, T}(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}})
    maximum_dist = T(0.0)
    for dim in 1:length(p)
        d = max_dist_dim(rect, p, dim)
        maximum_dist += d*d
    end
    return maximum_dist
end
function max_euclidean{N, T}(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}})
    sqrt(max_euclideansq(rect, p))
end

# Functions that return both minimum and maximum for convenience
@inline function minmax_dist_dim{N, T}(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}, dim::Int)
    minimum_d = min_dist_dim(rect, p, dim)
    maximum_d = max_dist_dim(rect, p, dim)
    return minimum_d, maximum_d
end

@inline function minmax_euclideansq{N, T}(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}})
    minimum_dist = min_euclideansq(rect, p)
    maximum_dist = max_euclideansq(rect, p)
    return minimum_dist, maximum_dist
end

function minmax_euclidean{N, T}(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}})
    minimumsq, maximumsq = minmax_euclideansq(rect, p)
    return sqrt(minimumsq), sqrt(maximumsq)
end


text:
module TestFEM
    using FactCheck
    using FEM

    include("test_geomesh.jl")
    include("test_materials.jl")
    include("test_elements.jl")
    include("test_interpolators.jl")
    include("test_feproblem.jl")
    #include("test_vtkexport.jl")
end


text:
function compute_derivatives(obj::FieldAdditionalEquations, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "computeDerivatives", Vector{RealFieldElement}, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_name(obj::FieldAdditionalEquations)
    return jcall(obj, "getName", JString, ())
end

function init(obj::FieldAdditionalEquations, arg0::FieldSpacecraftState, arg1::FieldAbsoluteDate)
    return jcall(obj, "init", void, (FieldSpacecraftState, FieldAbsoluteDate), arg0, arg1)
end


text:
struct Control
    iter::Int ## total number of iterations (incl. warmup)
    warmup::Int ## number of warmup iterations
    chains::Int ## number of chains
    save_warmup::Bool ## store and return warmup iterations in Fit
    common_delta::Bool ## common log-HRs across classes
    common_beta::Bool ## common fixed effects across classes
    adapt_delta::Float64 ## target mean acceptance rate for NUTS transitions
end

Control(iter, chains) = Control(iter=iter, chains=chains)

function Control(;iter=2000, warmup=Int(floor(iter/2)), chains=4, save_warmup=false,
                 common_delta=true, common_beta=true, adapt_delta=0.8)
    @argcheck iter > warmup
    Control(iter, Int(floor(iter/2)), chains, save_warmup, common_delta, common_beta, adapt_delta)
end

struct Chain{Tpar<:AbstractParam, Tprior<:AbstractPrior}
    theta::Vector{Tpar}
    accept::Tpar
    loglik_p::Matrix{Float64} ## pointwise log-likelihood
    prior::Tprior
    init::Tpar
end

struct Fit{Tpar<:AbstractParam, Tprior<:AbstractPrior} <: AbstractFit
    control::Control
    chains::Vector{Chain{Tpar, Tprior}}
end

function show(io::IO, x::Data)
    println(io, "Subjects: $(length(x.T))")
    println(io, "Repeated measurements: $(sum(x.n))")
    println(io, "Events: $(sum(x.E))")
end

function show(io::IO, x::Chain{T,S}) where {T,S}
    cm = countmap([par.G for par in x.theta])
    cm0 = countmap([par.G-par.G0 for par in x.theta])

    println(io, "Number of groups: $(collect(zip(keys(cm), values(cm))))")
    println(io, "Number of non-empty groups: $(collect(zip(keys(cm0), values(cm0))))")
end

function show(io::IO, x::Fit{T}) where T <: AbstractParam
    println(io, "Number of chains: $(x.control.chains)")
    println(io, "Iterations: $(x.control.iter)")
    for i in eachindex(x.chains)
        println(io, "Chain $(i):")
        show(io, x.chains[i])
        println(io, "\n")
    end
end

merge(a...) = foldr(merge, a)
merge(a::Chain{T,S}) where {T,S} = a

function merge(a::Chain{T,S}, b::Chain{T,S}) where {T <: AbstractParam, S <: AbstractPrior}
    Chain{T,S}(vcat(a.theta, b.theta), a.accept, vcat(a.loglik_p, b.loglik_p), a.prior, a.init)
end

size(x::Fit) = (length(x.chains), length(first(x.chains)))
slice(x::Chain, t) = x.theta[t]

length(x::Fit) = length(x.chains)
length(x::Chain) = length(x.theta)

function get_samples(x::Fit)
    K, M = size(x)
    res = Matrix{typeof(first(x.chains).init)}(undef, M, K)

    for k in 1:K
        res[:,k] = get_samples(x.chains[k])
    end
    res
end

get_samples(x::Chain) = x.theta

get_loglik_p(x::Chain) = x.loglik_p
get_loglik_p(x::Fit) = get_loglik_p(merge(x.chains...))

get_acceptance_rates(x::Fit) = get_acceptance_rates(merge(x.chains...))
get_acceptance_rates(x::Chain) = x.accept[x.theta[1].G]

## Calculate class-membership probabilities for each subject
function class_prob(xi, Xp)
    cp = exp.(Xp * xi')  ## N x G
    cp ./ sum(cp, dims=2)[:,1] ## normalize each row
end

function build_cov_matrix(S::Vector{T}, L::Matrix{T}) where T
    SD = Diagonal(S)
    Array(Symmetric(SD*L*L'*SD))
end

build_cov_matrix(theta::Tp, g::Int) where Tp <: AbstractParam = buidl_cov_matrix(theta.S[g,:], theta.L[g,:,:])

function post_cind_prob(theta::ParamSpline, data)
    @unpack Y, times, n, Xl, Z, Xp, Xe, T, E, haz_basis, cumhaz_basis = data
    @unpack xi, beta, mu, sigma, w, delta, S, L = theta

    G, q = size(mu)
    N = length(data.T)
    cp = class_prob(xi, Xp)
    res = Matrix{Float64}(undef, N, G)

    for g in 1:G
        Sigma = build_cov_matrix(S[g,:], L[g,:,:])

        for i in 1:N
            nt = 1:n[i]
            lll_vec = loglik_p_long(beta[g,:], mu[g,:], Sigma, sigma, Y[i,1:n[i]], times[i,1:n[i]], n[i], Xl[i,1:n[i],:], Z[i,1:n[i],:])
            lle_vec = loglik_p_surv(T[i], E[i], Xe[i:i,:], haz_basis[i:i,:], cumhaz_basis[i:i,:], w[g,:], delta[g,:])
            res[i,g] = lll_vec + lle_vec + log(cp[g])
        end
    end
    res
end

function post_cind_prob(theta::Param, data)
    @unpack Y, times, n, Xl, Z, Xp, Xe, T, E = data
    @unpack xi, beta, mu, sigma, bl_shape, bl_scale, delta, S, L = theta

    G, q = size(mu)
    N = length(data.T)
    cp = class_prob(xi, Xp)
    res = Matrix{Float64}(undef, N, G)

    for g in 1:G
        Sigma = build_cov_matrix(S[g,:], L[g,:,:])

        for i in 1:N
            nt = 1:n[i]
            lll_vec = loglik_p_long(beta[g,:], mu[g,:], Sigma, sigma, Y[i,1:n[i]], times[i,1:n[i]], n[i], Xl[i,1:n[i],:], Z[i,1:n[i],:])
            lle_vec = loglik_p_surv(T[i], E[i], Xe[i:i,:], bl_shape[g], bl_scale[g], delta[g,:])
            res[i,g] = lll_vec + lle_vec + log(cp[g])
        end
    end
    res
end

## sample from discrete distribution with probabilities given by exp.(alpha)/sum(exp.(alpha)) using Gumbel-max trick
sample_discrete(alpha) = last(findmax(alpha + rand(Gumbel(), length(alpha))))

function sample_cind(theta::Tp, data::Data) where Tp <: AbstractParam
    alpha = post_cind_prob(theta, data)
    [sample_discrete(r) for r in eachrow(alpha)]
end

function sample_sigma(b, beta, mu, sigma, Xl, Z, Y, n, times, cind, prior)
    N = size(Xl)[1]
    q = size(Xl)[3]
    v = 0.0

    for i in 1:N
        g = cind[i]
        u = b[i,:] ##mu[g,:] + b[i,:]

        for j in 1:n[i]
            v += (Y[i,j] - dot(Z[i,j,:], u) - dot(Xl[i,j,:], beta[g,:]))^2
        end
    end
    v /= 2

    log_target(x::Float64) = logpdf(InverseGamma(sum(n)/2 - 1, v), exp(2*x)) + logpdf(prior, exp(x)) + x

    y, ac = metropolis(log_target, log(sigma), rand(Normal(log(sigma), 0.1)))

    exp(y), ac
end

function sample_mu_conj(Sigma, beta, mu, sigma, Xl, Z, Y, n, prior_mu)
    N = size(Y,1)
    q = length(mu)

    V0 = inv(cov(prior_mu))
    ll_mean_mu = zeros(Float64, q)
    ll_inv_cov_mu = zeros(Float64, q, q)

    for i in 1:N
        tt = 1:n[i]
        Zi = Z[i,tt,:]
        Xli = Xl[i,tt,:]
        C = inv(Array(Symmetric(Zi*Sigma*Zi' + sigma^2*diagm(0 => ones(Float64, n[i])))))
        ll_mean_mu += Zi'*C*(Y[i,tt] - Xli*beta)
        ll_inv_cov_mu += Zi'*C*Zi
    end

    post_cov = inv(Symmetric(V0 + ll_inv_cov_mu))
    post_mean = post_cov * (V0*mean(prior_mu) + ll_mean_mu)
    rand(MvNormal(post_mean, post_cov)), 1
end

function sample_beta_conj(theta, data, prior_beta)
    @unpack beta, mu, sigma, S, L, cind = theta
    @unpack Xl, Z, Y, n = data

    N = size(Y,1)
    l = size(beta,2)
    G, q = size(mu)

    V0 = inv(cov(prior_beta))

    ll_mean_beta = zeros(Float64, l)
    ll_inv_cov_beta = zeros(Float64, l, l)

    Sigma = zeros(eltype(beta), G, q, q)
    for g in 1:G
        Sigma[g,:,:] = build_cov_matrix(S[g,:], L[g,:,:])
    end

    for i in 1:N
        tt = 1:n[i]
        Zi = Z[i,tt,:]
        Xli = Xl[i,tt,:]
        g = cind[i]
        C = inv(Array(Symmetric(Zi*Sigma[g,:,:]*Zi' + sigma^2*diagm(0 => ones(Float64, n[i])))))
        ll_mean_beta += Xli'*C*(Y[i,tt] - Zi*mu[g,:])
        ll_inv_cov_beta += Xli'*C*Xli
    end

    post_cov = inv(Symmetric(V0 + ll_inv_cov_beta))
    post_mean = post_cov * (V0*mean(prior_beta) + ll_mean_beta)
    rand(MvNormal(post_mean, post_cov)), 1
end

function sample_mu_beta_conj(Sigma, beta, mu, sigma, Xl, Z, Y, n, prior_mu, prior_beta)
    N = size(Y,1)
    q = length(mu)
    l = length(beta)

    mb0 = vcat(mean(prior_mu), mean(prior_beta))
    V0 = zeros(Float64, q+l, q+l)
    ## assuming independent priors for mu and beta
    V0[1:q, 1:q] = inv(cov(prior_mu))
    V0[(q+1):end, (q+1):end] = inv(cov(prior_beta))

    ll_mean_mu = zeros(Float64, q)
    ll_mean_beta = zeros(Float64, l)
    ll_inv_cov_mu = zeros(Float64, q, q)
    ll_inv_cov_beta = zeros(Float64, l, l)

    for i in 1:N
        tt = 1:n[i]
        Zi = Z[i,tt,:]
        Xli = Xl[i,tt,:]
        C = inv(Array(Symmetric(Zi*Sigma*Zi' + sigma^2*diagm(0 => ones(Float64, n[i])))))
        ll_mean_mu += Zi'*C*(Y[i,tt] - Xli*beta)
        ll_mean_beta += Xli'*C*(Y[i,tt] - Zi*mu)
        ll_inv_cov_mu += Zi'*C*Zi
        ll_inv_cov_beta += Xli'*C*Xli
    end

    ll_inv_cov = zeros(Float64, q+l, q+l)
    ll_inv_cov[1:q, 1:q] = ll_inv_cov_mu
    ll_inv_cov[(q+1):end, (q+1):end] = ll_inv_cov_beta

    post_cov = inv(Symmetric(V0 + ll_inv_cov))
    post_mean = post_cov * (V0*mb0 + vcat(ll_mean_mu, ll_mean_beta))

    mu_beta = rand(MvNormal(post_mean, post_cov))
    mu_beta[1:q], mu_beta[(q+1):end], 1
end

function inv_cov_SL(S, L)
    ##try
        invSL = inv(LowerTriangular(Diagonal(S)*L))
        return invSL'invSL
    ##catch ex
    ##    @info S, L
    ##    error("singular")
    ##end
    ##return diagm(0 => ones(Float64, length(S)))
end

function sample_b_conj(theta, data)
    @unpack G, b, beta, mu, sigma, cind, S, L = theta
    @unpack T, E, Z, Xe, Xl, Z, Y, n, times = data

    N = length(T)
    q = size(b,2)
    z = Matrix{Float64}(undef, N, q)
    invSigma = [inv_cov_SL(S[g,:], L[g,:,:]) for g in 1:G]
    z0 = rand(MvNormal(diagm(0 => ones(eltype(b), q))), N)

    for i in 1:N
        g = cind[i]
        tt = 1:n[i]
        Zi = Z[i,tt,:]
        ll_mean = Zi'*(Y[i,tt] - Xl[i,tt,:]*beta[g,:])
        ll_inv_cov = Zi'*Zi / sigma^2

        V = Symmetric(invSigma[g] + ll_inv_cov)
        V_chol = cholesky!(V; check=false)

        if !issuccess(V_chol)
            z[i,:] = b[i,:]
        else
            VL = inv(V_chol.L)
            post_mean = VL'*VL*(invSigma[g]*mu[g,:] + ll_mean / sigma^2)
            z[i,:] = post_mean + VL'*z0[:,i]
        end
    end
    z, 1
end

function update_lmm_params_common!(theta::T, accept::T, data::Data, prior) where T <: AbstractParam
    G, q = size(theta.mu)

    theta.beta[1,:], ac = sample_beta_conj(theta, data, prior.beta_mv)

    for g in 1:G
        theta.beta[g,:] = theta.beta[1,:]
    end
    accept.beta .+= ac

    for g in 1:G
        sel = theta.cind .== g

        if sum(sel) > 0
            data_g = slice_data(data, sel)

            Sigma = build_cov_matrix(theta.S[g,:], theta.L[g,:,:])

            theta.mu[g,:], ac = sample_mu_conj(Sigma, theta.beta[g,:], theta.mu[g,:], theta.sigma,
                                               data_g.Xl, data_g.Z, data_g.Y, data_g.n, prior.mu_mv)
            accept.mu[g,:] .+= ac
        else
            theta.mu[g,:] = rand(prior.mu_mv)
            accept.mu[g,:] .+= 1
        end
    end

    theta.b, ac = sample_b_conj(theta, data)
    accept.b .+= ac

    theta.sigma, ac = sample_sigma(theta.b, theta.beta, theta.mu, theta.sigma, data.Xl, data.Z, data.Y, data.n, data.times, theta.cind, prior.sigma)
    accept.sigma += ac
end

function update_lmm_params_mix!(theta::T, accept::T, data::Data, prior) where T <: AbstractParam
    G, q = size(theta.mu)

    for g in 1:G
        sel = theta.cind .== g

        if sum(sel) > 0
            data_g = slice_data(data, sel)

            Sigma = build_cov_matrix(theta.S[g,:], theta.L[g,:,:])

            theta.mu[g,:], theta.beta[g,:], ac = sample_mu_beta_conj(Sigma, theta.beta[g,:], theta.mu[g,:], theta.sigma, data_g.Xl, data_g.Z, data_g.Y, data_g.n, prior.mu_mv, prior.beta_mv)
            accept.beta[g,:] .+= ac
            accept.mu[g,:] .+= ac
        else
            theta.mu[g,:] = rand(prior.mu_mv)
            accept.mu[g,:] .+= 1

            theta.beta[g,:] = rand(prior.beta_mv)
            accept.beta[g,:] .+= 1
        end
    end

    theta.b, ac = sample_b_conj(theta, data)
    accept.b .+= ac

    theta.sigma, ac = sample_sigma(theta.b, theta.beta, theta.mu, theta.sigma, data.Xl, data.Z, data.Y, data.n, data.times, theta.cind, prior.sigma)
    accept.sigma += ac
end

function transform_zero_sum(xi)
    G = size(xi,1)
    A = diagm(0 => ones(eltype(xi), G)) .- 1/G
    A*xi
end

function transform_ref_group(xi)
    G, p = size(xi)
    A = diagm(0 => ones(eltype(xi), G-1)) .- 1/G
    inv(A)*xi[1:(G-1),:]
end

function sample_xi_NUTS(sampler, xi, cind, Xp)
    G, p = size(xi)
    ##init = xi[1:(G-1),:]
    init = transform_ref_group(xi)
    res, ac = nuts_mlogit_step!(sampler, cind_to_matrix(cind, G), init)
    ##vcat(res.xi, zeros(eltype(res.xi), 1, p)), ac
    new_xi = transform_zero_sum(vcat(res.xi, zeros(eltype(res.xi), 1, p)))
    new_xi, ac
end

function update_cp_params!(theta::T, accept::T, data::Data, prior, sampler) where T <: AbstractParam
    G, q = size(theta.mu)

    if G == 1
        theta.cind .= 1
        return loglik(theta, data)
    else
        theta.cind = sample_cind(theta, data)
        accept.cind .+= 1

        theta.xi, ac = sample_xi_NUTS(sampler, theta.xi, theta.cind, data.Xp)
        ##theta.xi, ac = sample_xi(theta.xi, theta.cind, data.Xp)
        accept.xi .+= ac
    end
end

function cind_to_matrix(cind::Vector{Int}, G::Int)
    Y = zeros(Int, length(cind), G)
    for i in 1:length(cind)
        Y[i, cind[i]] = 1
    end
    Y
end

function init_NUTS_samplers(warmup, G, theta, data, prior, cp=Dict{Int, NUTSGibbsSampler}(),
                            surv=Dict{Int, NUTSGibbsSampler}(),
                            cov=Dict{Int, NUTSGibbsSampler}(); common_delta=false)
    if !haskey(surv, G)
        if G > 1
            cp[G] = nuts_mlogit_init(cind_to_matrix(theta.cind, G), data.Xp, prior.xi, warmup)
            @debug "Initialized NUTS sampler for class membership regression (G = $(G))"
        end

        if common_delta
            surv[G] = surv_model_init2(theta, data, prior, warmup)
        else
            surv[G] = surv_model_init(theta, data, prior, warmup)
        end
        @debug "Initialized NUTS sampler for survival submodel (G = $(G))"

        ##@info "Initializing NUTS sampler for LMM submodel"
        ##lmm[G] = NUTS_JLCM_lmm_init(theta, data, prior)

        cov[G] = cov_model_init(theta, data, prior, warmup)
        @debug "Initialized NUTS sampler for covariance model (G = $(G))"

        ##@info "Initializing NUTS sampler for JLCM model (G = $(G))"
        ##cov[G] = NUTS_JLCM_init(theta, data, prior)
    end

    cp, surv, cov
end

function sampling(M::Int, data::Data, par_type::Type{Tp}, prior::Tprior=Prior(data), prop::Tprior=prior;
                  chains=1, G=1, warmup=Int(floor(M/2)), common_delta=false, common_beta=false,
                  return_warmup=false) where {Tp <: AbstractParam, Tprior <: AbstractPrior}
    control = Control(M, warmup, chains, return_warmup, common_delta, common_beta, 0.8)
    sampling(control, data, par_type, prior, prop; G=G)
end

function sampling(control::Control, data::Data, par_type::Type{Tp}, prior::Tprior=Prior(data), prop::Tprior=prior; G=1) where {Tp <: AbstractParam, Tprior <: AbstractPrior}

    sampling(control, data, [Tp(data,G) for k in 1:control.chains], prior, prop)
end

function sampling(M, data::Data, x::T; warmup=floor(Int, M/2), from_init=true) where T <: AbstractFit
    c = x.control
    ctrl = Control(M, warmup, c.chains, c.save_warmup, c.common_delta, c.common_beta, c.adapt_delta)

    if from_init
        init = [xx.init for xx in x.chains]
    else
        init = [xx.theta[end] for xx in x.chains]
    end

    sampling(ctrl, data, init, x.chains[1].prior)
end

function sampling(control::Control, data::Data, init::Vector{Tp}, prior::Tprior=Prior(data), prop::Tprior=prior) where {Tp <: AbstractParam, Tprior <: AbstractPrior}

    @assert length(init) == control.chains

    chains = progress_pmap(1:control.chains) do k
        @info "Chain $(k) of $(control.chains)"
        sampling(control, data, init[k]; prior=prior, prop=prop)
    end

    Fit(control, chains)
end

function sampling(control::Control, data::Data, init::Tp=Param(data, 1);
                  prior::Tpr=Prior(), prop::Tpr=prior) where {Tp <: AbstractParam, Tpr <: AbstractPrior}
    G = init.G
    M = control.iter
    warmup = control.warmup

    offset = control.save_warmup ? 0 : warmup

    @info "Taking $(M) samples including $(warmup) warmup iterations in $(G)-class model"

    loglik_ar = Matrix{Float64}(undef, M-offset, length(data.T))
    theta_stor = Vector{Tp}(undef, M-offset)

    theta = deepcopy(init)

    if offset == 0
        theta_stor[1] = deepcopy(theta)
        loglik_ar[1,:] = loglik_p(theta, data)
    end

    accept = Tp(data, G; set_zero=true)

    cp_sampler, surv_sampler, cov_sampler = init_NUTS_samplers(warmup, G, theta, data, prior; common_delta=control.common_delta)
    ac = 0.0

    ##@showprogress 1 "Sampling..."
    for t in 2:M
        if G > 1
            update_cp_params!(theta, accept, data, prior, cp_sampler[G])
        end

        if control.common_beta
            update_lmm_params_common!(theta, accept, data, prior)
        else
            update_lmm_params_mix!(theta, accept, data, prior)
        end

        theta, ac = cov_model_step!(cov_sampler[G], theta)
        accept.S .+= ac
        accept.L .+= ac
        accept.sigma += ac

        if control.common_delta
            theta, ac = surv_model_step2!(surv_sampler[G], theta)
        else
            theta, ac = surv_model_step!(surv_sampler[G], theta)
        end

        accept.delta .+= ac
        accept.bl_shape .+= ac
        accept.bl_scale .+= ac

        theta.G0 = sum([sum(theta.cind .== g) == 0 for g in 1:G])

        if t > offset
            ll = loglik_p(theta, data)
            loglik_ar[t-offset,:] = ll
            theta_stor[t-offset] = deepcopy(theta)
        end
    end

    Chain{Tp,Tpr}(theta_stor, accept, loglik_ar, prior, init)
end


text:
using ..Fed: STATS

"""
    QDiffStaticPayloadSerde(chunksize, is_client)

Applies both quantization and diff-deduplication to the weights before
serialization in order to compress the payload. Quantization is a lossy process,
thus a lost of information is to be expected after deserialization.
"""
mutable struct QDiffStaticPayloadSerde{T <: Unsigned} <: PayloadSerde
    chunksize::Int
    is_client::Bool
    quantizer::Quantizer{T}
    original_data::Vector{T}

    QDiffStaticPayloadSerde{T}(chunksize, is_client) where T <: Unsigned = new(
        chunksize,
        is_client,
        Quantizer{T}(-1.0f0, 1.0f0),
        Vector{Vector{T}}(undef, 0)
    )
end

"""
    serialize_payload(::QDiffPayloadSerde, weights)

Serializes `weights` with the `QuantizedDedupPayloadSerde` where quantization
and diff-deduplication of values are applied before serialization.
"""
function serialize_payload(
    p::QDiffStaticPayloadSerde{T},
    weights::Vector{Float32}
)::Vector{UInt8} where T <: Unsigned
    # quantize weights
    qweights = [quantize(p.quantizer, w) for w in weights]
    STATS.base.req_data = qweights

    if p.is_client
        patched_file = diff(p.original_data, qweights, sizeof(T))
    else
        p.original_data = qweights
        patched_file = Patch{T}(UnitRange[], qweights)
    end

    payload = QDiffPayload{T}(patched_file, p.quantizer.minval, p.quantizer.maxval)

    return pack(payload)
end

"""
    deserialize_payload(::QuantizedDedupPayloadSerde, data, from)

Deserializes `data` with the `QuantizedDedupPayloadSerde` where regeneration of
deduplicated values and dequantizatoin are applied after deserialization.
"""
function deserialize_payload(
    p::QDiffStaticPayloadSerde{T},
    data::Vector{UInt8},
    from::String
)::Vector{Float32} where T <: Unsigned
    payload = unpack(data)

    if p.is_client
        qweights = payload.data.data
        p.original_data = qweights
    else
        delta_earned = sum(length.(payload.data.range)) * sizeof(T) - (8 * length(payload.data.range))
        STATS.network.delta_earned += delta_earned

        qweights = patch(p.original_data, payload.data)
        push!(STATS.base.res_data, qweights)
    end

    # dequantize weights
    weights = [dequantize(p.quantizer, w) for w in qweights]

    return weights
end


text:
# This file is part of GenericMRRR.jl, released under the MIT "Expat" license

# The methods in this file are translated from
# the dqds routines in LAPACK.
# LAPACK authors are (obviously) not responsible for translation errors.

# LAPACK is released under a BSD license, and is
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.

module DQDS
using LinearAlgebra
using ..GenericMRRR: safemin

# Some of the dqds implementation here assumes over/underflow control and
# NaN handling as in IEEE-754.
# Here such code is gated with a trait, to allow for generalization.

abstract type FPModel end
struct IEEEFPModel <: FPModel end
struct UnknownFPModel <: FPModel end

# Ideally the default should be unknown, but that would make things annoying
# for Quadmath, ArbFloats, etc.
# FPModel(::Type{T}) where T = UnknownFPModel()
FPModel(::Type{T}) where T = IEEEFPModel()

# Open questions: DoubleFloats
FPModel(::Type{T}) where T <: Base.IEEEFloat = IEEEFPModel()
FPModel(::Type{BigFloat}) = IEEEFPModel()

# The bulk of the code works on a two-fold interleaving for data locality.
# We use (q,e) to update (qq,ee) and vice versa - this is the ping-pong scheme.
# This should be storage-equivalent to the array used in the LAPACK routines.
struct QDtet{T}
  q::T
  qq::T
  e::T
  ee::T
end

# We handle the ping-pong alternation with value type params.
# Values correspond to pp in LAPACK routines.

# The compiler might handle these cleanly.
@inline function update!(zv::AbstractVector{QDtet{T}}, j, ::Val{0}, qnew, enew) where T
    @inbounds zv[j] = QDtet(zv[j].q, qnew, zv[j].e, enew)
end
@inline function update!(zv::AbstractVector{QDtet{T}}, j, ::Val{1}, qnew, enew) where T
    @inbounds zv[j] = QDtet(qnew, zv[j].qq, enew, zv[j].ee)
end
@inline function update0!(zv::AbstractVector{QDtet{T}}, j, qnew, enew) where T
    @inbounds zv[j] = QDtet(zv[j].q, qnew, zv[j].e, enew)
end
@inline function update1!(zv::AbstractVector{QDtet{T}}, j, qnew, enew) where T
    @inbounds zv[j] = QDtet(qnew, zv[j].qq, enew, zv[j].ee)
end
@inline function _save_eig!(zv::AbstractVector{QDtet{T}}, j, w) where T
    @inbounds zv[j] = QDtet(w, zv[j].qq, zv[j].e, zv[j].ee)
end

mutable struct DQDS_State{T}
    nfail::Int
    iter::Int
    ndiv::Int
    ttype::Int
    dmin1::T
    dmin2::T
    dn::T
    dn1::T
    dn2::T
    δσ::T
    g::T
    τ::T
    function DQDS_State{T}(iter, ndiv) where T
        z0 = zero(T)
        new{T}(0,iter,ndiv,0,z0,z0,z0,z0,z0,z0,z0,z0)
    end
end

function _dqds_tol(n, ::Type{T}) where T
    log(T(n)) * 4 * eps(T)
end

# Use dqds algorithm to compute all eigenvalues of a tridiagonal SPD matrix A
# to high relative accuracy.
# Elements of the LU factorization of a tridiagonal similar to A must be
# interleaved into the qd vector Z on entry.
# Translation of DLASQ2 from LAPACK
#
# On algorithm failure DLASQ2 returns with a flag but I found no recovery logic,
# so here we just throw.
function _dqds_eigvals!(n, Z::AbstractVector{T}) where T
    z0 = zero(T)
    safmin = safemin(T)
    tol = 100 * eps(T)
    tol2 = tol^2

    simple_result(tr1, tr2) = (tr1, tr2, 0, 0.0, 0.0)

    if length(Z) < 2*n
        throw(ArgumentError("Z must be have at least 2n elements"))
    end
    if n == 0
        return simple_result(z0, z0)
    end
    if n == 1
        if Z[1] < 0
            throw(ArgumentError("Z must be PSD"))
        end
        return simple_result(Z[1], Z[1])
    end
    if n == 2
        if Z[2] < 0 || Z[3] < 0
            throw(ArgumentError("Z must be PSD"))
        end
        Ztrace = Z[1]+Z[2]+Z[3]
        w1, w2 = _qd_2x2(Z[1],Z[2],Z[3],tol2)
        Z[1] = w1
        Z[2] = w2
        esum = w1 + w2
        return simple_result(Ztrace, esum)
    end

    # check for negative data, compute sums of q's and e's
    # Note: paper says q,e must be positive but DLASQ2 allows zeros.
    Z[2n] = z0
    emin = Z[2]
    qmax = z0
    zmax = z0
    d = z0
    e = z0
    for k in 1:2:2*(n-1)
        if Z[k] < 0
            throw(ArgumentError("Z must be PSD: diag elt $k"))
        elseif Z[k+1] < 0
            throw(ArgumentError("Z must be PSD: off-diag elt $k"))
        end
        d += Z[k]
        e += Z[k+1]
        qmax = max(qmax, Z[k])
        emin = min(emin, Z[k+1])
        zmax = max(qmax, zmax, Z[k+1])
    end
    if Z[2*n-1] < 0
            throw(ArgumentError("Z must be PSD: diag elt 2n-1"))
    end
    d += Z[2*n-1]
    qmax = max(qmax, Z[2*n-1])
    zmax = max(qmax, zmax)

    Ztrace = d + e

    # check for diagonality
    if e == 0
        for k in 2:n
            Z[k] = Z[2*k-1]
        end
        sort!(view(Z,1:n), rev=true)
        Z[2*n-1] = d
        return simple_result(Ztrace, Ztrace)
    end

    if Ztrace == 0
        Z[2*n-1] = z0
        return simple_result(Ztrace, z0)
    end

    # rearrange data for locality
    # Flip the qd array if warranted.
    cbias = T(3/2)
    zv = Vector{QDtet{T}}(undef,n)
    if cbias * Z[1] < Z[2*n-1]
        for k in 1:n-1
            zv[k] = QDtet(Z[2(n+1-k)-1],z0,Z[2*(n-k)],z0)
        end
        zv[n] = QDtet(Z[1],z0,z0,z0)
    else
        for k in 1:n
            zv[k] = QDtet(Z[2k-1],z0,Z[2k],z0)
        end
    end

    i0 = 1
    n0 = n

    # initial split checking via dqd and Li's test
    for pp in 0:1
        qmax, emin = _split_check0!(zv, i0, n0, tol2, Val(pp))
    end

    # set up for and then loop over calls to _dqds_3 (good dqds steps)

    state = DQDS_State{T}(2,2*(n0-i0))

    array_done = false
    σ = z0
    for iwhila in 1:n+1
        if n0 < 1
            array_done = true
            break
        end
        # E[n0] holds the value of σ when submatrix i0:n0
        # splits from the rest, but is negated as a flag.
        state.δσ = z0 # reset for this iteration
        if n0 == n
            σ = z0
        else
            σ = -zv[n0].e
        end
        if σ < 0
            # This is a logic err, so panic is appropriate.
            throw(ErrorException("split marked by positive value in E"))
        end

        # find last unreduced submatrix's top index i0, qmax, and emin.
        # Find Geršgorin-type bound if Q's are much greater than E's.

        emax = z0
        if n0 > i0
            emin = abs(zv[n0-1].e)
        else
            emin = z0
        end
        qmin = zv[n0].q
        qmax = qmin
        early = false
        i4x = 0
        for i in n0:-1:2
            i4x = i
            if zv[i-1].e <= 0
                early = true
                break
            end
            if qmin >= 4*emax
                qmin = min(qmin, zv[i].q)
                emax = max(emax, zv[i-1].e)
            end
            qmax = max(qmax, zv[i-1].q + zv[i-1].e)
            emin = min(emin, zv[i-1].e)
        end
        if !early
            i4x = 1
        end
        i0 = i4x
        pp = 0
        if n0 - i0 > 1
            dee = zv[i0].q
            deemin = dee
            kmin = i0
            for i in i0:n0-1
                dee = zv[i+1].q * (dee / (dee + zv[i+1].e))
                if dee <= deemin
                    deemin = dee
                    kmin = i
                end
            end
            # flip if appropriate
            if ((kmin-i0)*2 < n0-kmin) && (deemin <= zv[n0].q * (1/T(2)))
                ipn4 = 4*(i0+n0)
                pp = 2
                ip = i0+n0
                for j in i0:(i0+n0-1) >> 1
                    zd0 = QDtet(zv[ip-j].q, zv[ip-j].qq, zv[ip-j-1].e, zv[ip-j-1].ee)
                    zv[ip-j-1] = QDtet(zv[ip-j-1].q, zv[ip-j-1].qq, zv[j].e, zv[j].ee)
                    zv[ip-j] = QDtet(zv[j].q, zv[j].qq, zv[ip-j].e, zv[ip-j].ee)
                    zv[j] = zd0
                end
            end
        end
        # put -(initial shift) into dmin
        # (use lower bound on Geršgorin if positive)
        dmin = -max(z0, qmin-2*sqrt(qmin)*sqrt(emax))

        # now i0:n0 is unreduced
        # use pp ∈ [0,1] for ping-pong; 2 indicates that Z was flipped
        # so tests for deflation in DLASQ3 should not be performed.

        nbig = 100 * (n0-i0+1)
        submatrix_done = false
        for iwhilb in 1:nbig
            if i0 > n0
                submatrix_done = true
                break
            end
            # While submatrix is unfinished take a good dqds step.
            n0,pp,dmin,σ = _dqds_3!(zv,i0,n0,pp,qmax,σ,dmin,state)

            pp = 1 - pp
            # When emin is very small check for splits
            if pp == 0 && (n0-i0 >= 3)
                if zv[n0].ee <= tol2 * qmax || zv[n0].e < tol2 * σ
                    qmax, emin, i0 = _split_check!(zv,i0,n0,tol2,σ)
                end

            end
        end # iwhilb
        if !submatrix_done
            # Iteration limit reached:
            # Restore shift σ and place new d's and e's in qd array.
            # This might need to be done for serveral blocks.
            i1 = i0
            n1 = n0
            blocks_remain = true
            while blocks_remain
                tempq = zv[i0].q
                qnew = zv[i0].q + σ
                zv[i0] = QDtet(qnew,zv[i0].qq,zv[i0].e,zv[i0].ee)
                for k in i0+1:n0
                    tempe = zv[k-1].e
                    enew = zv[k-1].e * (tempq / zv[k-1].q)
                    tempq = zv[k].q
                    qnew = tempq + σ + tempe - enew
                    zv[k-1] = QDtet(zv[k-1].q,zv[k-1].qq,enew,zv[k-1].ee)
                    zv[k] = QDtet(qnew,zv[k].qq,zv[k].e,zv[k].ee)
                end
                # prepare to repeat on previous block if any
                if i1 > 1
                    n1 = i1-1
                    while (i1 >= 2) && (zv[i1-1].e >= 0)
                        i1 -= 1
                    end
                    σ = -zv[n1].e
                else
                    blocks_remain = false
                end
            end
            # convergence failure; store what we have
            for k in 1:n
                Z[2k-1] = zv[k].q
                # Only block 1..n0 is unfinished.
                # Remainder of e's must be negligible, but junk might be
                # left in workspace.
                if k < n0
                    Z[2k] = zv[k].e
                else
                    Z[2k] = z0
                end
            end
            throw(ErrorException("inner loop failed to converge"))
            # info = 2
            return simple_result(Ztrace, z0)
        end
    end # iwhila
    if !array_done
        throw(ErrorException("outer loop failed"))
        # info = 3
        # return
    end

    # move q's to the front
    for k in 1:n
        Z[k] = zv[k].q
    end
    # sort and compute sum of eigenvalues
    sort!(view(Z,1:n),rev=true)
    esum = z0
    for k in n:-1:1
        esum += Z[k]
    end
    # return trace, sum, and performance monitors
    return Ztrace, esum, state.iter, state.ndiv / n^2, 100 * state.nfail / state.iter
end

for (pp,qf,ef,qqf,eef) in
     ((0, :q, :e, :qq, :ee),
      (1, :qq, :ee, :q, :e))
    @eval begin
        function _split_check0!(zv::AbstractVector{QDtet{T}}, i0, n0, tol2, ::Val{$pp}) where T
            z0 = T(0)
            safmin = safemin(T)
            d = zv[n0].$qf
            for i in (n0-1):-1:i0
                if zv[i].$ef <= tol2*d
                    d = zv[i].$qf
                    update!(zv,i,Val(1-$pp),d,-z0)
                else
                    d = zv[i].$qf * (d / (d + zv[i].$ef))
                end
            end

            # dqd maps Z to ZZ plus Li's test.
            emin = zv[i0+1].$qf # Z[4*i0+pp+1]
            d = zv[i0].$qf # Z[4*i0+pp-3]
            for i in i0:n0-1
                #for i4 in 4*i0+pp:4:4*(n0-1)+pp
                qnew = d + zv[i].$ef
                if zv[i].$ef < tol2*d
                    update!(zv,i,Val(1-$pp), zv[i].q, -z0)
                    qnew = d
                    enew = z0
                    d = zv[i+1].$qf
                elseif ((safmin * zv[i+1].$qf < qnew)
                        && (safmin * qnew < zv[i+1].$qf))
                    temp = zv[i+1].$qf / qnew
                    enew = zv[i].$ef * temp
                    d *= temp
                else
                    enew = zv[i+1].$qf * (zv[i].$ef / qnew)
                    d - zv[i+1].$qf * (d / qnew)
                end
                emin = min(emin, enew)
                update!(zv,i,Val($pp), qnew, enew)
            end
            update!(zv,n0,Val($pp), d, z0)
            # find qmax
            qmax = zv[i0].$qqf
            for i in i0+1:n0
                qmax = max(qmax, zv[i0].$qqf)
            end
            return qmax, emin
        end
    end # eval
end # pp loop

# compute eigenvalues of 2x2 sym. PSD matrix in qd form
@inline function _qd_2x2(q1::T, e1::T, q2::T, tol2) where T
    if q2 > q1
        q2, q1 = q1, q2
    end
    t = ((q1 - q2) + e1) * (1/T(2))
    if e1 > q2*tol2 && t != 0
        s = q2 * (e1 / t)
        if s <= t
            s = q2 * (e1 / (t * (1 + sqrt(1 + s/t))))
        else
            s = q2 * (e1 / (t + sqrt(t) * sqrt(t + s)))
        end
        t = q1 + (s + e1)
        q2 = q2 * (q1 / t)
        q1 = t
    end
    return q1, q2
end

function _split_check!(zv::AbstractArray{QDtet{T}},i0,n0,tol2,σ) where T
    z0 = T(0)
    splt = i0-1
    qmax = zv[i0].q
    emin = zv[i0].e
    old_emin = zv[i0].ee
    for i in i0:(n0-3)
        if zv[i].ee <= tol2 * zv[i].q || zv[i].e <= tol2 * σ
            #  zv[i].e = -σ
            zv[i] = QDtet(zv[i].q, zv[i].qq, -σ, zv[i].ee)
            splt = i
            qmax = z0
            emin = zv[i+1].e
            old_emin = zv[i+1].ee
        else
            qmax = max(qmax, zv[i+1].q)
            emin = min(emin, zv[i].e)
            old_emin = min(old_emin, zv[i].ee)
        end
    end
    zv[n0] = QDtet(zv[n0].q, zv[n0].qq, emin, old_emin)
    i0 = splt + 1
    return qmax, emin, i0
end

# n0,pp,dmin,σ = _dqds_3!(Z,i0,n0,pp,qmax,σ,dmin,state)
# Check for deflation, compute shift, and invoke dqds.
# In case of failure, change shifts and try again.
# translation of DLASQ3 from LAPACK
# Documentation for dlasq3 suggests it doesn't use old σ or dmin,
# but that's wrong.
# Deflation tests should be skipped if pp==2.
# That seems to depend on a confusing index in dlasq3; the logic is slightly
# different here but appears to work.
function _dqds_3!(zv::AbstractArray{QDtet{T}},i0,n0in,pp,qmax,σ,dmin,state) where T
    z0 = zero(T)
    myeps = eps(T)
    tol = 100myeps
    tol2 = tol^2
    cbias = T(3)/2

    n0 = n0in
    while true
        # check for deflation
        if n0 < i0
            return n0,pp,dmin,σ
        end

        if n0 == i0
            # deflate 1
            if pp==0
                w1 = zv[n0].q + σ
            elseif pp==1
                w1 = zv[n0].qq + σ
            else
                throw(ErrorException("simple deflation w/ pp=2: logic error"))
            end
            _save_eig!(zv, n0, w1)
            n0 -= 1
            continue
        end

        if n0 == i0+1
            deflate2 = true
        else
            if pp==0
                e1 = zv[n0-1].e
                ee1 = zv[n0-1].ee
                q1 = zv[n0-1].q
                w1 = σ + zv[n0].q
                e2 = zv[n0-2].e
                ee2 = zv[n0-2].ee
                q2 = zv[n0-2].q
            else
                e1 = zv[n0-1].ee
                ee1 = zv[n0-1].e
                q1 = zv[n0-1].qq
                w1 = σ + zv[n0].qq
                e2 = zv[n0-2].ee
                ee2 = zv[n0-2].e
                q2 = zv[n0-2].qq
            end
            if pp==2
                deflate2 = false
            else
                # if E[n0-1] is negligible, deflate 1
                if (e1 <= tol2 * w1) || (ee1 <= tol2 * q1)
                    _save_eig!(zv, n0, w1)
                    n0 -= 1
                    continue
                end
                deflate2 = (e2 <= tol2 * σ) || (ee2 <= tol2 * q2)
            end
        end
        if !deflate2
            break
        end

        if pp==0
            w1, w2 = _qd_2x2(zv[n0-1].q, zv[n0-1].e, zv[n0].q, tol2)
        else
            w1, w2 = _qd_2x2(zv[n0-1].qq, zv[n0-1].ee, zv[n0].qq, tol2)
        end
        w1 += σ
        w2 += σ
        _save_eig!(zv,n0,w1)
        _save_eig!(zv,n0-1,w2)
        n0 -= 2
        continue
    end # deflation loop

    if pp == 2
        pp = 0
    end
    # reverse qd array if warranted
    if dmin <= 0 || n0 < n0in
        flipped, xdmin2, xqmax = _maybe_flip!(zv, i0, n0, Val(pp))
        if flipped
            qmax = max( qmax, xqmax)
            state.dmin2 = min( state.dmin2, xdmin2)
            dmin = -z0
        end
    end # if reversing

    # choose a shift
    _dqds_4!(state, zv, i0, n0, pp, n0in, dmin)

    underflag = false
    while true
        if n0 > i0+1
            dmin = _dqds_5!(state,zv,i0,n0,pp,σ,myeps)
        end
        state.ndiv += n0-i0+2
        state.iter += 1
        if dmin >= 0 && state.dmin1 >= 0
            # success
            break
        end
        etest = (pp == 0) ? zv[n0-1].ee : zv[n0-1].e
        if   ((dmin < 0) && (state.dmin1 > 0)
              && (etest < tol * (σ + state.dn1))
              && (abs(state.dn) < tol*σ)
              )
            # convergence hidden by negative dn
            if pp==0
                zv[n0] = QDtet(zv[n0].q, z0, zv[n0].e, zv[n0].ee)
            else
                zv[n0] = QDtet(z0, zv[n0].qq, zv[n0].e, zv[n0].ee)
            end
            dmin = z0
            break
        elseif dmin < 0
            # τ too big; select another and try again
            state.nfail += 1
            if state.ttype < -22
                # failed twice; play it safe.
                state.τ = z0
            elseif state.dmin1 > 0
                # late failure: gives excellent shifts
                state.τ = (state.τ + dmin) * (1 - 2*myeps)
                state.ttype -= 11
            else
                # early failure
                state.τ /= 4
                state.ttype -= 12
            end
            continue
        elseif isnan(dmin)
            if state.τ == 0
                underflag = true
                break
            else
                state.τ = z0
                continue
            end
        else
            underflag = true
            break
        end
    end # dqds loop
    if underflag
        dmin = _dqd_6!(state, zv,i0,n0,Val(pp))
        state.ndiv += n0-i0+2
        state.iter += 1
        state.τ = z0
    end

    if state.τ < σ
        state.δσ += state.τ
        t = σ + state.δσ
        state.δσ -= (t - σ)
    else
        t = σ + state.τ
        state.δσ = σ - (t - state.τ) + state.δσ
    end
    σ = t
    return n0,pp,dmin,σ
end

function _maybe_flip!(zv::AbstractVector{QDtet{T}}, i0, n0, pp) where T
    z0 = T(0)
    cbias = T(3)/T(2)
    if pp==0
        flipping = cbias * zv[i0].q < zv[n0].q
    else
        flipping = cbias * zv[i0].qq < zv[n0].qq
    end
    if flipping
        ip = i0+n0
        for j in i0:(i0+n0-1) >> 1
            zd0 = QDtet(zv[ip-j].q, zv[ip-j].qq, zv[ip-j-1].e, zv[ip-j-1].ee)
            zv[ip-j-1] = QDtet(zv[ip-j-1].q, zv[ip-j-1].qq, zv[j].e, zv[j].ee)
            zv[ip-j] = QDtet(zv[j].q, zv[j].qq, zv[ip-j].e, zv[ip-j].ee)
            zv[j] = zd0
        end
        if n0-i0 <= 4
            zv[n0] = QDtet(zv[n0].q, zv[n0].qq, zv[i0].e, zv[i0].ee)
        end
        if pp==0
            dmin2 = zv[n0].e # Z[4*n0+pp-1]
            qmax = max( zv[i0].q, zv[i0+1].q)
        else
            dmin2 = zv[n0].ee # Z[4*n0+pp-1]
            qmax = max( zv[i0].qq, zv[i0+1].qq)
        end
        e1 = min( zv[n0].e, zv[i0].e, zv[i0+1].e)
        ee1 = min( zv[n0].ee, zv[i0].ee, zv[i0+1].ee)
        zv[n0] = QDtet(zv[n0].q, zv[n0].qq, e1, ee1)
    else
        dmin2, qmax = z0, z0
    end
    return flipping, dmin2, qmax
end

# compute an approximation to the smallest eigenvalue using values of d
# from previous transform.
# modifies `state`: sets τ, ttype;  updates g
# translation of DLASQ4 from LAPACK

# NOTE: Fernando (note, Feb 2020) claims that one should assign `state.τ = s`
# before all of the unusual return statements.
# I haven't found the rationale for these unusual branches, so I'm puzzled.
function _dqds_4!(state, zv::AbstractArray{QDtet{T}}, i0, n0, pp, n0in, dmin) where T
    z0 = zero(T)
    const1 = 0.563
    const2 = 1.01
    const3 = 1.04
    if dmin <= 0
        state.τ = -dmin
        state.ttype = -1
        return
    end
    if pp == 0
        q0 =  zv[n0].q
        q1 = zv[n0-1].q
        q2 = zv[n0-2].q
        e1 = zv[n0-1].e
        e2 = zv[n0-2].e
        if n0-i0>2
            q3 = zv[n0-3].q
            e3 = zv[n0-3].e
        end
        ee1 = zv[n0-1].ee
        ee2 = zv[n0-2].ee
    else
        q0 =  zv[n0].qq
        q1 = zv[n0-1].qq
        q2 = zv[n0-2].qq
        e1 = zv[n0-1].ee
        e2 = zv[n0-2].ee
        if n0-i0>2
            q3 = zv[n0-3].qq
            e3 = zv[n0-3].ee
        end
        ee1 = zv[n0-1].e
        ee2 = zv[n0-2].e
    end

    if n0in == n0
        # no eigenvalues deflated
        if dmin == state.dn || dmin == state.dn1
            b1 = sqrt(q0) * sqrt(e1)
            b2 = sqrt(q1) * sqrt(e2)
            a2 = q1 + e1

            # cases 2 and 3
            if dmin == state.dn && state.dmin1 == state.dn1
                gap2 = state.dmin2 - a2 - state.dmin2 / 4
                if gap2 > 0 && gap2 > b2
                    gap1 = a2 - state.dn - (b2 / gap2) * b2
                else
                    gap1 = a2 - state.dn - (b1 + b2)
                end
                if gap1 > 0 && gap1 > b1
                    s = max(state.dn - (b1 / gap1) * b1, dmin * (1/T(2)))
                    state.ttype = -2
                else
                    s = z0
                    if state.dn > b1
                        s = state.dn - b1
                    end
                    if a2 > (b1 + b2)
                        s = min(s, a2 - (b1+b2))
                    end
                    s = max( s, dmin / 3)
                    state.ttype = -3
                end
            else
                # case 4
                state.ttype = -4
                s = dmin / 4
                if dmin == state.dn
                    γ = state.dn
                    a2 = z0
                    if e1 > q1
                        return
                    end
                    b2 = e1 / q1
                    npp = n0 - 2
                else
                    γ = state.dn1
                    if pp==0
                        if zv[n0-1].ee > zv[n0-1].qq
                            return
                        end
                        a2 = zv[n0-1].ee / zv[n0-1].qq
                        if e2 > q2
                            return
                        end
                    else
                        if zv[n0-1].e > zv[n0-1].q
                            return
                        end
                        a2 = zv[n0-1].e / zv[n0-1].q
                        if e2 > q2
                            return
                        end
                    end
                    b2 = e2 / q2
                    npp = n0-3
                end
                # approximate contribution to norm squared from i < nn-1
                a2 += b2
                if pp==0
                    @inbounds for i in npp:-1:i0
                        if b2 == 0
                            break
                        end
                        b1 = b2
                        if zv[i].ee > zv[i].qq
                            return
                        end
                        b2 *= zv[i].ee / zv[i].qq
                        a2 += b2
                        if 100*max(b1,b2) < a2 || const1 < a2
                            break
                        end
                    end
                else
                    @inbounds for i in npp:-1:i0
                        if b2 == 0
                            break
                        end
                        b1 = b2
                        if zv[i].e > zv[i].q
                            return
                        end
                        b2 *= zv[i].e / zv[i].q
                        a2 += b2
                        if 100*max(b1,b2) < a2 || const1 < a2
                            break
                        end
                    end
                end
                a2 *= const3
                # Rayleigh quotient residual bound
                if a2 < const1
                    s = γ * (1 - sqrt(a2)) / (1 + a2)
                end
            end
        elseif dmin == state.dn2
            # case 5
            state.ttype = -5
            s = dmin / 4
            # compute contribution to norm squared from i > nn-2
            if pp == 0
                b1 = zv[n0].qq
                b2 = zv[n0-1].qq
                γ = state.dn2
                if zv[n0-2].ee > b2 || zv[n0-1].ee > b1
                    return
                end
                a2 = (zv[n0-2].ee / b2) * (1 + zv[n0-1].ee / b1)
            else
                b1 = zv[n0].q
                b2 = zv[n0-1].q
                γ = state.dn2
                if zv[n0-2].e > b2 || zv[n0-1].e > b1
                    return
                end
                a2 = (zv[n0-2].e / b2) * (1 + zv[n0-1].e / b1)
            end
            # approximate contribution from i < nn-2
            if n0-i0>2
                b2 = e3 / q3
                a2 += b2
                if pp == 0
                  @inbounds for i in n0-4:-1:i0
                    if b2 == 0
                        break
                    end
                    b1 = b2
                    if zv[i].e > zv[i].q
                        return
                    end
                    b2 *= zv[i].e / zv[i].q
                    a2 += b2
                    if 100 * max(b2, b1) < a2 || const1 < a2
                        break
                    end
                  end
                else
                  @inbounds for i in n0-4:-1:i0
                    if b2 == 0
                        break
                    end
                    b1 = b2
                    if zv[i].ee > zv[i].qq
                        return
                    end
                    b2 *= zv[i].ee / zv[i].qq
                    a2 += b2
                    if 100 * max(b2, b1) < a2 || const1 < a2
                        break
                    end
                  end
                end
                a2 *= const3
            end
            if a2 < const1
                s = γ * (1 - sqrt(a2)) / (1 + a2)
            end
        else
            # case 6, no information
            if state.ttype == -6
                state.g += (1 - state.g) / T(3)
            elseif state.ttype == -18
                state.g = 1 / T(12)
            else
                state.g = 1 / T(4)
            end
            s = state.g * dmin
            state.ttype = -6
        end
    elseif n0in == n0+1
        # one eigval deflated; use dmin1, dn1
        if state.dmin1 == state.dn1 && state.dmin2 == state.dn2
            # cases 7 and 8
            state.ttype = -7
            s = state.dmin1 * (1/T(3))
            if e1 > q1
                return
            end
            b1 = e1 / q1
            b2 = b1
            if b2 != 0
                if pp == 0
                    @inbounds for i in (n0-2):-1:i0

                        a2 = b1
                        if zv[i].e > zv[i].q
                            return
                        end
                        b1 *= zv[i].e / zv[i].q
                        b2 += b1
                        if 100 * max(b1,a2) < b2
                            break
                        end
                    end
                else
                    @inbounds for i in (n0-2):-1:i0
                        a2 = b1
                        if zv[i].ee > zv[i].qq
                            return
                        end
                        b1 *= zv[i].ee / zv[i].qq
                        b2 += b1
                        if 100 * max(b1,a2) < b2
                            break
                        end
                    end
                end
            end
            b2 = sqrt(const3 * b2)
            a2 = state.dmin1 / (1 + b2^2)
            gap2 = state.dmin2 * (1/T(2)) - a2
            if gap2 > 0 && gap2 > b2*a2
                s = max(s, a2 * (1-const2*a2*(b2/gap2)*b2))
            else
                s = max(s, a2*(1-const2*b2))
                state.ttype = -8
            end
        else
            # case 9
            s = state.dmin1 / 4
            if state.dmin1 == state.dn1
                s = state.dmin1 * (1/T(2))
            end
            state.ttype = -9
        end
    elseif n0in == n0+2
        # two eigvales deflated; use dmin2, dn2
        if state.dmin1 == state.dn1 && state.dmin2 == state.dn2
            # cases 10 and 11
            state.ttype = -10
            s = state.dmin1 / 3
            if e1 > q1
                return
            end
            b1 = e1 / q1
            b2 = b1
            if b2 != 0
                if pp == 0
                  @inbounds for i in (n0-2):-1:i0
                    if zv[i].e > zv[i].q
                        return
                    end
                    b1 *= zv[i].e / zv[i].q
                    b2 += b1
                    if 100 * b1 < b2
                        break
                    end
                  end
                else
                  @inbounds for i in (n0-2):-1:i0
                    if zv[i].ee > zv[i].qq
                        return
                    end
                    b1 *= zv[i].ee / zv[i].qq
                    b2 += b1
                    if 100 * b1 < b2
                        break
                    end
                  end
                end
            end
            b2 = sqrt(const3 * b2)
            a2 = state.dmin2 / (1 + b2^2)
            gap2 = q1 + e2 - sqrt(q2) * sqrt(e2) - a2
            if gap2 > 0 && gap2 > b2*a2
                s = max(s, a2 * (1-const2*a2*(b2/gap2)*b2))
            else
                s = max(s, a2*(1-const2*b2))
            end
        else
            # case 11
            s = state.dmin2 / 4
            state.ttype = -11
        end
    elseif n0in > n0+2
        # case 12: more than 2 eigvals deflated; no information
        s = z0
        state.ttype = -12
    end
    state.τ = s
end

# compute one dqds transform in ping-pong form
# "translation" of DLASQ5 from LAPACK
function _dqds_5!(state,zv::AbstractArray{QDtet{T}},i0,n0,pp,σ,ϵ) where T
    __dqds_5!(state,zv,i0,n0,Val(pp),σ,ϵ,FPModel(T))
end

# I could repeat myself, or I could spend a few hours trying to optimize some metacode...
for (pp, qf, ef, update) in
    ((0, :q, :e, :update0!),
     (1, :qq, :ee, :update1!))
 @eval begin
  function __dqds_5!(state,zv::AbstractArray{QDtet{T}},i0,n0,::Val{$pp},σ,ϵ,::IEEEFPModel) where T
    z0 = zero(T)
    dthresh = ϵ * (σ + state.τ)
    negl_τ = state.τ < dthresh * (1/T(2))
    if negl_τ
        state.τ = z0
    end
    @inbounds begin
        # we set d's to 0 if they are small enough
        # if branches are costly, consider dispatching on negl_τ
        emin = zv[i0+1].$qf
        d = zv[i0].$qf - state.τ
        dmin = d
        state.dmin1 = -zv[i0].$qf
        @inbounds for j in i0:(n0-3)
            qnew = d + zv[j].$ef
            t = zv[j+1].$qf / qnew
            d = d * t - state.τ
            if negl_τ && (d < dthresh)
                d = z0
            end
            dmin = min( dmin, d)
            enew = zv[j].$ef * t
            emin = min( enew, emin)
            $update(zv,j,qnew,enew)
        end

        # unroll last two steps
        state.dn2 = d
        state.dmin2 = dmin
        j = (n0-2)
        qnew = state.dn2 + zv[j].$ef
        enew = zv[j+1].$qf * (zv[j].$ef / qnew)
        state.dn1 = zv[j+1].$qf * (state.dn2 / qnew) - state.τ
        $update(zv,j,qnew,enew)
        dmin = min(dmin, state.dn1)

        state.dmin1 = dmin
        j += 1
        qnew = state.dn1 + zv[j].$ef
        enew = zv[j+1].$qf * (zv[j].$ef / qnew)
        state.dn = zv[j+1].$qf * (state.dn1 / qnew) - state.τ
        $update(zv,j,qnew,enew)
        dmin = min(dmin, state.dn)

        update!(zv,n0,Val($pp),state.dn,emin)
    end
      return dmin
  end
# compute one dqd (zero shift) transform in ping-pong form,
# with over/underflow handling.
# Translation of DLASQ6 from LAPACK
  function _dqd_6!(state, zv::AbstractArray{QDtet{T}},i0,n0,::Val{$pp}) where T
    z0 = zero(T)
    safmin = safemin(T)
    emin = zv[i0+1].$qf
    d = zv[i0].$qf
    dmin = d

    for j in i0:(n0-3)
        qnew = d + zv[j].$ef
        if qnew == 0
            enew = z0
            d = zv[j+1].$qf
            dmin = d
            emin = z0
        elseif (safmin * zv[j+1].$qf < qnew) && (safmin * qnew < zv[j+1].$qf)
            t = zv[j+1].$qf / qnew
            d *= t
            enew = zv[j].$ef * t
        else
            enew = zv[j+1].$qf * (zv[j].$ef / qnew)
            d = zv[j+1].$qf * (d / qnew)
        end
        update!(zv,j,Val($pp),qnew,enew)
        dmin = min(dmin, d)
        emin = min(emin, enew)
    end

    # unroll last two steps
    state.dn2 = d
    state.dmin2 = dmin
    j = (n0-2)
    qnew = state.dn2 + zv[j].$ef
    if qnew == 0
        enew = z0
        state.dn1 = zv[j+1].$qf
        dmin = state.dn1
        emin = z0
    elseif (safmin * zv[j+1].$qf < qnew) && (safmin * qnew < zv[j+1].$qf)
        t = zv[j+1].$qf / qnew
        enew *= t
        state.dn1 = state.dn2 * t
    else
        enew = zv[j+1].$qf * (zv[j].$ef / qnew)
        state.dn1 = zv[j+1].$qf * (state.dn2 / qnew)
    end
    update!(zv,j,Val($pp),qnew,enew)
    dmin = min(dmin, state.dn1)

    state.dmin1 = dmin
    j += 1
    qnew = state.dn1 + zv[j].$ef
    if qnew == 0
        enew = z0
        state.dn = zv[j+1].$qf
        dmin = state.dn
        emin = z0
    elseif (safmin * zv[j+1].$qf < qnew) && (safmin * qnew < zv[j+1].$qf)
        t = zv[j+1].$qf / qnew
        enew = zv[j].$ef * t
        state.dn = state.dn1 * t
    else
        enew = zv[j+1].$qf * (zv[j].$ef / qnew)
        state.dn = zv[j+1].$qf * (state.dn1 / qnew)
    end
    update!(zv,j,Val($pp),qnew,enew)
    dmin = min(dmin, state.dn)

    update!(zv,n0,Val($pp),state.dn,emin)

    return dmin
  end # function
 end # eval block
end # pp loop

end # module