text:
"""
"""
struct CellDof{DS} <: CellDatum
  cell_dof::AbstractArray{<:Union{Dof,AbstractArray{<:Dof}}}
  trian::Triangulation
  domain_style::DS
end

get_data(f::CellDof) = f.cell_dof
get_triangulation(f::CellDof) = f.trian
DomainStyle(::Type{CellDof{DS}}) where DS = DS()

function change_domain(a::CellDof,::ReferenceDomain,::PhysicalDomain)
  @notimplemented
end

function change_domain(a::CellDof,::PhysicalDomain,::ReferenceDomain)
  @notimplemented
end

function change_domain(a::CellDof,target_trian::Triangulation,target_domain::DomainStyle)
  @notimplemented
end

function change_domain(a::CellDof,target_trian::RestrictedTriangulation,target_domain::DomainStyle)
  @notimplementedif DomainStyle(a) != target_domain
  trian_a = get_triangulation(a)
  @notimplementedif ! have_compatible_domains(trian_a,get_background_triangulation(target_trian))
  cell_dof = get_data(a)
  tcell_to_cell = get_cell_to_bgcell(target_trian)
  tcell_dof = lazy_map(Reindex(cell_dof),tcell_to_cell)
  CellDof(tcell_dof,target_trian,DomainStyle(a))
end

# Evaluation of CellDof

(a::CellDof)(f) = evaluate(a,f)

function evaluate!(cache,s::CellDof,f::CellField)

  trian_f = get_triangulation(f)
  trian_s = get_triangulation(s)

  if trian_f !== trian_s
    @unreachable """\n
    A CellDof can only be evaluated on a CellField defined on the same Triangulation.
    """
  end

  b = change_domain(f,s.domain_style)
  lazy_map(evaluate,get_data(s),get_data(b))
end

function evaluate!(cache, ::CellField, ::CellDof)
  @unreachable """\n
  CellField (f) objects cannot be evaluated at CellDof (s) objects.
  However, CellDofs objects can be evaluated at CellField objects.
  Did you mean evaluate(f,s) instead of evaluate(s,f), i.e.
  f(s) instead of s(f)?
  """
end


text:
module Adaptation
export Adaptation

using LinearAlgebra: LinearAlgebra
using Statistics: Statistics
using Parameters: @unpack, @pack!

using ..AdvancedHMC: DEBUG, AbstractScalarOrVec

abstract type AbstractAdaptor end
function getM⁻¹ end
function getϵ end
function adapt! end
function reset! end
function initialize! end
function finalize! end
export AbstractAdaptor, adapt!, initialize!, finalize!, reset!, getϵ, getM⁻¹

struct NoAdaptation <: AbstractAdaptor end
export NoAdaptation
include("stepsize.jl")
export StepSizeAdaptor, NesterovDualAveraging
include("massmatrix.jl")
export MassMatrixAdaptor, UnitMassMatrix, WelfordVar, WelfordCov

##
## Composite adaptors
## TODO: generalise this to a list of adaptors
##

struct NaiveHMCAdaptor{M<:MassMatrixAdaptor, Tssa<:StepSizeAdaptor} <: AbstractAdaptor
    pc  :: M
    ssa :: Tssa
end
Base.show(io::IO, a::NaiveHMCAdaptor) = print(io, "NaiveHMCAdaptor(pc=$(a.pc), ssa=$(a.ssa))")

getM⁻¹(ca::NaiveHMCAdaptor) = getM⁻¹(ca.pc)
getϵ(ca::NaiveHMCAdaptor) = getϵ(ca.ssa)

# TODO: implement consensus adaptor
function adapt!(
    nca::NaiveHMCAdaptor,
    θ::AbstractVecOrMat{<:AbstractFloat},
    α::AbstractScalarOrVec{<:AbstractFloat}
)
    adapt!(nca.ssa, θ, α)
    adapt!(nca.pc, θ, α)
end
function reset!(aca::NaiveHMCAdaptor)
    reset!(aca.ssa)
    reset!(aca.pc)
end
initialize!(adaptor::NaiveHMCAdaptor, n_adapts::Int) = nothing
finalize!(aca::NaiveHMCAdaptor) = finalize!(aca.ssa)

include("stan_adaptor.jl")
export NaiveHMCAdaptor, StanHMCAdaptor

end # module

ext: jl


text:
# data is generated by test/gendata.jl

# load test data
countsint = Matrix{Int}(CSV.read(joinpath(testdir,"data","countsint.csv.gz"), DataFrame))
covarsdf = CSV.read(joinpath(testdir,"data","covarsdf.csv.gz"), DataFrame)
transform!(covarsdf, :cat => categorical, renamecols=false)
counts = convert(SparseMatrixCSC{Float64,Int}, countsint)

# construct equivalent covars matrix so we can show how that api works too
covars = ModelMatrix(ModelFrame(@formula(y ~ x + z + cat + y), covarsdf)).m[:,2:end]

n, d = size(counts)
p = size(covars, 2)
ncats = length(levels(covarsdf.cat))

# used for testing predict
newcovars = covars[1:10,:]

# used for testing srproj and fit(CIR)
projdir = size(covars,2)

γdistrom = 1.0

coefsRdistrom = sparse(Matrix{Float64}(CSV.read(joinpath(testdir,"data","dmr_coefsRdistrom.csv.gz"), DataFrame)))
zRdistrom = Matrix{Float64}(CSV.read(joinpath(testdir,"data","dmr_zRdistrom.csv.gz"), DataFrame))
z1Rdistrom = Matrix{Float64}(CSV.read(joinpath(testdir,"data","dmr_z1Rdistrom.csv.gz"), DataFrame))
predictRdistrom = Matrix{Float64}(CSV.read(joinpath(testdir,"data","dmr_predictRdistrom.csv.gz"), DataFrame))


text:
export iauFk425
"""
Convert B1950.0 FK4 star catalog data to J2000.0 FK5.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

This function converts a star's catalog data from the old FK4
(Bessel-Newcomb) system to the later IAU 1976 FK5 (Fricke) system.

Given: (all B1950.0, FK4)
   - r1950,d1950    double   B1950.0 RA,Dec (rad)
   - dr1950,dd1950  double   B1950.0 proper motions (rad/trop.yr)
   - p1950          double   parallax (arcsec)
   - v1950          double   radial velocity (km/s, +ve = moving away)

Returned: (all J2000.0, FK5)
   - r2000,d2000    double   J2000.0 RA,Dec (rad)
   - dr2000,dd2000  double   J2000.0 proper motions (rad/Jul.yr)
   - p2000          double   parallax (arcsec)
   - v2000          double   radial velocity (km/s, +ve = moving away)

Notes:

   1. The proper motions in RA are dRA/dt rather than cos(Dec)*dRA/dt,
   and are per year rather than per century.

   2. The conversion is somewhat complicated, for several reasons:

   . Change of standard epoch from B1950.0 to J2000.0.

   . An intermediate transition date of 1984 January 1.0 TT.

   . A change of precession model.

   . Change of time unit for proper motion (tropical to Julian).

   . FK4 positions include the E-terms of aberration, to simplify
      the hand computation of annual aberration.  FK5 positions
      assume a rigorous aberration computation based on the Earth's
      barycentric velocity.

   . The E-terms also affect proper motions, and in particular cause
      objects at large distances to exhibit fictitious proper
      motions.

   The algorithm is based on Smith et al. (1989) and Yallop et al.
   (1989), which presented a matrix method due to Standish (1982) as
   developed by Aoki et al. (1983), using Kinoshita's development of
   Andoyer's post-Newcomb precession.  The numerical constants from
   Seidelmann (1992) are used canonically.

   3. Conversion from B1950.0 FK4 to J2000.0 FK5 only is provided for.
   Conversions for different epochs and equinoxes would require
   additional treatment for precession, proper motion and E-terms.

   4. In the FK4 catalog the proper motions of stars within 10 degrees
   of the poles do not embody differential E-terms effects and
   should, strictly speaking, be handled in a different manner from
   stars outside these regions.  However, given the general lack of
   homogeneity of the star data available for routine astrometry,
   the difficulties of handling positions that may have been
   determined from astrometric fields spanning the polar and non-
   polar regions, the likelihood that the differential E-terms
   effect was not taken into account when allowing for proper motion
   in past astrometry, and the undesirability of a discontinuity in
   the algorithm, the decision has been made in this SOFA algorithm
   to include the effects of differential E-terms on the proper
   motions for all stars, whether polar or not.  At epoch J2000.0,
   and measuring "on the sky" rather than in terms of RA change, the
   errors resulting from this simplification are less than
   1 milliarcsecond in position and 1 milliarcsecond per century in
   proper motion.

Called:
   - iauAnp       normalize angle into range 0 to 2pi
   - iauPv2s      pv-vector to spherical coordinates
   - iauPdp       scalar product of two p-vectors
   - iauPvmpv     pv-vector minus pv_vector
   - iauPvppv     pv-vector plus pv_vector
   - iauS2pv      spherical coordinates to pv-vector
   - iauSxp       multiply p-vector by scalar

References:

   - Aoki, S. et al., 1983, "Conversion matrix of epoch B1950.0
   FK4-based positions of stars to epoch J2000.0 positions in
   accordance with the new IAU resolutions".  Astron.Astrophys.
   128, 263-267.

   - Seidelmann, P.K. (ed), 1992, "Explanatory Supplement to the
   Astronomical Almanac", ISBN 0-935702-68-7.

   - Smith, C.A. et al., 1989, "The transformation of astrometric
   catalog systems to the equinox J2000.0".  Astron.J. 97, 265.

   - Standish, E.M., 1982, "Conversion of positions and proper motions
   from B1950.0 to the IAU system at J2000.0".  Astron.Astrophys.,
   115, 1, 20-22.

   - Yallop, B.D. et al., 1989, "Transformation of mean star places
   from FK4 B1950.0 to FK5 J2000.0 using matrices in 6-space".
   Astron.J. 97, 274.

This revision:   2018 December 5

SOFA release 2019-07-22

Copyright (C) 2019 IAU SOFA Board.  See notes at end.
"""

# void iauFk425(double r1950, double d1950,
#               double dr1950, double dd1950,
#               double p1950, double v1950,
#               double *r2000, double *d2000,
#               double *dr2000, double *dd2000,
#               double *p2000, double *v2000)

function iauFk425(r1950::Real, d1950::Real, dr1950::Real, dd1950::Real,
                  p1950::Real, v1950::Real)
   ref_r2000 = Ref{Float64}(0.0)
   ref_d2000 = Ref{Float64}(0.0)
   ref_dr2000 = Ref{Float64}(0.0)
   ref_dd2000 = Ref{Float64}(0.0)
   ref_p2000 = Ref{Float64}(0.0)
   ref_v2000 = Ref{Float64}(0.0)

   ccall((:iauFk425, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, ),
         convert(Float64, r1950), convert(Float64, d1950),
         convert(Float64, dr1950), convert(Float64, dd1950),
         convert(Float64, p1950), convert(Float64, v1950),
         ref_r2000, ref_d2000,
         ref_dr2000, ref_dd2000,
         ref_p2000, ref_v2000)

   return ref_r2000[], ref_d2000[], ref_dr2000[], ref_dd2000[], ref_p2000[], ref_v2000[]
end


text:
module ConstraintDomains

# NOTE - Reexport.jl
using Reexport

@reexport using PatternFolds

# Imports
using Dictionaries

# Functions exports
export AbstractDomain, ContinuousDomain, DiscreteDomain, RangeDomain, DiscreteSet
export domain
export domain_size
export get_domain
export add!, delete!
export merge_domains, intersect_domains

# Includes
include("utils.jl")
include("common.jl")
include("continuous.jl")
include("discrete.jl")
include("general.jl")

end


text:
# JIT compilation of Julia code to PTX

export cufunction

struct CompilerContext
    # core invocation
    f::Core.Function
    tt::DataType
    cap::VersionNumber
    kernel::Bool

    # optional properties
    alias::Union{Nothing,String}
    minthreads::Union{Nothing,CuDim}
    maxthreads::Union{Nothing,CuDim}
    blocks_per_sm::Union{Nothing,Integer}
    maxregs::Union{Nothing,Integer}

    # hacks
    inner_f::Union{Nothing,Core.Function}

    CompilerContext(f, tt, cap, kernel; inner_f=nothing, alias=nothing,
                    minthreads=nothing, maxthreads=nothing, blocks_per_sm=nothing, maxregs=nothing) =
        new(f, tt, cap, kernel, alias, minthreads, maxthreads, blocks_per_sm, maxregs, inner_f)
end

struct CompilerError <: Exception
    ctx::CompilerContext
    message::String
    meta::Dict
end

compiler_error(ctx::CompilerContext, message="unknown error"; kwargs...) =
    throw(CompilerError(ctx, message, kwargs))

function Base.showerror(io::IO, err::CompilerError)
    ctx = err.ctx
    fn = typeof(something(ctx.inner_f, ctx.f)).name.mt.name
    args = join(ctx.tt.parameters, ", ")
    print(io, "could not compile $fn($args) for GPU; $(err.message)\n")
    if haskey(err.meta, :errors) && isa(err.meta[:errors], Vector{UnsupportedIRError})
        for suberr in err.meta[:errors]
            print(io, "\nReason: ")
            Base.showerror(io, suberr, suberr.compiletrace)
            println(io)
        end
        print(io, "\nTry inspecting generated code with the @device_code_... macros")
    else
        for (key,val) in err.meta
            print(io, "\n- $key = $val")
        end
    end
end

#
# main code generation functions
#

function module_setup(mod::LLVM.Module)
    # NOTE: NVPTX::TargetMachine's data layout doesn't match the NVPTX user guide,
    #       so we specify it ourselves
    if Int === Int64
        triple!(mod, "nvptx64-nvidia-cuda")
        datalayout!(mod, "e-p:64:64:64-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64-f32:32:32-f64:64:64-v16:16:16-v32:32:32-v64:64:64-v128:128:128-n16:32:64")
    else
        triple!(mod, "nvptx-nvidia-cuda")
        datalayout!(mod, "e-p:32:32:32-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64-f32:32:32-f64:64:64-v16:16:16-v32:32:32-v64:64:64-v128:128:128-n16:32:64")
    end

    # add debug info metadata
    push!(metadata(mod), "llvm.module.flags",
         MDNode([ConstantInt(Int32(1), JuliaContext()),    # llvm::Module::Error
                 MDString("Debug Info Version"),
                 ConstantInt(DEBUG_METADATA_VERSION(), JuliaContext())]))
end

# make function names safe for PTX
safe_fn(fn::String) = replace(fn, r"[^aA-zZ0-9_]"=>"_")
safe_fn(f::Core.Function) = safe_fn(String(typeof(f).name.mt.name))
safe_fn(f::LLVM.Function) = safe_fn(LLVM.name(f))

function raise_exception(insblock::BasicBlock, ex::Value)
    fun = LLVM.parent(insblock)
    mod = LLVM.parent(fun)
    ctx = context(mod)

    builder = Builder(ctx)
    position!(builder, insblock)

    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    call!(builder, trap)
end

# maintain our own "global unique" suffix for disambiguating kernels
globalUnique = 0

function irgen(ctx::CompilerContext)
    # get the method instance
    isa(ctx.f, Core.Builtin) && compiler_error(ctx, "function is not a generic function")
    world = typemax(UInt)
    meth = which(ctx.f, ctx.tt)
    sig_tt = Tuple{typeof(ctx.f), ctx.tt.parameters...}
    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig_tt, meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti)
    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)

    # set-up the compiler interface
    function hook_module_setup(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        module_setup(LLVM.Module(ref))
    end
    function hook_raise_exception(insblock::Ptr{Cvoid}, ex::Ptr{Cvoid})
        insblock = convert(LLVM.API.LLVMValueRef, insblock)
        ex = convert(LLVM.API.LLVMValueRef, ex)
        raise_exception(BasicBlock(insblock), Value(ex))
    end
    dependencies = Vector{LLVM.Module}()
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        push!(dependencies, LLVM.Module(ref))
    end
    params = Base.CodegenParams(cached=false,
                                track_allocations=false,
                                code_coverage=false,
                                static_alloc=false,
                                prefer_specsig=true,
                                module_setup=hook_module_setup,
                                module_activation=hook_module_activation,
                                raise_exception=hook_raise_exception)

    # get the code
    mod = let
        ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                    (Any, UInt, Bool, Bool, Base.CodegenParams),
                    linfo, world, #=wrapper=#false, #=optimize=#false, params)
        if ref == C_NULL
            compiler_error(ctx, "the Julia compiler could not generate LLVM IR")
        end

        llvmf = LLVM.Function(ref)
        LLVM.parent(llvmf)
    end

    # the main module should contain a single jfptr_ function definition,
    # e.g. jlcall_kernel_vadd_62977
    definitions = filter(f->!isdeclaration(f), functions(mod))
    wrapper = let
        fs = collect(filter(f->startswith(LLVM.name(f), "jfptr_"), definitions))
        @assert length(fs) == 1
        fs[1]
    end

    # the jlcall wrapper function should point us to the actual entry-point,
    # e.g. julia_kernel_vadd_62984
    entry_tag = let
        m = match(r"jfptr_(.+)_\d+", LLVM.name(wrapper))
        @assert m != nothing
        m.captures[1]
    end
    unsafe_delete!(mod, wrapper)
    entry = let
        re = Regex("julia_$(entry_tag)_\\d+")
        llvmcall_re = Regex("julia_$(entry_tag)_\\d+u\\d+")
        fs = collect(filter(f->occursin(re, LLVM.name(f)) &&
                               !occursin(llvmcall_re, LLVM.name(f)), definitions))
        if length(fs) != 1
            compiler_error(f, tt, cap, "could not find single entry-point";
                           entry=>entry_tag, available=>[LLVM.name.(definitions)])
        end
        fs[1]
    end

    # link in dependent modules
    link!.(Ref(mod), dependencies)

    # clean up incompatibilities
    for llvmf in functions(mod)
        # only occurs in debug builds
        delete!(function_attributes(llvmf), EnumAttribute("sspreq", 0, JuliaContext()))

        # make function names safe for ptxas
        # (LLVM ought to do this, see eg. D17738 and D19126), but fails
        # TODO: fix all globals?
        llvmfn = LLVM.name(llvmf)
        if !isdeclaration(llvmf)
            llvmfn′ = safe_fn(llvmf)
            if llvmfn != llvmfn′
                LLVM.name!(llvmf, llvmfn′)
            end
        end
    end

    # rename the entry point
    llvmfn = replace(LLVM.name(entry), r"_\d+$"=>"")
    ## add a friendlier alias
    alias = something(ctx.alias, String(typeof(something(ctx.inner_f, ctx.f)).name.mt.name))
    if startswith(alias, '#')
        alias = "anonymous"
    else
        alias = safe_fn(alias)
    end
    llvmfn = replace(llvmfn, r"_.+" => "_$alias")
    ## append a global unique counter
    global globalUnique
    globalUnique += 1
    llvmfn *= "_$globalUnique"
    LLVM.name!(entry, llvmfn)

    return mod, entry
end

# promote a function to a kernel
# FIXME: sig vs tt (code_llvm vs cufunction)
function promote_kernel!(ctx::CompilerContext, mod::LLVM.Module, entry_f::LLVM.Function)
    kernel = wrap_entry!(ctx, mod, entry_f)

    # property annotations
    # TODO: belongs in irgen? doesn't maxntidx doesn't appear in ptx code?

    annotations = LLVM.Value[kernel]

    ## kernel metadata
    append!(annotations, [MDString("kernel"), ConstantInt(Int32(1), JuliaContext())])

    ## expected CTA sizes
    for (typ,vals) in (:req=>ctx.minthreads, :max=>ctx.maxthreads)
        if vals != nothing
            bounds = CUDAdrv.CuDim3(vals)
            for dim in (:x, :y, :z)
                bound = getfield(bounds, dim)
                append!(annotations, [MDString("$(typ)ntid$(dim)"),
                                      ConstantInt(Int32(bound), JuliaContext())])
            end
        end
    end

    if ctx.blocks_per_sm != nothing
        append!(annotations, [MDString("minctasm"),
                              ConstantInt(Int32(ctx.blocks_per_sm), JuliaContext())])
    end

    if ctx.maxregs != nothing
        append!(annotations, [MDString("maxnreg"),
                              ConstantInt(Int32(ctx.maxregs), JuliaContext())])
    end

    push!(metadata(mod), "nvvm.annotations", MDNode(annotations))

    return kernel
end

# generate a kernel wrapper to fix & improve argument passing
function wrap_entry!(ctx::CompilerContext, mod::LLVM.Module, entry_f::LLVM.Function)
    entry_ft = eltype(llvmtype(entry_f))
    @assert return_type(entry_ft) == LLVM.VoidType(JuliaContext())

    # filter out ghost types, which don't occur in the LLVM function signatures
    sig = Base.signature_type(ctx.f, ctx.tt)
    julia_types = filter(dt->!isghosttype(dt), sig.parameters)

    # generate the wrapper function type & definition
    global globalUnique
    function wrapper_type(julia_t, codegen_t)
        if !isbitstype(julia_t)
            # don't pass jl_value_t by value; it's an opaque structure
            return codegen_t
        elseif isa(codegen_t, LLVM.PointerType) && !(julia_t <: Ptr)
            # we didn't specify a pointer, but codegen passes one anyway.
            # make the wrapper accept the underlying value instead.
            return eltype(codegen_t)
        else
            return codegen_t
        end
    end
    wrapper_types = LLVM.LLVMType[wrapper_type(julia_t, codegen_t)
                                  for (julia_t, codegen_t)
                                  in zip(julia_types, parameters(entry_ft))]
    wrapper_fn = replace(LLVM.name(entry_f), r"^.+?_"=>"ptxcall_") # change the CC tag
    wrapper_ft = LLVM.FunctionType(LLVM.VoidType(JuliaContext()), wrapper_types)
    wrapper_f = LLVM.Function(mod, wrapper_fn, wrapper_ft)

    # emit IR performing the "conversions"
    Builder(JuliaContext()) do builder
        entry = BasicBlock(wrapper_f, "entry", JuliaContext())
        position!(builder, entry)

        wrapper_args = Vector{LLVM.Value}()

        # perform argument conversions
        codegen_types = parameters(entry_ft)
        wrapper_params = parameters(wrapper_f)
        for (julia_t, codegen_t, wrapper_t, wrapper_param) in
            zip(julia_types, codegen_types, wrapper_types, wrapper_params)
            if codegen_t != wrapper_t
                # the wrapper argument doesn't match the kernel parameter type.
                # this only happens when codegen wants to pass a pointer.
                @assert isa(codegen_t, LLVM.PointerType)
                @assert eltype(codegen_t) == wrapper_t

                # copy the argument value to a stack slot, and reference it.
                ptr = alloca!(builder, wrapper_t)
                if LLVM.addrspace(codegen_t) != 0
                    ptr = addrspacecast!(builder, ptr, codegen_t)
                end
                store!(builder, wrapper_param, ptr)
                push!(wrapper_args, ptr)

                # Julia marks parameters as TBAA immutable;
                # this is incompatible with us storing to a stack slot, so clear TBAA
                # TODO: tag with alternative information (eg. TBAA, or invariant groups)
                entry_params = collect(parameters(entry_f))
                candidate_uses = []
                for param in entry_params
                    append!(candidate_uses, collect(uses(param)))
                end
                while !isempty(candidate_uses)
                    usepair = popfirst!(candidate_uses)
                    inst = user(usepair)

                    md = metadata(inst)
                    if haskey(md, LLVM.MD_tbaa)
                        delete!(md, LLVM.MD_tbaa)
                    end

                    # follow along certain pointer operations
                    if isa(inst, LLVM.GetElementPtrInst) ||
                       isa(inst, LLVM.BitCastInst) ||
                       isa(inst, LLVM.AddrSpaceCastInst)
                        append!(candidate_uses, collect(uses(inst)))
                    end
                end
            else
                push!(wrapper_args, wrapper_param)
            end
        end

        call!(builder, entry_f, wrapper_args)

        ret!(builder)
    end

    # early-inline the original entry function into the wrapper
    push!(function_attributes(entry_f), EnumAttribute("alwaysinline", 0, JuliaContext()))
    linkage!(entry_f, LLVM.API.LLVMInternalLinkage)
    ModulePassManager() do pm
        always_inliner!(pm)
        run!(pm, mod)
    end

    return wrapper_f
end

const libdevices = Dict{String, LLVM.Module}()
function link_libdevice!(ctx::CompilerContext, mod::LLVM.Module)
    CUDAnative.configured || return

    # find libdevice
    path = if isa(libdevice, Dict)
        # select the most recent & compatible library
        vers = keys(CUDAnative.libdevice)
        compat_vers = Set(ver for ver in vers if ver <= ctx.cap)
        isempty(compat_vers) && error("No compatible CUDA device library available")
        ver = maximum(compat_vers)
        path = libdevice[ver]
    else
        libdevice
    end

    # load the library, once
    if !haskey(libdevices, path)
        open(path) do io
            libdevice_mod = parse(LLVM.Module, read(io), JuliaContext())
            name!(libdevice_mod, "libdevice")
            libdevices[path] = libdevice_mod
        end
    end
    libdevice_mod = LLVM.Module(libdevices[path])

    # override libdevice's triple and datalayout to avoid warnings
    triple!(libdevice_mod, triple(mod))
    datalayout!(libdevice_mod, datalayout(mod))

    # 1. save list of external functions
    exports = map(LLVM.name, functions(mod))
    filter!(fn->!haskey(functions(libdevice_mod), fn), exports)

    # 2. link with libdevice
    link!(mod, libdevice_mod)

    ModulePassManager() do pm
        # 3. internalize all functions not in list from (1)
        internalize!(pm, exports)

        # 4. eliminate all unused internal functions
        #
        # this isn't necessary, as we do the same in optimize! to inline kernel wrappers,
        # but it results _much_ smaller modules which are easier to handle on optimize=false
        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)

        # 5. run NVVMReflect pass
        push!(metadata(mod), "nvvm-reflect-ftz",
              MDNode([ConstantInt(Int32(1), JuliaContext())]))

        # 6. run standard optimization pipeline
        #
        #    see `optimize!`

        run!(pm, mod)
    end
end

function machine(cap::VersionNumber, triple::String)
    InitializeNVPTXTarget()
    InitializeNVPTXTargetInfo()
    t = Target(triple)

    InitializeNVPTXTargetMC()
    cpu = "sm_$(cap.major)$(cap.minor)"
    if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
        # in the case of CUDA 9, we use sync intrinsics from PTX ISA 6.0+
        feat = "+ptx60"
    else
        feat = ""
    end
    tm = TargetMachine(t, triple, cpu, feat)

    return tm
end

# Optimize a bitcode module according to a certain device capability.
function optimize!(ctx::CompilerContext, mod::LLVM.Module, entry::LLVM.Function)
    tm = machine(ctx.cap, triple(mod))

    # GPU code is _very_ sensitive to register pressure and local memory usage,
    # so forcibly inline every function definition into the entry point
    # and internalize all other functions (enabling ABI-breaking optimizations).
    # FIXME: this is too coarse. use a proper inliner tuned for GPUs
    ModulePassManager() do pm
        no_inline = EnumAttribute("noinline", 0, JuliaContext())
        always_inline = EnumAttribute("alwaysinline", 0, JuliaContext())
        for f in filter(f->f!=entry && !isdeclaration(f), functions(mod))
            attrs = function_attributes(f)
            if !(no_inline in collect(attrs))
                push!(attrs, always_inline)
            end
            linkage!(f, LLVM.API.LLVMInternalLinkage)
        end
        always_inliner!(pm)
        run!(pm, mod)
    end

    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        ccall(:jl_add_optimization_passes, Cvoid,
              (LLVM.API.LLVMPassManagerRef, Cint),
              LLVM.ref(pm), Base.JLOptions().opt_level)

        # NVPTX's target machine info enables runtime unrolling,
        # but Julia's pass sequence only invokes the simple unroller.
        loop_unroll!(pm)
        instruction_combining!(pm)  # clean-up redundancy
        licm!(pm)                   # the inner runtime check might be outer loop invariant

        cfgsimplification!(pm)

        # CUDAnative's compiler internalizes non-inlined child functions, making it possible
        # to rewrite them (whereas the Julia JIT caches those functions exactly);
        # this opens up some more optimization opportunities
        dead_arg_elimination!(pm)   # parent doesn't use return value --> ret void

        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)

        # FIXME: see #195
        PassManagerBuilder() do pmb
            populate!(pm, pmb)
        end

        run!(pm, mod)
    end
end

function mcgen(ctx::CompilerContext, mod::LLVM.Module, f::LLVM.Function)
    tm = machine(ctx.cap, triple(mod))

    InitializeNVPTXAsmPrinter()
    return String(emit(tm, mod, LLVM.API.LLVMAssemblyFile))
end

# Compile a function to PTX, returning the assembly and an entry point.
# Not to be used directly, see `cufunction` instead.
function compile_function(ctx::CompilerContext)
    ## high-level code generation (Julia AST)

    @debug "(Re)compiling function" ctx

    validate_invocation(ctx)

    ## low-level code generation (LLVM IR)

    mod, entry = irgen(ctx)

    if ctx.kernel
        entry = promote_kernel!(ctx, mod, entry)
    end

    @trace("Module entry point: ", LLVM.name(entry))

    # link libdevice, if it might be necessary
    # TODO: should be more find-grained -- only matching functions actually in this libdevice
    if any(f->isdeclaration(f) && intrinsic_id(f)==0, functions(mod))
        link_libdevice!(ctx, mod)
    end

    # optimize the IR
    optimize!(ctx, mod, entry)

    # make sure any non-isbits arguments are unused
    real_arg_i = 0
    sig = Base.signature_type(ctx.f, ctx.tt)
    for (arg_i,dt) in enumerate(sig.parameters)
        isghosttype(dt) && continue
        real_arg_i += 1

        if !isbitstype(dt)
            param = parameters(entry)[real_arg_i]
            if !isempty(uses(param))
                compiler_error(ctx, "passing and using non-bitstype argument";
                               argument=arg_i, argument_type=dt)
            end
        end
    end

    # validate generated IR
    errors = unique(validate_ir(mod))
    if !isempty(errors)
        compiler_error(ctx, "unsupported LLVM IR"; errors=errors)
    end

    ## machine code generation (PTX assembly)

    module_asm = mcgen(ctx, mod, entry)

    return module_asm, LLVM.name(entry)
end

# (::CompilerContext)
const compile_hook = Ref{Union{Nothing,Function}}(nothing)

# Main entry point for compiling a Julia function + argtypes to a callable CUDA function
function cufunction(dev::CuDevice, @nospecialize(f), @nospecialize(tt); kwargs...)
    CUDAnative.configured || error("CUDAnative.jl has not been configured; cannot JIT code.")
    @assert isa(f, Core.Function)

    # select a capability level
    dev_cap = capability(dev)
    compat_caps = filter(cap -> cap <= dev_cap, target_support)
    isempty(compat_caps) &&
        error("Device capability v$dev_cap not supported by available toolchain")
    cap = maximum(compat_caps)

    ctx = CompilerContext(f, tt, cap, true; kwargs...)

    if compile_hook[] != nothing
        global globalUnique
        previous_globalUnique = globalUnique
        compile_hook[](ctx)
        globalUnique = previous_globalUnique
    end

    (module_asm, module_entry) = compile_function(ctx)

    # enable debug options based on Julia's debug setting
    jit_options = Dict{CUDAdrv.CUjit_option,Any}()
    if Base.JLOptions().debug_level == 1
        jit_options[CUDAdrv.GENERATE_LINE_INFO] = true
    elseif Base.JLOptions().debug_level >= 2
        jit_options[CUDAdrv.GENERATE_DEBUG_INFO] = true
    end
    cuda_mod = CuModule(module_asm, jit_options)
    cuda_fun = CuFunction(cuda_mod, module_entry)

    @debug begin
        attr = attributes(cuda_fun)
        bin_ver = VersionNumber(divrem(attr[CUDAdrv.FUNC_ATTRIBUTE_BINARY_VERSION],10)...)
        ptx_ver = VersionNumber(divrem(attr[CUDAdrv.FUNC_ATTRIBUTE_PTX_VERSION],10)...)
        regs = attr[CUDAdrv.FUNC_ATTRIBUTE_NUM_REGS]
        local_mem = attr[CUDAdrv.FUNC_ATTRIBUTE_LOCAL_SIZE_BYTES]
        shared_mem = attr[CUDAdrv.FUNC_ATTRIBUTE_SHARED_SIZE_BYTES]
        constant_mem = attr[CUDAdrv.FUNC_ATTRIBUTE_CONST_SIZE_BYTES]
        """Compiled $f to PTX $ptx_ver for SM $bin_ver using $regs registers.
           Memory usage: $local_mem B local, $shared_mem B shared, $constant_mem B constant"""
    end

    return cuda_fun, cuda_mod
end

function init_jit()
    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("--nvptx-fma-level=1")
end


text:
# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-pinfo
# https://searchcode.com/file/21350024/blacs/BLACS-openmpi/SRC/MPI/blacs_pinfo_.c
# input : nothing
# output: rank, num of processes
function pinfo()
    mypnum, nprocs = zeros(ScaInt,1), zeros(ScaInt,1)
    ccall((:blacs_pinfo_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{ScaInt}),
        mypnum, nprocs)
    return mypnum[1], nprocs[1]
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-get
# https://searchcode.com/file/21350069/blacs/BLACS-openmpi/SRC/MPI/blacs_get_.c
# input : BLACS context ( MPI communicator )
#         BLACS indicator for BLACS contxt which should be returned
# output: BLACS context
function get(icontxt::ScaInt, what::ScaInt)
    val = zeros(ScaInt,1)
    ccall((:blacs_get_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
        Ref(icontxt), Ref(what), val)
    return val[1]
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-gridinit
# https://searchcode.com/file/21350091/blacs/BLACS-openmpi/SRC/MPI/blacs_gridinit_.c
# input : BLACS context
#         BLACS layout; how to map processes to BLACS grid
#         num of rows and cols in the process grid
# output: BLACS context
function gridinit(icontxt::ScaInt, layout::Char, nprow::ScaInt, npcol::ScaInt)
    ocontxt = Array{ScaInt,1}([icontxt])
    ccall((:blacs_gridinit_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{Char}, Ptr{ScaInt}, Ptr{ScaInt}),
        ocontxt, f_pchar(layout), Ref(nprow), Ref(npcol))
    return ocontxt[1]
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-gridinfo
# https://searchcode.com/file/21350044/blacs/BLACS-openmpi/SRC/MPI/blacs_gridinfo_.c
# input : BLACS context
# output: num of process rows and cols in the current process id
#         row and col coordinate in the current process id
function gridinfo(icontxt::ScaInt)
    nprow = zeros(ScaInt,1)
    npcol = zeros(ScaInt,1)
    myrow = zeros(ScaInt,1)
    mycol = zeros(ScaInt,1)
    ccall((:blacs_gridinfo_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
        Ref(icontxt), nprow, npcol, myrow, mycol)
    return nprow[1], npcol[1], myrow[1], mycol[1]
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-pnum
# https://searchcode.com/file/21350042/blacs/BLACS-openmpi/SRC/MPI/blacs_pnum_.c
# input : BLACS context
#         row and col coordinate in a process id
# output: process number
function pnum(icontxt::ScaInt, prow::ScaInt, pcol::ScaInt)
    return ccall((:blacs_pnum_, libscalapack), ScaInt,
        (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
        Ref(icontxt), Ref(prow), Ref(pcol))
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-pcoord
# https://searchcode.com/file/21350056/blacs/BLACS-openmpi/SRC/MPI/blacs_pcoord_.c
# input : BLACS context, process number
# output: row and col coordinate in the pnum process id
function pcoord(icontxt::ScaInt, pnum::ScaInt)
    prow = zeros(ScaInt,1)
    pcol = zeros(ScaInt,1)
    ccall((:blacs_pcoord_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
        Ref(icontxt), Ref(pnum), prow, pcol)
    return prow[1], pcol[1]
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-barrier
# https://searchcode.com/file/21350053/blacs/BLACS-openmpi/SRC/MPI/blacs_barrier_.c
# input : BLACS context
#         BLACS parameter for row/col/all grid to participate barrier
# output: nothing
function barrier(icontxt::ScaInt, scope::Char)
    ccall((:blacs_barrier_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{Char}),
        Ref(icontxt), f_pchar(scope))
end

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-gridexit
# https://searchcode.com/file/21350036/blacs/BLACS-openmpi/SRC/MPI/blacs_gridexit_.c
# input : BLACS context
# output: nothing
gridexit(icontxt::ScaInt) = ccall((:blacs_gridexit_, libscalapack), Nothing, (Ptr{ScaInt},), Ref(icontxt))

# https://software.intel.com/en-us/mkl-developer-reference-c-blacs-exit
# https://searchcode.com/file/21349846/blacs/BLACS-openmpi/SRC/MPI/blacs_exit_.c
# input : BLACS flag whether continue message passing or not after done
# output: nothing
exit(continue_::ScaInt = 0) = ccall((:blacs_exit_, libscalapack), Nothing, (Ptr{ScaInt},), Ref(continue_))


text:
function test_annotation_row(recording, id, span; custom...)
    row = @compat (; recording, id, span)
    row_with_custom = (; row..., custom...)

    @test Arrow.Table(Legolas.tobuffer([row], Legolas.Schema("onda.annotation@1"); validate=true)) isa Arrow.Table
    @test Arrow.Table(Legolas.tobuffer([row_with_custom], Legolas.Schema("onda.annotation@1"); validate=true)) isa Arrow.Table

    # intended normalization of input fields for constructor
    recording::UUID = recording isa UUID ? recording : UUID(recording)
    id::UUID = id isa UUID ? id : UUID(id)
    span::TimeSpan = TimeSpan(span)
    norm_row = @compat (; recording, id, span)
    norm_row_with_custom = (; norm_row..., custom...)

    @test has_rows(Annotation(row), norm_row)
    @test has_rows(Annotation(row_with_custom), norm_row_with_custom)
    @test has_rows(Annotation(Annotation(row)), norm_row)
    @test has_rows(Annotation(Annotation(row_with_custom)), norm_row_with_custom)
    @test has_rows(Annotation(Tables.Row(row)), norm_row)
    @test has_rows(Annotation(Tables.Row(row_with_custom)), norm_row_with_custom)
    @test has_rows(Annotation(; row...), norm_row)
    @test has_rows(Annotation(; row..., custom...), norm_row_with_custom)
end

@testset "`Annotation` construction/access" begin
    custom = (a="test", b=1, c=[2.0, 3.0])
    test_annotation_row(UInt128(uuid4()), UInt128(uuid4()), (start=Nanosecond(1), stop=Nanosecond(100)); custom...)
    test_annotation_row(uuid4(), uuid4(), TimeSpan(Nanosecond(1), Nanosecond(100)); custom...)
end

@testset "`onda.annotation` validation" begin
    template = (recording=uuid4(), id=uuid4(), span=TimeSpan(0, 1), custom=1234)
    @test Annotation(template) isa Annotation
    good = [template, Tables.rowmerge(template; id=uuid4()), Tables.rowmerge(template; id=uuid4())]
    @test validate_annotations(good) === good
    @test_throws ArgumentError validate_annotations(vcat(good, template))
    @test_throws ArgumentError validate_annotations([template, template, template])
    @test_throws ArgumentError validate_annotations((x=[1, 2, 3], y=["lol", "bad", "table"]))
end

@testset "`merge_overlapping_annotations`" begin
    recs = (uuid4(), uuid4(), uuid4())
    sources = [#= 1 =#  Annotation(recording=recs[1], id=uuid4(), span=TimeSpan(0, 100)),
               #= 2 =#  Annotation(recording=recs[2], id=uuid4(), span=TimeSpan(55, 100)),
               #= 3 =#  Annotation(recording=recs[1], id=uuid4(), span=TimeSpan(34, 76)),
               #= 4 =#  Annotation(recording=recs[1], id=uuid4(), span=TimeSpan(120, 176)),
               #= 5 =#  Annotation(recording=recs[2], id=uuid4(), span=TimeSpan(67, 95)),
               #= 6 =#  Annotation(recording=recs[2], id=uuid4(), span=TimeSpan(15, 170)),
               #= 7 =#  Annotation(recording=recs[1], id=uuid4(), span=TimeSpan(43, 89)),
               #= 8 =#  Annotation(recording=recs[3], id=uuid4(), span=TimeSpan(0, 50)),
               #= 9 =#  Annotation(recording=recs[2], id=uuid4(), span=TimeSpan(2, 10)),
               #= 10 =# Annotation(recording=recs[1], id=uuid4(), span=TimeSpan(111, 140)),
               #= 11 =# Annotation(recording=recs[3], id=uuid4(), span=TimeSpan(60, 100)),
               #= 12 =# Annotation(recording=recs[3], id=uuid4(), span=TimeSpan(23, 80)),
               #= 13 =# Annotation(recording=recs[3], id=uuid4(), span=TimeSpan(100, 110)),
               #= 14 =# Annotation(recording=recs[1], id=uuid4(), span=TimeSpan(200, 300))]
    merged = Tables.columns(merge_overlapping_annotations(sources))
    @test Tables.columnnames(merged) == (:recording, :id, :span, :from)
    sources_id = [row.id for row in sources]
    @test !any(in(id, sources_id) for id in merged.id)
    merged = @compat Set(Tables.rowtable((; merged.recording, merged.span, merged.from)))
    expected = Set([(recording=recs[1], span=TimeSpan(0, 100), from=[sources[1].id, sources[3].id, sources[7].id]),
                    (recording=recs[1], span=TimeSpan(111, 176), from=[sources[10].id, sources[4].id]),
                    (recording=recs[1], span=TimeSpan(200, 300), from=[sources[14].id]),
                    (recording=recs[2], span=TimeSpan(15, 170), from=[sources[6].id, sources[2].id, sources[5].id]),
                    (recording=recs[2], span=TimeSpan(2, 10), from=[sources[9].id]),
                    (recording=recs[3], span=TimeSpan(100, 110), from=[sources[13].id]),
                    (recording=recs[3], span=TimeSpan(0, 100), from=[sources[8].id, sources[12].id, sources[11].id])])
    @test expected == merged
end


text:
#!/usr/bin/env julia
using XAM,Statistics,GenomicFeatures

bamName = "PC-HLA-EG4.bam"

bamIn = open(BAM.Reader,bamName, index="$bamName.bai")
targetChr="chr1"

numAligns = 0
numAligns+=1

allquals=Int[]
for record in eachoverlap(bamIn,targetChr,1:100000)
	numAligns+=1
	push!(allquals,BAM.mappingquality(record))
end
close(bamIn)

println("Num aligns: ",numAligns)
println("Mean quality: ",mean(allquals))

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
function StaticRanges.similar_type(
    ia::IndicesArray;
    axes_type::Type=axes_type(A),
    parent_type::Type=parent_type(A)
   )
    return _similar_type(ia, parent_type, axes_type)
end

function _similar_type(
    a::IndicesArray;
    new_axes_type::Type=axes_type(a),
    new_parent_type::Type=parent_type(a)
   )
    return IndicesArray{eltype(new_parent_type),ndims(new_parent_type),new_parent_type,new_axes_type}
end

#=
function Base.similar(
    a::IndicesArray{T},
    eltype::Type=T,
    axs=axes(a)
   ) where {T}
    return IndicesArray(similar(parent(a), eltype, map(length, axs)), _drop_empty(axs))
end
=#

function Base.similar(
    a::AbstractArray{T},
    eltype::Type,
    axs::Tuple{Vararg{<:AbstractIndex}}
   ) where {T}
    return rebuild(a, similar(a, eltype, map(length, axs)), _drop_empty(axs))
end

function Base.similar(
    ::Type{A},
    eltype::Type,
    axs::Tuple{Vararg{<:AbstractIndex}}
   ) where {A<:AbstractArray}
    return rebuild(A, similar(A, eltype, map(length, axs)), _drop_empty(axs))
end

function Base.similar(
    ::Type{A},
    axs::Tuple{Vararg{<:AbstractIndex}}
   ) where {A<:AbstractArray}
    return rebuild(A, similar(A, map(length, axs)), _drop_empty(axs))
end


text:
using Test
using QBase

@testset "./src/channels.jl" begin

@testset "repalacer_channel()" begin
    @testset "simple qubit examples" begin
        ρ = State([1 0;0 0])
        σ = State([0 0;0 1])

        r = replacer_channel(ρ, σ, 0.5)

        @test r isa State
        @test r == [0.5 0;0 0.5]
    end
    @testset "errors" begin
        ρ = State([1 0;0 0])
        σ = State([0 0;0 1])

        @test_throws DomainError replacer_channel(ρ, σ, -0.1)
        @test_throws DomainError replacer_channel(ρ, σ, 1.1)
        @test_throws DomainError replacer_channel(
            State([1 0;0 0]), State([0 0 0;0 1 0;0 0 0]), 0.5
        )
    end
end

@testset "depolarizing_channel()" begin
    @testset "simple State examples" begin
        ρ0 = State([1 0;0 0])

        ρ0_out1 = depolarizing_channel(ρ0, 1)

        @test ρ0_out1 isa State
        @test ρ0_out1 == ρ0

        ρ0_out2 = depolarizing_channel(ρ0, 0)

        @test ρ0_out2 isa State
        @test ρ0_out2 == [1/2 0;0 1/2]

        ρ0_out3 = depolarizing_channel(ρ0, 0.5)

        @test ρ0_out3 isa State
        @test ρ0_out3 == [3/4 0;0 1/4]
    end

    @testset "simple State examples" begin
        ρ = State([0.5 0 0.5;0 0 0;0.5 0 0.5])

        ρ_out1 = depolarizing_channel(ρ, 1)

        @test ρ_out1 isa State
        @test ρ_out1 == ρ

        ρ_out2 = depolarizing_channel(ρ, 0.5)

        @test ρ_out2 isa State
        @test ρ_out2 == [(1/4+1/6) 0 1/4;0 1/6 0;1/4 0 (1/4+1/6)]
    end

    @testset "errors" begin
        ρ = State([1 0;0 0])

        @test_throws DomainError depolarizing_channel(ρ, 1.1)
        @test_throws DomainError depolarizing_channel(ρ, -0.1)
    end
end

@testset "erasure_channel()" begin
    @testset "simple State examples" begin
        ρ = State([1 0;0 0])

        ρ_out1 = erasure_channel(ρ, 1)

        @test ρ_out1 isa State
        @test ρ_out1 == [1 0 0;0 0 0;0 0 0]

        ρ_out2 = erasure_channel(ρ, 0)

        @test ρ_out2 isa State
        @test ρ_out2 == [0 0 0;0 0 0;0 0 1]

        ρ_out3 = erasure_channel(ρ, 0.5)

        @test ρ_out3 isa State
        @test ρ_out3 == [0.5 0 0;0 0 0;0 0 0.5]
    end

    @testset "simple State examples" begin
        ρ = State([0.5 0 0.5;0 0 0;0.5 0 0.5])

        ρ_out1 = erasure_channel(ρ, 1)

        @test ρ_out1 isa State
        @test ρ_out1 == [0.5 0 0.5 0;0 0 0 0;0.5 0 0.5 0;0 0 0 0]

        ρ_out2 = erasure_channel(ρ, 0.5)

        @test ρ_out2 isa State
        @test ρ_out2 == [1/4 0 1/4 0;0 0 0 0;1/4 0 1/4 0;0 0 0 1/2]
    end

    @testset "errors" begin
        ρ = State([1 0;0 0])

        @test_throws DomainError erasure_channel(ρ, 1.1)
        @test_throws DomainError erasure_channel(ρ, -0.1)
    end
end

end


text:
using DateParsing
using Base.Test

include("perf.jl")


text:
"""
    update_ista!(B, resid, grad, X, Y, Z, norms, lambda, reg, stepsize)

Updates coefficient estimates in place for each ISTA iteration when `X` and
`Z` are both standardized.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- grad = 2d array of floats consisting of the gradient
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response observations
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = `nothing`
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- stepsize = 1d array consisting of a float; step size of updates

# Value

None; updates coefficients in place

"""
function update_ista!(B::AbstractArray{Float64,2},
                      resid::AbstractArray{Float64,2},
                      grad::AbstractArray{Float64,2},
                      X::AbstractArray{Float64,2},
                      Y::AbstractArray{Float64,2},
                      Z::AbstractArray{Float64,2},
                      norms::Nothing, lambda::Float64, reg::BitArray{2},
                      stepsize::AbstractArray{Float64,1})

    # Update gradient
    calc_grad!(grad, X, Z, resid)

    # Cycle through all the coefficients to perform assignments
    for j = 1:size(B,2), i = 1:size(B,1)
        b2update = B[i,j] - stepsize[1] * grad[i,j] # L2 update
        b2sign = sign(b2update) # L2 update sign

        # Update the current coefficient
        b = reg[i,j] ? prox(B[i,j], grad[i,j], b2sign, lambda, norms,
                            stepsize[1]) : b2update

        # Update coefficient estimate in the current copy
        if b != B[i,j]
            B[i,j] = b
        end
    end

    # Update residuals
    calc_resid!(resid, X, Y, Z, B)
end

"""
    update_ista!(B, resid, grad, X, Y, Z, norms, lambda, reg, stepsize)

Updates coefficient estimates in place for each ISTA iteration when `X` and
`Z` are not standardized.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- grad = 2d array of floats consisting of the gradient
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response observations
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- stepsize = 1d array consisting of a float; step size of updates

# Value

None; updates coefficients in place

"""
function update_ista!(B::AbstractArray{Float64,2},
                      resid::AbstractArray{Float64,2},
                      grad::AbstractArray{Float64,2},
                      X::AbstractArray{Float64,2},
                      Y::AbstractArray{Float64,2},
                      Z::AbstractArray{Float64,2},
                      norms::AbstractArray{Float64,2},
                      lambda::Float64, reg::BitArray{2},
                      stepsize::AbstractArray{Float64,1})

    # Update gradient
    calc_grad!(grad, X, Z, resid)

    # Cycle through all the coefficients to perform assignments
    for j = 1:size(B,2), i = 1:size(B,1)
        b2update = B[i,j] - stepsize[1] * grad[i,j]./norms[i,j] # L2 updates
        b2sign = sign(b2update) # L2 update signs

        # Update the current coefficient
        b = reg[i,j] ? prox(B[i,j], grad[i,j], b2sign, lambda, norms[i,j],
                            stepsize[1]) : b2update

        # Update coefficient estimate in the current copy
        if b != B[i,j]
            B[i,j] = b
        end
    end

    # Update residuals
    calc_resid!(resid, X, Y, Z, B)
end

"""
    ista!(X, Y, Z, lambda, B, regXidx, regZidx, reg, norms;
          isVerbose, stepsize, thresh, maxiter)

Performs ISTA with fixed step size.

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- lambda = lambda penalty, a floating scalar
- B = 2d array of floats consisting of starting coefficient estimates
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`

# Keyword arguments

- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates. Defaults to `0.01`.
- thresh = threshold at which the coefficients are considered to have
  converged, a floating scalar. Defaults to `10^(-7)`.
- maxiter = maximum number of update iterations. Defaults to `10^10`.

# Value

None; updates coefficients in place

# Some notes

Convergence is determined as when the log ratio of the current and previous
criteria is less than the threshold `thresh`.

The default method for choosing the fixed step size for `ista!` is to use the
reciprocal of the product of the maximum eigenvalues of `X*transpose(X)`
and `Z*transpose(Z)`. This is computed by the `mlmnet` function when `ista!`
is passed into the `fun` argument and `setStepSize` is set to `true`. If
`setStepSize` is set to `false`, the value of the `stepsize` argument will be
used as the fixed step size. Note that obtaining the eigenvalues when `X`
and/or `Z` are very large may exceed computational limitations.

"""
function ista!(X::AbstractArray{Float64,2}, Y::AbstractArray{Float64,2},
               Z::AbstractArray{Float64,2}, lambda::Float64,
               B::AbstractArray{Float64,2},
               regXidx::AbstractArray{Int64,1},
               regZidx::AbstractArray{Int64,1}, reg::BitArray{2}, norms;
               isVerbose::Bool=true, stepsize::Float64=0.01,
               thresh::Float64=10.0^(-7), maxiter::Int=10^10)

    # Pre-allocating arrays to store residuals using and gradient.
    resid = Array{Float64}(undef, size(Y))
    grad = Array{Float64}(undef, size(B))

    # Calculate residuals
    calc_resid!(resid, X, Y, Z, B)

    # Denominators of criterion
    crit_denom = [size(X,1)*size(Z,1), size(X,2)*size(Z,2)]
    # Placeholder to store the old criterion
    oldcrit = 1.0
    # Calculate the current criterion
    crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

    stepsize = [stepsize]
    iter = 0
    # Iterate until coefficients converge or maximum iterations have been
    # reached.
    while (abs(log(crit/oldcrit)) > thresh) && (iter < maxiter)
        # Store the current criterion
        oldcrit = crit
        # Update the coefficient estimates
        update_ista!(B, resid, grad, X, Y, Z, norms, lambda, reg, stepsize)
        # Calculate the criterion after updating
        crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

        iter += 1 # Increment the number of iterations
        # Warning message if coefficient estimates do not converge.
        if iter==(maxiter)
            println("Estimates did not converge.")
        end
    end

    println_verbose(string("Criterion: ", crit), isVerbose)
    println_verbose(string("Number of iterations: ", iter), isVerbose)
end


text:
abstract type AbstractHydroFormulation <: AbstractDeviceFormulation end
abstract type AbstractHydroDispatchFormulation <: AbstractHydroFormulation end
abstract type AbstractHydroUnitCommitment <: AbstractHydroFormulation end
abstract type AbstractHydroReservoirFormulation <: AbstractHydroDispatchFormulation end
struct HydroDispatchRunOfRiver <: AbstractHydroDispatchFormulation end
struct HydroDispatchReservoirFlow <: AbstractHydroReservoirFormulation end
struct HydroDispatchReservoirStorage <: AbstractHydroReservoirFormulation end
#=
# Commenting out all Unit Commitment formulations as all Hydro UC
# formulations are currently not supported
struct HydroCommitmentRunOfRiver <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirFlow <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirStorage <: AbstractHydroUnitCommitment end
=#
########################### Hydro generation variables #################################
function activepower_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
) where {H <: PSY.HydroGen}
    add_variable(
        psi_container,
        devices,
        variable_name(ACTIVE_POWER, H),
        false,
        :nodal_balance_active;
        lb_value = d -> PSY.get_activepowerlimits(d).min,
        ub_value = d -> PSY.get_activepowerlimits(d).max,
        init_value = d -> PSY.get_activepower(d),
    )

    return
end

function reactivepower_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
) where {H <: PSY.HydroGen}
    add_variable(
        psi_container,
        devices,
        variable_name(REACTIVE_POWER, H),
        false,
        :nodal_balance_reactive;
        ub_value = d -> PSY.get_reactivepowerlimits(d).max,
        lb_value = d -> PSY.get_reactivepowerlimits(d).min,
        init_value = d -> PSY.get_reactivepower(d),
    )

    return
end

function energy_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
) where {H <: PSY.HydroGen}
    add_variable(
        psi_container,
        devices,
        variable_name(ENERGY, H),
        false;
        ub_value = d -> PSY.get_storage_capacity(d),
        lb_value = d -> 0.0,
        init_value = d -> PSY.get_initial_storage(d),
    )

    return
end

#=
function inflow_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
) where {H <: PSY.HydroGen}
    add_variable(
        psi_container,
        devices,
        variable_name(INFLOW, H),
        false;
        ub_value = d -> PSY.get_inflow(d),
        lb_value = d -> 0.0,
    )

    return
end
=#

function spillage_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
) where {H <: PSY.HydroGen}
    add_variable(
        psi_container,
        devices,
        variable_name(SPILLAGE, H),
        false;
        lb_value = d -> 0.0,
    )

    return
end

"""
This function add the variables for power generation commitment to the model
"""
#=
function commitment_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
) where {H <: PSY.HydroGen}
    time_steps = model_time_steps(psi_container)
    var_names = [variable_name(ON, H), variable_name(START, H), variable_name(STOP, H)]

    for v in var_names
        add_variable(psi_container, devices, v, true)
    end

    return
end

# All Hydro UC formulations are currently not supported
### Constraints for Thermal Generation without commitment variables ####
"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitment_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H,D},
    system_formulation::Type{S},
    feedforward::Union{Nothing,AbstractAffectFeedForward},
) where {H<:PSY.HydroGen,D<:AbstractHydroUnitCommitment,S<:PM.AbstractPowerModel}
    device_commitment(
        psi_container,
        get_initial_conditions(psi_container, ICKey(DeviceStatus, H)),
        constraint_name(COMMITMENT, H),
        (variable_name(START, H), variable_name(STOP, H), variable_name(ON, H)),
    )

    return
end
=#
####################################### Reactive Power Constraints #########################
function reactivepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, D},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroGen, D <: AbstractHydroDispatchFormulation}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limits = PSY.get_reactivepowerlimits(d)
        name = PSY.get_name(d)
        range_data = DeviceRange(name, limits)
        constraint_data[ix] = range_data
    end

    device_range(
        psi_container,
        constraint_data,
        constraint_name(REACTIVE_RANGE, H),
        variable_name(REACTIVE_POWER, H),
    )
    return
end

function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, <:AbstractHydroDispatchFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroGen}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    if !parameters && !use_forecast_data
        constraint_data = Vector{DeviceRange}(undef, length(devices))
        for (ix, d) in enumerate(devices)
            name = PSY.get_name(d)
            ub = PSY.get_activepower(d)
            limits = (min = 0.0, max = ub)
            range_data = DeviceRange(name, limits)
            add_device_services!(range_data, d, model)
            constraint_data[ix] = range_data
        end
        device_range(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE_RANGE, H),
            variable_name(ACTIVE_POWER, H),
        )
        return
    end

    forecast_label = "get_rating"
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, x -> PSY.get_rating(x), ts_vector)
        add_device_services!(timeseries_data, d, model)
        constraint_data[ix] = timeseries_data
    end

    if parameters
        device_timeseries_param_ub(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE, H),
            UpdateRef{H}(ACTIVE_POWER, forecast_label),
            variable_name(ACTIVE_POWER, H),
        )
    else
        device_timeseries_ub(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE, H),
            variable_name(ACTIVE_POWER, H),
        )
    end
    return
end

function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, <:AbstractHydroReservoirFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroGen}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_activepowerlimits(d)
        range_data = DeviceRange(name, limits)
        add_device_services!(range_data, d, model)
        constraint_data[ix] = range_data
    end
    device_range(
        psi_container,
        constraint_data,
        constraint_name(ACTIVE_RANGE, H),
        variable_name(ACTIVE_POWER, H),
    )
    return
end

#=
# All Hydro UC formulations are currently not supported
function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H,<:AbstractHydroUnitCommitment},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing,AbstractAffectFeedForward},
) where {H<:PSY.HydroGen}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    ts_data_active, constraint_data = get_time_series(
        psi_container,
        devices,
        model,
        x -> PSY.get_activepowerlimits(x),
    )

    if !parameters && !use_forecast_data
        device_semicontinuousrange(
            psi_container,
            constraint_data,
            constraint_name(ACTIVE_RANGE, H),
            variable_name(ACTIVE_POWER, H),
            variable_name(ON, H),
        )
        return
    end

    if parameters
        device_timeseries_ub_bigM(
            psi_container,
            ts_data_active,
            constraint_name(ACTIVE_RANGE, H),
            variable_name(ACTIVE_POWER, H),
            UpdateRef{H}(ON, "get_rating"),
            variable_name(ON, H),
        )
    else
        device_timeseries_ub_bin(
            psi_container,
            ts_data_active,
            constraint_name(ACTIVE_RANGE, H),
            variable_name(ACTIVE_POWER, H),
            variable_name(ON, H),
        )
    end

    return
end
=#
######################## Inflow constraints ############################
#=
# TODO: Determine if this is useful for ROR formulation ?
function inflow_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H,<:AbstractHydroDispatchFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing,AbstractAffectFeedForward},
) where {H<:PSY.HydroGen}

    return
end

function inflow_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{PSY.HydroEnergyReservoir,HydroDispatchReservoirStorage},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing,AbstractAffectFeedForward},
) where {H<:PSY.HydroGen}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    ts_data_inflow, constraint_data = _get_inflow_time_series(
        psi_container,
        devices,
        model,
        x -> (min = 0.0, max = PSY.get_inflow(x)),
    )

    if !parameters && !use_forecast_data
        device_range(
            psi_container,
            constraint_data,
            constraint_name(INFLOW_RANGE, H),
            variable_name(INFLOW, H),
        )
        return
    end

    if parameters
        device_timeseries_param_ub(psi_container,
                            ts_data_inflow,
                            constraint_name(INFLOW_RANGE, H),
                            UpdateRef{H}(INFLOW_RANGE, "get_inflow"),
                            variable_name(INFLOW, H))
    else
        device_timeseries_ub(
            psi_container,
            ts_data_inflow,
            constraint_name(INFLOW_RANGE, H),
            variable_name(INFLOW, H),
        )
    end

    return
end
=#
######################## Energy balance constraints ############################

function energy_balance_constraint!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, HydroDispatchReservoirStorage},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroEnergyReservoir}
    key = ICKey(EnergyLevel, H)
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)

    if !has_initial_conditions(psi_container.initial_conditions, key)
        throw(IS.DataFormatError("Initial Conditions for $(H) Energy Constraints not in the model"))
    end

    forecast_label = "get_inflow"
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, x -> PSY.get_rating(x), ts_vector)
        add_device_services!(timeseries_data, d, model)
        constraint_data[ix] = timeseries_data
    end

    if parameters
        energy_balance_external_input_param(
            psi_container,
            get_initial_conditions(psi_container, key),
            constraint_data,
            constraint_name(ENERGY_CAPACITY, H),
            (
                variable_name(SPILLAGE, H),
                variable_name(ACTIVE_POWER, H),
                variable_name(ENERGY, H),
            ),
            UpdateRef{H}(INFLOW, forecast_label),
        )
    else
        energy_balance_external_input(
            psi_container,
            get_initial_conditions(psi_container, key),
            constraint_data,
            constraint_name(ENERGY_CAPACITY, H),
            (
                variable_name(SPILLAGE, H),
                variable_name(ACTIVE_POWER, H),
                variable_name(ENERGY, H),
            ),
        )
    end
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::Type{<:AbstractHydroUnitCommitment},
) where {H <: PSY.HydroGen}
    status_init(psi_container.initial_conditions, devices)
    output_init(psi_container.initial_conditions, devices)
    duration_init(psi_container.initial_conditions, devices)

    return
end

function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::Type{D},
) where {H <: PSY.HydroGen, D <: AbstractHydroDispatchFormulation}
    output_init.initial_conditions_container(psi_container, devices)

    return
end

########################## Addition of to the nodal balances ###############################
function nodal_expression!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    system_formulation::Type{<:PM.AbstractPowerModel},
) where {H <: PSY.HydroGen}
    nodal_expression!(psi_container, devices, PM.AbstractActivePowerModel)
    # Commented out since PF = 1.0 is the assumtion for RoR Hydro
    #=
     parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)
    if parameters
     forecast_label = "get_rating"
     peak_value_function = x -> PSY.get_rating(x) * sin(acos(PSY.get_powerfactor(x)))
    else
     forecast_label = "get_rating"
     peak_value_function = x -> PSY.get_reactivepower(x)
    end
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
     ts_vector = get_time_series(psi_container, d, forecast_label)
     timeseries_data = DeviceTimeSeries(d, peak_value_function, ts_vector)
     constraint_data[ix] = timeseries_data
    end

    if parameters
     include_parameters(
         psi_container,
         constraint_data,
         UpdateRef{R}(REACTIVE_POWER, forecast_label),
         :nodal_balance_active,
     )
     return
    else
     for t in model_time_steps(psi_container)
         for device in constraint_data
             add_to_expression!(
                 psi_container.expressions[:nodal_balance_reactive],
                 device.bus_number,
                 t,
                 device.multiplier * device.timeseries[t],
             )
         end
     end
    end
    return
    =#
end

function nodal_expression!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    system_formulation::Type{<:PM.AbstractActivePowerModel},
) where {H <: PSY.HydroGen}
    parameters = model_has_parameters(psi_container)
    use_forecast_data = model_uses_forecasts(psi_container)
    if use_forecast_data
        forecast_label = "get_rating"
        peak_value_function = x -> PSY.get_rating(x)
    else
        forecast_label = ""
        peak_value_function = x -> PSY.get_activepower(x)
    end
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, peak_value_function, ts_vector)
        constraint_data[ix] = timeseries_data
    end

    if parameters
        include_parameters(
            psi_container,
            constraint_data,
            UpdateRef{H}(ACTIVE_POWER, forecast_label),
            :nodal_balance_active,
        )
        return
    else
        for t in model_time_steps(psi_container)
            for device in constraint_data
                add_to_expression!(
                    psi_container.expressions[:nodal_balance_active],
                    device.bus_number,
                    t,
                    device.multiplier * device.timeseries[t],
                )
            end
        end
    end
    return
end

##################################### Hydro generation cost ############################
function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.HydroEnergyReservoir},
    device_formulation::Type{D},
    system_formulation::Type{<:PM.AbstractPowerModel},
) where {D <: AbstractHydroFormulation}
    add_to_cost(
        psi_container,
        devices,
        variable_name(ACTIVE_POWER, PSY.HydroEnergyReservoir),
        :fixed,
        -1.0,
    )

    return
end

function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::Type{D},
    system_formulation::Type{<:PM.AbstractPowerModel},
) where {D <: AbstractHydroFormulation, H <: PSY.HydroGen}

    return
end

##################################### Water/Energy Limit Constraint ############################
function energy_limit_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, <:AbstractHydroDispatchFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::IntegralLimitFF,
) where {H <: PSY.HydroGen}
    return
end

function energy_limit_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::DeviceModel{H, <:AbstractHydroDispatchFormulation},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroGen}

    forecast_label = "get_storage_capacity"
    constraint_data = Vector{DeviceTimeSeries}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, forecast_label)
        timeseries_data = DeviceTimeSeries(d, x -> PSY.get_rating(x), ts_vector)
        add_device_services!(timeseries_data, d, model)
        constraint_data[ix] = timeseries_data
    end

    if model_has_parameters(psi_container)
        device_energy_limit_param_ub(
            psi_container,
            constraint_data,
            constraint_name(ENERGY_LIMIT, H),
            UpdateRef{H}(ENERGY_BUDGET, forecast_label),
            variable_name(ACTIVE_POWER, H),
        )
    else
        device_energy_limit_ub(
            psi_container,
            constraint_data,
            constraint_name(ENERGY_LIMIT),
            variable_name(ACTIVE_POWER, H),
        )
    end
end

function device_energy_limit_param_ub(
    psi_container::PSIContainer,
    energy_limit_data::Vector{DeviceTimeSeries},
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    variable = get_variable(psi_container, var_name)
    set_name = (r.name for r in energy_limit_data)
    constraint = add_cons_container!(psi_container, cons_name, set_name)
    container = add_param_container!(psi_container, param_reference, set_name, 1)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    for data in energy_limit_data
        name = data.name
        multiplier[name, 1] = data.multiplier
        param[name, 1] = PJ.add_parameter(psi_container.JuMPmodel, sum(data.timeseries))
        constraint[name] = JuMP.@constraint(
            psi_container.JuMPmodel,
            sum([variable[name, t] for t in time_steps]) <= multiplier[name, 1] * param[name, 1]
        )
    end

    return
end

function device_energy_limit_ub(
    psi_container::PSIContainer,
    energy_limit_data::Vector{DeviceTimeSeries},
    cons_name::Symbol,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    variable = get_variable(psi_container, var_name)
    names = (r.name for r in energy_limit_data)
    constraint = add_cons_container!(psi_container, cons_name, names)

    for data in energy_limit_data
        name = data.name
        forecast = data.timeseries
        multiplier = data.multiplier
        constraint[name] = JuMP.@constraint(
            psi_container.JuMPmodel,
            sum([variable[name, t] for t in time_steps]) <= multiplier * sum(forecast)
        )
    end

    return
end


text:
include("MITgcmDatas.jl")
include("days_Parse_Selector.jl")

using DataFrames
typealias Tsize Tuple{Int64, Vararg{Int64}}

function bread(singleFile :: ASCIIString, m :: MITgcmDatas;
               varoffset :: Int64 = 1,
               varsize :: Tsize = (m.nx, m.ny, m.nz) )
  open(singleFile) do fid
    seek(fid,8*(varoffset))
    x = read(fid, Float64, varsize)
    map(ntoh,x)
  end
end

function bread( m :: MITgcmDatas, varfile :: ASCIIString;
                selector :: Function = days ,
                varoffset :: Int64 = 0,
                whether_mean :: Bool = false,
                varsize :: Tsize = (m.nx, m.ny, m.nz) )
  (files, timesteps) = filestimes(m, varfile)
  (fs,ts) = map((files, timesteps)) do x
    x[selector(timesteps, m)]
  end
  println("Reading $(length(ts)) files")

  if (whether_mean)
    result = zeros(varsize...)
    for ii = 1:length(ts)
      result[:,:,:] +=  bread("$(m.Dir)/$(fs[ii])",m, varoffset = varoffset, varsize = varsize)
      println("Filename is $(m.Dir)/$(fs[ii]), with offsize as $varoffset")
    end
    return result/length(ts)
  else
    result = zeros(varsize..., length(ts))
    for ii = 1:length(ts)
      result[:,:,:,ii] = bread("$(m.Dir)/$(fs[ii])",m, varoffset = varoffset, varsize = varsize)
    end
    return result
  end

end

function bread( m :: MITgcmDatas,
                FileDataOffset :: DataFrame, varname :: ASCIIString;
                selector :: Function =  exactTimeStep,
                whether_mean :: Bool = true
                )
  df = FileDataOffset[FileDataOffset[:varname] .==  varname, :]
  @assert (size(df,1) == 1)
  varoffset = df[1, :varoffset]
  varsize = df[1, :varsize]
  varfile = df[1, :varfile]
  return bread( m , varfile ; selector = selector ,
                varoffset = varoffset,
                whether_mean = whether_mean,
                varsize = varsize)
end

export bread


text:
# using GeometryTypes
# using GLVisualize, Colors, GLWindow
# w = glscreen()
# @async renderloop(w)
#
#
# """
# Returns a tuple of points and indices which can be rendered as line segments.
# You can pass an optional `step` plus `stepscale`, to scale the linesegment at
# every `step` by `stepscale`.
# Example output for `step=3`:
# ```
# .              .
# |    .    .    |
# |____|____|____|
# ```
# """
# function tick_axis{N,T}(a::Point{N,T}, b::Point{N,T}, tick_dir::Vec{N,T}, step, step_long=1, step_long_scale=1.)
#     indices = Cuint[]
#     points  = Point{N,T}[]
#     tick_axis!(points, indices, a, b, tick_dir, step, step_long, step_long_scale)
# end
#
# """
# inplace version of `tick_axis`. Improves speed for e.g. animations
# """
# function tick_axis!{N,T}(
#         points, indices,
#         a::Point{N,T}, b::Point{N,T}, tick_dir::Vec{N,T}, step, step_long=1, step_long_scale=1.
#     )
#     ab = Vec{N,T}(b-a)
#     len = norm(ab)
#
#     dir = (ab/len)*step
#     # generate grid point edges like:
#     # :    :    :    : (.) <- optional ending point
#     step_count_high = ceil(Int, len / step)+1
#     step_count_low = floor(Int, len / step)+1
#
#     endpoint = step_count_high-step_count_low #0 or 1
#     if step_count_high*2 > length(points) # resize if necessary
#         resize!(points, step_count_low*2+endpoint)
#         resize!(indices, step_count_low*4+(endpoint*2)) # we can assume, that indices are wrong as well
#     end
#     points[end] = b # insert end point
#     if endpoint==0
#         indices[end], indices[end-1] = length(points)-2, length(points)-1# 0 indexed endpoint
#     end
#     current_position, k = a, 0
#     for i=1:2:(step_count_low*2)
#         current_position = Vec{N,T}(a)+(dir*k)
#         points[i] = current_position
#         scaling = (k%step_long==0) ? step_long_scale : 1.
#         points[i+1] = current_position+(tick_dir*scaling)
#         k+=1
#     end
#     # connect the edges
#     # .    .    .    .
#     #|____|____|____|
#     i, base_index = 1, 0
#     for count=1:(step_count_high-1)
#         # connect tick_dir
#         indices[i  ] = base_index
#         indices[i+1] = base_index+1
#         # connect to next segment
#         indices[i+2] = base_index
#         indices[i+3] = base_index+2
#         i+=4;base_index+=2
#     end
#     points, indices
# end
#
# function annoteded_axis(a, b, tick_dir, richtext, dim)
#     s = 0.1
#     step_long=5; step_long_scale=2f0
#     points, indices = tick_axis(a, b, tick_dir, s, step_long, step_long_scale)
#     i = 2
#     while i <= length(points)
#         p = points[i]
#         ppadd = p + tick_dir*step_long_scale*1.5f0
#         x = string(p[dim])
#         pos = GLVisualize.calc_position(
#             x, Point2f0(0), Vec2f0(0.003),
#             richtext.defaults[:font], richtext.text.atlas
#         )
#         pos3d = map(x-> ppadd+Point3f0(x, 0), pos)
#         insert!(richtext, x, richtext.cursor, Dict(
#             :position => pos3d,
#             :scale => Vec2f0(0.003),
#             :color => RGBA{Float32}(0,0,0,1)
#         ))
#         richtext.cursor += length(x)
#         i += 2*step_long
#     end
#     _view(visualize(points, :linesegment, indices=indices), camera=:perspective)
# end
# function grid_axis(cube, richtext)
#     o = origin(cube)
#     w = widths(cube)
#     c_axis = ntuple(i->unit(Point3f0, i), 3)
#     axis = ntuple(i->w[i].*c_axis[i], 3)
#     xpoint = o + axis[1]
#     ypoint = o + axis[2]
#     xypoint = o + axis[1] + axis[2]
#     origins = Point3f0[xpoint, xpoint, xypoint]
#     targets = Point3f0[xpoint+axis[3], xypoint, o+axis[2]]
#     tickdirs = map([c_axis[1], c_axis[1], c_axis[2]]) do x
#         Vec3f0(0.03*x)
#     end
#     for i=1:3
#         annoteded_axis(origins[i], targets[i], tickdirs[i], richtext, i)
#     end
#     _view(visualize(cube, :grid))
# end
# bb = AABB{Float32}(Vec3f0(-0.4, 0.5, -0.7), Vec3f0(2, 3, 4))
# N = 1000
# _view(visualize(
#     randstring(N),
#     scale=fill(Vec2f0(0.1), N),
#     color=fill(RGBA{Float32}(0,0,0,0), N),
#     rotation=fill(Vec3f0(0,0,1), N),
#     position=fill(Point3f0(0), N),
#     billboard=true
# ), camera=:perspective)
#
#
# richtext = GLVisualize.RichText(GLVisualize.Text(renderlist(w)[1]))
#
# grid_axis(bb, richtext)
# _view(visualize(GLNormalMesh(Sphere(Point3f0(0.5,2,2), 0.2f0))))
#
#
# # richtext.cursor
# #
# # insert!(richtext, x, 0:1, Dict(
# #     :position => Point3f0[Point3f0(0,i*0.8,0) for i=0:(length(x)-1)],
# #     :scale => Vec2f0(0.05),
# #     :rotation => Vec3f0(pi,0,4*pi)
# # ))
#
# # #
# # #
# # # immutable Glyph{N,T,C<:Colorant}
# # #     position::Point{N,T}
# # #     scale::Vec{2,T}
# # #     uv_width::Vec{4,T}
# # #     offset::Vec{2,T}
# # #     rotation::Vec{N,T}
# # #     color::C
# # # end
# # #
# # #
# # # type TextBuffer
# # #     text::ArrayOfStructs{Glyph{N,T}}
# # #     atlas
# # # end
# # # length(t::TextBuffer) = length(t.text)
# # # setindex!(t::TextBuffer, value::Glyph, index) = t.text[index] = value
# # # function setindex!(t::TextBuffer, value::Tuple, index)
# # #     char, font, position, scale = value
# # #     t.text[index] = Glyph(
# # #         position,
# # #         scale,
# # #         get_uv_offsetwidth!(t.atlas, char, font),
# # #         get_scale!(t.atlas, char, font, scale),
# # #         glyph_bearing!(t.atlas, char, font, scale)
# # #     )
# # # end
# # # function splice!(text::TextBuffer, to_insert::TextTypes, start_index=1)
# # #     len = offset+length(to_insert)
# # #     if len <= length(text)
# # #         resize!(text, len)
# # #     end
# # #     i = start_index
# # #     for char in to_insert
# # #         text[i] = (char, font, position, scale)
# # #         position = calc_position(position, start_pos, atlas, char, font, scale)
# # #         i+=1
# # #     end
# # #     positions, uvs, scales
# # # end
# # #
# # # function tick_labels(t_start::Vec{N,T}, t_end::Vec{N,T}, tick_dir::Vec{N,T}, step::T, long_step=5, precision)
# # #     # do text only for long step
# # #     step_count = floor(Int, (b-a) / step / long_step)
# # #     positions = Point{N,T}[]
# # #     sizehint!(positions, step_count*precision) # should be roughly that many digitis, but can't say for sure
# # #     tick = a
# # #     for k=1:step_count
# # #         tick += step*long_step
# # #         str = Showoff.format_fixed_scientific(tick, precision, false)
# # #         for char in str
# # #
# # #
# # #
# # #
# # # end
# # #
# # # points, indices = gen_grid(Vec3f0(0,0,1), Vec3f0(0,0.1,0), linspace(0f0,3f0, 20), 5, 1.5)
# # # empty!(w.renderlist)
# # # _view(visualize(points, :linesegment, indices=indices), camera=:perspective)
# # # map(Int, indices)


text:
module GaussianProcesses_

export GPClassifier

import MLJBase

using CategoricalArrays

import ..GaussianProcesses # strange lazy-loading syntax

const GP = GaussianProcesses

# here T is target type:
const CD{T,C} = MLJBase.CategoricalDecoder{Int,false,T,1,UInt32,C}
const GPClassifierFitResultType{T} =
    Tuple{GP.GPE,     # TODO: make this a concrete type for ensembling efficiency
          Union{CD{T,CategoricalValue{T,UInt32}},
                CD{T,CategoricalString{UInt32}}}}

mutable struct GPClassifier{T, M<:GP.Mean, K<:GP.Kernel} <: MLJBase.Deterministic{GPClassifierFitResultType{T}}
    target_type::Type{T} # target is CategoricalArray{target_type}
    mean::M
    kernel::K
end

function GPClassifier(
    ; target_type=Int
    , mean=GP.MeanZero()
    , kernel=GP.SE(0.0,1.0)) # binary

    model = GPClassifier(
        target_type
        , mean
        , kernel)

    message = MLJBase.clean!(model)
    isempty(message) || @warn message

    return model
end

# function MLJBase.clean!

function MLJBase.fit(model::GPClassifier{T2,M,K}
            , verbosity::Int
            , X
            , y::CategoricalVector{T}) where {T,T2,M,K}

    Xmatrix = MLJBase.matrix(X)

    T == T2 || throw(ErrorException("Type, $T, of target incompatible "*
                                    "with type, $T2, of $model."))

    decoder = MLJBase.CategoricalDecoder(y, Int)
    y_plain = MLJBase.transform(decoder, y)

    if VERSION < v"1.0"
        XT = collect(transpose(Xmatrix))
        yP = convert(Vector{Float64}, y_plain)
        gp = GP.GPE(XT
                  , yP
                  , model.mean
                  , model.kernel)

        GP.fit!(gp, XT, yP)
    else
        gp = GP.GPE(transpose(Xmatrix)
                  , y_plain
                  , model.mean
                  , model.kernel)
        GP.fit!(gp, transpose(Xmatrix), y_plain)
    end

    fitresult = (gp, decoder)

    cache = nothing
    report = nothing

    return fitresult, cache, report
end

function MLJBase.predict(model::GPClassifier{T}
                       , fitresult
                       , Xnew) where T

    Xmatrix = MLJBase.matrix(Xnew)

    gp, decoder = fitresult

    nlevels = length(decoder.pool.levels)
    pred = GP.predict_y(gp, transpose(Xmatrix))[1] # Float
    # rounding with clamping between 1 and nlevels
    pred_rc = clamp.(round.(Int, pred), 1, nlevels)

    return MLJBase.inverse_transform(decoder, pred_rc)
end

# metadata:
MLJBase.load_path(::Type{<:GPClassifier}) = "MLJModels.GaussianProcesses_.GPClassifier" # lazy-loaded from MLJ
MLJBase.package_name(::Type{<:GPClassifier}) = "GaussianProcesses"
MLJBase.package_uuid(::Type{<:GPClassifier}) = "891a1506-143c-57d2-908e-e1f8e92e6de9"
MLJBase.package_url(::Type{<:GPClassifier}) = "https://github.com/STOR-i/GaussianProcesses.jl"
MLJBase.is_pure_julia(::Type{<:GPClassifier}) = true
MLJBase.input_scitypes(::Type{<:GPClassifier}) = MLJBase.Continuous
MLJBase.target_scitype(::Type{<:GPClassifier}) = MLJBase.FiniteOrderedFactor
MLJBase.input_is_multivariate(::Type{<:GPClassifier}) = true

end # module


text:
using DashBootstrapComponents

button_group = dbc_buttongroup(
    [
        dbc_button("First"),
        dbc_button("Second"),
        dbc_dropdownmenu(
            [dbc_dropdownmenuitem("Item 1"), dbc_dropdownmenuitem("Item 2")],
            label="Dropdown",
            group=true,
        ),
    ],
    vertical=true,
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module cudadevrt

using CUDA, BenchmarkTools, Random

const threads = 256
#simple add matrix and vector kernel
function kernel_add_mat_vec(m, x1, x2, y)
    # one block per column
    offset = (blockIdx().x-1) * m
    @inbounds xtmp = x2[blockIdx().x]
    for i = threadIdx().x : blockDim().x : m
        @inbounds y[offset + i] = x1[offset + i] + xtmp
    end
    return
end

function add!(y, x1, x2)
    m, n = size(x1)
    @cuda blocks = n, 1 threads = threads kernel_add_mat_vec(m, x1, x2, y)
end

function main()
    Random.seed!(1)
    m, n = 3072, 1536    # 256 multiplier
    x1 = cu(randn(Float32, (m, n)) .+ Float32(0.5))
    x2 = cu(randn(Float32, (1, n)) .+ Float32(0.5))
    y1 = similar(x1)

    results = @benchmark CUDA.@sync add!($y1, $x1, $x2)

    # BenchmarkTools captures inputs, JuliaCI/BenchmarkTools.jl#127, so forcibly free them
    CUDA.unsafe_free!(x1)
    CUDA.unsafe_free!(x2)
    CUDA.unsafe_free!(y1)

    return results
end

end

cudadevrt.main()


text:
# All tests for 3D parallel PUMI meshes go here

# MPI Declarations
comm = MPI.COMM_WORLD
comm_world = MPI.MPI_COMM_WORLD
comm_self = MPI.COMM_SELF
my_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

opts = Dict{ASCIIString, Any}()
opts["order"] = 1
opts["dimensions"] = 3
opts["use_DG"] = true
opts["Tsbp"] = Float64
opts["Tmsh"] = Complex128
opts["operator_type"] = "SBPOmega"
opts["smb_name"] = "../src/mesh_files/tet8cube_np2.smb"
opts["dmg_name"] = ".null"

opts["numBC"] = 3
opts["BC1"] = [ 0, 1, 2, 3]
opts["BC1_name"] = "ExpBC"
opts["BC2"] = [4]
opts["BC2_name"] = "noPenetrationBC"
opts["BC3"] = [10]
opts["BC3_name"] = "noPenetrationBC"

opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
opts["jac_type"] = 2
opts["jac_method"] = 2
opts["run_type"] = 5

Tsbp = opts["Tsbp"]
Tmsh = opts["Tmsh"]
shape_type = 2
sbp = getTetSBPOmega(degree=opts["order"], Tsbp=Tsbp)
ref_verts = sbp.vtx
face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
topo = ElementTopology{3}(face_verts)
sbpface = TetFace{Tsbp}(opts["order"], sbp.cub, ref_verts)
mesh = PumiMeshDG3{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"], sbp,
                         opts, sbpface, topo; dofpernode=1,
                         coloring_distance=opts["coloring_distance"],
                         shape_type=shape_type)

orig_vert_coords = deepcopy(mesh.vert_coords)

@testset "--- Check if Unique wall coordinates are being computed correctly across all ranks ---" begin

  geom_faces = opts["BC2"]
  wallCoords = FFD.getGlobalUniqueWallCorrdsArray(mesh, geom_faces)
  if my_rank == 0
    @test norm(wallCoords -[0.0 0.0 0.0 0.0 0.0 0.0
                            0.0 0.5 0.0 0.5 1.0 1.0
                            0.5 1.0 1.0 0.5 1.0 0.5] ) < 1e-14
  end

  if my_rank == 1
    @test norm(wallCoords -[0.0 0.0 0.0
                            0.0 0.5 1.0
                            0.0 0.0 0.0]) < 1e-14
  end

end # End facts
MPI.Barrier(comm)

@testset "--- Checking FFD on 3D parallel DG Pumi meshes ---" begin

  @testset "Check control point manipulation with nonlinear mapping on full mesh" begin

    ndim = mesh.dim
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]
    offset = [0., 0., 0.5] # No offset in the X & Y direction
    full_geom = true

    # Create a mapping object using nonlinear mapping
    map, box = initializeFFD(mesh, sbp, order, nControlPts, offset, full_geom)
    @fact norm(box.box_bounds - [0.0 0.0 -0.5
                                 1.0 1.0 1.5]) < 1e-14

    # Rigid body rotation
    theta = -20*pi/180  # Rotate wall coordinates by 20 degrees
    rotMat = [cos(theta) -sin(theta) 0
              sin(theta) cos(theta)  0
              0          0           1] # Rotation matrix
    # Rotate the control points
    for k = 1:map.nctl[3]
      for j = 1:map.nctl[2]
        for i = 1:map.nctl[1]
          map.cp_xyz[:,i,j,k] = rotMat*map.cp_xyz[:,i,j,k]
        end
      end
    end

    # Rigid body translation
    map.cp_xyz[1,:,:,:] += 0.2
    map.cp_xyz[2,:,:,:] += 0.3
    map.cp_xyz[3,:,:,:] += 0.5

    vertices = evalVolume(map, mesh)
    commitToPumi(map, mesh, sbp, vertices, opts)

    fname = string("./testvalues/translation_plus_rotation_DG_3D_tet8cube_rank", my_rank,".dat")
    # f = open(fname, "w")
    # for i = 1:length(mesh.vert_coords)
    #   println(f, mesh.vert_coords[i])
    # end
    # close(f)

    test_values = readdlm(fname)
    for i = 1:length(test_values)
      err = abs.(test_values[i] - mesh.vert_coords[i])
      @test  err  < 1e-14
    end

  end # End context("Check control point manipulation with nonlinear mapping on full mesh")

  MPI.Barrier(comm)
  # Reset the coordinates and mesh to the original value
  for i = 1:mesh.numEl
    update_coords(mesh, i, real(orig_vert_coords[:,:,i]))
  end
  commit_coords(mesh, sbp, opts)

  @testset "Check control point manipulation on a geometry face" begin

    ndim = mesh.dim
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]
    offset = [0.5, 0.5, 0.5] # No offset in the X & Y direction
    geom_faces = opts["BC2"]

    # Create a mapping object using nonlinear mapping
    map, box = initializeFFD(mesh, sbp, order, nControlPts, offset, false, geom_faces)

    # Rigid body rotation
    theta = -2*pi/180  # Rotate wall coordinates by 2 degrees
    rotMat = [cos(theta) -sin(theta) 0
              sin(theta) cos(theta)  0
              0          0           1] # Rotation matrix
    # Rotate the control points
    for k = 1:map.nctl[3]
      for j = 1:map.nctl[2]
        for i = 1:map.nctl[1]
          map.cp_xyz[:,i,j,k] = rotMat*map.cp_xyz[:,i,j,k]
        end
      end
    end

    # Rigid body translation
    map.cp_xyz[1,:,:,:] += 0.02
    map.cp_xyz[2,:,:,:] += 0.03
    map.cp_xyz[3,:,:,:] += 0.05

    vertices = evalSurface(map, mesh)
    commitToPumi(map, mesh, sbp, vertices, opts)

    fname = string("./testvalues/translation_plus_rotation_DG_3D_tet8cube_face4_rank", my_rank, ".dat")
    # f = open(fname, "w")
    # for i = 1:length(mesh.vert_coords)
    #   println(f, mesh.vert_coords[i])
    # end
    # close(f)

    test_values = readdlm(fname)
    for i = 1:length(test_values)
      err = abs.(test_values[i] - mesh.vert_coords[i])
      @test  err  < 1e-14
    end

  end # End context("Check control point manipulation on a geometry face")

  MPI.Barrier(comm)
  # Reset the coordinates and mesh to the original value
  for i = 1:mesh.numEl
    update_coords(mesh, i, real(orig_vert_coords[:,:,i]))
  end
  commit_coords(mesh, sbp, opts)

  @testset "--- Checking evaldXdControlPointProduct for 3D DG Mesh ---" begin

    ndim = mesh.dim
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]
    offset = [0.5, 0.5, 0.5] # No offset in the X & Y direction
    geom_faces = opts["BC2"]

    # Create a mapping object using nonlinear mapping
    map, box = initializeFFD(mesh, sbp, order, nControlPts, offset, false, geom_faces)

    fill!(map.work, 0.0)

    # Create seed vector
    # - Get original wall coordinates
    orig_wallCoords = FFD.getGlobalUniqueWallCorrdsArray(mesh, geom_faces)
    Xs_bar = ones(3, size(orig_wallCoords,2))
    evaldXdControlPointProduct(map, mesh, vec(Xs_bar))

    fname = "./testvalues/evaldXdControlPointProduct_tet8cube.dat"

    test_values = readdlm(fname)
    @test ( length(map.work) )== length(test_values)
    for i = 1:length(map.work)
      err = abs.(test_values[i] - map.work[i])
      @test  err  < 1e-14
    end

  end # End context("--- Checking evaldXdControlPointProduct for 2D DG Mesh ---")

end # End facts("--- Checking FFD on 3D serial DG Pumi meshes ---")


text:
export AbstractLandmarkSelection,
       LandmarkBasedRepresentation,
       RandomLandmarkSelection,
       EvenlySpacedLandmarkSelection,
       select_landmarks,
       BresenhamLandmarkSelection,
       MS3
using StatsBase
"""
```julia
abstract type AbstractLandmarkSelection end
```
Abstract type that defines how to sample data points. Types that inherint from `AbstractLandmarkSelection` has to implements the following interface:

```julia
select_landmarks{L<:AbstractLandmarkSelection}(c::L, X)
```
The `select_landmarks`function returns an array with the indices of the sampled points.
# Arguments

- ```c::T<:AbstractLandmarkSelecion```. The landmark selection type.
- ```d::D<:DataAccessor```.  The [`DataAccessor`](ref) type.
- ```X```. The data. The data to be sampled.
"""
abstract type AbstractLandmarkSelection end

"""
```
struct RandomLandmarkSelection <: LandmarkSelection
```
`Random` random samples `n` data points from a dataset.

"""

struct RandomLandmarkSelection <: AbstractLandmarkSelection
end
"""
```julia
select_landmarks(c::RandomLandmarkSelection,d::T,n::Integer, X)
```
The function returns `n`random points according to `RandomLandmarkSelection`
# Arguments
- c::RandomLandmarkSelection.
- n::Integer. The number of data points to sample.
- X. The data to be sampled.
"""
function select_landmarks(c::RandomLandmarkSelection,n::Integer, X)
  return  StatsBase.sample(1:number_of_patterns(X), n, replace = false)
end

"""
```
struct EvenlySpacedLandmarkSelection <: AbstractLandmarkSelection
```
The `EvenlySpacedLandmarkSelection` selection method selects  `n` evenly spaced points  from a dataset.

"""
struct EvenlySpacedLandmarkSelection  <: AbstractLandmarkSelection
end
"""
```
select_landmarks(c::EvenlySpacedLandmarkSelection,n::Integer, X)
```

"""
function select_landmarks(c::EvenlySpacedLandmarkSelection, n::Integer, X)
    m = number_of_patterns(X)
    return collect(1:round(Integer,floor(m/n)):m)[1:n]
end

struct BresenhamLandmarkSelection  <: AbstractLandmarkSelection
end
function select_landmarks(c::BresenhamLandmarkSelection, m::Integer, X)
    n = number_of_patterns(X)
    return round.(Integer, [(i-1)*n//m + n//(2*m) for i=1:m])
end

"""
```
struct MS3 <: AbstractLandmarkSelection
    proportion::Float64
    sim::Function
end
```
The `MS3` selection method selects  `m`
NYSTROM SAMPLING DEPENDS ON THE EIGENSPECTRUM SHAPE OF THE DATA
"""

struct MS3 <: AbstractLandmarkSelection
    proportion::Float64
    sim::Function
end
function select_landmarks(c::MS3, m::Integer, X)
   cant = number_of_patterns(X)
   points  = [rand(1:cant); rand(1:cant)]
    while (length(points)<m)
       T_candidates  = setdiff(1:cant, points)
       T = rand(T_candidates, round(Integer,length(T_candidates)*c.proportion))
       min_simmilarity = Inf
       min_point = 0
        for t in T
           simmilarity = 0
            for p in points
                simmilarity = simmilarity + c.sim(get_element(X,p), get_element(X,t))^2
            end
            if simmilarity < min_simmilarity
                min_simmilarity = simmilarity
                min_point = t
            end
        end
        push!(points, min_point)
    end
    return points
end