text:
# Routines for an unstructured mesh that is contained in a hierarchical, rectangularly
# partitioned mesh data structure
mutable struct HierarchicalRectangularlyPartitionedMesh
    name::String
    rect::Rectangle_2D
    mesh::Ref{UnstructuredMesh_2D}
    parent::Ref{HierarchicalRectangularlyPartitionedMesh}
    children::Vector{Ref{HierarchicalRectangularlyPartitionedMesh}}
end

function HierarchicalRectangularlyPartitionedMesh(;
        name::String = "DefaultName",
        rect::Rectangle_2D = Rectangle_2D(),
        mesh::Ref{UnstructuredMesh_2D} = Ref{UnstructuredMesh_2D}(),
        parent::Ref{HierarchicalRectangularlyPartitionedMesh}
              = Ref{HierarchicalRectangularlyPartitionedMesh}(),
        children::Vector{Ref{HierarchicalRectangularlyPartitionedMesh}}
                       = Ref{HierarchicalRectangularlyPartitionedMesh}[]
    )
    this = HierarchicalRectangularlyPartitionedMesh(name, rect, mesh, parent, children)
    # If this mesh has a parent, add this mesh to the parent's children vector
    if isassigned(parent)
        push!(parent[].children, Ref(this))
    end
    return this
end

Base.broadcastable(HRPM::HierarchicalRectangularlyPartitionedMesh) = Ref(HRPM)

# Add the boundary edges to each mesh in the HRPM
# Type-stable
function add_boundary_edges!(HRPM::HierarchicalRectangularlyPartitionedMesh)
    apply_function_recursively_to_HRPM_meshes!(add_boundary_edges, HRPM)
end

# Add the connectivity to each mesh in the HRPM
# Type-stable
function add_connectivity!(HRPM::HierarchicalRectangularlyPartitionedMesh)
    apply_function_recursively_to_HRPM_meshes!(add_connectivity, HRPM)
end

# Add the edges to each mesh in the HRPM
# Type-stable
function add_edges!(HRPM::HierarchicalRectangularlyPartitionedMesh)
    apply_function_recursively_to_HRPM_meshes!(add_edges, HRPM)
end

# Add every field to each mesh in the HRPM
# Type-stable
function add_everything!(HRPM::HierarchicalRectangularlyPartitionedMesh)
    apply_function_recursively_to_HRPM_meshes!(add_everything, HRPM)
end

# Add materialized edges to each mesh in the HRPM
# Type-stable
function add_materialized_edges!(HRPM::HierarchicalRectangularlyPartitionedMesh)
    apply_function_recursively_to_HRPM_meshes!(add_materialized_edges, HRPM)
end

# Add materialized faces to each mesh in the HRPM
# Type-stable
function add_materialized_faces!(HRPM::HierarchicalRectangularlyPartitionedMesh)
    apply_function_recursively_to_HRPM_meshes!(add_materialized_faces, HRPM)
end

# Apply a function, f, to each of the meshes in the HRPM
function apply_function_recursively_to_HRPM_meshes!(f::Function,
                                                    HRPM::HierarchicalRectangularlyPartitionedMesh)
    nchildren = length(HRPM.children)
    if isassigned(HRPM.mesh)
        HRPM.mesh[] = f(HRPM.mesh[])
    elseif 0 < nchildren
        for ichild = 1:nchildren
            apply_function_recursively_to_HRPM_meshes!(f, HRPM.children[ichild][])
        end
    end
    return nothing
end

# Apply a mutating function, f!, to each of the meshes in the HRPM
function apply_mutating_function_recursively_to_HRPM_meshes!(f::Function,
                                                    HRPM::HierarchicalRectangularlyPartitionedMesh)
    nchildren = length(HRPM.children)
    if isassigned(HRPM.mesh)
        f(HRPM.mesh[])
    elseif 0 < nchildren
        for ichild = 1:nchildren
            apply_mutating_function_recursively_to_HRPM_meshes!(f, HRPM.children[ichild][])
        end
    end
    return nothing
end

# Return the axis-aligned bounding box of the HRPM
function boundingbox(HRPM::HierarchicalRectangularlyPartitionedMesh)
    if HRPM.rect !== Rectangle_2D()
        return HRPM.rect
    elseif isassigned(HRPM.mesh)
        HRPM.rect = boundingbox(HRPM.mesh[], boundary_shape="Rectangle")
        return HRPM.rect
    elseif 0 < length(HRPM.children)
        HRPM.rect = mapreduce(x->boundingbox(x[]), union, HRPM.children)
        return HRPM.rect
    else
        @error "Something went wrong"
        return Rectangle_2D()
    end
end

# # Fill a statically sized, mutable vector (coord) with the necessary indices to navigate from the
# # root HRPM, through the children, to the base mesh, and return the face ID to which the point p
# # may be found in.
# # Example:
# # For an HRPM with 4 levels:
# # [1, 2, 1, 16]
# # denotes HRPM.children[1][].children[2][].children[1][].mesh[].faces[16] contains p
# # If the face is found, return true. Otherwise, return false
# function find_face!(p::Point_2D, coord::MVector{N, UInt32},
#                     HRPM::HierarchicalRectangularlyPartitionedMesh)
#     in_rect = p ∈  HRPM.rect
#     if !in_rect
#         return false
#     elseif in_rect && (0 < length(HRPM.children))
#         for (i, child) in enumerate(HRPM.children)
#             in_child = find_face!(p, coord, child[])
#             if in_child
#                 coord[findfirst(x->x==0, coord)] = U(i)
#                 return true
#             end
#         end
#         return false
#     elseif in_rect && isassigned(HRPM.mesh)
#         face = find_face(p, HRPM.mesh[]::UnstructuredMesh_2D{F, U})
#         coord[findfirst(x->x==0, coord)] = face
#         reverse!(coord)
#         return face == 0 ? false : true
#     end
#     return false
# end
#
# # Get the intersection algorithm that will be used for l ∩ HRPM
# # Not type-stable
# function get_intersection_algorithm(HRPM::HierarchicalRectangularlyPartitionedMesh)
#     if isassigned(HRPM.mesh)
#         return get_intersection_algorithm(HRPM.mesh[])
#     else
#         return get_intersection_algorithm(HRPM.children[1][])
#     end
# end
#
# # Check if the HRPM's mesh faces are materialized
# # Not type-stable
# function has_materialized_faces(HRPM::HierarchicalRectangularlyPartitionedMesh)
#     if isassigned(HRPM.mesh)
#         if length(HRPM.mesh[].materialized_faces) !== 0
#             return true
#         else
#             return false
#         end
#     else
#         return has_materialized_faces(HRPM.children[1][])
#     end
# end
#
# # Intersection a line with the HRPM. Returns a vector of points, ordered by distance from
# # the line's start point
# # Not type-stable
# function intersect(l::LineSegment_2D{F},
#                    HRPM::HierarchicalRectangularlyPartitionedMesh{F, U}) where {F <: AbstractFloat,
#                                                                                U <: Unsigned}
#     # An array to hold all of the intersection points
#     intersection_points = Point_2D{F}[]
#     nchildren = length(HRPM.children)
#     if 0 < (l ∩ HRPM.rect)[1]
#         if isassigned(HRPM.mesh)
#             append!(intersection_points, l ∩ HRPM.mesh[]::UnstructuredMesh_2D{F,U})
#         elseif 0 < nchildren
#             for ichild = 1:nchildren
#                 append!(intersection_points,
#                         l ∩ HRPM.children[ichild][]::HierarchicalRectangularlyPartitionedMesh{F, U})
#             end
#             return sort_points(l.points[1], intersection_points)
#         end
#     end
#     return intersection_points
# end
#
# Return the height of the HRPM (number of edges between this node and the leaf)
#         A        H = 2
#        / \
#       /   \
#      B     C     H = 1
#     / \   / \
#    D   E F   G   H = 0
# Not type-stable
function node_height(HRPM::HierarchicalRectangularlyPartitionedMesh)
    if length(HRPM.children) === 0
        return 0
    elseif 0 < length(HRPM.children)
        return node_height(HRPM.children[1][]) + 1
    else
        @error "Something went wrong"
        return -100
    end
end

# Get the level (distance from current node to root + 1) of the HRPM
#         A        L = 1
#        / \
#       /   \
#      B     C     L = 2
#     / \   / \
#    D   E F   G   L = 3
# Not type-stable
function node_level(HRPM::HierarchicalRectangularlyPartitionedMesh; current_level::Int64=1)
    if isassigned(HRPM.parent)
        return node_level(HRPM.parent[]; current_level = current_level + 1)
    else
        return current_level
    end
end

# Partition a mesh into an HRPM based upon the names of its face sets.
# Must contain face sets of the form "GRID_LN_X_Y" where N,X,Y are integers
# N is the level of the node and X,Y are indices of the mesh's location in a rectangular
# grid
function partition_rectangularly(mesh::UnstructuredMesh_2D)
    @info "Converting $(typeof(mesh)) into HierarchicalRectangularlyPartitionedMesh"
    # Extract set names, grid names, and max level
    set_names, grid_names, max_level = process_partition_rectangularly_input(mesh)

    # Create a tree to store grid relationships.
    root = create_HRPM_tree(mesh, grid_names, max_level)

    # Construct the leaf meshes
    leaf_meshes = create_HRPM_leaf_meshes(mesh, grid_names, max_level)

    # Construct the mesh hierarchy
    HRPM = create_HRPM(root, leaf_meshes)
    return HRPM
end

# How to display the HRPM in the REPL
function Base.show(io::IO, HRPM::HierarchicalRectangularlyPartitionedMesh; relative_offset::Int64 = 0)
    println(io, "HierarchicalRectangularlyPartitionedMesh")
    println(io, "  ├─ Name      : $(HRPM.name)")
    println(io, "  ├─ Bounding Box : $(HRPM.rect)")
    println(io, "  ├─ Mesh      : $(isassigned(HRPM.mesh))")
    println(io, "  ├─ Children  : $(length(HRPM.children))")
    if isassigned(HRPM.parent)
        println(io, "  └─ Parent    : $(HRPM.parent[].name)")
    else
        println(io, "  └─ Parent    : None")
    end
end

# Functions used in partition_rectangularly
# ------------------------------------------------------------------------------------------------------
function attach_HRPM_children!(HRPM::HierarchicalRectangularlyPartitionedMesh,
                               tree::Tree,
                               leaf_meshes::Vector{<:UnstructuredMesh_2D})
    for child in tree.children
        name = child[].data
        child_mesh = HierarchicalRectangularlyPartitionedMesh(name = name,
                                                              parent = Ref(HRPM) )
        for leaf_mesh in leaf_meshes
            if name == leaf_mesh.name
                child_mesh.mesh[] = leaf_mesh
            end
        end
        attach_HRPM_children!(child_mesh, child[], leaf_meshes)
    end
    return nothing
end

# Construct the HRPM
function create_HRPM(tree::Tree, leaf_meshes::Vector{UnstructuredMesh_2D})
    # Construct the HRPM from the top down
    root = HierarchicalRectangularlyPartitionedMesh( name = tree.data )
    attach_HRPM_children!(root, tree, leaf_meshes)
    # Add the rectangles
    boundingbox(root)
    return root
end

# Construct the leaf meshes
function create_HRPM_leaf_meshes(mesh::UnstructuredMesh_2D, grid_names::Vector{String}, max_level::Int64)
    # Generate the leaf meshes (The smallest spatially)
    leaf_meshes = UnstructuredMesh_2D[]
    leaf_names = String[]
    for name in grid_names
        level = parse(Int64, name[7])
        if level == max_level
            push!(leaf_names, name)
        end
    end
    for name in leaf_names
        push!(leaf_meshes, submesh(name, mesh))
    end
    # remove grid levels
    for leaf_mesh in leaf_meshes
        for name in keys(leaf_mesh.face_sets)
            if occursin("GRID_", uppercase(name))
                delete!(leaf_mesh.face_sets, name)
            end
        end
    end
    return leaf_meshes
end

# Create a tree to store grid relationships.
function create_HRPM_tree(mesh::UnstructuredMesh_2D, grid_names::Vector{String}, max_level::Int64)
    root = Tree( data = mesh.name )
    current_nodes = Tree[]
    next_nodes = Tree[]
    old_grid_names = copy(grid_names)
    new_grid_names = copy(grid_names)
    # Do first level
    for grid_name in old_grid_names
        grid_level = parse(Int64, grid_name[7])
        if grid_level === 1
            # Add to appropriate node (root)
            push!(next_nodes, Tree(data = grid_name; parent=Ref(root)))
            filter!(x->x ≠ grid_name, new_grid_names)
        end
    end
    # Do all other levels:
    for level in 2:max_level
        old_grid_names = copy(new_grid_names)
        current_nodes = next_nodes
        next_nodes = []
        for grid_name in old_grid_names
            grid_level = parse(Int64, grid_name[7])
            if grid_level == level
                # find the parent for this grid
                grid_faces = mesh.face_sets[grid_name]
                for node in current_nodes
                    node_faces = mesh.face_sets[node.data]
                    if grid_faces ⊆ node_faces
                        push!(next_nodes, Tree(data = grid_name, parent=Ref(node)))
                        filter!(x->x ≠ grid_name, new_grid_names)
                        break
                    end
                end
            end
        end
    end
    return root
end

# Us this the last child in the parent's list of children?
# offset determines if the nth-parent is the last child
function is_last_child(HRPM::HierarchicalRectangularlyPartitionedMesh; relative_offset::Int64=0)
    if !isassigned(HRPM.parent)
        return true
    end
    if relative_offset > 0
        return is_last_child(HRPM.parent[]; relative_offset=relative_offset-1)
    else
        nsiblings = length(HRPM.parent[].children) - 1
        return (HRPM.parent[].children[nsiblings + 1][] == HRPM)
    end
end

# Extract set names, grid names, and max level
function process_partition_rectangularly_input(mesh::UnstructuredMesh_2D)
    set_names = collect(keys(mesh.face_sets))
    grid_names = copy(set_names)
    for set_name in set_names
        if !occursin("GRID_", uppercase(set_name))
            filter!(x->x ≠ set_name, grid_names)
        end
    end
    if length(grid_names) === 0
        @error "No grid face sets in mesh"
    end

    # Get the number of grid levels
    max_level = 0
    for grid_name in grid_names
        level = parse(Int64, grid_name[7])
        if max_level < level
            max_level = level
        end
    end

    return set_names, grid_names, max_level
end

# Plot
# -------------------------------------------------------------------------------------------------
if enable_visualization
    function linesegments!(HRPM::HierarchicalRectangularlyPartitionedMesh)
        apply_mutating_function_recursively_to_HRPM_meshes!(linesegments!, HRPM)
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct EvictionPolicy <: SwaggerModel

    function EvictionPolicy(;)
        o = new()
        o
    end
end # type EvictionPolicy

const _property_map_EvictionPolicy = Dict{Symbol,Symbol}()
const _property_types_EvictionPolicy = Dict{Symbol,String}()
Base.propertynames(::Type{ EvictionPolicy }) = collect(keys(_property_map_EvictionPolicy))
Swagger.property_type(::Type{ EvictionPolicy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_EvictionPolicy[name]))}
Swagger.field_name(::Type{ EvictionPolicy }, property_name::Symbol) =  _property_map_EvictionPolicy[property_name]

function check_required(o::EvictionPolicy)
    true
end

function validate_property(::Type{ EvictionPolicy }, name::Symbol, val)
end

ext: jl


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements sample function
# 3. Exports Type

# using Random

"""
    LinSysVecRowRandCyclic <: LinSysVecRowSelect

A mutable structure with a field to store a cycling order. Randomly specifies a cycling
order the equations of a linear system. Once this ordering is exhausted by the solver,
a new random ordering is specified. This process is repeated

# Fields
- `order::Union{Vector{Int64},Nothing}`

Calling `LinSysVecRowOneRandCyclic()` defaults to setting `order` to `nothing`. The `sample`
function will handle the re-initialization of the fields once the system is provided.
"""
mutable struct LinSysVecRowRandCyclic <: LinSysVecRowSelect
    order::Union{Vector{Int64},Nothing}
end
LinSysVecRowRandCyclic() = LinSysVecRowRandCyclic(nothing)

# Common sample interface for linear systems
function sample(
    type::LinSysVecRowRandCyclic,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)

    mod_ind = mod(iter, 1:length(b))
    if mod_ind == 1
        type.order = randperm(length(b))
    end

    eqn_ind = type.order[mod_ind]
    return A[eqn_ind, :], b[eqn_ind]
end

#export LinSysVecRowRandCyclic


text:
function TRL(vks,uks,Tmat,k,
             pbits,nbits,jocc_p,jocc_n,
             SPEs,ppinfo,nninfo,bis,bfs,block_tasks,
             p_NiNfs,n_NiNfs,Mps,delMs,Vpn,
              eval_jj,oPP,oNN,oPNu,oPNd,Jidxs,
             tdims,num_ev,num_history,lm,ls,en,tol,to,doubleLanczos=false)
    mdim = tdims[end]
    TF=[false]

    lnJ = ls
    Jvs = [zeros(Float64,mdim) for i=1:lnJ]
    Jvret = [zeros(Float64,mdim)]
    Jmat = zeros(Float64,lnJ,lnJ)
    Jvs[1] .= vks[1]
    Jtol = 1.e-6
    JTF = [false];beta_J = 1.0
    if doubleLanczos
        Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,
                 pbits,nbits,tdims,eval_jj,
                 Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to)
        vks[1] .= Jvs[1]
        for i=1:lnJ;Jvs[i] .=0.0;end
    end
    elit = 1
    while TF[1]==false
        for it = k:lm-1
            vk =vks[it]; vkp1 =vks[it+1]
            @timeit to "operate H" begin
                operate_H!(vk,vkp1,
                           pbits,nbits,
                           jocc_p,jocc_n,SPEs,
                           ppinfo,nninfo,
                           tdims,bis,bfs,block_tasks,
                           p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to)
            end
            talpha =  dot(vk,vkp1)
            Tmat[it,it] = talpha
            diagonalize_T!(it,num_ev,Tmat,en,num_history,TF,tol)
            #print_vec("En TR(d)L $it ",en[1])
            if TF[1];elit=it;break;end
            axpy!(-talpha,vk,vkp1)
            svks = @views vks[1:it-1]
            @timeit to "ReORTH" ReORTH(it,vkp1,svks)
            tbeta = sqrt(dot(vkp1,vkp1))
            tmp = 1.0/tbeta
            vkp1 .*= tmp
            Tmat[it+1,it] = tbeta; Tmat[it,it+1] = tbeta
            if doubleLanczos
                if tbeta < Jtol;TF[1]=true;elit=it;break;end

                @timeit to "JJ lanczos" begin
                    for i =2:length(Jvs);Jvs[i] .=0.0;end
                    Jmat .= 0.0;JTF[1] = false
                    Jvs[1] .= vkp1
                    Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,
                             pbits,nbits,tdims,eval_jj,
                             Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to)
                    vkp1 .= Jvs[1]
                    if TF[1];elit=it;break;end
                end
            end
        end
        if TF[1] == false
            @timeit to "Restart" begin
                ThickRestart(vks,uks,Tmat,lm,ls)
            end
        end
        k = ls+1
    end
    return elit
end

function Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,
                  pbits,nbits,tdims,eval_jj,
                  Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to)
    mdim = tdims[end]
    lnJ = length(Jvs)
    eljit=1; k = 1; inow=k
    while JTF[1] == false
        for it = k:lnJ-1
            inow = it
            vk = Jvs[it]; vkp1 = Jvs[it+1]
            operate_J!(vk,vkp1,pbits,nbits,tdims,
                       Jidxs,oPP,oNN,oPNu,oPNd,beta_J)
            axpy!(eval_jj,vk,vkp1)## vkp1 .+= eval_jj .* vk
            Jmat[it,it] = dot(vk,vkp1)
            teval = eigvals(@views Jmat[1:it,1:it])[1]
            if abs(teval-eval_jj) < Jtol
                eljit=it;JTF[1] = true;break
            end
            if JTF[1];eljit=it;break;end
            ReORTH(it,vkp1,Jvs)
            beta = sqrt(dot(vkp1,vkp1))
            if beta < 1.e-4;eljit=it;TF[1]=true;JTF[1]=true;break;end
            vkp1 .*= 1.0/beta
            Jmat[it+1,it] = beta; Jmat[it,it+1] = beta
            eljit = it
        end
        if JTF[1]==false
            ThickRestart_J(Jvs,Jvret,Jmat,eljit+1,1,eval_jj,Jtol)
            k=2
        end
    end
    if inow > k
        ThickRestart_J(Jvs,Jvret,Jmat,eljit+1,1,eval_jj,Jtol)
    end
    return nothing
end

function operate_H!(wf,twf,
                    pbits,nbits,
                    jocc_p,jocc_n,
                    SPEs,ppinfo,nninfo,
                    tdims,bis,bfs,block_tasks,
                    p_NiNfs,n_NiNfs,Vpn,
                    Mps,delMs,
                    to=nothing)
    #lblock = length(pbits)
    #@inbounds @qthreads for bi=1:lblock
    @inbounds @threads for bi in block_tasks
        if bi==0; continue;end #empty job
        ret = [0,0,0]; ret2 = [0,0,0]
        idim = tdims[bi]
        l_Np = length(pbits[bi])
        l_Nn = length(nbits[bi])
        offset = idim -l_Nn

        #@timeit to "pn1" begin
        @inbounds for (bfidx,bf) in enumerate(bfs[bi])
            bisearch!(delMs,Mps[bf]-Mps[bi],ret) #!! bf = bfs[bi][j]
            Vs = Vpn[ret[1]]
            fdim = tdims[bf]
            l_Nn_f = length(nbits[bf])
            p_NiNf = p_NiNfs[bi][bfidx]
            n_NiNf = n_NiNfs[bi][bfidx]
            off_f = fdim-l_Nn_f
            @inbounds for (nth,V) in enumerate(Vs)
                Npifs = p_NiNf[nth]; Nnifs = n_NiNf[nth]
                @inbounds for Npif in Npifs
                    tMi = offset+ Npif.i .* l_Nn
                    tMf = off_f + Npif.f .* l_Nn_f
                    phase_p = Npif.phase
                    @inbounds for Nnif in Nnifs
                        Mi = tMi + Nnif.i; Mf = tMf + Nnif.f
                        phase_n = Nnif.phase
                        twf[Mi] += ifelse(phase_p!=phase_n,-V,V) .* wf[Mf]
                    end
                end
            end
        end
        #end
        #@timeit to "pn2" begin
        @inbounds for j = 1:length(bis[bi])-1    #### tbf = bi !!!
            tbi = bis[bi][j]
            bisearch!(delMs,Mps[bi]-Mps[tbi],ret)
            bisearch_ord!(bfs[tbi],bi,ret2)
            fdim=tdims[tbi]
            l_Nn_i=length(nbits[tbi])
            p_NiNf = p_NiNfs[tbi][ret2[1]]
            n_NiNf = n_NiNfs[tbi][ret2[1]]
            off_f = fdim - l_Nn_i
            @inbounds for (nth,V) in enumerate(Vpn[ret[1]])
                Npifs = p_NiNf[nth]
                Nnifs = n_NiNf[nth]
                @inbounds for Npif in Npifs
                    tMi = off_f  + Npif.i *l_Nn_i #idim <-> fdim
                    tMf = offset + Npif.f*l_Nn
                    phase_p = Npif.phase
                    @inbounds for Nnif in Nnifs
                        Mi = tMi + Nnif.i; Mf = tMf + Nnif.f
                        phase_n = Nnif.phase
                        twf[Mf] += ifelse(phase_p!=phase_n,-V,V) .* wf[Mi]
                    end
                end
            end
        end
        #end
        #@timeit to "pp/nn" begin
        ### pp/nn interaction
        @inbounds for (Npi,tinfo) in enumerate(ppinfo[bi])
            tMi = offset + Npi*l_Nn
            @inbounds for (jj,tmp) in enumerate(tinfo)
                tMf = offset + l_Nn*tmp.f
                fac = tmp.coef
                @inbounds @simd for nidx = 1:l_Nn
                    Mi = tMi + nidx; Mf = tMf + nidx
                    twf[Mf] += fac .* wf[Mi]
                    twf[Mi] += fac .* wf[Mf]
                end
            end
        end
        @inbounds for (Nni,tinfo) in enumerate(nninfo[bi])
            tMi = offset + Nni
            @inbounds for (jj,tmp) in enumerate(tinfo)
                tMf = offset + tmp.f
                fac = tmp.coef
                @inbounds @simd for pidx = 1:l_Np
                    Mi = tMi + pidx*l_Nn; Mf = tMf + pidx*l_Nn
                    twf[Mf] += fac .* wf[Mi]
                    twf[Mi] += fac .* wf[Mf]
                end
            end
        end
        #end
        #@timeit to "1b" begin
        ### one-body operator
        @inbounds for pidx = 1:l_Np
            tMi = offset + pidx*l_Nn
            @inbounds for nidx =1:l_Nn
                Mi = tMi + nidx
                twf[Mi] += (dot(SPEs[1],jocc_p[bi][pidx])+
                            dot(SPEs[2],jocc_n[bi][nidx])) .* wf[Mi]
            end
        end
        #end
    end
    return nothing
end

function diagonalize_T!(k::Int64,num_ev::Int64,
                        Tmat,
                        en::Array{Array{Float64,1}},
                        num_history::Int64,
                        TF::Array{Bool,1},
                        tol::Float64) where{FA<:Array{Float64,1},
                                            FA2<:Array{Float64,2}}
    for ith = num_history:-1:2
        en[ith] .= en[ith-1]
    end
    n = minimum([num_ev,k])
    @views en[1][1:n] .= eigvals(@views Tmat[1:k,1:k] )[1:n]
    if all( en[2]-en[1] .< tol) ; TF[1]=true; end
    return nothing
end

function TRBL(q,vks,uks,Tmat,Beta_H,pbits,nbits,jocc_p,jocc_n,SPEs,
               ppinfo,nninfo,bis,bfs,block_tasks,
               p_NiNfs,n_NiNfs,Mps,delMs,Vpn,tdims,
               eval_jj,oPP,oNN,oPNu,oPNd,Jidxs,
               num_ev,num_history,lm,ls_sub,en,tol,to,
               doubleLanczos=false)
    ls = q*ls_sub
    mdim = tdims[end]
    TF=[false]
    elit = 1
    inow = 1; itmin = 1; itmax = div(lm,q)-1
    rescount = 0
    Vt = zeros(Float64,q,mdim)
    R  = zeros(Float64,q,q)
    Beta_J = zeros(Float64,q,q)

    lnJ = 20 # maximum([q*4,10])
    Jvs = [zeros(Float64,q,mdim) for i=1:lnJ]
    Jvret = [zeros(Float64,q,mdim) ]
    Jmat = zeros(Float64,lnJ*q,lnJ*q)
    JTF = [false];Jtol=1.e-6
    U = zeros(Float64,mdim,lnJ*q)
    Mat = zeros(Float64,mdim,q)
    if doubleLanczos
        @timeit to "JJ lanczos" begin
            Jvs[1] .= vks[1]
            bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,
                          pbits,nbits,tdims,eval_jj,
                          Jidxs,oPP,oNN,oPNu,oPNd,to,Beta_H,U,Mat)
            vks[1] .= Jvs[1];Beta_J .= 0.0;JTF[1]=false
            for i=1:lnJ;Jvs[i] .= 0.0;end
            V = vks[1]
        end
    end
    while TF[1]==false
        for it = itmin:itmax
            inow = it
            V  = vks[it]
            HV = vks[it+1]
            @timeit to "bl_operateH" begin
                bl_operate_H!(q,V,HV,pbits,nbits,
                              jocc_p,jocc_n,SPEs,ppinfo,nninfo,
                              tdims,bis,bfs,block_tasks,
                              p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to)
            end
            BLAS.gemm!('N','T',1.0,V,HV,0.0,R) ##mul!(R,V,HV')
            if issymmetric(R) == false
                @inbounds for i=1:q;for j=i:q
                    R[i,j] = R[j,i]
                end;end
            end
            @views Tmat[q*it-q+1:q*it,q*it-q+1:q*it] .= R
            diagonalize_T!(it*q,num_ev,Tmat,en,num_history,TF,tol)
            #print_vec("En TR(d)BL $it ",en[1])
            BLAS.gemm!('N','N',-1.0,R,V,1.0,HV)#mul!(HV,R,V,-1.0,1.0)
            s_vks = @views vks[1:it-1]
            @timeit to "ReORTH" bl_ReORTH(q,it,mdim,s_vks,HV,Vt,R)
            bl_QR!(HV',Beta_H,mdim,q)#myQR!(HV',Beta_H,mdim,q)
            if doubleLanczos
                tnorm = norm(Diagonal(Beta_H),Inf)
                if tnorm < Jtol
                    println("Hbn norm $tnorm");TF[1]=true;elit=it;break
                end
                @timeit to "JJ lanczos" begin
                    for i=2:lnJ; Jvs[i] .= 0.0; end
                    JTF[1] = false; Jmat .= 0.0;Jvs[1] .= HV
                    Vt .= 0.0; R .= 0.0
                    bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,
                                  tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,to,Beta_H,U,Mat)
                    HV .= Jvs[1]; Beta_J .=0.0
                end
            end
            add_bl_T!(q,it,Tmat,Beta_H)
            if TF[1];elit = it;break;end
        end
        if TF[1] == false
            @timeit to "bl_Restart" begin
                bl_ThickRestart(q,vks,uks,Beta_H,Tmat,inow,ls_sub,mdim,Vt)
            end
            itmin = ls_sub + 1
        end
    end
    return elit
end

function bl_ReORTH(q,it,mdim,vks,HV,Vt,R)
    #LinearAlgebra.jl:  mul!(C, A, B, α, β) -> C: ABalpha + Cbeta
    #BLAS.gemm!(tA, tB, alpha, A, B,beta,C)   C := alpha*A*B + beta*C
    # :tA,tB 'N': normal 'T': transpose
    @inbounds for l = it-1:-1:1
        Vt .= vks[l]
        BLAS.gemm!('N','T', 1.0,HV,Vt,0.0,R) #mul!(R,HV,Vt')
        BLAS.gemm!('N','N',-1.0,R,Vt,1.0,HV) #mul!(HV,R,Vt,-1.0,1.0)
    end
    R .= 0.0
    return nothing
end

# wf:V, twf:HV
function bl_operate_H!(q,wf,twf,pbits,nbits,jocc_p,jocc_n,
                       SPEs,ppinfo,nninfo,tdims,bis,bfs,block_tasks,
                       p_NiNfs,n_NiNfs,Vpn, Mps,delMs,to=nothing)
    #lblock = length(pbits)
    #@inbounds @qthreads for bi=1:lblock
    @inbounds @threads for bi in block_tasks
        if bi==0; continue;end #empty job
        ret = [0,0,0]; ret2 = [0,0,0]
        idim = tdims[bi]
        l_Np = length(pbits[bi])
        l_Nn = length(nbits[bi])
        offset = idim -l_Nn
        #@timeit to "pn1" begin
        @inbounds for (bfidx,bf) in enumerate(bfs[bi])
            bisearch!(delMs,Mps[bf]-Mps[bi],ret) #!! bf = bfs[bi][j]
            Vs = Vpn[ret[1]]
            fdim = tdims[bf]
            l_Nn_f = length(nbits[bf])
            p_NiNf = p_NiNfs[bi][bfidx]
            n_NiNf = n_NiNfs[bi][bfidx]
            off_f = fdim-l_Nn_f
            @inbounds for (nth,V) in enumerate(Vs)
                Npifs = p_NiNf[nth]; Nnifs = n_NiNf[nth]
                @inbounds for Npif in Npifs
                    tMi = offset+ Npif.i .* l_Nn
                    tMf = off_f + Npif.f .* l_Nn_f
                    phase_p = Npif.phase
                    @inbounds for Nnif in Nnifs
                        Mi = tMi + Nnif.i; Mf = tMf + Nnif.f
                        phase_n = Nnif.phase
                        coeff = ifelse(phase_p!=phase_n,-V,V)
                        w_f = @views twf[:,Mi]
                        w_i = @views  wf[:,Mf]
                        @inbounds @simd for b=1:q
                            w_f[b] += coeff .* w_i[b]
                        end
                    end
                end
            end
        end
        #end
        #@timeit to "pn2" begin
        @inbounds for j = 1:length(bis[bi])-1    #### tbf = bi !!!
            tbi = bis[bi][j]
            bisearch!(delMs,Mps[bi]-Mps[tbi],ret)
            bisearch_ord!(bfs[tbi],bi,ret2)
            fdim=tdims[tbi]
            l_Nn_i=length(nbits[tbi])
            p_NiNf = p_NiNfs[tbi][ret2[1]]
            n_NiNf = n_NiNfs[tbi][ret2[1]]
            off_f = fdim - l_Nn_i
            @inbounds for (nth,V) in enumerate(Vpn[ret[1]])
                Npifs = p_NiNf[nth]
                Nnifs = n_NiNf[nth]
                @inbounds for Npif in Npifs
                    tMi = off_f  + Npif.i *l_Nn_i #idim <-> fdim
                    tMf = offset + Npif.f*l_Nn
                    phase_p = Npif.phase
                    @inbounds for Nnif in Nnifs
                        Mi = tMi + Nnif.i; Mf = tMf + Nnif.f
                        phase_n = Nnif.phase
                        coeff = ifelse(phase_p!=phase_n,-V,V)
                        w_f = @views twf[:,Mf]
                        w_i = @views  wf[:,Mi]
                        @inbounds @simd for b=1:q
                            w_f[b] += coeff .* w_i[b]
                        end
                    end
                end
            end
        end
        #end
        #@timeit to "pp/nn" begin
        ### pp/nn interaction
        @inbounds for (Npi,tinfo) in enumerate(ppinfo[bi])
            tMi = offset + Npi*l_Nn
            @inbounds for (jj,tmp) in enumerate(tinfo)
                tMf = offset + l_Nn*tmp.f
                fac = tmp.coef
                @inbounds @simd for nidx = 1:l_Nn
                    Mi = tMi + nidx; Mf = tMf + nidx
                    w_f1 = @views twf[:,Mf]
                    w_i1 = @views  wf[:,Mi]
                    w_f2 = @views twf[:,Mi]
                    w_i2 = @views  wf[:,Mf]
                    @inbounds @simd for b=1:q
                        w_f1[b] += fac .* w_i1[b]
                        w_f2[b] += fac .* w_i2[b]
                    end
                end
            end
        end
        @inbounds for (Nni,tinfo) in enumerate(nninfo[bi])
            tMi = offset + Nni
            @inbounds for (jj,tmp) in enumerate(tinfo)
                tMf = offset + tmp.f
                fac = tmp.coef
                @inbounds @simd for pidx = 1:l_Np
                    Mi = tMi + pidx*l_Nn; Mf = tMf + pidx*l_Nn
                    w_f1 = @views twf[:,Mf]
                    w_i1 = @views  wf[:,Mi]
                    w_f2 = @views twf[:,Mi]
                    w_i2 = @views  wf[:,Mf]
                    @inbounds @simd for b=1:q
                        w_f1[b] += fac .* w_i1[b]
                        w_f2[b] += fac .* w_i2[b]
                    end
                end
            end
        end
        #end
        #@timeit to "1b" begin
        ### one-body operator
        @inbounds for pidx = 1:l_Np
            tMi = offset + pidx*l_Nn
            @inbounds for nidx =1:l_Nn
                Mi = tMi + nidx
                coeff = (dot(SPEs[1],jocc_p[bi][pidx])+
                         dot(SPEs[2],jocc_n[bi][nidx]))
                w_f = @views twf[:,Mi]
                w_i = @views  wf[:,Mi]
                @inbounds @simd for b=1:q
                    w_f[b] += coeff .* w_i[b]
                end

            end
        end
        #end
    end
    return nothing
end

function Jcompress(q,vks,Tmat,inow,ls_sub,mdim,R,bnout,U,Mat,Beta_J;use_LAPACK=false)
    lm = q*inow; ls = q*ls_sub
    vals = [0.0]
    vecs = [0.0;0.0]
    if use_LAPACK # used for only for debug (you need "l_diag.F90" by S.Yoshida)
        M = Tmat[1:lm,1:lm]
        ccall((:diagonalize_double_,"l_diag.so"),Nothing,
              (Ref{Int64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Int64}),
              lm,M,vals,vecs,ls)
    else
        M = @views Tmat[1:lm,1:lm]
        vals,vecs = eigen(M) # eigen(Symmetric(M,:L))
    end
    # x = maximum([ 1.0-sum(vecs[1:lm,i].^2) for i=1:ls])
    # if x > 1.e-6
    #     println("WARNING: JJ block norm diff. in jj_refine ", x)
    # end

    tv = @views vecs[1:ls,1:ls]
    BLAS.gemm!('T','N',1.0,tv,bnout,0.0,R)# mul!(R,tv,bnout)
    bnout .= R

    Tmat .= 0.0; U .= 0.0
    for i=1:inow
        bv = vks[i]
        for b=1:q
            j = q*(i-1)+b
            u = @views U[:,j]
            u .= @views bv[b,:]
        end
    end
    tU = @views U[:,1:lm]
    tv = @views vecs[1:lm,1:ls]
    BLAS.gemm!('T','T',1.0,tv,tU,0.0,vks[1])
    Beta_J .= @views vecs[1:ls,1:ls]

    #println("V1 compressed ",vks[1][1,:])
    #println("V2 compressed ",vks[1][2,:])
    return nothing
end

function bl_ThickRestart(q,vks,uks,R,Tmat,inow,ls_sub,mdim,Vt)
    lm = q*inow
    ls = q*ls_sub
    r = zeros(Float64,q,ls)
    vals,vecs = eigen(@views Tmat[1:lm,1:lm])
    Tmat .= 0.0
    for k = 1:ls; Tmat[k,k] = vals[k]; end
    tv = @views vecs[lm-q+1:lm,1:ls]
    mul!(r,R,tv) #BLAS.gemm!('N','N',1.0,R,tv,0.0,r)
    Tmat[ls+1:ls+q,1:ls] .= r
    Tmat[1:ls,ls+1:ls+q] .= r'
    for bi = 1:ls_sub
        uk = uks[bi] #.= 0.0
        uk .= 0.0
        for b=1:q
            k = q*(bi-1) + b
            for j=1:inow
                vk = vks[j]
                for bj = 1:q
                    v = @views vk[bj,:]
                    idx = q*(j-1) +bj
                    fac = vecs[idx,k]
                    @inbounds @simd for m=1:mdim
                        uk[b,m] += fac .* v[m]
                    end
                end
            end
        end
    end
    for i = 1:ls_sub
        vks[i] .= uks[i]
    end
    vks[ls_sub+1] .= vks[inow+1]
    for i= ls_sub+2:length(vks)
        vks[i] .= 0.0
    end
    for i=1:ls_sub+1
        V = vks[i]
        if i>1
            for j = i-1:-1:1
                Vt .= vks[j]
                mul!(R,V,Vt')
                mul!(V,R,Vt,-1.0,1.0)
            end
            bl_QR!(V',R,mdim,q)
        end
    end
    return nothing
end

function bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,
                       pbits,nbits,tdims,eval_jj,
                       Jidxs,oPP,oNN,oPNu,oPNd,to,bnout,
                       U,Mat;verbose=false)
    mdim = tdims[end]
    lnJ = length(Jvs)
    itmin = 1;itmax = lnJ-1; inow = 0
    rescount = 0
    while JTF[1] == false
        for it = itmin:itmax
            inow = it
            V = Jvs[it]; JV = Jvs[it+1]
            bl_operate_J!(q,V,JV,pbits,nbits,tdims,
                          Jidxs,oPP,oNN,oPNu,oPNd)
            axpy!(eval_jj,V,JV)##JV .+= eval_jj .* V
            BLAS.gemm!('N','T',1.0,V,JV,0.0,R)
            @inbounds for j=1:q;for i=1:j
                R[j,i] = R[i,j]
            end;end
            @views Jmat[q*it-q+1:q*it,q*it-q+1:q*it] .= R
            tJmat = @views Jmat[1:it*q,1:it*q]

            jeval = eigvals(tJmat)[1:q]
            if verbose; print_vec("$it jeval=> ",jeval);end
            if jeval[1] - eval_jj < -Jtol
                if verbose;println("neg JJ @$it jeval");end
                JTF[1] = true;break
            end
            if all( abs.(jeval.-eval_jj) .< Jtol)
                if verbose;println("J converged @$it");end
                JTF[1] = true;break
            end
            BLAS.gemm!('N','N',-1.0,R,V,1.0,JV) #mul!(JV,R,V,-1.0,1.0)
            s_vs = @views Jvs[1:it-1]
            bl_ReORTH(q,it,mdim,s_vs,JV,Vt,R)
            bl_QR!(JV',Beta_J,mdim,q)

            tnorm = norm(Diagonal(Beta_J),Inf)
            add_bl_T!(q,it,Jmat,Beta_J)
            if tnorm < Jtol
                println("Jbn norm $tnorm")
                JTF[1]=true;break
            end
        end
        if rescount == 20;println("JJ not converged");return nothing;end
        if JTF[1]==false
            tJmat = @views Jmat[1:inow*q,1:inow*q]
            bl_ThickRestart(q,Jvs,Jvret,Beta_J,Jmat,
                            inow,1,mdim,Vt)
            itmin = 2;rescount += 1
        end
    end
    if inow > itmin
        #@timeit to "Restart"
        Jcompress(q,Jvs,Jmat,inow,1,mdim,R,bnout,U,Mat,Beta_J)
    end
    return nothing
end

function bl_operate_J!(q,Rvec,Jv,
                       pbits,nbits,tdims,
                       Jidxs,
                       oPP,oNN,oPNu,oPNd,
                       beta_J=1.0)
    lblock=length(pbits)
    @inbounds @threads for bi in Jidxs
        if bi==0;continue;end
        idim = tdims[bi]
        lNn = length(nbits[bi])
        lNp = length(pbits[bi])
        opPP = oPP[bi]
        opNN = oNN[bi]
        offset = idim-lNn
        @inbounds for tmp in opPP
            Npi =tmp.Mi; Npf=tmp.Mf; fac=tmp.fac .* beta_J
            tMi = offset + Npi*lNn
            tMf = offset + Npf*lNn
            @inbounds for nidx = 1:lNn
                Mi = tMi+nidx; Mf = tMf+nidx
                w_f1 = @views Jv[:,Mf]
                w_i1 = @views Rvec[:,Mi]
                w_f2 = @views Jv[:,Mi]
                w_i2 = @views Rvec[:,Mf]
                @inbounds @simd for b=1:q
                    w_f1[b] += fac .* w_i1[b]
                    w_f2[b] += fac .* w_i2[b]
                end
            end
        end
        @inbounds for tmp in opNN #nn
            Nni =tmp.Mi; Nnf=tmp.Mf; fac=tmp.fac .* beta_J
            tMi = offset + Nni
            tMf = offset + Nnf
            @inbounds for pidx = 1:lNp
                Mi = tMi+pidx*lNn; Mf = tMf+pidx*lNn
                w_f1 = @views Jv[:,Mf]
                w_i1 = @views Rvec[:,Mi]
                w_f2 = @views Jv[:,Mi]
                w_i2 = @views Rvec[:,Mf]
                @inbounds @simd for b=1:q
                    w_f1[b] += fac .* w_i1[b]
                    w_f2[b] += fac .* w_i2[b]
                end
            end
        end
    end
    @inbounds @threads for bi = 2:lblock
        operator = oPNd[bi]
        bf = bi-1
        l_Nn_i = length(nbits[bi])
        l_Nn_f = length(nbits[bf])
        off_i = tdims[bi] - l_Nn_i
        off_f = tdims[bf] - l_Nn_f
        @inbounds for top in operator
            pj = top.pjump
            nj = top.njump
            fac =top.fac .* beta_J
            @inbounds for tmp_p in pj
                phase_p=tmp_p.phase
                tMi = off_i + tmp_p.i * l_Nn_i
                tMf = off_f + tmp_p.f * l_Nn_f
                @inbounds for tmp_n in nj
                    phase_n=tmp_n.phase
                    Mi = tMi + tmp_n.i
                    Mf = tMf + tmp_n.f
                    coeff = ifelse(phase_p!=phase_n,-fac,fac)
                    w_f = @views Jv[:,Mf]; w_i = @views Rvec[:,Mi]
                    @inbounds @simd for b=1:q
                        w_f[b] += coeff .* w_i[b]
                    end
                end
            end
        end
    end
    @inbounds @threads for bi = 1:lblock-1
        operator = oPNu[bi]
        bf = bi+1
        l_Nn_i = length(nbits[bi])
        l_Nn_f = length(nbits[bf])
        off_i = tdims[bi] - l_Nn_i
        off_f = tdims[bf] - l_Nn_f
        @inbounds for top in operator
            pj  = top.pjump
            nj  = top.njump
            fac = top.fac .* beta_J
            @inbounds for tmp_p in pj
                phase_p=tmp_p.phase
                tMi = off_i + tmp_p.i * l_Nn_i
                tMf = off_f + tmp_p.f * l_Nn_f
                @inbounds for tmp_n in nj
                    Nni = tmp_n.i; Nnf = tmp_n.f; phase_n=tmp_n.phase
                    Mi = tMi + tmp_n.i
                    Mf = tMf + tmp_n.f
                    coeff = ifelse(phase_p!=phase_n,-fac,fac)
                    w_f = @views Jv[:,Mf]; w_i = @views Rvec[:,Mi]
                    @inbounds @simd for b=1:q
                        w_f[b] += coeff .* w_i[b]
                    end
                end
            end
        end
    end
    return nothing
end


text:
# Common Functions

### Dimensionality Checks

function check_dims(X::AbstractMatrix; dims::Integer)
    dims ∈ (1, 2) || throw(ArgumentError("dims should be 1 or 2 (got $dims)"))

    n = size(X, dims)
    p = size(X, dims == 1 ? 2 : 1)

    return (n, p)
end

### Data/Parameter Validation

function validate_priors(π::AbstractVector{T}) where T
    total = zero(T)

    for (k, πₖ) in enumerate(π)
        if !(zero(T) < πₖ < one(T))
            throw(DomainError(πₖ, "class prior probability at class index $(k) must be " *
                                  "in the open interval (0,1)"))
        end
        total += πₖ
    end

    if !isapprox(total, one(T))
        throw(ArgumentError("class priors vector π must sum to 1 (got $(total))"))
    end

    return length(π)
end

### Class Calculations

"""
    class_counts!(nₘ, y)

Overwrite vector `nₘ` with counts of each index in `y`
"""
function class_counts!(nₘ::Vector{T}, y::Vector{<:Integer}) where {T<:Integer}
    m = length(nₘ)

    nₘ .= zero(T)

    for i = 1:length(y)
        yᵢ = y[i]
        1 ≤ yᵢ ≤ m || throw(BoundsError(nₘ, yᵢ))
        nₘ[yᵢ] += 1
    end

    return nₘ
end

"""
    _class_statistics!(M, nₘ, X, y)

Backend for `class_statistics!` - assumes `dims=2`.
"""
function _class_statistics!(M::AbstractMatrix, nₘ::Vector{<:Integer}, X::AbstractMatrix,
                            y::Vector{<:Integer})
    p, n = size(X)
    p₂, m = size(M)
    n₂ = length(y)
    m₂ = length(nₘ)

    p == p₂ || throw(DimensionMismatch("predictor count mismatch between M and X"))
    n == n₂ || throw(DimensionMismatch("observation count mismatch between y and X"))
    m == m₂ || throw(DimensionMismatch("class count mismatch between nₘ and M"))

    T = eltype(nₘ)

    M .= zero(eltype(M))
    nₘ .= zero(T)  # track counts to ensure an observation for each class

    for i = 1:n
        yᵢ = y[i]
        1 ≤ yᵢ ≤ m || throw(BoundsError(M, (1, yᵢ)))
        nₘ[yᵢ] += one(T)
        @inbounds for j = 1:p
            M[j, yᵢ] += X[j, i]
        end
    end

    all(nₘ .≥ 1) || error("must have at least one observation per class")
    broadcast!(/, M, M, transpose(nₘ))

    return (M, nₘ)
end

"""
    class_statistics!(M, nₘ, X, y; dims=1)

Overwrites matrix `M` with class centroids from `X` based on class indexes from `y`. Use
`dims=1` for row-based observations and `dims=2` for column-based observations.
"""
function class_statistics!(M::AbstractMatrix, nₘ::Vector{<:Integer}, X::AbstractMatrix,
                           y::Vector{<:Integer}; dims::Integer=1)
    n, p = check_dims(X, dims=dims)
    m, p₂ = check_dims(M, dims=dims)
    n₂ = length(y)
    m₂ = length(nₘ)

    altdims = dims == 1 ? 2 : 1

    p == p₂ || throw(DimensionMismatch("predictor count along dimension $(altdims) of X " *
                                       "must match dimension $(altdims) of M (got $(p) " *
                                       "and $(p₂))"))
    n == n₂ || throw(DimensionMismatch("observation count along length of y must match " *
                                       "dimension $(dims) of X (got $(n) and $(n₂))"))
    m == m₂ || throw(DimensionMismatch("class count along length of nₘ must match " *
                                       "dimension $(dims) of M (got $(m₂) and $(m))"))

    if dims == 1
        _class_statistics!(transpose(M), nₘ, transpose(X), y)
        return (M, nₘ)
    else dims == 2
        return _class_statistics!(M, nₘ, X, y)
    end
end

### Regularization

"""
    regularize!(Σ₁, Σ₂, λ)
"""
function regularize!(Σ₁::AbstractMatrix{T}, Σ₂::AbstractMatrix{T}, λ::T) where {T}
    (m = size(Σ₁, 1)) == size(Σ₂, 1) || throw(DimensionMismatch(""))
    (n = size(Σ₁, 2)) == size(Σ₂, 2) || throw(DimensionMismatch(""))

    m == n || throw(DimensionMismatch("matrices Σ₁ and Σ₂ must be square"))

    0 ≤ λ ≤ 1 || throw(DomainError(λ, "λ must be in the interval [0,1]"))

    for ij in eachindex(Σ₁)
        Σ₁[ij] = (1-λ)*Σ₁[ij] + λ*Σ₂[ij]
    end

    return Σ₁
end

@inline regularize(x, y, γ) = (1-γ)*x + γ*y

"""
    regularize!(Σ, γ)

Shrink `Σ` matrix towards the average eigenvalue multiplied by the identity matrix.
"""
function regularize!(Σ::AbstractMatrix{T}, γ::T) where {T}
    (p = size(Σ, 1)) == size(Σ, 2) || throw(DimensionMismatch("Σ must be square"))

    0 ≤ γ ≤ 1 || throw(DomainError(γ, "γ must be in the interval [0,1]"))

    a =  γ*tr(Σ)/p  # Average eigenvalue scaled by γ

    broadcast!(*, Σ, Σ, 1 - γ)
    for i = 1:p
        Σ[i, i] += a
    end

    return Σ
end


text:
import Tokenize.Lexers: peekchar, readchar, iswhitespace, emit, emit_error,  accept_batch, eof

const EmptyWS = Tokens.EMPTY_WS
const SemiColonWS = Tokens.SEMICOLON_WS
const NewLineWS = Tokens.NEWLINE_WS
const WS = Tokens.WS
const EmptyWSToken = RawToken(EmptyWS, (0, 0), (0, 0), -1, -1)

mutable struct Closer
    newline::Bool
    semicolon::Bool
    tuple::Bool
    comma::Bool
    paren::Bool
    brace::Bool
    inmacro::Bool
    insquare::Bool
    inref::Bool
    inwhere::Bool
    square::Bool
    block::Bool
    ifop::Bool
    range::Bool
    ws::Bool
    wsop::Bool
    unary::Bool
    for_generator::Bool
    precedence::Int
end
Closer() = Closer(true, true, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, -1)

mutable struct ParseState
    l::Lexer{Base.GenericIOBuffer{Array{UInt8,1}},RawToken}
    done::Bool # Remove this
    lt::RawToken
    t::RawToken
    nt::RawToken
    nnt::RawToken
    lws::RawToken
    ws::RawToken
    nws::RawToken
    nnws::RawToken
    closer::Closer
    errored::Bool
end
function ParseState(str::Union{IOBuffer,String})
    ps = ParseState(tokenize(str, RawToken), false, RawToken(), RawToken(), RawToken(), RawToken(), RawToken(), RawToken(), RawToken(), RawToken(), Closer(), false)
    return next(next(ps))
end

function ParseState(str::Union{IOBuffer,String}, loc::Int)
    ps = ParseState(str)
    prevpos = position(ps)
    while ps.nt.startbyte < loc
        next(ps)
        prevpos = loop_check(ps, prevpos)
    end
    return ps
end

function Base.show(io::IO, ps::ParseState)
    println(io, "ParseState at $(position(ps.l.io))")
    println(io, "last    : ", kindof(ps.lt), " ($(ps.lt))", "    ($(wstype(ps.lws)))")
    println(io, "current : ", kindof(ps.t), " ($(ps.t))", "    ($(wstype(ps.ws)))")
    println(io, "next    : ", kindof(ps.nt), " ($(ps.nt))", "    ($(wstype(ps.nws)))")
end
peekchar(ps::ParseState) = peekchar(ps.l)
if !applicable(Base.peek, Tuple{IOBuffer, Char})
function _peek(s::IO, ::Type{T}) where T
    mark(s)
    try read(s, T)::T
    finally
        reset(s)
    end
end
peekchar(io) = _peek(io, Char)
else
peekchar(io) = peek(io, Char)
end

wstype(t::AbstractToken) = kindof(t) == EmptyWS ? "empty" :
                    kindof(t) == NewLineWS ? "ws w/ newline" :
                    kindof(t) == SemiColonWS ? "ws w/ semicolon" : "ws"

function next(ps::ParseState)
    #  shift old tokens
    ps.lt = ps.t
    ps.t = ps.nt
    ps.nt = ps.nnt
    ps.lws = ps.ws
    ps.ws = ps.nws
    ps.nws = ps.nnws

    ps.nnt = Tokenize.Lexers.next_token(ps.l)

    # combines whitespace, comments and semicolons
    if iswhitespace(peekchar(ps.l)) || peekchar(ps.l) == '#' || peekchar(ps.l) == ';'
        ps.nnws = lex_ws_comment(ps.l, readchar(ps.l))
    else
        ps.nnws = EmptyWSToken
    end

    return ps
end

function Base.seek(ps::ParseState, offset)
    seek(ps.l, offset)
    next(next(ps))
end

Base.position(ps::ParseState) = ps.nt.startbyte

current_line(ps::ParseState) = ps.nt.startpos[1]
"""
    lex_ws_comment(l::Lexer, c)

Having hit an initial whitespace/comment/semicolon continues collecting similar
`Chars` until they end. Returns a WS token with an indication of newlines/ semicolons. Indicating a semicolons takes precedence over line breaks as the former is equivalent to the former in most cases.
"""
function read_ws_comment(l, c::Char)
    newline = c == '\n'
    semicolon = c == ';'
    if c == '#'
        newline = read_comment(l)
    else
        newline, semicolon = read_ws(l, newline, semicolon)
    end
    while iswhitespace(peekchar(l)) || peekchar(l) == '#' || peekchar(l) == ';'
        c = readchar(l)
        if c == '#'
            read_comment(l)
            newline = newline || peekchar(l) == '\n'
            semicolon = semicolon || peekchar(l) == ';'
        elseif c == ';'
            semicolon = true
        else
            newline |= c == '\n'
            newline, semicolon = read_ws(l, newline, semicolon)
        end
    end
    return newline, semicolon
end

function lex_ws_comment(l::Lexer, c::Char)
    newline, semicolon = read_ws_comment(l, c)
    return emit(l, semicolon ? SemiColonWS :
                   newline ? NewLineWS : WS)
end

function read_ws(l, newline, semicolon)
    while iswhitespace(peekchar(l))
        c = readchar(l)
        c == '\n' && (newline = true)
        c == ';' && (semicolon = true)
    end
    return newline, semicolon
end

function read_comment(l)
    if peekchar(l) != '='
        while true
            pc = peekchar(l)
            if pc == '\n' || eof(pc)
                return true
            end
            readchar(l)
        end
    else
        c = readchar(l) # consume the '='
        n_start, n_end = 1, 0
        while true
            if eof(c)
                return false
            end
            nc = readchar(l)
            if c == '#' && nc == '='
                n_start += 1
            elseif c == '=' && nc == '#'
                n_end += 1
            end
            if n_start == n_end
                return false
            end
            c = nc
        end
    end
end

# Functions relating to tokens
isemptyws(t::AbstractToken) = kindof(t) == EmptyWS
isnewlinews(t::AbstractToken) = kindof(t) === NewLineWS
isendoflinews(t::AbstractToken) = kindof(t) == SemiColonWS || kindof(t) == NewLineWS
@inline val(token::AbstractToken, ps::ParseState) = String(ps.l.io.data[token.startbyte + 1:token.endbyte + 1])
both_symbol_and_op(t::AbstractToken) = kindof(t) === Tokens.WHERE || kindof(t) === Tokens.IN || kindof(t) === Tokens.ISA
isprefixableliteral(t::AbstractToken) = (kindof(t) === Tokens.STRING || kindof(t) === Tokens.TRIPLE_STRING || kindof(t) === Tokens.CMD || kindof(t) === Tokens.TRIPLE_CMD)
isassignmentop(t::AbstractToken) = Tokens.begin_assignments < kindof(t) < Tokens.end_assignments

isidentifier(t::AbstractToken) = kindof(t) === Tokens.IDENTIFIER
isliteral(t::AbstractToken) = Tokens.begin_literal < kindof(t) < Tokens.end_literal
isbool(t::AbstractToken) =  Tokens.TRUE ≤ kindof(t) ≤ Tokens.FALSE
iscomma(t::AbstractToken) =  kindof(t) === Tokens.COMMA
iscolon(t::AbstractToken) =  kindof(t) === Tokens.COLON
iskeyword(t::AbstractToken) = Tokens.iskeyword(kindof(t))
isinstance(t::AbstractToken) = isidentifier(t) || isliteral(t) || isbool(t) || iskeyword(t)
ispunctuation(t::AbstractToken) = iscomma(t) || kindof(t) === Tokens.END || Tokens.LSQUARE ≤ kindof(t) ≤ Tokens.RPAREN || kindof(t) === Tokens.AT_SIGN
isnumberliteral(t::AbstractToken) = isliteral(t) && (Tokens.INTEGER ≤ kindof(t) ≤ Tokens.FLOAT)


text:
using MLJLinearModels
using BenchmarkTools, Random, LinearAlgebra
DO_COMPARISONS = false; include("../testutils.jl")

n, p = 50_000, 500
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.5)

# ======================== #
# ELEMENTARIES
#
# > apply_X!        ✅ Sept 13, 2019
# > apply_Xt!       ✅ Sept 13, 2019
# > get_residuals!  ✅ Sept 13, 2019
# ======================== #

# ------------
# No intercept

R.allocate(n, p)

Xθ = similar(y)
Xtθ = similar(θ)
r = R.SCRATCH_N[]

# Sept 13, 2019 :: 7.83 ms (0 allocations: 0 bytes)
@btime R.apply_X!($Xθ, $X, $θ);

# Sept 13, 2019 :: 7.94 ms (0 allocations: 0 bytes) [pretty much only apply_X!]
@btime R.get_residuals!($r, $X, $θ, $y);

# Sept 13, 2019 :: 7.9 ms (0 allocations: 0 bytes)
@btime R.apply_Xt!($Xtθ, $X, $y);

# --------------
# With intercept

R.allocate(n, p+1)

Xθ = similar(y)
Xtθ = similar(θ1)
r = R.SCRATCH_N[]

# Sept 13, 2019 :: 7.949 ms (1 allocation: 48 bytes) -- alloc for the view
@btime R.apply_X!($Xθ, $X, $θ1);

# Sept 13, 2019 :: 7.955 ms (1 allocation: 48 bytes) -- mostly apply_X!
@btime R.get_residuals!($r, $X, $θ1, $y);

# Sept 13, 2019 :: 7.883 ms (1 allocation: 48 bytes)
@btime R.apply_Xt!($Xtθ, $X, $y);


text:
# implementation of the GPUCompiler interfaces for generating PTX code

## target

export PTXCompilerTarget

Base.@kwdef struct PTXCompilerTarget <: AbstractCompilerTarget
    cap::VersionNumber
    ptx::VersionNumber = v"6.0" # for compatibility with older versions of CUDA.jl

    # codegen quirks
    ## can we emit debug info in the PTX assembly?
    debuginfo::Bool = false
    ## can exceptions use `exit` (which doesn't kill the GPU), or should they use `trap`?
    exitable::Bool = false

    # optional properties
    minthreads::Union{Nothing,Int,NTuple{<:Any,Int}} = nothing
    maxthreads::Union{Nothing,Int,NTuple{<:Any,Int}} = nothing
    blocks_per_sm::Union{Nothing,Int} = nothing
    maxregs::Union{Nothing,Int} = nothing
end

source_code(::PTXCompilerTarget) = "ptx"

llvm_triple(::PTXCompilerTarget) = Int===Int64 ? "nvptx64-nvidia-cuda" : "nvptx-nvidia-cuda"

function llvm_machine(target::PTXCompilerTarget)
    triple = llvm_triple(target)
    t = Target(triple=triple)

    tm = TargetMachine(t, triple, "sm_$(target.cap.major)$(target.cap.minor)",
                       "+ptx$(target.ptx.major)$(target.ptx.minor)")
    asm_verbosity!(tm, true)

    return tm
end

# the default datalayout does not match the one in the NVPTX user guide
llvm_datalayout(::PTXCompilerTarget) = Int===Int64 ?
    "e-p:64:64:64-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64"*
     "-f32:32:32-f64:64:64-v16:16:16-v32:32:32-v64:64:64-v128:128:128-n16:32:64" :
    "e-p:32:32:32-i1:8:8-i8:8:8-i16:16:16-i32:32:32-i64:64:64"*
     "-f32:32:32-f64:64:64-v16:16:16-v32:32:32-v64:64:64-v128:128:128-n16:32:64"

## job

function Base.show(io::IO, job::CompilerJob{PTXCompilerTarget})
    print(io, "PTX CompilerJob of ", job.source)
    print(io, " for sm_$(job.target.cap.major)$(job.target.cap.minor)")

    job.target.minthreads !== nothing && print(io, ", minthreads=$(job.target.minthreads)")
    job.target.maxthreads !== nothing && print(io, ", maxthreads=$(job.target.maxthreads)")
    job.target.blocks_per_sm !== nothing && print(io, ", blocks_per_sm=$(job.target.blocks_per_sm)")
    job.target.maxregs !== nothing && print(io, ", maxregs=$(job.target.maxregs)")
end

const ptx_intrinsics = ("vprintf", "__assertfail", "malloc", "free")
isintrinsic(::CompilerJob{PTXCompilerTarget}, fn::String) = in(fn, ptx_intrinsics)

# XXX: the debuginfo part should be handled by GPUCompiler as it applies to all back-ends.
runtime_slug(job::CompilerJob{PTXCompilerTarget}) =
    "ptx-sm_$(job.target.cap.major)$(job.target.cap.minor)" *
       "-debuginfo=$(Int(llvm_debug_info(job)))" *
       "-exitable=$(job.target.exitable)"

function process_module!(job::CompilerJob{PTXCompilerTarget}, mod::LLVM.Module)
    # calling convention
    if LLVM.version() >= v"8"
        for f in functions(mod)
            # JuliaGPU/GPUCompiler.jl#97
            #callconv!(f, LLVM.API.LLVMPTXDeviceCallConv)
        end
    end
end

function process_entry!(job::CompilerJob{PTXCompilerTarget}, mod::LLVM.Module, entry::LLVM.Function)
    ctx = context(mod)

    if job.source.kernel
        # work around bad byval codegen (JuliaGPU/GPUCompiler.jl#92)
        entry = lower_byval(job, mod, entry)

        # property annotations
        annotations = LLVM.Value[entry]

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

        ## expected CTA sizes
        if job.target.minthreads !== nothing
            for (dim, name) in enumerate([:x, :y, :z])
                bound = dim <= length(job.target.minthreads) ? job.target.minthreads[dim] : 1
                append!(annotations, [MDString("reqntid$name"),
                                    ConstantInt(Int32(bound), ctx)])
            end
        end
        if job.target.maxthreads !== nothing
            for (dim, name) in enumerate([:x, :y, :z])
                bound = dim <= length(job.target.maxthreads) ? job.target.maxthreads[dim] : 1
                append!(annotations, [MDString("maxntid$name"),
                                    ConstantInt(Int32(bound), ctx)])
            end
        end

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

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

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

        if LLVM.version() >= v"8"
            # calling convention
            callconv!(entry, LLVM.API.LLVMPTXKernelCallConv)
        end
    else
        # we can't look up device functions using the CUDA APIs, so alias them to a global
        gv = GlobalVariable(mod, llvmtype(entry), LLVM.name(entry) * "_slot")
        initializer!(gv, entry)
        linkage!(gv, LLVM.API.LLVMLinkOnceODRLinkage)
        set_used!(mod, gv)
    end

    return entry
end

function add_lowering_passes!(job::CompilerJob{PTXCompilerTarget}, pm::LLVM.PassManager)
    add!(pm, FunctionPass("HideUnreachable", hide_unreachable!))
    add!(pm, ModulePass("HideTrap", hide_trap!))
end

function optimize_module!(job::CompilerJob{PTXCompilerTarget}, mod::LLVM.Module)
    tm = llvm_machine(job.target)
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)

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

        # the above loop unroll pass might have unrolled regular, non-runtime nested loops.
        # that code still needs to be optimized (arguably, multiple unroll passes should be
        # scheduled by the Julia optimizer). do so here, instead of re-optimizing entirely.
        early_csemem_ssa!(pm) # TODO: gvn instead? see NVPTXTargetMachine.cpp::addEarlyCSEOrGVNPass
        dead_store_elimination!(pm)

        cfgsimplification!(pm)

        # get rid of the internalized functions; now possible unused
        global_dce!(pm)

        run!(pm, mod)
    end
end

function llvm_debug_info(job::CompilerJob{PTXCompilerTarget})
    # allow overriding the debug info from CUDA.jl
    if job.target.debuginfo
        invoke(llvm_debug_info, Tuple{CompilerJob}, job)
    else
        LLVM.API.LLVMDebugEmissionKindNoDebug
    end
end

## LLVM passes

# HACK: this pass removes `unreachable` information from LLVM
#
# `ptxas` is buggy and cannot deal with thread-divergent control flow in the presence of
# shared memory (see JuliaGPU/CUDAnative.jl#4). avoid that by rewriting control flow to fall
# through any other block. this is semantically invalid, but the code is unreachable anyhow
# (and we expect it to be preceded by eg. a noreturn function, or a trap).
#
# TODO: can LLVM do this with structured CFGs? It seems to have some support, but seemingly
#       only to prevent introducing non-structureness during optimization (ie. the front-end
#       is still responsible for generating structured control flow).
function hide_unreachable!(fun::LLVM.Function)
    job = current_job::CompilerJob
    ctx = context(fun)
    changed = false
    @timeit_debug to "hide unreachable" begin

    # remove `noreturn` attributes
    #
    # when calling a `noreturn` function, LLVM places an `unreachable` after the call.
    # this leads to an early `ret` from the function.
    attrs = function_attributes(fun)
    delete!(attrs, EnumAttribute("noreturn", 0, ctx))

    # build a map of basic block predecessors
    predecessors = Dict(bb => Set{LLVM.BasicBlock}() for bb in blocks(fun))
    @timeit_debug to "predecessors" for bb in blocks(fun)
        insts = instructions(bb)
        if !isempty(insts)
            inst = last(insts)
            if isterminator(inst)
                for bb′ in successors(inst)
                    push!(predecessors[bb′], bb)
                end
            end
        end
    end

    # scan for unreachable terminators and alternative successors
    worklist = Pair{LLVM.BasicBlock, Union{Nothing,LLVM.BasicBlock}}[]
    @timeit_debug to "find" for bb in blocks(fun)
        unreachable = terminator(bb)
        if isa(unreachable, LLVM.UnreachableInst)
            unsafe_delete!(bb, unreachable)
            changed = true

            try
                terminator(bb)
                # the basic-block is still terminated properly, nothing to do
                # (this can happen with `ret; unreachable`)
                # TODO: `unreachable; unreachable`
            catch ex
                isa(ex, UndefRefError) || rethrow(ex)
                let builder = Builder(ctx)
                    position!(builder, bb)

                    # find the strict predecessors to this block
                    preds = collect(predecessors[bb])

                    # find a fallthrough block: recursively look at predecessors
                    # and find a successor that branches to any other block
                    fallthrough = nothing
                    while !isempty(preds)
                        # find an alternative successor
                        for pred in preds, succ in successors(terminator(pred))
                            if succ != bb
                                fallthrough = succ
                                break
                            end
                        end
                        fallthrough === nothing || break

                        # recurse upwards
                        old_preds = copy(preds)
                        empty!(preds)
                        for pred in old_preds
                            append!(preds, predecessors[pred])
                        end
                    end
                    push!(worklist, bb => fallthrough)

                    dispose(builder)
                end
            end
        end
    end

    # apply the pending terminator rewrites
    @timeit_debug to "replace" if !isempty(worklist)
        let builder = Builder(ctx)
            for (bb, fallthrough) in worklist
                position!(builder, bb)
                if fallthrough !== nothing
                    br!(builder, fallthrough)
                else
                    # couldn't find any other successor. this happens with functions
                    # that only contain a single block, or when the block is dead.
                    ft = eltype(llvmtype(fun))
                    if return_type(ft) == LLVM.VoidType(ctx)
                        # even though returning can lead to invalid control flow,
                        # it mostly happens with functions that just throw,
                        # and leaving the unreachable there would make the optimizer
                        # place another after the call.
                        ret!(builder)
                    else
                        unreachable!(builder)
                    end
                end
            end
        end
    end

    end
    return changed
end

# HACK: this pass removes calls to `trap` and replaces them with inline assembly
#
# if LLVM knows we're trapping, code is marked `unreachable` (see `hide_unreachable!`).
function hide_trap!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "hide trap" begin

    # inline assembly to exit a thread, hiding control flow from LLVM
    exit_ft = LLVM.FunctionType(LLVM.VoidType(ctx))
    exit = if job.target.exitable
        InlineAsm(exit_ft, "exit;", "", true)
    else
        InlineAsm(exit_ft, "trap;", "", true)
    end

    if haskey(functions(mod), "llvm.trap")
        trap = functions(mod)["llvm.trap"]

        for use in uses(trap)
            val = user(use)
            if isa(val, LLVM.CallInst)
                let builder = Builder(ctx)
                    position!(builder, val)
                    call!(builder, exit)
                    dispose(builder)
                end
                unsafe_delete!(LLVM.parent(val), val)
                changed = true
            end
        end
    end

    end
    return changed
end


text:
export CRNN_Batch_TCS
export DNN_Batch_TCS
export RNN_Batch_TCS

"""
    DNN_Batch_TCS(p::Variable{Array{T}},
                  seqlabels::Vector,
                  inputlens;
                  background::Int=1,
                  foreground::Int=2,
                  weight=1.0) where T

a batch of concatenated input sequence is processed by neural networks into `p`

# Inputs
`p`         : 2-D Variable, probability or weighted probability\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`inputlens` : records each input sequence's length, like [20,17,...]\n
`weight`    : weight for TCS loss

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│TCSLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function DNN_Batch_TCS(p::Variable{Array{T}},
                       seqlabels::Vector,
                       inputlens;
                       background::Int=1,
                       foreground::Int=2,
                       weight=1.0) where T
    batchsize = length(seqlabels)
    loglikely = zeros(T, batchsize)
    I, F = indexbounds(inputlens)
    r = zero(ᵛ(p))

    Threads.@threads for b = 1:batchsize
        span = I[b]:F[b]
        r[:,span], loglikely[b] = TCS(p.value[:,span], seqlabels[b], background=background, foreground=foreground)
        loglikely[b] /= length(seqlabels[b]) * 3 + 1
    end

    if p.backprop
        function DNN_Batch_TCS_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, DNN_Batch_TCS_Backward)
    end
    return sum(loglikely)/batchsize
end

"""
    RNN_Batch_TCS(p::Variable{Array{T}},
                  seqlabels::Vector,
                  inputlens;
                  background::Int=1,
                  foreground::Int=2,
                  weight=1.0) where T

a batch of padded input sequence is processed by neural networks into `p`

# Inputs
`p`         : 3-D Variable with shape (featdims,timesteps,batchsize), probability or weighted probability\n
`seqlabels` : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`inputlens` : records each input sequence's length, like [20,17,...]\n
`weight`    : weight for TCS loss

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│TCSLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function RNN_Batch_TCS(p::Variable{Array{T}},
                       seqlabels::Vector,
                       inputlens;
                       background::Int=1,
                       foreground::Int=2,
                       weight=1.0) where T
    batchsize = length(seqlabels)
    loglikely = zeros(T, batchsize)
    r = zero(ᵛ(p))

    Threads.@threads for b = 1:batchsize
        Tᵇ = inputlens[b]
        Lᵇ = length(seqlabels[b])
        r[:,1:Tᵇ,b], loglikely[b] = TCS(p.value[:,1:Tᵇ,b], seqlabels[b], background=background, foreground=foreground)
        loglikely[b] /= Lᵇ * 3 + 1
    end

    if p.backprop
        function RNN_Batch_TCS_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, RNN_Batch_TCS_Backward)
    end
    return sum(loglikely)/batchsize
end

"""
    CRNN_Batch_TCS(p::Variable{Array{T}},
                   seqlabels::Vector;
                   background::Int=1,
                   foreground::Int=2,
                   weight=1.0) where T

a batch of padded input sequence is processed by neural networks into `p`

# Main Inputs
`p`            : 3-D Variable with shape (featdims,timesteps,batchsize), probability or weighted probability\n
`seqlabels`    : a batch of sequential labels, like [[i,j,k],[x,y],...]\n
`weight`       : weight for TCS loss

# Structure
                   ┌───┐
                   │ │ │
                   │ W ├──►──┐
                   │ │ │     │
                   └───┘     │
    ┌───┐          ┌───┐     │     ┌───┐
    │ │ │ softmax  │ │ │   ┌─┴─┐   │ │ │   ┌───────┐
    │ X ├─────────►│ Y ├──►│ × ├──►│ P ├──►│TCSLOSS│◄── (seqLabel)
    │ │ │          │ │ │   └───┘   │ │ │   └───┬───┘
    └─┬─┘          └─┬─┘           └───┘       │
      │              │+                        ▼
    ┌─┴─┐            ▼                       ┌─┴─┐
    │ │ │          ┌─┴─┐ -                   │ │ │
    │ δ │◄─────────┤ - │─────────◄───────────┤ r │
    │ │ │          └───┘                     │ │ │
    └───┘                                    └───┘
"""
function CRNN_Batch_TCS(p::Variable{Array{T}},
                        seqlabels::Vector;
                        background::Int=1,
                        foreground::Int=2,
                        weight=1.0) where T
    featdims, timesteps, batchsize = size(p)
    loglikely = zeros(T, batchsize)
    r = zero(ᵛ(p))

    Threads.@threads for b = 1:batchsize
        r[:,:,b], loglikely[b] = TCS(p.value[:,:,b], seqlabels[b], background=background, foreground=foreground)
        loglikely[b] /= length(seqlabels[b]) * 3 + 1
    end

    if p.backprop
        function CRNN_Batch_TCS_Backward()
            if need2computeδ!(p)
                if weight==1.0
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T))
                else
                    δ(p) .-= r ./ (ᵛ(p) .+ eps(T)) .* weight
                end
            end
        end
        push!(graph.backward, CRNN_Batch_TCS_Backward)
    end
    return sum(loglikely)/batchsize
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
# Largely inspired by CPLEXCP.jl.
# Main differences:
# - no support for quadratic functions (for now?)
# - no support for adding int/bool constraint on variables (only at creation)
# TODOs:
# - objective
# - var/cons name
# - parameters: number of threads, time limit, raw
# - looks like JaCoP is as immutable as possible: looks like it's impossible to change the name of a variable, e.g.

@enum(VariableType, CONTINUOUS, BINARY, INTEGER, SET, CIRCUIT)

mutable struct VariableInfo
    index::MOI.VariableIndex
    variable::Variable
    name::String # TODO: remove me?
    type::VariableType

    # Variable bounds: either integers or floats. Needed for is_valid on VariableIndex
    # and bound constraints.
    lb::Union{Nothing, Integer, Float64}
    ub::Union{Nothing, Integer, Float64}
end

function VariableInfo(index::MOI.VariableIndex, variable::Variable)
    return VariableInfo(
        index,
        variable,
        "",
        INTEGER,
        nothing,
        nothing,
    )
end

mutable struct ConstraintInfo
    # Only necessary information to access an existing constraint, delete it when needed.
    index::MOI.ConstraintIndex
    constraint::Union{Constraint, Nothing}
    f::Union{MOI.AbstractScalarFunction, MOI.AbstractVectorFunction}
    set::MOI.AbstractSet
    name::String
end

function ConstraintInfo(
    index::MOI.ConstraintIndex,
    constraint::Union{Constraint, Nothing},
    f::Union{MOI.AbstractScalarFunction, MOI.AbstractVectorFunction},
    set::MOI.AbstractSet,
)
    return ConstraintInfo(index, constraint, f, set, "")
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    # The low-level JaCoP store.
    inner::Store

    # The model name.
    name::String

    # A mapping from the MOI.VariableIndex to the CPLEX variable object.
    # VariableInfo also stores some additional fields like the type of variable.
    variable_info::CleverDicts.CleverDict{MOI.VariableIndex, VariableInfo}

    # A mapping from the MOI.ConstraintIndex to the CPLEX variable object.
    # VariableInfo also stores some additional fields like the type of variable.
    constraint_info::Dict{MOI.ConstraintIndex, ConstraintInfo}

    # # Memorise the objective sense and the function separately, as the Concert
    # # API forces to give both at the same time.
    # objective_sense::MOI.OptimizationSense
    # objective_function::Union{Nothing, MOI.AbstractScalarFunction}
    # objective_function_cp::Union{Nothing, NumExpr}
    # objective_cp::Union{Nothing, IloObjective}

    # # Cache parts of a solution.
    # cached_solution_state::Union{Nothing, Bool}

    # # Mappings from variable and constraint names to their indices. These are
    # # lazily built on-demand, so most of the time, they are `nothing`.
    # # The solver's functionality is not useful in this case, as it can only
    # # handle integer variables. Moreover, bound constraints do not have names
    # # for the solver.
    # name_to_variable::Union{Nothing, Dict{String, MOI.VariableIndex}}
    # name_to_constraint::Union{Nothing, Dict{String, MOI.ConstraintIndex}}

    """
        Optimizer()

    Create a new Optimizer object.
    """
    function Optimizer()
        model = new()
        model.inner = Store()

        model.variable_info =
            CleverDicts.CleverDict{MOI.VariableIndex, VariableInfo}()
        model.constraint_info = Dict{MOI.ConstraintIndex, ConstraintInfo}()

        # model.objective_sense = MOI.FEASIBILITY_SENSE
        # model.objective_function = nothing
        # model.objective_function_cp = nothing
        # model.objective_cp = nothing

        # model.cached_solution_state = nothing
        # model.callback_state = CB_NONE

        MOI.empty!(model)
        return model
    end
end

Base.show(io::IO, model::Optimizer) = show(io, model.inner)

function MOI.empty!(model::Optimizer)
    model.inner = Store()
    model.name = ""
    empty!(model.variable_info)
    empty!(model.constraint_info)

    # model.objective_sense = MOI.FEASIBILITY_SENSE
    # model.objective_function = nothing
    # model.objective_function_cp = nothing
    # model.objective_cp = nothing

    # model.cached_solution_state = nothing
    # model.name_to_variable = nothing
    # model.name_to_constraint = nothing
    return
end

function MOI.is_empty(model::Optimizer)
    !isempty(model.name) && return false
    !isempty(model.variable_info) && return false
    !isempty(model.constraint_info) && return false
    # model.objective_sense != MOI.FEASIBILITY_SENSE && return false
    # model.objective_function !== nothing && return false
    # model.objective_function_cp !== nothing && return false
    # model.objective_cp !== nothing && return false
    # model.name_to_variable !== nothing && return false
    # model.name_to_constraint !== nothing && return false
    # model.cached_solution_state !== nothing && return false
    return true
end

MOI.get(::Optimizer, ::MOI.SolverName) = "JaCoP"

## Types of objectives and constraints that are supported.

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{F},
) where {
    F <: Union{
        MOI.VariableIndex,
        MOI.ScalarAffineFunction{Float64},
    },
}
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{F},
) where {
    T <: Union{Int32, Float64},
    F <: Union{
        MOI.EqualTo{T},
        MOI.LessThan{T},
        MOI.GreaterThan{T},
        MOI.Interval{T},
    },
}
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{T}},
    ::Type{F},
) where {
    T <: Union{Int32, Float64},
    F <:
    Union{MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T}},
    # No interval!
}
    return true
end

MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{MOI.VariableIndex}) = true
function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintName,
    ::Type{<:MOI.ConstraintIndex},
)
    return true
end

MOI.supports(::Optimizer, ::MOI.Name) = true
# MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = true
# MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
# MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
# MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    return MOI.Utilities.default_copy_to(dest, src)
end

function MOI.get(::Optimizer, ::MOI.ListOfVariableAttributesSet)
    return MOI.AbstractVariableAttribute[MOI.VariableName()]
end

function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
    attributes = Any[MOI.ObjectiveSense()]
    typ = MOI.get(model, MOI.ObjectiveFunctionType())
    if typ !== nothing
        push!(attributes, MOI.ObjectiveFunction{typ}())
    end
    if MOI.get(model, MOI.Name()) != ""
        push!(attributes, MOI.Name())
    end
    return attributes
end

function MOI.get(::Optimizer, ::MOI.ListOfConstraintAttributesSet)
    return MOI.AbstractConstraintAttribute[MOI.ConstraintName()]
end


text:
# Tests for Atomistic.jl

using Atomistic
using PeriodicTable
using Test
using Unitful
using UnitfulAtomic

@testset "Atomistic.jl" begin
    @testset "Unit Tests" begin
        include("unit_tests_ab_integration.jl")
        include("unit_tests_nbs_integration.jl")
        include("unit_tests_dftk_integration.jl")
    end
    @testset "Integration Tests" begin
        include("integration_test_ab_nbs.jl")
        include("integration_test_ab_nbs_ip.jl")
        include("integration_test_ab_nbs_dftk.jl")
    end
end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
""""
Enhanced string literals

String literals with Swift-like format
C and Python like formatting
LaTex, Emoji, HTML, and Unicode names

Copyright 2016-2018 Gandalf Software, Inc., Scott P. Jones
Portions originally from code in Julia, copyright Julia contributors
Licensed under MIT License, see LICENSE.md
"""
module StringLiterals

using Format
#export FormatSpec, FormatExpr, printfmt, printfmtln, format, generate_formatter
#export pyfmt, cfmt, fmt, fmt_default, fmt_default!, reset!, defaultSpec

using StrTables, LaTeX_Entities, HTML_Entities, Unicode_Entities, Emoji_Entities

export @f_str, @F_str, @sinterpolate, @pr_str, @PR_str, @pr, @PR
export s_unescape_string, s_escape_string, s_print_unescaped, s_print_escaped

@static if VERSION < v"0.7.0-DEV"
    const _parse = parse
    const _ParseError = ParseError
    _sprint(f, s) = sprint(endof(s), f, s)
    _sprint(f, s, c) = sprint(endof(s), f, s, c)
    const outhex = hex
else
    const _parse = Meta.parse
    const _ParseError = Base.Meta.ParseError
    _sprint(f, s) = sprint(f, s; sizehint=lastindex(s))
    _sprint(f, s, c) = sprint(f, s, c; sizehint=lastindex(s))
    outhex(c, p=1) = string(UInt32(c), base=16, pad=p)
end

@static if VERSION >= v"1.0.0"
    const next = iterate
    start(iter) = 1
    done(iter, state) = iterate(iter, state) === nothing
end

"""
String macro with more Swift-like syntax, plus support for emojis and LaTeX names
"""
macro f_str(str) ; s_interp_parse(false, str) ; end

"""
String macro with more Swift-like syntax, plus support for emojis and LaTeX names, also legacy
"""
macro F_str(str) ; s_interp_parse(true, str) ; end

"""
String macros that calls print directly
"""
macro pr_str(str) ; s_print(false, str) ; end
macro PR_str(str) ; s_print(true, str) ; end

"""
Interpolates one or more strings using more Swift-like syntax
julia> x = "World"; @sinterpolate "Hello \\(x)"
"Hello World"
"""
macro sinterpolate(args...) ; s_interp_parse(false, args...) ; end

throw_arg_err(msg) = throw(ArgumentError(msg))
throw_arg_err(msg, s) = throw_arg_err(string(msg, repr(s)))

"""
Handle Unicode character constant, of form \\u{<hexdigits>}
"""
function s_parse_unicode(io, s,  i)
    done(s,i) &&
        throw_arg_err("Incomplete \\u{...} in ", s)
    c, i = next(s, i)
    c != '{' &&
        throw_arg_err("\\u missing opening { in ", s)
    done(s,i) &&
        throw_arg_err("Incomplete \\u{...} in ", s)
    beg = i
    c, i = next(s, i)
    n::UInt32 = 0
    k = 0
    while c != '}'
        done(s, i) &&
            throw_arg_err("\\u{ missing closing } in ", s)
        (k += 1) > 6 &&
            throw_arg_err("Unicode constant too long in ", s)
        n = n<<4 + c - ('0' <= c <= '9' ? '0' :
                        'a' <= c <= 'f' ? 'a' - 10 :
                        'A' <= c <= 'F' ? 'A' - 10 :
                        throw_arg_err("\\u missing closing } in ", s))
        c, i = next(s,i)
    end
    k == 0 &&
        throw_arg_err("\\u{} has no hex digits in ", s)
    ((0x0d800 <= n <= 0x0dfff) || n > 0x10ffff) &&
        throw_arg_err("Invalid Unicode character constant ", s[beg-3:i-1])
    print(io, Char(n))
    i
end

"""
Handle Emoji character, of form \\:name:
"""
function s_parse_emoji(io, s,  i)
    beg = i # start location
    c, i = next(s, i)
    while c != ':'
        done(s, i) &&
            throw_arg_err("\\: missing closing : in ", s)
        c, i = next(s, i)
    end
    emojistr = lookupname(Emoji_Entities.default, s[beg:i-2])
    emojistr == "" &&
        throw_arg_err("Invalid Emoji name in ", s)
    print(io, emojistr)
    i
end

"""
Handle LaTeX character/string, of form \\<name>
"""
function s_parse_latex(io, s,  i)
    beg = i # start location
    c, i = next(s, i)
    while c != '>'
        done(s, i) &&
            throw_arg_err("\\< missing closing > in ", s)
        c, i = next(s, i)
    end
    latexstr = lookupname(LaTeX_Entities.default, s[beg:i-2])
    latexstr == "" &&
        throw_arg_err("Invalid LaTeX name in ", s)
    print(io, latexstr)
    i
end

"""
Handle HTML character/string, of form \\&name;
"""
function s_parse_html(io, s,  i)
    beg = i # start location
    c, i = next(s, i)
    while c != ';'
        done(s, i) &&
            throw_arg_err("\\& missing ending ; in ", s)
        c, i = next(s, i)
    end
    htmlstr = lookupname(HTML_Entities.default, s[beg:i-2])
    htmlstr == "" &&
        throw_arg_err("Invalid HTML name in ", s)
    print(io, htmlstr)
    i
end

"""
Handle Unicode name, of form \\N{name}, from Python
"""
function s_parse_uniname(io, s,  i)
    done(s, i) &&
        throw_arg_err("\\N incomplete in ", s)
    c, i = next(s, i)
    c != '{' &&
        throw_arg_err("\\N missing initial { in ", s)
    done(s, i) &&
        throw_arg_err("\\N{ incomplete in ", s)
    beg = i # start location
    c, i = next(s, i)
    while c != '}'
        done(s, i) &&
            throw_arg_err("\\N{ missing closing } in ", s)
        c, i = next(s, i)
    end
    unistr = lookupname(Unicode_Entities.default, s[beg:i-2])
    unistr == "" &&
        throw_arg_err("Invalid Unicode name in ", s)
    print(io, unistr)
    i
end

"""
String interpolation parsing, allow legacy \$, \\xHH, \\uHHHH, \\UHHHHHHHH
"""
s_print_unescaped_legacy(io, s::AbstractString) = s_print_unescaped(io, s, true)

"""
String interpolation parsing
Based on code resurrected from Julia base:
https://github.com/JuliaLang/julia/blob/deab8eabd7089e2699a8f3a9598177b62cbb1733/base/string.jl
"""
function s_print_unescaped(io, s::AbstractString, flg::Bool=false)
    i = start(s)
    while !done(s,i)
        c, i = next(s,i)
        if !done(s,i) && c == '\\'
            c, i = next(s,i)
            if c == 'u'
                i = flg ? s_parse_unicode_legacy(io, s, i, c) : s_parse_unicode(io, s, i)
            elseif c == ':'     # Emoji
                i = s_parse_emoji(io, s, i)
            elseif c == '&'     # HTML
                i = s_parse_html(io, s, i)
            elseif c == '<'     # LaTeX
                i = s_parse_latex(io, s, i)
            elseif c == 'N'     # Unicode name
                i = s_parse_uniname(io, s, i)
            elseif flg && c == 'U'
                i = s_parse_unicode_legacy(io, s, i, c)
            elseif flg && c == 'x' # hex byte
                i = s_parse_hex_legacy(io, s, i)
            else
                c = (c == '0' ? '\0' :
                     c == '$' ? '$'  :
                     c == '"' ? '"'  :
                     c == '\'' ? '\'' :
                     c == '\\' ? '\\' :
                     c == 'a' ? '\a' :
                     c == 'b' ? '\b' :
                     c == 't' ? '\t' :
                     c == 'n' ? '\n' :
                     c == 'v' ? '\v' :
                     c == 'f' ? '\f' :
                     c == 'r' ? '\r' :
                     c == 'e' ? '\e' :
                     throw_arg_err(string("Invalid \\",c," sequence in "), s))
                write(io, UInt8(c))
            end
        else
            print(io, c)
        end
    end
end

function s_parse_unicode_legacy(io, s, i, c)
    done(s, i) &&
        throw_arg_err(string("\\", c, " used with no following hex digits"))
    beg = i
    m = (c == 'u' ? 4 : 8)
    if s[i] == '{'
        m == 4 || throw_arg_err("{ only allowed with \\u")
        return s_parse_unicode(io, s, i)
    end
    n = k = 0
    while (k += 1) <= m && !done(s,i)
        c, j = next(s,i)
        n = '0' <= c <= '9' ? n<<4 + c-'0' :
            'a' <= c <= 'f' ? n<<4 + c-'a'+10 :
            'A' <= c <= 'F' ? n<<4 + c-'A'+10 : break
        i = j
    end
    k == 1 &&
        throw_arg_err(string("\\", c, " used with no following hex digits"))
    ((0x0d800 <= n <= 0x0dfff) || n > 0x10ffff) &&
        throw_arg_err("Invalid Unicode character constant ", s[beg-2:i-1])
    print(io, Char(n))
    i
end

hexerr() = throw_arg_err("\\x used with no following hex digits")
function s_parse_hex_legacy(io, s, i)
    done(s,i) && hexerr()
    c, i = next(s,i)
    n = '0' <= c <= '9' ? c-'0' :
        'a' <= c <= 'f' ? c-'a'+10 :
        'A' <= c <= 'F' ? c-'A'+10 : hexerr()
    if done(s,i)
        write(io, UInt8(n))
        return i
    end
    c, j = next(s,i)
    n = '0' <= c <= '9' ? n<<4 + c-'0' :
        'a' <= c <= 'f' ? n<<4 + c-'a'+10 :
        'A' <= c <= 'F' ? n<<4 + c-'A'+10 : (j = i; n)
    write(io, UInt8(n))
    j
end

s_unescape_string(s::AbstractString) = _sprint(s_print_unescaped, s)

function s_print_escaped(io, s::AbstractString, esc::Union{AbstractString, Char})
    i = start(s)
    while !done(s,i)
        c, i = next(s, i)
        c == '\0'       ? print(io, "\\0") :
        c == '\e'       ? print(io, "\\e") :
        c == '\\'       ? print(io, "\\\\") :
        c in esc        ? print(io, '\\', c) :
        '\a' <= c <= '\r' ? print(io, '\\', "abtnvfr"[Int(c)-6]) :
        isprint(c)      ? print(io, c) :
                          print(io, "\\u{", outhex(c), "}")
    end
end

s_escape_string(s::AbstractString) = _sprint(s_print_escaped, s, '\"')

const ByteStr = String

s_print(flg::Bool, s::AbstractString) = s_print(flg, s, flg ? s_unescape_str : s_unescape_legacy)
function s_print(flg::Bool, s::AbstractString, unescape::Function)
    sx = s_interp_parse_vec(flg, s, unescape)
    (length(sx) == 1 && isa(sx[1], String)
     ? Expr(:call, :print, sx[1])
     : Expr(:call, :print, sx...))
end

function s_interp_parse(flg::Bool, s::AbstractString, unescape::Function, p::Function)
    sx = s_interp_parse_vec(flg, s, unescape)
    length(sx) == 1 && isa(sx[1], String) ? sx[1] : Expr(:call, :sprint, p, sx...)
end

function s_interp_parse_vec(flg::Bool, s::AbstractString, unescape::Function)
    sx = []
    i = j = start(s)
    while !done(s, j)
        c, k = next(s, j)
        if c == '\\' && !done(s, k)
            if s[k] == '('
                # Handle interpolation
                isempty(s[i:j-1]) ||
                    push!(sx, unescape(s[i:j-1]))
                ex, j = _parse(s, k, greedy=false)
                isa(ex, Expr) && (ex.head === :continue) &&
                    throw(_ParseError("Incomplete expression"))
                push!(sx, esc(ex))
                i = j
            elseif s[k] == '%'
                # Move past \\, c should point to '%'
                c, k = next(s, k)
                done(s, k) && throw(_ParseError("Incomplete % expression"))
                # Handle interpolation
                if !isempty(s[i:j-1])
                    push!(sx, unescape(s[i:j-1]))
                end
                if s[k] == '('
                    # Need to find end to parse to
                    _, j = _parse(s, k, greedy=false)
                    # This is a bit hacky, and probably doesn't perform as well as it could,
                    # but it works! Same below.
                    str = string("(StringLiterals.fmt", s[k:j-1], ')')
                else
                    # Move past %, c should point to letter
                    beg = k
                    while true
                        c, k = next(s, k)
                        done(s, k) &&
                            throw(_ParseError("Incomplete % expression"))
                        s[k] == '(' && break
                    end
                    _, j = _parse(s, k, greedy=false)
                    str = string("(StringLiterals.cfmt(\"", s[beg-1:k-1], "\",", s[k+1:j-1], ')')
                end
                ex, _ = _parse(str, 1, greedy=false)
                isa(ex, Expr) && (ex.head === :continue) &&
                    throw(_ParseError("Incomplete expression"))
                push!(sx, esc(ex))
                i = j
            elseif s[k] == '{'
                # Move past \\, c should point to '{'
                c, k = next(s, k)
                done(s, k) &&
                    throw(_ParseError("Incomplete {...} Python format expression"))
                # Handle interpolation
                isempty(s[i:j-1]) ||
                    push!(sx, unescape(s[i:j-1]))
                beg = k # start location
                c, k = next(s, k)
                while c != '}'
                    done(s, k) &&
                        throw_arg_err("\\{ missing closing } in ", s)
                    c, k = next(s, k)
                end
                done(s, k) &&
                    throw(_ParseError("Missing (expr) in Python format expression"))
                c, k = next(s, k)
                c != '(' &&
                    throw(_ParseError(string("Missing (expr) in Python format expression: ", c)))
                # Need to find end to parse to
                _, j = _parse(s, k-1, greedy=false)
                # This is a bit hacky, and probably doesn't perform as well as it could,
                # but it works! Same below.
                str = string("(StringLiterals.pyfmt(\"", s[beg:k-3], "\",", s[k:j-1], ')')
                ex, _ = _parse(str, 1, greedy=false)
                isa(ex, Expr) && (ex.head === :continue) &&
                    throw(_ParseError("Incomplete expression"))
                push!(sx, esc(ex))
                i = j
            elseif flg && s[k] == '$'
                isempty(s[i:j-1]) ||
                    push!(sx, unescape(s[i:j-1]))
                i = k
                # Move past \\, c should point to '$'
                c, j = next(s, k)
            else
                j = k
            end
        elseif flg && c == '$'
            isempty(s[i:j-1]) ||
                push!(sx, unescape(s[i:j-1]))
            ex, j = _parse(s, k, greedy=false)
            isa(ex,Expr) && ex.head === :continue &&
                throw(_ParseError("incomplete expression"))
            push!(sx, esc(ex))
            i = j
        else
            j = k
        end
    end
    isempty(s[i:end]) ||
        push!(sx, unescape(s[i:j-1]))
    sx
end

function s_unescape_str(s)
    s = s_unescape_string(s)
    isvalid(String, s) ? s : throw_arg_err("Invalid UTF-8 sequence")
end
function s_unescape_legacy(s)
    s = _sprint(s_print_unescaped_legacy, s)
    isvalid(String, s) ? s : throw_arg_err("Invalid UTF-8 sequence")
end

s_interp_parse(flg::Bool, s::AbstractString, u::Function) = s_interp_parse(flg, s, u, print)
s_interp_parse(flg::Bool, s::AbstractString) =
    s_interp_parse(flg, s, flg ? s_unescape_legacy : s_unescape_str)

end # module


text:
function lab_21(r)
    go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    horizontal = 0
    vertical = 0

    border_ver = false
    border_hor = false
    side = Ost

    while !isborder(r, Nord)
        while move_if_possible!(r, side)
            horizontal, border_hor = find_special!(r, horizontal, border_hor, Nord)
            vertical, border_ver = find_special!(r, vertical, border_ver, West)
        end

        horizontal, border_hor = find_special!(r, horizontal, border_hor, Nord)
        vertical, border_ver = find_special!(r, vertical, border_ver, West)

        side = inverse_side(side)

        move!(r, North)
    end

    return horizontal, vertical
end

function find_special!(r::Robot, my_ans::Int, border_now::Bool, side::HorizonSide)
    if isborder(r, side)
        border_now = true
    end
    if !isborder(r, side) && border_now
        border_now = false
        my_ans += 1
    end
    return my_ans, border_now
end


text:
using ..Internal

abstract FunctionSystem

@mustimplement evaluate{T <: Number}(basis::FunctionSystem, n::Integer, x::T)

"The base type for systems of (orthogonal) polynomials"
abstract PolynomialSystem <: FunctionSystem

@mustimplement recurrence_coeff(basis::PolynomialSystem, i::Integer)

@mustimplement issymmetric(basis::PolynomialSystem)

function rc_array(basis::PolynomialSystem, n::Integer)
  arr = [recurrence_coeff(basis, k) for k in 0:n]
  flip_tuple_array(arr)
end

function rc_array_monic(basis::PolynomialSystem, n::Integer)
  r = rc_array(basis::PolynomialSystem, n::Integer)
  # extract columns
  a = -r[1]
  b = r[3][2:end]
  c = r[2][:]

  # convert to monic polynomials
  α = a ./ c
  β = b ./ (c[1:end-1] .* c[2:end])
  return α, β
end

"Evaluate polynomials at `x` up to degree `n-1`"
function evaluate{T <: Number}(basis::PolynomialSystem, n::Integer, x::T)
  y = zeros(T, n+1)
  y[1] = 1
  for k=0:n-1
    a, b, c = recurrence_coeff(basis, k)
    y[k+2] = (a + b * x) * y[k+1]
    if k>0
      y[k+2] = y[k+2] - c * y[k]
    end
  end
  return y
end

"The stochastic Hermite polynomials"
immutable HermitePolynomials <: PolynomialSystem; end

recurrence_coeff(::HermitePolynomials, k::Integer) =
  0, 1, k

issymmetric(::HermitePolynomials) = true

"The Legendre polynomials"
immutable LegendrePolynomials <: PolynomialSystem; end

recurrence_coeff(::LegendrePolynomials, k::Integer) =
  0, (2*k+1) // (k+1), k // (k+1)

issymmetric(::LegendrePolynomials) = true

"The Laguerre polynomials"
immutable LaguerrePolynomials{T<:Real} <: PolynomialSystem
  α::T
end

recurrence_coeff(L::LaguerrePolynomials, k::Integer) =
  (2k + 1 + L.α) / (k+1),
  -1 / (k+1),
  (k + L.α) / (k+1)

issymmetric(::LaguerrePolynomials) = false

immutable ChebyshevTPolynomials <: PolynomialSystem; end

recurrence_coeff(basis::ChebyshevTPolynomials, k::Integer) =
  0, (k==0 ? 1 : 2), 1

issymmetric(::ChebyshevTPolynomials) = true

immutable ChebyshevUPolynomials <: PolynomialSystem; end

recurrence_coeff(::ChebyshevUPolynomials, k::Integer) =
  0, 2, 1

issymmetric(::ChebyshevUPolynomials) = true

import Alea.Quadrature.gauss_rule
function gauss_rule(basis::PolynomialSystem, n::Integer)
  α, β = rc_array_monic(basis, n-1)

  gauss_rule((α, β))
end

export
  PolynomialSystem,
  evaluate,
  recurrence_coeff,
  issymmetric,
  HermitePolynomials,
  LegendrePolynomials,
  LaguerrePolynomials,
  ChebyshevTPolynomials,
  ChebyshevUPolynomials