text:
const CEPS = Ptr{Cvoid}

struct SlepcEPS
    ptr::Ref{CEPS}

    SlepcEPS() = new(Ref{Ptr{Cvoid}}())
end

# allows us to pass SlepcEPS objects directly into CEPS ccall signatures
Base.cconvert(::Type{CEPS}, eps::SlepcEPS) = eps.ptr[]

"""
    EPSCreate(comm::MPI.Comm, eps::SlepcEPS)

Wrapper for EPSCreate
"""
function EPSCreate(comm::MPI.Comm, eps::SlepcEPS)
    error = ccall((:EPSCreate, libslepc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CEPS}), comm, eps.ptr)
    @assert iszero(error)
end

function EPSCreate(comm::MPI.Comm = MPI.COMM_WORLD)
    eps = SlepcEPS()
    EPSCreate(comm, eps)
    return eps
end

"""
    EPSSetOperators(eps::SlepcEPS, A::PetscMat, B::PetscMat)

Wrapper for EPSSetOperators with two matrices
"""
function EPSSetOperators(eps::SlepcEPS, A::PetscMat, B::PetscMat)
    error = ccall((:EPSSetOperators, libslepc), PetscErrorCode, (CEPS, CMat, CMat), eps, A, B)
    @assert iszero(error)
end

"""
    EPSSetOperators(eps::SlepcEPS, A::PetscMat)

Wrapper for EPSSetOperators with only one matrix
"""
function EPSSetOperators(eps::SlepcEPS, A::PetscMat)
    error = ccall((:EPSSetOperators, libslepc), PetscErrorCode, (CEPS, CMat, CMat), eps, A, C_NULL)
    @assert iszero(error)
end

"""
    EPSSetUp(eps::SlepcEPS)

Wrapper for EPSSetup
"""
function EPSSetUp(eps::SlepcEPS)
    error = ccall((:EPSSetUp, libslepc), PetscErrorCode, (CEPS,), eps)
    @assert iszero(error)
end

"""
    EPSSetFromOptions(eps::SlepcEPS)

Wrapper for EPSSetFromOptions
"""
function EPSSetFromOptions(eps::SlepcEPS)
    error = ccall((:EPSSetFromOptions, libslepc), PetscErrorCode, (CEPS,), eps)
    @assert iszero(error)
end

"""
    EPSGetOperators(eps::SlepcEPS)

Wrapper for EPSGetOperators
"""
function EPSGetOperators(eps::SlepcEPS)
    A = PetscMat()
    B = PetscMat()

    error = ccall((:EPSGetOperators, libslepc), PetscErrorCode, (CEPS, Ptr{CMat}, Ptr{CMat}), eps, A.ptr, B.ptr)
    @assert iszero(error)

    return A, B
end

"""
    EPSSolve(eps::SlepcEPS)

Wrapper for EPSSolve
"""
function EPSSolve(eps::SlepcEPS)
    error = ccall((:EPSSolve, libslepc), PetscErrorCode, (CEPS,), eps)
    @assert iszero(error)
end

"""
    EPSGetConverged(eps::SlepcEPS)

Wrapper for `EPSGetConverged` : return the number of converged eigenvalues.
"""
function EPSGetConverged(eps::SlepcEPS)
    nconv = Ref{PetscInt}(0)

    error = ccall((:EPSGetConverged, libslepc), PetscErrorCode, (CEPS, Ref{PetscInt}), eps, nconv)
    @assert iszero(error)

    return nconv[]
end

"""
    EPSGetEigenvalue(eps::SlepcEPS, ieig)

Wrapper for EPSGetEigenvalue . SLEPc 0-based indexing is used : `0 < ieig < EPSGetConverged-1` ()

A tuple (real, imag) is returned.
"""
function EPSGetEigenvalue(eps::SlepcEPS, ieig)
    eigr = Ref{PetscScalar}()
    eigi = Ref{PetscScalar}()

    error = ccall((:EPSGetEigenvalue, libslepc),
                   PetscErrorCode,
                   (CEPS, PetscInt, Ref{PetscScalar}, Ref{PetscScalar}),
                   eps, PetscInt(ieig), eigr, eigi
    )
    @assert iszero(error)

    return eigr[], eigi[]
end

"""
    EPSGetTolerances(eps::SlepcEPS)

Wrapper for `EPSGetTolerances`
"""
function EPSGetTolerances(eps::SlepcEPS)
    tol = Ref{PetscReal}(0)
    maxits = Ref{PetscInt}(0)

    error = ccall((:EPSGetTolerances, libslepc), PetscErrorCode, (CEPS, Ref{PetscReal}, Ref{PetscInt}), eps, tol, maxits)
    @assert iszero(error)

    return tol[], maxits[]
end

"""
    EPSGetEigenpair(eps::SlepcEPS, ieig, vecr::PetscVec, veci::PetscVec)

Wrapper for `EPSGetEigenpair`. SLEPc 0-based indexing is used : `0 < ieig < EPSGetConverged-1`
"""
function EPSGetEigenpair(eps::SlepcEPS, ieig, vecr::PetscVec, veci::PetscVec)
    eigr = Ref{PetscScalar}(0.)
    eigi = Ref{PetscScalar}(0.)

    error = ccall((:EPSGetEigenpair, libslepc),
                    PetscErrorCode,
                    (CEPS, PetscInt, Ref{PetscScalar}, Ref{PetscScalar}, CVec, CVec),
                    eps, PetscInt(ieig), eigr, eigi, vecr, veci
    )
    @assert iszero(error)

    return eigr[], eigi[], vecr, veci
end

"""
    EPSGetEigenpair(eps::SlepcEPS, mat::PetscMat, ieig)

Wrapper for EPSGetEigenpair without providing pre-allocated vec. but providing one of the operator matrix
"""
function EPSGetEigenpair(eps::SlepcEPS, mat::PetscMat, ieig)
    eigr = Ref{PetscScalar}()
    eigi = Ref{PetscScalar}()
    vecr, veci = MatCreateVecs(mat)

    return EPSGetEigenpair(eps, ieig, vecr, veci)
end
"""
    EPSGetEigenpair(eps::SlepcEPS, ieig)

Wrapper for EPSGetEigenpair without providing pre-allocated vec. nor one of the operator matrix
"""
function EPSGetEigenpair(eps::SlepcEPS, ieig)
    A, B = EPSGetOperators(eps)
    vecr, veci = MatCreateVecs(A)

    return EPSGetEigenpair(eps, ieig, vecr, veci)
end

"""
    EPSView(eps::SlepcEPS, view::PetscViewer = C_NULL)

Wrapper for EPSView
"""
function EPSView(eps::SlepcEPS, viewer::PetscViewer = C_NULL)
    error = ccall((:EPSView, libslepc), PetscErrorCode, (CEPS, PetscViewer), eps, viewer)
    @assert iszero(error)
end

"""
    EPSDestroy(eps::SlepcEPS)

Wrapper for EPSDestroy
"""
function EPSDestroy(eps::SlepcEPS)
    error = ccall((:EPSDestroy, libslepc), PetscErrorCode, (Ptr{CEPS},), eps.ptr)
    @assert iszero(error)
end


text:
#=
The edit distance between two strings refers to the minimum number of character insertions, deletions, and substitutions required to change one string to the other. For example, the edit distance between “kitten” and “sitting” is three: substitute the “k” for “s”, substitute the “e” for “i”, and append a “g”.

Given two strings, compute the edit distance between them.
=#

using Test
include("Solutions/problem31_edit_distance.jl")

@test editdistance("kitten", "sitting") == 3
@test editdistance("cat", "dog") == 3

ext: jl

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
#=##############################################################################
# DESCRIPTION
    Vehicle type handling all defined geometries and their properties.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

################################################################################
# ABSTRACT VEHICLE TYPE
################################################################################
"""
    `AbstractVehicle{N, M, R}`

Type handling all geometries and subsystems that define a flight vehicle.

`N` indicates the number of tilting systems in this vehicle, while and `M`
indicates the number of rotor systems. `R` is a Real type.

Implementations must have the following properties:
    * `V::Array{R, 1}`          : Current velocity vector of the vehicle in the
                                    global coordinate system.
    * `W::Array{R, 1}`          : Current angular velocity vector of the vehicle
                                    in the global coordinate system.
"""
abstract type AbstractVehicle{N, M, R} end

##### FUNCTIONS REQUIRED IN IMPLEMENTATIONS ####################################
"""
    `get_ntltsys(self::AbstractVehicle)`
Return number of tilting systems.
"""
function get_ntltsys(self::AbstractVehicle)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `get_nrtrsys(self::AbstractVehicle)`
Return number of rotor systems.
"""
function get_nrtrsys(self::AbstractVehicle)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `add_dV(self::AbstractVehicle, dV)`
Increments the velocity of the vehicle by `dV`.
"""
function add_dV(self::AbstractVehicle, dV)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `add_dW(self::AbstractVehicle, dV)`
Increments the angular velocity of the vehicle by `dW`.
"""
function add_dW(self::AbstractVehicle, dW)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `set_V(self::AbstractVehicle, V)`
Set current vehicle velocity to `V`.
"""
function set_V(self::AbstractVehicle, V)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `set_W(self::AbstractVehicle, W)`
Set current vehicle angular velocity to `W`.
"""
function set_W(self::AbstractVehicle, W)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `nextstep_kinematic(self::AbstractVehicle, angles)`
Tilts every tilting system of this vehicle into its corresponding new angle,
where `angles[i]=[Ax, Ay, Az]` is the new angle of the i-th tilting system (in
degrees).
"""
function tilt_systems(self::AbstractVehicle, angles)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `nextstep_kinematic(self::AbstractVehicle, dt::Real)`
Translates and rotates the vehicle in a time step `dt` according to current
linear and angular velocity.
"""
function nextstep_kinematic(self::AbstractVehicle, dt::Real)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `precalculations(self, pfield, t, dt)`
Precalculations before calling the solver and before shedding trialing wake.
"""
function precalculations(self::AbstractVehicle, Vinf::Function,
                                pfield::vpm.ParticleField, t::Real, dt::Real)
    error("$(typeof(self)) has no implementation yet!")
end

"""
Shed VPM wake. If `unsteady_shedcrit<=0` it sheds wake due to loading
distribution, otherwise it will shed wake due to unsteady loading.
"""
function shed_wake(self::AbstractVehicle, Vinf::Function,
                            pfield::vpm.ParticleField, dt::Real; t=0.0,
                            unsteady_shedcrit=-1.0, p_per_step=1,
                            sigmafactor=1.0, overwrite_sigma=nothing)
    error("$(typeof(self)) has no implementation yet!")
end

# """
# Solves the aerodynamics of this vehicle.
# """
# function solve(self::AbstractVehicle, Vinf::Function, pfield::vpm.ParticleField,
#                 wake_coupled::Bool, vpm_solver::String,  t::Real, dt::Real,
#                                                       rlx::Real, sigma::Real)
#     error("$(typeof(self)) has no implementation yet!")
# end

"""
Returns a function that generates an array of particles representing the surface
of the vehicle as a collection of vortex particles.
"""
function generate_static_particle_fun(self::AbstractVehicle, sigma::Real)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `save_vtk(self::AbstractVehicle, filename; path=nothing, num=nothing,
optargs...)`

Output VTK files with vehicle geometry and solution fields.
"""
function save_vtk(self::AbstractVehicle, filename;
                        path=nothing, num=nothing, optargs...)
    error("$(typeof(self)) has no implementation yet!")
end

##### COMMON FUNCTIONS  ########################################################
##### COMMON INTERNAL FUNCTIONS  ###############################################
##### END OF ABSTRACT VEHICLE ##################################################

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
struct TrainingSampleResult
    sample::TrainingSample
    outputs::Vector{SingleNetworkOutput}
end

struct TrainingSampleResultSVG
    result::TrainingSampleResult
end

svg(x::TrainingSampleResult) = TrainingSampleResultSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::TrainingSampleResultSVG)
    sample = x.result.sample
    initial_grid = getinferencegrid(
        sample.image_state,
        sample.targets[1],
        "INITIAL INFERENCE",
        "INPUT"
    )
    recurrent_grids = getrecurrentinfgrid(sample)
    g = Grid(
        [initial_grid, recurrent_grids...],
        :,
        1,
        header=(text="TRAINING SAMPLE"),
        line_thickness=0,
        margin=15,
    )
    show(io, m, g)
end


text:
module HelloX
using Example
using UnicodePlots

greet() = print("Hello World!")

function julia_main()
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function moonbar()
    try
        print("\x1b[?25l") # hide cursor
        for i in 1:3
            for 🌝 ∈ '🌑':'🌘'
                sleep(0.1)
                print(🌝)
                print('\r')
            end
        end
    finally
        print("\x1b[?25h") # unhide cursor
    end
end

function real_main()
    @show ARGS
    @show Base.PROGRAM_FILE
    println("Hello World from HelloX.jl")
    println(Example.hello("World from Example.jl"))
    plt = lineplot([cos, sin], -π/2, 2π)
    plt = lineplot!(plt, -0.5, .2, name = "line")
    show(plt)
    println() # newline
    moonbar()
    return
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end # module


text:
using Base.Order
function get_colptr_LUT{T,Ti}(A::SparseMatrixCSC{T,Ti}, i0::Integer, i1::Integer)
    i0 = convert(Ti, i0)
    i1 = convert(Ti, i1)
    if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end
    r1 = Int(A.colptr[i1])
    r2 = Int(A.colptr[i1+1]-1)

    i = (r1 > r2) ? r1 : searchsortedfirst(A.rowval, i0, r1, r2, Forward)

    return i
end

function get_colptrs(K::SparseMatrixCSC, fp::FEProblem)
    #colptrs = zeros(Int, length(K.nzval))
    colptrs = Int[]
    z = 1
    for section in fp.sections
        z = get_colptrs(K, z, colptrs, section, fp.nodes)
    end
    return colptrs
end

function get_colptrs(K::SparseMatrixCSC, z::Int, colptrs::Vector{Int},
                     section::FESection, nodes::Vector{FENode2})
    for element in section.elements
               dof1_n = 0
        for vertex1 in element.vertices
            for dof1 in nodes[vertex1].dofs
                dof1_n += 1
                if dof1.active
                    dof2_n = 0
                    for vertex2 in element.vertices
                        for dof2 in nodes[vertex2].dofs
                            dof2_n += 1
                            if dof2.active
                                push!(colptrs, get_colptr_LUT(K, dof1.eq_n, dof2.eq_n))
                            end
                        end
                    end
                end
            end
        end
    end
    return z
end

#Ke = stiffness(element, nodes, mat)

#=using Base.Order

immutable FEMSParseMatrix{T, Ti}
    K::SparseMatrixCSC{T,Ti}
    eqn2nzval::Vector{Int}
end

@inline get_mat(K::FEMSParseMatrix) = K.K
@inline addval(K::FEMSParseMatrix, eqn::Int, val) = K.K.nzval[K.eqn2nzval[eqn]] += val

function gen_sparse(fp::FEProblem)
    dof_rows = Int[]
    dof_cols = Int[]
    for section in fp.sections
        gen_sparse(section, fp.nodes, dof_rows, dof_cols)
    end

    K = Base.sparse(dof_rows, dof_cols, ones(length(dof_rows)), fp.n_eqs, fp.n_eqs)

    for section in fp.sections
        gen_eq2nzval(section, fp.nodes, dof_rows, dof_cols)
    end

    eqn2nzval =
    for element in section.elements
        for (dof1, i) in activedofs(element, nodes)
            for (dof2, j) in activedofs(element, nodes)
                idx = get_colptr_LUT(K, dof1.eq_n, dof2.eq_n)
                push!(colptrs, g)
                #z += 1
            end
        end
    end
    return z
end

function
    #colptrs = zeros(Int, length(K.nzval))
    colptrs = Int[]
    z = 1
    for section in fp.sections
        z = get_colptrs(K, z, colptrs, section, fp.nodes)
    end
    return colptrs
end

function gen_sparse(section::FESection, nodes::Vector{FENode2},
                   dof_rows::Vector{Int}, dof_cols::Vector{Int})
    for element in section.elements
        for (dof1, _) in activedofs(element, nodes)
            for (dof2, _) in activedofs(element, nodes)
                push!(dof_rows, dof1.eq_n)
                push!(dof_cols, dof2.eq_n)
            end
        end
    end
end

function get_nzval{T,Ti}(A::SparseMatrixCSC{T,Ti}, i0::Integer, i1::Integer)
    i0 = convert(Ti, i0)
    i1 = convert(Ti, i1)
    if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end
    r1 = Int(A.colptr[i1])
    r2 = Int(A.colptr[i1+1]-1)

    i = (r1 > r2) ? r1 : searchsortedfirst(A.rowval, i0, r1, r2, Forward)

    return i
end

=#


text:
# Geographic phases

"""
    PhaseGeog(model, name, evlon, evlat, depth, time, dtdd, inc, takeoff)
    PhaseGeog(..., pierce, lon, lat, radius)

Construct a `Phase3Geog, which represents a single event-station path and one single phase
arrival between two geographic points.

`PhaseGeog` objects have the following fields:

- `model`: Name of model used to calculate phase properties
- `name`: Name of phase
- `evlon`, `evlat`, `depth`: Event longitude and latitude (in °) and depth (in km)
- `stlon`, `stlat`: Station longitude and latitude in °
- `delta`: Epicentral distance in °
- `time`: Travel time in s
- `dtdd`: Horizontal slowness (ray parameter) in s/°
- `inc`: Incidence angle at receiver, measured from downwards in °
- `takeoff`: Takeoff angle at source, measured from downwards in °

The following fields are filled if the ray path has been calculated via `path`:

- `lon`, `lat: `Vector`s of geographic coodinates in ° along the path
- `radius`: `Vector` of radii in km along the path
"""
struct PhaseGeog{T} <: AbstractPhase
    model::String
    name::String
    evlon::T
    evlat::T
    stlon::T
    stlat::T
    depth::T
    delta::T
    time::T
    dtdd::T
    inc::T
    takeoff::T
    pierce::Vector{T}
    lon::Vector{T}
    lat::Vector{T}
    radius::Vector{T}
end
PhaseGeog(model, name, evlon, evlat, stlon, stlat, depth, delta, time, dtdd, inc, takeoff) =
    PhaseGeog{TauPyFloat}(model, name, evlon, evlat, stlon, stlat, depth, delta, time,
                          dtdd, inc, takeoff, [], [], [], [])
PhaseGeog(p::PhaseGeog, pierce, lon, lat, radius) =
    PhaseGeog{TauPyFloat}(p.model, p.name, p.evlon, p.evlat, p.evlon, p.evlat, p.depth,
                          p.delta, p.time, p.dtdd, p.inc, p.takeoff, pierce, lon, lat, radius)

"""
    path(event_lon, event_lat, depth, station_lon, station_lat, phase="ttall"; model="$DEFAULT_MODEL", cache=true) -> Vector{PhaseGeog}

Create a set of `PhaseGeogs` which contain the computed ray path for a set of `phase`s
from event at (`event_lon`, `event_lat`)° and `depth` km deep, recorded at a station
at (`station_lon`, `station_lat`)°.  Optionally specify the model (one of
$(AVAILABLE_MODELS)).

If no arrivals are found for the geometry, an empty `Vector{Phase}` is returned.

Set `cache` to `false` to always recompute the ray path and ignore the module cache.
"""
function path(event_lon, event_lat, depth, station_lon, station_lat, phase="ttall";
              model=DEFAULT_MODEL, cache=true)
    phase = phase isa AbstractString ? [phase] : phase
    _call_taup(:get_ray_paths_geo, model, depth, event_lat, event_lon,
               station_lat, station_lon, phase; cache=cache)
end

"""
    travel_time(event_lon, event_lat, depth, station_lon, station_lat, phase="ttall"; model="$DEFAULT_MODEL", cache=true) -> Vector{PhaseGeog}

Return a `Vector` of `PhaseGeog`s, given an event `depth` km deep located at
(`event_lon`, `event_lat`)°, recorded at a station at (`station_lon`, `station_lat`)°.
Optionally specify the model (one of: $(AVAILABLE_MODELS)).

If no arrivals are found for the geometry, an empty `Vector{Phase}` is returned.

Set `cache` to `false` to always recompute the ray path and ignore the module cache.
"""
function travel_time(event_lon, event_lat, depth, station_lon, station_lat, phase="ttall";
                     model=DEFAULT_MODEL, cache=true)
    phase = phase isa AbstractString ? [phase] : phase
    _call_taup(:get_travel_times_geo, model, depth, event_lat, event_lon, station_lat, station_lon, phase; cache=cache)
end

"""Helper function which takes `obspy.taup.Arrivals` and returns `PhaseGeog`s.
Arguments should be the same as those passed to osbpy.taup calls, prefixed
by `model`."""
function _phases_from_arrivals(arr, model, depth, event_lat, event_lon,
                               station_lat, station_lon, phase)
    p = Vector{PhaseGeog{TauPyFloat}}()
    for a in arr
        name = a.name
        delta = a.distance
        depth = a.source_depth
        time = a.time
        dtdd = deg2rad(a.ray_param)
        inc = a.incident_angle
        takeoff = a.takeoff_angle
        path = a.path
        lon, lat, radius = if path === nothing
            TauPyFloat[], TauPyFloat[], TauPyFloat[]
        else
            pathlist = path.tolist()
            [pp[6] for pp in pathlist], [pp[5] for pp in pathlist],
                RADIUS[model] .- [pp[4] for pp in pathlist]
        end
        push!(p, PhaseGeog{TauPyFloat}(model, name, event_lon, event_lat, station_lon,
                                       station_lat, depth, delta, time, dtdd, inc, takeoff,
                                       [], lon, lat, radius))
    end
    p
end

max_stars_repo_name: anowacki/TauPy


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
using CartesianMesh
# using Revise
using HDGElasticity

function allequal(u,v)
    return all(u .== v)
end

function allapprox(u,v)
    return all(u .≈ v)
end

function allapprox(u,v,atol)
    @assert length(u) == length(v)
    return all([isapprox(u[i],v[i],atol=atol) for i in 1:length(u)])
end

function distance_function(coords,xc)
    return coords[1,:] .- xc
end

x0 = [0.,0.]
widths = [2.,1.]
nelements = [2,1]
mesh = UniformMesh(x0,widths,nelements)
basis = TensorProductBasis(2,1)
quad1d = ImplicitDomainQuadrature.ReferenceQuadratureRule(2)
poly = InterpolatingPolynomial(1,basis)
NF = HDGElasticity.number_of_basis_functions(basis)
coords = HDGElasticity.nodal_coordinates(mesh,basis)
xc = 0.75
coeffs = reshape(distance_function(coords,xc),NF,:)
cellsign = HDGElasticity.cell_signatures(coeffs,poly)
quads = HDGElasticity.element_quadratures(cellsign,coeffs,poly,quad1d)

py = quad1d.points

p1,w1 = ImplicitDomainQuadrature.transform(quad1d,0.50,1.0)
testp1 = [p1[1]  p1[2]  p1[1]  p1[2]
          py[1]  py[1]  py[2]  py[2]]
testw1 = kron(quad1d.weights,w1)
@test allapprox(testp1,quads[1,1].points)
@test allapprox(testw1,quads[1,1].weights)

p2,w2 = ImplicitDomainQuadrature.transform(quad1d,-1.,0.5)
testp2 = [p2[1]  p2[2]  p2[1]  p2[2]
          py[1]  py[1]  py[2]  py[2]]
testw2 = kron(quad1d.weights,w2)
@test allapprox(testp2,quads[2,1].points)
@test allapprox(testw2,quads[2,1].weights)

qp3 = tensor_product_quadrature(2,2)
@test allequal(qp3.points,quads[1,2].points)
@test allequal(qp3.weights,quads[1,2].weights)
@test !isassigned(quads,2,2)

cellsign[1] = +1
HDGElasticity.element_quadratures!(quads,cellsign,coeffs,poly,quad1d)
@test allequal(qp3.points,quads[1,1].points)

update!(poly,coeffs[:,1])
facemaps = HDGElasticity.reference_cell_facemaps(2)
facequads = reshape([Vector{QuadratureRule{1}}(undef,4) for i = 1:4],2,2)
HDGElasticity.update_face_quadratures!(facequads[1,1],poly,facemaps,+1,quad1d)

testp,testw = ImplicitDomainQuadrature.transform(quad1d,0.5,1.0)
@test allapprox(facequads[1,1][1].points,testp)
@test allapprox(facequads[1,1][1].weights,testw)
@test allapprox(facequads[1,1][3].points,testp)
@test allapprox(facequads[1,1][3].weights,testw)
@test allapprox(facequads[1,1][2].points,quad1d.points)
@test allapprox(facequads[1,1][2].weights,quad1d.weights)
@test length(facequads[1,1][4]) == 0

HDGElasticity.update_face_quadratures!(facequads[2,1],poly,facemaps,-1,quad1d)

p2,w2 = ImplicitDomainQuadrature.transform(quad1d,-1.,0.5)
@test allequal(facequads[2,1][1].points,p2)
@test allequal(facequads[2,1][1].weights,w2)
@test length(facequads[2,1][2]) == 0
@test allequal(facequads[2,1][3].points,p2)
@test allequal(facequads[2,1][3].weights,w2)

@test allapprox(facequads[2,1][4].points,quad1d.points)
@test allapprox(facequads[2,1][4].weights,quad1d.weights)

@test all([!isassigned(facequads[2,2],i) for i = 1:4])

cellsign = HDGElasticity.cell_signatures(coeffs,poly)
facequads = reshape([Vector{QuadratureRule{1}}(undef,4) for i = 1:4],2,2)

HDGElasticity.face_quadratures!(facequads,cellsign,coeffs,poly,facemaps,
    quad1d)

@test allequal(facequads[1,1][1].points,p1)
@test allequal(facequads[1,1][1].weights,w1)
@test allequal(facequads[1,1][3].points,p1)
@test allequal(facequads[1,1][3].weights,w1)
@test allapprox(facequads[1,1][2].points,quad1d.points)
@test allapprox(facequads[1,1][2].weights,quad1d.weights)
@test length(facequads[1,1][4]) == 0
@test allequal(facequads[2,1][1].points,p2)
@test allequal(facequads[2,1][1].weights,w2)
@test allequal(facequads[2,1][3].points,p2)
@test allequal(facequads[2,1][3].weights,w2)
@test allapprox(facequads[2,1][4].points,quad1d.points)
@test allapprox(facequads[2,1][4].weights,quad1d.weights)
# @test all([!isassigned(facequads[2,2],i) for i = 1:4])
@test all([length(fq) == 0 for fq in facequads[2,2]])

facequads = HDGElasticity.face_quadratures(cellsign,coeffs,poly,facemaps,quad1d)

@test allequal(facequads[1,1][1].points,p1)
@test allequal(facequads[1,1][1].weights,w1)
@test allequal(facequads[1,1][3].points,p1)
@test allequal(facequads[1,1][3].weights,w1)
@test allapprox(facequads[1,1][2].points,quad1d.points)
@test allapprox(facequads[1,1][2].weights,quad1d.weights)
@test length(facequads[1,1][4]) == 0
@test allequal(facequads[2,1][1].points,p2)
@test allequal(facequads[2,1][1].weights,w2)
@test allequal(facequads[2,1][3].points,p2)
@test allequal(facequads[2,1][3].weights,w2)
@test allapprox(facequads[2,1][4].points,quad1d.points)
@test allapprox(facequads[2,1][4].weights,quad1d.weights)
# @test all([!isassigned(facequads[2,2],i) for i = 1:4])
@test all([length(fq) == 0 for fq in facequads[2,2]])

basis1d = TensorProductBasis(1,1)
quad1d = tensor_product_quadrature(1,2)
icoeffs = HDGElasticity.interface_coefficients(cellsign,coeffs,poly,basis1d,quad1d)
testicoeffs = [0.5,-1.0,0.5,1.0]
@test allapprox(icoeffs[1],testicoeffs)
@test !isassigned(icoeffs,2)

vbasis = TensorProductBasis(2,4)
sbasis = TensorProductBasis(1,4)
poly = InterpolatingPolynomial(1,vbasis)
quad1d = ImplicitDomainQuadrature.ReferenceQuadratureRule(5)
iquad = tensor_product_quadrature(1,5)
vtpq = tensor_product_quadrature(2,5)
coords = HDGElasticity.nodal_coordinates(mesh,vbasis)
NF = HDGElasticity.number_of_basis_functions(vbasis)
xc = 0.75
coeffs = reshape(distance_function(coords,xc),NF,:)

facemaps = HDGElasticity.reference_cell_facemaps(2)
cellsign = HDGElasticity.cell_signatures(coeffs,poly)

vquads = HDGElasticity.element_quadratures(cellsign,coeffs,poly,quad1d)
fquads = HDGElasticity.face_quadratures(cellsign,coeffs,poly,facemaps,quad1d)
icoeffs = HDGElasticity.interface_coefficients(cellsign,coeffs,poly,sbasis,iquad)

dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
cellmap = HDGElasticity.CellMap(dgmesh.domain[1])
imap = InterpolatingPolynomial(2,sbasis)
inormals = HDGElasticity.interface_normals(dgmesh,icoeffs,imap,
    coeffs,poly,iquad.points)
testnormals = zeros(2,5)
testnormals[1,:] .= 1.0
@test allapprox(testnormals,inormals[1],1e-15)
@test !isassigned(inormals,2)
ftpq = repeat([iquad],4)
fnormals = HDGElasticity.reference_normals()

ufs = HDGElasticity.UniformFunctionSpace(vbasis,sbasis,vtpq,ftpq,vquads,
    fquads,fnormals,icoeffs,iquad,imap,inormals)
@test size(ufs.isactiveface) == (2,2)
@test allequal(ufs.isactiveface[1,1],[1,1,1,0])
@test allequal(ufs.isactiveface[2,1],[1,0,1,1])
@test allequal(ufs.isactiveface[1,2],[1,1,1,1])
@test allequal(ufs.isactiveface[2,2],[0,0,0,0])
@test ufs.dofsperelement == 10

ufs = HDGElasticity.UniformFunctionSpace(dgmesh,4,5,coeffs,poly)
@test ufs.dofsperelement == 10

mesh = UniformMesh([0.,0.],[4.,1.],[2,1])
coords = HDGElasticity.nodal_coordinates(mesh,vbasis)

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end

testn = [1.,1.]/sqrt(2.)
coeffs = reshape(plane_distance_function(coords,testn,[0.5,0.]),NF,:)
dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,4,5,coeffs,poly)
@test sum(ufs.vquads[2,1].weights) ≈ 0.25
@test sum(ufs.vquads[1,1].weights) ≈ 3.75
@test sum(ufs.vquads[1,2].weights) ≈ 4.0
@test !isassigned(ufs.vquads,2,2)
@test ufs.dofsperelement == 10

testnormals = repeat(testn,inner=(1,5))
@test allapprox(testnormals,ufs.inormals[1])
@test !isassigned(ufs.inormals,2)


text:
module SciMLTutorials

using Weave, Pkg, InteractiveUtils, IJulia

repo_directory = joinpath(@__DIR__,"..")
cssfile = joinpath(@__DIR__, "..", "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")

function weave_file(folder,file,build_list=(:script,:html,:pdf,:github,:notebook); kwargs...)
  tmp = joinpath(repo_directory,"tutorials",folder,file)
  Pkg.activate(dirname(tmp))
  Pkg.instantiate()
  args = Dict{Symbol,String}(:folder=>folder,:file=>file)
  if :script ∈ build_list
    println("Building Script")
    dir = joinpath(repo_directory,"script",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "script"
    tangle(tmp;out_path=dir)
  end
  if :html ∈ build_list
    println("Building HTML")
    dir = joinpath(repo_directory,"html",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "html"
    weave(tmp,doctype = "md2html",out_path=dir,args=args; fig_ext=".svg", css=cssfile, kwargs...)
  end
  if :pdf ∈ build_list
    println("Building PDF")
    dir = joinpath(repo_directory,"pdf",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "pdf"
    try
      weave(tmp,doctype="md2pdf",out_path=dir,args=args; template=latexfile, kwargs...)
    catch ex
      @warn "PDF generation failed" exception=(ex, catch_backtrace())
    end
  end
  if :github ∈ build_list
    println("Building Github Markdown")
    dir = joinpath(repo_directory,"markdown",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "github"
    weave(tmp,doctype = "github",out_path=dir,args=args; kwargs...)
  end
  if :notebook ∈ build_list
    println("Building Notebook")
    dir = joinpath(repo_directory,"notebook",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "notebook"
    Weave.convert_doc(tmp,joinpath(dir,file[1:end-4]*".ipynb"))
  end
end

function weave_all()
  for folder in readdir(joinpath(repo_directory,"tutorials"))
    folder == "test.jmd" && continue
    weave_folder(folder)
  end
end

function weave_folder(folder)
  for file in readdir(joinpath(repo_directory,"tutorials",folder))
    println("Building $(joinpath(folder,file)))")
    try
      weave_file(folder,file)
    catch
    end
  end
end

function tutorial_footer(folder=nothing, file=nothing; remove_homedir=true)
    display("text/markdown", """
    ## Appendix

     This tutorial is part of the SciMLTutorials.jl repository, found at: <https://github.com/SciMLTutorials/SciMLTutorials.jl>.
     For more information on doing scientific machine learning (SciML) with open source software, check out <https://sciml.ai/>.
    """)
    if folder !== nothing && file !== nothing
        display("text/markdown", """
        To locally run this tutorial, do the following commands:
        ```
        using SciMLTutorials
        SciMLTutorials.weave_file("$folder","$file")
        ```
        """)
    end
    display("text/markdown", "Computer Information:")
    vinfo = sprint(InteractiveUtils.versioninfo)
    display("text/markdown",  """
    ```
    $(vinfo)
    ```
    """)

    ctx = Pkg.API.Context()
    pkgs = Pkg.Display.status(Pkg.API.Context(), use_as_api=true);
    projfile = ctx.env.project_file
    remove_homedir && (projfile = replace(projfile, homedir() => "~"))

    display("text/markdown","""
    Package Information:
    """)

    md = ""
    md *= "```\nStatus `$(projfile)`\n"

    for pkg in pkgs
        if !isnothing(pkg.old) && pkg.old.ver !== nothing
          md *= "[$(string(pkg.uuid))] $(string(pkg.name)) $(string(pkg.old.ver))\n"
        else
          md *= "[$(string(pkg.uuid))] $(string(pkg.name))\n"
        end
    end
    md *= "```"
    display("text/markdown", md)
end

function open_notebooks()
  Base.eval(Main, Meta.parse("import IJulia"))
  path = joinpath(repo_directory,"notebook")
  IJulia.notebook(;dir=path)
end

end


text:
using Pkg
Pkg.activate("..")
using Profile, TimerOutputs, StatProfilerHTML
include("peps_util.jl")

Nx  = 4
Ny  = 4
mdim = 3
io = open("prof_$(string(Nx))_$mdim.txt", "w+")
logger = SimpleLogger(io)
global_logger(logger)
J  = 1.0
sites = siteinds("S=1/2",Nx*Ny)
println("Beginning A")
A = checkerboardPEPS(sites, Nx, Ny, mindim=mdim)
cA = cuPEPS(A)
H  = makeCuH_XXZ(Nx, Ny, J)
@info "Built cA and H"
Ls = buildLs(cA, H; mindim=mdim, maxdim=mdim)
@info "Built first Ls"
Rs = buildRs(cA, H; mindim=mdim, maxdim=mdim)
@info "Built first Rs"
cA, Ls, Rs = rightwardSweep(cA, Ls, Rs, H; mindim=mdim, maxdim=mdim)
cA, Ls, Rs = leftwardSweep(cA, Ls, Rs, H; mindim=mdim, maxdim=mdim)
Profile.init(n=10_000_000)
@profile doSweeps(cA, Ls, Rs, H; mindim=mdim, maxdim=mdim)
statprofilehtml()


text:
Base.hash(x::ComponentArray, h::UInt) = hash(keys(x), hash(getdata(x), h))

Base.:(==)(x::ComponentArray, y::ComponentArray) = getdata(x)==getdata(y) && getaxes(x)==getaxes(y)
Base.:(==)(x::ComponentArray, y::AbstractArray) = getdata(x)==y && keys(x)==keys(y) # For equality with LabelledArrays
Base.:(==)(x::AbstractArray, y::ComponentArray) = y==x

Base.keys(x::ComponentVector) = keys(indexmap(getaxes(x)[1]))

Base.haskey(x::ComponentVector, s::Symbol) = haskey(indexmap(getaxes(x)[1]), s)

Base.propertynames(x::ComponentVector) = propertynames(indexmap(getaxes(x)[1]))

# Property access for ComponentVectors goes through _get/_setindex
@inline Base.getproperty(x::ComponentVector, s::Symbol) = _getindex(Base.maybeview, x, Val(s))
@inline Base.getproperty(x::ComponentVector, s::Val) = _getindex(Base.maybeview, x, s)

@inline Base.setproperty!(x::ComponentVector, s::Symbol, v) = _setindex!(x, v, Val(s))
@inline Base.setproperty!(x::ComponentVector, s::Val, v) = _setindex!(x, v, s)


text:
using JuMP
using Gurobi
using MathProgBase

include("functions.jl")

# solveDEA: solving dea problems by a smart algorithm
function solveDEA(m::JuMP.Model ; incrementSize = 100 , Tol = 10.0^-6 , lpUB = Inf , extremeValueSetFlag = 0)

  # for i in varStartValue
  #   println("Dictionary Key : $(i[1]) , Value : $(i[2])")
  # end
  #
  # for i in varEndValue
  #   println("Dictionary Key : $(i[1]) , Value : $(i[2])")
  # end
  #
  # for i in varGroup
  #   println("Dictionary Key : $(i[1]) , Group : $(i[2])")
  # end

  # parameter
  # m: model variable of JuMP which is entered by JuMP's scripts
  # incrementSize: smallLP increment size when its size under lpUP for resampling
  # Tol: the accuracy for solving DEA problem
  # lpUB: the limit size of LP
  # extremeValueSetFlag: first sampling take extremeValueSet or not (default =  turn off)

  # data preprocessing
  # getting the data of full-size model which is entered by users
  fullsizeObjSense, fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeNumVar, fullsizeNumConst = getFullsizeModelData(m)

  # println(fullsizeObjSense) # no adjust Min
  # println(size(fullsizeObjCoeff)) # slice <50005,>
  # println(size(fullsizeConstLB)) # slice <19, >
  # println(size(fullsizeConstUB)) # slice <19, >
  # println(size(fullsizeConstMatrix)) # slice <19, 50005>
  # println(size(fullsizeVarLB)) # slice <50005,>
  # println(size(fullsizeVarUB)) # slice <50005,>

  denseFullsizeConstMatrix = full(fullsizeConstMatrix)

  # initial varibles related with data model
  iterations = 0
  dimension, scale = size(denseFullsizeConstMatrix)
  lpSize = 0
  additiveSize = 0
  precision = round(Int,round(log(10, 1/Tol) ,1) )

  # user parameters simple exception catch
  lpUB = (lpUB > scale) ? scale : lpUB
  incrementSize = (incrementSize > lpUB) ? lpUB : incrementSize
  increment = [incrementSize , incrementSize]

  # create dictionary for fulldatamatrix mapping
  fulldatamatrix_dict = Dict([denseFullsizeConstMatrix[1:end,i] => i for i in 1:length(denseFullsizeConstMatrix[1,:])])

  # create critical component and data component from fulldatamatrix
  dataConstMatrix_dict = 0
  criticalConstMatrix_dict = 0

  if varibleInitFlag == false
    dataConstMatrix_dict = find(fullsizeObjCoeff.==0)
    criticalConstMatrix = denseFullsizeConstMatrix[:,find(fullsizeObjCoeff.!=0)]
  else
    for i in varGroup
      if contains(i[2],"nc")
        tmp_dict = collect(varStartValue[i[1]]:varEndValue[i[1]])
        if dataConstMatrix_dict == 0
          dataConstMatrix_dict = tmp_dict
        else
          dataConstMatrix_dict = [dataConstMatrix_dict;tmp_dict]
        end
      else
        tmp_dict = collect(varStartValue[i[1]]:varEndValue[i[1]])
        if criticalConstMatrix_dict == 0
          criticalConstMatrix_dict = tmp_dict
        else
          criticalConstMatrix_dict = [criticalConstMatrix_dict;tmp_dict]
        end
      end
    end
    dataConstMatrix = denseFullsizeConstMatrix[:,dataConstMatrix_dict]
    criticalConstMatrix = denseFullsizeConstMatrix[:,criticalConstMatrix_dict]
    #dataConstMatrix_dict = collect(1:50000)
    #criticalConstMatrix = denseFullsizeConstMatrix[:, collect(50001:50005)]
  end

  # create map for fulldatamatrix mapping to dataConstMatrix
  dataConstMatrix_inv_dict = zeros(Int64, scale)
  for i in 1:length(dataConstMatrix_dict)
    dataConstMatrix_inv_dict[dataConstMatrix_dict[i] ] = i
  end

  # first sampling
  smallLP = criticalConstMatrix[1:end,:]
  critical_size = length(Array(criticalConstMatrix[1,:]))
  lpSize = lpSize + critical_size

  # create map for whether dataConstMatrix used or not
  dataConstMatrix_map = zeros(length(dataConstMatrix[1,:]))
  # record dataConstMatrix map sequence
  dataConstMatrix_map_seq = zeros(lpUB - critical_size)

  if (extremeValueSetFlag == 0)
    # create small LP (Critical data + first datamatrix random data in incrementSize)
    randNum = rand(1:length(dataConstMatrix[1,:]) , incrementSize)
    smallLP = [smallLP  dataConstMatrix[:,randNum[:]] ]
    dataConstMatrix_map[randNum] = 1
    dataConstMatrix_map_seq[1 : incrementSize] = randNum
    lpSize = lpSize + incrementSize
  else
    # create small LP (Critical data + extremeValueSet)
    extremeValueSet = zeros(Int64, dimension)
    for i in 1 : dimension
      extremeValueSet[i] = (fullsizeConstUB[i] == Inf) ? findmin(dataConstMatrix[i,:])[2] : findmax(dataConstMatrix[i,:])[2]
    end
    extremeValueSet =  sort(extremeValueSet)
    smallLP = [smallLP  dataConstMatrix[:,extremeValueSet[:]] ]
    dataConstMatrix_map[extremeValueSet] = 1
    dataConstMatrix_map_seq[1 : dimension] = extremeValueSet
    lpSize = lpSize + dimension
  end

  # generate other small condition matrix
  smallVarLB = zeros(lpSize)
  smallVarUB = zeros(lpSize)
  smallObjCoeff = zeros(lpSize)

  for i in 1 : lpSize
    smallVarLB[i] = fullsizeVarLB[fulldatamatrix_dict[smallLP[:,i]]]
    smallVarUB[i] = fullsizeVarUB[fulldatamatrix_dict[smallLP[:,i]]]
    smallObjCoeff[i] = fullsizeObjCoeff[fulldatamatrix_dict[smallLP[:,i]]]
  end

  smallConstLB = fullsizeConstLB
  smallConstUB = fullsizeConstUB
  smallSense = fullsizeObjSense

  # definition and initialization solving Status
  optStatus = false   # the optimality status of model
  iterations = 0      # the iterations that algorithm uses to solve model

  # JuMP solver initial define
  ### 20170924 need to add other solver that can be used in JuMP
  if solver_Outputflag_check(string(m.solver)) == "GurobiSolver"# check Output flag option for user catch
    setsolver(m , GurobiSolver(OutputFlag = 0))
  end

  lm = MathProgBase.LinearQuadraticModel(m.solver)  # model variable of MathProgBase, we use this to build small model
  stat = :Default   # model status
  m.objVal = NaN      # objective value of full-size model
  m.colVal = fill(NaN, fullsizeNumVar)    # values of all variables of full-size model
  m.linconstrDuals = Array(Float64, 0)    # dual values of full-size model

  # loop until model being solved to optimal or algorithm processing over 1000 iterations
  while optStatus == false
    iterations = iterations + 1
    if iterations > 1000
      warn("Iterations > 1000")
      stat = :Unsolved
      break
    end
    optStatus = true

    # load problems
    lm = MathProgBase.LinearQuadraticModel(m.solver)
    MathProgBase.loadproblem!(lm, smallLP, smallVarLB, smallVarUB, smallObjCoeff, smallConstLB, smallConstUB, smallSense)

    MathProgBase.optimize!(lm)
    stat = MathProgBase.status(lm)

    # get constduals
    constDuals = MathProgBase.getconstrduals(lm)
    # create check optimal array
    checkopt_arr = constDuals.'* fullsizeConstMatrix

    # addictive size calculate
    if lpSize < lpUB # there's room to increase sample size
      additiveSize = increment[min(iterations, length(increment))]
      if lpSize + additiveSize > lpUB
        additiveSize = lpUB - lpSize
      end
    else
        # if the lpSize is already equal to lpUB, then let additiveSize to be 1
        additiveSize = 1
    end

    # check optimal function
    infeasibleIdx = Array(Int, additiveSize)
    infeasibleValue = Array(Float64, additiveSize)

    # create notset for data index which didn't use in dataConstMatrix
    notset = find(dataConstMatrix_map.==0)
    for i in 1:length(notset)
      notset[i] = dataConstMatrix_dict[notset[i] ]
    end

    if smallSense == :Min
      if additiveSize == 1
        infeasibleValue, infeasibleIdx = findmax(checkopt_arr[notset])
        infeasibleIdx = dataConstMatrix_inv_dict[notset[infeasibleIdx] ]
        dataConstMatrix_map[infeasibleIdx] = 1
        optStatus = infeasibleValue[1] <= Tol
      else
        infeasibleValue, infeasibleIdx = largestn(checkopt_arr[notset], additiveSize)
        for i in 1:length(infeasibleIdx)
          infeasibleIdx[i] = dataConstMatrix_inv_dict[notset[infeasibleIdx[i]]]
        end
        dataConstMatrix_map[infeasibleIdx] = 1
        optStatus = findmax(infeasibleValue)[1] <= Tol
      end
    else
      if additiveSize == 1
        infeasibleValue, infeasibleIdx = findmin(checkopt_arr[notset])
        infeasibleIdx = dataConstMatrix_inv_dict[notset[infeasibleIdx] ]
        dataConstMatrix_map[infeasibleIdx] = 1
        optStatus = infeasibleValue[1] >= -Tol
      else
        infeasibleValue, infeasibleIdx = smallestn(checkopt_arr[notset], additiveSize)
        infeasibleIdx = dataConstMatrix_inv_dict[notset[infeasibleIdx] ]
        for i in 1:length(infeasibleIdx)
          infeasibleIdx[i] = dataConstMatrix_inv_dict[notset[infeasibleIdx[i]]]
        end
        optStatus = findmin(infeasibleValue)[1] >= -Tol
      end
    end

    # resampling and dynamic resampling
    if optStatus == false
      if additiveSize == 1
        # create delete index
        critical_length = length(Array(criticalConstMatrix[1,:]))

        deleteIdx = rand(1 : lpUB - critical_length)
        # marked used and not optimal data for "2"
        dataConstMatrix_map[dataConstMatrix_map_seq[deleteIdx]] = 2
        # swap now element at delete index
        dataConstMatrix_map_seq[deleteIdx] = infeasibleIdx
        # swap dataConstMatrix
        smallLP[:,critical_length + deleteIdx] = dataConstMatrix[:,infeasibleIdx]

        # swap other small condition matrix
        smallVarLB[critical_length + deleteIdx] = fullsizeVarLB[fulldatamatrix_dict[smallLP[:,critical_length + deleteIdx]]]
        smallVarUB[critical_length + deleteIdx] = fullsizeVarUB[fulldatamatrix_dict[smallLP[:,critical_length + deleteIdx]]]
        smallObjCoeff[critical_length + deleteIdx] = fullsizeObjCoeff[fulldatamatrix_dict[smallLP[:,critical_length + deleteIdx]]]
      else
        # increment data with incrementSize into smallLP
        smallLP = [smallLP  dataConstMatrix[:,infeasibleIdx] ]

        dataConstMatrix_map_seq[lpSize - critical_size + 1 : lpSize - critical_size + additiveSize] = infeasibleIdx

        # increment data with incrementSize into other small condition matrix
        tmp_smallVarLB = zeros(additiveSize)
        tmp_smallVarUB = zeros(additiveSize)
        tmp_smallObjCoeff = zeros(additiveSize)
        for i in lpSize + 1 : lpSize + additiveSize
          tmp_smallVarLB[i - lpSize] = fullsizeVarLB[fulldatamatrix_dict[smallLP[:,i]]]
          tmp_smallVarUB[i - lpSize] = fullsizeVarUB[fulldatamatrix_dict[smallLP[:,i]]]
          tmp_smallObjCoeff[i - lpSize] = fullsizeObjCoeff[fulldatamatrix_dict[smallLP[:,i]]]
        end

        smallVarLB = [smallVarLB ; tmp_smallVarLB]
        smallVarUB = [smallVarUB ; tmp_smallVarUB]
        smallObjCoeff = [smallObjCoeff ; tmp_smallObjCoeff]

        lpSize = lpSize + additiveSize

      end
    end
  end

  # Extracting solution
  if stat == :Optimal
    # recording objective value
    m.objVal = round(MathProgBase.getobjval(lm), precision)

    # recording clambda (DEA solution)
    m.colVal = fill(0.0, fullsizeNumVar)
    solution = MathProgBase.getsolution(lm)
    for i in 1 : length(solution)
      m.colVal[fulldatamatrix_dict[smallLP[:,i] ] ] = round(solution[i], precision)
    end

    # recording dual value
    m.linconstrDuals = try
        round(MathProgBase.getconstrduals(lm)[1:fullsizeNumConst], precision)
    catch
        fill(NaN, fullsizeNumConst)
    end

    println("iterations : $iterations")

  #   duals = m.linconstrDuals
  #
  #   # recording slack value
  #   slack = zeros(Float64,5)
  #   for k in 1:dimension
  #     for i in 1:length(find(m.colVal.!=0))
  #       slack[k] += m.colVal[find(m.colVal.!=0)[i]] * denseFullsizeConstMatrix[k,find(m.colVal.!=0)[i]]
  #     end
  #     if smallConstLB[k] == -Inf
  #       slack[k] = round(abs(smallConstUB[k] - slack[k]), precision)
  #     else
  #       slack[k] = round(abs(slack[k] - smallConstLB[k]), precision)
  #     end
  #   end
  end
  # return duals,slack

end


text:
# function NormalDifference(df, group, var; ind_summary = mean, ylabel = "Median ...", xyfont = font(18, "Bookman Light"))
#     res_plt, res_group, res_individual = mean_sem_scatter(df, group, var; ind_summary = ind_summary)
#     res_test = test_difference(res_individual, group, var; normality = true)
#     pooledSD = sqrt(((res_group.Central[2])^2 + (res_group.Central[2])^2)/2)
#     res_effect = (res_group.Central[2] - res_group.Central[2]) / pooledSD
#     # if res_effect < 0.5
#     #     res_effect = ((res_test.nx * res_test.ny) - res_test.U) / (res_test.nx * res_test.ny)
#     # end
#     e = round(res_effect; digits =2)
#     add_info!(res_plt, res_group, pvalue(res_test), e; normality = true)
#     xprop = ("Group", xyfont)
#     yprop = (ylabel, xyfont)
#     plot!(xaxis = xprop, yaxis = yprop)
#     return (plt = res_plt,
#         test = res_test,
#         effect = res_effect,
#         groupdf = res_group,
#         individual_df = res_individual)
# end

"""
    `Difference(df, group, var; ind_summary = mean, ylabel = "Median ...", xyfont = font(18, "Bookman Light"), ylims = nothing)`
    A - plot a scatter with median ± CI of the variable var per each value of group
    B - perform MannWhitneyUTest and annotate the results
    C - Calculate the effect size if p < 0.05
    C - return a named tuple with the plot, MannWhitneyUTest, effect size
        a Dataframe with the group median and CI of var,
        and a Dataframe with the individual mouse mean value of var
"""

function Difference(df, group, var; ind_summary = mean, ylabel = "Median ...", xyfont = font(18, "Bookman Light"), ylims = nothing)
    res_plt, res_group, res_individual = median_ci_scatter(df, group, var; ind_summary = ind_summary)
    res_test = MWU_test(res_individual, group, var)
    res_effect = res_test.U/(res_test.nx * res_test.ny)
    if res_effect < 0.5
        res_effect = ((res_test.nx * res_test.ny) - res_test.U) / (res_test.nx * res_test.ny)
    end
    e = round(res_effect; digits =2)
    add_info!(res_plt, res_group, pvalue(res_test), e; ylims = ylims)
    xprop = ("Group", xyfont)
    yprop = (ylabel, xyfont)
    plot!(xaxis = xprop, yaxis = yprop)
    return (plt = res_plt,
        test = res_test,
        effect = res_effect,
        groupdf = res_group,
        individual_df = res_individual)
end

"""
    `median_ci_scatter(df, grouping, var)`
    A - calculate the mean of the variable var for each animal using individual summary
    B - calculate the median and ci for each group using group_summary
    C - Plots the value as a scatter plot ± CI
"""

function median_ci_scatter(df, grouping, var; ind_summary = mean)
    df1 = individual_summary(df, grouping, var; summary = ind_summary)
    df2 = group_summary(df1, grouping, var; normality = false)
    if typeof(grouping) <: AbstractVector && sizeof(grouping) > 1
        df2[!,:xaxis] = [join(x,"_") for x in eachrow(df2[:, grouping])]
        xaxis = :xaxis
    else
        xaxis = grouping
    end
    try
        cases = levels(df[:,grouping])
        df2[!,:color] = [val == cases[1] ? OGCol1 : OGCol2 for val in df2[:,xaxis]]
    catch
        df2[!,:color] .= :gray75
    end
    plt = @df df2 scatter(1:nrow(df2),:Central,
        yerror = :ERR,
        xlims = (0.5, nrow(df2) + 0.5),
        xticks = (1:nrow(df2),string.(cols(xaxis))),
        color = :color,
        legend = false)
    return plt, df2, df1
end

function add_info!(plt, df, p, e; normality = false, ylims = nothing)
    plt_lim = maximum(df.Central .+ last.(df.ERR))
    plt_span = plt_lim/10
    ref, span = plt_lim + plt_span, plt_span/2
    add_bar!(plt, ref, span)
    if p >= 0.05
        pval = "N.S."
    else
        round(p, digits = 2) == 0 ? pval = round(p, digits = 3) : pval = round(p, digits = 2)
        add_effect!(plt, ref+span/2, span/2, e)
    end
    if normality
        pmessage = "T-test, p = $(pval)"
    else
        pmessage = "Mann-Whitney U test, p = $(pval)"
    end
    add_pvalue!(plt, ref+span/2, span/2, pmessage)
    isnothing(ylims) && (ylims = (0,plt_lim + 2plt_span))
    yaxis!(ylims = ylims)
    return plt
end

function add_bar!(plt, ref, span)
    plot!(plt,[1,1],[ref,ref+span], linecolor = :black)
    plot!(plt,[2,2],[ref,ref+span], linecolor = :black)
    plot!(plt,[1,2],[ref+span/2,ref+span/2], linecolor = :black)
end

function add_pvalue!(plt, ref, span, p)
    if isa(p, String)
        message = p
    else
        message = p < 0.01 ? "p < 0.01" : "p < 0.05"
    end
    annotate!(plt,[(1.5,ref + span,
        Plots.text(message,
        8, :center))])
    return plt
end

function add_effect!(plt, ref, span, e)
    message = "Effect size = $e"
    annotate!(plt,[(1.5,ref - span,
        Plots.text(message,
        8, :center))])
    return plt
end

"""
    `function_analysis(df,variable, f; grouping = nothing, step =0.05, calc = :basic,
            color = [:auto], linestyle = [:auto])`
    Apply the function f over the vaariable var per each value of grouping and
    plots the result over the variable var
"""

function function_analysis(df,var, f; grouping = nothing, step =0.05, calc = :basic,
        color = [:auto], linestyle = [:auto])
    subgroups = isnothing(grouping) ? [:MouseID] : vcat(:MouseID,grouping)
    xaxis = range(extrema(df[:, var])..., step = step)
    dd1 = combine(groupby(df,subgroups), var => (t-> f(t,xaxis = xaxis)) => AsTable)
    rename!(dd1, Dict(:Xaxis => var))
    sort!(dd1,[:MouseID,var])
    if calc == :bootstrapping
        dd2 = combine(groupby(dd1,grouping)) do dd3
            group_summary(dd3,var,:fy; normality = false)
        end
        dd2[!,:low] = [x[1] for x in dd2.ERR]
        dd2[!,:up] = [x[2] for x in dd2.ERR]
    elseif calc == :quantiles
        dd2 = combine(groupby(dd1,[grouping,var]), :fy =>(t-> (Central = mean(t),
        low= abs(mean(t) - quantile(t,0.25),
        up = abs(quantile(t,0.975)-mean(t))),
        # ERR = (abs(mean(t) - quantile(t,0.25)) + abs(quantile(t,0.975)-mean(t)))/2,
        SEM = sem(t))) => AsTable)
    elseif calc == :basic
        dd2 = combine(groupby(dd1,[grouping,var]), :fy =>(t-> (Central = mean(t),up = sem(t), low = sem(t))) => AsTable)
    end
    sort!(dd2,var)

    plt = @df dd2 plot(cols(var),:Central, ribbon = (:low, :up), group = cols(grouping), linecolor = :auto, color = color, linestyle = linestyle)
    return plt, dd2
end

function mediansurvival_analysis(streakdf,variable, grouping; plt = plot())
    dd1 = combine(groupby(streakdf,[:MouseID,grouping]), variable => median => variable)
    # dd2 = combine(groupby(dd1,grouping), variable => (t-> (Mean = mean(t),Sem = sem(t))) => AsTable)
    dd2 = group_summary(dd1,grouping,variable; normality = false)
    # dd2[!,:low] = [x[1] for x in dd2.ERR]
    # dd2[!,:up] = [x[2] for x in dd2.ERR]
    isa(dd2[:, grouping], CategoricalArray) && (dd2[!,grouping] = string.(dd2[!,grouping]))
    println(isa(dd2[:, grouping], CategoricalArray))
    @df dd2 scatter!(plt,cols(grouping), :Central, yerror = :ERR,
        # xlims = (-0.25,2.25), xlabel = "Group",
        ylabel = "Median survival time", label = "")
    return plt
end


text:
_length(iter) = _length(iter, Base.IteratorSize(iter))
_length(iter, ::Base.HasLength) = length(iter)
_length(iter, ::Base.HasShape)  = length(iter)
_length(iter, ::Base.IsInfinite) = Inf
_length(iter, ::Base.SizeUnknown) = "NA"

_length_str(iter) = _length_str(iter, Base.IteratorSize(iter))
_length_str(iter, ::Base.HasLength) = string(length(iter), ", ")
_length_str(iter, ::Base.HasShape)  = string(length(iter), ", ")
_length_str(iter, ::Base.IsInfinite) = ""
_length_str(iter, ::Base.SizeUnknown) = ""

Base.IteratorEltype(::Type{T}) where {T<:DataIterator} = Base.HasEltype()
Base.eltype(::Type{DataIterator{E,T}}) where {E,T} = E

# There is no contract that says these methods will work
# It may be that some DataIterator subtypes do not support getindex
# and/or don't support collect
getobs(A::AbstractDataIterator) = map(getobs,collect(A))
getobs(A::AbstractDataIterator, i) = getobs(A[i])
getobs(A::AbstractDataIterator{T}) where {T<:Tuple} = map(x->map(getobs,x), collect(A))
getobs(A::AbstractDataIterator{T}, i::Integer) where {T<:Tuple} = getobs.(A[i])

DataSubset(data::T, indices, obsdim) where {T<:DataIterator} =
    throw(MethodError(DataSubset, (data,indices,obsdim)))

# To avoid overflow when infinite
_next_idx(iter, idx) = _next_idx(Base.IteratorSize(iter), idx)
_next_idx(::Base.IteratorSize, idx) = idx + 1
_next_idx(::Base.IsInfinite, idx) = 1

# --------------------------------------------------------------------
# ObsIterator

function Base.show(io::IO, iter::ObsIterator{E,T}) where {E,T}
    if get(io, :compact, false)
        print(io, typeof(iter).name.name, "{", E, ",", T, "} with " , _length(iter), " observations")
    else
        print(io, summary(iter), "\n Iterator providing ", _length(iter), " observations")
    end
end

# --------------------------------------------------------------------
# BatchIterator

function Base.show(io::IO, iter::BatchIterator{E,T}) where {E,T}
    if get(io, :compact, false)
        print(io, typeof(iter).name.name, "{", E, ",", T, "} with " , _length(iter), " batches")
    else
        print(io, summary(iter), "\n Iterator providing ", _length(iter), " batches of size ", batchsize(iter))
    end
end

# --------------------------------------------------------------------

"""
    RandomObs(data, [count], [obsdim])

Description
============

Create an iterator that generates `count` randomly sampled
observations from `data`. In the case `count` is not provided,
it will generate random samples indefinitely.

Arguments
==========

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`nobs`](@ref) (see Details for more information).

- **`count`** : Optional. The number of randomly sampled
    observations that the iterator will generate before stopping.
    If omitted, the iterator will generate randomly sampled
    observations forever.

- **`obsdim`** : Optional. If it makes sense for the type of
    `data`, `obsdim` can be used to specify which dimension of
    `data` denotes the observations. It can be specified in a
    type-stable manner as a positional argument (see
    `?LearnBase.ObsDim`), or more conveniently as a smart keyword
    argument.

Details
========

For `RandomObs` to work on some data structure, the type of the
given variable `data` must implement the data container
interface. See `?DataSubset` for more info.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Tom Breloff (Github: https://github.com/tbreloff)

Examples
=========

```julia
X, Y = MLDataUtils.load_iris()

# go over 500 randomly sampled observations in X
i = 0
for x in RandomObs(X, 500) # also: RandomObs(X, count = 500)
    @assert typeof(x) <: SubArray{Float64,1}
    @assert length(x) == 4
    i += 1
end
@assert i == 500

# if no count it provided the iterator will generate samples forever
for x in RandomObs(X)
    # this loop will never stop unless break is used
    if true; break; end
end

# also works for multiple data arguments (e.g. labeled data)
for (x,y) in RandomObs((X,Y), count = 100)
    @assert typeof(x) <: SubArray{Float64,1}
    @assert typeof(y) <: String
end
```

see also
=========

[`BalancedObs`](@ref), [`RandomBatches`](@ref),
[`ObsView`](@ref), [`BatchView`](@ref), [`shuffleobs`](@ref),
[`DataSubset`](@ref), [`BufferGetObs`](@ref)
"""
struct RandomObs{E,T,O,I} <: ObsIterator{E,T}
    data::T
    count::Int
    obsdim::O
end

function RandomObs(data::T, count::Int, obsdim::O) where {T,O}
    count > 0 || throw(ArgumentError("count has to be greater than 0"))
    E = typeof(datasubset(data, 1, obsdim))
    RandomObs{E,T,O,Base.HasLength}(data, count, obsdim)
end

function RandomObs(data::T, obsdim::O) where {T,O}
    E = typeof(datasubset(data, 1, obsdim))
    RandomObs{E,T,O,Base.IsInfinite}(data, 1337, obsdim)
end

RandomObs(data::T, count::Int; obsdim = default_obsdim(data)) where {T} =
    RandomObs(data, count, convert(LearnBase.ObsDimension,obsdim))

# convenience constructor.
RandomObs(data, ::Nothing, obsdim) = RandomObs(data, obsdim)
function RandomObs(data; count = nothing, obsdim = default_obsdim(data))
    RandomObs(data, count, convert(LearnBase.ObsDimension,obsdim))
end

function Base.iterate(iter::RandomObs, idx = 1)
    idx > _length(iter) && return nothing
    (datasubset(iter.data, rand(1:nobs(iter.data, iter.obsdim)), iter.obsdim),
     _next_idx(iter,idx))
end

Base.eltype(::Type{RandomObs{E,T,O,I}}) where {E,T,O,I} = E
Base.IteratorSize(::Type{RandomObs{E,T,O,I}}) where {E,T,O,I} = I()
Base.length(iter::RandomObs{E,T,O,Base.HasLength}) where {E,T,O} = iter.count
nobs(iter::RandomObs) = nobs(iter.data, iter.obsdim)

function Base.summary(iter::RandomObs)
    io = IOBuffer()
    print(io, typeof(iter).name.name, "(")
    Base.showarg(io, iter.data, false)
    print(io, ", ", _length_str(iter))
    print(io, replace(string(iter.obsdim), "LearnBase." => ""))
    print(io, ')')
    first(readlines(seek(io,0)))
end

# --------------------------------------------------------------------

"""
    BalancedObs([f], data, [count], [obsdim])

Description
============

Create an iterator that generates `count` randomly sampled
observations from `data`. In the case `count` is not provided,
it will generate random samples indefinitely.

In contrast to [`RandomObs`](@ref), `BalancedObs` expects `data`
to be a labeled data container. It uses the label distribution of
`data` to make sure every label has an equal probability to be
sampled from.

Arguments
==========

- **`f`** : Optional. A function that should be applied to each
    observation individually in order to extract or compute the
    target for that observation. This function is only used once
    during construction to determine which label each observation
    belongs to.

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref),
    [`nobs`](@ref), and optionally [`gettargets`](@ref) (see
    Details for more information).

- **`count`** : Optional. The number of randomly sampled
    observations that the iterator will generate before stopping.
    If omitted, the iterator will generate randomly sampled
    observations forever.

- **`obsdim`** : Optional. If it makes sense for the type of
    `data`, `obsdim` can be used to specify which dimension of
    `data` denotes the observations. It can be specified in a
    type-stable manner as a positional argument (see
    `?LearnBase.ObsDim`), or more conveniently as a smart keyword
    argument.

Details
========

For `BalancedObs` to work on some data structure, the type of the
given variable `data` must implement the labeled data container
interface. See `?DataSubset` for more info.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
=========

```julia
# load first 55 observations of the iris data as example
# - 50 observations for "setosa"
# -  5 observations for "versicolor"
X, Y = MLDataUtils.load_iris(55)

# go over 100 balanced samples observations in X
num_versicolor = 0
for (x,y) in BalancedObs((X,Y), 100) # also: BalancedObs((X,Y), count = 100)
    @assert typeof(x) <: SubArray{Float64,1}
    @assert length(x) == 4
    # count how many times we sample a versicolor observation
    num_versicolor += y[] == "versicolor" ? 1 : 0
end
println(num_versicolor) # around 50

# if no count it provided the iterator will generate samples forever
for (x,y) in BalancedObs((X,Y))
    # this loop will never stop unless break is used
    if true; break; end
end
```

see also
=========

[`RandomObs`](@ref), [`targets`](@ref), [`RandomBatches`](@ref),
[`ObsView`](@ref), [`BatchView`](@ref), [`shuffleobs`](@ref),
[`DataSubset`](@ref), [`BufferGetObs`](@ref)
"""
struct BalancedObs{E,T,L,O,I} <: ObsIterator{E,T}
    data::T
    labelmap::L
    count::Int
    obsdim::O
end

function BalancedObs(f::Function, data, count::Int, obsdim)
    count > 0 || throw(ArgumentError("count has to be greater than 0"))
    lm = labelmap(eachtarget(f, data, obsdim))
    E = typeof(datasubset(data, 1, obsdim))
    BalancedObs{E,typeof(data),typeof(lm),typeof(obsdim),Base.HasLength}(data, lm, count, obsdim)
end
BalancedObs(data, count::Int, obsdim) = BalancedObs(identity, data, count, obsdim)

function BalancedObs(f::Function, data, obsdim)
    lm = labelmap(eachtarget(f, data, obsdim))
    E = typeof(datasubset(data, 1, obsdim))
    BalancedObs{E,typeof(data),typeof(lm),typeof(obsdim),Base.IsInfinite}(data, lm, 1337, obsdim)
end
BalancedObs(data, obsdim) = BalancedObs(identity, data, obsdim)

BalancedObs(data, count::Int; obsdim = default_obsdim(data)) =
    BalancedObs(data, count, convert(LearnBase.ObsDimension,obsdim))
BalancedObs(f::Function, data, count::Int; obsdim = default_obsdim(data)) =
    BalancedObs(f, data, count, convert(LearnBase.ObsDimension,obsdim))

# convenience constructor.
BalancedObs(data, ::Nothing, obsdim) = BalancedObs(data, obsdim)
BalancedObs(f::Function, data, ::Nothing, obsdim) = BalancedObs(f, data, obsdim)
function BalancedObs(data; count = nothing, obsdim = default_obsdim(data))
    BalancedObs(data, count, convert(LearnBase.ObsDimension,obsdim))
end
function BalancedObs(f::Function, data; count = nothing, obsdim = default_obsdim(data))
    BalancedObs(f, data, count, convert(LearnBase.ObsDimension,obsdim))
end

function Base.iterate(iter::BalancedObs, idx = 1)
    idx > _length(iter) && return nothing
    # uniformly random select a label
    obsidx = rand(iter.labelmap).second
    # uniformly random select observation from that label
    (datasubset(iter.data, rand(obsidx), iter.obsdim),
     _next_idx(iter,idx))
end

Base.eltype(::Type{BalancedObs{E,T,L,O,I}}) where {E,T,L,O,I} = E
Base.IteratorSize(::Type{BalancedObs{E,T,L,O,I}}) where {E,T,L,O,I} = I()
Base.length(iter::BalancedObs{E,T,L,O,Base.HasLength}) where {E,T,L,O} = iter.count
nobs(iter::BalancedObs) = nobs(iter.data, iter.obsdim)

function Base.summary(iter::BalancedObs)
    io = IOBuffer()
    print(io, typeof(iter).name.name, '(')
    Base.showarg(io, iter.data, false)
    print(io, ", ", _length_str(iter))
    print(io, replace(string(iter.obsdim), "LearnBase." => ""))
    print(io, ')')
    first(readlines(seek(io,0)))
end

# --------------------------------------------------------------------

"""
    RandomBatches(data, [size], [count], [obsdim])

Description
============

Create an iterator that generates `count` randomly sampled
batches from `data` with a batch-size of `size` .
In the case `count` is not provided, it will generate random
batches indefinitely.

Arguments
==========

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`nobs`](@ref) (see Details for more information).

- **`size`** : Optional. The batch-size of each batch.
    I.e. the number of randomly sampled observations in each batch

- **`count`** : Optional. The number of randomly sampled batches
    that the iterator will generate before stopping. If omitted,
    the iterator will generate randomly sampled batches forever.

- **`obsdim`** : Optional. If it makes sense for the type of
    `data`, `obsdim` can be used to specify which dimension of
    `data` denotes the observations. It can be specified in a
    type-stable manner as a positional argument (see
    `?LearnBase.ObsDim`), or more conveniently as a smart keyword
    argument.

Details
========

For `RandomBatches` to work on some data structure, the type of
given variable `data` must implement the data container
interface. See `?DataSubset` for more info.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)
- Tom Breloff (Github: https://github.com/tbreloff)

Examples
=========

```julia
X, Y = MLDataUtils.load_iris()

# go over 500 randomly sampled batches of batchsize 10
i = 0
for x in RandomBatches(X, 10, 500) # also: RandomObs(X, size = 10, count = 500)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert size(x) == (4,10)
    i += 1
end
@assert i == 500

# if no count it provided the iterator will generate samples forever
for x in RandomBatches(X, 10)
    # this loop will never stop unless break is used
    if true; break; end
end

# also works for multiple data arguments (e.g. labeled data)
for (x,y) in RandomBatches((X,Y), 10, 500)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
end
```

see also
=========

[`RandomObs`](@ref), [`BatchView`](@ref), [`ObsView`](@ref),
[`shuffleobs`](@ref), [`DataSubset`](@ref), [`BufferGetObs`](@ref)
"""
struct RandomBatches{E,T,O,I} <: BatchIterator{E,T}
    data::T
    size::Int
    count::Int
    obsdim::O
end

function RandomBatches(data::T, size::Int, count::Int, obsdim::O) where {T,O}
    size  > 0 || throw(ArgumentError("size has to be greater than 0"))
    count > 0 || throw(ArgumentError("count has to be greater than 0"))
    E = typeof(datasubset(data, rand(1:nobs(data,obsdim), size), obsdim))
    RandomBatches{E,T,O,Base.HasLength}(data, size, count, obsdim)
end

function RandomBatches(data::T, size::Int, obsdim::O) where {T,O}
    size > 0 || throw(ArgumentError("size has to be greater than 0"))
    E = typeof(datasubset(data, rand(1:nobs(data,obsdim), size), obsdim))
    RandomBatches{E,T,O,Base.IsInfinite}(data, size, 1337, obsdim)
end

RandomBatches(data::T, size::Int; obsdim = default_obsdim(data)) where {T} =
    RandomBatches(data, size, convert(LearnBase.ObsDimension,obsdim))

RandomBatches(data::T, size::Int, count::Int; obsdim = default_obsdim(data)) where {T} =
    RandomBatches(data, size, count, convert(LearnBase.ObsDimension,obsdim))

# convenience constructor.
RandomBatches(data::T, size::Int, ::Nothing, obsdim) where {T} =
    RandomBatches(data, size, obsdim)

function RandomBatches(data::T; size::Int = -1, count = nothing, obsdim = default_obsdim(data)) where T
    nobsdim = convert(LearnBase.ObsDimension,obsdim)
    nsize = size < 0 ? default_batch_size(data, nobsdim)::Int : size
    RandomBatches(data, nsize, count, nobsdim)
end

function Base.iterate(iter::RandomBatches, idx = 1)
    idx > _length(iter) && return nothing
    # maybe use StatsBase.sample instead of rand in order to avoid
    # replacement. That said I would like to avoid keyword arguments
    # and currently sample needs "replace" to be specified as such
    indices = rand(1:nobs(iter.data, iter.obsdim), iter.size)
    (datasubset(iter.data, indices, iter.obsdim), _next_idx(iter, idx))
end

Base.eltype(::Type{RandomBatches{E,T,O,I}}) where {E,T,O,I} = E
Base.IteratorSize(::Type{RandomBatches{E,T,O,I}}) where {E,T,O,I} = I()
Base.length(iter::RandomBatches{E,T,O,Base.HasLength}) where {E,T,O} = iter.count
nobs(iter::RandomBatches) = nobs(iter.data, iter.obsdim)
batchsize(iter::RandomBatches) = iter.size

function Base.summary(iter::RandomBatches)
    io = IOBuffer()
    print(io, typeof(iter).name.name, "(")
    Base.showarg(io, iter.data, false)
    print(io, ", ")
    print(io, iter.size, ", ")
    print(io, _length_str(iter))
    print(io, replace(string(iter.obsdim), "LearnBase." => ""))
    print(io, ')')
    first(readlines(seek(io,0)))
end

# --------------------------------------------------------------------

"""
    BufferGetObs(iterator, [buffer])

A stateful iterator that stores the output of
`iterate(iterator,state)` into `buffer` using `getobs!(buffer,
...)`. Depending on the type of data provided by `iterator` this
may be more memory efficient than `getobs(...)`. In the case of
array data, for example, this allows for cache-efficient
processing of each element without allocating a temporary array.

Note that not all types of data support buffering, because it is
the author's choice to opt-in and implement a custom `getobs!`.
For those types that do not provide a custom `getobs!`, the
buffer will be ignored and the result of `getobs(...)` returned.

see [`eachobs`](@ref) and [`eachbatch`](@ref) for concrete examples.
"""
struct BufferGetObs{TElem,TIter}
    iter::TIter
    buffer::TElem
end

function BufferGetObs(iter::T) where T
    buffer = _getobs_eltype(iter, eltype(T))
    BufferGetObs{typeof(buffer),T}(iter, buffer)
end

_getobs_eltype(iter, ::Type) = getobs(first(iter))
_getobs_eltype(iter, ::Type{T}) where {T<:Tuple} = map(getobs, first(iter))

function Base.show(io::IO, iter::BufferGetObs{E,T}) where {E,T}
    if get(io, :compact, false)
        print(io, typeof(iter).name.name, "{", E, ",", T, "} with " , _length(iter), " elements")
    else
        print(io, summary(iter), "\n Iterator providing ", _length(iter), " elements")
    end
end

# TODO: simplify this mess
function Base.iterate(b::BufferGetObs)
    inner = Base.iterate(b.iter)
    inner == nothing && return nothing
    subset, nidx = inner
    (getobs!(b.buffer, subset), nidx)
end
function Base.iterate(b::BufferGetObs{T}) where T<:Tuple
    inner = Base.iterate(b.iter)
    inner == nothing && return nothing
    subset, nidx = inner
    (map(getobs!, b.buffer, subset), nidx)
end
function Base.iterate(b::BufferGetObs, idx)
    inner = Base.iterate(b.iter, idx)
    inner == nothing && return nothing
    subset, nidx = inner
    (getobs!(b.buffer, subset), nidx)
end
function Base.iterate(b::BufferGetObs{T}, idx) where T<:Tuple
    inner = Base.iterate(b.iter, idx)
    inner == nothing && return nothing
    subset, nidx = inner
    (map(getobs!, b.buffer, subset), nidx)
end

Base.eltype(::Type{BufferGetObs{E,T}}) where {E,T} = E
Base.IteratorSize(::Type{BufferGetObs{E,T}}) where {E,T} = Base.IteratorSize(T)
Base.length(b::BufferGetObs) = length(b.iter)
Base.size(b::BufferGetObs, I...) = size(b.iter, I...)
nobs(b::BufferGetObs) = nobs(b.iter)
batchsize(b::BufferGetObs) = batchsize(b.iter)

function Base.summary(b::BufferGetObs)
    io = IOBuffer()
    print(io, typeof(b).name.name, "(")
    Base.showarg(io, b.iter, false)
    print(io, ", ")
    Base.showarg(io, b.buffer, false)
    print(io, ')')
    first(readlines(seek(io,0)))
end

# --------------------------------------------------------------------

"""
    eachobs(data, [obsdim])

Iterate over `data` one observation at a time. If supported by
the type of `data`, a buffer will be preallocated and reused for
memory efficiency.

IMPORTANT: Avoid using `collect`, because in general each
iteration could return the same object with mutated values.
If that behaviour is undesired use `obsview` instead.

```julia
X = rand(4,100)
for x in eachobs(X)
    # loop entered 100 times
    @assert typeof(x) <: Vector{Float64}
    @assert size(x) == (4,)
end
```

In the case of arrays it is assumed that the observations are
represented by the last array dimension. This can be overwritten.

```julia
# This time flip the dimensions of the matrix
X = rand(100,4)
A = eachobs(X, obsdim=1)
# The behaviour remains the same as before
@assert eltype(A) <: Array{Float64,1}
@assert length(A) == 100
```

Multiple variables are supported (e.g. for labeled data)

```julia
for (x,y) in eachobs((X,Y))
    # ...
end
```

Note that internally `eachobs(data, obsdim)` maps to
`BufferGetObs(obsview(data, obsdim))`.

```julia
@assert typeof(eachobs(X)) <: BufferGetObs
@assert typeof(eachobs(X).iter) <: ObsView
```

This means that the following code:

```julia
for obs in eachobs(data, obsdim)
    # ...
end
```

is roughly equivalent to:

```julia
obs = getobs(data, 1, obsdim) # use first element to preallocate buffer
for _ in obsview(data, obsdim)
    getobs!(obs, _) # reuse buffer each iteration
    # ...
end
```

see [`BufferGetObs`](@ref), [`obsview`](@ref), and
[`getobs!`](@ref) for more info. also see [`eachbatch`](@ref) for
a mini-batch version.
"""
eachobs(data, obsdim) = BufferGetObs(ObsView(data, obsdim))
eachobs(data; obsdim = default_obsdim(data)) = eachobs(data, convert(LearnBase.ObsDimension,obsdim))

# --------------------------------------------------------------------

"""
    eachbatch(data, [size], [count], [obsdim])

Iterate over `data` one batch at a time. If supported by the type
of `data`, a buffer will be preallocated and reused for memory
efficiency.

IMPORTANT: Avoid using `collect`, because in general each
iteration could return the same object with mutated values.
If that behaviour is undesired use `BatchView` instead.

The (constant) batch-size can be either provided directly using
`size` or indirectly using `count`, which derives `size` based on
`nobs`. In the case that the size of the dataset is not dividable
by the specified (or inferred) `size`, the remaining observations
will be ignored.

```julia
X = rand(4,150)
for x in eachbatch(X, size = 10) # or: eachbatch(X, count = 15)
    # loop entered 15 times
    @assert typeof(x) <: Matrix{Float64}
    @assert size(x) == (4,10)
end
```

In the case of arrays it is assumed that the observations are
represented by the last array dimension. This can be overwritten.

```julia
# This time flip the dimensions of the matrix
X = rand(150,4)
A = eachbatch(X, size = 10, obsdim = 1)
# The behaviour remains the same as before
@assert eltype(A) <: Array{Float64,2}
@assert length(A) == 15
```

Multiple variables are supported (e.g. for labeled data)

```julia
for (x,y) in eachbatch((X,Y))
    # ...
end
```

Note that internally `eachbatch(data, ...)` maps to
`BufferGetObs(batchview(data, ...))`.

```julia
@assert typeof(eachbatch(X)) <: BufferGetObs
@assert typeof(eachbatch(X).iter) <: BatchView
```

This means that the following code:

```julia
for batch in eachbatch(data, batchsize, obsdim)
    # ...
end
```

is roughly equivalent to:

```julia
batch = getobs(data, collect(1:batchsize), obsdim) # use first element to preallocate buffer
for _ in batchview(data, batchsize, obsdim)
    getobs!(batch, _) # reuse buffer each iteration
    # ...
end
```

see [`BufferGetObs`](@ref), [`batchview`](@ref), and
[`getobs!`](@ref) for more info. also see [`eachobs`](@ref) for a
single-observation version.
"""
function eachbatch(data; size = -1, maxsize = -1, count = -1, obsdim = default_obsdim(data))
    maxsize != -1 && size != -1 && throw(ArgumentError("Providing both \"size\" and \"maxsize\" is not supported"))
    if maxsize != -1
        # set upto to true in order to allow a flexible batch size
        BufferGetObs(BatchView(data, maxsize, count, convert(LearnBase.ObsDimension,obsdim), true))
    else
        # force given batch size
        BufferGetObs(BatchView(data, size, count, convert(LearnBase.ObsDimension,obsdim)))
    end
end

eachbatch(data, obsdim::T) where {T<:Union{Tuple,ObsDimension}} =
    BufferGetObs(BatchView(data, -1, -1, obsdim))

eachbatch(data, size::Int, obsdim::T = default_obsdim(data), upto::Bool = false) where {T<:Union{Tuple,ObsDimension}} =
    BufferGetObs(BatchView(data, size, -1, obsdim, upto))

eachbatch(data, size::Int, count::Int, obsdim::T = default_obsdim(data), upto::Bool = false) where {T<:Union{Tuple,ObsDimension}} =
    BufferGetObs(BatchView(data, size, count, obsdim, upto))


text:
using YahooFinance

get_symbols("GOOG", "2018-12-20", "2019-12-15")


text:
module GAPLib

export GAPData, GAPAgent, GAPJob, loadGAP

const data_path = joinpath(pkgdir(GAPLib), "data")

using ZipFile

include("Data.jl")
include("Loader.jl")

end


text:
import MitosisStochasticDiffEq as MSDE
using StochasticDiffEq
using Mitosis
using Test, Random
using Statistics

K = 200_000
Random.seed!(1234)

# define SDE function
f(u,p,t) = p[1]*u .+ p[2]
g(u,p,t) = p[3]

# time span
tstart = 0.0
tend = 1.0
dt = 0.01
trange = tstart:dt:tend

# initial condition
u0 = 1.1

# set true model parameters
par = [-0.2, 0.1, 0.9]

# set of linear parameters Eq.~(2.2)
plin = copy(par)
sdekernel = MSDE.SDEKernel(f,g,trange,par)

@testset "ensemble sampling tests" begin
  samples1 = MSDE.sample(sdekernel, u0, K, save_noise=false)
  samples2 = [MSDE.sample(sdekernel, u0, save_noise=false)[1] for _ in 1:K]

  @test isapprox(mean(samples1), mean(samples2), rtol=1e-2)
  @test isapprox(cov(samples1), cov(samples2), rtol=1e-2)
end

# initial values for ODE
mynames = (:logscale, :μ, :Σ);
myvalues = [0.0, 1.5, 0.1];
NT = NamedTuple{mynames}(myvalues)

sdetildekernel = MSDE.SDEKernel(f,g,trange,plin)
message, backward = MSDE.backwardfilter(sdetildekernel, NT)

sdetildekernel2 = MSDE.SDEKernel(Mitosis.AffineMap(plin[1], plin[2]), Mitosis.ConstantMap(plin[3]), trange, plin)
message2, backward2 = MSDE.backwardfilter(sdetildekernel2, NT)

@testset "ensemble backward tests" begin
  @test isapprox(backward, backward2, rtol=1e-15)
  @test isapprox(Array(message.sol), Array(message2.sol), rtol=1e-15)
end

x0 = 1.34
ll0 = randn()

samples1 = []
for k=1:K
  push!(samples1, MSDE.forwardguiding(sdekernel, message, (x0, ll0))[1][1])
end

samples2 = MSDE.forwardguiding(sdekernel, message, (x0, ll0), numtraj=K)[1][1]

@testset "ensemble guiding tests" begin
  @test isapprox(mean(samples1), mean(samples2), rtol=1e-2)
  @test isapprox(cov(samples1), cov(samples2), rtol=1e-2)
end


text:
if Sys.islinux()
    max_watches = parse(Int, chomp(read(`sysctl -n fs.inotify.max_user_watches`, String)))
    if max_watches <= 8192 && !isfile(joinpath(@__DIR__, "user_watches"))
        default_watches = 65536
        msg = """
Revise needs to be notified by your filesystem about changes to your code. Your operating
system imposes a limit of $max_watches watched files and directories. Exceeding this limit
can cause errors, `ERROR: start_watching (File Monitor): no space left on device (ENOSPC)`.
To prevent this from happening, Revise will try to increase the number to $default_watches;
however, if you encounter the above error, consider increasing the number even further.
See the documentation for details.

Changing the number of watches is an operation that requires `sudo` privileges.
"""
        # Because output is redirected to a file, we ensure the message is printed on the screen
        run(pipeline(`echo $msg`, stdout="/dev/tty"))
        # Check whether the user has sudo privileges (will prompt for password)
        has_sudo = try success(`sudo -v`) catch err false end
        if has_sudo
            run(pipeline(`echo $default_watches`, `sudo tee -a /proc/sys/fs/inotify/max_user_watches`))
        else
            print(msg)  # dumps msg to the log file
            @warn "You lack sudo privileges, consider contacting your system administrator."
        end
    end
end


text:
using YaoBlocks.AD, YaoBlocks
using Test

@testset "patches" begin
    include("patches.jl")
end

@testset "outerproduct_and_projection" begin
    include("outerproduct_and_projection.jl")
end

@testset "NoParams" begin
    include("NoParams.jl")
end

@testset "apply_back" begin
    include("apply_back.jl")
end

@testset "mat_back" begin
    include("mat_back.jl")
end

@testset "specializes" begin
    include("specializes.jl")
end

@testset "chainrules_patch" begin
    include("chainrules_patch.jl")
end

max_stars_repo_name: QuantumBFS/YaoStdCircuit

max_issues_repo_name: QuantumBFS/YaoStdCircuit


text:
transpose(A::GBVecOrMat) = Transpose(A)

@kwargs function Base.copy!(C::GBMatrix, A::LinearAlgebra.Transpose{<:Any, <:GBMatrix})
    libgb.GrB_transpose(C, mask, accum, A.parent, desc)
end

@kwargs function Base.copy(A::LinearAlgebra.Transpose{<:Any, <:GBMatrix})
    C = similar(A.parent, size(A.parent,2), size(A.parent, 1))
    copy!(C, A)
    return C
end

function _handletranspose(
    A::GBVecMatOrTrans,
    desc::Union{Descriptor, Nothing} = nothing,
    B::Union{GBVecMatOrTrans, Nothing} = nothing
)
    if A isa Transpose
        desc = desc + Descriptors.T0
        A = A.parent
    end
    if B isa Transpose
        desc = desc + Descriptors.T1
        B = B.parent
    end
    @show desc
    return A, desc, B
end


text:
@enum Severity::UInt8 Asymptomatic=1 Mild Severe Critical # 2bits

@enum HealthState::UInt8 Healthy Incubating Infectious MildSymptoms SevereSymptoms CriticalSymptoms Recovered Dead # 3 bits

@enum FreedomState::UInt8 Free HomeQuarantine HomeTreatment Hospitalized Released # 3 bits

@enum DetectionStatus::UInt8 Undetected UnderObservation TestPending Detected #2 bits

@enum ContactKind::UInt8 NoContact=0 HouseholdContact HospitalContact FriendshipContact SporadicContact ConstantKernelContact OutsideContact # 3 bits

@enum DetectionKind::UInt8 NoDetection=0 OutsideQuarantineDetction=1 FromQuarantineDetection FromTrackingDetection

@enum TrackingKind::UInt8 begin
  NotTracked=0
  MildCaseTracked #unused yet
  HospitalTracked #unused yet
  QuarantineTracked #unused yet
  ClassicalTracked
  PhoneTracked
end


text:
#!/usr/local/bin/julia

# ./Mesh_test.jl

push!(LOAD_PATH, "/home/matt/Documents/julia/3dprinting/")

using Mesh
using STL

n = Net()

v1 = vertex!(n, 0.0, 0.0, 0.0)
v2 = vertex!(n, 0.0, 1.0, 0.0)
v3 = vertex!(n, 3.0, 0.0, 0.0)
v4 = vertex!(n, 0.0, 0.0, 1.0)
face!(n, v1, v2, v3)
face!(n, v3, v4, v2)
face!(n, v1, v4, v3)
face!(n, v4, v1, v2)

STL_ASCII(n, "tri.stl")


text:
@inline logdet_coef(y::AbstractFixedSizeVector, ::Any, σ::AbstractFixedSizeVector) = One()

function univariate_normal_quote(
    M::Int, ::Type{T}, yisvec::Bool,
    @nospecialize(μisvec::Union{Bool,Nothing}), @nospecialize(σisvec::Union{Bool,Nothing}),
    (track_y, track_μ, track_σ)::NTuple{3,Bool},
    (inity, initμ, initσ)::NTuple{3,Bool},
    partial::Bool, calclogdet::Bool = false
) where {T}
    if M > 1
        pre_quote = quote
            qf = SIMDPirates.vbroadcast($(VectorizationBase.pick_vector(M,T)), zero($T))
        end
    else
        pre_quote = quote qf = zero($T) end
    end
    if partial && track_σ
        push!(pre_quote.args, :(σ = ∂canonicalize_Σ(σin)))
    else
        push!(pre_quote.args, :(σ = canonicalize_Σ(σin)))
    end
    return_expr = quote end
    if yisvec # not allowed to be nothing, must be bool
        yexpr = :(y[i])
    else
        yexpr = :y
    end
    if μisvec === nothing
        δexpr = yexpr
        #if not, then it is a bool
    elseif μisvec # == true
        δexpr = :($yexpr - μ[i])
    else # μisvec == false
        δexpr = :($yexpr - μ)
    end
    # add target
    # @show M, T, yisvec, μisvec, σisvec
    # @show track_y, track_μ, track_σ
    # @show inity, initμ, initσ
    # @show partial, calclogdet
    if σisvec === nothing
        # push!(pre_quote.args, :(σ⁻¹ = One() ))
        loop_expr = quote
            δ = $δexpr
            δσ⁻² = δ
            qf = vfmadd(δ, δ, qf)
        end
    elseif σisvec# == true
        push!(pre_quote.args, :(σ⁻² = precision(σ)))
        loop_expr = quote
            δ = $δexpr
            σ⁻²ᵢ = σ⁻²[i]
            δσ⁻² = δ * σ⁻²ᵢ
            qf = vfmadd(δσ⁻², δ, qf)
        end
        if calclogdet && track_σ
            push!(pre_quote.args, :(nld = vbroadcast(Vec{$W,$T},zero($T))))
            push!(pre_quote.args, :(ldv = nlogdet(σ)))
            push!(loop_expr, :(nld = vadd(nld, ldv[i])))
        end
    elseif partial && (track_y || track_μ) #σisvec == false, and we're calculating ∂y or ∂μ
        push!(pre_quote.args, :(σ⁻² = precision(σ)))
        if track_σ
            loop_expr = quote
                δ = $δexpr
                δσ⁻² = δ * σ⁻²
                δ²σ⁻² = δσ⁻² * δ
                qf = vadd(δ²σ⁻², qf)
            end
        else
            loop_expr = quote
                δ = $δexpr
                δσ⁻² = δ * σ⁻²
                qf = vfmadd(δσ⁻², δ, qf)
            end
        end
    else #σisvec == false
        # we do not need to keep track of δσ⁻², so we multiply at the end of the reduction.
        loop_expr = quote
            δ = $δexpr
            qf = vfmadd(δ, δ, qf)
        end
        push!(return_expr.args, :(qf = vmul(qf, precision(σ))))
    end
    if !σisvec && calclogdet && track_σ
        if M == 1
            push!(return_expr.args, :(nld = nlogdet(σ)))
        else
            push!(return_expr.args, :(nld = vmul($(T(M)), nlogdet(σ))))
        end
    end
    if partial
        if track_y
            if yisvec
                if inity
                    push!(loop_expr.args, :(∂y[i]  = -δσ⁻²))
                else
                    push!(loop_expr.args, :(∂y[i] -=  δσ⁻²))
                end
            else
                push!(pre_quote.args, :(∂ys  = zero($T)))
                push!(loop_expr.args, :(∂ys -= δσ⁻²))
                push!(return_expr.args, inity ? :(∂y[] = ∂ys) : :(∂y[] += ∂ys))
            end
        end
        if track_μ
            if μisvec# == true
                if initμ
                    push!(loop_expr.args, :(∂μ[i]  = δσ⁻²))
                else
                    push!(loop_expr.args, :(∂μ[i] += δσ⁻²))
                end
            else#if μisvec == false
                if yisvec
                    push!(pre_quote.args, :(∂μs = zero($T)))
                    push!(loop_expr.args, :(∂μs += δσ⁻²))
                else # both yisvec && μisvec
                    push!(return_expr.args, :(∂μs = Base.FastMath.sub_fast(∂y)))
                end
                push!(return_expr.args, initμ ? :(∂μ[] = ∂μs) : :(∂μ[] += ∂μs))
            end
        end
        if track_σ
            eq = initσ ? :(=) : :(+=)
            if σisvec# == true
                push!(pre_quote.args, :(∂k∂σ  = ∂k∂i(σ)))
                if calclogdet
                    push!(pre_quote.args, :( ∂ld∂σ = ∂ld∂i(σ) ) )
                    set∂σ = :(∂k∂σ[i] * δσ⁻² - ∂ld∂σ[i])
                else
                    set∂σ = :(∂k∂σ[i] * δσ⁻²)
                end
                push!(loop_expr.args, Expr(eq, :(∂σ[i]), set∂σ))
            else#if σisvec == false
                if calclogdet
                    if yisvec
                        push!(return_expr.args, Expr(eq, :(∂σ[]), :(∂k∂i(σ) * vsum(qf) - $(T(M)) * ∂ld∂i(σ))))
                    else
                        push!(return_expr.args, Expr(eq, :(∂σ[]), :(∂k∂i(σ) * vsum(qf) - ∂ld∂i(σ))))
                    end
                else
                    push!(return_expr.args, Expr(eq, :(∂σ[]), :(∂k∂i(σ) * vsum(qf))))
                end
            end
        end
    end
    if calclogdet && track_σ
        if yisvec && !σisvec
            retex = :(tadd(vmul($(T(-0.5)), qf), nld))
        else # either both isvec, or neither
            retex = :(vfmadd($(T(-0.5)), qf, nld))
        end
    else
        retex = :(vmul($(T(-0.5)), qf))
    end
    q = if yisvec
        loop_expr = quote
            LoopVectorization.@vvectorize $T 4 for i ∈ 1:$M
                $loop_expr
            end
        end
        quote
            # $(Expr(:meta,:inline))
            $pre_quote
            $(macroexpand(ProbabilityDistributions, loop_expr))
            $return_expr
            $retex
        end
    else
        quote
            # $(Expr(:meta,:inline))
            @fastmath begin
                $pre_quote
                $loop_expr
                $return_expr
            end
            $retex
        end
    end
    if (yisvec && !σisvec) || (4VectorizationBase.pick_vector_width(T) <= M)
        pushfirst!(q.args, Expr(:meta, :inline)) # inline Target...
    end
    simplify_expr(q)
end
function alloc_univariate_normal_quote(M, S, @nospecialize(T), (track_y, track_μ, track_σ), sp = true)
    N = length(S.parameters)
    X = Vector{Int}(undef, N)
    X[1] = 1
    for n in 2:N
        X[n] = X[n-1] * (S.parameters[n-1])::Int
    end
    P = Tuple{X...}
    q = quote end
    if track_y
        if yisvec
            if sp
                push!(q.args, :((sp,∂y) = PtrArray{$S,$T,$N,$P}(sp)))
            else
                push!(q.args, :(∂y = FixedSizeArray{$S,$T,$N,$P}(undef) ))
            end
        else
            push!(q.args, :(∂y = Ref{$T}()))
        end
    end
    if track_μ
        if μisvec == true
            if sp
                push!(q.args, :((sp,∂μ) = PtrArray{$S,$T,$N,$P}(sp)))
            else
                push!(q.args, :(∂μ = FixedSizeArray{$S,$T,$N,$P}(undef) ))
            end
        elseif μisvec == false
            push!(q.args, :(∂μ = Ref{$T}()))
        end
    end
    if track_σ
        if σisvec == true
            if sp
                push!(q.args, :((sp,∂σ) = PtrArray{$S,$T,$N,$P}(sp)))
            else
                push!(q.args, :(∂σ = FixedSizeArray{$S,$T,$N,$P}(undef) ))
            end
        elseif σisvec == false
            push!(q.args, :(∂σ = Ref{$T}()))
        end
    end
    q
end

@noinline function univariate_normal_length(SV::Core.SimpleVector, RV::Core.SimpleVector)
    PaddedMatrices.isdense(SV, RV) || throw("Non-unit stride not yet supported.")
    # because isdense passed, last(RV) must equal prod(SV[1:N-1])
    last(RV)::Int * last(SV)::Int
end

@inline Normal(y::T) where {T <: Real} = Base.FastMath.mul_fast(T(-0.5), Base.FastMath.abs2_fast(y))
@inline Normal(::Val{(true,)}, y::Real) = Normal(y)
@inline Normal(::Val{(false,)}, y::T) where {T <: Real} = zero(T)

@generated function Normal(::Val{(true,)}, y::AbstractFixedSizeArray{S,T,N,X,L}) where {S,T,N,X,L}
    nrows = first(S.parameters)::Int
    if PaddedMatrices.isdense(S.parameters, X.parameters)
        M = N == 1 ? nrows : L
        W = VectorizationBase.pick_vector_width(M, T)
        q = quote
            out = vbroadcast(Vec{$W,$T}, zero($T))
            @vvectorize $T 4 for m ∈ 1:$M
                yₘ = y[m]
                out = vmuladd(yₘ, yₘ, out)
            end
            vmul($(T(-0.5)), out)
        end
        (4VectorizationBase.pick_vector_width(T) <= M) && pushfirst!(q.args, Expr(:meta, :inline))
        q
    else
        P = (X.parameters[2])::Int
        W = VectorizationBase.pick_vector_width(nrows, T)
        return quote
            out = vbroadcast(Vec{$W,$T}, zero($T))
            ind = 0
            Base.Cartesian.@nloops $(N-1) i j -> 1:size(y,j+1) begin
                @vvectorize $T 4 for i_0 ∈ 1:$nrows
                    yᵢ = y[ind + i_0]
                    out = vmuladd(yᵢ, yᵢ, out)
                end
                ind += $P
            end
            vmul($(T(-0.5)), out)
        end
    end
end
@generated function ∂Normal!(∂y::AbstractFixedSizeArray{S,T,N,X,L}, y::AbstractFixedSizeArray{S,T,N,X,L}) where {S,T,N,X,L}
    nrows = first(S.parameters)::Int
    if PaddedMatrices.isdense(S.parameters, X.parameters)
        M = N == 1 ? nrows : L
        W = VectorizationBase.pick_vector_width(M, T)
        q = quote
            out = vbroadcast(Vec{$W,$T}, zero($T))
            @vvectorize $T 4 for m ∈ 1:$M
                yₘ = y[m]
                ∂y[m] = -yₘ
                out = vmuladd(yₘ, yₘ, out)
            end
            vmul($(T(-0.5)), out)
        end
        (4VectorizationBase.pick_vector_width(T) <= M) && pushfirst!(q.args, Expr(:meta, :inline))
        return q
    else
        P = (X.parameters[2])::Int
        W = VectorizationBase.pick_vector_width(nrows, T)
        return quote
            out = vbroadcast(Vec{$W,$T}, zero($T))
            ind = 0
            Base.Cartesian.@nloops $(N-1) i j -> 1:size(y,j+1) begin
                @vvectorize $T 4 for i_0 ∈ 1:$nrows
                    yᵢ = y[ind + i_0]
                    ∂y[ind + i_0] = -yᵢ
                    out = vmuladd(yᵢ, yᵢ, out)
                end
                ind += $P
            end
            vmul($(T(-0.5)), out)
        end
    end
end

@inline function ∂Normal!(∂y::U, y::T) where {T <: Real, U}
    t = Base.FastMath.mul_fast(T(-0.5), Base.FastMath.abs2_fast(y))
    if isinitialized(U)
        ∂y[] = Base.FastMath.sub_fast(∂y[], y)
    else
        ∂y[] = Base.FastMath.sub_fast(y)
    end
    t
end
@inline ∂Normal!(::Nothing, y::T) where {T <: Real} = zero(T)
@generated function Normal_kernel(x::AbstractFixedSizeArray{S,T,N,X,L}) where {S,T,N,X,L}
    P = univariate_normal_length(S.parameters, X.parameters)
    V = VectorizationBase.pick_vector(P,T)
    q = quote
        # $(Expr(:meta,:inline))
        target = vbroadcast($V, zero($T))
        @vvectorize $T 4 for p ∈ 1:$P
            xₚ = x[p]
            target = vmuladd(xₚ, xₚ, target)
        end
        vmul(target, vbroadcast($V, $(T(-0.5))))
    end
    (4VectorizationBase.pick_vector_width(T) <= P) && pushfirst!(q.args, Expr(:meta, :inline))
    q
end
@inline function Normal_kernel(::Val{track}, x::AbstractFixedSizeArray{S,T}) where {S,T,track}
    first(track) ? Normal_kernel(x) : Zero()
end

# @eval loop to specify normals
for yisvec ∈ (true,false)
    if yisvec
        args1 = [:(y::AbstractFixedSizeArray{S,T,N,R,L})]
        M = :(univariate_normal_length(S.parameters, R.parameters))
        whereparams = [:S,:N,:R,:L]
        σisvec = :(σin <: AbstractFixedSizeArray)
    else
        args1 = [:(y::Union{T,RealFloat{<:Any,T}})]
        M = 1
        whereparams = Symbol[]
        σisvec = false
    end
    ∂args1 = [:(∂y::∂YN)]; ∂whereparams1 = [:∂YN, :∂ΣN]
    for μ ∈ (false,true)
        if μ
            args2 = push!(copy(args1), yisvec ? :(μ::Union{T,Int,<:AbstractFixedSizeArray{S,T,N,R,L}}) : :(μ::Union{T,Int}))
            track_μ = :(track[2])
            initμ = :(!isinitialized(∂MN))
            ∂args2 = push!(copy(∂args1), :(∂μ::∂MN));
            ∂whereparams2 = push!(copy(∂whereparams1), :∂MN)
            μisvec = :(μ <: AbstractFixedSizeArray)
            ∂track_μ = :(∂MN !== Nothing)
        else
            args2 = args1
            ∂args2 = ∂args1
            ∂whereparams2 = ∂whereparams1
            track_μ = false
            μisvec = nothing
            initμ = false
            ∂track_μ = false
        end
        ∂args3 = push!(copy(∂args2), :(∂σ::∂ΣN))
        args3 = push!(copy(args2), yisvec ? :(σin::Union{T,Int,RealFloat{<:Any,T},Precision,<:AbstractFixedSizeArray{S,T,N,R,L}}) : :(σin::Union{T,Int,RealFloat{<:Any,T},Precision{T}}))
        # @show μ, args3, μisvec
        for calclogdet ∈ (true,false)
            n = calclogdet ? :Normal : :Normal_kernel
            ∂n = calclogdet ? :∂Normal! : :∂Normal_kernel!
            @eval @generated function $n(::Val{track}, $(args3...)) where {track, $(whereparams...), T <: Union{Float32,Float64}}
                univariate_normal_quote(
                    $M, T, $yisvec, $μisvec, $σisvec,
                    (first(track), $track_μ, last(track)),
                    (false, false, false),
                    false, $calclogdet
                )
            end
            @eval @generated function $n($(args3...)) where {$(whereparams...), T <: Union{Float32,Float64}}
                univariate_normal_quote(
                    $M, T, $yisvec, $μisvec, $σisvec,
                    ((true, $μ, true)),
                    (false, false, false),
                    false, $calclogdet
                )
            end

            @eval @generated function $∂n($(∂args3...), $(args3...)) where {$(∂whereparams2...), $(whereparams...), T <: Union{Float32,Float64}}
                univariate_normal_quote(
                    $M, T, $yisvec, $μisvec, $σisvec,
                    (∂YN !== Nothing, $∂track_μ, ∂ΣN !== Nothing),
                    (!isinitialized(∂y), $initμ, !isinitialized(∂σ)),
                    true, $calclogdet
                )
            end
        end
    end
end


text:
# This file includes the buffer test set

@testset "BufferTestSet" begin
    @info "Running BufferTestSet ..."

    # Simple Buffer construction
    buf = Buffer(5)
    @test eltype(buf) == Float64
    @test length(buf) == 5
    @test size(buf) == (5,)
    @test mode(buf) == Cyclic
    @test buf.index == 1
    @test isempty(buf)
    @test buf.state == :empty
    @test size(buf) == (5,)
    @test isa(buf, AbstractArray)

    # Buffer data length
    buf = Buffer(5)
    @test datalength(buf) == 5
    buf = Buffer(3, 10)
    @test datalength(buf) == 10

    # Writing values into Buffers
    buf = Buffer(5)
    write!(buf, 1.)
    @test !isempty(buf)
    @test !isfull(buf)
    @test buf.index == 2
    @test buf.state == :nonempty

    # Reading from buffers
    val = read(buf)
    @test val == 1.
    @test buf.index == 2
    @test !isempty(buf)

    # More on buffer construction
    buf = Buffer{Fifo}(Float64, 2, 5)
    buf = Buffer{Fifo}(Float64, 5)
    buf = Buffer{Fifo}(5)
    buf = Buffer{Normal}(5)
    buf = Buffer(5)

    # # Filling buffers
    # buf = Buffer{Cyclic}(5)
    # fill!(buf, 1.)
    # @test outbuf(buf) == ones(5)
    # buf = Buffer{Normal}(2,5)
    # fill!(buf, [1, 1])
    # @test buf.data == ones(2, 5)

    # Writing into Buffers with different modes
    for bufmode in [Normal, Lifo, Fifo]
        buf = Buffer{bufmode}(2, 5)
        for item in 1 : 5
            write!(buf, [item, item])
        end
        @test outbuf(buf) == [5. 4. 3. 2. 1.; 5. 4. 3. 2. 1.]
        @test isfull(buf)
        @test buf.index == 6
        @test_throws Exception write!(buf, [1., 2.])  # When full, data cannot be written into buffers.
    end

    buf = Buffer{Cyclic}(2, 5)
    for item in 1 : 5
        write!(buf, [item, item])
    end
    @test outbuf(buf) == [5. 4. 3. 2. 1.; 5. 4. 3. 2. 1.]
    @test isfull(buf)
    @test buf.index == 6
    temp = outbuf(buf)
    write!(buf, [6., 6.])  # When full, data can be written into Cyclic buffers.
    @test outbuf(buf) == hcat([6., 6.], temp[:, 1:end-1])

    # Reading from Buffers with different modes
    for bufmode in [Normal, Cyclic]
        buf = Buffer{bufmode}(5)
        foreach(item -> write!(buf, item), 1 : 5)
        for i = 1 : 5
            @test read(buf) == 5
        end
        @test !isempty(buf)
    end

    buf = Buffer{Fifo}(5)
    foreach(item -> write!(buf, item), 1 : 5)
    for i = 1 : 5
        @test read(buf) == i
    end
    @test isempty(buf)
    @test_throws Exception read(buf)  # When buffer is empty, no more reads.

    buf = Buffer{Lifo}(5)
    foreach(item -> write!(buf, item), 1 : 5)
    vals = collect(5:-1:1)
    for i = 1 : 5
        @test read(buf) == vals[i]
    end
    @test isempty(buf)
    @test_throws Exception read(buf)  # When buffer is empty, no more reads.

    @info "Done BufferTestSet."
end  # testset