text:
using FreeTypeAbstraction: height_insensitive_boundingbox

"""
Calculates the exact boundingbox of a Scene/Plot, without considering any transformation
"""
raw_boundingbox(x::Atomic) = data_limits(x)

rootparent(x) = rootparent(parent(x))
rootparent(x::Scene) = x

# function raw_boundingbox(x::Annotations)
#     bb = raw_boundingbox(x.plots)
#     inv(modelmatrix(rootparent(x))) * bb
# end

raw_boundingbox(x::Combined) = raw_boundingbox(x.plots)
boundingbox(x) = raw_boundingbox(x)

function combined_modelmatrix(x)
    m = Mat4f0(I)
    while true
        m = modelmatrix(x) * m
        if parent(x) !== nothing && parent(x) isa Combined
            x = parent(x)
        else
            break
        end
    end
    return m
end

function modelmatrix(x)
    t = transformation(x)
    transformationmatrix(t.translation[], t.scale[], t.rotation[])
end

function boundingbox(x::Atomic)
    bb = raw_boundingbox(x)
    return combined_modelmatrix(x) * bb
end

boundingbox(scene::Scene) = raw_boundingbox(scene)
function raw_boundingbox(scene::Scene)
    if scene[OldAxis] !== nothing
        return raw_boundingbox(scene[OldAxis])
    elseif scene.limits[] !== automatic
        return scene_limits(scene)
    elseif cameracontrols(scene) == EmptyCamera()
        # Empty camera means this is a parent scene that itself doesn't display anything
        return raw_boundingbox(scene.children)
    else
        plots = plots_from_camera(scene)
        children = filter(scene.children) do child
            child.camera == scene.camera
        end
        return raw_boundingbox([plots; children])
    end
end

function raw_boundingbox(plots::Vector)
    isempty(plots) && return FRect3D()
    plot_idx = iterate(plots)
    bb = FRect3D()
    while plot_idx !== nothing
        plot, idx = plot_idx
        plot_idx = iterate(plots, idx)
        # isvisible(plot) || continue
        bb2 = boundingbox(plot)
        isfinite_rect(bb) || (bb = bb2)
        isfinite_rect(bb2) || continue
        bb = union(bb, bb2)
    end
    return bb
end

function project_widths(matrix, vec)
    pr = project(matrix, vec)
    zero = project(matrix, zeros(typeof(vec)))
    return pr - zero
end

function rotate_bbox(bb::FRect3D, rot)
    points = decompose(Point3f0, bb)
    FRect3D(Ref(rot) .* points)
end

function boundingbox(x::Text, text::String, position::VecTypes)

    glyphlayout = x._glyphlayout[]

    pos = to_ndim(Point3f0, x.position[], 0)
    rot = convert_attribute(x.rotation[], key"rotation"())

    if x.space[] == :data
        data_text_boundingbox(x[1][], glyphlayout, rot, pos)
    elseif x.space[] == :screen
        data_limits(x)
    end
end

function data_text_boundingbox(string::String, glyphlayout::Glyphlayout, rotation::Quaternion, position::Point3f0)
    bb = FRect3D()
    glyphorigins, glyphbbs = glyphlayout.origins, glyphlayout.bboxes

    for (char, charo, glyphbb) in zip(string, glyphorigins, glyphbbs)
        # ignore line breaks
        # char in ('\r', '\n') && continue

        charbb = rotate_bbox(FRect3D(glyphbb), rotation) + charo + position
        if !isfinite_rect(bb)
            bb = charbb
        else
            bb = union(bb, charbb)
        end
    end

    bb
end

function boundingbox(x::Text, texts::AbstractArray, positions::AbstractArray)

    layouts = x._glyphlayout[]
    rotations = convert_attribute(x.rotation[], key"rotation"())
    positions_3d = to_ndim.(Point3f0, positions, 0)

    if x.space[] == :data
        bb = FRect3D()
        broadcast_foreach(texts, positions_3d, layouts, rotations) do text, pos, layout, rot
            if !isfinite_rect(bb)
                bb = data_text_boundingbox(text, layout, rot, pos)
            else
                bb = union(bb, data_text_boundingbox(text, layout, rot, pos))
            end
        end
        bb
    elseif x.space[] == :screen
        data_limits(x)
    end
end

function boundingbox(x::Text)

    # TODO: this is only necessary because of the "text in text" recipe weirdness
    if !isempty(x.plots)
        # the "real" text is inside the outer text
        x = x.plots[1]
    end

    boundingbox(x, to_value(x[1]), to_value(x[:position]))
end

#     if x.space[] == :data
#         if x[1][] isa AbstractArray
#             bb = FRect3D()
#             for
#         else
#             boundingbox(x, to_value(x[1]))
#         end
#     elseif x.space[] == :screen
#         data_limits(x)
#     end
# end

# boundingbox(x::Text) = boundingbox(x, to_value(x[1]))

# function boundingbox(
#         text::String, position, textsize;
#         font = "default", align = (:left, :bottom), rotation = 0.0
#     )
#     return boundingbox(
#         text, position, textsize,
#         to_font(font), to_align(align), to_rotation(rotation)
#     )
# end

"""
Calculate an approximation of a tight rectangle around a 2D rectangle rotated by `angle` radians.
This is not perfect but works well enough. Check an A vs X to see the difference.
"""
function rotatedrect(rect::Rect{2}, angle)
    ox, oy = rect.origin
    wx, wy = rect.widths
    points = @SMatrix([
        ox oy;
        ox oy+wy;
        ox+wx oy;
        ox+wx oy+wy;
    ])
    mrot = @SMatrix([
        cos(angle) -sin(angle);
        sin(angle) cos(angle);
    ])
    rotated = mrot * points'

    rmins = minimum(rotated, dims = 2)
    rmaxs = maximum(rotated, dims = 2)

    return Rect2D(rmins..., (rmaxs .- rmins)...)
end

function quaternion_to_2d_angle(quat)
    # this assumes that the quaternion was calculated from a simple 2d rotation as well
    return 2acos(quat[4]) * (signbit(quat[1]) ? -1 : 1)
end

function boundingbox(
        text::String, position, textsize, font,
        align, rotation, model, justification, lineheight;
        # use the font's ascenders and descenders for the bounding box
        # this means that a string's boundingbox doesn't change in the vertical
        # dimension when characters change (for example numbers during an animation)
        # this is not wanted in most cases because of the jitter it creates when
        # the boundingbox slightly changes size in each frame (in MakieLayout mostly)
        use_vertical_dimensions_from_font = true
    )
    atlas = get_texture_atlas()
    N = length(text)
    ctext_state1 = iterate(text)
    ctext_state1 === nothing && return FRect3D()

    # call the layouting algorithm to find out where all the glyphs end up
    # this is kind of a doubling, maybe it could be avoided if at creation all
    # positions would be populated in the text object, but that seems convoluted
    if position isa VecTypes
        glyph_offsets = layout_text(text, textsize, font, align,
            rotation, model, justification, lineheight)
        position = [go .+ to_ndim(Point3f0, position, 0) for go in glyph_offsets]
    end

    bbox = Ref(FRect3D())            # use ref to avoid julia#15276
    ctext_state = Ref{Union{Nothing,typeof(ctext_state1)}}(ctext_state1)

    let position = position
        broadcast_foreach(1:N, rotation, font, textsize) do i, rotation, font, scale
            c, text_state = ctext_state[]::typeof(ctext_state1)
            ctext_state[] = iterate(text, text_state)

            if !(c in ('\r', '\n'))
                bb_unitspace = if use_vertical_dimensions_from_font
                    height_insensitive_boundingbox(
                        FreeTypeAbstraction.get_extent(font, c), font)
                else
                    inkboundingbox(FreeTypeAbstraction.get_extent(font, c))
                end

                scaled_bb = bb_unitspace * scale

                # TODO this only works in 2d
                rot_2d_radians = quaternion_to_2d_angle(rotation)
                rotated_bb = rotatedrect(scaled_bb, rot_2d_radians)

                # bb = rectdiv(bb, 1.5)
                shifted_bb = FRect3D(rotated_bb) + position[i]
                bboxc = bbox[]
                if !isfinite_rect(bboxc)
                    bbox[] = shifted_bb
                else
                    bbox[] = union(bboxc, shifted_bb)
                end
            end
        end
    end
    return bbox[]
end


text:
dataref_spec =
[
	(:data, AbstractString, nothing),
	(:field, Union{AbstractString, Vector, Dict{Any,Any}}, nothing),
	(:fields, Union{AbstractString, Vector, Dict{Any,Any}}, nothing),
	(:sort, Union{Bool, Dict{Any,Any}}, nothing)
]

primitivefactory(:VegaDataRef, dataref_spec)

ext: jl

max_stars_repo_name: jpthor/vega


text:
__precompile__()

module SWSs

import Base: filter!, filter, getindex
using Base.Dates

using AutoHashEquals
using Geodesy

import SphericalGeom

export
    AbstractSWS,
    SWS,
    SWSTime,
    azimuth,
    backazimuth,
    midpoint,
    read_sheba,
    read_sheba_time,
    surf_dist,
    write_sheba

#=
    Types
=#
abstract type AbstractSWS end

"""
    SWS

Type containing information about a shear wave splitting measurement.
"""
@auto_hash_equals immutable SWS <: AbstractSWS
    evlo::Float64 # Degrees
    evla::Float64
    evdp::Float64 # km below sea level
    stlo::Float64
    stla::Float64
    stdp::Float64 # km below sea level
    stnm::String
    phi::Float64 # degrees
    dphi::Float64
    dt::Float64 # s
    ddt::Float64
    spol::Float64
    dspol::Float64
    ex::Float64 # m relative to
    ey::Float64
    ez::Float64
    sx::Float64
    sy::Float64
    sz::Float64
    dist::Float64 # Distance in m from source to receiver, straight line
    incidence::Float64 # Degrees from the downward direction
end

"""
    SWS(evlo, evla, evdp, stlo, stla, stdp, stnm, phi, dphi, dt, ddt, spol, dspol, origin) -> ::SWSobs

Return an `SWS` instance with extra fields filled in.
`origin` is a Geodesy.LLA instance and must be given
"""
function SWS(evlo, evla, evdp, stlo, stla, stdp, stnm, phi, dphi, dt, ddt, spol, dspol, origin)
    e = ENU(LLA(evla, evlo, -evdp*1000.), origin, wgs84)
    s = ENU(LLA(stla, stlo, -stdp*1000.), origin, wgs84)
    dist = distance(e, s)
    incidence = incidence_angle(e, s)
    SWS(evlo, evla, evdp, stlo, stla, stdp, stnm, phi, dphi, dt, ddt, spol, dspol,
           e.e, e.n, e.u, s.e, s.n, s.u, dist, incidence)
end

"Allow indexing into arrays of SWS by using the symbol form of the fieldname"
getindex(A::Array{SWS}, s::Symbol) = Array{typeof(getfield(A[1], s))}([getfield(a, s) for a in A])
getindex(t::SWS, s::Symbol) = getfield(t, s)
"Allow indexing into arrays of SWS obs by using a boolean array"
function getindex(A::Array{SWS}, b::BitArray)
    @assert length(b) == length(A) "Indexing array does not have correct dimensions"
    out = Array{SWS,1}()
    for i in 1:length(A)
        if b[i] push!(out, A[i]) end
    end
    out
end

"""
    SWSTime

Type containing a shear wave splitting measurement as well as the date and time of
the event.

Access the usual SWS field with `.sws` and the date with `.time`.
"""
immutable SWSTime <: AbstractSWS
    sws::SWS
    time::DateTime
end

"""
    getindex(s::Union{SWSTime,Array{SWSTime}}, k::Symbol) -> v

Return an array of values `v` for the fields named `k` in an array of `SWSTime`s.

Use `:time` for `k` to get the `DateTime`s for each measurement; otherwise, use `:sws`
to get an array of `SWS`s; or finally use any fieldname in the `SWS` type to get that.
"""
getindex(s::Array{SWSTime}, k::Symbol) = k == :time ? [ss.time for ss in s] :
                                             k == :sws ? [ss.sws for ss in s] :
                                             [getfield(ss.sws, k) for ss in s]
getindex(s::SWSTime, k::Symbol) = k == :time ? s.time :
                                      k == :sws ? s.sws : getfield(s.sws, k)

#=
    Functions on AbstractSWSs
=#

"Return the azimuth from the station to the receiver"
azimuth(s::AbstractSWS) = SphericalGeom.azimuth(s[:evlo], s[:evla], s[:stlo], s[:stla], true)
"Return the backazimuth from the receiver to the station"
backazimuth(s::AbstractSWS) = SphericalGeom.azimuth(s[:stlo], s[:stla], s[:evlo], s[:evla], true)

"Return the surface distance between the event and station"
surf_dist(s::Union{AbstractSWS,Array{<:AbstractSWS}}) = sqrt.((s[:ex] .- s[:sx]).^2 + (s[:ey] .- s[:sy]).^2)
"Return the incidence angle, measured away from downwards, from the event to the station"
incidence_angle(e::ENU, s::ENU) = rad2deg(atan2(s.u - e.u, sqrt((s.e - e.e).^2 + (s.n - e.n).^2)))

"""
    filter!(a::Array{SWS}, args...; kwargs...) -> a
Remove from the array of SWS measurements `a` any measurements which do not match
the filters specified by the `args` and `kwargs`.

    filter(a::Array{SWS}, args...; kwargs...) -> a_copied

Return a copy of `a` without non-matching measurements.

### Arguments

Arguments must be tuples containing a function and two values, the upper and lower limits.
Function must accept one argument (a SWS) and return one value.

```julia
julia> horizontal_dist(s::SWS) = sqrt((s.ex - s.sx)^2 + (s.ey - s.sy)^2);

julia> filter!(splits, (horizontal_dist, 0, 20_000)); # Measurements with epicentral distance less than 20 km
```

### Keyword arguments

Each keyword argument should corresond to
a field in the SWS type, and to it should be passed a tuple of two values, the lower
and upper bound (inclusive) that field may take.

The exception is the String field `stnm`, which accepts only a regular expression.

```julia
julia> filter!(splits, evdp=(-Inf,100)); # Keep only events shallower than 100 km
```
"""
function filter!(a::Array{<:AbstractSWS}, args...; kwargs...)
    # args
    for arg in args
        length(arg) == 3 && typeof(arg[1]) <: Function ||
            throw(ArgumentError("Arguments must be length-3 tuples"))
    end
    # kwargs
    field_types = Dict{Symbol,DataType}()
    for (f,t) in zip(fieldnames(SWS), SWS.types)
        field_types[f] = t
    end
    for (k,v) in kwargs
        k in keys(field_types) && field_types[k] == String &&
            !(typeof(v) <: Union{Regex,String}) &&
            throw(ArgumentError("Value for $k must be a String or Regex"))
    end
    delete_indices = Int[]
    fields_done = false
    for (i,s) in enumerate(a)
        fields_done = false
        for (f,v1,v2) in args
            fields_done && continue
            if !(v1 <= f(s) <= v2)
                push!(delete_indices, i)
                fields_done = true
            end
        end
        for (k,v) in kwargs
            fields_done && continue
            if haskey(field_types, k) && field_types[k] == String
                if !ismatch(Regex(v), s[k])
                    push!(delete_indices, i)
                    fields_done = true
                end
            else
                if !(v[1] <= s[k] <= v[end])
                    push!(delete_indices, i)
                    fields_done = true
                end
            end
        end
    end
    deleteat!(a, delete_indices)
    a
end

filter(a::Array{<:AbstractSWS}, args...; kwargs...) = filter!(deepcopy(a), args...; kwargs...)

@doc (@doc filter!) filter

"""
    midpoint(s::Union{SWS,Array{SWS}}) -> x, y, z

Return the midpoint of a straight line between the source and receiver in cartesian coordinates.
"""
midpoint(s::Union{AbstractSWS,Array{<:AbstractSWS}}) =
    (s[:sx] + s[:ex])/2, (s[:sy] + s[:ey])/2, (s[:sz] + s[:ez])/2

"""
    midpoint_lonlat(s::Union{SWS,Array{SWS}})

Return the midpoint of a straight line between the source and receiver
in (longitude, latitude, depth) coordinates.  This assumes that distances are small
and longitude and latitude can be treated as cartesian.
"""
midpoint_lonlat(s::Union{SWS,Array{SWS}}) =
    (s[:evlo] + s[:stlo])/2, (s[:evla] + s[:evla])/2 + (s[:evdp] - s[:stel]/1000)/2

# TODO: Profile this and speed it up: it's probably really slow
"""
    read_sheba(file, origin, stdps=Dict{String,<:Real}) -> s::Vector{SWS}

Read the shear wave splitting measurements `s` from a `file` written by the SHEBA [1]
program.  `origin` is the `Gedoesy.LLA` coordinate of the centre of the UR grid.

To set the depth of stations, supply `stdps` (a `Dict{String,<:Real}`) which
contains the station names as keys and their depths in m as the values.

##### References

1. Wüstefeld, A., Al-Harrasi, O., Verdon, J., Wookey, J. and Kendall, J-M. (2010).
   A strategy for automated analysis of passive microseismic data to image seismic
   anisotropy and fracture characteristics.  Geophysical Prospecting, 58, 755–773.
   doi:10.1111/j.1365-2478.2010.00891.x (available at https://github.com/jwookey/sheba)
"""
function read_sheba(file, origin, stdps=Dict{String,Float64}())
    d = readdlm(file)
    d = d[.!ismatch.(r"^[#%]", string.(d[:,1])), :]
    n = size(d, 1)
    s = Array{SWS}(n)
    for i in 1:n
        elon = d[i,4]
        elat = d[i,3]
        slon = d[i,6]
        slat = d[i,5]
        edep = d[i,7]
        sta = d[i,20]
        phi = d[i,11]
        dphi = d[i,12]
        dt = d[i,13]
        ddt = d[i,14]
        spol = d[i,15]
        dspol = d[i,16]
        stdp = haskey(stdps, sta) ? stdps[sta] : 0.0
        s[i] = SWS(elon, elat, edep, slon, slat, stdp, sta, phi, dphi, dt, ddt,
            spol, dspol, origin)
    end
    s
end

"""
    read_sheba_time(file, args...; comment=r\"^[#%]\") -> st::Vector{SWSTime}

Read results from a SHEBA-formatted file, including the date and time of the
events listed therein, returning the values `st`.

Optionally specify the Regex used to determine which lines are comments with `comment`.
"""
function read_sheba_time(file, args...; comment=r"^[#%]")
    s = read_sheba(file, args...)
    n = length(s)
    st = Vector{SWSTime}(n)
    data = readdlm(file)
    data = data[.!ismatch.(comment, string.(data[:,1])), :]
    size(data, 1) == n || error("Length of splits read ($n) does "*
        "not match the number of non-comment lines ($(size(data, 1))) in file '$file'")
    for i in 1:n
        y = data[i,1]÷1000
        d = data[i,1]%1000
        h = data[i,2]÷100
        m = data[i,2]%100
        date = Date(y) + Day(d)
        st[i] = SWSTime(s[i], DateTime(year(date), month(date), day(date), h, m))
    end
    st
end

"""
    write_sheba(s::Array{<:AbstractSWS}, file)

Write the shear wave splitting measurements `s` to a file `file` in the format used
by the SHEBA [1] program.

##### References

1. Wüstefeld, A., Al-Harrasi, O., Verdon, J., Wookey, J. and Kendall, J-M. (2010).
   A strategy for automated analysis of passive microseismic data to image seismic
   anisotropy and fracture characteristics.  Geophysical Prospecting, 58, 755–773.
   doi:10.1111/j.1365-2478.2010.00891.x (available at https://github.com/jwookey/sheba)
"""
function write_sheba(s::Array{<:AbstractSWS}, file)
    open(file, "w") do f
        println(f, "%  DATE TIME    EVLA    EVLO    STLA    STLO    EVDP    DIST     AZI     BAZ    FAST   DFAST    TLAG   DTLAG    SPOL   DSPOL    WBEG    WEND  STAT %  DATE TIME EIGORIG EIGCORR      Q     SNR   NDF   STAT SNR")
        for ss in s
            datetime = if typeof(ss) <: SWSTime
                t = ss[:time]
                @sprintf("%04d%03d %02d%02d", year(t), _dayofyear(t), hour(t), minute(t))
            else
                "DATE TIME"
            end
            println(f, datetime, " ", ss[:evla], " ", ss[:evlo], " ", ss[:stla], " ",
                ss[:stlo], " ", ss[:evdp], " ",
                SphericalGeom.delta(ss[:evlo], ss[:evla], ss[:stlo], ss[:stla], true), " ",
                azimuth(ss), " ", backazimuth(ss), " ",
                ss[:phi], " ", ss[:dphi], " ", ss[:dt], " ", ss[:ddt], " ", ss[:spol], " ",
                ss[:dspol], " WBEG WEND % ", ss[:stnm], " DATE TIME EIGORIG EIGCORR Q SNR NDF ",
                ss[:stnm], " SNR")
        end
    end
end

"""
    _dayofyear(d::DateTime) -> doy::Int

Return the day of the year of a `DateTime` `d`."""
_dayofyear(d::DateTime) = Dates.value(Date(year(d), month(d), day(d))
                                       - Date(year(d)) + Day(1))

end # module

max_stars_repo_name: anowacki/assorted-julia-modules

max_issues_repo_name: anowacki/assorted-julia-modules

max_forks_repo_name: anowacki/assorted-julia-modules


text:
####
# Matrix memory layout traits
#
# if MemoryLayout(A) returns BandedColumnMajor, you must override
# pointer and leadingdimension
# in addition to the banded matrix interface
####

abstract type AbstractBlockBandedLayout <: MemoryLayout end

struct BandedBlockBandedColumnMajor <: AbstractBlockBandedLayout end
struct BandedBlockBandedRowMajor <: AbstractBlockBandedLayout end
struct BlockBandedColumnMajor <: AbstractBlockBandedLayout end
struct BlockBandedRowMajor <: AbstractBlockBandedLayout end

transposelayout(::BandedBlockBandedColumnMajor) = BandedBlockBandedRowMajor()
transposelayout(::BandedBlockBandedRowMajor) = BandedBlockBandedColumnMajor()
transposelayout(::BlockBandedColumnMajor) = BlockBandedRowMajor()
transposelayout(::BlockBandedRowMajor) = BlockBandedColumnMajor()
conjlayout(::Type{<:Complex}, M::AbstractBlockBandedLayout) = ConjLayout(M)

# Here we override broadcasting for banded matrices.
# The design is to to exploit the broadcast machinery so that
# banded matrices that conform to the banded matrix interface but are not
# <: AbstractBandedMatrix can get access to fast copyto!, lmul!, rmul!, axpy!, etc.
# using broadcast variants (B .= A, B .= 2.0 .* A, etc.)

abstract type AbstractBlockBandedStyle <: AbstractArrayStyle{2} end
struct BandedBlockBandedStyle <: AbstractBlockBandedStyle end
abstract type AbstractBlockSkylineStyle <: AbstractBlockBandedStyle end
struct BlockSkylineStyle <: AbstractBlockSkylineStyle end
struct BlockBandedStyle <: AbstractBlockSkylineStyle end
BlockSkylineStyle(::Val{2}) = BlockBandedStyle()
BlockBandedStyle(::Val{2}) = BlockBandedStyle()
BandedBlockBandedStyle(::Val{2}) = BandedBlockBandedStyle()
BroadcastStyle(::DefaultArrayStyle{2}, ::AbstractBlockBandedStyle) = DefaultArrayStyle{2}()
BroadcastStyle(::AbstractBlockBandedStyle, ::DefaultArrayStyle{2}) = DefaultArrayStyle{2}()

BroadcastStyle(::BlockSkylineStyle, ::BandedBlockBandedStyle) = BlockSkylineStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::BlockSkylineStyle) = BlockSkylineStyle()
BroadcastStyle(::BlockSkylineStyle, ::BlockBandedStyle) = BlockSkylineStyle()
BroadcastStyle(::BlockBandedStyle, ::BlockSkylineStyle) = BlockSkylineStyle()

BroadcastStyle(::BlockBandedStyle, ::BandedBlockBandedStyle) = BlockBandedStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::BlockBandedStyle) = BlockBandedStyle()

####
# Default to standard Array broadcast
#
# This is because, for example, exp.(B) is not banded
####

copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle}) =
   copyto!(dest, Broadcasted{DefaultArrayStyle{2}}(bc.f, bc.args, bc.axes))

similar(bc::Broadcasted{<:AbstractBlockBandedStyle}, ::Type{T}) where T =
    similar(Broadcasted{DefaultArrayStyle{2}}(bc.f, bc.args, bc.axes), T)

##
# copyto!
##

function checkblocks(A, B)
    Arows, Acols = cumulsizes(A)
    Brows, Bcols = cumulsizes(B)
    if Acols ≠ Bcols || Arows ≠ Brows
        throw(DimensionMismatch("*"))
    end
end

function _blockbanded_copyto!(dest::AbstractMatrix{T}, src::AbstractMatrix) where T
    @boundscheck checkblocks(dest, src)

    dl, du = colblockbandwidths(dest)
    sl, su = colblockbandwidths(src)
    M,N = nblocks(src)
    # Source matrix must fit within bands of destination matrix
    all(dl .≥ min.(sl,Ref(M-1))) && all(du .≥ min.(su,Ref(N-1))) || throw(BandError(dest))

    for J = 1:N
        for K = max(1,J-du[J]):min(J-su[J]-1,M)
            view(dest,Block(K),Block(J)) .= zero(T)
        end
        for K = max(1,J-su[J]):min(J+sl[J],M)
            view(dest,Block(K),Block(J)) .= view(src,Block(K),Block(J))
        end
        for K = max(1,J+sl[J]+1):min(J+dl[J],M)
            view(dest,Block(K),Block(J)) .= zero(T)
        end
    end
    dest
end

function blockbanded_copyto!(dest::AbstractMatrix, src::AbstractMatrix)
    if isblockbanded(dest)
        _blockbanded_copyto!(dest, src)
    else
        _blockbanded_copyto!(PseudoBlockArray(dest, blocksizes(src).block_sizes), src)
    end
    dest
end

copyto!(dest::AbstractMatrix, src::AbstractBlockBandedMatrix) = blockbanded_copyto!(dest, src)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(identity)})
    (A,) = bc.args
    dest ≡ A && return dest
    blockbanded_copyto!(dest, A)
end

##
# lmul!/rmul!
##

######
# extra marrix routines
#####

function blockbanded_fill!(B::AbstractMatrix{T}, x) where T
    x == zero(T) || throw(BandError(B))

    M,N = nblocks(B)
    for J = 1:N, K = blockcolrange(B,J)
        fill!(view(B,K,Block(J)), x)
    end
    B
end

function blockbanded_rmul!(B::AbstractMatrix{T}, x::Number) where T
    x == zero(T) || throw(BandError(B))

    M,N = nblocks(B)
    for J = 1:N, K = blockcolrange(B,J)
        rmul!(view(B,K,Block(J)), x)
    end
    B
end

function blockbanded_lmul!(x::Number, B::AbstractMatrix{T}) where T
    x == zero(T) || throw(BandError(B))

    M,N = nblocks(B)
    for J = 1:N, K = blockcolrange(B,J)
        lmul!(x, view(B,K,Block(J)))
    end
    B
end

function blockbanded_fill!(B::AbstractBlockBandedMatrix, x::Number)
    x == zero(T) || throw(BandError(B))

    fill!(B.data, x)
end

function blockbanded_rmul!(A::AbstractBlockBandedMatrix, x::Number)
    rmul!(A.data, x)
    A
end

function blockbanded_lmul!(x::Number, A::AbstractBlockBandedMatrix)
    lmul!(x, A.data)
    A
end

lmul!(x::Number, A::AbstractBlockBandedMatrix) = blockbanded_lmul!(x, A)
rmul!(A::AbstractBlockBandedMatrix, x::Number) = blockbanded_rmul!(A, x)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(-), <:Tuple{<:AbstractMatrix}})
    A, = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_lmul!(-1, dest)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(-), <:Tuple{<:AbstractMatrix}}, ::Type{T}) where T =
    similar(bc.args[1], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:Number,<:AbstractMatrix}})
    α,A = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_lmul!(α, dest)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:Number,<:AbstractMatrix}}, ::Type{T}) where T =
    similar(bc.args[2], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:AbstractMatrix,<:Number}})
    A,α = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_rmul!(dest, α)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:AbstractMatrix,<:Number}}, ::Type{T}) where T =
    similar(bc.args[1], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(/), <:Tuple{<:AbstractMatrix,<:Number}})
    A,α = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_rmul!(dest, inv(α))
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(/), <:Tuple{<:AbstractMatrix,<:Number}}, ::Type{T}) where T =
    similar(bc.args[1], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(\), <:Tuple{<:Number,<:AbstractMatrix}})
    α,A = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_lmul!(inv(α), dest)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(\), <:Tuple{<:Number,<:AbstractMatrix}}, ::Type{T}) where T =
    similar(bc.args[2], T)

##
# axpy!
##

# these are the routines of the banded interface of other AbstractMatrices
function blockbanded_axpy!(a, X::AbstractMatrix, Y::AbstractMatrix)
    size(X) == size(Y) || throw(DimensionMismatch())

    for J=Block(1):Block(nblocks(X,2)), K=blockcolrange(X,J)
        view(Y,K,J) .= a .* view(X,K,J) .+ view(Y,K,J)
    end
    Y
end

function _combine_blocksizes(A, B)
    A_size = blocksizes(A).block_sizes
    A_size == blocksizes(B).block_sizes || throw(DimensionMismatch())
    A_size
end

_combine_blocksizes(::Diagonal, B) = blocksizes(B)
_combine_blocksizes(A, ::Diagonal) = blocksizes(A)

function _combined_blocksizes(A, B)
    blocksizes(A) == blocksizes(B) || throw(DimensionMismatch("Block sizes do not agree"))
    (A,B)
end

_combined_blocksizes(A::Diagonal, B) = PseudoBlockArray(A, blocksizes(B).block_sizes), B
_combined_blocksizes(A, B::Diagonal) = A, PseudoBlockArray(B, blocksizes(A).block_sizes)

for op in (:+, :-)
    @eval begin
        function similar(bc::Broadcasted{<:BlockSkylineStyle, <:Any, typeof($op), <:Tuple{<:AbstractMatrix,<:AbstractMatrix}}, ::Type{T}) where T
            A,B = bc.args
            Al,Au = colblockbandwidths(A)
            Bl,Bu = colblockbandwidths(B)
            BlockSkylineMatrix{T}(undef, _combine_blocksizes(A,B), (max.(Al,Bl), max.(Au,Bu)))
        end

        function similar(bc::Broadcasted{<:BlockBandedStyle, <:Any, typeof($op), <:Tuple{<:AbstractMatrix,<:AbstractMatrix}}, ::Type{T}) where T
            A,B = bc.args
            Al,Au = blockbandwidths(A)
            Bl,Bu = blockbandwidths(B)
            BlockBandedMatrix{T}(undef, _combine_blocksizes(A,B), (max(Al,Bl), max(Au,Bu)))
        end

        function similar(bc::Broadcasted{<:BandedBlockBandedStyle, <:Any, typeof($op), <:Tuple{<:AbstractMatrix,<:AbstractMatrix}}, ::Type{T}) where T
            A,B = bc.args
            Al,Au = blockbandwidths(A)
            Bl,Bu = blockbandwidths(B)
            Aλ,Aμ = subblockbandwidths(A)
            Bλ,Bμ = subblockbandwidths(B)

            BandedBlockBandedMatrix{T}(undef, _combine_blocksizes(A,B), (max(Al,Bl), max(Au,Bu)), (max(Aλ,Bλ), max(Aμ,Bμ)))
        end
    end
end

for op in (:+, :-)
    @eval function copyto!(C::AbstractArray{T}, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof($op),
                                                                <:Tuple{<:AbstractMatrix,<:AbstractMatrix}}) where T
        A,B = _combined_blocksizes(bc.args...)
        A_l,A_u = blockbandwidths(A)
        B_l,B_u = blockbandwidths(B)
        C_l,C_u = blockbandwidths(C)

        size(A) == size(B) == size(C) || throw(DimensionMismatch())
        N,M = nblocks(C)

        for J̃ = 1:M
            J = Block(J̃)
            for K = Block.(max(1,J̃-min(A_u,B_u)):min(N,J̃+min(A_l,B_l)))
                view(C,K,J) .= $op.(view(A,K,J), view(B,K,J))
            end
            for K = Block.(max(1,J̃-A_u):min(N,J̃-B_u-1))
                view(C,K,J) .= view(A,K,J)
            end
            for K = Block.(max(1,J̃-B_u):min(N,J̃-A_u-1))
                view(C,K,J) .= $op.(view(B,K,J))
            end
            for K = Block.(max(1,J̃+B_l+1):min(N,J̃+A_u))
                view(C,K,J) .= view(A,K,J)
            end
            for K = Block.(max(1,J̃+A_l+1):min(N,J̃+B_u))
                view(C,K,J) .= $op.(view(B,K,J))
            end
            for K = Block.(max(J̃-C_u,1):min(J̃-max(A_u,B_u)-1,N))
                view(C,K,J) .= zero(T)
            end
            for K = Block.(max(J̃+max(A_l,B_l)+1,1):min(J̃+C_u,N))
                view(C,K,J) .= zero(T)
            end
        end

        C
    end
end

function copyto!(dest::AbstractArray{T}, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(+),
                                                        <:Tuple{<:Broadcasted{<:AbstractBlockBandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                                                        <:AbstractMatrix}}) where T
    αA,B = bc.args
    α,A = αA.args
    dest ≡ B || blockbanded_copyto!(dest, B)
    blockbanded_axpy!(α, A, dest)
end

function similar(bc::Broadcasted{BlockBandedStyle, <:Any, typeof(+),
                        <:Tuple{<:Broadcasted{<:AbstractBlockBandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                        <:AbstractMatrix}}, ::Type{T}) where T
    αA,B = bc.args
    α,A = αA.args
    Al,Au = blockbandwidths(A)
    Bl,Bu = blockbandwidths(B)
    BlockBandedMatrix{T}(undef, _combine_blocksizes(A,B), (max(Al,Bl), max(Au,Bu)))
end

function similar(bc::Broadcasted{BandedBlockBandedStyle, <:Any, typeof(+),
                        <:Tuple{<:Broadcasted{<:AbstractBlockBandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                        <:AbstractMatrix}}, ::Type{T}) where T
    αA,B = bc.args
    α,A = αA.args
    Al,Au = blockbandwidths(A)
    Bl,Bu = blockbandwidths(B)
    Aλ,Aμ = subblockbandwidths(A)
    Bλ,Bμ = subblockbandwidths(B)

    BandedBlockBandedMatrix{T}(undef, _combine_blocksizes(A,B), (max(Al,Bl), max(Au,Bu)), (max(Aλ,Bλ), max(Aμ,Bμ)))
end

####
# Special case for Diagonal multiplicartion
####


text:
### State estimation tests
path = abspath(joinpath(dirname(Base.find_package("JuliaGrid")), ".."), "test/estimation/")
case14dc = string(path, "case14dc.xlsx")
estimation_incdc = string(path, "estimation_incdc.xlsx")
manousakis2010 = string(path, "manousakis2010.xlsx")
case14seobsbad = string(path, "case14seobsbad.xlsx")

@testset "DC State Estimation" begin
    runse("case14se.xlsx", "dc"; save = string(path, "save.h5"))
    runse("case14se.xlsx", "dc"; save = string(path, "save.xlsx"))
        rm(string(path, "save.h5"))
        rm(string(path, "save.xlsx"))
        @test !("save.h5" in cd(readdir, path))
        @test !("save.xlsx" in cd(readdir, path))

    runse(estimation_incdc, "dc"; save = string(path, "save.h5"))
    runse(estimation_incdc, "dc"; save = string(path, "save.xlsx"))
            rm(string(path, "save.h5"))
            rm(string(path, "save.xlsx"))
            @test !("save.h5" in cd(readdir, path))
            @test !("save.xlsx" in cd(readdir, path))

    runse(estimation_incdc, "dc")
        fieldsres = (:main, :flow, :estimate, :error, :baddata, :observability)
        @test all(fieldsres == fieldnames(JuliaGrid.StateEstimationDC))

    pf, = runpf(case14dc, "dc")
    se, = runse(case14dc, "dc", "main", "flow", "estimate", "error")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - se.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - se.main[:, 3])) < accuracy
        @test maximum(abs.(pf.flow[:, 4] - se.flow[:, 4])) < accuracy
        @test maximum(abs.(pf.flow[:, 5] - se.flow[:, 5])) < accuracy
        @test se.error[1] < 1e-10
        @test se.error[2] < 1e-10
        @test se.error[3] < 1e-10

    lav, = runse(case14dc, "dc", "lav", "main")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - lav.main[:, 2])) < accuracy

    data = runmg(estimation_incdc; runflow = 0, pmuset = "complete", pmuvariance = ["complete" 1e-30], legacyset = ["Pij" 0])
        results, measurements, = runse(data, "dc", "estimate"; bad = ["pass" 3])
        @test maximum(abs.(results.estimate[:, 8] - results.estimate[:, 10])) < 1e-8
        @test all(measurements.pmuVoltage[:, 7] .== 1)
        @test all(measurements.pmuVoltage[:, 6] .== 1e-30)
        @test all(measurements.pmuVoltage[:, 7] .== 1)
        @test all(measurements.pmuCurrent .== 0)
        @test all(measurements.legacyCurrent .== 0)
        @test all(measurements.legacyInjection .== 0)
        @test all(measurements.legacyVoltage .== 0)

    results, = runse(manousakis2010, "dc", "observe", "estimate")
        islands = [[5; 6; 11; 12; 13], [10], [14], [4; 7; 8; 9], [1], [2; 3]]
        for (k, i) in enumerate(islands)
            for j in results.observability
                if issubset(i, j)
                    islands[k] = []
                end
            end
        end
        idx = findall(x->x==3, results.estimate[:, 2])
        @test (all(isempty.(islands)))
        @test all(results.estimate[idx, 4] .== 3) && all(results.estimate[idx, 5] .== [9; 10])

    results, = runse(manousakis2010, "dc", "observe"; observe = ["pivot" 1e-8 "Pij" 1e4])
        idx = findall(x->x==3, results.estimate[:, 2])
        @test all(results.estimate[idx, 4] .== 1) && all(results.estimate[idx, 5] .== [16; 17])

    results, = runse(manousakis2010, "dc", "observe"; observe = ["pivot" 1e-8 "Ti" 1e4])
        idx = findall(x->x==3, results.estimate[:, 2])
        @test all(results.estimate[idx, 4] .== 8) && all(results.estimate[idx, 5] .== [10; 14])

    results, = runse(manousakis2010, "dc"; observe = ["islandBP" "islandMax" 500 "islandBreak" 15 "islandStopping" 1 "islandTreshold" 1e10 "restoreBP" "restoreMax" 500])
        islands = [[5; 6; 11; 12; 13], [10], [14], [4; 7; 8; 9], [1], [2; 3]]
        for (k, i) in enumerate(islands)
            for j in results.observability
                if issubset(i, j)
                islands[k] = []
                end
            end
        end
        idx = findall(x->x==3, results.estimate[:, 2])
        @test (all(isempty.(islands)))
        @test all(results.estimate[idx, 4] .== 3) && all(results.estimate[idx, 5] .== [9; 10])

    results, = runse("case14se.xlsx", "dc", "observe", "bad", "main"; observe = ["pivot" 1e-8 "Ti" 1e-4 "Pi" 1e-4], bad = ["pass" 4])
        Ti = results.main[:, 2]
        results, = runse(case14seobsbad, "dc", "observe", "bad", "main"; observe = ["pivot" 1e-8 "Ti" 1e-4 "Pi" 1e-4], bad = ["pass" 2])
        Tinew = results.main[:, 2]
        @test maximum(abs.(Ti - Tinew)) < 1.0

    data = runmg("case_ACTIVSg70k.h5"; legacyset = "complete", legacyvariance = ["Pij" 1e-2 "Pi" 1e-3])
        measure, = data
        measure.legacyInjection[105, 2] = 3000; measure.legacyInjection[105, 4] = 1
        measure.legacyInjection[5, 2] = 1000; measure.legacyInjection[5, 4] = 1
        results, = runse(data, "dc"; bad = ["pass" 2 "treshold" 2])
        @test results.baddata[1, 3] == 3 && results.baddata[1, 4] == 105
        @test results.baddata[2, 3] == 3 && results.baddata[2, 4] == 5
end

@testset "PMU State Estimation" begin
    runse("case14se.xlsx", "pmu"; save = string(path, "save.h5"))
    runse("case14se.xlsx", "pmu"; save = string(path, "save.xlsx"))
        rm(string(path, "save.h5"))
        rm(string(path, "save.xlsx"))
        @test !("save.h5" in cd(readdir, path))
        @test !("save.xlsx" in cd(readdir, path))

    runse("case14se.xlsx", "dc")
        fieldsres = (:main, :flow, :estimate, :error, :baddata, :observability)
        @test all(fieldsres == fieldnames(JuliaGrid.StateEstimationPMU))

    pf, = runpf("case30.h5", "nr")
    data = runmg("case30.h5"; runflow = 1, pmuset = "complete", pmuvariance = ["complete" 1e-120])
    se, = runse(data, "pmu", "main", "flow", "estimate", "error")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - se.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - se.main[:, 3])) < accuracy
        @test se.error[1] < 1e-10
        @test se.error[2] < 1e-10

    lav, = runse(data, "pmu", "lav")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - lav.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - lav.main[:, 3])) < accuracy
end

@testset "Nonlinear State Estimation" begin
    pf, = runpf("case14.h5", "nr")
    data = runmg("case14.h5"; runflow = 1, pmuset = "complete", pmuvariance = ["complete" 1e-30], legacyset = "complete", legacyvariance = ["complete" 1e-30])
    se, = runse(data, "nonlinear", "main", "flow", "estimate", "error")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - se.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - se.main[:, 3])) < accuracy

    pf, = runpf("case300.h5", "nr")
    data = runmg("case300.h5"; runflow = 1, legacyset = ["Pij" "all" "Pi" "all" "Qij" "all" "Qi" "all"], legacyvariance = ["complete" 1e-120])
    se, = runse(data, "nonlinear")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - se.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - se.main[:, 3])) < accuracy

    data = runmg("case300.h5"; runflow = 1, pmuset = ["Ti" "all" "Vi" "all" "Iij" "all"], pmuvariance = ["complete" 1e-30])
    se, = runse(data, "nonlinear")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - se.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - se.main[:, 3])) < accuracy

    pf, = runpf("case_ieee30.h5", "nr")
    data = runmg("case_ieee30.h5"; runflow = 1, pmuset = ["Ti" "all" "Vi" "all" "Dij" 10], pmuvariance = ["complete" 1e-120])
    se, = runse(data, "nonlinear")
        accuracy = 1e-8
        @test maximum(abs.(pf.main[:, 2] - se.main[:, 2])) < accuracy
        @test maximum(abs.(pf.main[:, 3] - se.main[:, 3])) < accuracy
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
module StencilInterpolations

using CUDA
using VectorizationBase: Vec, vsum, vadd

include("utils.jl")
include("gather.jl")
include("scatter.jl")
include("kernels.jl")
include("bilinear.jl")
include("trilinear.jl")

export scattering, scattering!, gathering!
export grid2particle, grid2particle!
export lerp, ndlinear, random_particles

end # module


text:
# utils_TaylorN.jl: N-variables Taylor expansions
#
# Last modification: 2014.06.07
#
# Luis Benet & David P. Sanders
# UNAM
#

## Default values for the maximum degree of polynomials (MAXORDER) and
##  number of variables considered (NUMVARS)
const MAXORDER = [4]
const NUMVARS = [2]

## Hash tables
"""
`generateIndicesTable`: generates the dictionary `indicesTable`.
NUMVARS: number of variables for the polynomial expansions <--> lNiv
MAXORDER: maximum degree of polynomials <--> nPart
`generatePosTable`: generates the dictionary `posTable`; is the "inverse" of `indicesTable`
"""
function generateIndicesTable()
    numVars::Int = NUMVARS[end]
    maxDeg::Int = MAXORDER[end]
    DDic = Dict{Int, Array{Int,1}}()
    #
    if numVars==1
        info( string("Using `Taylor` constructor rather than `TaylorN` is\n",
            "MUCH faster for 1-variable expansions.\n") )
        for k = 0:maxDeg
            DDic[k+1] = [k]
        end
        return DDic
    end
    #
    nCoefTot = binomial( numVars+maxDeg, maxDeg)
    iindices = zeros(Int, numVars)
    pos = 0
    for kDeg = 0:maxDeg
        iV = numVars
        for iz = 0:kDeg
            iindices[end] = iz
            @inbounds pos, DDic[pos] = pos2indices!( iV, kDeg, iindices, pos, DDic)
        end
    end
    return DDic
end
function pos2indices!(iV::Int, kDeg::Int, iIndices::Array{Int,1}, pos::Int, dict::Dict{Int,Array{Int,1}})
    #
    jVar = iV-1
    kDegNext = kDeg - iIndices[iV]
    if jVar > 1
        for jDeg=0:kDegNext
            iIndices[jVar] = jDeg
            @inbounds pos, dict[pos] = pos2indices!( jVar, kDegNext, iIndices, pos, dict)
        end
    else
        iIndices[1] = kDegNext
        pos += 1
        @inbounds dict[pos] = iIndices[1:end]
    end
    return pos, iIndices[1:end]
end

const indicesTable = [ generateIndicesTable() ]
const sizeTable = [ length( indicesTable[end]) ]

function generatePosTable()
    DDic = Dict{Array{Int,1}, Int}()
    for i = 1:sizeTable[end]
        v = indicesTable[end][i]
        DDic[ v ] = i
    end
    return DDic
end
const posTable = [ generatePosTable() ]

## Functions to obtain the number of homogenous coefficients of given degree
# """Returns the number of homogeneous coefficients of degree k for NUMVARS:
#     binomial(k+NUMVARS-1,k)"""
# function numHomogCoefK(k::Int)
#     k == 0 && return 1
#     return binomial( k + NUMVARS[end] - 1, k )
# end
"""Returns the position of the first homogeneous coefficient of degree k for NUMVARS:
    binomial(k+NUMVARS-1,k-1)+1"""
function posHomogCoefK(k::Int)
    k == 0 && return 1
    return binomial( k + NUMVARS[end]-1, k-1 ) + 1
end
"""Generates the auxiliary table of the initial and final position for the coefficients of
homogeneous polynomial of given degree (`posAuxTable`)"""
function generateAuxTable()
    DDic = Dict{Int,(Int, Int)}()
    for k=0:MAXORDER[end]
        posI = posHomogCoefK(k)
        posF = posHomogCoefK(k+1)-1
        DDic[k] = ( posI, posF )
    end
    return DDic
end
const posAuxTable = [ generateAuxTable() ]

## Utilities to get/set MAXORDER and NUMVARS; they reset the indicesTable
get_maxOrder() = MAXORDER[end]
function set_maxOrder(n::Int)
    @assert n >= 0
    MAXORDER[end] = n
    info(string("MAXORDER is now ", n, "; hash tables regenerated.\n"))
    indicesTable[end] = generateIndicesTable()
    sizeTable[end] = length( indicesTable[end] )
    posTable[end] = generatePosTable()
    posAuxTable[end] = generateAuxTable()
    return n
end
#
get_numVars() = NUMVARS[end]
function set_numVars(n::Int)
    @assert n > 0
    NUMVARS[end] = n
    info(string("NUMVARS is now ", n, "; hash tables regenerated.\n"))
    indicesTable[end] = generateIndicesTable()
    sizeTable[end] = length( indicesTable[end] )
    posTable[end] = generatePosTable()
    posAuxTable[end] = generateAuxTable()
    return n
end

## TaylorN Constructors ##
immutable TaylorN{T<:Number} <: AbstractSeries{T,NUMVARS[end]}
   coeffs :: Array{T,1}
   order :: Int
   numVars :: Int
   function TaylorN(coeffs::Array{T,1}, order::Int, numVars::Int)
        @assert order <= MAXORDER[end]
        lencoef = length(coeffs)
        nCoefTot = sizeTable[end]
        @assert lencoef <= nCoefTot
        v = zeros(T, nCoefTot)
        #@inbounds v[1:lencoef] = coeffs[1:lencoef]
        for i = 1:lencoef
            @inbounds v[i] = coeffs[i]
        end
        new(v, order, NUMVARS[end])
   end
end
TaylorN{T<:Number}(x::TaylorN{T}, order::Int) = TaylorN{T}(x.coeffs, order, NUMVARS[end])
TaylorN{T<:Number}(x::TaylorN{T}) = TaylorN{T}(x.coeffs, x.order, NUMVARS[end])
#
TaylorN{T<:Number}(coeffs::Array{T,1}, order::Int) = TaylorN{T}(coeffs, order, NUMVARS[end])
TaylorN{T<:Number}(coeffs::Array{T,1}) = TaylorN{T}(coeffs, MAXORDER[end], NUMVARS[end])
#
TaylorN{T<:Number}(x::T, order::Int) = TaylorN{T}([x], order, NUMVARS[end])
TaylorN{T<:Number}(x::T) = TaylorN{T}([x], 0, NUMVARS[end])

## Type, length ##
eltype{T<:Number}(::TaylorN{T}) = T
length(a::TaylorN) = length( a.coeffs )
get_numVars(x::TaylorN) = x.numVars
get_maxOrder(x::TaylorN) = x.order

## Conversion and promotion rules ##
convert{T<:Number}(::Type{TaylorN{T}}, a::TaylorN) = TaylorN(convert(Array{T,1}, a.coeffs), a.order)
convert{T<:Number, S<:Number}(::Type{TaylorN{T}}, b::Array{S,1}) = TaylorN(convert(Array{T,1},b))
convert{T<:Number}(::Type{TaylorN{T}}, b::Number) = TaylorN([convert(T,b)], 0)
promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{TaylorN{S}}) = TaylorN{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{Array{S,1}}) = TaylorN{promote_type(T, S)}
promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{S}) = TaylorN{promote_type(T, S)}

## Auxiliary functions ##
# function firstnonzero{T<:Number}(a::Taylor{T})
#     order = a.order
#     nonzero::Int = order+1
#     z = zero(T)
#     for i = 1:order+1
#         if a.coeffs[i] != z
#             nonzero = i-1
#             break
#         end
#     end
#     nonzero
# end
function fixshape{T<:Number, S<:Number}(a::TaylorN{T}, b::TaylorN{S})
    @assert a.numVars == b.numVars
    order = max(a.order, b.order)
    a1, b1 = promote(a, b)
    return TaylorN(a1, order), TaylorN(b1, order), order
end

## real, imag, conj and ctranspose ##
for f in (:real, :imag, :conj)
    @eval ($f){T<:Number}(a::TaylorN{T}) = TaylorN(($f)(a.coeffs), a.order)
end
ctranspose{T<:Number}(a::TaylorN{T}) = conj(a)

## zero and one ##
zero{T<:Number}(a::TaylorN{T}) = TaylorN(zero(T), a.order)
one{T<:Number}(a::TaylorN{T}) = TaylorN(one(T), a.order)

## Equality ##
function ==(a::TaylorN, b::TaylorN)
    a1, b1, order = fixshape(a, b)
    return a1.coeffs == b1.coeffs
end

## Addition and substraction ##
for f in (:+, :-)
    @eval begin
        function ($f)(a::TaylorN, b::TaylorN)
            a1, b1, order = fixshape(a, b)
            v = ($f)(a1.coeffs, b1.coeffs)
            return TaylorN(v, order)
        end
        ($f)(a::TaylorN) = TaylorN(($f)(a.coeffs), a.order)
    end
end

## Multiplication ##
function *(a::TaylorN, b::TaylorN)
    a1, b1, order = fixshape(a, b)
    T = eltype(a1)
    nCoefTot = sizeTable[end]
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = a1.coeffs[1] * b1.coeffs[1]   ## 0th-order
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = mulHomogCoefN( k, a1.coeffs, b1.coeffs )
    end
    TaylorN(coeffs, order)
end
function mulHomogCoefN{T<:Number}( k::Int, ac::Array{T,1}, bc::Array{T,1} )
    k==0 && return ac[1] * bc[1]
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    coeffs = zeros(T, numCoefk)
    for ka = 0:k
        kb = k-ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds coeffs[pos] += ac[ja] * bc[jb]
            end
        end
    end
    return coeffs
end

## Division ##
function /(a::TaylorN, b::TaylorN)
    a1, b1, order = fixshape(a, b)
    #!?orddivfact, cdivfact = divfactorization(a1, b1) # order and coefficient of first factorized term
    @assert b1.coeffs[1] != zero(b1.coeffs[1])
    cLHopital = a1.coeffs[1] / b1.coeffs[1]
    T = typeof(cLHopital)
    v1 = convert(Array{T,1}, a1.coeffs)
    v2 = convert(Array{T,1}, b1.coeffs)
    nCoefTot = sizeTable[end]
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = cLHopital
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = divHomogCoefN( k, v1, v2, coeffs, 0 )
    end
    TaylorN(coeffs, order)
end
# function divfactorization(a1::Taylor, b1::Taylor)
#     # order of first factorized term; a1 and b1 are assumed to be of the same order (length)
#     a1nz = firstnonzero(a1)
#     b1nz = firstnonzero(b1)
#     orddivfact = min(a1nz, b1nz)
#     if orddivfact > a1.order
#         orddivfact = a1.order
#     end
#     cdivfact = a1.coeffs[orddivfact+1] / b1.coeffs[orddivfact+1]
#     aux = abs2(cdivfact)
#     # Is the polynomial factorizable?
#     if isinf(aux) || isnan(aux)
#         info("Order k=$(orddivfact) => coeff[$(orddivfact+1)]=$(cdivfact)")
#         error("Division does not define a Taylor polynomial\n",
#             " or its first non-zero coefficient is Inf/NaN.\n")
#     ##else orddivfact>0
#     ##    warn("Factorizing the polynomial.\n",
#     ##        "The last k=$(orddivfact) Taylor coefficients ARE SET to 0.\n")
#     end
#     return orddivfact, cdivfact
# end
# Homogeneous coefficient for the division
function divHomogCoefN{T<:Number}( k::Int, ac::Array{T,1}, bc::Array{T,1},
    coeffs::Array{T,1}, ordLHopital::Int)
    #
    k == ordLHopital && return ac[ordLHopital+1] / bc[ordLHopital+1]
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    cc = zeros(T, numCoefk)
    @inbounds cc = mulHomogCoefN(k, coeffs, bc)
    @inbounds cc[1:end] = (ac[posI:posF]-cc[1:end]) / bc[ordLHopital+1]
    cc
end

## Division functions: rem and mod
for op in (:mod, :rem)
    @eval begin
        function ($op){T<:Real}(a::TaylorN{T}, x::Real)
            coeffs = a.coeffs
            coeffs[1] = ($op)(a.coeffs[1], x)
            return TaylorN( coeffs, a.order )
        end
    end
end
function mod2pi{T<:Real}(a::TaylorN{T})
    coeffs = a.coeffs
    coeffs[1] = mod2pi( a.coeffs[1] )
    return TaylorN( coeffs, a.order )
end

## Int power ##
function ^(a::TaylorN, n::Integer)
    uno = one(a)
    n < 0 && return uno / a^(-n)
    n == 0 && return uno
    if n%2 == 0     # even power
        n == 2 && return square(a)
        pow = div(n, 2)
        return square( a^pow )
    else            # odd power
        n == 1 && return a
        pow = div(n-1, 2)
        return a*square( a^pow )
    end
end
## Rational power ##
^(a::TaylorN,x::Rational) = a^(x.num/x.den)
## Real power ##
function ^(a::TaylorN, x::Real)
    uno = one(a)
    x == zero(x) && return uno
    x == 0.5 && return sqrt(a)
    order = a.order
    # First non-zero coefficient
    # ...
    @assert a.coeffs[1] != zero(a.coeffs[1])
    aux = ( a.coeffs[1] )^x
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    nCoefTot = sizeTable[end]
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = aux
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = powHomogCoefN(k, v, x, coeffs, 0)
    end
    TaylorN(coeffs,order)
end
function powHomogCoefN{T<:Number}(k::Int, ac::Array{T,1}, x::Real,
    coeffs::Array{T,1}, knull::Int)
    #
    k == knull && return (ac[ordLHopital+1])^x
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    cc = zeros(T, numCoefk)
    for ka = 0:k-1
        kb = k-ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds cc[pos] += ( x*kb-ka ) * ac[jb] * coeffs[ja]
            end
        end
    end
    @inbounds cc[1:end] = cc[1:end] / ( k * ac[1] )
    cc
end
^{T<:Number,S<:Number}(a::TaylorN{T}, x::Complex{S}) = exp( x*log(a) )
^(a::TaylorN, b::TaylorN) = exp( b*log(a) )

## Square ##
function square{T<:Number}(a::TaylorN{T})
    order = a.order
    nCoefTot = sizeTable[end]
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = a.coeffs[1]^2
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = squareHomogCoefN( k, a.coeffs )
    end
    TaylorN(coeffs,order)
end
# Homogeneous coefficients for square
function squareHomogCoefN{T<:Number}(k::Int, ac::Array{T,1})
    k == 0 && return (ac[1])^2
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    coeffs = zeros(T, numCoefk)
    kodd = k%2
    kend = div(k - 2 + kodd, 2)
    for ka = 0:kend
        kb = k - ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds coeffs[pos] += ac[ja] * ac[jb]
            end
        end
    end
    @inbounds coeffs[1:end] = 2coeffs[1:end]
    #
    if kodd == 0
        ka = div(k,2)
        @inbounds posaI, posaF = posAuxTable[end][ka]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posaI:posaF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds coeffs[pos] += ac[ja] * ac[jb]
            end
        end
    end
    coeffs
end

## sqrt ##
function sqrt(a::TaylorN)
    order = a.order
    nCoefTot = sizeTable[end]
    aux = sqrt(a.coeffs[1])
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = aux
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = sqrtHomogCoefN( k, v, coeffs, 0 )
    end
    TaylorN(coeffs,order)
end
# Homogeneous coefficients for sqrt
function sqrtHomogCoefN{T<:Number}(k::Int, ac::Array{T,1}, coeffs::Array{T,1}, knull::Int)
    k == 0 && return sqrt( ac[1] )
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    cc = zeros(T, numCoefk)
    kodd = k%2
    kend = div(k - 2 + kodd, 2)
    for ka = 1:kend
        kb = k - ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds cc[pos] += coeffs[ja] * coeffs[jb]
            end
        end
    end
    @inbounds cc[1:end] = 2cc[1:end]
    #
    if kodd == 0
        ka = div(k,2)
        @inbounds posaI, posaF = posAuxTable[end][ka]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posaI:posaF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds cc[pos] += coeffs[ja] * coeffs[jb]
            end
        end
    end
    @inbounds cc[1:end] = (ac[posI:posF] - cc[1:end]) / (2coeffs[1])
    cc
end

## exp ##
function exp(a::TaylorN)
    order = a.order
    nCoefTot = sizeTable[end]
    aux = exp(a.coeffs[1])
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = aux
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = expHomogCoefN( k, v, coeffs )
    end
    TaylorN(coeffs,order)
end
# Homogeneous coefficients for exp
function expHomogCoefN{T<:Number}(k::Int, ac::Array{T,1}, coeffs::Array{T,1})
    k == 0 && return exp( ac[1] )
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    cc = zeros(T, numCoefk)
    for ka = 0:k-1
        kb = k - ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds cc[pos] += kb * coeffs[ja] * ac[jb]
            end
        end
    end
    @inbounds cc[1:end] = cc[1:end] / k
    cc
end

## log ##
function log(a::TaylorN)
    order = a.order
    nCoefTot = sizeTable[end]
    aux = log(a.coeffs[1])
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    coeffs = zeros(T, nCoefTot)
    coeffs[1] = aux
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds coeffs[posI:posF] = logHomogCoefN( k, v, coeffs )
    end
    TaylorN(coeffs, order)
end
# Homogeneous coefficients for square
function logHomogCoefN{T<:Number}(k::Int, ac::Array{T,1}, coeffs::Array{T,1})
    k == 0 && return log( ac[1] )
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    cc = zeros(T, numCoefk)
    for ka = 1:k-1
        kb = k - ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds cc[pos] += ka * coeffs[ja] * ac[jb]
            end
        end
    end
    @inbounds cc[1:end] = ( ac[posI:posF] - cc[1:end] / k ) / ac[1]
    cc
end

## sin and cos ##
sin(a::TaylorN) = sincos(a)[1]
cos(a::TaylorN) = sincos(a)[2]
function sincos(a::TaylorN)
    order = a.order
    nCoefTot = sizeTable[end]
    aux = sin(a.coeffs[1])
    T = typeof(aux)
    v = convert(Array{T,1}, a.coeffs)
    sincoeffs = zeros(T,nCoefTot)
    coscoeffs = zeros(T,nCoefTot)
    sincoeffs[1] = aux
    coscoeffs[1] = cos( a.coeffs[1] )
    for k = 1:order
        @inbounds posI, posF = posAuxTable[end][k]
        @inbounds sincoeffs[posI:posF], coscoeffs[posI:posF] =
            sincosHomogCoefN(k, v, sincoeffs, coscoeffs)
    end
    return TaylorN( sincoeffs, order ), TaylorN( coscoeffs, order )
end
# Homogeneous coefficients for sin and cos
function sincosHomogCoefN{T<:Number}(k::Int, ac::Array{T,1}, scoeffs::Array{T,1}, ccoeffs::Array{T,1})
    k == 0 && return sin( ac[1] ), cos( ac[1] )
    @inbounds posI, posF = posAuxTable[end][k]
    numCoefk = posF - posI + 1
    sv = zeros(T, numCoefk)
    cv = zeros(T, numCoefk)
    for ka = 0:k-1
        kb = k - ka
        @inbounds posaI, posaF = posAuxTable[end][ka]
        @inbounds posbI, posbF = posAuxTable[end][kb]
        for ja = posaI:posaF
            @inbounds inda = indicesTable[end][ja]
            for jb = posbI:posbF
                @inbounds indb = indicesTable[end][jb]
                @inbounds pos = posTable[end][ inda + indb ]
                pos = pos - posI + 1
                @inbounds begin
                    x = kb * ac[jb]
                    sv[pos] += x * ccoeffs[ja]
                    cv[pos] -= x * scoeffs[ja]
                end
            end
        end
    end
    @inbounds sv[1:end] = sv[1:end] / k
    @inbounds cv[1:end] = cv[1:end] / k
    sv, cv
end

## Differentiation ##
"""Partial differentiation of a TaylorN series with respect to the r-th variable"""
function diffTaylor{T<:Number}(a::TaylorN{T}, r::Int)
    @assert 1 <= r <= NUMVARS[end]
    order = a.order
    nCoefTot = sizeTable[end]
    coeffs = zeros(T, nCoefTot)
    iIndices = zeros(Int, NUMVARS[end])
    for pos = 1:nCoefTot
        @inbounds iIndices[1:end] = indicesTable[end][pos]
        @inbounds n = iIndices[r]
        n == 0 && continue
        @inbounds iIndices[r] -= 1
        @inbounds posI = posTable[end][ iIndices ]
        @inbounds coeffs[posI] = n * a.coeffs[pos]
    end
    return TaylorN( coeffs, order )
end
diffTaylor(a::TaylorN) = diffTaylor(a, 1)

## TO BE DONE: Integration...

## Evaluates a Taylor polynomial on a given point ##
## NEEDS REVISION since yields only approx results
function evalTaylor{T<:Number,S<:Number}(a::TaylorN{T}, vals::Array{S,1} )
    numVars = NUMVARS[end]
    @assert length(vals) == numVars
    R = promote_type(T,S)
    suma = zero(R)
    nCoefTot = sizeTable[end]
    iIndices = zeros(Int, numVars)
    for pos = nCoefTot:-1:1
        @inbounds iIndices[1:end] = indicesTable[end][pos]
        @inbounds val = a.coeffs[pos]
        val == zero(R) && continue
        for k = 1:numVars
            @inbounds val *= vals[k]^iIndices[k]
        end
        suma += val
    end
    suma
end
evalTaylor{T<:Number}(a::TaylorN{T}) = evalTaylor(a, zeros(T, a.numVars) )


text:
function operator_pass(x, state)
    if x isa CSTParser.BinaryOpCall || x isa CSTParser.BinarySyntaxOpCall
        if CSTParser.precedence(x.op) in (8,13,14,16)
            ensure_no_space_after(x.arg1, state, state.offset)
            ensure_no_space_after(x.op, state, state.offset + x.arg1.fullspan)
        else
            ensure_single_space_after(x.arg1, state, state.offset)
            ensure_single_space_after(x.op, state, state.offset + x.arg1.fullspan)
        end
    elseif x isa CSTParser.WhereOpCall
        ensure_single_space_after(x.op, state, state.offset + x.arg1.fullspan)
        n = length(x.args)
        offset = state.offset + x.arg1.fullspan + x.op.fullspan
        for (i,a) in enumerate(x.args)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x isa CSTParser.EXPR{CSTParser.ColonOpCall}
        offset = state.offset
        n = length(x.args)
        for (i,a) in enumerate(x.args)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x isa CSTParser.EXPR{CSTParser.ChainOpCall} || x isa CSTParser.EXPR{CSTParser.Comparison}
        offset = state.offset
        n = length(x.args)
        for (i,a) in enumerate(x.args)
            if i != n
                ensure_single_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function tuple_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.TupleH}
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if a isa CSTParser.PUNCTUATION && a.kind == Tokens.COMMA && i !=n && !(x.args[i+1] isa CSTParser.PUNCTUATION)
                ensure_single_space_after(a, state, offset)
            elseif i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function curly_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.Curly}
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function call_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.Call}
        offset = state.offset + x.args[1].fullspan
        n = length(x)
        for (i, a) in enumerate(x)
            i == 1 && continue
            if a isa CSTParser.PUNCTUATION && a.kind == Tokens.COMMA
                ensure_single_space_after(a, state, offset)
            # elseif a isa CSTParser.EXPR{CSTParser.Parameters}
            elseif i != n && !(x.args[i + 1] isa CSTParser.EXPR{CSTParser.Parameters})
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x isa CSTParser.EXPR{CSTParser.Kw}
        ensure_single_space_after(x.args[1], state, state.offset)
        ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
    end
end


text:
function build_checkerboard(mc::AbstractDQMC{CBGeneric})
  l = mc.l

  l.groups = UnitRange[]
  edges_used = zeros(Int64, l.n_bonds)
  l.checkerboard = zeros(3, l.n_bonds)
  group_start = 1
  group_end = 1

  while minimum(edges_used) == 0
    sites_used = zeros(Int64, l.sites)

    for id in 1:l.n_bonds
      src, trg, typ = l.bonds[id,1:3]

      if edges_used[id] == 1 continue end
      if sites_used[src] == 1 continue end
      if sites_used[trg] == 1 continue end

      edges_used[id] = 1
      sites_used[src] = 1
      sites_used[trg] = 1

      l.checkerboard[:, group_end] = [src, trg, id]
      group_end += 1
    end
    push!(l.groups, group_start:group_end-1)
    group_start = group_end
  end
  l.n_groups = length(l.groups)
end

# helper to cutoff numerical zeros
rem_eff_zeros!(X::AbstractArray) = map!(e -> abs.(e) < 1e-15 ? zero(e) : e,X,X)

function init_checkerboard_matrices(mc::AbstractDQMC{CBGeneric})
  l = mc.l
  p = mc.p
  H = heltype(mc)

  println("Initializing hopping exponentials (Checkerboard, generic)")

  build_checkerboard(mc)

  n_groups = l.n_groups
  eT_half = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2) # group, spin (up, down), flavor (x, y)
  eT_half_inv = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2)
  eT = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2)
  eT_inv = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2)

  U = [0.0, 1.0]
  R = [1.0, 0.0]
  UR = [1.0, 1.0]
  DR = [1.0, -1.0]
  UU = [0.0, 2.0]
  RR = [2.0, 0.0]

  for f in 1:2
    for s in 1:2
      for (g, gr) in enumerate(l.groups)
        # Build hopping matrix of individual chkr group
        T = zeros(H, l.sites, l.sites)

        for i in gr
          src = l.checkerboard[1,i]
          trg = l.checkerboard[2,i]
          bond = l.checkerboard[3,i]
          v = l.bond_vecs[bond,:]

          # nn
          if v == U && p.hoppings != "none"
            T[trg, src] = T[src, trg] += -l.t[2,f]
          elseif v == R
            T[trg, src] = T[src, trg] += -l.t[1,f]

          #Nnn
          elseif (v == UR || v == DR) && p.Nhoppings != "none"
            T[trg, src] = T[src, trg] += -l.tN[1,f]

          #NNnn
          elseif v == UU && p.NNhoppings != "none"
            T[trg, src] = T[src, trg] += -l.tNN[2,f]
          elseif v == RR
            T[trg, src] = T[src, trg] += -l.tNN[1,f]
          end
        end

        eT_half[g,s,f] = sparse(rem_eff_zeros!(exp(- 0.5 * p.delta_tau * T)))
        eT_half_inv[g,s,f] = sparse(rem_eff_zeros!(exp(0.5 * p.delta_tau * T)))
        eT[g,s,f] = sparse(rem_eff_zeros!(exp(- p.delta_tau * T)))
        eT_inv[g,s,f] = sparse(rem_eff_zeros!(exp(p.delta_tau * T)))
      end
    end
  end

  if p.opdim == 3
    l.chkr_hop_half = [cat(eT_half[g,1,1], eT_half[g,2,2], eT_half[g,2,1], eT_half[g,1,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_half_inv = [cat(eT_half_inv[g,1,1], eT_half_inv[g,2,2], eT_half_inv[g,2,1], eT_half_inv[g,1,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop = [cat(eT[g,1,1], eT[g,2,2], eT[g,2,1], eT[g,1,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_inv = [cat(eT_inv[g,1,1], eT_inv[g,2,2], eT_inv[g,2,1], eT_inv[g,1,2], dims=(1,2)) for g in 1:n_groups]

  else # O(2) and O(1) model
    l.chkr_hop_half = [cat(eT_half[g,1,1], eT_half[g,2,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_half_inv = [cat(eT_half_inv[g,1,1], eT_half_inv[g,2,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop = [cat(eT[g,1,1], eT[g,2,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_inv = [cat(eT_inv[g,1,1], eT_inv[g,2,2], dims=(1,2)) for g in 1:n_groups]
  end
  l.chkr_hop_half_dagger = adjoint.(l.chkr_hop_half)
  l.chkr_hop_dagger = adjoint.(l.chkr_hop)

  muv = vcat(fill(p.mu1,l.sites),fill(p.mu2,l.sites))
  muv = repeat(muv, outer=[Int(p.flv/2)])
  l.chkr_mu_half = sparse(Diagonal(exp.(-0.5*p.delta_tau * -muv)))
  l.chkr_mu_half_inv = sparse(Diagonal(exp.(0.5*p.delta_tau * -muv)))
  l.chkr_mu = sparse(Diagonal(exp.(-p.delta_tau * -muv)))
  l.chkr_mu_inv = sparse(Diagonal(exp.(p.delta_tau * -muv)))

  fold_chkr_grps(mc)

  hop_mat_exp_chkr = foldl(*,l.chkr_hop_half) * sqrt.(l.chkr_mu)
  r = effreldiff(l.hopping_matrix_exp,hop_mat_exp_chkr)
  r[findall(x -> x == zero(x), hop_mat_exp_chkr)] .= 0
  println("Checkerboard (generic) - exact (abs):\t\t", maximum(absdiff(l.hopping_matrix_exp,hop_mat_exp_chkr)))
end

function fold_chkr_grps(mc::AbstractDQMC{CBGeneric})
  # ### Heuristic folding (merging) of chkr groups.
  #
  # We often have many sparse chkr groups which stay sparse, when folded (multiplied together).
  # It is advantageous to reduce the number of sparse matrices that we multiply all the time in
  # multiply_B_*! functions.
  #
  # The applied heuristic is as follows:
  # Fold chkr matrices until a SPARSITY_LIMIT is reached. Save the result as one "folded" matrix.
  # Continue in the same spirit with remaining chkr matrices until we visited them all.
  #
  # Massive(!) speed improvement!!
  l = mc.l
  H = heltype(mc)
  flv = mc.p.flv
  N = mc.l.sites
  SPARSITY_LIMIT = mc.p.sparsity_limit

  @show SPARSITY_LIMIT

  l.chkr_hop_half_folded = SparseMatrixCSC{H, Int64}[]
  cur = SparseMatrixCSC{H}(I, flv*N, flv*N)
  l.folded = UnitRange{Int64}[]
  fstart = 2
  fstop = 1

  for i in 2:l.n_groups
    if i != 2 && (count(!iszero, cur)/length(cur)) > SPARSITY_LIMIT
      push!(l.chkr_hop_half_folded, cur)
      push!(l.folded, fstart:fstop)
      fstart = i
      fstop = i
      cur = SparseMatrixCSC{H}(I, flv*N, flv*N)
    else
      fstop += 1
    end

    cur = l.chkr_hop_half[i] * cur
  end

  push!(l.chkr_hop_half_folded, cur)
  push!(l.folded, fstart:fstop)

  l.n_folded = length(l.folded)

  l.chkr_hop_half_inv_folded = Vector{SparseMatrixCSC{H, Int64}}(undef, l.n_folded)
  l.chkr_hop_half_dagger_folded = Vector{SparseMatrixCSC{H, Int64}}(undef, l.n_folded)
  l.chkr_hop_half_folded_rev = Vector{SparseMatrixCSC{H, Int64}}(undef, l.n_folded)
  l.chkr_hop_half_inv_folded_rev = Vector{SparseMatrixCSC{H, Int64}}(undef, l.n_folded)
  l.chkr_hop_half_dagger_folded_rev = Vector{SparseMatrixCSC{H, Int64}}(undef, l.n_folded)

  for (i, rng) in enumerate(l.folded)
    cur_rev = SparseMatrixCSC{H}(I, flv*N, flv*N)
    cur_inv = SparseMatrixCSC{H}(I, flv*N, flv*N)
    cur_inv_rev = SparseMatrixCSC{H}(I, flv*N, flv*N)
    cur_dagger = SparseMatrixCSC{H}(I, flv*N, flv*N)
    cur_dagger_rev = SparseMatrixCSC{H}(I, flv*N, flv*N)
    for k in rng
      cur_rev = cur_rev * l.chkr_hop_half[k]
      cur_inv = l.chkr_hop_half_inv[k] * cur_inv
      cur_inv_rev = cur_inv_rev * l.chkr_hop_half_inv[k]
      cur_dagger = l.chkr_hop_half_dagger[k] * cur_dagger
      cur_dagger_rev = cur_dagger_rev * l.chkr_hop_half_dagger[k]
    end
    l.chkr_hop_half_folded_rev[i] = cur_rev
    l.chkr_hop_half_inv_folded[i] = cur_inv
    l.chkr_hop_half_inv_folded_rev[i] = cur_inv_rev
    l.chkr_hop_half_dagger_folded[i] = cur_dagger
    l.chkr_hop_half_dagger_folded_rev[i] = cur_dagger_rev
  end

  @show l.n_groups
  @show l.n_folded
  nothing
end

#### WITH ARTIFICIAL B-FIELD

function init_checkerboard_matrices_Bfield(mc::AbstractDQMC{CBGeneric})
  l = mc.l
  p = mc.p
  H = heltype(mc)

  println("Initializing hopping exponentials (Bfield, Checkerboard, generic)")

  build_checkerboard(mc)

  n_groups = l.n_groups
  eT_half = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2) # group, spin (up, down), flavor (x, y)
  eT_half_inv = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2)
  eT = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2)
  eT_inv = Array{SparseMatrixCSC{H, Int}, 3}(undef, n_groups,2,2)

  B = zeros(2,2) # rowidx = spin up,down, colidx = flavor
  if p.Bfield
    B[1,1] = B[2,2] = 2 * pi / l.sites
    B[1,2] = B[2,1] = - 2 * pi / l.sites
  end

  U = [0.0, 1.0]
  R = [1.0, 0.0]
  UR = [1.0, 1.0]
  DR = [1.0, -1.0]
  UU = [0.0, 2.0]
  RR = [2.0, 0.0]

  for f in 1:2
    for s in 1:2
      for (g, gr) in enumerate(l.groups)
        # Build hopping matrix of individual chkr group
        T = zeros(H, l.sites, l.sites)

        for i in gr
          src = l.checkerboard[1,i]
          trg = l.checkerboard[2,i]
          bond = l.checkerboard[3,i]
          v = l.bond_vecs[bond,:]

          if v == U && p.hoppings != "none"
            T[trg, src] += - exp(im * l.peirls[s,f][trg,src]) * l.t[2,f]
            T[src, trg] += - exp(im * l.peirls[s,f][src,trg]) * l.t[2,f]
          elseif v == R && p.hoppings != "none"
            T[trg, src] += - exp(im * l.peirls[s,f][trg,src]) * l.t[1,f]
            T[src, trg] += - exp(im * l.peirls[s,f][src,trg]) * l.t[1,f]

          #Nnn
          elseif (v == UR || v == DR) && p.Nhoppings != "none"
            T[trg, src] += - exp(im * l.peirls[s,f][trg,src]) * l.tN[1,f]
            T[src, trg] += - exp(im * l.peirls[s,f][src,trg]) * l.tN[1,f]

          #NNnn
          elseif v == UU  && p.NNhoppings != "none"
            T[trg, src] += - exp(im * l.peirls[s,f][trg,src]) * l.tNN[2,f]
            T[src, trg] += - exp(im * l.peirls[s,f][src,trg]) * l.tNN[2,f]
          elseif v == RR && p.NNhoppings != "none"
            T[trg, src] += - exp(im * l.peirls[s,f][trg,src]) * l.tNN[1,f]
            T[src, trg] += - exp(im * l.peirls[s,f][src,trg]) * l.tNN[1,f]
          end
        end

        eT_half[g,s,f] = sparse(rem_eff_zeros!(exp(- 0.5 * p.delta_tau * T)))
        eT_half_inv[g,s,f] = sparse(rem_eff_zeros!(exp(0.5 * p.delta_tau * T)))
        eT[g,s,f] = sparse(rem_eff_zeros!(exp(- p.delta_tau * T)))
        eT_inv[g,s,f] = sparse(rem_eff_zeros!(exp(p.delta_tau * T)))
      end
    end
  end

  if p.opdim == 3
    l.chkr_hop_half = [cat(eT_half[g,1,1], eT_half[g,2,2], eT_half[g,2,1], eT_half[g,1,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_half_inv = [cat(eT_half_inv[g,1,1], eT_half_inv[g,2,2], eT_half_inv[g,2,1], eT_half_inv[g,1,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop = [cat(eT[g,1,1], eT[g,2,2], eT[g,2,1], eT[g,1,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_inv = [cat(eT_inv[g,1,1], eT_inv[g,2,2], eT_inv[g,2,1], eT_inv[g,1,2], dims=(1,2)) for g in 1:n_groups]

  else # O(2) and O(1) model
    l.chkr_hop_half = [cat(eT_half[g,1,1], eT_half[g,2,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_half_inv = [cat(eT_half_inv[g,1,1], eT_half_inv[g,2,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop = [cat(eT[g,1,1], eT[g,2,2], dims=(1,2)) for g in 1:n_groups]
    l.chkr_hop_inv = [cat(eT_inv[g,1,1], eT_inv[g,2,2], dims=(1,2)) for g in 1:n_groups]
  end
  l.chkr_hop_half_dagger = adjoint.(l.chkr_hop_half)
  l.chkr_hop_dagger = adjoint.(l.chkr_hop)

  muv = vcat(fill(p.mu1,l.sites),fill(p.mu2,l.sites))
  muv = repeat(muv, outer=[Int(p.flv/2)])
  l.chkr_mu_half = sparse(Diagonal(exp.(-0.5*p.delta_tau * -muv)))
  l.chkr_mu_half_inv = sparse(Diagonal(exp.(0.5*p.delta_tau * -muv)))
  l.chkr_mu = sparse(Diagonal(exp.(-p.delta_tau * -muv)))
  l.chkr_mu_inv = sparse(Diagonal(exp.(p.delta_tau * -muv)))

  fold_chkr_grps(mc)

  hop_mat_exp_chkr = foldl(*,l.chkr_hop_half) * sqrt.(l.chkr_mu)
  r = effreldiff(l.hopping_matrix_exp,hop_mat_exp_chkr)
  r[findall(x -> x == zero(x), hop_mat_exp_chkr)] .= 0
  println("Checkerboard (Bfield, generic) - exact (abs):\t\t", maximum(absdiff(l.hopping_matrix_exp,hop_mat_exp_chkr)))
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2019 Michael Reed

function orbit(E,f::Function,bi::Matrix{Float64},orb::Int=0,depth::Int=1,incr::Int=384)
    x,N,N2,orbit,bis = real_orb(E,f,bi,orb,depth,incr)
    # trim graph
    d=1.07
    xl = bi[1:2]
    yl = [minimum([d*minimum(N[:,2]),0]),maximum([d*maximum(N[:,2]),0])]
    # plot background lines
    plt = UnicodePlots.lineplot(x[:],N[:,1],xlim=xl,ylim=yl,color=:magenta,name="y=real(z)")
    UnicodePlots.lineplot!(plt,x[:],N[:,2],color=:blue,name="ϕ(real(z))")
    # plot orbit cobweb path
    UnicodePlots.scatterplot!(plt,orbit[:,1],orbit[:,2],color=:red,name="(zₙ,ϕ(zₙ))")
    # set legend
    leg = ["ϕ^{$x}(z)" for x ∈ 2:depth]
    # plot f^2,f^3,f^4,...
    for h ∈ 3:depth+1
        UnicodePlots.lineplot!(plt,x,N[:,h],color=isodd(h) ? :yellow : :green,name=leg[h-2])
    end
    if ~(orb == 0)
        UnicodePlots.lineplot!(plt,collect(range(bi[1],stop=bi[2],length=length(N2))),N2[:],color=:cyan,name="ϕ(z_{0:$orb})")
        funt = ", IC: z₀ = $(bis[3]), n∈0:$orb"
    else
        funt = ""
    end
    # set title
    UnicodePlots.title!(plt,"z ↦ $E$funt")
    return plt
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ExplicitStructure(connectivities)

A data structure that encodes all incidence relations
of a mesh explicitly as a list of `connectivities`.

See also `[`Connectivity`](@ref).
"""
struct ExplicitStructure{C<:Connectivity} <: TopologicalStructure
  connectivities::Vector{C}
end

connectivities(s::ExplicitStructure) = s.connectivities

==(s1::ExplicitStructure, s2::ExplicitStructure) =
  Set(s1.connectivities) == Set(s2.connectivities)


text:
export HTTPServer

struct HTTPServer
    io::Base.IOServer
    handler::Function
    # options::Dict{Symbol,Any}
end

function createpipelinehandler(opts)
    opts = Dict{Symbol, Any}(opts)
    fn = function(socket::IO)
        connection = Connection(socket)
        allowedopts = [
            :maxheader,
            :maxuri,
            :maxmethod,
            :maxversion,
            :maxbody,
            :maxtrailer,
            :keepcomments,
            :strict
        ]
        processopts = Vector{Pair{Symbol,Any}}()
        for opt in allowedopts
            if haskey(opts, opt)
                processopts[opt] = opts[opt]
            end
        end
        processrequest!(connection; processopts...)
        # req = Request(connection)
        # res = Response(connection)
        next(connection)

        # TODO: make req + res
    end
end

function HTTPServer(io::Base.IOServer, p::Pipeline; options...)

    fallthroughhandler = function(conn)
        # TODO: content negotation
        throw(HTTPError(404; message="Not Found"))
    end
    errorhandler = function(error, conn) # TODO: req, res
        shouldclose = false
        if !isa(error, HTTPError)
            error = HTTPError(500)
        end

        status_code = error.code
        status_message = STATUS_CODES[status_code]
        safe = error.safe
        message = safe ? error.message : status_message
        shoudclose = error.shouldclose
        contentlength = length(message)

        headers = Headers(
            "Content-Length" => string(contentlength),
            "Content-Type" => "text/plain;charset=utf-8"
        )

        writeresponse!(conn, status_code, status_message, headers, IOBuffer(message))

        if shouldclose
            close(conn.socket)
        end

        fin(conn)
    end
    p = pipe(p, fallthroughhandler)
    p = recover(p, errorhandler)
    p = pipe(createpipelinehandler(options), p)
    handler = function(socket::IO)
        Distributed.@spawn begin
            println("running handler")
            run(p, socket)
            true
        end
    end

    HTTPServer(io, handler)
end

function Base.run(server::HTTPServer)
    while isopen(server.io)
        println("ran")
        sock = accept(server.io)
        println("new socket")
        begin
            while isopen(sock)
                res = resolve_future(server.handler(sock))
                println("open $(isopen(sock))")
            end
            println("socket closed $(isopen(sock))")
        end

        println("Server open? $(isopen(server.io))")
    end
end


text:
module Kegg

# include -
include("Include.jl")

# export methods -
export get_reactions_for_ec_number
export get_ec_number_for_gene
export get_sequence_for_gene
export get_sequence_for_protein
export get_organism_codes
export build_metabolite_matching_table

# export types -
export Reaction
export Sequence
export Organism

end # module


text:
module TestFromComposite

using Test
using Tables
using MLJBase
using ..Models
using ..TestUtilities
using CategoricalArrays
using StableRNGs
using Parameters
rng = StableRNG(616161)

ridge_model = FooBarRegressor(lambda=0.1)
selector_model = FeatureSelector()

## FROM_NETWORK_PREPROCESS

# supervised:
Xs = source(nothing)
ys = source(nothing)
z  = log(ys)
stand = UnivariateStandardizer()
standM = machine(stand, z)
u = transform(standM, z)
hot = OneHotEncoder()
hotM = machine(hot, Xs)
W = transform(hotM, Xs)
knn = KNNRegressor()
knnM = machine(knn, W, u)
oak = DecisionTreeRegressor()
oakM = machine(oak, W, u)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)

mach_ex = :(machine(Deterministic(), Xs, ys; predict=yhat))

## TESTING `from_network_preprocess`

ex = Meta.parse(
    "begin
         mutable struct CompositeX
             knn_rgs=knn
             one_hot_enc=hot
         end
         target_scitype=AbstractVector{<:Continuous}
         input_scitype=Table(Continuous,Multiclass)
     end")
mach_, modeltype_ex, struct_ex, no_fields, dic =
    MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex)
eval(Parameters.with_kw(struct_ex, TestFromComposite, false))
@test supertype(CompositeX) == DeterministicComposite
composite = CompositeX()
@test composite.knn_rgs == knn
@test composite.one_hot_enc == hot
@test dic[:target_scitype] == :(AbstractVector{<:Continuous})
@test dic[:input_scitype] == :(Table(Continuous, Multiclass))

ex = Meta.parse(
    "begin
         mutable struct Composite4 <: ProbabilisticComposite
             knn_rgs=knn
             one_hot_enc=hot
         end
     end")
mach_, modeltype_ex, struct_ex =
    MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex)
eval(Parameters.with_kw(struct_ex, TestFromComposite, false))
@test supertype(Composite4) == ProbabilisticComposite

ex = Meta.parse(
    "mutable struct Composite2
        knn_rgs=knn
        one_hot_enc=hot
     end")
mach_, modeltype_ex, struct_ex, no_fields, dic =
    MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex)
eval(Parameters.with_kw(struct_ex, TestFromComposite, false))
composite = Composite2()
@test composite.knn_rgs == knn
@test composite.one_hot_enc == hot

ex = Meta.parse(
    "begin
         mutable struct Composite6 <: Probabilistic
             knn_rgs=knn
             one_hot_enc=hot
         end
     end")
@test_logs((:warn, r"New composite"),
           MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex))

ex = Meta.parse(
    "begin
         mutable struct Composite20
             knn_rgs=knn
             one_hot_enc=hot
         end
         target_scitype == Continuous
     end")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex))

ex = Meta.parse(
    "begin
         mutable struct Composite20
             knn_rgs=knn
             one_hot_enc=hot
         end
         Continuous
     end")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex))

ex = Meta.parse(
    "begin
         mutable struct Composite20
             knn_rgs=knn
             one_hot_enc=hot
         end
         43 = Continuous
     end")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex))

ex = Meta.parse(
    "begin
         mutable struct Composite7 < Probabilistic
             knn_rgs=knn
             one_hot_enc=hot
         end
     end")
@test_throws(ArgumentError,
           MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex))

@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestFromComposite, knn, ex))

ex = Meta.parse(
    "begin
         Composite3(
             knn_rgs=knn,
             one_hot_enc=hot)
     end")
@test_throws(ArgumentError,
             MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex))

ex = Meta.parse(
    "begin
         mutable struct Composite8
             knn_rgs::KNNRegressor=knn
             one_hot_enc=hot
         end
     end")
mach_, modeltype_ex, struct_ex =
    MLJBase.from_network_preprocess(TestFromComposite, mach_ex, ex)
eval(Parameters.with_kw(struct_ex, TestFromComposite, false))
VERSION ≥ v"1.3.0-" &&
    @test fieldtypes(Composite8) == (KNNRegressor, Any)

# test that you cannot leave "default" component models unspecified:
modeltype_ex = :Composite9
struct_ex = :(mutable struct Composite9 <: DeterministicComposite
              knn_rgs::KNNRegressor
              one_hot_enc = hot
              end)
@test_logs (:error, r"Problem instantiating") begin
    @test_throws Exception begin
        MLJBase.from_network_(TestFromComposite,
                              mach_ex, modeltype_ex,
                              struct_ex, false, Dict{Symbol,Any}())
    end
end

## TEST MACRO-EXPORTED  NETWORKS
# (CANNOT WRAP IN @testset)

# some actual data:
N = 10
X = MLJBase.table(rand(N, 3))
y = rand(N)
w = rand(N)

# supervised with sample weights:
ws = source()
knnM = machine(knn, W, u, ws)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)

@from_network machine(Deterministic(), Xs, ys, ws; predict=yhat) begin
    mutable struct CompositeX1
        knn_rgs=knn
        one_hot_enc=hot
    end
    supports_weights = true
    target_scitype = AbstractVector{<:Continuous}
end
model = CompositeX1()
@test supports_weights(model)
@test target_scitype(model) == AbstractVector{<:Continuous}
predict(fit!(machine(model, X, y, w), verbosity=-1), X);

# unsupervised:
@from_network machine(Unsupervised(), Xs; transform=W) begin
    mutable struct CompositeX2
        one_hot_enc=hot
    end
end
model = CompositeX2()
transform(fit!(machine(model, X), verbosity=-1), X)

# second supervised test:
fea = FeatureSelector()
feaM = machine(fea, Xs)
G = transform(feaM, Xs)
hotM = machine(hot, G)
H = transform(hotM, G)
elm = DecisionTreeClassifier()
elmM = machine(elm, H, ys)
yhat = predict(elmM, H)

@from_network machine(Probabilistic(), Xs, ys; predict=yhat) begin
    mutable struct CompositeX3
        selector=fea
        one_hot=hot
        tree=elm
    end
end
model = CompositeX3()
y = coerce(y, Multiclass)
@test predict(fit!(machine(model, X, y), verbosity=-1), X) isa
    AbstractVector{<:UnivariateFinite}

# yet more examples:
x1 = map(n -> mod(n,3), rand(rng, UInt8, 100)) |> categorical;
x2 = randn(rng, 100);
X = (x1=x1, x2=x2);
y = x2.^2;

Xs = source(X)
ys = source(y)
z = log(ys)
stand = UnivariateStandardizer()
standM = machine(stand, z)
u = transform(standM, z)
hot = OneHotEncoder()
hotM = machine(hot, Xs)
W = transform(hotM, Xs)
knn = KNNRegressor()
knnM = machine(knn, W, u)
oak = DecisionTreeRegressor()
oakM = machine(oak, W, u)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)

mach = machine(Deterministic(), Xs, ys; predict=yhat)

@from_network mach begin
    mutable struct Composite10
        knn_rgs::KNNRegressor=knn
        one_hot_enc=hot
    end
end

model_ = Composite10()

mach = machine(model_, X, y)

@test_model_sequence(fit_only!(mach),
                     [(:train, model_), (:train, stand), (:train, hot),
                      (:train, knn), (:train, oak)],
                     [(:train, model_), (:train, hot), (:train, stand),
                      (:train, knn), (:train, oak)],
                     [(:train, model_), (:train, stand), (:train, hot),
                      (:train, oak), (:train, knn)],
                     [(:train, model_), (:train, hot), (:train, stand),
                      (:train, oak), (:train, knn)])

model_.knn_rgs.K = 55
knn = model_.knn_rgs
@test_model_sequence(fit_only!(mach),
                     [(:update, model_), (:skip, stand), (:skip, hot),
                      (:update, knn), (:skip, oak)],
                     [(:update, model_), (:skip, hot), (:skip, stand),
                      (:update, knn), (:skip, oak)],
                     [(:update, model_), (:skip, stand), (:skip, hot),
                      (:skip, oak), (:update, knn)],
                     [(:update, model_), (:skip, hot), (:skip, stand),
                      (:skip, oak), (:update, knn)])

@test MLJBase.tree(mach.fitresult.predict).arg1.arg1.arg1.arg1.model.K == 55

# check data anonymity:
@test all(x->(x===nothing),
          [s.data for s in sources(mach.fitresult.predict)])

multistand = Standardizer()
multistandM = machine(multistand, W)
W2 = transform(multistandM, W)

mach = machine(Unsupervised(), Xs; transform=W2)

@from_network mach begin
    mutable struct MyTransformer
        one_hot=hot
    end
end

model_ = MyTransformer()

mach = machine(model_, X)
@test_model_sequence fit_only!(mach) [(:train, model_),
                                      (:train, hot), (:train, multistand)]

model_.one_hot.drop_last=true
hot = model_.one_hot
@test_model_sequence fit_only!(mach) [(:update, model_),
                                      (:update, hot), (:train, multistand)]

# check nested fitted_params:
FP = MLJBase.fitted_params(mach)
@test keys(FP) == (:one_hot, :machines, :fitted_params_given_machine)
@test Set(FP.one_hot.fitresult.all_features) == Set(keys(X))

# check data anomynity:
@test all(x->(x===nothing),
          [s.data for s in sources(mach.fitresult.transform)])

transform(mach, X);

## TEST MACRO-EXPORTED SUPERVISED NETWORK WITH SAMPLE WEIGHTS

rng = StableRNG(56161)
N = 500
X = (x = rand(rng, 3N), );
y = categorical(rand(rng, "abc", 3N));
# define class weights :a, :b, :c in ration 2:4:1
w = map(y) do η
    if η == 'a'
        return 2
    elseif η == 'b'
        return 4
    else
        return 1
    end
end;
Xs = source(X)
ys = source(y)
ws = source(w)

standM = machine(Standardizer(), Xs)
W = transform(standM, Xs)

rgs = ConstantClassifier() # supports weights
rgsM = machine(rgs, W, ys, ws)
yhat = predict(rgsM, W)

fit!(yhat, verbosity=0)
fit!(yhat, rows=1:div(N,2), verbosity=0)
yhat(rows=1:div(N,2));

mach = machine(Probabilistic(), Xs, ys, ws; predict=yhat)

@from_network mach begin
    mutable struct MyComposite
        regressor=rgs
    end
    supports_weights=true
end

my_composite = MyComposite()
@test MLJBase.supports_weights(my_composite)
mach = fit!(machine(my_composite, X, y), verbosity=0)
Xnew = selectrows(X, 1:div(N,2))
predict(mach, Xnew)[1]
posterior = predict(mach, Xnew)[1]

# "posterior" is roughly uniform:
@test abs(pdf(posterior, 'b')/(pdf(posterior, 'a'))  - 1) < 0.15
@test abs(pdf(posterior, 'b')/(pdf(posterior, 'c'))  - 1) < 0.15

# now add weights:
mach = fit!(machine(my_composite, X, y, w), rows=1:div(N,2), verbosity=0)
posterior = predict(mach, Xnew)[1]

# "posterior" is skewed appropriately in weighted case:
@test abs(pdf(posterior, 'b')/(2*pdf(posterior, 'a'))  - 1) < 0.15
@test abs(pdf(posterior, 'b')/(4*pdf(posterior, 'c'))  - 1) < 0.19

# composite with no fields:
mach = machine(Probabilistic(), Xs, ys, ws; predict=yhat)
@from_network mach begin
    struct CompositeWithNoFields
    end
end
composite_with_no_fields = CompositeWithNoFields()
mach = fit!(machine(composite_with_no_fields, X, y), verbosity=0)

## EXPORTING A TRANSFORMER WITH PREDICT AND TRANSFORM

# A dummy clustering model:
mutable struct DummyClusterer <: Unsupervised
    n::Int
end
DummyClusterer(; n=3) = DummyClusterer(n)
function MLJBase.fit(model::DummyClusterer, verbosity::Int, X)
    Xmatrix = Tables.matrix(X)
    n = min(size(Xmatrix, 2), model.n)
    centres = Xmatrix[1:n, :]
    levels = categorical(1:n)
    report = (centres=centres,)
    fitresult = levels
    return fitresult, nothing, report
end
MLJBase.transform(model::DummyClusterer, fitresult, Xnew) =
    selectcols(Xnew, 1:length(fitresult))
MLJBase.predict(model::DummyClusterer, fitresult, Xnew) =
    [fill(fitresult[1], nrows(Xnew))...]

N = 20
X = (a = rand(N), b = categorical(rand("FM", N)))

Xs = source(X)
W = transform(machine(OneHotEncoder(), Xs), Xs)
clust = DummyClusterer(n=2)
m = machine(clust, W)
yhat = predict(m, W)
Wout = transform(m, W)
foo = first(yhat)
mach = machine(Unsupervised(), Xs;
               predict=yhat,
               transform=Wout,
               report=(foo=foo,))

@from_network mach begin
    mutable struct WrappedClusterer
        clusterer::Unsupervised = clust
    end
    input_scitype = Table(Continuous,Multiclass)
end

model = WrappedClusterer()
mach = fit!(machine(model, X), verbosity=0)
fit!(yhat)
@test predict(mach, X) == yhat()
@test transform(mach, X).a ≈ Wout().a
rep = report(mach)
@test rep.foo == yhat() |> first

## EXPORTING A STATIC LEARNING NETWORK (NO TRAINING ARGUMENTS)

age = [23, 45, 34, 25, 67]
X = (age = age,
     gender = categorical(['m', 'm', 'f', 'm', 'f']))

struct MyStaticTransformer <: Static
    ftr::Symbol
end

MLJBase.transform(transf::MyStaticTransformer, verbosity, X) =
    selectcols(X, transf.ftr)

Xs = source()
W = transform(machine(MyStaticTransformer(:age)), Xs)
Z = 2*W

@from_network machine(Static(), Xs; transform=Z) begin
    struct NoTraining
    end
end

mach = fit!(machine(NoTraining()), verbosity=0)
@test transform(mach, X) == 2*X.age

## TESTINGS A STACK AND IN PARTICULAR FITTED_PARAMS

folds(data, nfolds) =
    partition(1:nrows(data), (1/nfolds for i in 1:(nfolds-1))...);

model1 = RidgeRegressor()
model2 = KNNRegressor(K=1)
judge = KNNRegressor(K=1)

X = source()
y = source()

folds(X::AbstractNode, nfolds) = node(XX->folds(XX, nfolds), X)
MLJBase.restrict(X::AbstractNode, f::AbstractNode, i) =
    node((XX, ff) -> restrict(XX, ff, i), X, f);
MLJBase.corestrict(X::AbstractNode, f::AbstractNode, i) =
    node((XX, ff) -> corestrict(XX, ff, i), X, f);

f = folds(X, 3)

m11 = machine(model1, corestrict(X, f, 1), corestrict(y, f, 1))
m12 = machine(model1, corestrict(X, f, 2), corestrict(y, f, 2))
m13 = machine(model1, corestrict(X, f, 3), corestrict(y, f, 3))

y11 = predict(m11, restrict(X, f, 1));
y12 = predict(m12, restrict(X, f, 2));
y13 = predict(m13, restrict(X, f, 3));

m21 = machine(model2, corestrict(X, f, 1), corestrict(y, f, 1))
m22 = machine(model2, corestrict(X, f, 2), corestrict(y, f, 2))
m23 = machine(model2, corestrict(X, f, 3), corestrict(y, f, 3))

y21 = predict(m21, restrict(X, f, 1));
y22 = predict(m22, restrict(X, f, 2));
y23 = predict(m23, restrict(X, f, 3));

y1_oos = vcat(y11, y12, y13);
y2_oos = vcat(y21, y22, y23);

X_oos = MLJBase.table(hcat(y1_oos, y2_oos))

m_judge = machine(judge, X_oos, y)

m1 = machine(model1, X, y)
m2 = machine(model2, X, y)

y1 = predict(m1, X);
y2 = predict(m2, X);

X_judge = MLJBase.table(hcat(y1, y2))
yhat = predict(m_judge, X_judge)

@from_network machine(Deterministic(), X, y; predict=yhat) begin
    mutable struct MyStack
        regressor1=model1
        regressor2=model2
        judge=judge
    end
end

my_stack = MyStack()
X, y = make_regression(18, 2)
mach = machine(my_stack, X, y)
fit!(mach, verbosity=0)

fp = fitted_params(mach)
@test keys(fp.judge) == (:tree,)
@test length(fp.regressor1) == 4
@test length(fp.regressor2) == 4
@test keys(fp.regressor1[1]) == (:coefficients, :intercept)
@test keys(fp.regressor2[1]) == (:tree,)

## ISSUE #377

target_stand = Standardizer()
stand = Standardizer()
rgs = KNNRegressor()

X = source()
y = source()

mach1 = machine(target_stand, y)
z = transform(mach1, y)
mach2 = machine(stand, X)
W = transform(mach2, X)
mach3 = machine(rgs, W, z)
zhat = predict(mach3, W)
yhat = inverse_transform(mach1, zhat)

@from_network machine(Deterministic(), X, y; predict=yhat) begin
    mutable struct CompositeA
        rgs=rgs
        stand=stand
        target=target_stand
    end
end

X, y = make_regression(20, 2);
model = CompositeA(stand=stand, target=stand)
mach = machine(model, X, y)
@test_logs((:error, r"The hyper"),
           (:error, r"Problem"),
           (:info, r"Running"),
           (:info, r"Type checks okay"),
           @test_throws(ArgumentError,
                        fit!(mach, verbosity=-1)))

end

true


text:
using JDF
using CSV
using DataFrames

#generate data
n=20_000
df0=DataFrame(v=repeat(vcat("abc"),n));
allowmissing(df0)
df0.v = convert(Vector{Union{Missing,String}},df0.v)
df0.v[end] = missing

#write file
csvfile = raw"C:\tmp\afile.csv"
CSV.write(csvfile ,df0)

#read file
csvSep=','
df =  CSV.read(csvfile, DataFrame,threaded=true, delim=csvSep, pool=0.05,strict=true, lazystrings=false);

#save jdf
jdffi = raw"C:\tmp\df.jdf"
jdffile = JDF.savejdf(jdffi, df)

#load jdf
dfloaded = JDF.loadjdf(jdffi)

df.v #<- this one is pooled as expected
dfloaded.v #<- not pooled anymore

df.v

using DataAPI
levels(df.v)
describe(df.v)
DataAPI.refpool(df.v)
refvalue(df.v)
DataAPI.refarray(df.v)
DataAPI.defaultarray(df.v)

Base.summarysize(df)/1024/1024
Base.summarysize(dfloaded)/1024/1024


text:
using Test

const tests = [
    "core",
    "conversions",
    "graph-theory",
    "matrices",
    "operations",
    "properties",
    "queries",
    "traits",
    "utils"
]

for test in tests
    @testset "$test" begin include("$test.jl") end
end


text:
using DashBootstrapComponents, DashHtmlComponents

card = dbc_card(
    [
        dbc_cardimg(src="/static/images/placeholder286x180.png", top=true),
        dbc_cardbody([
            html_h4("Card title", className="card-title"),
            html_p(
                "Some quick example text to build on the card title and " *
                "make up the bulk of the card's content.",
                className="card-text",
            ),
            dbc_button("Go somewhere", color="primary"),
        ]),
    ],
    style=Dict("width" => "18rem"),
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components