text:
##==============================================================================
## LEGACY SUPPORT FOR ZMQ IN CAESAR
##==============================================================================

export listSolvekeys

export _evalType

# not sure if and where this is still being used
function _evalType(pt::String)::Type
  @error "_evalType has been deprecated, use DFG serialization methods instead."
  try
    getfield(Main, Symbol(pt))
  catch ex
    io = IOBuffer()
    showerror(io, ex, catch_backtrace())
    err = String(take!(io))
    error("_evalType: Unable to locate factor/distribution type '$pt' in main context (e.g. do a using on all your factor libraries). Please check that this factor type is loaded into main. Stack trace = $err")
  end
end

##==============================================================================
## TODO deprecated
##==============================================================================

# TODO Replace with _CalcFactorParametric API
function (s::MsgPrior{<:ParametricTypes})(X1::AbstractVector{T};
                        userdata::Union{Nothing,FactorMetadata}=nothing) where T<:Real

  if isa(s.Z, Normal)
    meas = s.Z.μ
    σ = s.Z.σ
    #TODO confirm signs
    res = meas - X1[1]
    return (res./σ) .^ 2

  elseif isa(s.Z, MvNormal)
    meas = mean(s.Z)
    iΣ = invcov(s.Z)
    #TODO confirm math : Σ^(1/2)*X
    res = meas .- X1
    return res' * iΣ * res

  else
    #this should not happen
    @error("$s not suported, please use non-parametric")
  end                    #
end

# maybe replace X with a type.
function (s::Prior{<:ParametricTypes})(X1::AbstractVector{T};
                    userdata::Union{Nothing,FactorMetadata}=nothing) where T <: Real

  if isa(s.Z, Normal)
    meas = s.Z.μ
    σ = s.Z.σ
    #TODO confirm signs
    res = meas - X1[1]
    return (res./σ) .^ 2

  elseif isa(s.Z, MvNormal)
    meas = mean(s.Z)
    iΣ = invcov(s.Z)
    #TODO confirm math : Σ^(1/2)*X
    res = meas .- X1
    return res' * iΣ * res # + 2*log(1/(  sqrt(det(Σ)*(2pi)^k) )) ## cancel ×1/2 in calling function ## k = dim(μ)
  else
    #this should not happen
    @error("$s not suported, please use non-parametric")
  end
end

# parametric specific functor
# TODO Replace with _CalcFactorParametric API
function (s::LinearRelative{N,<:ParametricTypes})(
                                X1::AbstractArray{<:Real},
                                X2::AbstractArray{<:Real};
                                userdata::Union{Nothing,FactorMetadata}=nothing ) where N
  #
  # can I change userdata to a keyword arg, DF, No will be resolved with consolidation, #467
  # FIXME, replace if with dispatch
  if isa(s.Z, Normal)
    meas = mean(s.Z)
    σ = std(s.Z)
    # res = similar(X2)
    res = meas - (X2[1] - X1[1])
    return (res/σ) .^ 2

  elseif isa(s.Z, MvNormal)
    meas = mean(s.Z)
    iΣ = invcov(s.Z)
    #TODO confirm math : Σ^(1/2)*X
    res = meas .- (X2 .- X1)
    return res' * iΣ * res

  else
    #this should not happen
    @error("$s not supported, please use non-parametric")
  end
end

# parametric specific functor
# TODO Replace with _CalcFactorParametric API
function (s::EuclidDistance{<:ParametricTypes})(X1::AbstractArray{<:Real},
                                                X2::AbstractArray{<:Real};
                                                userdata::Union{Nothing,FactorMetadata}=nothing )
  #
  # can I change userdata to a keyword arg, DF, No will be resolved with consolidation
  #
  # FIXME, replace if with dispatch
  if isa(s.Z, Normal)
    meas = mean(s.Z)
    σ = std(s.Z)
    # res = similar(X2)
    res = meas - norm(X2 - X1)
    res *= res
    return res/(σ^2)

  elseif isa(s.Z, MvNormal)
    meas = mean(s.Z)
    iΣ = invcov(s.Z)
    #TODO confirm math : Σ^(1/2)*X
    res = meas .- (X2 .- X1)
    return res' * iΣ * res

  else
    #this should not happen
    @error("$s not suported, please use non-parametric")
  end
end

"""
$(TYPEDEF)

TODO TO BE DEPRECATED
"""
mutable struct PotProd
    Xi::Symbol # Int
    prev::Array{Float64,2}
    product::Array{Float64,2}
    potentials::Array{BallTreeDensity,1}
    potentialfac::Vector{Symbol}
end

"""
$(TYPEDEF)

TODO TO BE DEPRECATED
"""
mutable struct CliqGibbsMC
    prods::Array{PotProd,1}
    lbls::Vector{Symbol}
    CliqGibbsMC() = new()
    CliqGibbsMC(a,b) = new(a,b)
end

"""
$(TYPEDEF)

TODO TO BE DEPRECATED
"""
mutable struct DebugCliqMCMC
  mcmc::Union{Nothing, Array{CliqGibbsMC,1}}
  outmsg::LikelihoodMessage
  outmsglbls::Dict{Symbol, Symbol} # Int
  priorprods::Vector{CliqGibbsMC}
  DebugCliqMCMC() = new()
  DebugCliqMCMC(a,b,c,d) = new(a,b,c,d)
end

##==============================================================================
## Deprecate code below before v0.22
##==============================================================================

"""
(cf::CalcFactor)( res::AbstractVector{<:Real}, meas..., params...)

Default fallback for the standard factor calculation interface, as in `cf.factor(residual, noise_process, parameters)`,
where factors are either library standard or user out-of-library factor definitions.  See documentation for
more details and tutorials on using your own factors (designed to be as easy as possible).

Notes
- These residual calculations use used to find non-Gaussian / multimodal (incl. PPE) and conventional Gaussian estimates.
- `cf.legacyMeas == (measparams[1:cf._measCount]...,)`

Example
```julia
# TBD
```
"""
function (cf::CalcFactor)(measparams... ) #where {T<:FunctorInferenceType,M,P<:Tuple,X<:AbstractVector} {T,M,P,X}
  #
  # NOTE this is a legacy interface
  res = zeros(size(cf._legacyMeas,1))
  cf.factor(res, cf.metadata, cf._sampleIdx, cf._legacyMeas, cf._legacyParams...)
end

@deprecate prodmultipleonefullpartials(w...;kw...) prodmultiplefullpartials(w...;kw...)

# """
#     $(SIGNATURES)

# Multiply a single full and several partial dimension constraints.

# DevNotes
# - FIXME consolidate partial and full product AMP API, relates to #1010
# - TODO -- reuse memory rather than rand here
# """
# function prodmultipleonefullpartials( dens::Vector{BallTreeDensity},
#                                       partials::Dict{Int, Vector{BallTreeDensity}},
#                                       Ndims::Int,
#                                       N::Int,
#                                       manis::Tuple  )
#   #

#   # TODO -- should this be [1.0] or ones(Ndims)
#   denspts = getPoints(dens[1])
#   pGM = deepcopy(denspts)

#   for (dimnum,pp) in partials
#     push!(pp, AMP.manikde!(pGM[dimnum:dimnum,:], (manis[dimnum],) ))
#   end

#   # do each partial dimension individually
#   for (dimnum,pp) in partials
#     pGM[dimnum,:] = AMP.manifoldProduct(pp, (manis[dimnum],), Niter=1) |> getPoints
#   end
#   return pGM
# end

"""
    $(SIGNATURES)

Multiply different dimensions from partial constraints individually.

DevNotes
- FIXME Integrate with `manifoldProduct`, see #1010
"""
function productpartials!(pGM::Array{Float64,2},
                          dummy::BallTreeDensity,
                          partials::Dict{Int, Vector{BallTreeDensity}},
                          manis::Tuple  )
  #
  @warn "productpartials! is being deprecated without direct replacement, see future versions of AMP.manifoldProduct instead."
  # do each partial dimension individually
  for (dimnum,pp) in partials
    pGM[dimnum,:] = AMP.manifoldProduct(pp, (manis[dimnum],), Niter=1) |> getPoints
  end
  nothing
end

@deprecate freshSamples!(w...;kw...) sampleFactor!(w...;kw...)
@deprecate freshSamples(w...;kw...) sampleFactor(w...;kw...)

@deprecate numericSolutionCCW!(w...;kw...) _solveCCWNumeric!(w...;kw...)

@deprecate solveFactorMeasurements( dfg::AbstractDFG,fctsym::Symbol,solveKey::Symbol=:default;retries::Int=3 ) approxDeconv(dfg,fctsym,solveKey,retries=retries)

##==============================================================================
## Deprecate code below before v0.21
##==============================================================================

# see DFG #590
@deprecate extractdistribution(x) convert(SamplableBelief, x)

"""
$(TYPEDEF)

MUST BE DEPRECATED IN FAVOR OF EuclidDistance
"""
mutable struct Ranged <: AbstractRelativeRoots
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    Ranged() = new()
    Ranged(x...) = new(x[1], x[2], x[3])
end
"""
$(TYPEDEF)
"""
mutable struct PackedRanged <: PackedInferenceType
    Zij::Array{Float64,1}
    Cov::Array{Float64,1}
    W::Array{Float64,1}
    PackedRanged() = new()
    PackedRanged(x...) = new(x[1], x[2], x[3])
end
function convert(::Type{Ranged}, r::PackedRanged)
  return Ranged(r.Zij, r.Cov, r.W)
end
function convert(::Type{PackedRanged}, r::Ranged)
  return PackedRanged(r.Zij, r.Cov, r.W)
end
function (ra::Ranged)(res::AbstractVector{<:Real},
    userdata::FactorMetadata,
    idx::Int,
    meas::Tuple{<:AbstractArray{<:Real,2}},
    p1::AbstractArray{<:Real},
    l1::AbstractArray{<:Real})

  res[1] = meas[1][1,idx] - abs(l1[1,idx] - p1[1,idx])
  nothing
end
function getSample(ra::Ranged, N::Int=1)
  @warn("`::Ranged` is being deprecated in favor of `::EuclidDistance`")
  ret = zeros(1,N)
  for i in 1:N
    ret[1,i] = ra.Cov[1]*randn()+ra.Zij[1]
  end
  # rand(Distributions.Normal(odo.Zij[1],odo.Cov[1]), N)'
  return (ret,)
end

function getSample( cfo::CalcFactor,
                    N::Int=1)
  #
  if !hasfield(typeof(cfo.factor), :specialSampler)
    @warn "`getSample(::MyFactor, ::Int)` API is being deprecated, use `getSample(cf::CalcFactor{<:$(typeof(cfo.factor).name)}, ::Int=1) = cf.factor...` instead.  Similarly for residual factor calculations, see IIF #467."
    getSample(cfo.factor, N)
  else
    @warn "`myfactor.specialSampler` API is being deprecated, use `getSample(cf::CalcFactor{<:$(typeof(cfo.factor).name)}, ::Int=1) = cf.metadata.___` instead.  Similarly for residual factor calculations, see IIF #467."
    cfo.factor.specialSampler(cfo.factor, N, cfo.metadata, cfo.metadata.fullvariables...)
  end
end

function freshSamples(usrfnc::T,
                      N::Int,
                      fmd::FactorMetadata,
                      vnd::Vector=[]) where { T <: FunctorInferenceType }
  #
  @warn("freshSamples(::FunctorInferenceType, ::Int, ::FactorMetadata) has been deprecated, use freshSamples!(::CommonConvWrapper, ::Int) instead.")
  if !hasfield(T, :specialSampler)
    getSample(usrfnc, N)
  else
    usrfnc.specialSampler(usrfnc, N, fmd, vnd...)
  end
end

function calcZDim(usrfnc::T,
                  Xi::Vector{<:DFGVariable},
                  fmd::FactorMetadata=FactorMetadata(Xi, getLabel.(Xi), Vector{Matrix{Float64}}(), :null, nothing) ) where {T <: FunctorInferenceType}
  #
  # # zdim = T != GenericMarginal ? size(getSample(usrfnc, 2)[1],1) : 0
  zdim = if T != GenericMarginal
    error("this version of calcZDim has been made obsolete, use calcZDim(::CalcFactor) instead.")
  #   vnds = Xi # (x->getSolverData(x)).(Xi)
  #   smpls = freshSamples(usrfnc, 2, fmd, vnds)[1]
  #   size(smpls,1)
  else
    0
  end
  return zdim
end

# # TODO make in-place memory version
# function getSample( s::Mixture{N_,F,S,T},
#                     N::Int=1,
#                     special...;
#                     kw...  ) where {N_,F<:FunctorInferenceType,S,T}
#   #
#   # TODO consolidate #927, case if mechanics has a special sampler
#   # TODO slight bit of waste in computation, but easiest way to ensure special tricks in s.mechanics::F are included
#   ## example case is old FluxModelsPose2Pose2 requiring velocity
#   smplLambda = hasfield(typeof(s.mechanics), :specialSampler) ? ()->s.specialSampler(s.mechanics, N, special...; kw...)[1] : ()->getSample(s.mechanics, N)[1]
#   smpls = smplLambda()
#     # smpls = Array{Float64,2}(undef,s.dims,N)
#   #out memory should be right size first
#   length(s.labels) != N ? resize!(s, N) : nothing
#   s.labels .= rand(s.diversity, N)
#   for i in 1:N
#     mixComponent = s.components[s.labels[i]]
#     smpls[:,i] = rand(mixComponent,1)
#   end

#   # TODO only does first element of meas::Tuple at this stage, see #1099
#   (smpls, ) # s.labels
# end

# # should not be called in case of Prior
# (s::Mixture)( res::AbstractArray{<:Real},
#               fmd::FactorMetadata,
#               idx::Int,
#               meas::Tuple,
#               X... ) = s.mechanics(res, fmd, idx, meas, X...)
#

# function (s::LinearRelative)( res::AbstractArray{<:Real},
#                               userdata::FactorMetadata,
#                               idx::Int,
#                               meas::Tuple,
#                               X1::AbstractArray{<:Real,2},
#                               X2::AbstractArray{<:Real,2}  )
#   #
#   res[:] = meas[1][:,idx] - (X2[:,idx] - X1[:,idx])
#   nothing
# end

# function (s::EuclidDistance)( res::AbstractArray{<:Real},
#                               fmd::FactorMetadata,
#                               idx::Int,
#                               meas::Tuple,
#                               X1::AbstractArray{<:Real,2},
#                               X2::AbstractArray{<:Real,2}  )
#   #
#   res[1] = meas[1][1,idx] - norm(X2[:,idx] - X1[:,idx])
#   res[1] ^= 2
#   return res[1]
# end

function PriorSphere1(mu::Array{Float64}, cov::Array{Float64,2}, W::Vector{Float64})
  @warn "PriorSphere1(mu,cov,W) is deprecated in favor of PriorSphere1(T(...)) -- use for example PriorSphere1(MvNormal(mu, cov))"
  PriorSphere1(MvNormal(mu[:], cov))
end

@deprecate fetchMsgUpThis(cliq::TreeClique) IIF.getMessageBuffer(cliq).upTx

#


text:
# dicom2nifti.jl
# Simple DICOM to NIfTI converter

# This script walks through the current working directory looking for DICOM
# files and sorts them into volumes. It then creates directories corresponding
# to each protocol, and saves nii files named by series number. It has
# (barely) been tested with DICOM files from a Siemens scanner. It may not
# work with other scanners. It may give incorrect information. Test thoroughly
# and use at your own risk.

using DICOM, NIfTI, ArgParse

function parse_commandline()
    s = ArgParseSettings()
    s.description = "DICOM to NIfTI converter"

    @add_arg_table s begin
        "--extension"
            help = "assume files with given extension are DICOM files. Set to blank for no extension."
            default = ".dcm"
        "--sphinx"
            help = "correct for sphinx position (for monkeys)"
            action = :store_true
        "--axes"
            help = """orientation of tensor in the resulting volume. This is a three-character label
            that determines how the data is stored, and how it is displayed in applications that do
            not make use of the orientation matrix in the NIfTI header. Valid characters are:
            L (= left), R (= right), A (= anterior), P (= posterior), I (= inferior), and
            S (= superior). If not specified, the orientation in the DICOM file is preserved."""
            default = nothing
        "dicomdir"
            help = "directory of DICOM files"
            required = true
        "niftidir"
            help = "directory in which to put NIfTI files"
            required = true
    end

    parse_args(s)
end

function walk(fn::Function, path::String)
    for fname in readdir(path)
        if fname == "." || fname == ".."
            continue
        end
        fpath = joinpath(path, fname)
        if isdir(lstat(fpath))
            walk(fn, fpath)
        else
            fn(fpath)
        end
    end
end

function findtag(d, tag, from, what)
    tag = lookup(d, tag)
    if tag == false
        error("$what missing from $from")
    end
    tag.data
end

function main()
    cmd = parse_commandline()

    if cmd["axes"] == nothing
        axes = nothing
    else
        axes = uppercase(cmd["axes"])
        if length(axes) != 3 || (!('R' in axes) && !('L' in axes)) ||
                (!('A' in axes) && !('P' in axes)) ||
                (!('S' in axes) && !('I' in axes))
            error("Invalid axes provided")
        end
    end

    transform = cmd["sphinx"] ? [1 0 0; 0 0 1; 0 -1 0] : [-1 0 0; 0 -1 0; 0 0 1]

    # Load all DICOMs into an array, indexed by series number
    dicoms = Dict{Int, Any}()
    walk(cmd["dicomdir"]) do fpath
        if !endswith(fpath, cmd["extension"])
            return
        end

        local d
        f = open(fpath, "r")
        d = DICOM.dcm_parse(f)

        series_number_index = findfirst((x) -> x.tag == (0x0020, 0x0011), d)
        if series_number_index == 0
            println("No series number found for $fpath; ignoring")
            return
        end
        series_number = d[series_number_index].data[1]

        dicom_arr = get(dicoms, series_number, nothing)
        if dicom_arr == nothing
            dicoms[series_number] = dicom_arr = {}
        end
        push!(dicom_arr, d)
    end

    # Sort DICOMs
    for (series_number, slices) in dicoms
        d = slices[1]

        protocol_name = findtag(d, (0x0018, 0x1030), series_number, "Protocol Name")[1]
        orientation = findtag(d, (0x0020, 0x0037), series_number,
            "Image Orientation (Patient)")
        pixel_spacing = findtag(d, (0x0028, 0x0030), series_number, "Pixel Spacing")
        slice_thickness = findtag(d, (0x0018, 0x0050), series_number, "Slice Thickness")[1]
        time_step = lookup(d, (0x0018, 0x0080)) # This is the TR, but may not be the time step
        phase_encoding_direction = lookup(d, (0x0018, 0x1312))

        # Convert orentation to RAS
        orientation = transform*reshape(orientation, (3, 2))

        # Determine Z as cross product of X and Y
        orientation = hcat(orientation, cross(orientation[:, 1], orientation[:, 2]))

        # Scale by pixel size
        orientation = orientation*[pixel_spacing[1] 0 0; 0 pixel_spacing[2] 0; 0 0 slice_thickness]

        # Find Z coordinates of each slice
        image_positions = transform*hcat(
            [float32(findtag(d, (0x0020, 0x0032), series_number, "Image Position (Patient)"))
            for d in slices]...)
        z_coords = [dot(image_positions[:, i], orientation[:, 3])
            for i = 1:size(image_positions, 2)]

        # Sort Z coordinates
        p = sortperm(z_coords)

        # Concatenate slices along Z axis according to increasing Z coordinate
        raw = cat(3,
            [findtag(slices[i], (0x7FE0, 0x0010), series_number, "Pixel Data")[1] for i in p]...)

        orientation = float32(hcat(orientation, image_positions[:, p[1]]))

        dim_info = phase_encoding_direction.data[1] == "ROW" ? (1, 2, 3) :
                phase_encoding_direction.data[1] == "COL" ? (2, 1, 3) :
                (0, 0, 0)

        voxel_size = tuple(pixel_spacing..., slice_thickness)

        # Permute according to axes specified on command line
        if axes != nothing
            # Determine permutation of current volume to RAS
            ras = Array(Int, 3)
            sign = Array(Bool, 3)
            rg = [1:3;]
            abs_orientation = abs(orientation)
            for i = 1:3
                idx = findmax(abs_orientation[i, rg])[2]
                ras[i] = splice!(rg, idx)
                sign[i] = orientation[i, ras[i]] >= 0
            end

            # Determine permutation of RAS to user-preferred orientation
            specified_axes = Int[(c in "RL" ? 1 :
                                  c in "AP" ? 2 : 3)
                                 for c in axes]
            specified_sign = Bool[c in "RAS" for c in axes]
            output_axes = ras[specified_axes]

            # Permute
            raw = permutedims(raw, output_axes)
            orientation[:, 1:3] = orientation[:, output_axes]
            voxel_size = voxel_size[output_axes]
            dim_info = dim_info[output_axes]

            # Flip signs and dimensions
            flip_sign = find(sign[specified_axes] .!= specified_sign)
            orientation[:, flip_sign] = -orientation[:, flip_sign]
            for dim in flip_sign
                raw = flipdim(raw, dim)
                orientation[dim, 4] += (sign[specified_axes[dim]] ? 1 : -1) *
                    (size(raw, dim) - 1)*voxel_size[dim]
            end
        end

        # Create a directory for each protocol
        protocol_dir = joinpath(cmd["niftidir"], replace(protocol_name, "/", ""))
        if !isdir(protocol_dir)
            mkdir(protocol_dir)
        end

        # Write NIfTI volumes
        ni = NIVolume(raw; voxel_size=voxel_size,
            orientation=orientation, dim_info=dim_info,
            time_step=time_step != false && !isempty(time_step.data) ? time_step.data[1] : 0f0)
        niwrite(joinpath(protocol_dir, "$(series_number).nii"), ni)
    end
end

main()

ext: jl


text:
# This function is very important for the computation of Floquet multipliers and checks that the eigensolvers compute the eigenvalues with largest modulus instead of their default behaviour which is with largest real part. If this option is not properly set, bifurcations of periodic orbits will be missed.
function checkFloquetOptions(eigls::AbstractEigenSolver)
	if eigls isa DefaultEig
		return @set eigls.which = abs
	elseif eigls isa EigArpack
		return setproperties(eigls; which = :LM, by = abs)
	elseif eigls isa EigArnoldiMethod
		return setproperties(eigls; which = ArnoldiMethod.LM(), by = abs)
	end
	if eigls isa EigKrylovKit
		return @set eigls.which = :LM
	end
end

# see https://discourse.julialang.org/t/uniform-scaling-inplace-addition-with-matrix/59928/5

####################################################################################################
# Computation of Floquet Coefficients for periodic orbit problems

"""
	floquet = FloquetQaD(eigsolver::AbstractEigenSolver)

This composite type implements the computation of the eigenvalues of the monodromy matrix in the case of periodic orbits problems (based on the Shooting method or Finite Differences (Trapeze method)), also called the Floquet multipliers. The method, dubbed Quick and Dirty (QaD), is not numerically very precise for large / small Floquet exponents when the number of time sections is large because of many matrix products. It allows, nevertheless, to detect bifurcations. The arguments are as follows:
- `eigsolver::AbstractEigenSolver` solver used to compute the eigenvalues.

If `eigsolver == DefaultEig()`, then the monodromy matrix is formed and its eigenvalues are computed. Otherwise, a Matrix-Free version of the monodromy is used.

!!! danger "Floquet multipliers computation"
    The computation of Floquet multipliers is necessary for the detection of bifurcations of periodic orbits (which is done by analyzing the Floquet exponents obtained from the Floquet multipliers). Hence, the eigensolver `eigsolver` needs to compute the eigenvalues with largest modulus (and not with largest real part which is their default behavior). This can be done by changing the option `which = :LM` of `eigsolver`. Nevertheless, note that for most implemented eigensolvers in the current Package, the proper option is set.
"""
struct FloquetQaD{E <: AbstractEigenSolver } <: AbstractFloquetSolver
	eigsolver::E
	function FloquetQaD(eigls::AbstractEigenSolver)
		eigls2 = checkFloquetOptions(eigls)
		return new{typeof(eigls2)}(eigls2)
	end
	FloquetQaD(eigls::FloquetQaD) = eigls
end
geteigenvector(eig::FloquetQaD, vecs, n::Union{Int, Array{Int64,1}}) = geteigenvector(eig.eigsolver, vecs, n)

function (fl::FloquetQaD)(J, nev; kwargs...)
	if fl.eigsolver isa DefaultEig
		# we build the monodromy matrix and compute the spectrum
		monodromy = MonodromyQaD(J)
	else
		# we use a Matrix Free version
		monodromy = dx -> MonodromyQaD(J, dx)
	end
	vals, vecs, cv, info = fl.eigsolver(monodromy, nev)
	Inf in vals && @warn "Detecting infinite eigenvalue during the computation of Floquet coefficients"

	# the `vals` should be sorted by largest modulus, but we need the log of them sorted this way
	logvals = log.(complex.(vals))
	I = sortperm(logvals, by = real, rev = true)
	return logvals[I], geteigenvector(fl.eigsolver, vecs, I), cv, info
end
####################################################################################################
# ShootingProblem
# Matrix free monodromy operator
function MonodromyQaD(JacSH::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}, du::AbstractVector) where {Tpb <: ShootingProblem, Tjacpb, Torbitguess, Tp}
	sh = JacSH.pb
	x = JacSH.x
	p = JacSH.par

	# period of the cycle
	T = getPeriod(sh, x)

	# extract parameters
	M = getMeshSize(sh)
	N = div(length(x) - 1, M)

	# extract the time slices
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)

	out = copy(du)

	for ii in 1:M
		# call the jacobian of the flow
		@views out .= evolve(sh.flow, Val(:SerialdFlow), xc[:, ii], p, out, sh.ds[ii] * T).du
	end
	return out
end

# Compute the monodromy matrix at `x` explicitely, not suitable for large systems
function MonodromyQaD(JacSH::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}) where {Tpb <: ShootingProblem, Tjacpb, Torbitguess, Tp}
	sh = JacSH.pb
	x = JacSH.x
	p = JacSH.par

	# period of the cycle
	T = getPeriod(sh, x)

	# extract parameters
	M = getMeshSize(sh)
	N = div(length(x) - 1, M)

	Mono = zeros(N, N)

	# extract the time slices
	xv = @view x[1:end-1]
	xc = reshape(xv, N, M)
	du = zeros(N)

	for ii in 1:N
		du[ii] = 1
		# call jacobian of the flow
		@views Mono[:, ii] .= evolve(sh.flow, Val(:SerialdFlow), xc[:, 1], p, du, T).du
		du[ii] = 0
	end

	return Mono
end

# Compute the monodromy matrix at `x` explicitely, not suitable for large systems
# it is based on a matrix expression of the Jacobian of the shooting functional. We thus
# just extract the blocks needed to compute the monodromy
function MonodromyQaD(JacSH::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}) where {Tpb <: ShootingProblem, Tjacpb <: AbstractMatrix, Torbitguess, Tp}
	J = JacSH.jacpb
	sh = JacSH.pb
	M = getMeshSize(sh)
	N = div(length(JacSH.x) - 1, M)
	mono = copy(J[1:N, 1:N])
	if M == 1
		return mono + I
	end
	tmp = similar(mono)
	r = N
	for ii = 1:M-1
		# mono .= J[r+1:r+N, r+1:r+N] * mono
		@views mul!(tmp, J[r+1:r+N, r+1:r+N], mono)
		mono .= tmp
		r += N
	end
	return mono
end

# This function is used to reconstruct the spatio-temporal eigenvector of the shooting functional sh
# at position x from the Floquet eigenvector ζ
@views function MonodromyQaD(::Val{:ExtractEigenVector}, sh::ShootingProblem, x::AbstractVector, par, ζ::AbstractVector)

	# period of the cycle
	T = getPeriod(sh, x)

	# extract parameters
	M = getMeshSize(sh)
	N = div(length(x) - 1, M)

	# extract the time slices
	xv = x[1:end-1]
	xc = reshape(xv, N, M)

	out = evolve(sh.flow, Val(:SerialdFlow), xc[:, 1], par, ζ, sh.ds[1] * T).du
	out_a = [copy(out)]

	for ii in 2:M
		# call the jacobian of the flow
		out .= evolve(sh.flow, Val(:SerialdFlow), xc[:, ii], par, out, sh.ds[ii] * T).du
		push!(out_a, copy(out))
	end
	return out_a
end
####################################################################################################
# PoincareShooting

# matrix free evaluation of monodromy operator
function MonodromyQaD(JacSH::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}, dx_bar::AbstractVector) where {Tpb <: PoincareShootingProblem, Tjacpb, Torbitguess, Tp}
	psh = JacSH.pb
	x_bar = JacSH.x
	p = JacSH.par

	M = getMeshSize(psh)
	Nm1 = div(length(x_bar), M)

	# reshape the period orbit guess into a Matrix
	x_barc = reshape(x_bar, Nm1, M)
	@assert length(dx_bar) == Nm1 "Please provide the right dimension to your matrix-free eigensolver, it must be $Nm1."

	xc = similar(x_bar, Nm1 + 1)
	outbar = copy(dx_bar)
	outc = similar(dx_bar, Nm1 + 1)

	for ii in 1:M
		E!(psh.section,  xc,  view(x_barc, :, ii), ii)
		dE!(psh.section, outc, outbar, ii)
		outc .= diffPoincareMap(psh, xc, p, outc, ii)
		# check to <outc, normals[ii]> = 0
		# println("--> ii=$ii, <out, normali> = ", dot(outc, sh.section.normals[ii]))
		dR!(psh.section, outbar, outc, ii)
	end
	return outbar

end

# matrix based formulation of monodromy operator, not suitable for large systems
# it is based on a matrix expression of the Jacobian of the shooting functional. We thus
# just extract the blocks needed to compute the monodromy
function MonodromyQaD(JacSH::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}) where {Tpb <: PoincareShootingProblem, Tjacpb <: AbstractMatrix, Torbitguess, Tp}
	J = JacSH.jacpb
	sh = JacSH.pb
	T = eltype(J)

	M = getMeshSize(sh)
	Nj = length(JacSH.x)
	N = div(Nj, M)

	if M == 1
		return I-J
	end

	mono = copy(J[N+1:2N, 1:N])
	tmp = similar(mono)
	r1 = mod(2N, Nj)
	r2 = N
	for ii = 1:M-1
		# mono .= J[r+1:r+N, r+1:r+N] * mono
		@views mul!(tmp, J[r1+1:r1+N, r2+1:r2+N], mono)
		mono .= tmp
		r1 = mod(r1 + N, Nj)
		r2 += N
	end
	# the structure of the functional imposes to take into account the sign
	sgn = iseven(M) ? one(T) : -one(T)
	mono .*= sgn
	return mono
end

# This function is used to reconstruct the spatio-temporal eigenvector of the shooting functional sh
# at position x from the Floquet eigenvector ζ
@views function MonodromyQaD(::Val{:ExtractEigenVector}, psh::PoincareShootingProblem, x_bar::AbstractVector, p, ζ::AbstractVector)
	#  ζ is of size (N-1)
	M = getMeshSize(psh)
	Nm1 = length(ζ)
	dx = similar(x_bar, length(ζ) + 1)

	x_barc = reshape(x_bar, Nm1, M)
	xc = similar(x_bar, Nm1 + 1)
	dx_bar = similar(x_bar, Nm1)
	outbar = copy(dx_bar)
	outc = similar(dx_bar, Nm1 + 1)
	out_a = typeof(xc)[]

	for ii in 1:M
		E!(psh.section,  xc,  view(x_barc, :, ii), ii)
		dE!(psh.section, outc, outbar, ii)
		outc .= diffPoincareMap(psh, xc, p, outc, ii)
		# check to <outc, normals[ii]> = 0
		# println("--> ii=$ii, <out, normali> = ", dot(outc, sh.section.normals[ii]))
		dR!(psh.section, outbar, outc, ii)
		push!(out_a, copy(outbar))
	end
	return out_a
end
####################################################################################################
# PeriodicOrbitTrapProblem

# Matrix-Free version of the monodromy operator
function MonodromyQaD(JacFW::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}, du::AbstractVector) where {Tpb <: PeriodicOrbitTrapProblem, Tjacpb, Torbitguess, Tp}
	poPb = JacFW.pb
	u0 = JacFW.x
	par = JacFW.par

	# extraction of various constants
	M, N = size(poPb)

	# period of the cycle
	T = extractPeriodFDTrap(poPb, u0)

	# time step
	h =  T * getTimeStep(poPb, 1)
	Typeh = typeof(h)

	out = copy(du)

	u0c = getTimeSlices(u0, N, M)

	@views out .= out .+ h/2 .* apply(poPb.J(u0c[:, M-1], par), out)
	# res = (I - h/2 * poPb.J(u0c[:, 1])) \ out
	@views res, _ = poPb.linsolver(poPb.J(u0c[:, 1], par), out; a₀ = convert(Typeh, 1), a₁ = -h/2)
	out .= res

	for ii in 2:M-1
		h =  T * getTimeStep(poPb, ii)
		@views out .= out .+ h/2 .* apply(poPb.J(u0c[:, ii-1], par), out)
		# res = (I - h/2 * poPb.J(u0c[:, ii])) \ out
		@views res, _ = poPb.linsolver(poPb.J(u0c[:, ii], par), out; a₀ = convert(Typeh, 1), a₁ = -h/2)
		out .= res
	end

	return out
end

# This function is used to reconstruct the spatio-temporal eigenvector of the Trapezoid functional
# at position x from the Floquet eigenvector ζ
function MonodromyQaD(::Val{:ExtractEigenVector}, poPb::PeriodicOrbitTrapProblem, u0::AbstractVector, par, ζ::AbstractVector)
	# extraction of various constants
	M = poPb.M
	N = poPb.N

	# period of the cycle
	T = extractPeriodFDTrap(poPb, u0)

	# time step
	h =  T * getTimeStep(poPb, 1)
	Typeh = typeof(h)

	out = copy(ζ)

	u0c = getTimeSlices(u0, N, M)

	@views out .= out .+ h/2 .* apply(poPb.J(u0c[:, M-1], par), out)
	# res = (I - h/2 * poPb.J(u0c[:, 1])) \ out
	@views res, _ = poPb.linsolver(poPb.J(u0c[:, 1], par), out; a₀ = convert(Typeh, 1), a₁ = -h/2)
	out .= res
	out_a = [copy(out)]
	# push!(out_a, copy(out))

	for ii in 2:M-1
		h =  T * getTimeStep(poPb, ii)
		@views out .= out .+ h/2 .* apply(poPb.J(u0c[:, ii-1], par), out)
		# res = (I - h/2 * poPb.J(u0c[:, ii])) \ out
		@views res, _ = poPb.linsolver(poPb.J(u0c[:, ii], par), out; a₀ = convert(Typeh, 1), a₁ = -h/2)
		out .= res
		push!(out_a, copy(out))
	end
	push!(out_a, copy(ζ))

	return out_a
end

# Compute the monodromy matrix at `u0` explicitely, not suitable for large systems
function MonodromyQaD(JacFW::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp})  where {Tpb <: PeriodicOrbitTrapProblem, Tjacpb, Torbitguess, Tp}

	poPb = JacFW.pb
	u0 = JacFW.x
	par = JacFW.par

	# extraction of various constants
	M, N = size(poPb)

	# period of the cycle
	T = extractPeriodFDTrap(poPb, u0)

	# time step
	h =  T * getTimeStep(poPb, 1)

	u0c = getTimeSlices(u0, N, M)

	@views mono = Array(I - h/2 * (poPb.J(u0c[:, 1], par))) \ Array(I + h/2 * poPb.J(u0c[:, M-1], par))
	temp = similar(mono)

	for ii in 2:M-1
		# for some reason, the next line is faster than doing (I - h/2 * (poPb.J(u0c[:, ii]))) \ ...
		# also I - h/2 .* J seems to hurt (a little) the performances
		h =  T * getTimeStep(poPb, ii)
		@views temp = Array(I - h/2 * (poPb.J(u0c[:, ii], par))) \ Array(I + h/2 * poPb.J(u0c[:, ii-1], par))
		mono .= temp * mono
	end
	return mono
end
####################################################################################################
# simplified version of
# Fairgrieve, Thomas F., and Allan D. Jepson. “O. K. Floquet Multipliers.” SIAM Journal on Numerical Analysis 28, no. 5 (October 1991): 1446–62. https://doi.org/10.1137/0728075.
struct FloquetLUColl{E <: AbstractEigenSolver, Tb} <: AbstractFloquetSolver
	eigsolver::E		# not really needed
	B::Tb
	"""
	## Arguments
	- `eigls` an eigensolver
	- `ntot` total number of unknowns (without countinig the period)
	- `n` space dimension
	"""
	function FloquetLUColl(eigls::AbstractEigenSolver, ntot::Int, n::Int)
		eigls2 = checkFloquetOptions(eigls)
		# build the mass matrix
		B = zeros(ntot, ntot)
		B[end-n+1:end, end-n+1:end] .= I(n)
		return new{typeof(eigls2), typeof(B)}(eigls2, B)
	end
	FloquetLUColl(eigls::FloquetLUColl) = eigls
end

# based on Fairgrieve, Thomas F., and Allan D. Jepson. “O. K. Floquet Multipliers.” SIAM Journal on Numerical Analysis 28, no. 5 (October 1991): 1446–62. https://doi.org/10.1137/0728075.
@views function (fl::FloquetLUColl)(JacColl::FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}, nev; kwargs...) where {Tpb <: PeriodicOrbitOCollProblem, Tjacpb <: AbstractMatrix, Torbitguess, Tp}
	prob = JacColl.pb
	_J = JacColl.jacpb
	n, m, Ntst = size(prob)
	J = _J[1:end-1, 1:end-1]
	# case of v(0)
	J[end-n+1:end, 1:n] .= I(n)
	# case of v(1)
	J[end-n+1:end, end-n+1:end] .= -I(n)
	# solve generalized eigenvalue problem
	values, vecs = eigen(J, fl.B)
	# remove infinite eigenvalues
	ind = isinf.(values)
	indvalid = ind .== false
	vals = values[indvalid]
	# this are the Floquet multipliers
	μ = @. Complex(1 / (1 + vals))
	return log.(μ), Complex.(vecs[indvalid, :]), true
end


text:
# This file is part of GenericSchur.jl, released under the MIT "Expat" license
# Portions derived from LAPACK, see below.

# Schur decomposition (QZ algorithm) for generalized eigen-problems

# Developer notes:
# We don't currently implement GeneralizedHessenberg types so anyone else
# using internal functions is doomed to confusion.

using LinearAlgebra: givensAlgorithm

# decomposition is A = Q S Z', B = Q  Tmat Z'

function schur!(A::StridedMatrix{T}, B::StridedMatrix{T}; kwargs...) where {T<:Complex}
    ggschur!(A, B; kwargs...)
end

# ggschur! is similar to LAPACK zgges:
# Q is Vsl, Z is Vsr
# Tmat overwrites B
# S overwrites A

function ggschur!(A::StridedMatrix{T}, B::StridedMatrix{T};
                  wantQ::Bool=true, wantZ::Bool=true,
                  scale::Bool=true,
                  kwargs...) where T <: Complex

    n = checksquare(A)
    nb = checksquare(B)
    if nb != n
        throw(ArgumentError("matrices must have the same sizes"))
    end

    if scale
        scaleA, cscale, anrm = _scale!(A)
        scaleB, cscaleb, bnrm = _scale!(B)
    else
        scaleA = false
        scaleB = false
    end
    # maybe balance here

    bqr = qr!(B)
    # apply bqr.Q to A
    lmul!(bqr.Q', A)
    if wantQ
        Q = Matrix(bqr.Q)
    else
        Q = Matrix{T}(undef, 0, 0)
    end
    if wantZ
        Z = Matrix{T}(I, n, n)
    else
        Z =  Matrix{T}(undef, 0, 0)
    end

    # materializing R may waste memory; can we rely on storage in modified B?
    A, B, Q, Z = _hessenberg!(A, bqr.R, Q, Z)
    α, β, S, Tmat, Q, Z = _gqz!(A, B, Q, Z, true; kwargs...)

    # TODO if balancing, unbalance Q, Z

    if scaleA
        safescale!(S, cscale, anrm)
        # checkme: is this always correct?
        α = diag(S)
    end
    if scaleB
        safescale!(Tmat, cscaleb, bnrm)
        # checkme: is this always correct?
        β = diag(Tmat)
    end

    return GeneralizedSchur(S, Tmat, α, β, Q, Z)
end

# B temporarily loses triangularity, so is not specially typed.
# compare to zgghrd
function _hessenberg!(A::StridedMatrix{T}, B::StridedMatrix{T}, Q, Z;
                      ilo=1, ihi=size(A,1)
                      ) where T
    n = checksquare(A)
    wantQ = !isempty(Q)
    wantZ = !isempty(Z)

    triu!(B)

    for jc = ilo:ihi-2
        for jr = ihi:-1:jc+2
            # rotate rows jr-1,jr to null A[jr,jc]
            Gq,r = givens(A,jr-1,jr,jc)
            lmul!(Gq, A)
            lmul!(Gq, B)
            if wantQ
                rmul!(Q, Gq')
            end
            # rotate cols jr,jr-1 to null B[jr,jr-1]
            Gz,r = givens(B',jr,jr-1,jr)
            rmul!(A,Gz')
            rmul!(B,Gz')
            if wantZ
                rmul!(Z, Gz')
            end
        end
    end

    return triu!(A,-1), triu!(B), Q, Z
end

# single-shift QZ algo
#
# translated from LAPACK's zhgeqz
# LAPACK is released under a BSD license, and is
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
#
# H is Hessenberg
# B is upper-triangular
# This is an internal routine so we don't check.
function _gqz!(H::StridedMatrix{T}, B::StridedMatrix{T}, Q, Z, wantSchur;
               debug = false,
               ilo=1, ihi=size(H,1), maxiter=100*(ihi-ilo+1)
               ) where {T <: Complex}
    n = checksquare(H)
    wantQ = !isempty(Q)
    wantZ = !isempty(Z)

    if ilo > 1 || ihi < n
        @warn "almost surely not working for nontrivial ilo,ihi"
    end
    α = Vector{T}(undef, n)
    β = Vector{T}(undef, n)

    RT = real(T)
    ulp = eps(RT)
    safmin = safemin(RT)
    anorm = norm(view(H,ilo:ihi,ilo:ihi),2) # Frobenius
    bnorm = norm(view(B,ilo:ihi,ilo:ihi),2) # Frobenius
    atol = max(safmin, ulp * anorm)
    btol = max(safmin, ulp * bnorm)
    ascale = one(RT) / max(safmin, anorm)
    bscale = one(RT) / max(safmin, bnorm)
    half = 1 / RT(2)

    # TODO: if ihi < n, deal with ihi+1:n

    ifirst = ilo
    ilast = ihi

    if wantSchur
        ifirstm, ilastm = 1,n
    else
        ifirstm, ilastm = ilo,ihi
    end

    shiftcount = 0
    eshift = zero(T)

    niter = 0
    while true
        niter += 1
        if niter > maxiter
            @warn "convergence failure; factorization incomplete at ilast=$ilast"
            break
        end

        # warning: _istart, ifirst, ilast are as in LAPACK, not GLA
        # note zlartg(f,g,c,s,r) -> (c,s,r) = givensAlgorithm(f,g)

        # Split H if possible
        # Use 2 tests:
        # 1: H[j,j-1] == 0 || j=ilo
        # 2: B[j,j] = 0

        null_Btail = false # branch 50 in zhgeqz
        deflate_H = false  # branch 60 in zhgeqz
        ns_Bblock_found = false # branch 70 in zhgeqz

        if ilast == ilo
            deflate_H = true
        elseif abs1(H[ilast,ilast-1]) < atol
            H[ilast, ilast-1] = 0
            deflate_H = true
        end
        # For now, we follow LAPACK quite closely.
        # br is a branch indicator corresponding to goto targets
        br = deflate_H ? :br60 : :br_none
        if !deflate_H
            # deflate_H is only tentatively false at this point
            if abs(B[ilast,ilast]) <= btol
                B[ilast,ilast] = 0
                null_Btail = true
                br = :br50
            else
                # general case
                for j=ilast-1:-1:ilo
                    # Test 1
                    split_test1 = false
                    if j==ilo
                        split_test1 = true
                    else
                        if abs1(H[j,j-1]) <= atol
                            H[j,j-1] = 0
                            split_test1 = true
                        end
                    end
                    # Test 2
                    split_test2 = abs(B[j,j]) < btol
                    if split_test2
                        B[j,j] = 0
                        # Test 1a: check for 2 consec. small subdiags in H
                        split_test1a = !split_test1 &&
                            (abs1(H[j,j-1]) * (ascale * abs1(H[j+1,j]))
                             <= abs1(H[j,j]) * (ascale * atol))
                        # if both tests pass, split 1x1 block off top
                        # if remaining leading diagonal elt vanishes, iterate
                        if split_test1 || split_test1a
                            @mydebug println("splitting 1x1 at top (j=$j)")
                            for jch=j:ilast-1
                                G,r = givens(H[jch,jch],H[jch+1,jch],jch,jch+1)
                                H[jch,jch] = r
                                H[jch+1,jch] = 0
                                lmul!(G,view(H,:,jch+1:ilastm))
                                lmul!(G,view(B,:,jch+1:ilastm))
                                if wantQ
                                    rmul!(Q, G')
                                end
                                if split_test1a
                                    H[jch,jch-1] *= G.c
                                end
                                split_test1a = false
                                if abs1(B[jch+1, jch+1]) > btol
                                    if jch+1 >= ilast
                                        deflate_H = true
                                        br = :br60
                                        break
                                    else
                                        ifirst = jch+1
                                        ns_Bblock_found = true
                                        br = :br70
                                        break
                                    end
                                end
                                B[jch+1, jch+1] = 0
                            end # jch loop
                            if !ns_Bblock_found && !deflate_H
                                null_Btail = true
                                br = :br50
                            end
                            break
                        else
                            # only test 2 passed: chase 0 to B[ilast,ilast]
                            # then process as above
                            @mydebug println("chasing diag 0 in B from $j to $ilast")
                            for jch=j:ilast-1
                                G,r = givens(B[jch,jch+1],B[jch+1,jch+1],
                                             jch,jch+1)
                                B[jch,jch+1] = r
                                B[jch+1,jch+1] = 0
                                if jch < ilastm-1
                                    lmul!(G,view(B,:,jch+2:ilastm))
                                end
                                lmul!(G,view(H,:,jch-1:ilastm))
                                if wantQ
                                    rmul!(Q,G')
                                end
                                G,r = givens(conj(H[jch+1,jch]),conj(H[jch+1,jch-1]),
                                             jch,jch-1)
                                H[jch+1,jch] = r'
                                H[jch+1,jch-1] = 0
                                rmul!(view(H,ifirstm:jch,:),G')
                                rmul!(view(B,ifirstm:jch+1,:),G')
                                if wantZ
                                    rmul!(Z,G')
                                end
                            end # jch loop
                            null_Btail = true
                            br = :br50
                        end # if split_test1 || split_test1a
                    elseif split_test1 # but not split_test2, i.e. tiny B[j,j]
                        ifirst = j
                        ns_Bblock_found = true
                        br = :br70
                        break
                    end # if B[j,j] tiny
                    # neither test passed; try next j
                end # j loop
                @mydebug println("after j loop branch $br")
                if br == :br_none
                    error("dev error: drop-through assumed impossible")
                end
            end # B[ilast,ilast] is/is-not tiny branches
            if null_Btail # br == :br50
                # B[ilast,ilast] is 0; clear H[ilast,ilast-1] to split
                G,r = givens(conj(H[ilast,ilast]),conj(H[ilast,ilast-1]),ilast,ilast-1)
                H[ilast,ilast] = r'
                H[ilast,ilast-1] = 0
                rmul!(view(H,ifirstm:ilast-1,:),G')
                rmul!(view(B,ifirstm:ilast-1,:),G')
                if wantZ
                    rmul!(Z,G')
                end
                deflate_H = true
                br = :br60
            end
        end # if trivial split
        if deflate_H
            # H[ilast, ilast-1] == 0: standardize B, set α,β
            absb = abs(B[ilast,ilast])
            @mydebug println("deflation ilast=$ilast, |B_tail|=$absb")
            if absb > safmin
                signbc = conj(B[ilast,ilast] / absb)
                B[ilast,ilast] = absb
                if wantSchur
                    B[ifirstm:ilast-1,ilast] .*= signbc
                    H[ifirstm:ilast,ilast] .*= signbc
                else
                    H[ilast,ilast] *= signbc
                end
                if wantZ
                    Z[:,ilast] .*= signbc
                end
            else
                B[ilast,ilast] = 0
            end
            α[ilast] = H[ilast,ilast]
            β[ilast] = B[ilast,ilast]
            ilast -= 1
            if ilast < ilo
                # we are done; exit outer loop
                break
            end
            shiftcount = 0
            eshift = zero(T)
            if !wantSchur
                ilastm = ilast
                if ifirstm > ilast
                    ifirstm = ilo
                end
            end
            continue # iteration loop
        end # deflation branch

        @assert ns_Bblock_found

        # QZ step
        # This iteration involves block ifirst:ilast, assumed nonempty
        if !wantSchur
            ifirstm = ifirst
        end

        # Compute shift
        shiftcount += 1
        # At this point ifirst < ilast and diagonal elts of B in that block
        # are all nontrivial.
        if shiftcount % 10 !=  0
            # Wilkinson shift, i.e. eigenvalue of last 2x2 block of
            # A B⁻¹ nearest to last elt

            # factor B as U D where U has unit diagonal, compute A D⁻¹ U⁻¹
            b11 = (bscale * B[ilast-1,ilast-1])
            b22 = (bscale * B[ilast,ilast])
            u12 = (bscale * B[ilast-1,ilast]) / b22
            ad11 = (ascale * H[ilast-1,ilast-1]) / b11
            ad21 = (ascale * H[ilast,ilast-1]) / b11
            ad12 = (ascale * H[ilast-1,ilast]) / b22
            ad22 = (ascale * H[ilast,ilast]) / b22
            abi22 = ad22 - u12 * ad21
            t1 = half * (ad11 + abi22)
            rtdisc = sqrt(t1^2 + ad12 * ad21 - ad11 * ad22)
            t2 = real(t1 - abi22) * real(rtdisc) + imag(t1-abi22) * imag(rtdisc)
            shift =  (t2 <= 0) ? (t1 + rtdisc) : (t1 - rtdisc)
            if isnan(shift)
                shift = zero(T)
                @warn "logic error in generalized Schur, check results"
                @mydebug println("NaN encountered bscale=$bscale " *
                                 "b11=$(B[ilast-1,ilast-1]) " *
                                 "b22=$(B[ilast,ilast])")
            end
        else
            # exceptional shift
            # "Chosen for no particularly good reason" (LAPACK)
            eshift += (ascale * H[ilast, ilast-1]) / (bscale * B[ilast-1, ilast-1])
            if isnan(eshift)
                @warn "logic error in generalized Schur, check results"
                @mydebug println("NaN encountered bscale=$bscale " *
                                 "b11=$(B[ilast-1,ilast-1]) ")
                eshift = zero(T)
            end
            shift = eshift
        end

        # check for two consecutive small subdiagonals
        local f1
        gotit = false
        for j=ilast-1:-1:ifirst+1
            _istart = j
            f1 = ascale * H[j,j] - shift * (bscale * B[j,j])
            t1 = abs1(f1)
            t2 = ascale * abs1(H[j+1,j])
            tr = max(t1,t2)
            if tr < one(RT) && tr != zero(RT)
                t1 /= tr
                t2 /= tr
            end
            if abs1(H[j,j-1]) * t2 <= t1 * atol
                gotit = true
                break
            end
        end
        if !gotit
            _istart = ifirst
            f1 = (ascale * H[ifirst, ifirst]
                  - shift * (bscale * B[ifirst, ifirst]))
        end

        # qz sweep
        @mydebug println("QZ sweep range $(_istart):$ilast shift=$shift")
        g2 = ascale*H[_istart+1, _istart]
        c,s,t3 = givensAlgorithm(f1, g2)
        for j=_istart:ilast-1
            if j > _istart
                    c,s,r = givensAlgorithm(H[j,j-1], H[j+1,j-1])
                    H[j,j-1] = r
                    H[j+1,j-1] = 0
            end
            G = Givens(j,j+1,T(c),s)
            lmul!(G, view(H,:,j:ilastm))
            lmul!(G, view(B,:,j:ilastm))
            if wantQ
                rmul!(Q, G')
            end
            c,s,r = givensAlgorithm(B[j+1,j+1], B[j+1,j])
            B[j+1,j+1] = r
            B[j+1,j] = 0
            G = Givens(j,j+1,T(c),s)
            rmul!(view(H,ifirstm:min(j+2,ilast),:),G)
            rmul!(view(B,ifirstm:j,:),G)
            if wantZ
                rmul!(Z, G)
            end
        end
    end # iteration loop
    @mydebug println("ggschur! done in $niter iters")
    # TODO if ilo > 1, deal with 1:ilo-1

    return α, β, H, B, Q, Z
end

function eigvecs(S::GeneralizedSchur, left=false)
    left && throw(ArgumentError("not implemented"))
    v = _geigvecs(S.S, S.T, S.Z)
    # CHECKME: Euclidean norm differs from LAPACK, so wait for upstream.
    # _enormalize!(v)
    return v
end

# right eigenvectors
# translated from LAPACK::ztgevc
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function _geigvecs(S::StridedMatrix{T}, P::StridedMatrix{T},
                    Z::StridedMatrix{T}=Matrix{T}(undef,0,0)
                    ) where {T <: Complex}
    n = size(S,1)
    RT = real(T)
    ulp = eps(RT)
    safmin = safemin(RT)
    smallnum = safmin * (n / ulp)
    bigx = one(RT) / smallnum
    bignum = one(RT) / (n * safmin)

    vectors = Matrix{T}(undef,n,n)
    v = zeros(T,n)
    if size(Z,1) > 0
        v2 = zeros(T,n)
    end

    # We use the 1-norms of the strictly upper part of S, P columns
    # to avoid overflow
    anorm = abs1(S[1,1])
    bnorm = abs1(P[1,1])
    snorms = zeros(RT,n)
    pnorms = zeros(RT,n)
    @inbounds for j=2:n
        for i=1:j-1
            snorms[j] += abs1(S[i,j])
            pnorms[j] += abs1(P[i,j])
        end
        anorm = max(anorm,snorms[j] + abs1(S[j,j]))
        bnorm = max(bnorm,pnorms[j] + abs1(P[j,j]))
    end
    ascale = one(RT) / max(anorm, safmin)
    bscale = one(RT) / max(bnorm, safmin)

    idx = n+1
    for ki=n:-1:1
        idx -= 1
        if abs1(S[ki,ki]) <= safmin && abs(real(P[ki,ki])) <= safmin
            # singular pencil; return unit eigenvector
            vectors[:,idx] .= zero(T)
            vectors[idx,idx] = one(T)
        else
            # compute coeffs a,b in (a A - b B) x = 0
            t = 1 / max(abs1(S[ki,ki]) * ascale,
                        abs(real(P[ki,ki])) * bscale, safmin)
            sα = (t * S[ki,ki]) * ascale
            sβ = (t * real(P[ki,ki])) * bscale
            acoeff = sβ * ascale
            bcoeff = sα * bscale
            # scale to avoid underflow

            lsa = abs(sβ) >= safmin && abs(acoeff) < smallnum
            lsb = abs1(sα) >= safmin && abs1(bcoeff) < smallnum
            s = one(RT)
            if lsa
                s = (smallnum / abs(sβ)) * min(anorm, bigx)
            end
            if lsb
                s = max(s, (smallnum / abs1(sα)) * min(bnorm, bigx))
            end
            if lsa || lsb
                s = min(s, 1 / (safmin * max( one(RT), abs(acoeff),
                                              abs1(bcoeff))))
                if lsa
                    acoeff = ascale * (s * sβ)
                else
                    acoeff *= s
                end
                if lsb
                    bcoeff = bscale * (s * sα)
                else
                    bcoeff *= s
                end
            end
            aca = abs(acoeff)
            acb = abs1(bcoeff)
            xmax = one(RT)
            v .= zero(T)
            dmin = max(ulp * aca * anorm, ulp * acb * bnorm, safmin)

            # triangular solve of (a A - b B) x = 0, columnwise
            # v[1:j-1] contains sums w
            # v[j+1:ki] contains x

            v[1:ki-1] .= acoeff * S[1:ki-1,ki] - bcoeff * P[1:ki-1,ki]
            v[ki] = one(T)

            for j=ki-1:-1:1
                # form x[j] = -v[j] / d
                # with scaling and perturbation
                d = acoeff * S[j,j] - bcoeff * P[j,j]
                if abs1(d) <= dmin
                    d = complex(dmin)
                end
                if abs1(d) < one(RT)
                    if abs1(v[j]) >= bignum * abs1(d)
                        t = 1 / abs1(v[j])
                        v[1:ki] *= t
                    end
                end
                v[j] = -v[j] / d
                if j > 1
                    # w = w + x[j] * (a S[:,j] - b P[:,j]) with scaling

                    if abs1(v[j]) > one(RT)
                        t = 1 / abs1(v[j])
                        if aca * snorms[j] + acb * pnorms[j] >= bignum * t
                            v[1:ki] *= t
                        end
                    end
                    ca = acoeff * v[j]
                    cb = bcoeff * v[j]
                    v[1:j-1] += ca * S[1:j-1,j] - cb * P[1:j-1,j]
                end
            end # for j (solve loop)
            if size(Z,1) > 0
                mul!(v2, Z, v)
                v, v2 = v2, v
                iend = n
            else
                iend = ki
            end

            xmax = zero(RT)
            for jr=1:iend
                xmax = max(xmax, abs1(v[jr]))
            end
            if xmax > safmin
                t = 1 / xmax
                vectors[1:iend,idx] .= t * v[1:iend]
            else
                iend = 0
            end
            vectors[iend+1:n,idx] .= zero(T)
        end # nonsingular branch
    end # index loop

    return vectors
end # function

"""
    canonicalize!(S::GeneralizedSchur)

rescale the fields of `S` so that `S.β` and the diagonal of `S.T` are real.
"""
function canonicalize!(S::GeneralizedSchur{Ty}) where Ty
    n = size(S.S,1)
    sf = safemin(real(Ty))
    for k=1:n
        scale = abs(S.T[k,k])
        if scale > sf
            t1 = conj(S.T[k,k] / scale)
            t2 = S.T[k,k] / scale
            S.T[k,k+1:n] .*= t1
            S.S[k,k:n] .*= t1
            if !(S.Q === nothing) && !isempty(S.Q)
                S.Q[:,k] .*= t2
            end
        else
            S.T[k,k] = zero(Ty)
        end
        S.α[k] = S.S[k,k]
        S.β[k] = S.T[k,k]
    end
    S
end


text:
using DynamicIterators
import DynamicIterators: dyniterate

struct Summed <: DynamicIterator
end

function dyniterate(::Summed, ::Nothing, y)
    y, y
end

function dyniterate(::Summed, i, y)
    i + y, i + y
end

@show collect(bind(1:5, Summed()))


text:
module Compiler

import ..Enzyme: Const, Active, Duplicated, DuplicatedNoNeed
import ..Enzyme: API, TypeTree, typetree, TypeAnalysis, FnTypeInfo

using LLVM, GPUCompiler, Libdl
import Enzyme_jll

import GPUCompiler: CompilerJob, FunctionSpec, codegen
using LLVM.Interop
import LLVM: Target, TargetMachine

# We have one global JIT and TM
const jit = Ref{OrcJIT}()
const tm  = Ref{TargetMachine}()

function __init__()
    opt_level = Base.JLOptions().opt_level
    if opt_level < 2
        optlevel = LLVM.API.LLVMCodeGenLevelNone
    elseif opt_level == 2
        optlevel = LLVM.API.LLVMCodeGenLevelDefault
    else
        optlevel = LLVM.API.LLVMCodeGenLevelAggressive
    end

    tm[] = GPUCompiler.JITTargetMachine(optlevel=optlevel)
    LLVM.asm_verbosity!(tm[], true)

    jit[] = OrcJIT(tm[]) # takes ownership of tm
    atexit() do
        dispose(jit[])
    end
end

# Define EnzymeTarget
Base.@kwdef struct EnzymeTarget <: AbstractCompilerTarget
end
GPUCompiler.llvm_triple(::EnzymeTarget) = Sys.MACHINE

# GPUCompiler.llvm_datalayout(::EnzymeTarget) =  nothing

function GPUCompiler.llvm_machine(::EnzymeTarget)
    return tm[]
end

module Runtime
    # the runtime library
    signal_exception() = return
    malloc(sz) = Base.Libc.malloc(sz)
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return
end

struct EnzymeCompilerParams <: AbstractCompilerParams end

## job

# TODO: We shouldn't blanket opt-out
GPUCompiler.check_invocation(job::CompilerJob{EnzymeTarget}, entry::LLVM.Function) = nothing

GPUCompiler.runtime_module(target::CompilerJob{EnzymeTarget}) = Runtime
GPUCompiler.isintrinsic(::CompilerJob{EnzymeTarget}, fn::String) = true
GPUCompiler.can_throw(::CompilerJob{EnzymeTarget}) = true

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
GPUCompiler.runtime_slug(job::CompilerJob{EnzymeTarget}) = "enzyme"

include("compiler/optimize.jl")
include("compiler/cassette.jl")

"""
Create the `FunctionSpec` pair, and lookup the primal return type.
"""
@inline function fspec(f::F, tt::TT) where {F, TT}
    # Entry for the cache look-up
    adjoint = FunctionSpec(f, tt, #=kernel=# false, #=name=# nothing)

    # primal function. Inferred here to get return type
    _tt = (tt.parameters...,)
    overdub_tt = Tuple{typeof(Compiler.CTX), F, map(eltype, _tt)...}
    primal = FunctionSpec(Cassette.overdub, overdub_tt, #=kernel=# false, #=name=# nothing)

    # can't return array since that's complicated.
    rt = Core.Compiler.return_type(Cassette.overdub, overdub_tt)
    if !(rt<:Union{AbstractFloat, Nothing})
        @error "Return type should be <:Union{Nothing, AbstractFloat}" rt adjoint primal
        error("Internal Enzyme Error")
    end
    return primal, adjoint, rt
end

function annotate!(mod)
    inactive = LLVM.StringAttribute("enzyme_inactive", "", context(mod))
    for inactivefn in ["jl_gc_queue_root"]
        fn = functions(mod)[inactivefn]
        push!(function_attributes(fn), inactive)
    end
end

function enzyme!(mod, primalf, adjoint, rt, split)
    ctx     = context(mod)
    rettype = convert(LLVMType, rt, ctx)
    dl      = string(LLVM.datalayout(mod))

    tt = [adjoint.tt.parameters...,]

    args_activity     = API.CDIFFE_TYPE[]
    uncacheable_args  = Bool[]
    args_typeInfo     = TypeTree[]
    args_known_values = API.IntList[]

    for T in tt
        if T <: Const
            push!(args_activity, API.DFT_CONSTANT)
        elseif T <: Active
            push!(args_activity, API.DFT_OUT_DIFF)
        elseif  T <: Duplicated
            push!(args_activity, API.DFT_DUP_ARG)
        elseif T <: DuplicatedNoNeed
            push!(args_activity, API.DFT_DUP_NONEED)
        else
            @assert("illegal annotation type")
        end
        typeTree = typetree(T, ctx, dl)
        push!(args_typeInfo, typeTree)
        if split
            push!(uncacheable_args, true)
        else
            push!(uncacheable_args, false)
        end
        push!(args_known_values, API.IntList())
    end

    # TODO ABI returned
    # The return of createprimal and gradient has this ABI
    #  It returns a struct containing the following values
    #     If requested, the original return value of the function
    #     If requested, the shadow return value of the function
    #     For each active (non duplicated) argument
    #       The adjoint of that argument

    if rt <: Integer
        retType = API.DFT_CONSTANT
    elseif rt <: AbstractFloat
        retType = API.DFT_OUT_DIFF
    elseif rt == Nothing
        retType = API.DFT_CONSTANT
    else
        error("What even is $rt")
    end

    TA = TypeAnalysis(triple(mod))
    global_AA = API.EnzymeGetGlobalAA(mod)
    retTT = typetree(rt, ctx, dl)

    typeInfo = FnTypeInfo(retTT, args_typeInfo, args_known_values)

    if split
        augmented = API.EnzymeCreateAugmentedPrimal(
            primalf, retType, args_activity, TA, global_AA, #=returnUsed=# true,
            typeInfo, uncacheable_args, #=forceAnonymousTape=# false, #=atomicAdd=# false, #=postOpt=# false)

        # 2. get new_primalf
        augmented_primalf = LLVM.Function(API.EnzymeExtractFunctionFromAugmentation(augmented))

        # TODOs:
        # 1. Handle mutable or !pointerfree arguments by introducing caching
        #     + specifically by setting uncacheable_args[i] = true
        # 2. Forward tape from augmented primalf to adjoint (as last arg)
        # 3. Make creation of augumented primalf vs joint forward and reverse optional

        tape = API.EnzymeExtractTapeTypeFromAugmentation(augmented)
        data = Array{Int64, 1}(undef, 3)
        existed = Array{UInt8, 1}(undef, 3)

        API.EnzymeExtractReturnInfo(augmented, data, existed)

        adjointf = LLVM.Function(API.EnzymeCreatePrimalAndGradient(
            primalf, retType, args_activity, TA, global_AA,
            #=returnValue=#false, #=dretUsed=#false, #=topLevel=#false,
            #=additionalArg=#tape, typeInfo,
            uncacheable_args, augmented, #=atomicAdd=#false, #=postOpt=#false))
    else
        adjointf = LLVM.Function(API.EnzymeCreatePrimalAndGradient(
            primalf, retType, args_activity, TA, global_AA,
            #=returnValue=#false, #=dretUsed=#false, #=topLevel=#true,
            #=additionalArg=#C_NULL, typeInfo,
            uncacheable_args, #=augmented=#C_NULL, #=atomicAdd=#false, #=postOpt=#false))
        augmented_primalf = nothing
    end

    API.EnzymeFreeGlobalAA(global_AA)
    return adjointf, augmented_primalf
end

include("compiler/thunk.jl")
include("compiler/reflection.jl")
# include("compiler/validation.jl")

end


text:
include("functions.jl")
include("steady.jl")
include("instantaneous.jl")


text:
using MatrixCompletion.Utilities.BatchUtils

@testset "$(format("BatchUtils: Tables Header"))" begin
  let
    tc = collect(1:100)
    batch = BatchFactory{SequentialScan}(size = 10)
    initialize(batch, tc)

    expected = [1:10, 11:20, 21:30, 31:40, 41:50, 51:60, 61:70, 71:80, 81:90, 91:100]
    ptr = 1
    while has_next(batch)
      @test tc[next(batch)] == collect(expected[ptr])
      ptr += 1
    end
  end

  let
    tc = collect(1:100)
    batch = BatchFactory{SequentialScan}(size = 64)
    initialize(batch, tc)
    expected = [1:64, 65:100]
    ptr = 1
    while has_next(batch)
      @test tc[next(batch)] == collect(expected[ptr])
      ptr += 1
    end
  end
end


text:
"""
    BCDiag

A structure of diagnostic properties of a Binary Classifier,
facilitates summary plots and tables.
"""
struct BCDiag
    n::Int                  # Number of observations
    n1::Int                 # Number of Class 1
    n0::Int                 # Number of Class 0
    baserate::Float64       # Class 1 incidence rate
    ks::Float64             # Maximum separation
    ksarg::Float64          # max sep occurs at this predictor value
    ksdep::Float64          # max sep occurs at this depth

    concordant::Int         # concordant pairs
    tied::Int               # tied pairs
    discordant::Int         # discordant pairs
    auc::Float64            # Area under ROC = (concordant + 0.5tied)/ total pairs
    gini::Float64           # Gini = 2C - 1

    grp::Vector             # group number [0-groups]
    depth::Vector           # Depth [0, 1]
    cdf1::Vector            # cum1 / n1 = True Positive Rate
    cdf0::Vector            # cum0 / n0 = False Positive Rate

    cntg::Vector            # count this group
    cnt1::Vector            # count class 1 this group
    prd1::Vector            # predicted class 1 this group
    rrObs::Vector           # response rate actual
    rrPrd::Vector           # response rate Predicted

    cumg::Vector            # cum count up to this group
    cum1::Vector            # cum count class 1 up to this group
    cpr1::Vector            # cum predicted class 1 up to this group
    crObs::Vector           # cum response rate actual
    crPrd::Vector           # cum response rate predicted
end

function Base.show(io::IO, ::MIME"text/plain", x::BCDiag)
    print(io,
        "Base rate: ", round(x.baserate, digits=4),
        "   n: ", x.n,
        "   n1: ", x.n1,
        "   n0: ", x.n0,
        "\nks:        ", round(x.ks, digits=4),
        "   occurs at value of ", x.ksarg,
        " depth of ", x.ksdep,
        "\nroc:       ", round(x.auc, digits=4),
        "   concordant pairs: ", x.concordant,
        "   tied pairs: ", x.tied,
        "   discordant pairs: ", x.discordant,
        "\nGini:      ", round(x.gini, digits=4),
    )
end

function Base.show(io::IO, x::BCDiag)
    print(io, round(x.ks, digits=4), "   ", round(x.auc, digits=4))
end

"""
    bcdiag(target, pred; groups = 100, rev = true, tie = 1e-6)

Perform diagnostics of a binary classifier.\\
`target` is a 2 level categorical variable, `pred` is probability of class 1.\\
`groups` is the number of bins to use for plotting/printing.\\
`rev` = true orders `pred` from high to low.\\
`tie` is the tolerance of `pred` where values are considered tied.\\

Returns a BCDiag struct which can be used for plotting or printing:
- `biasplot` is calibration plot of `target` response rate vs. `pred` response rate
- `ksplot` produces ksplot of cumulative distributions
- `rocplot` plots the Receiver Operating Characteristics curve
- `accuracyplot` plots the accuracy curve with adjustable utility
- `liftcurve` is the lift curve
- `cumliftsurve` is the cumulative lift surve
- `liftable` is the lift table as a DataFrame
- `cumliftable` is the cumulative lift table as a DataFrame
"""
function bcdiag(target::BitVector, pred::Vector; groups=100, rev=true, tie=1e-6)
    ks = kstest(target, pred; rev=rev)
    auroc = concordance(target, pred, tie)

    grpid = ranks(pred, groups=groups, rev=rev)
    frqtbl = freqtable(grpid, target)
    frq = frqtbl.array
    grp = NamedArrays.names(frqtbl, 1)

    n0, n1 = sum(frq, dims=1)
    n = n0 + n1

    cntg = vec(sum(frq, dims=2))            # bin count
    cnt1 = @view frq[:, 2]                  # bin count class 1
    cnt0 = @view frq[:, 1]                  # bin count class 0
    cumg = cumsum(cntg)                     # cum bin count
    cum1 = cumsum(cnt1)                     # cum class 1
    cum0 = cumsum(cnt0)                     # cum class 0

    depth = cumg ./ n                       # Depth including this group
    cdf1 = cum1 ./ n1                       # TPR: True Positive Rate
    cdf0 = cum0 ./ n0                       # FPR: False Positive Rate

    rrObs = cnt1 ./ cntg                    # actual response rate
    crObs = cum1 ./ cumg                    # cum actual response rate

    prd1 = Vector{Float64}(undef, length(grp))  # predicted class 1
    for (i, g) in enumerate(grp)
        prd1[i] = sum(pred[grpid .== g])
    end
    cpr1 = cumsum(prd1)                     # cum predicted class 1

    rrPrd = prd1 ./ cntg                    # predicted response rate
    crPrd = cpr1 ./ cumg                    # cum predicted response rate

    BCDiag(
        ks...,
        auroc...,
        grp, depth, cdf1, cdf0,
        cntg, cnt1, prd1, rrObs, rrPrd,
        cumg, cum1, cpr1, crObs, crPrd,
    )
end

function bcdiag(target::Vector, pred::Vector; groups=100, rev=true, tie=1e-6)
    uc = sort!(unique(target))
    length(uc) == 2 || error(ArgumentError("target must be 2 levels"))

    bcdiag(target .== uc[2], pred; rev=rev)
end

"""
    biasplot(x::BCDiag)

return bias calibration plot of `x` - actual response vs. predicted response
"""
function biasplot(x::BCDiag)
    dmin = min(minimum(x.rrPrd), minimum(x.rrObs))
    dmax = max(maximum(x.rrPrd), maximum(x.rrObs))

    plt = plot(
        size=(500, 500),
        aspect_ratio=:equal,
        legend=:bottomright,
        xguidefontsize=10,
        yguidefontsize=10,
        titlefontsize=11,
    )
    title!("Model Calibration Graph", titlefontsize=11)
    xlabel!("Predicted Response Rate", xguidefontsize=10)
    ylabel!("Observed Response Rate", yguidefontsize=10)

    plot!(x.rrPrd, x.rrObs, label=nothing, seriescolor=:blue)
    plot!([dmin, dmax], [dmin, dmax], label=nothing, width=0.5, seriescolor=:black)

    plt
end

"""
    ksplot(x::BCDiag)

return KS plot of `x` - CDF1 (True Positive) and CDF0 (False Positive) versus depth
"""
function ksplot(x::BCDiag)
    plt = plot(
        size=(500, 500),
        aspect_ratio=:equal,
        xlims=(0.0, 1.0),
        ylims=(0.0, 1.0),
        legend=:bottomright,
        xguidefontsize=10,
        yguidefontsize=10,
        titlefontsize=11,
    )
    title!("KS Graph", titlefontsize=11)
    xlabel!("Depth", xguidefontsize=10)
    ylabel!("Empirical CDF", yguidefontsize=10)

    plot!(
        [0.0; x.depth],
        [0.0; x.cdf1],
        label="Class 1",
        width=1.5,
        seriescolor=:blue,
    )

    plot!(
        [0.0; x.depth],
        [0.0; x.cdf0],
        label="Class 0",
        width=0.5,
        seriescolor=:black,
    )

    plot!(
        [0.0, x.baserate, 1.0],
        [0.0, 1.0, 1.0],
        label="Perfect Class 1",
        width=0.5,
        seriescolor=:blue,
        linestyle=:dash,
    )

    plot!(
        [0.0, x.baserate, 1.0],
        [0.0, 0.0, 1.0],
        label="Perfect Class 0",
        width=0.5,
        seriescolor=:black,
        linestyle=:dash,
    )
    ksidx = findlast(x.depth .<= x.ksdep)

    plot!(
        [x.depth[ksidx], x.depth[ksidx]],
        [x.cdf0[ksidx], x.cdf1[ksidx]],
        label=nothing,
        width=0.8,
        linestyle=:dash,
    )

    annotate!((
        0.03,
        0.97,
        text("Base rate = $(round(x.baserate, digits=4))", :black, :left, 9),
    ))
    annotate!((0.03, 0.93, text("Depth = $(round(x.ksdep, digits=4))", :black, :left, 9)))
    annotate!((0.03, 0.89, text("KSarg = $(round(x.ksarg, digits=4))", :black, :left, 9)))
    annotate!((0.03, 0.85, text("KS = $(round(x.ks, digits=4))", :black, :left, 9)))

    plt
end

"""
    rocplot(x::BCDiag)

return ROC plot of `x` - CDF1 (True Positive) vs. CDF0 (False Positive)
"""
function rocplot(x::BCDiag)
    plt = plot(
        size=(500, 500),
        aspect_ratio=:equal,
        xlims=(0.0, 1.0),
        ylims=(0.0, 1.0),
        legend=:bottomright,
        xguidefontsize=10,
        yguidefontsize=10,
        titlefontsize=11,
    )
    title!("ROC Graph", titlefontsize=11)
    xlabel!("False Positive Rate", xguidefontsize=10)
    ylabel!("True Positive Rate", yguidefontsize=10)

    plot!([0.0; x.cdf0], [0.0; x.cdf1], label="Model", seriescolor=:blue, width=1.5)
    plot!(
        [0.0, 1.0],
        [0.0, 1.0],
        label="Random",
        seriescolor=:black,
        width=0.3,
        linestyle=:dash,
    )

    annotate!((0.75, 0.22, text("auc = $(round(x.auc, digits=4))", :black, :left, 9)))
    annotate!((
        0.75,
        0.18,
        text("Gini = $(round(2. * x.auc - 1., digits=4))", :black, :left, 9),
    ))

    plt
end

"""
    accuracyplot(x::BCDiag; util=[1, 0, 0, 1])

Using `util` values for [TP, FN, FP, TN], produce accuracy plot and its [max, argmax, argdep].\\
Default `util` values of [1, 0, 0, 1] gives the standard accuracy value of (TP+TN)/N.
"""
function accuracyplot(x::BCDiag; util=[1, 0, 0, 1])
    p = x.baserate
    q = 1.0 - p

    # perfect classifier
    up = dot([p, 0.0, 0.0, q], util)                  # TP, FN, FP, TN of perfect classifier

    # 0.5 / 0.5 random classifier
    # ur5 = dot([0.5p, 0.5p, 0.5q, 0.5q], util)

    # baserate random classifier
    urb = dot([p * p, p * q, p * q, q * q], util)

    # model utility value [TP, FN, FP, TN]
    uc = [p .* x.cdf1 p .* (1.0 .- x.cdf1) q .* x.cdf0 q .* (1.0 .- x.cdf0)] * util
    ucmax = findmax(uc)

    plt = plot(
        size=(500, 500),
        xlims=(0.0, 1.0),
        legend=:best,
        xguidefontsize=10,
        yguidefontsize=10,
        titlefontsize=11,
    )
    title!("Accuracy Plot", titlefontsize=11)
    xlabel!("Depth", xguidefontsize=10)
    ylabel!("Utility", yguidefontsize=10)

    hline!([up], label="Perfect $(round(up, sigdigits=3))", linestyle=:dash)

    plot!(
        x.depth,
        uc,
        width=1.5,
        seriescolor=:blue,
        label="Model $(round(ucmax[1], sigdigits=3))",
    )
    vline!([x.depth[ucmax[2]]], label="Depth $(x.depth[ucmax[2]])")

    hline!([urb], label="Base Random $(round(urb, sigdigits=3))", linestyle=:dash)

    plt
end

"""
    liftcurve(x::BCDiag)

return lift curve plot of `x` - actual and predicted versus depth
"""
function liftcurve(x::BCDiag)
    plt = Plots.plot(
        size=(500, 500),
        xlims=(0.0, 1.0),
        ylims=(0.0, Inf), widen=true,
        legend=:topright,
        xguidefontsize=10,
        yguidefontsize=10,
        titlefontsize=11,
    )
    title!("Lift Curve", titlefontsize=11)
    xlabel!("Depth", xguidefontsize=10)
    ylabel!("Response Rate", yguidefontsize=10)

    plot!(x.depth, x.rrObs, label="Actual", width=1.0, seriescolor=:red)
    plot!(x.depth, x.rrPrd, label="Predicted", width=1.5, seriescolor=:blue)
    hline!([x.baserate], label="Base rate", linestyle=:dash)

    plt
end

"""
    cumliftcurve(x::BCDiag)

return cumulative lift curve plot of `x` - cumulative actual and predicted vs. depth
"""
function cumliftcurve(x::BCDiag)
    plt = Plots.plot(
        size=(500, 500),
        xlims=(0.0, 1.0),
        ylims=(0.0, Inf), widen=true,
        legend=:topright,
        xguidefontsize=10,
        yguidefontsize=10,
        titlefontsize=11,
        )
        title!("Cumulative Lift Curve", titlefontsize=11)
        xlabel!("Depth", xguidefontsize=10)
        ylabel!("Response Rate", yguidefontsize=10)

        plot!(x.depth, x.crObs, label="Actual", width=1.0, seriescolor=:red)
        plot!(x.depth, x.crPrd, label="Predicted", width=1.5, seriescolor=:blue)
        hline!([x.baserate], label="Base rate", linestyle=:dash)

    plt
end

"""
    liftable(x::BCDiag)

return lift table of `x` as a DataFrame
"""
function liftable(x::BCDiag)
    out = DataFrame(
        grp=x.grp,
        depth=x.depth,
        count=x.cntg,
        cntObs=x.cnt1,
        cntPrd=x.prd1,
        rrObs=x.rrObs,
        rrPred=x.rrPrd,
        liftObs=x.rrObs ./ x.baserate,
        liftPrd=x.rrPrd ./ x.baserate,
    )

    out
end

"""
    cumliftable(x::BCDiag)

return cumulative lift table of `x` as a DataFrame
"""
function cumliftable(x::BCDiag)
    out = DataFrame(
        grp=x.grp,
        depth=x.depth,
        count=x.cumg,
        cumObs=x.cum1,
        cumPrd=x.cpr1,
        crObs=x.cum1 ./ x.cumg,
        crPrd=x.cpr1 ./ x.cumg,
    )
    out.liftObs = out.crObs ./ x.baserate
    out.liftPrd = out.crPrd ./ x.baserate

    out
end


text:
export SequentialEnv

"""
    SequentialEnv(env)

Turn a simultaneous `env` into a sequential env.
"""
mutable struct SequentialEnv{E<:AbstractEnv} <: AbstractEnvWrapper
    env::E
    current_player_idx::Int
    actions::Vector{Any}
    function SequentialEnv(env::T) where T<:AbstractEnv
        @assert DynamicStyle(env) === SIMULTANEOUS "The SequentialEnv wrapper can only be applied to SIMULTANEOUS environments"
        new{T}(env, 1, Vector{Any}(undef, length(players(env))))
    end
end

RLBase.DynamicStyle(env::SequentialEnv) = SEQUENTIAL

RLBase.current_player(env::SequentialEnv) = env.current_player_idx

RLBase.legal_action_space(env::SequentialEnv) = legal_action_space(env, current_player(env))
RLBase.legal_action_space(env::SequentialEnv, player) = legal_action_space(env.env, player)

RLBase.action_space(env::SequentialEnv) = action_space(env, current_player(env))
RLBase.action_space(env::SequentialEnv, player) = action_space(env.env, player)

function RLBase.reset!(env::SequentialEnv)
    env.current_player_idx = 1
    reset!(env.env)
end

RLBase.reward(env::SequentialEnv) = reward(env, current_player(env))

RLBase.reward(env::SequentialEnv, player) = current_player(env) == 1 ? reward(env.env, player) : 0

function (env::SequentialEnv)(action)
    env.actions[env.current_player_idx] = action
    if env.current_player_idx == length(env.actions)
        env.env(env.actions)
        env.current_player_idx = 1
    else
        env.current_player_idx += 1
    end
end


text:
"""
Vibration analysis of free-floating thin plate

Free-vibration problem is solved for a homogeneous free-floating
(unsupported) square plate. This is the NAFEMS Benchmark, Test No. FV12.

The plate is discretized with shell elements. Because no displacements are
prevented, the structure has six rigid body modes (six zero vibration
frequencies).

The nonzero benchmark frequencies are (in hertz):

1.622, 2.360, 2.922, 4.190, 4.190,  7.356, 7.356, 7.668.
"""
module FV12_vibration_examples

using LinearAlgebra
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 = 8, visualize = true)
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.05*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
    l1 = collect(1:count(fens))
    # for i in [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);
    # @show sum(sum(M, dims = 1))/3,.

    # Solve
    OmegaShift = (0.5*2*pi)^2
    neigvs = 24
    evals, evecs, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    # @show nconv
    evals[:] = evals .- OmegaShift;
    fs = real(sqrt.(complex(evals)))/(2*pi)
    @info "Frequencies: $(round.(fs[7:10], digits=4))"

    # sol = eigen(Matrix(K+OmegaShift*M), Matrix(M))
    # evals = sol.values
    # evals[:] = evals .- OmegaShift;
    # fs = real(sqrt.(complex(evals)))/(2*pi)
    #     @show fs[9]
    #     evecs = sol.vectors

    # for index_ in 1:neigvs
    #     om = evals[index_]
    #     v = evecs[:, index_]
    #     @show index_, norm(K*v-om*M*v)
    #     @show (v'*K*v)/om
    #     @show (v'*M*v)
    # end

    # for n in dchi.dofnums[:, 6]
    #     # @show sqrt(K[n, n] / M[n, n]) / 2/pi
    #     if n != 0
    #         @show K[n, :]
    #     end
    # end

    # Visualization
    if !visualize
        return true
    end
    for ev in 1:10
        U = evecs[:, ev]
        scattersysvec!(dchi, (0.5*L)/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; title="$(ev)")
    end
    return true
end

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

    formul = FEMMShellQ4SRIModule
       # Report
    @info "FV12 free vibration, t/L=$(thickness/L), formulation=$(formul)"
    @info "Mesh: $n elements per side"

    tolerance = L/n/1000
    fens, fes = Q4block(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 = FESetShellQ4()
    accepttodelegate(fes, sfes)
    femm = formul.FEMMShellQ4SRI(IntegDomain(fes, GaussRule(2, 2), thickness), mater)
    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
    l1 = collect(1:count(fens))
    # for i in [6]
    #     setebc!(dchi, l1, true, i)
    # end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, dchi);
    # @show sum(sum(M, dims = 1))/3

    # Solve
    OmegaShift = (0.5*2*pi)^2
    neigvs = 24
    evals, evecs, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    # @show nconv
    evals[:] = evals .- OmegaShift;
    fs = real(sqrt.(complex(evals)))/(2*pi)
    @info "Frequencies: $(round.(fs[7:10], digits=4))"

    # sol = eigen(Matrix(K+OmegaShift*M), Matrix(M))
    # evals = sol.values
    # evals[:] = evals .- OmegaShift;
    # fs = real(sqrt.(complex(evals)))/(2*pi)
    #     @show fs[9]
    #     evecs = sol.vectors

    # for index_ in 1:neigvs
    #     om = evals[index_]
    #     v = evecs[:, index_]
    #     @show index_, norm(K*v-om*M*v)
    #     @show (v'*K*v)/om
    #     @show (v'*M*v)
    # end

    # for n in dchi.dofnums[:, 6]
    #     # @show sqrt(K[n, n] / M[n, n]) / 2/pi
    #     if n != 0
    #         @show K[n, :]
    #     end
    # end

    # Visualization
    if !visualize
        return true
    end
    for ev in 1:10
        U = evecs[:, ev]
        scattersysvec!(dchi, (0.5*L)/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; title="$(ev)")
    end
    return true
end

function single_q4sri()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.05*phun("m");
    L = 10.0*phun("m");

    lambdaij = sqrt.([13.49, 19.79, 24.43])
    fij(k) = lambdaij[k]^2/(2*pi*L^2)*sqrt(E*thickness^2/12/rho/(1-nu^2))

    @show fij.(1:3)
    # Mesh
    n = 8
    tolerance = L/n/1000
    fens, fes = Q4block(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 = FESetShellQ4()
    accepttodelegate(fes, sfes)
    femm = FEMMShellQ4SRIModule.FEMMShellQ4SRI(IntegDomain(fes, GaussRule(2, 2), thickness), mater)
    stiffness = FEMMShellQ4SRIModule.stiffness
    mass = FEMMShellQ4SRIModule.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
    l1 = collect(1:count(fens))
    # for i in [6]
    #     setebc!(dchi, l1, true, i)
    # end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    K = stiffness(femm, geom0, u0, Rfield0, dchi);
    M = mass(femm, geom0, dchi);
    # @show sum(sum(M, dims = 1))/3

    # Solve
    OmegaShift = (0.5*2*pi)^2
    neigvs = 10
    evals, evecs, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    # @show nconv
    evals[:] = evals .- OmegaShift;
    fs = real(sqrt.(complex(evals)))/(2*pi)
    @show fs

    # sol = eigen(Matrix(K+OmegaShift*M), Matrix(M))
    # evals = sol.values
    # evals[:] = evals .- OmegaShift;
    # fs = real(sqrt.(complex(evals)))/(2*pi)
    #     @show fs[9]
    #     evecs = sol.vectors

    # for index_ in 1:neigvs
    #     om = evals[index_]
    #     v = evecs[:, index_]
    #     @show index_, norm(K*v-om*M*v)
    #     @show (v'*K*v)/om
    #     @show (v'*M*v)
    # end

    # for n in dchi.dofnums[:, 6]
    #     # @show sqrt(K[n, n] / M[n, n]) / 2/pi
    #     if n != 0
    #         @show K[n, :]
    #     end
    # end

    # Visualization
    if !visualize
        return true
    end
    for ev in 7:neigvs
        U = evecs[:, ev]
        scattersysvec!(dchi, (0.5*L)/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; title="$(ev)")
    end
    return true
end

function test_convergence()
    formul = FEMMShellT3FFModule
    @info "FV12 free vibration, formulation=$(formul)"
    @info "1.622, 2.360, 2.922, 4.190, 4.190,  7.356, 7.356, 7.668."
    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:
__precompile__()
module LShapedSolvers

using TraitDispatch
using Parameters
using JuMP
using StochasticPrograms
using MathProgBase
using RecipesBase
using ProgressMeter

import Base: show, put!, wait, isready, take!, fetch
import StochasticPrograms: StructuredModel, optimsolver, optimize_structured!, fill_solution!
importall MathProgBase.SolverInterface

using Base.Order: Ordering, ReverseOrdering, Reverse
using DataStructures: PriorityQueue,enqueue!,dequeue!

export
    LShapedSolver,
    Crash,
    StructuredModel,
    optimsolver,
    optimize_structured!,
    fill_solution!,
    get_solution,
    get_objective_value

# Include files
include("LQSolver.jl")
include("subproblem.jl")
include("hyperplane.jl")
include("AbstractLShaped.jl")
include("regularization.jl")
include("parallel.jl")
include("solvers/LShaped.jl")
include("solvers/DLShaped.jl")
include("solvers/Regularized.jl")
include("solvers/DRegularized.jl")
include("solvers/TrustRegion.jl")
include("solvers/DTrustRegion.jl")
include("solvers/LevelSet.jl")
include("solvers/DLevelSet.jl")
include("crash.jl")
include("spinterface.jl")

end # module


text:
module Types
using Random

include("params.jl")
import .Params

abstract type Asset end

abstract type AssetManager end

abstract type Investor end

abstract type Order end

struct Equity <: Asset
    value::Array{Float64} # Value history of each stock
    beta::Array{Float64} # Stocks' betas
    vol::Array{Float64} # Stocks' volatilities
    impact::Array{Float64} # Symmetric price impact for 1 money unit bought/sold
    volume::Array{Float64} # Money value of the stock traded in a day
end

# NOTE: Mismatch between market value and what the index should be worth
# since the former determines the latter, while it should be vice versa.
# Implicit assumption is that there are other stocks not traded by our agents.
struct MarketIndex <: Asset
    value::Array{Float64} # Index, used to drive stock values
end

struct RetailInvestor <: Investor
    assets::Array{Float64} # N possible stakes in M funds and (1) cash
    horizon::Array{Int64} # Positive integer, draw from range, can use rand(S)
    threshold::Array{Float64} # ~N(0, 5)
end

struct EquityFund <: AssetManager
    holdings::Array{Float64} # Units of each stock held in the fund's portfolio
    stakes::Array{Float64} # Share of assets that each investor in the funds own
    value::Array{Float64} # Fund's value history
end

mutable struct BuyMarketOrder <: Order
    values::Array{Float64} # Money amount spent on each stock for one of .funds
    funds::Array{Int64} # Fund that a set of stocks (.values) gets bought for
end

mutable struct SellMarketOrder <: Order
    values::Array{Float64} # Money amount spent on each stock for one of .funds
    investors::Array{Int64} # Investor that stocks (.values) are liquidated for
end

end # module