text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/other.jl
# Defines isempty, empty, isnan, step, sign, abs, intersect, in.
#############################################################################

@inline function empty(x::MC{N,T}) where {N, T <: RelaxTag}
	MC{N,T}(Inf, -Inf, Interval{Float64}(Inf,-Inf), zeros(SVector{N,Float64}),
            zeros(SVector{N,Float64}), false)
end

interval_MC(x::MC{S,T}) where {S, T<:RelaxTag} = MC{S,T}(x.Intv)

@inline isnan(x::MC) = isnan(x.cc) || isnan(x.cv)

########### Defines differentiable step relaxations
@inline function cv_step(x::Float64, xL::Float64, xU::Float64)
	 xU <= 0.0 && (return 0.0, 0.0)
	 xL >= 0.0 && (return 1.0, 0.0)
	 x >= 0.0 ? ((x/xU)^2, 2.0*x/xU^2) : (0.0, 0.0)
end
@inline function cc_step(x::Float64, xL::Float64, xU::Float64)
	 xU <= 0.0 && (return 0.0, 0.0)
	 xL >= 0.0 && (return 1.0, 0.0)
	 x >= 0.0 ? (1.0, 0.0) : (1.0-(x/xL)^2, -2.0*x/xL^2)
end
@inline function cv_step_NS(x::Float64, xL::Float64, xU::Float64)
	 xU <= 0.0 && (return 0.0, 0.0)
	 xL >= 0.0 && (return 1.0, 0.0)
	 x > 0.0 ? (x/xU, 1.0/xU) : (0.0, 0.0)
end
@inline function cc_step_NS(x::Float64, xL::Float64, xU::Float64)
	 xU <= 0.0 && (return 0.0, 0.0)
	 xL >= 0.0 && (return 1.0, 0.0)
	 x >= 0.0 ? (1.0, 0.0) : ((1.0 - (x/xL)), (-x/xL))
end
@inline function step_kernel(x::MC{N, T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
	 xL = x.Intv.lo
	 xU = x.Intv.hi
	 midcc, cc_id = mid3(x.cc, x.cv, xU)
	 midcv, cv_id = mid3(x.cc, x.cv, xL)
	 cc, dcc = cc_step_NS(midcc, xL, xU)
	 cv, dcv = cv_step_NS(midcv, xL, xU)
	 cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
	 cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	 cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
	 return MC{N,T}(cv, cc, z, cv_grad, cc_grad, x.cnst)
end
@inline function step_kernel(x::MC{N, Diff}, z::Interval{Float64}) where N
	 xL = x.Intv.lo
	 xU = x.Intv.hi
	 midcc = mid3v(x.cc, x.cv, xU)
	 midcv = mid3v(x.cc, x.cv, xL)
	 cc, dcc = cc_step(midcc, xL, xU)
	 cv, dcv = cv_step(midcv, xL, xU)
	 cc1, gdcc1 = cc_step(x.cv, xL, xU)
	 cv1, gdcv1 = cv_step(x.cv, xL, xU)
	 cc2, gdcc2 = cc_step(x.cc, xL, xU)
	 cv2, gdcv2 = cv_step(x.cc, xL, xU)
	 cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
	 cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
	 return MC{N,Diff}(cv, cc, z, cv_grad, cc_grad, x.cnst)
end
@inline step(x::MC) = step_kernel(x, step(x.Intv))

########### Defines sign
@inline function sign_kernel(x::MC{N, T}, z::Interval{Float64}) where {N, T<:RelaxTag}
	zMC = -step(-x) + step(x)
	return MC{N,T}(zMC.cv, zMC.cc, z, zMC.cv_grad, zMC.cc_grad, zMC.cnst)
end
@inline sign(x::MC) = sign_kernel(x, sign(x.Intv))

@inline function cv_abs(x::Float64, xL::Float64, xU::Float64)
	xL >= 0.0 && (return x, 1.0)
	xU <= 0.0 && (return -x, -1.0)
	if x >= 0.0
		return xU*(x/xU)^(MC_DIFF_MU+1), (MC_DIFF_MU+1)*(x/xU)^MC_DIFF_MU
	end
	return -xL*(x/xL)^(MC_DIFF_MU+1), -(MC_DIFF_MU+1)*(x/xL)^MC_DIFF_MU
end
@inline cc_abs(x::Float64,xL::Float64,xU::Float64) = dline_seg(abs, sign, x, xL, xU)
@inline cv_abs_NS(x::Float64,xL::Float64,xU::Float64) = abs(x), sign(x)

@inline function abs_kernel(x::MC{N, T}, z::Interval{Float64}) where {N, T <: Union{NS,MV}}

     xL = x.Intv.lo
	 xU = x.Intv.hi
	 eps_min = mid3v(xL, x.Intv.hi, 0.0)
	 eps_max = (abs(xU) >= abs(xL)) ? xU : xL
	 midcc, cc_id = mid3(x.cc, x.cv, eps_max)
	 midcv, cv_id = mid3(x.cc, x.cv, eps_min)

	 cc, dcc = cc_abs(midcc, xL, xU)
	 cv, dcv = cv_abs_NS(midcv, xL, xU)
	 cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
	 cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	 cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
	 return MC{N,T}(cv, cc, z, cv_grad, cc_grad, x.cnst)
end
@inline function abs_kernel(x::MC{N, Diff}, z::Interval{Float64}) where N

    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_min, blank = mid3(xL, xU, 0.0)
    eps_max = (abs(xU) >= abs(xL)) ? xU : xL
    midcc = mid3v(x.cc, x.cv, eps_max)
    midcv = mid3v(x.cc, x.cv, eps_min)

    cc, dcc = cc_abs(midcc, xL, xU)
    cv, dcv = cv_abs(midcv, xL, xU)
    cc1, gdcc1 = cc_abs(x.cv, xL, xU)
    cv1, gdcv1 = cv_abs(x.cv, xL, xU)
    cc2, gdcc2 = cc_abs(x.cc, xL, xU)
    cv2, gdcv2 = cv_abs(x.cc, xL, xU)

    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N, Diff}(cv, cc, z, cv_grad, cc_grad, x.cnst)
end
@inline abs(x::MC) = abs_kernel(x, abs(x.Intv))

@inline abs2_kernel(x::MC, z::Interval{Float64}) = sqr_kernel(x, z)
@inline abs2(x::MC) = abs2_kernel(x, pow(x.Intv,2))

@inline function correct_intersect(x::MC{N,T}, cv::Float64, cc::Float64, Intv::Interval{Float64}, cv_grad::SVector{N,Float64},
	                               cc_grad::SVector{N,Float64}, cnst::Bool) where {N, T <: RelaxTag}

    if cv - cc < MC_INTERSECT_TOL
		if diam(Intv) > 3.0*MC_INTERSECT_TOL
			cv = max(cv - MC_INTERSECT_TOL, Intv.lo)
			cc = min(cc + MC_INTERSECT_TOL, Intv.hi)
			if cv === Intv.lo
				cv_grad  = zero(SVector{N,Float64})
			end
			if cc === Intv.hi
				cc_grad = zero(SVector{N,Float64})
			end
			return MC{N,T}(cv, cc, Intv, cv_grad, cc_grad, cnst)
		else
			MC_INTERSECT_NOOP_FALLBACK && (return x)
		end
	end
	return MC{N,T}(NaN, NaN, Intv, cv_grad, cc_grad, false)
end

"""
$(TYPEDSIGNATURES)

Intersects two McCormick objects by computing setting `cv = max(cv1, cv2)`
and `cc = min(cc1, cc2)` with appropriate subgradients. Interval are intersected in the
usual fashion. Note that in a typical reverse propagation scheme this may result in
`cv > cc` due to rounding error. This is addressed in the following manner:
- if `cv - cc < MC_INTERSECT_TOL` and `diam(x.Intv ∩ y.Intv) > 3*MC_INTERSECT_TOL`,
then `cv = max(cv - MC_INTERSECT_TOL, Intv.lo)` and `cc = min(cc - MC_INTERSECT_TOL, Intv.hi)`.
Subgradients set to zero if `cv == Intv.lo` or `cc == Intv.hi`.
- if `cv - cc < MC_INTERSECT_TOL` and `diam(x.Intv ∩ y.Intv) < 3*MC_INTERSECT_TOL` and
`MC_INTERSECT_NOOP_FALLBACK == true` then return `x`.
- else return `MC{N,T}(NaN, NaN, x.Intv ∩ y.Intv, ...)`
"""
@inline function intersect(x::MC{N,T}, y::MC{N,T}) where {N, T<:Union{NS, MV}}

     Intv = x.Intv ∩ y.Intv
	 isempty(Intv) && (return empty(x))

	 cnst1 = (x.cc < y.cc) ? x.cnst : y.cnst
	 cnst2 = (x.cv > y.cv) ? x.cnst : y.cnst

     cc = (x.cc < y.cc) ? x.cc : y.cc
     cc_grad = (x.cc < y.cc) ? x.cc_grad : y.cc_grad

     cv = (x.cv > y.cv) ? x.cv : y.cv
     cv_grad = (x.cv > y.cv) ? x.cv_grad : y.cv_grad

	 if cv <= cc
		 return MC{N,T}(cv, cc, Intv, cv_grad, cc_grad, cnst1 && cnst2)
	 end

	 return correct_intersect(x, cv, cc, Intv, cv_grad, cc_grad, cnst1 && cnst2)
end

"""
$(TYPEDSIGNATURES)

Intersects two `MC{N, Diff}` in a manner than preserves differentiability. Interval are intersected in the
usual fashion. Note that in a typical reverse propagation scheme this may result in
`cv > cc` due to rounding error. This is addressed in the following manner:
- if `cv - cc < MC_INTERSECT_TOL` and `diam(x.Intv ∩ y.Intv) > 3*MC_INTERSECT_TOL`,
then `cv = max(cv - MC_INTERSECT_TOL, Intv.lo)` and `cc = min(cc - MC_INTERSECT_TOL, Intv.hi)`.
Subgradients set to zero if `cv == Intv.lo` or `cc == Intv.hi`.
- if `cv - cc < MC_INTERSECT_TOL` and `diam(x.Intv ∩ y.Intv) < 3*MC_INTERSECT_TOL` and
`MC_INTERSECT_NOOP_FALLBACK == true` then return `x`.
- else return `MC{N,T}(NaN, NaN, x.Intv ∩ y.Intv, ...)`
"""
@inline function intersect(x::MC{N, Diff}, y::MC{N, Diff}) where N

	Intv = intersect(x.Intv, y.Intv)
	isempty(Intv) && (return empty(x))

    max_MC = x - max(x - y, 0.0)
    min_MC = y - max(y - x, 0.0)

	if max_MC.cv <= min_MC.cc
    	return MC{N, Diff}(max_MC.cv, min_MC.cc, Intv, max_MC.cv_grad, min_MC.cc_grad, x.cnst && y.cnst)
	end

	return correct_intersect(x, max_MC.cv, min_MC.cc, Intv, max_MC.cv_grad, min_MC.cc_grad, x.cnst && y.cnst)
end

@inline function intersect(x::MC{N, T}, y::Interval{Float64}) where {N, T<:Union{NS,MV}}
	cnst1 = x.cnst
	cnst2 = x.cnst
	if x.cv >= y.lo
  		cv = x.cv
  		cv_grad = x.cv_grad
	else
		cnst1 = true
  		cv = y.lo
  		cv_grad = zero(SVector{N,Float64})
	end
	if x.cc <= y.hi
  		cc = x.cc
  		cc_grad = x.cc_grad
	else
		cnst2 = true
  		cc = y.hi
  		cc_grad = zero(SVector{N,Float64})
	end
	if cv <= cc
		return MC{N, T}(cv, cc, intersect(x.Intv, y), cv_grad, cc_grad, cnst1 && cnst2)
	end
	return correct_intersect(x, cv, cc, intersect(x.Intv, y), cv_grad, cc_grad, cnst1 && cnst2)
end
@inline function intersect(x::MC{N, Diff}, y::Interval{Float64}) where N
     max_MC = x - max(x - y, 0.0)
     min_MC = y - max(y - x, 0.0)
	 if max_MC.cv <= min_MC.cc
		 MC{N, Diff}(max_MC.cv, min_MC.cc, intersect(x.Intv, y), max_MC.cv_grad, min_MC.cc_grad, x.cnst)
	 end
     return correct_intersect(x, max_MC.cv, min_MC.cc, intersect(x.Intv, y), max_MC.cv_grad, min_MC.cc_grad, x.cnst && y.cnst)
end

@inline function intersect(c::Float64, x::MC{N,T}) where {N, T<:RelaxTag}
	isempty(x) && (return empty(x))
	isnan(x) && (return nan(x))

	intv = intersect(x.Intv, Interval{Float64}(c))
	isempty(intv) && (return empty(x))

	cv = max(c, x.cv)
	cc = min(c, x.cc)

	cv_grad = (cv == c) ? zero(SVector{N,Float64}) : x.cv_grad
	cc_grad = (cc == c) ? zero(SVector{N,Float64}) : x.cc_grad
	if cv <= c <= cc
		return MC{N, T}(c, c, intv, cv_grad, cc_grad, true)
	end
	correct_intersect(x, cv, cc, intv, cv_grad, cc_grad, cnst1 && cnst2)
end
@inline intersect(x::MC{N,T}, c::Float64) where {N, T<:RelaxTag} = intersect(c, x)
@inline intersect(c::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = intersect(c, x)

@inline in(a::Int, x::MC) = in(a, x.Intv)
@inline in(a::T, x::MC) where T<:AbstractFloat = in(a, x.Intv)
@inline isempty(x::MC) = isempty(x.Intv)

"""
xlogx

The function `xlogx` is defined as `xlogx(x) = x*log(x)`.
"""
xlogx(x::Float64) = x*log(x)
xlogx_deriv(x::Float64) = log(x) + 1.0
xlogx_deriv2(x::Float64) = 1.0/x
cc_xlogx(x::Float64, xL::Float64, xU::Float64) = dline_seg(xlogx, xlogx_deriv, x, xL, xU)
cv_xlogx(x::Float64, xL::Float64, xU::Float64) = xlogx(x), xlogx_deriv(x)
function xlogx(x::Interval{Float64})
	min_pnt = one(Interval{Float64})/exp(one(Interval{Float64}))
	xlogx_xL = Interval(x.lo)*log(Interval(x.lo))
	xlogx_xU = Interval(x.hi)*log(Interval(x.hi))
	xlogx_min_bool = isdisjoint(min_pnt, x)
	if isdisjoint(min_pnt, x)
		min_val = min(xlogx_xL.lo, xlogx_xU.hi)
	else
		min_val = -min_pnt.hi
	end
	Interval{Float64}(min_val, max(xlogx_xL.lo, xlogx_xU.hi))
end
@inline function xlogx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_max = xlogx(xU) > xlogx(xL) ?  xU : xL
    eps_min = in(1.0/exp(1.0), x) ? 1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ?  xL : xU)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    cc, dcc = cc_xlogx(midcc, xL, xU)
    cv, dcv = cv_xlogx(midcv, xL, xU)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function xlogx_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_max = xlogx(xU) > xlogx(xL) ?  xU : xL
    eps_min = in(1.0/exp(1.0), x) ? 1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ?  xL : xU)
    midcc = mid3v(x.cv, x.cc, eps_max)
    midcv = mid3v(x.cv, x.cc, eps_min)
    cc, dcc = cc_xlogx(midcc, xL, xU)
    cv, dcv = cv_xlogx(midcv, xL, xU)
    gcc1, gdcc1 = cc_xlogx(x.cv, xL, xU)
    gcv1, gdcv1 = cv_xlogx(x.cv, xL, xU)
    gcc2, gdcc2 = cc_xlogx(x.cc, xL, xU)
    gcv2, gdcv2 = cv_xlogx(x.cc, xL, xU)
    cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
    cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
    return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
function xlogx(x::MC{N,T}) where {N, T <: RelaxTag}
	xlogx_kernel(x, xlogx(x.Intv))
end

"""
arh

The arrhenius function `arh` is defined as `arh(x) = exp(-k/x)`.
"""
arh(x, k) = exp(-k/x)
arh(x::Float64, k::MC) = arh_kernel(x, k, x.Intv)
arh_kernel(x::Float64, k::MC, z::Interval{Float64}) = exp(-k/x)
arh(x::Float64, k::Float64) = exp(-k/x)
arh_deriv(x::Float64, k::Float64) = k*exp(-k/x)/x^2
function arh_grad(g, x::Float64, k::Float64)
	g[1] = k*exp(-k/x)/x^2
	g[2] = -exp(-k/x)/x
	nothing
end
@inline function arh_env(x::Float64, y::Float64, k::Float64)
    (x^2)*(exp((-k/y) + (k/x)) - 1.0) - k*(y + x)
end
@inline function arh_envm(x::Float64, y::Float64, k::Float64)
    (x^2)*(exp((-k/y) + (k/x)) - 1.0) + k*(y + x)
end
function cv_arh(x::Float64, xL::Float64, xU::Float64, k::Float64, p::Float64)
	if k > 0.0
	    (xL >= k/2.0) && (return dline_seg(arh, arh_deriv, x, xL, xU, k)..., p)
	    (xU <= k/2.0) && (return arh(x, k), arh_deriv(x, k), p)
	    if p === Inf
	        p, flag = secant(0.0, k/2.0, 0.0, k/2.0, arh_env, xU, k)
	        flag && (p = golden_section(0.0, k/2.0, arh_env, xU, k))
	    end
	    (x <= p) && (return arh(x, k), arh_deriv(x, k), p)
	    return dline_seg(arh, arh_deriv, x, p, xU, k)..., p
	end
	return arh(x, k), arh_deriv(x, k), 0.0
end
function cc_arh(x::Float64, xL::Float64, xU::Float64, k::Float64, p::Float64)
	if k > 0.0
		(xL >= k/2.0) && (return arh(x,k), arh_deriv(x,k), p)
		(xU <= k/2.0) && (return dline_seg(arh, arh_deriv, x, xL, xU, k)..., p)
		if p === Inf
			p, flag = secant(k/2.0, xU, k/2.0, xU, arh_envm, xL, k)
			flag && (p = golden_section(k/2.0, xU, arh_envm, xL, k))
		end
		(x >= p) && (return arh(x,k), arh_deriv(x,k), p)
		return dline_seg(arh, arh_deriv, x, xL, p, k)..., p
	end
	return arh(x,k), arh_deriv(x,k), 0.0
end
function arh_kernel(x::MC{N,T}, k::Float64, z::Interval{Float64},
	                cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
	(k == 0.0) && return one(MC{N,T})
	in(0.0, x) && throw(DomainError(0.0))
	xL = x.Intv.lo
    xU = x.Intv.hi
	eps_min = xL
	eps_max = xU
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    cv, dcv, cv_p = cv_arh(midcv, xL, xU, k, cv_p)
    cc, dcc, cc_p = cc_arh(midcc, xL, xU, k, cc_p)
    cv_grad = mid_grad(x.cv_grad, x.cc_grad, cv_id)*dcv
    cc_grad = mid_grad(x.cv_grad, x.cc_grad, cc_id)*dcc
    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    return MC{N, T}(cv, cc, z, cv_grad, cc_grad, x.cnst), cv_p, cc_p
end
function arh_kernel(x::MC{N,Diff}, k::Float64, z::Interval{Float64},
	                cv_p::Float64, cc_p::Float64) where N
	(k == 0.0) && return one(MC{N,Diff})
	in(0.0, x) && throw(DomainError(0.0))
	xL = x.Intv.lo
	xU = x.Intv.hi
	eps_min = xL
	eps_max = xU
	midcv, cv_id = mid3(x.cv, x.cc, eps_min)
	midcc, cc_id = mid3(x.cv, x.cc, eps_max)
	cv, dcv, cv_p = cv_arh(midcv, xL, xU, k, cv_p)
	cc, dcc, cc_p = cc_arh(midcc, xL, xU, k, cc_p)
	gcv1, gdcv1, cv_p = cv_arh(x.cv, xL, xU, k, cv_p)
	gcc1, gdcc1, cc_p = cc_arh(x.cv, xL, xU, k, cc_p)
	gcv2, gdcv2, cv_p = cv_arh(x.cc, xL, xU, k, cv_p)
	gcc2, gdcc2, cc_p = cc_arh(x.cc, xL, xU, k, cc_p)
	cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
	cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
	return MC{N,Diff}(cv, cc, z, cv_grad, cc_grad, x.cnst), cv_p, cc_p
end
function arh_kernel(x::Float64, k::MC{N,T}, z::Interval{Float64},
	                cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
	exp(-k/x), 0.0, 0.0
end
function arh_kernel(x::Float64, k::MC{N,Diff}, z::Interval{Float64},
	                cv_p::Float64, cc_p::Float64) where N
	exp(-k/x), 0.0, 0.0
end
function arh(x::MC{N,T}, k::Float64) where {N, T <: RelaxTag}
	yMC, tp1, tp2 = arh_kernel(x, k, exp(-k/x.Intv), Inf, Inf)
	return yMC
end

"""
expax

The `expax` function is defined as `expax(x, a) = x*exp(a*x)`.

Form defined in Najman, Jaromił, Dominik Bongartz, and Alexander Mitsos.
"Relaxations of thermodynamic property and costing models in process engineering."
Computers & Chemical Engineering 130 (2019): 106571.
"""
xexpax(x, a) = x*exp(a*x)
xexpax(x::Float64, a::Float64) = x*exp(a*x)
function xexpax(x::Interval{Float64}, a::Float64)
	(a == 0.0) && return x
	xL = x.lo
	xU = x.hi
	fxL = xexpax(Interval(xL), Interval(a)).lo
	fxU = xexpax(Interval(xU), Interval(a)).hi
	zpnt = -1.0/a
	if a > 0.0
		yL = (xL <= zpnt <= xU) ? xexpax(zpnt, a) : ((fxL <= fxU) ? fxL : fxU)
		yU = (fxL >= fxU) ? fxL : fxU
	else
		yL = (fxL <= fxU) ? fxL : fxU
		yU = (xL <= zpnt <= xU) ? xexpax(zpnt, a) : ((fxL >= fxU) ? fxL : fxU)
	end
	Interval(yL, yU)
end
xexpax_deriv(x::Float64, a::Float64) = exp(a*x)*(a*x + 1.0)
function xexpax_grad(g, x::Float64, a::Float64)
	g[1] = exp(a*x)*(a*x + 1.0)
	g[2] = exp(a*x)*x^2
	nothing
end
@inline function xexpax_env(x::Float64, y::Float64, a::Float64)
    (y - x)*xexpax_deriv(x, a) - (xexpax(y, a) - xexpax(x, a))
end
@inline function xexpax_envm(x::Float64, y::Float64, a::Float64)
    (x - y) - (xexpax(x, a) - xexpax(y, a))/xexpax_deriv(x, a)
end
function cv_xexpax(x::Float64, xL::Float64, xU::Float64, a::Float64, p::Float64)
	if a > 0.0
	    (xU <= -2.0/a) && (return dline_seg(xexpax, xexpax_deriv, x, xL, xU, a)..., p)
	    (xL >= -2.0/a) && (return xexpax(x, a), xexpax_deriv(x, a), p)
	    if p === Inf
	        p, flag = secant(-2.0/a, xU, -2.0/a, xU, xexpax_env, xL, a)
	        flag && (p = golden_section(-2.0/a, xU, xexpax_env, xL, a))
	    end
	    (x <= p) && (return dline_seg(xexpax, xexpax_deriv, x, p, xL, a)..., p)
	    return xexpax(x, a), xexpax_deriv(x, a), p
	end
	(xL >= -2.0/a) && (return xexpax(x, a), xexpax_deriv(x, a), p)
	(xU <= -2.0/a) && (return dline_seg(xexpax, xexpax_deriv, x, xL, xU, a)..., p)
	if p === Inf
		p, flag = secant(-2.0/a, xU, -2.0/a, xU, xexpax_env, xL, a)
		flag && (p = golden_section(-2.0/a, xU, xexpax_env, xL, a))
	end
	(x >= p) && (return xexpax(x, a), xexpax_deriv(x, a), p)
	return dline_seg(xexpax, xexpax_deriv, x, xL, p, a)..., p
end
function cc_xexpax(x::Float64, xL::Float64, xU::Float64, a::Float64, p::Float64)
	if a > 0.0
		(xL >= -2.0/a) && (return dline_seg(xexpax, xexpax_deriv, x, xL, xU, a)..., p)
		(xU <= -2.0/a) && (return xexpax(x, a), xexpax_deriv(x, a), p)
		if p === Inf
			p, flag = secant(xL, -2.0/a, xL, -2.0/a, xexpax_env, xU, a)
			flag && (p = golden_section(xL, -2.0/a, xexpax_env, xU, a))
		end
		(x <= p) && (return xexpax(x,a), xexpax_deriv(x,a), p)
		return dline_seg(xexpax, xexpax_deriv, x, p, xU, a)..., p
	end
	(xL >= -2.0/a) && (return dline_seg(xexpax, xexpax_deriv, x, xL, xU, a)..., p)
	(xU <= -2.0/a) && (return xexpax(x,a), xexpax_deriv(x,a), p)
	if p === Inf
		p, flag = secant(xL, -2.0/a, xL, -2.0/a, xexpax_envm, xU, a)
		flag && (p = golden_section(xL, a/2.0, xexpax_envm, xU, a))
	end
	(x <= p) && (return xexpax(x, a), xexpax_deriv(x, a), p)
	return dline_seg(xexpax, xexpax_deriv, x, p, xU, a)..., p
end
function xexpax_kernel(x::Float64, a::MC{N,T}, z::Interval{Float64},
	                   cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
	x*exp(a*x), 0.0, 0.0
end
function xexpax_kernel(x::Float64, a::MC{N,Diff}, z::Interval{Float64},
	                   cv_p::Float64, cc_p::Float64) where N
	x*exp(a*x), 0.0, 0.0
end
function xexpax_kernel(x::MC{N,T}, a::Float64, z::Interval{Float64},
	                   cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
	(a == 0.0) && return one(MC{N,T})
	xL = x.Intv.lo
    xU = x.Intv.hi
	zpnt = -1.0/a
	fxL = xexpax(xL, a)
	fxU = xexpax(xU, a)
	if a > 0.0
		eps_min = (xL <= zpnt <= xU) ? zpnt : ((fxL <= fxU) ? xL : xU)
		eps_max = (fxL >= fxU) ? xL : xU
	else
		eps_min = (fxL <= fxU) ? xL : xU
		eps_max = (xL <= zpnt <= xU) ? zpnt : ((fxL >= fxU) ? xL : xU)
	end
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    cv, dcv, cv_p = cv_xexpax(midcv, xL, xU, a, cv_p)
    cc, dcc, cc_p = cc_xexpax(midcc, xL, xU, a, cc_p)
    cv_grad = mid_grad(x.cv_grad, x.cc_grad, cv_id)*dcv
    cc_grad = mid_grad(x.cv_grad, x.cc_grad, cc_id)*dcc
    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    return MC{N, T}(cv, cc, z, cv_grad, cc_grad, x.cnst), cv_p, cc_p
end
function xexpax_kernel(x::MC{N,Diff}, a::Float64, z::Interval{Float64},
	                   cv_p::Float64, cc_p::Float64) where N
	(a == 0.0) && return one(MC{N,Diff})
	xL = x.Intv.lo
	xU = x.Intv.hi
	zpnt = -1.0/a
	fxL = xexpax(xL, a)
	fxU = xexpax(xU, a)
	if a > 0.0
		eps_min = (xL <= zpnt <= xU) ? zpnt : ((fxL <= fxU) ? xL : xU)
		eps_max = (fxL >= fxU) ? xL : xU
	else
		eps_min = (fxL <= fxU) ? xL : xU
		eps_max = (xL <= zpnt <= xU) ? zpnt : ((fxL >= fxU) ? xL : xU)
	end
	midcv, cv_id = mid3(x.cv, x.cc, eps_min)
	midcc, cc_id = mid3(x.cv, x.cc, eps_max)
	cv, dcv, cv_p = cv_xexpax(midcv, xL, xU, a, cv_p)
	cc, dcc, cc_p = cc_xexpax(midcc, xL, xU, a, cc_p)
	gcv1, gdcv1, cv_p = cv_xexpax(x.cv, xL, xU, a, cv_p)
	gcc1, gdcc1, cc_p = cc_xexpax(x.cv, xL, xU, a, cc_p)
	gcv2, gdcv2, cv_p = cv_xexpax(x.cc, xL, xU, a, cv_p)
	gcc2, gdcc2, cc_p = cc_xexpax(x.cc, xL, xU, a, cc_p)
	cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad
	cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad
	return MC{N,Diff}(cv, cc, z, cv_grad, cc_grad, x.cnst), cv_p, cc_p
end
xexpax(x::Float64, a::MC) = x*exp(a*x)
function xexpax(x::MC{N,T}, a::Float64) where {N, T <: RelaxTag}
	intv_xexpax = xexpax(x.Intv, a)
	yMC, tp1, tp2 = xexpax_kernel(x, a, intv_xexpax, Inf, Inf)
	return yMC
end

"""
xabsx

The function `xabsx` is defined as `xabsx(x) = x*abs(x)`.
"""
xabsx(x::Float64) = x*abs(x)
xabsx_deriv(x::Float64) = 2*abs(x)
xabsx_deriv2(x::Float64) = 2*abs(x)/x
@inline function cc_xabsx(x::Float64, xL::Float64, xU::Float64)
	(xU <= 0.0) && (return xabsx(x), xabsx_deriv(x))
	(0.0 <= xL) && (return dline_seg(xabsx, xabsx_deriv, x, xL, xU))

	root_f(x, xL, xU) = (xabsx(xU) - xabsx(x)) - 2*(xU - x)*abs(x)
	root_df(x, xL, xU) = xabsx_deriv(x) - xU*xabsx_deriv2(x)
	inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU)
	flag && (inflection = golden_section(xL, xU, root_f, xL, xU))
	if x <= inflection
		return xabsx(x), xabsx_deriv(x)
	else
		return dline_seg(xabsx, xabsx_deriv, x, inflection, xU)
	end
end
@inline function cv_xabsx(x::Float64, xL::Float64, xU::Float64)
	(xU <= 0.0) && (return dline_seg(xabsx, xabsx_deriv, x, xL, xU))
	(0.0 <= xL) && (return xabsx(x), xabsx_deriv(x))

	root_f(x, xL, xU) = (xabsx(x) - xabsx(xL)) - 2*(x - xL)*abs(x)
	root_df(x, xL, xU) = -xabsx_deriv(x) + xL*xabsx_deriv2(x)
	inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU)
	flag && (inflection = golden_section(xL, xU, root_f, xL, xU))
	if x <= inflection
		return dline_seg(xabsx, xabsx_deriv, x, xL, inflection)
	else
		return xabsx(x), xabsx_deriv(x)
	end
end
function xabsx(x::Interval{Float64})
	xabsx_xL = Interval(x.lo)*abs(Interval(x.lo))
	xabsx_xU = Interval(x.hi)*abs(Interval(x.hi))
	Interval{Float64}(xabsx_xL.lo, xabsx_xU.hi)
end
@inline function xabsx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
	xL = x.Intv.lo
	xU = x.Intv.hi
	eps_max = xU
	eps_min = xL
	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
	cc, dcc = cc_xabsx(midcc, xL, xU)
	cv, dcv = cv_xabsx(midcv, xL, xU)
	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function xabsx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Diff}
	xL = x.Intv.lo
	xU = x.Intv.hi
	eps_min = xL
	eps_max = xU
	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
	cc, dcc = cc_xabsx(midcc, xL, xU)
	cv, dcv = cv_xabsx(midcv, xL, xU)
	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
function xabsx(x::MC{N,T}) where {N, T <: RelaxTag}
	xabsx_kernel(x, xabsx(x.Intv))
end


text:
using MechGlueRecursiveArrayTools
using Test

@testset "MechGlueRecursiveArrayTools.jl" begin
    # Write your tests here.
end

ext: jl


text:
using FileIO, ImageCore, Unitful, ImageAxes, ImageMetadata
using ImageAxes.AxisArrays
using Test, Base.CoreLogging

struct UnsupportedInt16 <: Signed
    val::Int16
end
Base.flipsign(x::UnsupportedInt16, y::UnsupportedInt16) = UnsupportedInt16(flipsign(x.val, y.val))
Base.Unsigned(x::UnsupportedInt16) = Unsigned(x.val)
Base.UInt16(x::UnsupportedInt16) = UInt16(x.val)
Base.Int(x::UnsupportedInt16) = Int(x.val)
Base.promote_rule(::Type{UnsupportedInt16}, ::Type{Int}) = Int

include("unu-make.jl")

@testset "NRRD" begin
    workdir = joinpath(tempdir(), "Images")
    writedir = joinpath(workdir, "write")
    if !isdir(workdir)
        mkdir(workdir)
    end
    if !isdir(writedir)
        mkdir(writedir)
    end

    @testset "Gray, raw" begin
        img = load(joinpath(dirname(@__FILE__), "io", "small.nrrd"))
        flush(stdout)
        @test eltype(img) == Float32
        @test ndims(img) == 3
        @test size(img) == (10,10,10)
        outname = joinpath(writedir, "small.nrrd")
        isfile(outname) && rm(outname)
        save(outname, img)
        imgc = load(outname)
        @test img == imgc
        @test axisnames(imgc) == axisnames(img)
        @test axisvalues(imgc) == axisvalues(img)
        # Check that FixedPoint types get properly encoded
        outname = joinpath(writedir, "small14.nrrd")
        img = rand(Gray{N2f14}, 5, 5, 2)
        save(outname, img)
        imgr = load(outname)
        @test eltype(imgr) == Gray{N2f14}
        @test imgr == img
    end

    @testset "Units" begin
        img = load(joinpath(dirname(@__FILE__), "io", "units.nhdr"))
        axsv = axisvalues(img)
        @test step(axsv[1]) == 0.1u"mm"
        @test step(axsv[2]) == 0.2u"mm"
        @test step(axsv[3]) == 1u"mm"
        axsn = axisnames(img)
        @test axsn == (:L, :P, :S)
        outname = joinpath(writedir, "units.nrrd")
        isfile(outname) && rm(outname)
        save(outname, img)
        imgc = load(outname)
        @test img == imgc
        axsv = axisvalues(imgc)
        @test step(axsv[1]) == 0.1u"mm"
        @test step(axsv[2]) == 0.2u"mm"
        @test step(axsv[3]) == 1u"mm"
        axsn = axisnames(imgc)
        @test axsn == (:L, :P, :S)
    end

    @testset "Compressed (gzip)" begin
        img = load(joinpath(dirname(@__FILE__), "io", "smallgz.nrrd"))
        @test ndims(img) == 3
        @test eltype(img) == Float32
        outname = joinpath(writedir, "smallgz.nrrd")
        isfile(outname) && rm(outname)
        save(outname, img)
        imgc = load(outname)
        @test img == imgc
    end

    @testset "Time is 4th dimension" begin
        img = load(joinpath(dirname(@__FILE__), "io", "small_time.nrrd"))
        @test timedim(img) == 4
        outname = joinpath(writedir, "small_time.nrrd")
        isfile(outname) && rm(outname)
        save(outname, img)
        imgc = load(outname)
        @test img == imgc
        @test axisnames(img)[4] == axisnames(imgc)[4] == :time
        @test axisvalues(img) == axisvalues(imgc)
        img = AxisArray(rand(N0f8, 6, 5, 3, 2), :x, :y, :z, :time)
        save(outname, img)
        imgr = load(outname)
        @test axisnames(imgr)[4] == :time
    end

    @testset "eltype" begin
        for (elty, name) in ((UInt8, "uint8"),
                             (UInt8, "u8"),
                             (Gray{N0f8}, "gray_u8"),
                             (RGB{N0f8}, "rgb_u8"),
                             (Gray{N0f16}, "gray_u16"),
                             (Gray{N4f12}, "gray_u12"),
                             (Gray{N4f12}, "gray_u12_hex"))
            fn = joinpath(dirname(@__FILE__), "io", "eltype_$name.nrrd")
            img = load(fn)
            @test eltype(img) == elty
            @test size(img) == (3,5)
            outname = joinpath(writedir, "elt.nrrd")
            save(outname, img)
            imgr = load(outname)
            @test eltype(imgr) == elty
        end
    end

    @testset "bool" begin
        # Bool is not one of the supported eltypes in the `type` field of
        # http://teem.sourceforge.net/nrrd/format.html.
        # Make sure we don't do something wrong
        fn = joinpath(writedir, "eltype_bool.nrrd")
        img = rand(Bool, 5, 6)
        save(fn, img)
        imgr = load(fn)
        @test imgr == img
        # We don't insist on eltype(imgr) == Bool
        img = UnsupportedInt16[1 2; 3 4]
        @info "The following error is expected"
        fn = joinpath(writedir, "eltype_unsupported.nrrd")
        ex = try save(fn, img) catch e e end
        @test ex.ex == ErrorException("integer type UnsupportedInt16 not supported")
    end

    @testset "Mmapped" begin
        fn = joinpath(writedir, "volume.nrrd")
        save(fn, zeros(UInt8,500,500,500))
        v = load(fn)
        @test size(v) == (500,500,500)
        @test all(x->x==0, v)
    end

    @testset "Endian mmapped" begin
        fn = joinpath(dirname(@__FILE__), "io", "bswap.nrrd")
        img = load(fn; mmap=true)
        @test img == [255 254; 1 2]
    end

    @testset "Fiji compatibility" begin
        for (name, ps) in (("fiji-stack-pixels-16bit.nrrd", (1,1,1)),
                           ("fiji-stack-microns-16bit.nrrd", (0.4μm, 0.4μm, 2.0μm)))
            fn = joinpath(dirname(@__FILE__), "io", name)
            v = load(fn)
            @test size(v) == (50,40,30)
            @test pixelspacing(v) == ps
        end
    end

    @testset "Header only" begin
        img = load(joinpath(dirname(@__FILE__), "io", "small.nrrd"))
        outname = joinpath(writedir, "small.nhdr")
        isfile(outname) && rm(outname)
        props = Dict("datafile"=>joinpath(writedir,"small.raw"))
        save(outname, img; props=props)
        imgc = load(outname)
        @test img == imgc
        @test axisnames(imgc) == axisnames(img)
        @test axisvalues(imgc) == axisvalues(img)
    end

    @testset "N0f16 ImageMeta" begin
        img = ImageMeta(rand(N0f16, 5, 7), info=false)
        outname = joinpath(writedir, "imagemeta.nhdr")
        save(outname, img)
        imgc = load(outname)
        @test img == imgc

        # With units
        img = AxisArray(img, (:y, :x), (1.0u"mm", 1.2u"mm"))
        save(outname, img)
        imgc = load(outname)
        @test img == imgc
        @test pixelspacing(imgc) == (1.0u"mm", 1.2u"mm")
        img = ImageMeta(img, info=false)
        save(outname, img)
        imgc = load(outname)
        @test img == imgc
        @test AxisArrays.HasAxes(img) isa AxisArrays.HasAxes{true}
        @test pixelspacing(imgc) == (1.0u"mm", 1.2u"mm")
    end

    @testset "Size warning" begin
        for (T, Traw) in ((N0f8, UInt8), (N0f16, UInt16))
            img = rand(T, 5, 7, 3)
            outname = joinpath(writedir, "imagetrunc.nhdr")
            props = Dict("datafile"=>joinpath(writedir, "imagetrunc.raw"))
            save(outname, img, props=props)
            # Create a too-small data array
            open(props["datafile"], "w") do io
                write(io, rand(Traw, 5, 7, 2))
            end
            # Test that it can be read but produces a warning
            logger = Test.TestLogger()
            with_logger(logger) do
                imgr = load(outname)
                @test size(imgr) == (5, 7, 2)
            end
            record = logger.logs[1]
            @test occursin("header indicates an array size (5, 7, 3), but the file size is consistent with at most (5, 7, 2)", record.message)
        end
    end

    @testset "Extended eltype" begin
        img = load(joinpath(dirname(@__FILE__), "io", "test_helix.nrrd"))
        @test axisnames(img)[end-2:end] == (:R, :A, :S)
    end

    GC.gc()  # to close any mmapped files
    try
        rm(workdir, recursive=true)
    catch
    end
end

include("readremote.jl")

nothing


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module HigherOrderFns

# This module provides higher order functions specialized for sparse arrays,
# particularly map[!]/broadcast[!] for SparseVectors and SparseMatrixCSCs at present.
import Base: map, map!, broadcast, copy, copyto!

using Base: front, tail, to_shape
using ..SparseArrays: SparseVector, SparseMatrixCSC, AbstractSparseVector, AbstractSparseMatrixCSC,
                      AbstractSparseMatrix, AbstractSparseArray, indtype, nnz, nzrange,
                      SparseVectorUnion, AdjOrTransSparseVectorUnion, nonzeroinds, nonzeros, rowvals, getcolptr
using Base.Broadcast: BroadcastStyle, Broadcasted, flatten
using LinearAlgebra

# This module is organized as follows:
# (0) Define BroadcastStyle rules and convenience types for dispatch
# (1) Define a common interface to SparseVectors and SparseMatrixCSCs sufficient for
#       map[!]/broadcast[!]'s purposes. The methods below are written against this interface.
# (2) Define entry points for map[!] (short children of _map_[not]zeropres!).
# (3) Define entry points for broadcast[!] (short children of _broadcast_[not]zeropres!).
# (4) Define _map_[not]zeropres! specialized for a single (input) sparse vector/matrix.
# (5) Define _map_[not]zeropres! specialized for a pair of (input) sparse vectors/matrices.
# (6) Define general _map_[not]zeropres! capable of handling >2 (input) sparse vectors/matrices.
# (7) Define _broadcast_[not]zeropres! specialized for a single (input) sparse vector/matrix.
# (8) Define _broadcast_[not]zeropres! specialized for a pair of (input) sparse vectors/matrices.
# (9) Define general _broadcast_[not]zeropres! capable of handling >2 (input) sparse vectors/matrices.
# (10) Define broadcast methods handling combinations of broadcast scalars and sparse vectors/matrices.
# (11) Define broadcast[!] methods handling combinations of scalars, sparse vectors/matrices,
#       structured matrices, and one- and two-dimensional Arrays.
# (12) Define map[!] methods handling combinations of sparse and structured matrices.

# (0) BroadcastStyle rules and convenience types for dispatch

SparseVecOrMat = Union{SparseVector,AbstractSparseMatrixCSC}

# broadcast container type promotion for combinations of sparse arrays and other types
struct SparseVecStyle <: Broadcast.AbstractArrayStyle{1} end
struct SparseMatStyle <: Broadcast.AbstractArrayStyle{2} end
Broadcast.BroadcastStyle(::Type{<:SparseVector}) = SparseVecStyle()
Broadcast.BroadcastStyle(::Type{<:AbstractSparseMatrixCSC}) = SparseMatStyle()
const SPVM = Union{SparseVecStyle,SparseMatStyle}

# SparseVecStyle handles 0-1 dimensions, SparseMatStyle 0-2 dimensions.
# SparseVecStyle promotes to SparseMatStyle for 2 dimensions.
# Fall back to DefaultArrayStyle for higher dimensionality.
SparseVecStyle(::Val{0}) = SparseVecStyle()
SparseVecStyle(::Val{1}) = SparseVecStyle()
SparseVecStyle(::Val{2}) = SparseMatStyle()
SparseVecStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()
SparseMatStyle(::Val{0}) = SparseMatStyle()
SparseMatStyle(::Val{1}) = SparseMatStyle()
SparseMatStyle(::Val{2}) = SparseMatStyle()
SparseMatStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()

Broadcast.BroadcastStyle(::SparseMatStyle, ::SparseVecStyle) = SparseMatStyle()

# Tuples promote to dense
Broadcast.BroadcastStyle(::SparseVecStyle, ::Broadcast.Style{Tuple}) = Broadcast.DefaultArrayStyle{1}()
Broadcast.BroadcastStyle(::SparseMatStyle, ::Broadcast.Style{Tuple}) = Broadcast.DefaultArrayStyle{2}()

struct PromoteToSparse <: Broadcast.AbstractArrayStyle{2} end
PromoteToSparse(::Val{0}) = PromoteToSparse()
PromoteToSparse(::Val{1}) = PromoteToSparse()
PromoteToSparse(::Val{2}) = PromoteToSparse()
PromoteToSparse(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()

const StructuredMatrix = Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}
Broadcast.BroadcastStyle(::Type{<:Adjoint{T,<:Union{SparseVector,SparseMatrixCSC}} where T}) = PromoteToSparse()
Broadcast.BroadcastStyle(::Type{<:Transpose{T,<:Union{SparseVector,SparseMatrixCSC}} where T}) = PromoteToSparse()

Broadcast.BroadcastStyle(s::SPVM, ::Broadcast.AbstractArrayStyle{0}) = s
Broadcast.BroadcastStyle(s::SPVM, ::Broadcast.DefaultArrayStyle{0}) = s
Broadcast.BroadcastStyle(::SPVM, ::Broadcast.DefaultArrayStyle{1}) = PromoteToSparse()
Broadcast.BroadcastStyle(::SPVM, ::Broadcast.DefaultArrayStyle{2}) = PromoteToSparse()

Broadcast.BroadcastStyle(::SPVM, ::LinearAlgebra.StructuredMatrixStyle{<:StructuredMatrix}) = PromoteToSparse()
Broadcast.BroadcastStyle(::PromoteToSparse, ::LinearAlgebra.StructuredMatrixStyle{<:StructuredMatrix}) = PromoteToSparse()

Broadcast.BroadcastStyle(::PromoteToSparse, ::SPVM) = PromoteToSparse()
Broadcast.BroadcastStyle(::PromoteToSparse, ::Broadcast.Style{Tuple}) = Broadcast.DefaultArrayStyle{2}()

# FIXME: currently sparse broadcasts are only well-tested on known array types, while any AbstractArray
# could report itself as a DefaultArrayStyle().
# See https://github.com/JuliaLang/julia/pull/23939#pullrequestreview-72075382 for more details
is_supported_sparse_broadcast() = true
is_supported_sparse_broadcast(::AbstractArray, rest...) = false
is_supported_sparse_broadcast(::AbstractSparseArray, rest...) = is_supported_sparse_broadcast(rest...)
is_supported_sparse_broadcast(::StructuredMatrix, rest...) = is_supported_sparse_broadcast(rest...)
is_supported_sparse_broadcast(::Array, rest...) = is_supported_sparse_broadcast(rest...)
is_supported_sparse_broadcast(t::Union{Transpose, Adjoint}, rest...) = is_supported_sparse_broadcast(t.parent, rest...)
is_supported_sparse_broadcast(x, rest...) = axes(x) === () && is_supported_sparse_broadcast(rest...)
is_supported_sparse_broadcast(x::Ref, rest...) = is_supported_sparse_broadcast(rest...)

can_skip_sparsification(f, rest...) = false
can_skip_sparsification(::typeof(*), ::SparseVectorUnion, ::AdjOrTransSparseVectorUnion) = true

# Dispatch on broadcast operations by number of arguments
const Broadcasted0{Style<:Union{Nothing,BroadcastStyle},Axes,F} =
    Broadcasted{Style,Axes,F,Tuple{}}
const SpBroadcasted1{Style<:SPVM,Axes,F,Args<:Tuple{SparseVecOrMat}} =
    Broadcasted{Style,Axes,F,Args}
const SpBroadcasted2{Style<:SPVM,Axes,F,Args<:Tuple{SparseVecOrMat,SparseVecOrMat}} =
    Broadcasted{Style,Axes,F,Args}

# (1) The definitions below provide a common interface to sparse vectors and matrices
# sufficient for the purposes of map[!]/broadcast[!]. This interface treats sparse vectors
# as n-by-one sparse matrices which, though technically incorrect, is how broacast[!] views
# sparse vectors in practice.
@inline numrows(A::SparseVector) = length(A)
@inline numrows(A::AbstractSparseMatrixCSC) = size(A, 1)
@inline numcols(A::SparseVector) = 1
@inline numcols(A::AbstractSparseMatrixCSC) = size(A, 2)
# numrows and numcols respectively yield size(A, 1) and size(A, 2), but avoid a branch
@inline columns(A::SparseVector) = 1
@inline columns(A::AbstractSparseMatrixCSC) = 1:size(A, 2)
@inline colrange(A::SparseVector, j) = 1:length(nonzeroinds(A))
@inline colrange(A::AbstractSparseMatrixCSC, j) = nzrange(A, j)
@inline colstartind(A::SparseVector, j) = one(indtype(A))
@inline colboundind(A::SparseVector, j) = convert(indtype(A), length(nonzeroinds(A)) + 1)
@inline colstartind(A::AbstractSparseMatrixCSC, j) = getcolptr(A)[j]
@inline colboundind(A::AbstractSparseMatrixCSC, j) = getcolptr(A)[j + 1]
@inline storedinds(A::SparseVector) = nonzeroinds(A)
@inline storedinds(A::AbstractSparseMatrixCSC) = rowvals(A)
@inline storedvals(A::SparseVecOrMat) = nonzeros(A)
@inline setcolptr!(A::SparseVector, j, val) = val
@inline setcolptr!(A::AbstractSparseMatrixCSC, j, val) = getcolptr(A)[j] = val
function trimstorage!(A::SparseVecOrMat, maxstored)
    resize!(storedinds(A), maxstored)
    resize!(storedvals(A), maxstored)
    return maxstored
end
function expandstorage!(A::SparseVecOrMat, maxstored)
    length(storedinds(A)) < maxstored && resize!(storedinds(A), maxstored)
    length(storedvals(A)) < maxstored && resize!(storedvals(A), maxstored)
    return maxstored
end

# (2) map[!] entry points
map(f::Tf, A::SparseVector) where {Tf} = _noshapecheck_map(f, A)
map(f::Tf, A::AbstractSparseMatrixCSC) where {Tf} = _noshapecheck_map(f, A)
map(f::Tf, A::AbstractSparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N}) where {Tf,N} =
    (_checksameshape(A, Bs...); _noshapecheck_map(f, A, Bs...))
map(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N} =
    (_checksameshape(A, Bs...); _noshapecheck_map(f, A, Bs...))
map!(f::Tf, C::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC, Bs::Vararg{SparseMatrixCSC,N}) where {Tf,N} =
    (_checksameshape(C, A, Bs...); _noshapecheck_map!(f, C, A, Bs...))
map!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N} =
    (_checksameshape(C, A, Bs...); _noshapecheck_map!(f, C, A, Bs...))
function _noshapecheck_map!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N}
    fofzeros = f(_zeros_eltypes(A, Bs...)...)
    fpreszeros = _iszero(fofzeros)
    return fpreszeros ? _map_zeropres!(f, C, A, Bs...) :
                        _map_notzeropres!(f, fofzeros, C, A, Bs...)
end
function _noshapecheck_map(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N}
    fofzeros = f(_zeros_eltypes(A, Bs...)...)
    fpreszeros = _iszero(fofzeros)
    maxnnzC = fpreszeros ? min(length(A), _sumnnzs(A, Bs...)) : length(A)
    entrytypeC = Base.Broadcast.combine_eltypes(f, (A, Bs...))
    indextypeC = _promote_indtype(A, Bs...)
    C = _allocres(size(A), indextypeC, entrytypeC, maxnnzC)
    return fpreszeros ? _map_zeropres!(f, C, A, Bs...) :
                        _map_notzeropres!(f, fofzeros, C, A, Bs...)
end
# (3) broadcast[!] entry points
copy(bc::SpBroadcasted1) = _noshapecheck_map(bc.f, bc.args[1])

@inline function copyto!(C::SparseVecOrMat, bc::Broadcasted0{Nothing})
    isempty(C) && return _finishempty!(C)
    f = bc.f
    fofnoargs = f()
    if _iszero(fofnoargs) # f() is zero, so empty C
        trimstorage!(C, 0)
        _finishempty!(C)
    else # f() is nonzero, so densify C and fill with independent calls to f()
        _densestructure!(C)
        storedvals(C)[1] = fofnoargs
        broadcast!(f, view(storedvals(C), 2:length(storedvals(C))))
    end
    return C
end

function _diffshape_broadcast(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N}
    fofzeros = f(_zeros_eltypes(A, Bs...)...)
    fpreszeros = _iszero(fofzeros)
    indextypeC = _promote_indtype(A, Bs...)
    entrytypeC = Base.Broadcast.combine_eltypes(f, (A, Bs...))
    shapeC = to_shape(Base.Broadcast.combine_axes(A, Bs...))
    maxnnzC = fpreszeros ? _checked_maxnnzbcres(shapeC, A, Bs...) : _densennz(shapeC)
    C = _allocres(shapeC, indextypeC, entrytypeC, maxnnzC)
    return fpreszeros ? _broadcast_zeropres!(f, C, A, Bs...) :
                        _broadcast_notzeropres!(f, fofzeros, C, A, Bs...)
end
# helper functions for map[!]/broadcast[!] entry points (and related methods below)
@inline _sumnnzs(A) = nnz(A)
@inline _sumnnzs(A, Bs...) = nnz(A) + _sumnnzs(Bs...)
@inline _iszero(x) = x == 0
@inline _iszero(x::Number) = Base.iszero(x)
@inline _iszero(x::AbstractArray) = Base.iszero(x)
@inline _zeros_eltypes(A) = (zero(eltype(A)),)
@inline _zeros_eltypes(A, Bs...) = (zero(eltype(A)), _zeros_eltypes(Bs...)...)
@inline _promote_indtype(A) = indtype(A)
@inline _promote_indtype(A, Bs...) = promote_type(indtype(A), _promote_indtype(Bs...))
@inline _aresameshape(A) = true
@inline _aresameshape(A, B) = size(A) == size(B)
@inline _aresameshape(A, B, Cs...) = _aresameshape(A, B) ? _aresameshape(B, Cs...) : false
@inline _checksameshape(As...) = _aresameshape(As...) || throw(DimensionMismatch("argument shapes must match"))
@inline _all_args_isa(t::Tuple{Any}, ::Type{T}) where T = isa(t[1], T)
@inline _all_args_isa(t::Tuple{Any,Vararg{Any}}, ::Type{T}) where T = isa(t[1], T) & _all_args_isa(tail(t), T)
@inline _all_args_isa(t::Tuple{Broadcasted}, ::Type{T}) where T = _all_args_isa(t[1].args, T)
@inline _all_args_isa(t::Tuple{Broadcasted,Vararg{Any}}, ::Type{T}) where T = _all_args_isa(t[1].args, T) & _all_args_isa(tail(t), T)
@inline _densennz(shape::NTuple{1}) = shape[1]
@inline _densennz(shape::NTuple{2}) = shape[1] * shape[2]
_maxnnzfrom(shape::NTuple{1}, A::SparseVector) = nnz(A) * div(shape[1], length(A))
_maxnnzfrom(shape::NTuple{2}, A::SparseVector) = nnz(A) * div(shape[1], length(A)) * shape[2]
_maxnnzfrom(shape::NTuple{2}, A::AbstractSparseMatrixCSC) = nnz(A) * div(shape[1], size(A, 1)) * div(shape[2], size(A, 2))
@inline _maxnnzfrom_each(shape, ::Tuple{}) = ()
@inline _maxnnzfrom_each(shape, As) = (_maxnnzfrom(shape, first(As)), _maxnnzfrom_each(shape, tail(As))...)
@inline _unchecked_maxnnzbcres(shape, As::Tuple) = min(_densennz(shape), sum(_maxnnzfrom_each(shape, As)))
@inline _unchecked_maxnnzbcres(shape, As...) = _unchecked_maxnnzbcres(shape, As)
@inline _checked_maxnnzbcres(shape::NTuple{1}, As...) = shape[1] != 0 ? _unchecked_maxnnzbcres(shape, As) : 0
@inline _checked_maxnnzbcres(shape::NTuple{2}, As...) = shape[1] != 0 && shape[2] != 0 ? _unchecked_maxnnzbcres(shape, As) : 0
@inline function _allocres(shape::NTuple{1}, indextype, entrytype, maxnnz)
    storedinds = Vector{indextype}(undef, maxnnz)
    storedvals = Vector{entrytype}(undef, maxnnz)
    return SparseVector(shape..., storedinds, storedvals)
end
@inline function _allocres(shape::NTuple{2}, indextype, entrytype, maxnnz)
    pointers = ones(indextype, shape[2] + 1)
    storedinds = Vector{indextype}(undef, maxnnz)
    storedvals = Vector{entrytype}(undef, maxnnz)
    return SparseMatrixCSC(shape..., pointers, storedinds, storedvals)
end

# (4) _map_zeropres!/_map_notzeropres! specialized for a single sparse vector/matrix
"Stores only the nonzero entries of `map(f, Array(A))` in `C`."
function _map_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat) where Tf
    spaceC::Int = min(length(storedinds(C)), length(storedvals(C)))
    Ck = 1
    @inbounds for j in columns(C)
        setcolptr!(C, j, Ck)
        for Ak in colrange(A, j)
            Cx = f(storedvals(A)[Ak])
            if !_iszero(Cx)
                Ck > spaceC && (spaceC = expandstorage!(C, Ck + nnz(A) - (Ak - 1)))
                storedinds(C)[Ck] = storedinds(A)[Ak]
                storedvals(C)[Ck] = Cx
                Ck += 1
            end
        end
    end
    @inbounds setcolptr!(C, numcols(C) + 1, Ck)
    trimstorage!(C, Ck - 1)
    return C
end
"""
Densifies `C`, storing `fillvalue` in place of each unstored entry in `A` and
`f(A[i])`/`f(A[i,j])` in place of each stored entry `A[i]`/`A[i,j]` in `A`.
"""
function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat) where Tf
    # Build dense matrix structure in C, expanding storage if necessary
    _densestructure!(C)
    # Populate values
    fill!(storedvals(C), fillvalue)
    @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
        for Ak in colrange(A, j)
            Cx = f(storedvals(A)[Ak])
            Cx != fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx)
        end
    end
    # NOTE: Combining the fill! above into the loop above to avoid multiple sweeps over /
    # nonsequential access of storedvals(C) does not appear to improve performance.
    return C
end
# helper functions for these methods and some of those below
@inline _densecoloffsets(A::SparseVector) = 0
@inline _densecoloffsets(A::AbstractSparseMatrixCSC) = 0:size(A, 1):(size(A, 1)*(size(A, 2) - 1))
function _densestructure!(A::SparseVector)
    expandstorage!(A, length(A))
    copyto!(nonzeroinds(A), 1:length(A))
    return A
end
function _densestructure!(A::AbstractSparseMatrixCSC)
    nnzA = size(A, 1) * size(A, 2)
    expandstorage!(A, nnzA)
    copyto!(getcolptr(A), 1:size(A, 1):(nnzA + 1))
    for k in _densecoloffsets(A)
        copyto!(rowvals(A), k + 1, 1:size(A, 1))
    end
    return A
end

# (5) _map_zeropres!/_map_notzeropres! specialized for a pair of sparse vectors/matrices
function _map_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf
    spaceC::Int = min(length(storedinds(C)), length(storedvals(C)))
    rowsentinelA = convert(indtype(A), numrows(C) + 1)
    rowsentinelB = convert(indtype(B), numrows(C) + 1)
    Ck = 1
    @inbounds for j in columns(C)
        setcolptr!(C, j, Ck)
        Ak, stopAk = colstartind(A, j), colboundind(A, j)
        Bk, stopBk = colstartind(B, j), colboundind(B, j)
        Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
        Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
        while true
            if Ai == Bi
                Ai == rowsentinelA && break # column complete
                Cx, Ci::indtype(C) = f(storedvals(A)[Ak], storedvals(B)[Bk]), Ai
                Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
            elseif Ai < Bi
                Cx, Ci = f(storedvals(A)[Ak], zero(eltype(B))), Ai
                Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
            else # Bi < Ai
                Cx, Ci = f(zero(eltype(A)), storedvals(B)[Bk]), Bi
                Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
            end
            # NOTE: The ordering of the conditional chain above impacts which matrices this
            # method performs best for. In the map situation (arguments have same shape, and
            # likely same or similar stored entry pattern), the Ai == Bi and termination
            # cases are equally or more likely than the Ai < Bi and Bi < Ai cases. Hence
            # the ordering of the conditional chain above differs from that in the
            # corresponding broadcast code (below).
            if !_iszero(Cx)
                Ck > spaceC && (spaceC = expandstorage!(C, Ck + (nnz(A) - (Ak - 1)) + (nnz(B) - (Bk - 1))))
                storedinds(C)[Ck] = Ci
                storedvals(C)[Ck] = Cx
                Ck += 1
            end
        end
    end
    @inbounds setcolptr!(C, numcols(C) + 1, Ck)
    trimstorage!(C, Ck - 1)
    return C
end
function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf
    # Build dense matrix structure in C, expanding storage if necessary
    _densestructure!(C)
    # Populate values
    fill!(storedvals(C), fillvalue)
    # NOTE: Combining this fill! into the loop below to avoid multiple sweeps over /
    # nonsequential access of storedvals(C) does not appear to improve performance.
    rowsentinelA = convert(indtype(A), numrows(A) + 1)
    rowsentinelB = convert(indtype(B), numrows(B) + 1)
    @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
        Ak, stopAk = colstartind(A, j), colboundind(A, j)
        Bk, stopBk = colstartind(B, j), colboundind(B, j)
        Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
        Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
        while true
            if Ai == Bi
                Ai == rowsentinelA && break # column complete
                Cx, Ci::indtype(C) = f(storedvals(A)[Ak], storedvals(B)[Bk]), Ai
                Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
            elseif Ai < Bi
                Cx, Ci = f(storedvals(A)[Ak], zero(eltype(B))), Ai
                Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
            else # Bi < Ai
                Cx, Ci = f(zero(eltype(A)), storedvals(B)[Bk]), Bi
                Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
            end
            Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
        end
    end
    return C
end

# (6) _map_zeropres!/_map_notzeropres! for more than two sparse matrices / vectors
function _map_zeropres!(f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N}
    spaceC::Int = min(length(storedinds(C)), length(storedvals(C)))
    rowsentinel = numrows(C) + 1
    Ck = 1
    stopks = _colstartind_all(1, As)
    @inbounds for j in columns(C)
        setcolptr!(C, j, Ck)
        ks = stopks
        stopks = _colboundind_all(j, As)
        rows = _rowforind_all(rowsentinel, ks, stopks, As)
        activerow = min(rows...)
        while activerow < rowsentinel
            vals, ks, rows = _fusedupdate_all(rowsentinel, activerow, rows, ks, stopks, As)
            Cx = f(vals...)
            if !_iszero(Cx)
                Ck > spaceC && (spaceC = expandstorage!(C, min(length(C), Ck + _sumnnzs(As...) - (sum(ks) - N))))
                storedinds(C)[Ck] = activerow
                storedvals(C)[Ck] = Cx
                Ck += 1
            end
            activerow = min(rows...)
        end
    end
    @inbounds setcolptr!(C, numcols(C) + 1, Ck)
    trimstorage!(C, Ck - 1)
    return C
end
function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N}
    # Build dense matrix structure in C, expanding storage if necessary
    _densestructure!(C)
    # Populate values
    fill!(storedvals(C), fillvalue)
    # NOTE: Combining this fill! into the loop below to avoid multiple sweeps over /
    # nonsequential access of nonzeros(C) does not appear to improve performance.
    rowsentinel = numrows(C) + 1
    stopks = _colstartind_all(1, As)
    @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
        ks = stopks
        stopks = _colboundind_all(j, As)
        rows = _rowforind_all(rowsentinel, ks, stopks, As)
        activerow = min(rows...)
        while activerow < rowsentinel
            vals, ks, rows = _fusedupdate_all(rowsentinel, activerow, rows, ks, stopks, As)
            Cx = f(vals...)
            Cx != fillvalue && (storedvals(C)[jo + activerow] = Cx)
            activerow = min(rows...)
        end
    end
    return C
end

# helper methods for map/map! methods just above
@inline _colstartind(j, A) = colstartind(A, j)
@inline _colstartind_all(j, ::Tuple{}) = ()
@inline _colstartind_all(j, As) = (
    _colstartind(j, first(As)),
    _colstartind_all(j, tail(As))...)
@inline _colboundind(j, A) = colboundind(A, j)
@inline _colboundind_all(j, ::Tuple{}) = ()
@inline _colboundind_all(j, As) = (
    _colboundind(j, first(As)),
    _colboundind_all(j, tail(As))...)
@inline _rowforind(rowsentinel, k, stopk, A) =
    k < stopk ? storedinds(A)[k] : convert(indtype(A), rowsentinel)
@inline _rowforind_all(rowsentinel, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
@inline _rowforind_all(rowsentinel, ks, stopks, As) = (
    _rowforind(rowsentinel, first(ks), first(stopks), first(As)),
    _rowforind_all(rowsentinel, tail(ks), tail(stopks), tail(As))...)

@inline function _fusedupdate(rowsentinel, activerow, row, k, stopk, A)
    # returns (val, nextk, nextrow)
    if row == activerow
        nextk = k + oneunit(k)
        (storedvals(A)[k], nextk, (nextk < stopk ? storedinds(A)[nextk] : oftype(row, rowsentinel)))
    else
        (zero(eltype(A)), k, row)
    end
end
@inline _fusedupdate_all(rowsentinel, activerow, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ((#=vals=#), (#=nextks=#), (#=nextrows=#))
@inline function _fusedupdate_all(rowsentinel, activerow, rows, ks, stopks, As)
    val, nextk, nextrow = _fusedupdate(rowsentinel, activerow, first(rows), first(ks), first(stopks), first(As))
    vals, nextks, nextrows = _fusedupdate_all(rowsentinel, activerow, tail(rows), tail(ks), tail(stopks), tail(As))
    return ((val, vals...), (nextk, nextks...), (nextrow, nextrows...))
end

# (7) _broadcast_zeropres!/_broadcast_notzeropres! specialized for a single (input) sparse vector/matrix
function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat) where Tf
    isempty(C) && return _finishempty!(C)
    spaceC::Int = min(length(storedinds(C)), length(storedvals(C)))
    # C and A cannot have the same shape, as we directed that case to map in broadcast's
    # entry point; here we need efficiently handle only heterogeneous C-A combinations where
    # one or both of C and A has at least one singleton dimension.
    #
    # We first divide the cases into two groups: those in which the input argument does not
    # expand vertically, and those in which the input argument expands vertically.
    #
    # Cases without vertical expansion
    Ck = 1
    if numrows(A) == numrows(C)
        @inbounds for j in columns(C)
            setcolptr!(C, j, Ck)
            bccolrangejA = numcols(A) == 1 ? colrange(A, 1) : colrange(A, j)
            for Ak in bccolrangejA
                Cx = f(storedvals(A)[Ak])
                if !_iszero(Cx)
                    Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A)))
                    storedinds(C)[Ck] = storedinds(A)[Ak]
                    storedvals(C)[Ck] = Cx
                    Ck += 1
                end
            end
        end
    # Cases with vertical expansion
    else # numrows(A) != numrows(C) (=> numrows(A) == 1)
        @inbounds for j in columns(C)
            setcolptr!(C, j, Ck)
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
            fofAx = f(Ax)
            # if fofAx is zero, then either A's jth column is empty, or A's jth column
            # contains a nonzero value x but f(Ax) is nonetheless zero, so we need store
            # nothing in C's jth column. if to the contrary fofAx is nonzero, then we must
            # densely populate C's jth column with fofAx.
            if !_iszero(fofAx)
                for Ci::indtype(C) in 1:numrows(C)
                    Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A)))
                    storedinds(C)[Ck] = Ci
                    storedvals(C)[Ck] = fofAx
                    Ck += 1
                end
            end
        end
    end
    @inbounds setcolptr!(C, numcols(C) + 1, Ck)
    trimstorage!(C, Ck - 1)
    return C
end
function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat) where Tf
    # For information on this code, see comments in similar code in _broadcast_zeropres! above
    # Build dense matrix structure in C, expanding storage if necessary
    _densestructure!(C)
    # Populate values
    fill!(storedvals(C), fillvalue)
    # Cases without vertical expansion
    if numrows(A) == numrows(C)
        @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
            bccolrangejA = numcols(A) == 1 ? colrange(A, 1) : colrange(A, j)
            for Ak in bccolrangejA
                Cx, Ci = f(storedvals(A)[Ak]), storedinds(A)[Ak]
                Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
            end
        end
    # Cases with vertical expansion
    else # numrows(A) != numrows(C) (=> numrows(A) == 1)
        svA, svC = storedvals(A), storedvals(C)
        @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Ax = Ak < stopAk ? svA[Ak] : zero(eltype(A))
            fofAx = f(Ax)
            if fofAx != fillvalue
                for i in (jo + 1):(jo + numrows(C))
                    svC[i] = fofAx
                end
            end
        end
    end
    return C
end

# (8) _broadcast_zeropres!/_broadcast_notzeropres! specialized for a pair of (input) sparse vectors/matrices
function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf
    isempty(C) && return _finishempty!(C)
    spaceC::Int = min(length(storedinds(C)), length(storedvals(C)))
    rowsentinelA = convert(indtype(A), numrows(C) + 1)
    rowsentinelB = convert(indtype(B), numrows(C) + 1)
    # C, A, and B cannot all have the same shape, as we directed that case to map in broadcast's
    # entry point; here we need efficiently handle only heterogeneous combinations of mats/vecs
    # with no singleton dimensions, one singleton dimension, and two singleton dimensions.
    # Cases involving objects with two singleton dimensions should be rare and optimizing
    # that case complicates the code appreciably, so we largely ignore that case's
    # performance below.
    #
    # We first divide the cases into two groups: those in which neither input argument
    # expands vertically, and those in which at least one argument expands vertically.
    #
    # NOTE: Placing the loops over columns outside the conditional chain segregating
    # argument shape combinations eliminates some code replication but unfortunately
    # hurts performance appreciably in some cases.
    #
    # Cases without vertical expansion
    Ck = 1
    if numrows(A) == numrows(B) == numrows(C)
        @inbounds for j in columns(C)
            setcolptr!(C, j, Ck)
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            # Restructuring this k/stopk code to avoid unnecessary colptr retrievals does
            # not improve performance signicantly. Leave in this less complex form.
            Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
            Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
            while true
                if Ai != Bi
                    if Ai < Bi
                        Cx, Ci = f(storedvals(A)[Ak], zero(eltype(B))), Ai
                        Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                    else # Ai > Bi
                        Cx, Ci = f(zero(eltype(A)), storedvals(B)[Bk]), Bi
                        Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                    end
                elseif #= Ai == Bi && =# Ai == rowsentinelA
                    break # column complete
                else #= Ai == Bi != rowsentinel =#
                    Cx, Ci::indtype(C) = f(storedvals(A)[Ak], storedvals(B)[Bk]), Ai
                    Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                    Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                end
                # NOTE: The ordering of the conditional chain above impacts which matrices
                # this method perform best for. In contrast to the map situation (arguments
                # have same shape, and likely same or similar stored entry pattern), where
                # the Ai == Bi and termination cases are equally or more likely than the
                # Ai < Bi and Bi < Ai cases, in the broadcast situation (arguments have
                # different shape, and likely largely disjoint expanded stored entry
                # pattern) the Ai < Bi and Bi < Ai cases are equally or more likely than the
                # Ai == Bi and termination cases. Hence the ordering of the conditional
                # chain above differs from that in the corresponding map code.
                if !_iszero(Cx)
                    Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A, B)))
                    storedinds(C)[Ck] = Ci
                    storedvals(C)[Ck] = Cx
                    Ck += 1
                end
            end
        end
    # Cases with vertical expansion
    elseif numrows(A) == numrows(B) == 1 # && numrows(C) != 1, vertically expand both A and B
        @inbounds for j in columns(C)
            setcolptr!(C, j, Ck)
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
            Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B))
            Cx = f(Ax, Bx)
            if !_iszero(Cx)
                for Ci::indtype(C) in 1:numrows(C)
                    Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A, B)))
                    storedinds(C)[Ck] = Ci
                    storedvals(C)[Ck] = Cx
                    Ck += 1
                end
            end
        end
    elseif numrows(A) == 1 # && numrows(B) == numrows(C) != 1 , vertically expand only A
        @inbounds for j in columns(C)
            setcolptr!(C, j, Ck)
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
            fvAzB = f(Ax, zero(eltype(B)))
            if _iszero(fvAzB)
                # either A's jth column is empty, or A's jth column contains a nonzero value
                # Ax but f(Ax, zero(eltype(B))) is nonetheless zero, so we can scan through
                # B's jth column without storing every entry in C's jth column
                while Bk < stopBk
                    Cx = f(Ax, storedvals(B)[Bk])
                    if !_iszero(Cx)
                        Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A, B)))
                        storedinds(C)[Ck] = storedinds(B)[Bk]
                        storedvals(C)[Ck] = Cx
                        Ck += 1
                    end
                    Bk += oneunit(Bk)
                end
            else
                # A's jth column is nonempty and f(Ax, zero(eltype(B))) is not zero, so
                # we must store (likely) every entry in C's jth column
                Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                for Ci::indtype(C) in 1:numrows(C)
                    if Bi == Ci
                        Cx = f(Ax, storedvals(B)[Bk])
                        Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                    else
                        Cx = fvAzB
                    end
                    if !_iszero(Cx)
                        Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A, B)))
                        storedinds(C)[Ck] = Ci
                        storedvals(C)[Ck] = Cx
                        Ck += 1
                    end
                end
            end
        end
    else # numrows(B) == 1 && numrows(A) == numrows(C) != 1, vertically expand only B
        @inbounds for j in columns(C)
            setcolptr!(C, j, Ck)
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B))
            fzAvB = f(zero(eltype(A)), Bx)
            if _iszero(fzAvB)
                # either B's jth column is empty, or B's jth column contains a nonzero value
                # Bx but f(zero(eltype(A)), Bx) is nonetheless zero, so we can scan through
                # A's jth column without storing every entry in C's jth column
                while Ak < stopAk
                    Cx = f(storedvals(A)[Ak], Bx)
                    if !_iszero(Cx)
                        Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A, B)))
                        storedinds(C)[Ck] = storedinds(A)[Ak]
                        storedvals(C)[Ck] = Cx
                        Ck += 1
                    end
                    Ak += oneunit(Ak)
                end
            else
                # B's jth column is nonempty and f(zero(eltype(A)), Bx) is not zero, so
                # we must store (likely) every entry in C's jth column
                Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                for Ci::indtype(C) in 1:numrows(C)
                    if Ai == Ci
                        Cx = f(storedvals(A)[Ak], Bx)
                        Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                    else
                        Cx = fzAvB
                    end
                    if !_iszero(Cx)
                        Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), A, B)))
                        storedinds(C)[Ck] = Ci
                        storedvals(C)[Ck] = Cx
                        Ck += 1
                    end
                end
            end
        end
    end
    @inbounds setcolptr!(C, numcols(C) + 1, Ck)
    trimstorage!(C, Ck - 1)
    return C
end
function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf
    # For information on this code, see comments in similar code in _broadcast_zeropres! above
    # Build dense matrix structure in C, expanding storage if necessary
    _densestructure!(C)
    # Populate values
    fill!(storedvals(C), fillvalue)
    rowsentinelA = convert(indtype(A), numrows(C) + 1)
    rowsentinelB = convert(indtype(B), numrows(C) + 1)
    # Cases without vertical expansion
    if numrows(A) == numrows(B) == numrows(C)
        @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
            Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
            while true
                if Ai < Bi
                    Cx, Ci = f(storedvals(A)[Ak], zero(eltype(B))), Ai
                    Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                elseif Ai > Bi
                    Cx, Ci = f(zero(eltype(A)), storedvals(B)[Bk]), Bi
                    Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                elseif #= Ai == Bi && =# Ai == rowsentinelA
                    break # column complete
                else #= Ai == Bi != rowsentinel =#
                    Cx, Ci::indtype(C) = f(storedvals(A)[Ak], storedvals(B)[Bk]), Ai
                    Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                    Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                end
                Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
            end
        end
    # Cases with vertical expansion
    elseif numrows(A) == numrows(B) == 1 # && numrows(C) != 1, vertically expand both A and B
        @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
            Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B))
            Cx = f(Ax, Bx)
            if Cx != fillvalue
                for Ck::Int in (jo + 1):(jo + numrows(C))
                    storedvals(C)[Ck] = Cx
                end
            end
        end
    elseif numrows(A) == 1 # && numrows(B) == numrows(C) != 1, vertically expand only A
        @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Ax = Ak < stopAk ? storedvals(A)[Ak] : zero(eltype(A))
            fvAzB = f(Ax, zero(eltype(B)))
            if fvAzB == fillvalue
                while Bk < stopBk
                    Cx = f(Ax, storedvals(B)[Bk])
                    Cx != fillvalue && (storedvals(C)[jo + storedinds(B)[Bk]] = Cx)
                    Bk += oneunit(Bk)
                end
            else
                Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                for Ci::indtype(C) in 1:numrows(C)
                    if Bi == Ci
                        Cx = f(Ax, storedvals(B)[Bk])
                        Bk += oneunit(Bk); Bi = Bk < stopBk ? storedinds(B)[Bk] : rowsentinelB
                    else
                        Cx = fvAzB
                    end
                    Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
                end
            end
        end
    else # numrows(B) == 1 && numrows(A) == numrows(C) != 1, vertically expand only B
        @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
            Ak, stopAk = numcols(A) == 1 ? (colstartind(A, 1), colboundind(A, 1)) : (colstartind(A, j), colboundind(A, j))
            Bk, stopBk = numcols(B) == 1 ? (colstartind(B, 1), colboundind(B, 1)) : (colstartind(B, j), colboundind(B, j))
            Bx = Bk < stopBk ? storedvals(B)[Bk] : zero(eltype(B))
            fzAvB = f(zero(eltype(A)), Bx)
            if fzAvB == fillvalue
                while Ak < stopAk
                    Cx = f(storedvals(A)[Ak], Bx)
                    Cx != fillvalue && (storedvals(C)[jo + storedinds(A)[Ak]] = Cx)
                    Ak += oneunit(Ak)
                end
            else
                Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                for Ci::indtype(C) in 1:numrows(C)
                    if Ai == Ci
                        Cx = f(storedvals(A)[Ak], Bx)
                        Ak += oneunit(Ak); Ai = Ak < stopAk ? storedinds(A)[Ak] : rowsentinelA
                    else
                        Cx = fzAvB
                    end
                    Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
                end
            end
        end
    end
    return C
end
_finishempty!(C::SparseVector) = C
_finishempty!(C::AbstractSparseMatrixCSC) = (fill!(getcolptr(C), 1); C)

# special case - vector outer product
_copy(f::typeof(*), x::SparseVectorUnion, y::AdjOrTransSparseVectorUnion) = _outer(x, y)
@inline _outer(x::SparseVectorUnion, y::Adjoint) = return _outer(conj, x, parent(y))
@inline _outer(x::SparseVectorUnion, y::Transpose) = return _outer(identity, x, parent(y))
function _outer(trans::Tf, x, y) where Tf
    nx = length(x)
    ny = length(y)
    rowvalx = nonzeroinds(x)
    rowvaly = nonzeroinds(y)
    nzvalsx = nonzeros(x)
    nzvalsy = nonzeros(y)
    nnzx = length(nzvalsx)
    nnzy = length(nzvalsy)

    nnzC = nnzx * nnzy
    Tv = typeof(oneunit(eltype(x)) * oneunit(eltype(y)))
    Ti = promote_type(indtype(x), indtype(y))
    colptrC = zeros(Ti, ny + 1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalsC = Vector{Tv}(undef, nnzC)

    idx = 0
    @inbounds colptrC[1] = 1
    @inbounds for jj = 1:nnzy
        yval = nzvalsy[jj]
        iszero(yval) && continue
        col = rowvaly[jj]
        yval = trans(yval)

        for ii = 1:nnzx
            xval = nzvalsx[ii]
            iszero(xval) && continue
            idx += 1
            colptrC[col+1] += 1
            rowvalC[idx] = rowvalx[ii]
            nzvalsC[idx] = xval * yval
        end
    end
    cumsum!(colptrC, colptrC)

    return SparseMatrixCSC(nx, ny, colptrC, rowvalC, nzvalsC)
end

# (9) _broadcast_zeropres!/_broadcast_notzeropres! for more than two (input) sparse vectors/matrices
function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N}
    isempty(C) && return _finishempty!(C)
    spaceC::Int = min(length(storedinds(C)), length(storedvals(C)))
    expandsverts = _expandsvert_all(C, As)
    expandshorzs = _expandshorz_all(C, As)
    rowsentinel = numrows(C) + 1
    Ck = 1
    @inbounds for j in columns(C)
        setcolptr!(C, j, Ck)
        ks = _startindforbccol_all(j, expandshorzs, As)
        stopks = _stopindforbccol_all(j, expandshorzs, As)
        # Neither fusing ks and stopks construction, nor restructuring them to avoid repeated
        # colptr lookups, improves performance significantly. So keep the less complex approach here.
        isemptys = _isemptycol_all(ks, stopks)
        defargs = _defargforcol_all(j, isemptys, expandsverts, ks, As)
        rows = _initrowforcol_all(j, rowsentinel, isemptys, expandsverts, ks, As)
        defaultCx = f(defargs...)
        activerow = min(rows...)
        if _iszero(defaultCx) # zero-preserving column scan
            while activerow < rowsentinel
                args, ks, rows = _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As)
                Cx = f(args...)
                if !_iszero(Cx)
                    Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), As)))
                    storedinds(C)[Ck] = activerow
                    storedvals(C)[Ck] = Cx
                    Ck += 1
                end
                activerow = min(rows...)
            end
        else # zero-non-preserving column scan
            for Ci in 1:numrows(C)
                if Ci == activerow
                    args, ks, rows = _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As)
                    Cx = f(args...)
                    activerow = min(rows...)
                else
                    Cx = defaultCx
                end
                if !_iszero(Cx)
                    Ck > spaceC && (spaceC = expandstorage!(C, _unchecked_maxnnzbcres(size(C), As)))
                    storedinds(C)[Ck] = Ci
                    storedvals(C)[Ck] = Cx
                    Ck += 1
                end
            end
        end
    end
    @inbounds setcolptr!(C, numcols(C) + 1, Ck)
    trimstorage!(C, Ck - 1)
    return C
end
function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N}
    isempty(C) && return _finishempty!(C)
    # Build dense matrix structure in C, expanding storage if necessary
    _densestructure!(C)
    # Populate values
    fill!(storedvals(C), fillvalue)
    expandsverts = _expandsvert_all(C, As)
    expandshorzs = _expandshorz_all(C, As)
    rowsentinel = numrows(C) + 1
    @inbounds for (j, jo) in zip(columns(C), _densecoloffsets(C))
        ks = _startindforbccol_all(j, expandshorzs, As)
        stopks = _stopindforbccol_all(j, expandshorzs, As)
        # Neither fusing ks and stopks construction, nor restructuring them to avoid repeated
        # colptr lookups, improves performance significantly. So keep the less complex approach here.
        isemptys = _isemptycol_all(ks, stopks)
        defargs = _defargforcol_all(j, isemptys, expandsverts, ks, As)
        rows = _initrowforcol_all(j, rowsentinel, isemptys, expandsverts, ks, As)
        defaultCx = f(defargs...)
        activerow = min(rows...)
        if defaultCx == fillvalue # fillvalue-preserving column scan
            while activerow < rowsentinel
                args, ks, rows = _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As)
                Cx = f(args...)
                Cx != fillvalue && (storedvals(C)[jo + activerow] = Cx)
                activerow = min(rows...)
            end
        else # fillvalue-non-preserving column scan
            for Ci in 1:numrows(C)
                if Ci == activerow
                    args, ks, rows = _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As)
                    Cx = f(args...)
                    activerow = min(rows...)
                else
                    Cx = defaultCx
                end
                Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx)
            end
        end
    end
    return C
end

# helper method for broadcast/broadcast! methods just above
@inline _expandsvert(C, A) = numrows(A) != numrows(C)
@inline _expandsvert_all(C, ::Tuple{}) = ()
@inline _expandsvert_all(C, As) = (_expandsvert(C, first(As)), _expandsvert_all(C, tail(As))...)
@inline _expandshorz(C, A) = numcols(A) != numcols(C)
@inline _expandshorz_all(C, ::Tuple{}) = ()
@inline _expandshorz_all(C, As) = (_expandshorz(C, first(As)), _expandshorz_all(C, tail(As))...)
@inline _startindforbccol(j, expandshorz, A) = expandshorz ? colstartind(A, 1) : colstartind(A, j)
@inline _startindforbccol_all(j, ::Tuple{}, ::Tuple{}) = ()
@inline _startindforbccol_all(j, expandshorzs, As) = (
    _startindforbccol(j, first(expandshorzs), first(As)),
    _startindforbccol_all(j, tail(expandshorzs), tail(As))...)
@inline _stopindforbccol(j, expandshorz, A) = expandshorz ? colboundind(A, 1) : colboundind(A, j)
@inline _stopindforbccol_all(j, ::Tuple{}, ::Tuple{}) = ()
@inline _stopindforbccol_all(j, expandshorzs, As) = (
    _stopindforbccol(j, first(expandshorzs), first(As)),
    _stopindforbccol_all(j, tail(expandshorzs), tail(As))...)
@inline _isemptycol(k, stopk) = k == stopk
@inline _isemptycol_all(::Tuple{}, ::Tuple{}) = ()
@inline _isemptycol_all(ks, stopks) = (
    _isemptycol(first(ks), first(stopks)),
    _isemptycol_all(tail(ks), tail(stopks))...)
@inline _initrowforcol(j, rowsentinel, isempty, expandsvert, k, A) =
    expandsvert || isempty ? convert(indtype(A), rowsentinel) : storedinds(A)[k]
@inline _initrowforcol_all(j, rowsentinel, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
@inline _initrowforcol_all(j, rowsentinel, isemptys, expandsverts, ks, As) = (
    _initrowforcol(j, rowsentinel, first(isemptys), first(expandsverts), first(ks), first(As)),
    _initrowforcol_all(j, rowsentinel, tail(isemptys), tail(expandsverts), tail(ks), tail(As))...)
@inline _defargforcol(j, isempty, expandsvert, k, A) =
    expandsvert && !isempty ? storedvals(A)[k] : zero(eltype(A))
@inline _defargforcol_all(j, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
@inline _defargforcol_all(j, isemptys, expandsverts, ks, As) = (
    _defargforcol(j, first(isemptys), first(expandsverts), first(ks), first(As)),
    _defargforcol_all(j, tail(isemptys), tail(expandsverts), tail(ks), tail(As))...)
@inline function _fusedupdatebc(rowsentinel, activerow, row, defarg, k, stopk, A)
    # returns (val, nextk, nextrow)
    if row == activerow
        nextk = k + oneunit(k)
        (storedvals(A)[k], nextk, (nextk < stopk ? storedinds(A)[nextk] : oftype(row, rowsentinel)))
    else
        (defarg, k, row)
    end
end
@inline _fusedupdatebc_all(rowsent, activerow, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ((#=vals=#), (#=nextks=#), (#=nextrows=#))
@inline function _fusedupdatebc_all(rowsentinel, activerow, rows, defargs, ks, stopks, As)
    val, nextk, nextrow = _fusedupdatebc(rowsentinel, activerow, first(rows), first(defargs), first(ks), first(stopks), first(As))
    vals, nextks, nextrows = _fusedupdatebc_all(rowsentinel, activerow, tail(rows), tail(defargs), tail(ks), tail(stopks), tail(As))
    return ((val, vals...), (nextk, nextks...), (nextrow, nextrows...))
end

# (10) broadcast over combinations of broadcast scalars and sparse vectors/matrices

# broadcast entry points for combinations of sparse arrays and other (scalar) types
@inline function copy(bc::Broadcasted{<:SPVM})
    bcf = flatten(bc)
    return _copy(bcf.f, bcf.args...)
end

_copy(f, args::SparseVector...) = _shapecheckbc(f, args...)
_copy(f, args::AbstractSparseMatrixCSC...) = _shapecheckbc(f, args...)
_copy(f, args::SparseVecOrMat...) = _diffshape_broadcast(f, args...)
# Otherwise, we incorporate scalars into the function and re-dispatch
function _copy(f, args...)
    parevalf, passedargstup = capturescalars(f, args)
    return _copy(parevalf, passedargstup...)
end
_copy(f) = throw(MethodError(_copy, (f,)))  # avoid method ambiguity

function _shapecheckbc(f, args...)
    _aresameshape(args...) ? _noshapecheck_map(f, args...) : _diffshape_broadcast(f, args...)
end

@inline function copyto!(dest::SparseVecOrMat, bc::Broadcasted{<:SPVM})
    if bc.f === identity && bc isa SpBroadcasted1 && Base.axes(dest) == (A = bc.args[1]; Base.axes(A))
        return copyto!(dest, A)
    end
    bcf = flatten(bc)
    As = map(arg->Base.unalias(dest, arg), bcf.args)
    return _copyto!(bcf.f, dest, As...)
end

@inline function _copyto!(f, dest, As::SparseVecOrMat...)
    _aresameshape(dest, As...) && return _noshapecheck_map!(f, dest, As...)
    Base.Broadcast.check_broadcast_axes(axes(dest), As...)
    fofzeros = f(_zeros_eltypes(As...)...)
    if _iszero(fofzeros)
        return _broadcast_zeropres!(f, dest, As...)
    else
        return _broadcast_notzeropres!(f, fofzeros, dest, As...)
    end
end

@inline function _copyto!(f, dest, args...)
    # args contains nothing but SparseVecOrMat and scalars
    # See below for capturescalars
    parevalf, passedsrcargstup = capturescalars(f, args)
    _copyto!(parevalf, dest, passedsrcargstup...)
end

# capturescalars takes a function (f) and a tuple of mixed sparse vectors/matrices and
# broadcast scalar arguments (mixedargs), and returns a function (parevalf, i.e. partially
# evaluated f) and a reduced argument tuple (passedargstup) containing only the sparse
# vectors/matrices in mixedargs in their original order, and such that the result of
# broadcast(parevalf, passedargstup...) is broadcast(f, mixedargs...)
@inline function capturescalars(f, mixedargs)
    let (passedsrcargstup, makeargs) = _capturescalars(mixedargs...)
        parevalf = (passed...) -> f(makeargs(passed...)...)
        return (parevalf, passedsrcargstup)
    end
end
# Work around losing Type{T}s as DataTypes within the tuple that makeargs creates
@inline capturescalars(f, mixedargs::Tuple{Ref{Type{T}}, Vararg{Any}}) where {T} =
    capturescalars((args...)->f(T, args...), Base.tail(mixedargs))
@inline capturescalars(f, mixedargs::Tuple{Ref{Type{T}}, Ref{Type{S}}, Vararg{Any}}) where {T, S} =
    # This definition is identical to the one above and necessary only for
    # avoiding method ambiguity.
    capturescalars((args...)->f(T, args...), Base.tail(mixedargs))
@inline capturescalars(f, mixedargs::Tuple{SparseVecOrMat, Ref{Type{T}}, Vararg{Any}}) where {T} =
    capturescalars((a1, args...)->f(a1, T, args...), (mixedargs[1], Base.tail(Base.tail(mixedargs))...))
@inline capturescalars(f, mixedargs::Tuple{Union{Ref,AbstractArray{<:Any,0}}, Ref{Type{T}}, Vararg{Any}}) where {T} =
    capturescalars((args...)->f(mixedargs[1], T, args...), Base.tail(Base.tail(mixedargs)))

nonscalararg(::SparseVecOrMat) = true
nonscalararg(::Any) = false
scalarwrappedarg(::Union{AbstractArray{<:Any,0},Ref}) = true
scalarwrappedarg(::Any) = false

@inline function _capturescalars()
    return (), () -> ()
end
@inline function _capturescalars(arg, mixedargs...)
    let (rest, f) = _capturescalars(mixedargs...)
        if nonscalararg(arg)
            return (arg, rest...), @inline function(head, tail...)
                (head, f(tail...)...)
            end # pass-through to broadcast
        elseif scalarwrappedarg(arg)
            return rest, @inline function(tail...)
                (arg[], f(tail...)...) # TODO: This can put a Type{T} in a tuple
            end # unwrap and add back scalararg after (in makeargs)
        else
            return rest, @inline function(tail...)
                (arg, f(tail...)...)
            end # add back scalararg after (in makeargs)
        end
    end
end
@inline function _capturescalars(arg) # this definition is just an optimization (to bottom out the recursion slightly sooner)
    if nonscalararg(arg)
        return (arg,), (head,) -> (head,) # pass-through
    elseif scalarwrappedarg(arg)
        return (), () -> (arg[],) # unwrap
    else
        return (), () -> (arg,) # add scalararg
    end
end

# NOTE: The following two method definitions work around #19096.
broadcast(f::Tf, ::Type{T}, A::AbstractSparseMatrixCSC) where {Tf,T} = broadcast(y -> f(T, y), A)
broadcast(f::Tf, A::AbstractSparseMatrixCSC, ::Type{T}) where {Tf,T} = broadcast(x -> f(x, T), A)

# (11) broadcast[!] over combinations of scalars, sparse vectors/matrices, structured matrices,
# and one- and two-dimensional Arrays (via promotion of structured matrices and Arrays)
#
# for combinations involving only scalars, sparse arrays, structured matrices, and dense
# vectors/matrices, promote all structured matrices and dense vectors/matrices to sparse
# and rebroadcast. otherwise, divert to generic AbstractArray broadcast code.

function copy(bc::Broadcasted{PromoteToSparse})
    bcf = flatten(bc)
    if can_skip_sparsification(bcf.f, bcf.args...)
        return _copy(bcf.f, bcf.args...)
    elseif is_supported_sparse_broadcast(bcf.args...)
        return _copy(bcf.f, map(_sparsifystructured, bcf.args)...)
    else
        return copy(convert(Broadcasted{Broadcast.DefaultArrayStyle{length(axes(bc))}}, bc))
    end
end

@inline function copyto!(dest::SparseVecOrMat, bc::Broadcasted{PromoteToSparse})
    bcf = flatten(bc)
    broadcast!(bcf.f, dest, map(_sparsifystructured, bcf.args)...)
end

_sparsifystructured(M::AbstractMatrix) = SparseMatrixCSC(M)
_sparsifystructured(V::AbstractVector) = SparseVector(V)
_sparsifystructured(M::AbstractSparseMatrix) = SparseMatrixCSC(M)
_sparsifystructured(V::AbstractSparseVector) = SparseVector(V)
_sparsifystructured(S::SparseVecOrMat) = S
_sparsifystructured(x) = x

# (12) map[!] over combinations of sparse and structured matrices
SparseOrStructuredMatrix = Union{SparseMatrixCSC,LinearAlgebra.StructuredMatrix}
map(f::Tf, A::SparseOrStructuredMatrix, Bs::Vararg{SparseOrStructuredMatrix,N}) where {Tf,N} =
    (_checksameshape(A, Bs...); _noshapecheck_map(f, _sparsifystructured(A), map(_sparsifystructured, Bs)...))
map!(f::Tf, C::AbstractSparseMatrixCSC, A::SparseOrStructuredMatrix, Bs::Vararg{SparseOrStructuredMatrix,N}) where {Tf,N} =
    (_checksameshape(C, A, Bs...); _noshapecheck_map!(f, C, _sparsifystructured(A), map(_sparsifystructured, Bs)...))

end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
"""

    Implements all three variants of the anomaly detectors based on k-nearest neighbors from
    From outliers to prototypes: Ordering data,
    Stefan Harmelinga and Guido Dornhegea and David Tax and Frank Meinecke and Klaus-Robert Muller, 2005

"""
struct KNNAnomaly{V<:Val}
    t::NNTree
    X::Matrix
    v::V
end

"""
    function KNNAnomaly(X::Matrix, v::Symbol, tree_type::Symbol = :BruteTree)

    create the k-nn anomaly detector with variant v::Symbol

"""
KNNAnomaly(X::Matrix, v::Symbol, tree_type::Symbol = :BruteTree) = KNNAnomaly(eval(tree_type)(X), X, Val(v))

"""
    kappa is the radius of a ball containing all points

"""
function StatsBase.predict(model::KNNAnomaly, x, k, v::V) where {V<:Val{:kappa}}
    inds, dists = NearestNeighbors.knn(model.t, x, k,true)
    map(d -> d[end],dists)
end

"""
    gamma is the average distance to all k-nearest neighbors

"""
function StatsBase.predict(model::KNNAnomaly, x, k, v::V) where {V<:Val{:gamma}}
    inds, dists = NearestNeighbors.knn(model.t, x, k)
    map(d -> Statistics.mean(d),dists)
end

"""
    delta is the average distance to all k-nearest neighbors

"""
function StatsBase.predict(model::KNNAnomaly, x, k, v::V) where {V<:Val{:delta}}
    inds, dists = NearestNeighbors.knn(model.t, x, k)
    map(i -> LinearAlgebra.norm(x[:,i[1]] - Statistics.mean(model.X[:,i[2]],dims=2)) ,enumerate(inds))
end

StatsBase.predict(model::KNNAnomaly, x, k, v::Symbol) = StatsBase.predict(model, x, k, Val(v))
StatsBase.predict(model::KNNAnomaly, x, k) = StatsBase.predict(model, x, k, model.v)


text:
# Entropy.jl

###### generalized entropy #####
"""
    gen_entropy(v, α)

    Compute the Generalized Entropy Index of a vector `v` at a specified parameter `α`.

# Examples
```julia
julia> using Inequality
julia> gen_entropy([8, 5, 1, 3, 5, 6, 7, 6, 3], 2)
0.09039256198347094
```
"""
gen_entropy(v::AbstractVector{<:Real}, α::Real)::Float64 = α == 1 ? theil(v) : (α == 0) ? (mld(v)) : (1/(α * (α - 1)) * sum(((v./Statistics.mean(v)).^α .- 1) * (1/length(v))) )

###### weighted generalized entropy #####
"""
    gen_entropy(v, w, α)

Compute the Generalized Entropy Index of a vector `v`, using weights given by a weight vector `w`
at a specified parameter `α`.

Weights must not be negative, missing or NaN. The weights and data vectors must have the same length.

# Examples
```julia
julia> using Inequality
julia> gen_entropy([8, 5, 1, 3, 5, 6, 7, 6, 3], collect(0.1:0.1:0.9), 2)
0.0709746654322026
```
"""
function gen_entropy(v::AbstractVector{<:Real},w::AbstractVector{<:Real}, α::Real)::Float64

    if α == 1
        return theil(v, w)
    elseif α == 0
        return mld(v, w)
    end

    checks_weights(v, w)

    return 1/(α * (α - 1)) * sum(((v./(sum(v .* w)/sum(w)) ) .^ α .-1) .* (w/sum(w)))

end

function gen_entropy(v::AbstractVector{<:Real},w::AbstractWeights, α::Real)::Float64

    if α == 1
        return theil(v, w)
    elseif α == 0
        return mld(v, w)
    end

    checks_weights(v, w)

    return 1/(α * (α - 1)) * sum(((v./(sum(v .* w)/w.sum) ) .^ α .-1) .* (w/w.sum))

end

"""
    wgen_entropy(v, w, α)

Compute the Generalized Entropy Index of `v` with weights `w` at a specified parameter `α`. See also [`gen_entropy`](gen_entropy)
"""
wgen_entropy(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, α::Real)::Float64 = gen_entropy(v,w,α)

wgen_entropy(v::AbstractVector{<:Real},  w::AbstractWeights, α::Real)::Float64 = gen_entropy(v,w,α)


text:
using AMQPClient

include("test_coverage.jl")
include("test_throughput.jl")
include("test_rpc.jl")

AMQPTestCoverage.runtests()
AMQPTestThroughput.runtests()
AMQPTestRPC.runtests()

if length(ARGS) > 0
    amqps_host = ARGS[1]
    virtualhost = ARGS[2]
    port = AMQPClient.AMQPS_DEFAULT_PORT

    login = ENV["AMQPPLAIN_LOGIN"]
    password = ENV["AMQPPLAIN_PASSWORD"]
    auth_params = Dict{String,Any}("MECHANISM"=>"AMQPLAIN", "LOGIN"=>login, "PASSWORD"=>password)

    AMQPTestCoverage.runtests(; host=amqps_host, port=AMQPClient.AMQPS_DEFAULT_PORT, virtualhost=virtualhost, amqps=amqps_configure(), auth_params=auth_params)
    AMQPTestThroughput.runtests(; host=amqps_host, port=AMQPClient.AMQPS_DEFAULT_PORT, tls=true)
    AMQPTestRPC.runtests(; host=amqps_host, port=AMQPClient.AMQPS_DEFAULT_PORT, amqps=amqps_configure())
end

exit(0)


text:
nrepeats_from_n(n) =
    if n < 2^22
        2^22 ÷ n
    else
        1
    end


text:
# FIXME: temporal patches to improve inferrability of Base methods
include("patches.jl")

# test utilities are extracted into a separate file for easier interactive testing from REPL
# i.e. julia -i test/interactive_utils.jl
include("interactive_utils.jl")

# stuff used across tests
using Test

const FIXTURE_DIR = normpath(@__DIR__, "fixtures")

const ERROR_REPORTS_FROM_SUM_OVER_STRING = let
    analyzer, = report_call(sum, (String,))
    @test !isempty(get_reports(analyzer))
    get_reports(analyzer)
end

function test_sum_over_string(ers)
    @test !isempty(ers)
    for target in ERROR_REPORTS_FROM_SUM_OVER_STRING
        @test any(ers) do er
            return er.msg == target.msg && er.sig == target.sig
        end
    end
end
test_sum_over_string(res::JET.VirtualProcessResult) = test_sum_over_string(res.inference_error_reports)
test_sum_over_string(analyzer::JET.AbstractAnalyzer)  = test_sum_over_string(get_reports(analyzer))

return