text:
# PART 5 – ASSESSMENT OF LOCAL METAL LOSS

using DataFrames
using CSV

# Determine Asessment Applicability
#Determine the assessment applicability
# @doc DesignCodeCriteria
# @doc MaterialToughness
# @doc CyclicService
# @doc Part5ComponentType
print("Begin -- Assessment Applicability and Component Type Checks\n")
creep_range = CreepRangeTemperature("Carbon Steel (UTS ≤ 414MPa (60 ksi))"; design_temperature=100.0, units="lbs-in-psi")
design = DesignCodeCriteria("ASME B31.3 Piping Code")
toughness = MaterialToughness("Certain")
cyclic = CyclicService(100, "Meets Part 14")
x = Part5ComponentType("Straight Section of Piping, Elbow or Bend - No Structural Attachments", vessel_orientation="horizontal", material="Carbon and Low Alloy Steels", D=24.75,Lss=120.0,H=0.0, NPS=10.0, design_temperature=100.0, units="lbs-in-psi")
part5_applicability = Part5AsessmentApplicability(x,design,toughness,cyclic,creep_range)

# For all assessments - determine the inspection data grid
M1 = [0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500]
M2 = [0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500]
M3 = [0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500]
M4 = [0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500]
M5 = [0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500]
M6 = [0.220, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500, 0.500]
CTPGrid = hcat(M6,M5,M4,M3,M2,M1) # build in descending order
CTPGrid = rotl90(CTPGrid) # rotate to correct orientation

# Level 1 fit for service
annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
# "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
# "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]
equipment_group = "piping" # "vessel", "tank"
flaw_location = "external" # "External","Internal"
pipe_code = "B31.8" # "B31.3", "B31.8"
FCA_string = "external"
metal_loss_categorization = "LTA" # "LTA" or "Groove-Like Flaw"
units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"
tmm_forcing = true
tmm = .165
tnom = 0.280 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.
trd = 0.280 # uniform thickness away from the local metal loss location established by thickness measurements at the time of the assessment.
FCAml = 0.05 # Future Corrosion Allowance applied to the region of metal loss.
FCA = 0.05 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).
LOSS = 0.045 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.
Do = 6.625 # Outside Diameter
D = Do - 2*(tnom) # Inside Dia.
P = 740.0 # internal design pressure.
S = 20000.0 # allowable stress.
E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.
MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.
F = 0.6 # B31.8 design factor
T = 1.0 # temperature factor B31.8
Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .
t = trd - LOSS - FCA # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.
tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).
spacings = 0.5 # spacings determine by visual inspection to adequately ccategorizse the corrosion -----------+ may add to CTP_Grid function for plotting purposes
# Flaw dimensions
s = 6.0 # longitudinal extent or length of the region of local metal loss based on future corroded thickness,
c = 2.0  # circumferential extent or length of the region of local metal loss (see Figure 5.2 and Figure 5.10), based on future corroded thickness, tc .
Ec = 1.0 # circumferential weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
El = 1.0 # longitudinal weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
RSFa = 0.9 # remaining strength factor - consult API 579 is go lower than 0.9

# For all assessments determine far enough from structural discontinuity
# Flaw-To-Major Structural Discontinuity Spacing
L1msd = [12.0] # distance to the nearest major structural discontinuity.
L2msd = [12.0] # distance to the nearest major structural discontinuity.
L3msd = [12.0] # distance to the nearest major structural discontinuity.
L4msd = [12.0] # distance to the nearest major structural discontinuity.
L5msd = [12.0] # distance to the nearest major structural discontinuity.
Lmsd = minimum([L1msd,L2msd,L3msd,L4msd,L5msd])[1]
if (Lmsd >= (1.8*(sqrt(D*(trd - LOSS - FCA)))))
    print("Satisfied - Flaw is located far enough from structural discontinuity\n")
    lmsd_satisfied = 1
else
    print("Not satisfied - Flaw is too close to a structural discontinuity - Conduct a level 3 assessment\n")
    lmsd_satisfied = 0
end

# Groove Like Flaw dimensions
gl = .05 # length of the Groove-Like Flaw based on future corroded condition.
gw = .4 # width of the Groove-Like Flaw based on future corroded condition.
gr = 0.1 # radius at the base of a Groove-Like Flaw based on future corroded condition.
β = 40.0 # see (Figure 5.4) :: orientation of the groove-like flaw with respect to the longitudinal axis or a parameter to compute an effective fracture toughness for a groove being evaluated as a crack-like flaw, as applicable.

FCAml_i = collect(0:0.001:.5)
RCA_out = zeros(19)
iterations_count = zeros(size(FCAml_i,1))
data_temp = FCAml_i[1]
i=1

# Iterate to find MAWPr maximum FCAml
for i = 1:size(FCAml_i,1)
    iterations_count[i] = i
    let RCA_out = RCA_out, data_temp = data_temp

    FCAml = FCAml_i[i]
    global part_5_lta_output = Part5LTALevel1(CTPGrid; remaining_life=remaining_life,tmm_forcing=tmm_forcing, tmm=tmm, annex2c_tmin_category=annex2c_tmin_category, pipe_code=pipe_code,equipment_group=equipment_group, flaw_location=flaw_location, FCA_string=FCA_string, metal_loss_categorization=metal_loss_categorization, units=units, Lmsd=Lmsd,tnom=tnom,
        trd=trd, FCA=FCA, FCAml=FCAml, LOSS=LOSS, Do=Do, D=D, P=P, S=S, E=E, MA=MA, Yb31=Yb31, t=t,tsl=tsl, F=F, T=T, spacings=spacings, s=s, c=c, El=El, Ec=Ec, RSFa=RSFa, gl=gl, gw=gw, gr=gr,β=β)

data_temp = part_5_lta_output[size(part_5_lta_output,1)+1:length(part_5_lta_output)]
RCA_out = append!(RCA_out , data_temp)

    end
end

max = Int64.(maximum(iterations_count))
if max == size(FCAml_i,1)
    RCA_out = reshape(RCA_out,19,max+1)
else
    RCA_out = reshape(RCA_out,19,max)
end
RCA_out = RCA_out[:,2:size(RCA_out,2)]

names_df = hcat(String.(part_5_lta_output[1:size(part_5_lta_output,1)]),RCA_out)

df = DataFrame(names_df)

CSV.write("C:/Users/Andrew.Bannerman/Desktop/MARS/15. PoC/3.24/PF-3001_FCA.csv", df;delim=',')

# account for variance in the MRWT
mrwt = tmm
upper = mrwt + .1
lower = mrwt - 0.1
x = collect(lower:0.0005:upper)
perc = (x .- mrwt) / mrwt *100
iterations_count = zeros(size(x,1))
data_out = zeros(19)
data_temp = x[1]

# Save datafrmaes

x_out = DataFrame(hcat(x))
perc_out = DataFrame(hcat(perc))

CSV.write("C:/Users/Andrew.Bannerman/Desktop/MARS/15. PoC/3.24/x_out.csv", x_out;delim=',')
CSV.write("C:/Users/Andrew.Bannerman/Desktop/MARS/15. PoC/3.24/perc_out.csv", perc_out;delim=',')

for i = 1:size(x,1)
    iterations_count[i] = i
    let data_out = data_out, data_temp = data_temp

    part_5_lta_output = Part5LTALevel1(CTPGrid; remaining_life=remaining_life,tmm_forcing=tmm_forcing, tmm=x[i], annex2c_tmin_category=annex2c_tmin_category, equipment_group=equipment_group, flaw_location=flaw_location, FCA_string=FCA_string, metal_loss_categorization=metal_loss_categorization, units=units, Lmsd=Lmsd, tnom=tnom,
        trd=trd, FCA=FCA, FCAml=FCAml, LOSS=LOSS, Do=Do, D=D, P=P, S=S, E=E, MA=MA, Yb31=Yb31, t=t,tsl=tsl, spacings=spacings, s=s, c=c, El=El, Ec=Ec, RSFa=RSFa, gl=gl, gw=gw, gr=gr,β=β)
data_temp = part_5_lta_output[size(part_5_lta_output,1)+1:length(part_5_lta_output)]
data_out = append!(data_out , data_temp)

    end
end

max = Int64.(maximum(iterations_count))
if max == size(x,1)
    data_out = reshape(data_out,19,max+1)
else
    data_out = reshape(data_out,19,max)
end
data_out = data_out[:,2:size(data_out,2)] # subet

names_df = hcat(String.(part_5_lta_output[1:19]),data_out)

df = DataFrame(names_df)

CSV.write("C:/Users/Andrew.Bannerman/Desktop/MARS/15. PoC/3.24/mawpr_thickness_sensitivty.csv", df;delim=',')


text:
using Pluto
using Test
using Markdown
using InteractiveUtils
using PlutoSerialization

module workspace100
module PlutoRunner
end
end
module var"workspace#200"
module PlutoRunner
end
end

@testset "PlutoSerialization.jl" begin
    @testset "workspace parsing" begin
        @test PlutoSerialization.ws_number(Module(:workspace)) === nothing
        @test PlutoSerialization.ws_number(Module(:workspace1)) == 1
        @test PlutoSerialization.ws_number(Module(Symbol("workspace#10"))) == 10
        @test PlutoSerialization.get_ws() == Dict(100 => workspace100, 200 => var"workspace#200")
    end

    @testset "roundtrip outside pluto" begin
        let buffer = IOBuffer()
            serialize(buffer, x -> x + 1)
            @test deserialize(seekstart(buffer))(1) == 2
        end
    end

    @testset "notebook" begin
        include("notebook.jl")
        s = Pluto.ServerSession()
        nb = Pluto.SessionActions.open(s, joinpath(@__DIR__, "notebook.jl"); run_async = false)
    end
end

ext: jl


text:
abstract type AbstractJaxTracer{T,N} <: AbstractJaxArray{T,N} end
#  o::PyObject
#  dims::NTuple{N,Int}
#end

const AbstractTracerVector{T} = AbstractJaxTracer{T,1}
const AbstractTracerMatrix{T} = AbstractJaxTracer{T,2}

Base.size(a::AbstractJaxTracer) = a.dims
Base.length(a::AbstractJaxTracer) = prod(a.dims)
Base.eltype(a::AbstractJaxTracer{T}) where {T} = T
Base.axes(a::AbstractJaxTracer) = map(Base.OneTo, size(a))

#wrap python properties
aval(a::AbstractJaxTracer) = a.o.aval

PyCall.PyObject(a::AbstractJaxTracer) = a.o

function Base.show(io::IO, ::MIME"text/plain", a::TA) where {TA<:AbstractJaxTracer}
  println(io, "$(_szstr(a.dims)) $TA:")
  print(io, " ", a.o)
end

function Base.show(io::IO, a::TA) where {TA<:AbstractJaxTracer}
  print(io, "$(_szstr(a.dims)) $TA")
end

for name in (:TracedArray, :ParallelTracerArray, :JVPTracerArray, :BatchTracerArray)
  @eval begin
    struct $name{T,N} <: AbstractJaxTracer{T,N}
      o::PyObject
      dims::NTuple{N,Int}
    end

    function $name(parr::PyObject)
      aval = parr.aval
      T = np_to_jl_type(aval.dtype)
      return $name{T,length(aval.shape)}(parr, reverse(aval.shape))
    end
    Base.convert(::Type{<:$name}, o::PyObject) = $name(o)
  end
end


text:
using LearnBase
import Base: cat, vcat, hcat

"""

"""
struct ArrayNode{A<:AbstractArray,C} <: AbstractNode
    data::A
    metadata::C
end

ArrayNode(data::AbstractMatrix) = ArrayNode(data, nothing)
# ArrayNode(data::AbstractNode, a...) = data

Flux.@functor ArrayNode

mapdata(f, x::ArrayNode) = ArrayNode(mapdata(f, x.data), x.metadata)

Base.ndims(x::ArrayNode) = Colon()
LearnBase.nobs(a::ArrayNode) = size(a.data, 2)
LearnBase.nobs(a::ArrayNode, ::Type{ObsDim.Last}) = nobs(a)

function reduce(::typeof(catobs), as::Vector{T}) where {T<:ArrayNode}
    data = reduce(catobs, [x.data for x in as])
    metadata = reduce(catobs, [a.metadata for a in as])
    ArrayNode(data, metadata)
end

# hcat and vcat only for ArrayNode
function Base.vcat(as::ArrayNode...)
    data = vcat([a.data for a in as]...)
    metadata = as[1].metadata == nothing ? nothing : as[1].metadata
    ArrayNode(data, metadata)
end

Base.hcat(as::ArrayNode...) = reduce(catobs, collect(as))

Base.getindex(x::ArrayNode, i::VecOrRange) = ArrayNode(subset(x.data, i), subset(x.metadata, i))

dsprint(io::IO, n::ArrayNode; pad=[], s="", tr=false) = paddedprint(io, "ArrayNode$(size(n.data))$(tr_repr(s, tr))")


text:
# Bar geometry summarizes data as vertical bars.
struct BarGeometry <: Gadfly.GeometryElement
    position::Symbol
    orientation::Symbol
    default_statistic::Gadfly.StatisticElement
    tag::Symbol
end
BarGeometry(; position=:stack, orientation=:vertical, tag=empty_tag) =
        BarGeometry(position, orientation, Gadfly.Stat.bar(position=position,
            orientation = orientation), tag)

element_aesthetics(geom::BarGeometry) = geom.orientation == :vertical ?
            [:xmin, :xmax, :y, :color, :alpha] : [:ymin, :ymax, :x, :color, :alpha]

default_statistic(geom::BarGeometry) = geom.default_statistic

"""
    Geom.bar[(; position=:stack, orientation=:vertical)]

Draw bars of height `y` centered at positions `x`, or from `xmin` to `xmax`.
If orientation is `:horizontal` switch x for y.  Optionally categorically
groups bars using the `color` aesthetic.  If `position` is `:stack` they will
be placed on top of each other;  if it is `:dodge` they will be placed side by
side.  For labelling dodged or stacked bar plots see [`Stat.dodge`](@ref).
If `position` is `:identity` then the bars can overlap (use the `alpha` aesthetic).
"""
const bar = BarGeometry

"""
    Geom.histogram[(; position=:stack, bincount=nothing, minbincount=3, maxbincount=150,
                    orientation=:vertical, density=false, limits=NamedTuple())]

Draw histograms from a series of observations in `x` or `y` optionally grouping
by `color`.  This geometry is equivalent to [`Geom.bar`](@ref) with
[`Stat.histogram`](@ref); see the latter for more information.
"""
histogram(; position=:stack, bincount=nothing,
            minbincount=3, maxbincount=150,
            orientation::Symbol=:vertical,
            density::Bool=false, limits::NamedTuple=NamedTuple(),
            tag::Symbol=empty_tag) =
    BarGeometry(position, orientation,
                Gadfly.Stat.histogram(bincount=bincount,
                                      minbincount=minbincount,
                                      maxbincount=maxbincount,
                                      position=position,
                                      orientation=orientation,
                                      density=density, limits=limits),
                tag)

# Render a single color bar chart
function render_bar(geom::BarGeometry,
                              theme::Gadfly.Theme,
                              aes::Gadfly.Aesthetics)
    rect = if geom.orientation == :horizontal
        XT = eltype(aes.x)
        xz = convert(XT, zero(XT))
        rectangle(min.(xz, aes.x),
            aes.ymin.*cy .- theme.bar_spacing/2,
            abs.(aes.x),
            (aes.ymax .- aes.ymin).*cy .- theme.bar_spacing, geom.tag)
    else
        YT = eltype(aes.y)
        yz = convert(YT, zero(YT))
        rectangle(aes.xmin.*cx .+ theme.bar_spacing/2,
            min.(yz, aes.y),
            (aes.xmax .- aes.xmin).*cx .- theme.bar_spacing,
            abs.(aes.y), geom.tag)
    end

    aes_alpha = eltype(aes.alpha) <: Int ? theme.alphas[aes.alpha] : aes.alpha
    sc =  if isa(theme.bar_highlight, Function)
        stroke(theme.bar_highlight(theme.default_color))
    elseif isa(theme.bar_highlight, Color)
        stroke(theme.bar_highlight)
    else
        stroke(nothing)
    end

    ctx = context()
    compose!(ctx, rect, fill(aes.color), sc, fillopacity(aes_alpha), svgclass("geometry"))
    return ctx
end

# Render a bar chart grouped by discrete colors and stacked.
function render_stacked_bar(geom::BarGeometry,
                                     theme::Gadfly.Theme,
                                     aes::Gadfly.Aesthetics,
                                     orientation::Symbol)
    # preserve factor orders of pooled data arrays
    if isa(aes.color, IndirectArray)
        idxs = sortperm(aes.color.index, rev=true)
    else
        idxs = 1:length(aes.color)
    end

    ctx = context()
    if orientation == :horizontal
        stack_height_dict = Dict()
        T = eltype(aes.x)
        z = convert(T, zero(T))
        for y in aes.ymin
            stack_height_dict[y] = z
        end
        stack_height = zeros(eltype(aes.x), length(idxs))

        for (i, j) in enumerate(idxs)
            if aes.x[j]>0
                stack_height[i] = stack_height_dict[aes.ymin[j]]
                stack_height_dict[aes.ymin[j]] += aes.x[j]
            else
                stack_height_dict[aes.ymin[j]] += aes.x[j]
                stack_height[i] = stack_height_dict[aes.ymin[j]]
            end
        end

        x0s = stack_height
        y0s = [aes.ymin[i]*cy + theme.bar_spacing/2 for i in idxs]
        widths = [abs(aes.x[i]) for i in idxs]
        heights = [(aes.ymax[i] - aes.ymin[i])*cy - theme.bar_spacing for i in idxs]
        compose!(ctx, rectangle(x0s, y0s, widths, heights, geom.tag))
    else
        stack_height_dict = Dict()
        T = eltype(aes.y)
        z = convert(T, zero(T))
        for x in aes.xmin
            stack_height_dict[x] = z
        end
        stack_height = zeros(eltype(aes.y), length(idxs))

        for (i, j) in enumerate(idxs)
            if aes.y[j]>0
                stack_height[i] = stack_height_dict[aes.xmin[j]]
                stack_height_dict[aes.xmin[j]] += aes.y[j]
            else
                stack_height_dict[aes.xmin[j]] += aes.y[j]
                stack_height[i] = stack_height_dict[aes.xmin[j]]
            end
        end

        x0s = [aes.xmin[i]*cx + theme.bar_spacing/2 for i in idxs]
        y0s = stack_height
        widths = [(aes.xmax[i] - aes.xmin[i])*cx - theme.bar_spacing for i in idxs]
        heights = [abs(aes.y[i]) for i in idxs]
        compose!(ctx, rectangle(x0s, y0s, widths, heights, geom.tag))
    end

    cs = [aes.color[i] for i in idxs]
    compose!(ctx, fill(cs), svgclass("geometry"))
    if isa(theme.bar_highlight, Function)
        compose!(ctx, stroke([theme.bar_highlight(c) for c in cs]))
    elseif isa(theme.bar_highlight, Color)
        compose!(ctx, stroke(theme.bar_highlight))
    else
        compose!(ctx, stroke(nothing))
    end
    return ctx
end

# Render a bar chart grouped by discrete colors and stuck next to each other.
function render_dodged_bar(geom::BarGeometry,
                                    theme::Gadfly.Theme,
                                    aes::Gadfly.Aesthetics,
                                    orientation::Symbol)
    # preserve factor orders of pooled data arrays
    if isa(aes.color, IndirectArray)
        idxs = sortperm(aes.color.index, rev=false)
    else
        idxs = 1:length(aes.color)
    end

    ctx = context()
    if orientation == :horizontal
        dodge_count = DefaultDict(() -> 0)
        for i in idxs
            dodge_count[aes.ymin[i]] += 1
        end

        dodge_height = Dict()
        dodge_pos_dict = DefaultDict(() -> 0cy)
        for i in idxs
            dodge_height[aes.ymin[i]] =
                ((aes.ymax[i] - aes.ymin[i]) / dodge_count[aes.ymin[i]]) * cy
            dodge_pos_dict[aes.ymin[i]] = aes.ymin[i]*cy
        end

        dodge_pos = Array{Measure}(undef, length(idxs))
        for (i, j) in enumerate(idxs)
            dodge_pos[i] = dodge_pos_dict[aes.ymin[j]] + theme.bar_spacing/2
            dodge_pos_dict[aes.ymin[j]] += dodge_height[aes.ymin[j]]
        end

        XT = eltype(aes.x)
        xz = convert(XT, zero(XT))

        aes_x = aes.x[idxs]
        compose!(ctx,
            rectangle([min(xz, x) for x in aes_x],
                dodge_pos,
                abs.(aes_x),
                [((aes.ymax[i] - aes.ymin[i])*cy - theme.bar_spacing) / dodge_count[aes.ymin[i]]
                 for i in idxs], geom.tag))
    else
        dodge_count = DefaultDict(() -> 0)
        for i in idxs
            dodge_count[aes.xmin[i]] += 1
        end

        dodge_width = Dict()
        dodge_pos_dict = DefaultDict(() -> 0cx)
        for i in idxs
            dodge_width[aes.xmin[i]] =
                ((aes.xmax[i] - aes.xmin[i]) / dodge_count[aes.xmin[i]]) * cx
            dodge_pos_dict[aes.xmin[i]] = aes.xmin[i]*cx
        end

        dodge_pos = Array{Measure}(undef, length(idxs))
        for (i, j) in enumerate(idxs)
            dodge_pos[i] = dodge_pos_dict[aes.xmin[j]] + theme.bar_spacing/2
            dodge_pos_dict[aes.xmin[j]] += dodge_width[aes.xmin[j]]
        end

        YT = eltype(aes.y)
        yz = convert(YT, zero(YT))

        aes_y = aes.y[idxs]
        compose!(ctx,
            rectangle(dodge_pos,
                [min(yz, y) for y in aes_y],
                [((aes.xmax[i] - aes.xmin[i])*cx - theme.bar_spacing) / dodge_count[aes.xmin[i]]
                 for i in idxs],
                abs.(aes_y), geom.tag))
    end

    cs = [aes.color[i] for i in idxs]
    compose!(ctx, fill(cs), svgclass("geometry"))
    if isa(theme.bar_highlight, Function)
        compose!(ctx, stroke([theme.bar_highlight(c) for c in cs]))
    elseif isa(theme.bar_highlight, Color)
        compose!(ctx, stroke(theme.bar_highlight))
    else
        compose!(ctx, stroke(nothing))
    end
    return ctx
end

# Render bar geometry.
#
# Args:
#   geom: bar geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::BarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    if geom.orientation == :horizontal
        Gadfly.assert_aesthetics_defined("BarGeometry", aes, :ymin, :ymax, :x)
        Gadfly.assert_aesthetics_equal_length("BarGeometry", aes, :ymin, :ymax, :x)
    elseif geom.orientation == :vertical
        Gadfly.assert_aesthetics_defined("BarGeometry", aes, :xmin, :xmax, :y)
        Gadfly.assert_aesthetics_equal_length("BarGeometry", aes, :xmin, :xmax, :y)
    else
        error("Orientation must be :horizontal or :vertical")
    end

    default_aes = Gadfly.Aesthetics()
    default_aes.color = RGBA{Float32}[theme.default_color]
    default_aes.alpha = Float64[theme.alphas[1]]
    aes = inherit(aes, default_aes)

    if geom.position==:identity || length(aes.color)==1
        ctx = render_bar(geom, theme, aes)
    elseif geom.position == :stack
        if geom.orientation == :horizontal
            for y in unique(aes.ymin)
               signs = map(sign, aes.x[aes.ymin.==y])
               in(1, signs) && in(-1, signs) &&
                       error("x values must be of the same sign for each unique y value")
            end
        elseif geom.orientation == :vertical
            for x in unique(aes.xmin)
               signs = map(sign, aes.y[aes.xmin.==x])
               in(1, signs) && in(-1, signs) &&
                       error("y values must be of the same sign for each unique x value")
            end
        end
        ctx = render_stacked_bar(geom, theme, aes, geom.orientation)
    elseif geom.position == :dodge
        ctx = render_dodged_bar(geom, theme, aes, geom.orientation)
    else
        error("$(geom.position) is not a valid position for the bar geometry.")
    end

    return compose!(context(), ctx,
        linewidth(theme.highlight_width),
        svgattribute("shape-rendering", "crispEdges"))
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 64af5a38-490a-11eb-34bc-37017a602974
using Revise, FFTResampling, Colors, PlutoUI, ImageShow

# ╔═╡ aec408dc-490a-11eb-0d5f-b167a9599994
md"#### Image Interpolation with FFT based sinc interpolation"

# ╔═╡ 18d15d82-50ee-11eb-1a86-abae54ad0050
begin
	N = 129
	x = range(-10.0, 10.0, length=33)[1:end-1]
	x_exact = range(-10.0, 10.0, length=513)[1:end-1]
	y = x'
	y_exact = x_exact'
	arr = sinc.(sqrt.(x .^2 .+ y .^2))
	arr_exact = sinc.(sqrt.(x_exact .^2 .+ y_exact .^2))
	arr_interp512 = resample(arr[1:end, 1:end], (512, 512));
end

# ╔═╡ 0a52fc98-76eb-11eb-2981-19c696ccd106
Gray.(arr)

# ╔═╡ fd0a8c24-504b-11eb-3e44-75fb8475f9da
md"### Slider to switch between images
"

# ╔═╡ 9245616e-490b-11eb-1648-c56a0155a111
img_all = cat(arr_exact, arr_interp512, dims=3);

# ╔═╡ 4c3e9134-490b-11eb-03a4-4d58d4d942da
md"
$(@bind index Slider(1:2))
"

# ╔═╡ 6a0cbf10-490b-11eb-10bd-359a3425da42
Gray.(img_all[:, :, index])#[100:300, 100:300, index])

# ╔═╡ Cell order:
# ╠═64af5a38-490a-11eb-34bc-37017a602974
# ╟─aec408dc-490a-11eb-0d5f-b167a9599994
# ╠═18d15d82-50ee-11eb-1a86-abae54ad0050
# ╠═0a52fc98-76eb-11eb-2981-19c696ccd106
# ╟─fd0a8c24-504b-11eb-3e44-75fb8475f9da
# ╠═9245616e-490b-11eb-1648-c56a0155a111
# ╟─4c3e9134-490b-11eb-03a4-4d58d4d942da
# ╠═6a0cbf10-490b-11eb-10bd-359a3425da42


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

module Broadcast

using .Base.Cartesian
using .Base: Indices, OneTo, linearindices, tail, to_shape,
            _msk_end, unsafe_bitgetindex, bitcache_chunks, bitcache_size, dumpbitcache,
            isoperator, promote_typejoin, unalias
import .Base: broadcast, broadcast!
export BroadcastStyle, broadcast_indices, broadcast_similar, broadcastable,
       broadcast_getindex, broadcast_setindex!, dotview, @__dot__

### Objects with customized broadcasting behavior should declare a BroadcastStyle

"""
`BroadcastStyle` is an abstract type and trait-function used to determine behavior of
objects under broadcasting. `BroadcastStyle(typeof(x))` returns the style associated
with `x`. To customize the broadcasting behavior of a type, one can declare a style
by defining a type/method pair

    struct MyContainerStyle <: BroadcastStyle end
    Base.BroadcastStyle(::Type{<:MyContainer}) = MyContainerStyle()

One then writes method(s) (at least [`broadcast_similar`](@ref)) operating on
`MyContainerStyle`. There are also several pre-defined subtypes of `BroadcastStyle`
that you may be able to leverage; see the
[Interfaces chapter](@ref man-interfaces-broadcasting) for more information.
"""
abstract type BroadcastStyle end

"""
`Broadcast.Style{C}()` defines a [`BroadcastStyle`](@ref) signaling through the type
parameter `C`. You can use this as an alternative to creating custom subtypes of `BroadcastStyle`,
for example

    Base.BroadcastStyle(::Type{<:MyContainer}) = Broadcast.Style{MyContainer}()

There is a pre-defined [`broadcast_similar`](@ref) method

    broadcast_similar(f, ::Style{C}, ::Type{ElType}, inds, args...) =
        similar(C, ElType, inds)

Naturally you can specialize this for your particular `C` (e.g., `MyContainer`).
"""
struct Style{T} <: BroadcastStyle end

BroadcastStyle(::Type{<:Tuple}) = Style{Tuple}()

struct Unknown <: BroadcastStyle end
BroadcastStyle(::Type{Union{}}) = Unknown()  # ambiguity resolution

"""
`Broadcast.AbstractArrayStyle{N} <: BroadcastStyle` is the abstract supertype for any style
associated with an `AbstractArray` type.
The `N` parameter is the dimensionality, which can be handy for AbstractArray types
that only support specific dimensionalities:

    struct SparseMatrixStyle <: Broadcast.AbstractArrayStyle{2} end
    Base.BroadcastStyle(::Type{<:SparseMatrixCSC}) = SparseMatrixStyle()

For AbstractArray types that support arbitrary dimensionality, `N` can be set to `Any`:

    struct MyArrayStyle <: Broadcast.AbstractArrayStyle{Any} end
    Base.BroadcastStyle(::Type{<:MyArray}) = MyArrayStyle()

In cases where you want to be able to mix multiple `AbstractArrayStyle`s and keep track
of dimensionality, your style needs to support a `Val` constructor:

    struct MyArrayStyleDim{N} <: Broadcast.AbstractArrayStyle{N} end
    (::Type{<:MyArrayStyleDim})(::Val{N}) where N = MyArrayStyleDim{N}()

Note that if two or more `AbstractArrayStyle` subtypes conflict, broadcasting machinery
will fall back to producing `Array`s. If this is undesirable, you may need to
define binary [`BroadcastStyle`](@ref) rules to control the output type.

See also [`Broadcast.DefaultArrayStyle`](@ref).
"""
abstract type AbstractArrayStyle{N} <: BroadcastStyle end

"""
`Broadcast.ArrayStyle{MyArrayType}()` is a [`BroadcastStyle`](@ref) indicating that an object
behaves as an array for broadcasting. It presents a simple way to construct
[`Broadcast.AbstractArrayStyle`](@ref)s for specific `AbstractArray` container types.
Broadcast styles created this way lose track of dimensionality; if keeping track is important
for your type, you should create your own custom [`Broadcast.AbstractArrayStyle`](@ref).
"""
struct ArrayStyle{A<:AbstractArray} <: AbstractArrayStyle{Any} end
(::Type{<:ArrayStyle{A}})(::Val) where A = A()

"""
`Broadcast.DefaultArrayStyle{N}()` is a [`BroadcastStyle`](@ref) indicating that an object
behaves as an `N`-dimensional array for broadcasting. Specifically, `DefaultArrayStyle` is
used for any
AbstractArray type that hasn't defined a specialized style, and in the absence of
overrides from other `broadcast` arguments the resulting output type is `Array`.
When there are multiple inputs to `broadcast`, `DefaultArrayStyle` "loses" to any other [`Broadcast.ArrayStyle`](@ref).
"""
struct DefaultArrayStyle{N} <: AbstractArrayStyle{N} end
(::Type{<:DefaultArrayStyle})(::Val{N}) where N = DefaultArrayStyle{N}()
const DefaultVectorStyle = DefaultArrayStyle{1}
const DefaultMatrixStyle = DefaultArrayStyle{2}
BroadcastStyle(::Type{<:AbstractArray{T,N}}) where {T,N} = DefaultArrayStyle{N}()
BroadcastStyle(::Type{<:Ref}) = DefaultArrayStyle{0}()
BroadcastStyle(::Type{T}) where {T} = DefaultArrayStyle{ndims(T)}()

# `ArrayConflict` is an internal type signaling that two or more different `AbstractArrayStyle`
# objects were supplied as arguments, and that no rule was defined for resolving the
# conflict. The resulting output is `Array`. While this is the same output type
# produced by `DefaultArrayStyle`, `ArrayConflict` "poisons" the BroadcastStyle so that
# 3 or more arguments still return an `ArrayConflict`.
struct ArrayConflict <: AbstractArrayStyle{Any} end

### Binary BroadcastStyle rules
"""
    BroadcastStyle(::Style1, ::Style2) = Style3()

Indicate how to resolve different `BroadcastStyle`s. For example,

    Broadcast.rule(::Primary, ::Secondary) = Primary()

would indicate that style `Primary` has precedence over `Secondary`.
You do not have to (and generally should not) define both argument orders.
The result does not have to be one of the input arguments, it could be a third type.

Please see the [Interfaces chapter](@ref man-interfaces-broadcasting) of the manual for
more information.
"""
BroadcastStyle(::S, ::S) where S<:BroadcastStyle = S() # homogeneous types preserved
# Fall back to Unknown. This is necessary to implement argument-swapping
BroadcastStyle(::BroadcastStyle, ::BroadcastStyle) = Unknown()
# Unknown loses to everything
BroadcastStyle(::Unknown, ::Unknown) = Unknown()
BroadcastStyle(::S, ::Unknown) where S<:BroadcastStyle = S()
# Precedence rules
BroadcastStyle(a::AbstractArrayStyle{0}, b::Style{Tuple}) = b
BroadcastStyle(a::AbstractArrayStyle, ::Style{Tuple})    = a
BroadcastStyle(::A, ::A) where A<:ArrayStyle             = A()
BroadcastStyle(::ArrayStyle, ::ArrayStyle)               = Unknown()
BroadcastStyle(::A, ::A) where A<:AbstractArrayStyle     = A()
Base.@pure function BroadcastStyle(a::A, b::B) where {A<:AbstractArrayStyle{M},B<:AbstractArrayStyle{N}} where {M,N}
    if Base.typename(A).wrapper == Base.typename(B).wrapper
        return A(_max(Val(M),Val(N)))
    end
    Unknown()
end
# Any specific array type beats DefaultArrayStyle
BroadcastStyle(a::AbstractArrayStyle{Any}, ::DefaultArrayStyle) = a
BroadcastStyle(a::AbstractArrayStyle{N}, ::DefaultArrayStyle{N}) where N = a
BroadcastStyle(a::AbstractArrayStyle{M}, ::DefaultArrayStyle{N}) where {M,N} =
    typeof(a)(_max(Val(M),Val(N)))

## Allocating the output container
"""
    broadcast_similar(f, ::BroadcastStyle, ::Type{ElType}, inds, As...)

Allocate an output object for [`broadcast`](@ref), appropriate for the indicated
[`Broadcast.BroadcastStyle`](@ref). `ElType` and `inds` specify the desired element type and indices of the
container.
`f` is the broadcast operation, and `As...` are the arguments supplied to `broadcast`.
"""
broadcast_similar(f, ::DefaultArrayStyle{N}, ::Type{ElType}, inds::Indices{N}, As...) where {N,ElType} =
    similar(Array{ElType}, inds)
broadcast_similar(f, ::DefaultArrayStyle{N}, ::Type{Bool}, inds::Indices{N}, As...) where N =
    similar(BitArray, inds)
# In cases of conflict we fall back on Array
broadcast_similar(f, ::ArrayConflict, ::Type{ElType}, inds::Indices, As...) where ElType =
    similar(Array{ElType}, inds)
broadcast_similar(f, ::ArrayConflict, ::Type{Bool}, inds::Indices, As...) =
    similar(BitArray, inds)

## Computing the result's indices. Most types probably won't need to specialize this.
broadcast_indices() = ()
broadcast_indices(::Type{T}) where T = ()
broadcast_indices(A) = broadcast_indices(combine_styles(A), A)
broadcast_indices(::Style{Tuple}, A) = (OneTo(length(A)),)
broadcast_indices(::DefaultArrayStyle{0}, A::Ref) = ()
broadcast_indices(::BroadcastStyle, A) = Base.axes(A)
"""
    Base.broadcast_indices(::SrcStyle, A)

Compute the indices for objects `A` with [`BroadcastStyle`](@ref) `SrcStyle`.
If needed, you can specialize this method for your styles.
You should only need to provide a custom implementation for non-AbstractArrayStyles.
"""
broadcast_indices

### End of methods that users will typically have to specialize ###

## Broadcasting utilities ##
# special cases defined for performance
broadcast(f, x::Number...) = f(x...)
@inline broadcast(f, t::NTuple{N,Any}, ts::Vararg{NTuple{N,Any}}) where {N} = map(f, t, ts...)

## logic for deciding the BroadcastStyle
# Dimensionality: computing max(M,N) in the type domain so we preserve inferrability
_max(V1::Val{Any}, V2::Val{Any}) = Val(Any)
_max(V1::Val{Any}, V2::Val{N}) where N = Val(Any)
_max(V1::Val{N}, V2::Val{Any}) where N = Val(Any)
_max(V1::Val, V2::Val) = __max(longest(ntuple(identity, V1), ntuple(identity, V2)))
__max(::NTuple{N,Bool}) where N = Val(N)
longest(t1::Tuple, t2::Tuple) = (true, longest(Base.tail(t1), Base.tail(t2))...)
longest(::Tuple{}, t2::Tuple) = (true, longest((), Base.tail(t2))...)
longest(t1::Tuple, ::Tuple{}) = (true, longest(Base.tail(t1), ())...)
longest(::Tuple{}, ::Tuple{}) = ()

# combine_styles operates on values (arbitrarily many)
combine_styles(c) = result_style(BroadcastStyle(typeof(c)))
combine_styles(c1, c2) = result_style(combine_styles(c1), combine_styles(c2))
@inline combine_styles(c1, c2, cs...) = result_style(combine_styles(c1), combine_styles(c2, cs...))

# result_style works on types (singletons and pairs), and leverages `BroadcastStyle`
result_style(s::BroadcastStyle) = s
result_style(s1::S, s2::S) where S<:BroadcastStyle = S()
# Test both orders so users typically only have to declare one order
result_style(s1, s2) = result_join(s1, s2, BroadcastStyle(s1, s2), BroadcastStyle(s2, s1))

# result_join is the final arbiter. Because `BroadcastStyle` for undeclared pairs results in Unknown,
# we defer to any case where the result of `BroadcastStyle` is known.
result_join(::Any, ::Any, ::Unknown, ::Unknown)   = Unknown()
result_join(::Any, ::Any, ::Unknown, s::BroadcastStyle) = s
result_join(::Any, ::Any, s::BroadcastStyle, ::Unknown) = s
# For AbstractArray types with specialized broadcasting and undefined precedence rules,
# we have to signal conflict. Because ArrayConflict is a subtype of AbstractArray,
# this will "poison" any future operations (if we instead returned `DefaultArrayStyle`, then for
# 3-array broadcasting the returned type would depend on argument order).
result_join(::AbstractArrayStyle, ::AbstractArrayStyle, ::Unknown, ::Unknown) =
    ArrayConflict()
# Fallbacks in case users define `rule` for both argument-orders (not recommended)
result_join(::Any, ::Any, ::S, ::S) where S<:BroadcastStyle = S()
@noinline function result_join(::S, ::T, ::U, ::V) where {S,T,U,V}
    error("""
conflicting broadcast rules defined
  Broadcast.BroadcastStyle(::$S, ::$T) = $U()
  Broadcast.BroadcastStyle(::$T, ::$S) = $V()
One of these should be undefined (and thus return Broadcast.Unknown).""")
end

# Indices utilities
combine_indices(A, B...) = broadcast_shape(broadcast_indices(A), combine_indices(B...))
combine_indices(A) = broadcast_indices(A)

# shape (i.e., tuple-of-indices) inputs
broadcast_shape(shape::Tuple) = shape
broadcast_shape(shape::Tuple, shape1::Tuple, shapes::Tuple...) = broadcast_shape(_bcs(shape, shape1), shapes...)
# _bcs consolidates two shapes into a single output shape
_bcs(::Tuple{}, ::Tuple{}) = ()
_bcs(::Tuple{}, newshape::Tuple) = (newshape[1], _bcs((), tail(newshape))...)
_bcs(shape::Tuple, ::Tuple{}) = (shape[1], _bcs(tail(shape), ())...)
function _bcs(shape::Tuple, newshape::Tuple)
    return (_bcs1(shape[1], newshape[1]), _bcs(tail(shape), tail(newshape))...)
end
# _bcs1 handles the logic for a single dimension
_bcs1(a::Integer, b::Integer) = a == 1 ? b : (b == 1 ? a : (a == b ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size"))))
_bcs1(a::Integer, b) = a == 1 ? b : (first(b) == 1 && last(b) == a ? b : throw(DimensionMismatch("arrays could not be broadcast to a common size")))
_bcs1(a, b::Integer) = _bcs1(b, a)
_bcs1(a, b) = _bcsm(b, a) ? b : (_bcsm(a, b) ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size")))
# _bcsm tests whether the second index is consistent with the first
_bcsm(a, b) = a == b || length(b) == 1
_bcsm(a, b::Number) = b == 1
_bcsm(a::Number, b::Number) = a == b || b == 1

## Check that all arguments are broadcast compatible with shape
# comparing one input against a shape
check_broadcast_shape(shp) = nothing
check_broadcast_shape(shp, ::Tuple{}) = nothing
check_broadcast_shape(::Tuple{}, ::Tuple{}) = nothing
check_broadcast_shape(::Tuple{}, Ashp::Tuple) = throw(DimensionMismatch("cannot broadcast array to have fewer dimensions"))
function check_broadcast_shape(shp, Ashp::Tuple)
    _bcsm(shp[1], Ashp[1]) || throw(DimensionMismatch("array could not be broadcast to match destination"))
    check_broadcast_shape(tail(shp), tail(Ashp))
end
check_broadcast_indices(shp, A) = check_broadcast_shape(shp, broadcast_indices(A))
# comparing many inputs
@inline function check_broadcast_indices(shp, A, As...)
    check_broadcast_indices(shp, A)
    check_broadcast_indices(shp, As...)
end

## Indexing manipulations

# newindex(I, keep, Idefault) replaces a CartesianIndex `I` with something that
# is appropriate for a particular broadcast array/scalar. `keep` is a
# NTuple{N,Bool}, where keep[d] == true means that one should preserve
# I[d]; if false, replace it with Idefault[d].
# If dot-broadcasting were already defined, this would be `ifelse.(keep, I, Idefault)`.
@inline newindex(I::CartesianIndex, keep, Idefault) = CartesianIndex(_newindex(I.I, keep, Idefault))
@inline _newindex(I, keep, Idefault) =
    (ifelse(keep[1], I[1], Idefault[1]), _newindex(tail(I), tail(keep), tail(Idefault))...)
@inline _newindex(I, keep::Tuple{}, Idefault) = ()  # truncate if keep is shorter than I

# newindexer(shape, A) generates `keep` and `Idefault` (for use by
# `newindex` above) for a particular array `A`, given the
# broadcast indices `shape`
# `keep` is equivalent to map(==, axes(A), shape) (but see #17126)
@inline newindexer(shape, A) = shapeindexer(shape, broadcast_indices(A))
@inline shapeindexer(shape, indsA::Tuple{}) = (), ()
@inline function shapeindexer(shape, indsA::Tuple)
    ind1 = indsA[1]
    keep, Idefault = shapeindexer(tail(shape), tail(indsA))
    (shape[1] == ind1, keep...), (first(ind1), Idefault...)
end

# Equivalent to map(x->newindexer(shape, x), As) (but see #17126)
map_newindexer(shape, ::Tuple{}) = (), ()
@inline function map_newindexer(shape, As)
    A1 = As[1]
    keeps, Idefaults = map_newindexer(shape, tail(As))
    keep, Idefault = newindexer(shape, A1)
    (keep, keeps...), (Idefault, Idefaults...)
end
@inline function map_newindexer(shape, A, Bs)
    keeps, Idefaults = map_newindexer(shape, Bs)
    keep, Idefault = newindexer(shape, A)
    (keep, keeps...), (Idefault, Idefaults...)
end

Base.@propagate_inbounds _broadcast_getindex(::Type{T}, I) where T = T
Base.@propagate_inbounds _broadcast_getindex(A, I) = _broadcast_getindex(combine_styles(A), A, I)
Base.@propagate_inbounds _broadcast_getindex(::DefaultArrayStyle{0}, A, I) = A[]
Base.@propagate_inbounds _broadcast_getindex(::Any, A, I) = A[I]
Base.@propagate_inbounds _broadcast_getindex(::Style{Tuple}, A::Tuple{Any}, I) = A[1]

## Broadcasting core
# nargs encodes the number of As arguments (which matches the number
# of keeps). The first two type parameters are to ensure specialization.
@generated function _broadcast!(f, B::AbstractArray, keeps::K, Idefaults::ID, A::AT, Bs::BT, ::Val{N}, iter) where {K,ID,AT,BT,N}
    nargs = N + 1
    quote
        $(Expr(:meta, :inline))
        # destructure the keeps and As tuples
        A_1 = A
        @nexprs $N i->(A_{i+1} = Bs[i])
        @nexprs $nargs i->(keep_i = keeps[i])
        @nexprs $nargs i->(Idefault_i = Idefaults[i])
        @simd for I in iter
            # reverse-broadcast the indices
            @nexprs $nargs i->(I_i = newindex(I, keep_i, Idefault_i))
            # extract array values
            @nexprs $nargs i->(@inbounds val_i = _broadcast_getindex(A_i, I_i))
            # call the function and store the result
            result = @ncall $nargs f val
            @inbounds B[I] = result
        end
        return B
    end
end

# For BitArray outputs, we cache the result in a "small" Vector{Bool},
# and then copy in chunks into the output
@generated function _broadcast!(f, B::BitArray, keeps::K, Idefaults::ID, A::AT, Bs::BT, ::Val{N}, iter) where {K,ID,AT,BT,N}
    nargs = N + 1
    quote
        $(Expr(:meta, :inline))
        # destructure the keeps and As tuples
        A_1 = A
        @nexprs $N i->(A_{i+1} = Bs[i])
        @nexprs $nargs i->(keep_i = keeps[i])
        @nexprs $nargs i->(Idefault_i = Idefaults[i])
        C = Vector{Bool}(undef, bitcache_size)
        Bc = B.chunks
        ind = 1
        cind = 1
        @simd for I in iter
            # reverse-broadcast the indices
            @nexprs $nargs i->(I_i = newindex(I, keep_i, Idefault_i))
            # extract array values
            @nexprs $nargs i->(@inbounds val_i = _broadcast_getindex(A_i, I_i))
            # call the function and store the result
            @inbounds C[ind] = @ncall $nargs f val
            ind += 1
            if ind > bitcache_size
                dumpbitcache(Bc, cind, C)
                cind += bitcache_chunks
                ind = 1
            end
        end
        if ind > 1
            @inbounds C[ind:bitcache_size] = false
            dumpbitcache(Bc, cind, C)
        end
        return B
    end
end

"""
    broadcastable(x)

Return either `x` or an object like `x` such that it supports `axes`, indexing, and its type supports `ndims`.

If `x` supports iteration, the returned value should have the same `axes` and indexing
behaviors as [`collect(x)`](@ref).

If `x` is not an `AbstractArray` but it supports `axes`, indexing, and its type supports
`ndims`, then `broadcastable(::typeof(x))` may be implemented to just return itself.
Further, if `x` defines its own [`BroadcastStyle`](@ref), then it must define its
`broadcastable` method to return itself for the custom style to have any effect.

# Examples
```jldoctest
julia> broadcastable([1,2,3]) # like `identity` since arrays already support axes and indexing
3-element Array{Int64,1}:
 1
 2
 3

julia> broadcastable(Int) # Types don't support axes, indexing, or iteration but are commonly used as scalars
Base.RefValue{Type{Int64}}(Int64)

julia> broadcastable("hello") # Strings break convention of matching iteration and act like a scalar instead
Base.RefValue{String}("hello")
```
"""
broadcastable(x::Union{Symbol,AbstractString,Function,UndefInitializer,Nothing,RoundingMode,Missing}) = Ref(x)
broadcastable(x::Ptr) = Ref{Ptr}(x) # Cannot use Ref(::Ptr) until ambiguous deprecation goes through
broadcastable(::Type{T}) where {T} = Ref{Type{T}}(T)
broadcastable(x::Union{AbstractArray,Number,Ref,Tuple}) = x
# In the future, default to collecting arguments. TODO: uncomment once deprecations are removed
# broadcastable(x) = collect(x)
# broadcastable(::Union{AbstractDict, NamedTuple}) = error("intentionally unimplemented to allow development in 1.x")

"""
    broadcast!(f, dest, As...)

Like [`broadcast`](@ref), but store the result of
`broadcast(f, As...)` in the `dest` array.
Note that `dest` is only used to store the result, and does not supply
arguments to `f` unless it is also listed in the `As`,
as in `broadcast!(f, A, A, B)` to perform `A[:] = broadcast(f, A, B)`.
"""
@inline function broadcast!(f::Tf, dest, As::Vararg{Any,N}) where {Tf,N}
    As′ = map(broadcastable, As)
    broadcast!(f, dest, combine_styles(As′...), As′...)
end
@inline broadcast!(f::Tf, dest, ::BroadcastStyle, As::Vararg{Any,N}) where {Tf,N} = broadcast!(f, dest, nothing, As...)

# Default behavior (separated out so that it can be called by users who want to extend broadcast!).
@inline function broadcast!(f, dest, ::Nothing, As::Vararg{Any, N}) where N
    if f isa typeof(identity) && N == 1
        A = As[1]
        if A isa AbstractArray && Base.axes(dest) == Base.axes(A)
            return copyto!(dest, A)
        end
    end
    _broadcast!(f, dest, As...)
    return dest
end

# Optimization for the case where all arguments are 0-dimensional
@inline function broadcast!(f, dest, ::AbstractArrayStyle{0}, As::Vararg{Any, N}) where N
    if dest isa AbstractArray
        if f isa typeof(identity) && N == 1
            return fill!(dest, As[1][])
        else
            @inbounds for I in eachindex(dest)
                dest[I] = f(map(getindex, As)...)
            end
            return dest
        end
    end
    _broadcast!(f, dest, As...)
    return dest
end

# For broadcasted assignments like `broadcast!(f, A, ..., A, ...)`, where `A`
# appears on both the LHS and the RHS of the `.=`, then we know we're only
# going to make one pass through the array, and even though `A` is aliasing
# against itself, the mutations won't affect the result as the indices on the
# LHS and RHS will always match. This is not true in general, but with the `.op=`
# syntax it's fairly common for an argument to be `===` a source.
broadcast_unalias(dest, src) = dest === src ? src : unalias(dest, src)

# This indirection allows size-dependent implementations.
@inline function _broadcast!(f, C, A, Bs::Vararg{Any,N}) where N
    shape = broadcast_indices(C)
    @boundscheck check_broadcast_indices(shape, A, Bs...)
    A′ = broadcast_unalias(C, A)
    Bs′ = map(B->broadcast_unalias(C, B), Bs)
    keeps, Idefaults = map_newindexer(shape, A′, Bs′)
    iter = CartesianIndices(shape)
    _broadcast!(f, C, keeps, Idefaults, A′, Bs′, Val(N), iter)
    return C
end

# broadcast with element type adjusted on-the-fly. This widens the element type of
# B as needed (allocating a new container and copying previously-computed values) to
# accommodate any incompatible new elements.
@generated function _broadcast!(f, B::AbstractArray, keeps::K, Idefaults::ID, As::AT, ::Val{nargs}, iter, st, count) where {K,ID,AT,nargs}
    quote
        $(Expr(:meta, :noinline))
        # destructure the keeps and As tuples
        @nexprs $nargs i->(A_i = As[i])
        @nexprs $nargs i->(keep_i = keeps[i])
        @nexprs $nargs i->(Idefault_i = Idefaults[i])
        while !done(iter, st)
            I, st = next(iter, st)
            # reverse-broadcast the indices
            @nexprs $nargs i->(I_i = newindex(I, keep_i, Idefault_i))
            # extract array values
            @nexprs $nargs i->(@inbounds val_i = _broadcast_getindex(A_i, I_i))
            # call the function
            V = @ncall $nargs f val
            # store the result
            if V isa eltype(B)
                @inbounds B[I] = V
            else
                # This element type doesn't fit in B. Allocate a new B with wider eltype,
                # copy over old values, and continue
                newB = Base.similar(B, promote_typejoin(eltype(B), typeof(V)))
                for II in Iterators.take(iter, count)
                    newB[II] = B[II]
                end
                newB[I] = V
                return _broadcast!(f, newB, keeps, Idefaults, As, Val(nargs), iter, st, count+1)
            end
            count += 1
        end
        return B
    end
end

maptoTuple(f) = Tuple{}
maptoTuple(f, a, b...) = Tuple{f(a), maptoTuple(f, b...).types...}

# An element type satisfying for all A:
# broadcast_getindex(
#     combine_styles(A),
#     A, broadcast_indices(A)
# )::_broadcast_getindex_eltype(A)
_broadcast_getindex_eltype(A) = _broadcast_getindex_eltype(combine_styles(A), A)
_broadcast_getindex_eltype(::BroadcastStyle, A) = eltype(A)  # Tuple, Array, etc.
_broadcast_getindex_eltype(::DefaultArrayStyle{0}, ::Ref{T}) where {T} = T

# Inferred eltype of result of broadcast(f, xs...)
combine_eltypes(f, A, As...) =
    Base._return_type(f, maptoTuple(_broadcast_getindex_eltype, A, As...))

"""
    broadcast(f, As...)

Broadcast the function `f` over the arrays, tuples, collections, `Ref`s and/or scalars `As`.

Broadcasting applies the function `f` over the elements of the container arguments and the
scalars themselves in `As`. Singleton and missing dimensions are expanded to match the
extents of the other arguments by virtually repeating the value. By default, only a limited
number of types are considered scalars, including `Number`s, `String`s, `Symbol`s, `Type`s,
`Function`s and some common singletons like `missing` and `nothing`. All other arguments are
iterated over or indexed into elementwise.

The resulting container type is established by the following rules:

 - If all the arguments are scalars or zero-dimensional arrays, it returns an unwrapped scalar.
 - If at least one argument is a tuple and all others are scalars or zero-dimensional arrays,
   it returns a tuple.
 - All other combinations of arguments default to returning an `Array`, but
   custom container types can define their own implementation and promotion-like
   rules to customize the result when they appear as arguments.

A special syntax exists for broadcasting: `f.(args...)` is equivalent to
`broadcast(f, args...)`, and nested `f.(g.(args...))` calls are fused into a
single broadcast loop.

# Examples
```jldoctest
julia> A = [1, 2, 3, 4, 5]
5-element Array{Int64,1}:
 1
 2
 3
 4
 5

julia> B = [1 2; 3 4; 5 6; 7 8; 9 10]
5×2 Array{Int64,2}:
 1   2
 3   4
 5   6
 7   8
 9  10

julia> broadcast(+, A, B)
5×2 Array{Int64,2}:
  2   3
  5   6
  8   9
 11  12
 14  15

julia> parse.(Int, ["1", "2"])
2-element Array{Int64,1}:
 1
 2

julia> abs.((1, -2))
(1, 2)

julia> broadcast(+, 1.0, (0, -2.0))
(1.0, -1.0)

julia> (+).([[0,2], [1,3]], Ref{Vector{Int}}([1,-1]))
2-element Array{Array{Int64,1},1}:
 [1, 1]
 [2, 2]

julia> string.(("one","two","three","four"), ": ", 1:4)
4-element Array{String,1}:
 "one: 1"
 "two: 2"
 "three: 3"
 "four: 4"

```
"""
@inline function broadcast(f, A, Bs...)
    A′ = broadcastable(A)
    Bs′ = map(broadcastable, Bs)
    broadcast(f, combine_styles(A′, Bs′...), nothing, nothing, A′, Bs′...)
end

# In the scalar case we unwrap the arguments and just call `f`
@inline broadcast(f, ::AbstractArrayStyle{0}, ::Nothing, ::Nothing, A, Bs...) = f(A[], map(getindex, Bs)...)

@inline broadcast(f, s::BroadcastStyle, ::Nothing, ::Nothing, A, Bs...) =
    broadcast(f, s, combine_eltypes(f, A, Bs...), combine_indices(A, Bs...), A, Bs...)

const NonleafHandlingTypes = Union{DefaultArrayStyle,ArrayConflict}

@inline function broadcast(f, s::NonleafHandlingTypes, ::Type{ElType}, inds::Indices, As...) where ElType
    if !Base.isconcretetype(ElType)
        return broadcast_nonleaf(f, s, ElType, inds, As...)
    end
    dest = broadcast_similar(f, s, ElType, inds, As...)
    broadcast!(f, dest, As...)
end

@inline function broadcast(f, s::BroadcastStyle, ::Type{ElType}, inds::Indices, As...) where ElType
    dest = broadcast_similar(f, s, ElType, inds, As...)
    broadcast!(f, dest, As...)
end

# When ElType is not concrete, use narrowing. Use the first element of each input to determine
# the starting output eltype; the _broadcast! method will widen `dest` as needed to
# accommodate later values.
function broadcast_nonleaf(f, s::NonleafHandlingTypes, ::Type{ElType}, shape::Indices, As...) where ElType
    nargs = length(As)
    iter = CartesianIndices(shape)
    if isempty(iter)
        return Base.similar(Array{ElType}, shape)
    end
    keeps, Idefaults = map_newindexer(shape, As)
    st = start(iter)
    I, st = next(iter, st)
    val = f([ _broadcast_getindex(As[i], newindex(I, keeps[i], Idefaults[i])) for i=1:nargs ]...)
    if val isa Bool
        dest = Base.similar(BitArray, shape)
    else
        dest = Base.similar(Array{typeof(val)}, shape)
    end
    dest[I] = val
    _broadcast!(f, dest, keeps, Idefaults, As, Val(nargs), iter, st, 1)
end

@inline broadcast(f, ::Style{Tuple}, ::Nothing, ::Nothing, A, Bs...) =
    tuplebroadcast(f, longest_tuple(A, Bs...), A, Bs...)
@inline tuplebroadcast(f, ::NTuple{N,Any}, As...) where {N} =
    ntuple(k -> f(tuplebroadcast_getargs(As, k)...), Val(N))
@inline tuplebroadcast(f, ::NTuple{N,Any}, ::Ref{Type{T}}, As...) where {N,T} =
    ntuple(k -> f(T, tuplebroadcast_getargs(As, k)...), Val(N))
longest_tuple(A::Tuple, B::Tuple, Bs...) = longest_tuple(_longest_tuple(A, B), Bs...)
longest_tuple(A, B::Tuple, Bs...) = longest_tuple(B, Bs...)
longest_tuple(A::Tuple, B, Bs...) = longest_tuple(A, Bs...)
longest_tuple(A, B, Bs...) = longest_tuple(Bs...)
longest_tuple(A::Tuple) = A
# Support only 1-tuples and N-tuples where there are no conflicts in N
_longest_tuple(A::Tuple{Any}, B::Tuple{Any}) = A
_longest_tuple(A::NTuple{N,Any}, B::NTuple{N,Any}) where N = A
_longest_tuple(A::NTuple{N,Any}, B::Tuple{Any}) where N = A
_longest_tuple(A::Tuple{Any}, B::NTuple{N,Any}) where N = B
@noinline _longest_tuple(A, B) =
    throw(DimensionMismatch("tuples $A and $B could not be broadcast to a common size"))

tuplebroadcast_getargs(::Tuple{}, k) = ()
@inline tuplebroadcast_getargs(As, k) =
    (_broadcast_getindex(first(As), k), tuplebroadcast_getargs(tail(As), k)...)

"""
    broadcast_getindex(A, inds...)

Equivalent to [`broadcast`](@ref)ing the `inds` arrays to a common size
and returning an array `[A[ks...] for ks in zip(indsb...)]` (where `indsb`
would be the broadcast `inds`). The shape of the output is equal to the shape of each
element of `indsb`.

# Examples
```jldoctest bc_getindex
julia> A = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> A[1:2, 1:2]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> broadcast_getindex(A, 1:2, 1:2)
2-element Array{Int64,1}:
 11
 22

julia> A[1:2, 2:-1:1]
2×2 Array{Int64,2}:
 12  11
 22  21

julia> broadcast_getindex(A, 1:2, 2:-1:1)
2-element Array{Int64,1}:
 12
 21
```
Because the indices are all vectors, these calls are like `[A[i[k], j[k]] for k = 1:2]`
where `i` and `j` are the two index vectors.

```jldoctest bc_getindex
julia> broadcast_getindex(A, 1:2, (1:2)')
2×2 Array{Int64,2}:
 11  12
 21  22

julia> broadcast_getindex(A, (1:2)', 1:2)
2×2 Array{Int64,2}:
 11  21
 12  22

julia> broadcast_getindex(A, [1 2 1; 1 2 2], [1, 2])
2×3 Array{Int64,2}:
 11  21  11
 12  22  22
```
"""
broadcast_getindex(src::AbstractArray, I::AbstractArray...) =
    broadcast_getindex!(Base.similar(Array{eltype(src)}, combine_indices(I...)),
                        src,
                        I...)

@generated function broadcast_getindex!(dest::AbstractArray, src::AbstractArray, I::AbstractArray...)
    N = length(I)
    Isplat = Expr[:(I[$d]) for d = 1:N]
    quote
        @nexprs $N d->(I_d = I[d])
        check_broadcast_indices(Base.axes(dest), $(Isplat...))  # unnecessary if this function is never called directly
        checkbounds(src, $(Isplat...))
        @nexprs $N d->(@nexprs $N k->(Ibcast_d_k = Base.axes(I_k, d) == OneTo(1)))
        @nloops $N i dest d->(@nexprs $N k->(j_d_k = Ibcast_d_k ? 1 : i_d)) begin
            @nexprs $N k->(@inbounds J_k = @nref $N I_k d->j_d_k)
            @inbounds (@nref $N dest i) = (@nref $N src J)
        end
        dest
    end
end

"""
    broadcast_setindex!(A, X, inds...)

Efficient element-by-element setting of the values of `A` in a pattern established by `inds`.
Equivalent to broadcasting the `X` and `inds` arrays to a common size, and then executing

    for (is, js) in zip(zip(indsb), eachindex(Xb))
        A[is...] = Xb[js...]
    end

where `Xb` and `indsb` are the broadcast `X` and `inds`.

See [`broadcast_getindex`](@ref) for examples of the treatment of `inds`.
"""
@generated function broadcast_setindex!(A::AbstractArray, x, I::AbstractArray...)
    N = length(I)
    Isplat = Expr[:(I[$d]) for d = 1:N]
    quote
        @nexprs $N d->(I_d = I[d])
        checkbounds(A, $(Isplat...))
        shape = combine_indices($(Isplat...))
        @nextract $N shape d->(length(shape) < d ? OneTo(1) : shape[d])
        @nexprs $N d->(@nexprs $N k->(Ibcast_d_k = Base.axes(I_k, d) == 1:1))
        if !isa(x, AbstractArray)
            xA = convert(eltype(A), x)
            @nloops $N i d->shape_d d->(@nexprs $N k->(j_d_k = Ibcast_d_k ? 1 : i_d)) begin
                @nexprs $N k->(@inbounds J_k = @nref $N I_k d->j_d_k)
                @inbounds (@nref $N A J) = xA
            end
        else
            X = x
            @nexprs $N d->(shapelen_d = length(shape_d))
            @ncall $N Base.setindex_shape_check X shapelen
            Xstate = start(X)
            @inbounds @nloops $N i d->shape_d d->(@nexprs $N k->(j_d_k = Ibcast_d_k ? 1 : i_d)) begin
                @nexprs $N k->(J_k = @nref $N I_k d->j_d_k)
                x_el, Xstate = next(X, Xstate)
                (@nref $N A J) = x_el
            end
        end
        A
    end
end

############################################################

# x[...] .= f.(y...) ---> broadcast!(f, dotview(x, ...), y...).
# The dotview function defaults to getindex, but we override it in
# a few cases to get the expected in-place behavior without affecting
# explicit calls to view.   (All of this can go away if slices
# are changed to generate views by default.)

Base.@propagate_inbounds dotview(args...) = Base.maybeview(args...)

############################################################
# The parser turns @. into a call to the __dot__ macro,
# which converts all function calls and assignments into
# broadcasting "dot" calls/assignments:

dottable(x) = false # avoid dotting spliced objects (e.g. view calls inserted by @view)
# don't add dots to dot operators
dottable(x::Symbol) = (!isoperator(x) || first(string(x)) != '.' || x === :..) && x !== :(:)
dottable(x::Expr) = x.head != :$
undot(x) = x
function undot(x::Expr)
    if x.head == :.=
        Expr(:(=), x.args...)
    elseif x.head == :block # occurs in for x=..., y=...
        Expr(:block, map(undot, x.args)...)
    else
        x
    end
end
__dot__(x) = x
function __dot__(x::Expr)
    dotargs = map(__dot__, x.args)
    if x.head == :call && dottable(x.args[1])
        Expr(:., dotargs[1], Expr(:tuple, dotargs[2:end]...))
    elseif x.head == :comparison
        Expr(:comparison, (iseven(i) && dottable(arg) && arg isa Symbol && isoperator(arg) ?
                               Symbol('.', arg) : arg for (i, arg) in pairs(dotargs))...)
    elseif x.head == :$
        x.args[1]
    elseif x.head == :let # don't add dots to `let x=...` assignments
        Expr(:let, undot(dotargs[1]), dotargs[2])
    elseif x.head == :for # don't add dots to for x=... assignments
        Expr(:for, undot(dotargs[1]), dotargs[2])
    elseif (x.head == :(=) || x.head == :function || x.head == :macro) &&
           Meta.isexpr(x.args[1], :call) # function or macro definition
        Expr(x.head, x.args[1], dotargs[2])
    else
        head = string(x.head)
        if last(head) == '=' && first(head) != '.'
            Expr(Symbol('.',head), dotargs...)
        else
            Expr(x.head, dotargs...)
        end
    end
end
"""
    @. expr

Convert every function call or operator in `expr` into a "dot call"
(e.g. convert `f(x)` to `f.(x)`), and convert every assignment in `expr`
to a "dot assignment" (e.g. convert `+=` to `.+=`).

If you want to *avoid* adding dots for selected function calls in
`expr`, splice those function calls in with `\$`.  For example,
`@. sqrt(abs(\$sort(x)))` is equivalent to `sqrt.(abs.(sort(x)))`
(no dot for `sort`).

(`@.` is equivalent to a call to `@__dot__`.)

# Examples
```jldoctest
julia> x = 1.0:3.0; y = similar(x);

julia> @. y = x + 3 * sin(x)
3-element Array{Float64,1}:
 3.5244129544236893
 4.727892280477045
 3.4233600241796016
```
"""
macro __dot__(x)
    esc(__dot__(x))
end

end # module

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function nlp_cvx_002_011(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)

    @objective(model, Min, (x - 3)^2 + (y - 2)^2)
    @constraint(model, 1 * x - 3 * y <= 3)
    @constraint(model, 1 * x - 5 * y <= 0)
    @constraint(model, 3 * x + 5 * y >= 15)
    @constraint(model, 7 * x + 2 * y >= 20)
    @constraint(model, 9 * x + 1 * y >= 20)
    @constraint(model, 3 * x + 7 * y >= 17)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 0, tol = objective_tol)
    return check_solution([x, y], [3, 2], tol = primal_tol)
end


text:
#!/usr/bin/env julia6
using MyStats, ReadGlobal, JLD2

function histindices(input::AbstractString,nbins::Integer)
    field = readfield(input)
    minf,maxf = min_max(field)
    ind = hist_indices(field,minf,maxf,nbins)
    bins = Bins(minf,maxf,nbins)
    @save "pdf_indices_$(input)_nbins$(nbins).jld2" ind bins
end

function main(f,ns)
    n = parse(Int,ns)
    histindices(f,n)
end

function main(f)
    n = 50
    histindices(f,n)
end

main(ARGS...)


text:
using DSGE
using HDF5, Base.Test
include("../util.jl")

path = dirname(@__FILE__)

# Initialize arguments to function
h5 = h5open("$path/../reference/kalman_filter_args.h5")
for arg in ["data", "lead", "a", "F", "b", "H", "var", "z0", "vz0"]
    eval(parse("$arg = read(h5, \"$arg\")"))
end
for arg in ["a", "b", "z0"]
    eval(parse("$arg = reshape(read(h5, \"$arg\"), length($arg), 1)"))
end

lead = round(Int,lead)

# Method with all arguments provided (9)
out_3 = kalman_filter(data, lead, a, F, b, H, var, z0, vz0)
out_9 = kalman_filter(data, lead, a, F, b, H, var, z0, vz0; allout=true)

h5 = h5open("$path/../reference/kalman_filter_out9.h5")
for out in [:L, :zend, :Pend, :pred, :vpred, :yprederror, :ystdprederror, :rmse,
            :rmsd, :filt, :vfilt]
    expect = read(h5, "$out")
    actual = out_9[out]

    if out == :zend
        expect = reshape(expect, length(expect), 1)
    end

    if ndims(expect) == 0
        @test_approx_eq expect actual
    else
        @test_matrix_approx_eq expect actual
    end
end
close(h5)

# Method with optional arguments omitted (7)
out_3 = kalman_filter(data, lead, a, F, b, H, var)
out_9 = kalman_filter(data, lead, a, F, b, H, var; allout=true)

h5 = h5open("$path/../reference/kalman_filter_out7.h5")
for out in [:L, :zend, :Pend, :pred, :vpred, :yprederror, :ystdprederror, :rmse,
            :rmsd, :filt, :vfilt]
    expect = read(h5, "$out")
    actual = out_9[out]

    if out == :zend
        expect = reshape(expect, length(expect), 1)
    end

    if out == :L
        @test_approx_eq_eps expect actual 1e-4
    elseif out ∈ [:Pend, :vpred, :vfilt]
        # These matrix entries are especially large, averaging 1e5, so we allow greater ϵ
        @test_matrix_approx_eq_eps expect actual 1e-1 1e-2
    else
        @test_matrix_approx_eq expect actual
    end
end
close(h5)

nothing


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

## Unary operators ##

"""
    conj!(A)

Transform an array to its complex conjugate in-place.

See also [`conj`](@ref).

# Examples
```jldoctest
julia> A = [1+im 2-im; 2+2im 3+im]
2×2 Array{Complex{Int64},2}:
 1+1im  2-1im
 2+2im  3+1im

julia> conj!(A);

julia> A
2×2 Array{Complex{Int64},2}:
 1-1im  2+1im
 2-2im  3-1im
```
"""
conj!(A::AbstractArray{<:Number}) = (@inbounds broadcast!(conj, A, A); A)

for f in (:-, :conj, :real, :imag)
    @eval ($f)(A::AbstractArray) = broadcast_preserving_zero_d($f, A)
end

## Binary arithmetic operators ##

for f in (:+, :-)
    @eval function ($f)(A::AbstractArray, B::AbstractArray)
        promote_shape(A, B) # check size compatibility
        broadcast_preserving_zero_d($f, A, B)
    end
end

function +(A::Array, Bs::Array...)
    for B in Bs
        promote_shape(A, B) # check size compatibility
    end
    broadcast_preserving_zero_d(+, A, Bs...)
end

for f in (:/, :\, :*)
    if f != :/
        @eval ($f)(A::Number, B::AbstractArray) = broadcast_preserving_zero_d($f, A, B)
    end
    if f != :\
        @eval ($f)(A::AbstractArray, B::Number) = broadcast_preserving_zero_d($f, A, B)
    end
end

## data movement ##

function reverse(A::Array{T}; dims::Integer) where T
    nd = ndims(A); d = dims
    1 ≤ d ≤ nd || throw(ArgumentError("dimension $d is not 1 ≤ $d ≤ $nd"))
    sd = size(A, d)
    if sd == 1 || isempty(A)
        return copy(A)
    end

    B = similar(A)

    nnd = 0
    for i = 1:nd
        nnd += Int(size(A,i)==1 || i==d)
    end
    if nnd==nd
        # reverse along the only non-singleton dimension
        for i = 1:sd
            B[i] = A[sd+1-i]
        end
        return B
    end

    d_in = size(A)
    leading = d_in[1:(d-1)]
    M = prod(leading)
    N = length(A)
    stride = M * sd

    if M==1
        for j = 0:stride:(N-stride)
            for i = 1:sd
                ri = sd+1-i
                B[j + ri] = A[j + i]
            end
        end
    else
        if allocatedinline(T) && M>200
            for i = 1:sd
                ri = sd+1-i
                for j=0:stride:(N-stride)
                    offs = j + 1 + (i-1)*M
                    boffs = j + 1 + (ri-1)*M
                    copyto!(B, boffs, A, offs, M)
                end
            end
        else
            for i = 1:sd
                ri = sd+1-i
                for j=0:stride:(N-stride)
                    offs = j + 1 + (i-1)*M
                    boffs = j + 1 + (ri-1)*M
                    for k=0:(M-1)
                        B[boffs + k] = A[offs + k]
                    end
                end
            end
        end
    end
    return B
end

"""
    rotl90(A)

Rotate matrix `A` left 90 degrees.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotl90(a)
2×2 Array{Int64,2}:
 2  4
 1  3
```
"""
function rotl90(A::AbstractMatrix)
    ind1, ind2 = axes(A)
    B = similar(A, (ind2,ind1))
    n = first(ind2)+last(ind2)
    for i=axes(A,1), j=ind2
        B[n-j,i] = A[i,j]
    end
    return B
end

"""
    rotr90(A)

Rotate matrix `A` right 90 degrees.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotr90(a)
2×2 Array{Int64,2}:
 3  1
 4  2
```
"""
function rotr90(A::AbstractMatrix)
    ind1, ind2 = axes(A)
    B = similar(A, (ind2,ind1))
    m = first(ind1)+last(ind1)
    for i=ind1, j=axes(A,2)
        B[j,m-i] = A[i,j]
    end
    return B
end
"""
    rot180(A)

Rotate matrix `A` 180 degrees.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rot180(a)
2×2 Array{Int64,2}:
 4  3
 2  1
```
"""
function rot180(A::AbstractMatrix)
    B = similar(A)
    ind1, ind2 = axes(A,1), axes(A,2)
    m, n = first(ind1)+last(ind1), first(ind2)+last(ind2)
    for j=ind2, i=ind1
        B[m-i,n-j] = A[i,j]
    end
    return B
end
"""
    rotl90(A, k)

Left-rotate matrix `A` 90 degrees counterclockwise an integer `k` number of times.
If `k` is a multiple of four (including zero), this is equivalent to a `copy`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotl90(a,1)
2×2 Array{Int64,2}:
 2  4
 1  3

julia> rotl90(a,2)
2×2 Array{Int64,2}:
 4  3
 2  1

julia> rotl90(a,3)
2×2 Array{Int64,2}:
 3  1
 4  2

julia> rotl90(a,4)
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
function rotl90(A::AbstractMatrix, k::Integer)
    k = mod(k, 4)
    k == 1 ? rotl90(A) :
    k == 2 ? rot180(A) :
    k == 3 ? rotr90(A) : copy(A)
end
"""
    rotr90(A, k)

Right-rotate matrix `A` 90 degrees clockwise an integer `k` number of times.
If `k` is a multiple of four (including zero), this is equivalent to a `copy`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rotr90(a,1)
2×2 Array{Int64,2}:
 3  1
 4  2

julia> rotr90(a,2)
2×2 Array{Int64,2}:
 4  3
 2  1

julia> rotr90(a,3)
2×2 Array{Int64,2}:
 2  4
 1  3

julia> rotr90(a,4)
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
rotr90(A::AbstractMatrix, k::Integer) = rotl90(A,-k)
"""
    rot180(A, k)

Rotate matrix `A` 180 degrees an integer `k` number of times.
If `k` is even, this is equivalent to a `copy`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> rot180(a,1)
2×2 Array{Int64,2}:
 4  3
 2  1

julia> rot180(a,2)
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
rot180(A::AbstractMatrix, k::Integer) = mod(k, 2) == 1 ? rot180(A) : copy(A)

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: milesfrain/julia

max_forks_repo_name: milesfrain/julia


text:
for (d, lattype) in enumerate((
        Chain,
        SquareLattice,
        CubicLattice
    ))
    constructor = d < 3 ? L -> lattype(L) : L -> lattype(d, L)

    @testset "$(d)D $lattype" begin
        for L in (3, 4)
            l = constructor(L)
            @test length(l) == L^d

            bonds = collect(neighbors(l, Val(true)))
            @test length(bonds) == 2*d*L^d
            @test allunique(bonds)
            for i in 0:length(l)-1
                # same source
                @test all(bonds[2d*i + 1][1] == x[1] for x in bonds[2d*i .+ (2:2d)])
                # different target
                @test allunique(x[2] for x in bonds[2d*i .+ (1:2d)])
            end

            reduced_bonds = collect(neighbors(l, Val(false)))
            @test length(reduced_bonds) == d*L^d
            # If directed is false, only one of (i, j) and (j, i)
            # should be kept.
            mirrored_bonds = [[trg, src] for (src, trg) in reduced_bonds]
            all_bonds = vcat(reduced_bonds, mirrored_bonds)
            @test sort(all_bonds) == sort(bonds)
        end
    end
end

using MonteCarlo: directed_norm

@testset "Lattice Iterators" begin
    # Set up a few test models/DQMCs
    m = HubbardModelAttractive(3, 3)
    dqmc1 = DQMC(m, beta=1.0)

    m = HubbardModelRepulsive(10, 1)
    dqmc2 = DQMC(m, beta=1.0)

    l = TriangularLattice(0, Lx=2, Ly=6)
    m = HubbardModelAttractive(l)
    dqmc3 = DQMC(m, beta=1.0)

    dqmcs = (dqmc1, dqmc2, dqmc3)

    @testset "Meta" begin
        for dqmc in dqmcs
            dirs = directions(dqmc)
            for i in 2:length(dirs)
                @test norm(dirs[i-1]) < norm(dirs[i]) + 1e-5
            end
        end

        dirs = [[cos(x), sin(x)] for x in range(0, 2pi-10eps(2pi), length=10)]
        norms = directed_norm.(dirs, 1e-6)
        for i in 2:length(norms)
            @test norms[i-1] < norms[i] < norms[i-1] + 1e-5
            @test 1.0 - 1e-5 < norms[i] < 1.0 + 1e-5
        end
    end

    @testset "EachSiteAndFlavor" begin
        for dqmc in dqmcs
            iter = EachSiteAndFlavor(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            Nflavors = MonteCarlo.nflavors(dqmc.model)
            @test collect(iter) == 1:Nsites*Nflavors
            @test length(iter) == Nsites*Nflavors
            @test eltype(iter) == Int64
            @test Base.IteratorSize(EachSiteAndFlavor) == Base.HasLength()
            @test Base.IteratorEltype(EachSiteAndFlavor) == Base.HasEltype()
        end
    end

    @testset "EachSite" begin
        for dqmc in dqmcs
            iter = EachSite(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            @test collect(iter) == 1:Nsites
            @test length(iter) == Nsites
            @test eltype(iter) == Int64
            @test Base.IteratorSize(EachSite) == Base.HasLength()
            @test Base.IteratorEltype(EachSite) == Base.HasEltype()
        end
    end

    @testset "OnSite" begin
        for dqmc in dqmcs
            iter = OnSite(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            @test collect(iter) == collect(zip(1:Nsites, 1:Nsites))
            @test length(iter) == Nsites
            @test eltype(iter) == Tuple{Int64, Int64}
            @test Base.IteratorSize(OnSite) == Base.HasLength()
            @test Base.IteratorEltype(OnSite) == Base.HasEltype()
        end
    end

    @testset "EachSitePair" begin
        for dqmc in dqmcs
            iter = EachSitePair(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            @test collect(iter) == [(i, j) for i in 1:Nsites for j in 1:Nsites]
            @test length(iter) == Nsites^2
            @test eltype(iter) == Tuple{Int64, Int64}
            @test Base.IteratorSize(EachSitePair) == Base.HasLength()
            @test Base.IteratorEltype(EachSitePair) == Base.HasEltype()
        end
    end

    @testset "EachSitePairByDistance" begin
        for dqmc in dqmcs
            iter = EachSitePairByDistance(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            @test length(iter) == Nsites^2
            @test eltype(iter) == NTuple{3, Int64}
            @test Base.IteratorSize(EachSitePairByDistance) == Base.HasLength()
            @test Base.IteratorEltype(EachSitePairByDistance) == Base.HasEltype()

            dirs = directions(dqmc)
            pos = MonteCarlo.positions(lattice(dqmc))
            wrap = MonteCarlo.generate_combinations(MonteCarlo.lattice_vectors(MonteCarlo.lattice(dqmc)))

            # Let's summarize these tests...
            check = true

            for (idx, src, trg) in iter
                _d = pos[src] - pos[trg]
                # Find lowest distance w/ periodic bounds
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check = check && (dirs[idx] ≈ d)
            end

            @test check

            _dirs = directions(iter, lattice(dqmc))
            check = true
            for (v, u) in zip(dirs, _dirs)
                check = check && (v ≈ u)
            end
            @test check
        end
    end

    @testset "EachLocalQuadByDistance" begin
        for dqmc in dqmcs
            iter = EachLocalQuadByDistance{6}(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            @test length(iter) == 6^2 * Nsites^2
            @test eltype(iter) == Tuple{Int64, UInt16, UInt16, UInt16, UInt16}
            @test Base.IteratorSize(EachLocalQuadByDistance) == Base.HasLength()
            @test Base.IteratorEltype(EachLocalQuadByDistance) == Base.HasEltype()

            dirs = directions(dqmc)
            pos = MonteCarlo.positions(lattice(dqmc))
            wrap = MonteCarlo.generate_combinations(MonteCarlo.lattice_vectors(MonteCarlo.lattice(dqmc)))

            # Let's summarize these tests...
            check12 = true
            check1 = true
            check2 = true

            idxs = CartesianIndices((length(dirs), 6, 6))

            for (lin, src1, trg1, src2, trg2) in iter
                idx12, idx1, idx2 = Tuple(idxs[lin])

                # src1 -- idx12 -- src2
                _d = pos[src1] - pos[src2]
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check12 = check12 && (dirs[idx12] ≈ d)

                # src1 -- idx1 -- trg1
                _d = pos[src1] - pos[trg1]
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check1 = check1 && (dirs[idx1] ≈ d)

                # src2 -- idx2 -- trg2
                _d = pos[src2] - pos[trg2]
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check2 = check2 && (dirs[idx2] ≈ d)
            end

            @test check12
            @test check1
            @test check2
        end
    end

    for dqmc in dqmcs
        @testset "EachLocalQuadBySyncedDistance" begin
            iter = EachLocalQuadBySyncedDistance{6}(dqmc, dqmc.model)
            Nsites = length(lattice(dqmc))
            # These are wrong on a lattice with a basis, because {6}
            # does not mean 6 surrounding sites, but rather 6 smallest directions (globally).
            # (A site directions may not apply to B sites)
            #@test length(iter) == 6^2 * Nsites^2
            #@test length(collect(iter)) == 6^2 * Nsites^2
            @test eltype(iter) == Tuple{Int64, UInt16, UInt16, UInt16, UInt16}
            @test Base.IteratorSize(EachLocalQuadBySyncedDistance) == Base.HasLength()
            @test Base.IteratorEltype(EachLocalQuadBySyncedDistance) == Base.HasEltype()

            dirs = directions(dqmc)
            pos = MonteCarlo.positions(lattice(dqmc))
            wrap = MonteCarlo.generate_combinations(MonteCarlo.lattice_vectors(MonteCarlo.lattice(dqmc)))

            # Let's summarize these tests...
            check12 = true
            check1 = true
            check2 = true

            idxs = CartesianIndices((length(dirs), 6))

            for (lin, src1, trg1, src2, trg2) in iter
                idx12, idx = Tuple(idxs[lin])

                # src1 -- idx12 -- src2
                _d = pos[src1] - pos[src2]
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check12 = check12 && (dirs[idx12] ≈ d)

                # src1 -- idx -- trg1
                _d = pos[src1] - pos[trg1]
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check1 = check1 && (dirs[idx] ≈ d)

                # src2 -- idx -- trg2
                _d = pos[src2] - pos[trg2]
                d = _d .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _d .+ v
                    if directed_norm(new_d, 1e-6) + 1e-6 < directed_norm(d, 1e-6)
                        d .= new_d
                    end
                end
                check2 = check2 && (dirs[idx] ≈ d)
            end

            @test check12
            @test check1
            @test check2
        end
    end

    @testset "Sum Wrapper" begin
        iter = EachSitePairByDistance(dqmc1, dqmc1.model)
        wrapped = Sum(iter)

        @test eltype(wrapped) == eltype(iter)
        @test length(wrapped) == length(iter)

        check = true
        vals = collect(iter)
        wals = collect(wrapped)
        for (v, w) in zip(vals, wals)
            check = check && (v == w)
        end
        @test check
    end

    @testset "Symmetry Wrapper" begin
        iter = EachSitePairByDistance(dqmc1, dqmc1.model)
        wrapped = ApplySymmetries{EachSitePairByDistance}([1], [0, 1, 1])(dqmc1, dqmc1.model)

        @test eltype(wrapped) == eltype(iter)
        @test length(wrapped) == length(iter)

        check = true
        vals = collect(iter)
        wals = collect(wrapped)
        for (v, w) in zip(vals, wals)
            check = check && (v == w)
        end
        @test check
    end

    @testset "SuperfluidDensity Wrapper" begin
        iter = EachSitePairByDistance(dqmc1, dqmc1.model)
        wrapped = MonteCarlo.SuperfluidDensity{EachSitePairByDistance}(
            1:3, [[1, 0], [0, 1], [1, 1]], [[0, 1], [1, 0], [1, -1]]
        )(dqmc1, dqmc1.model)

        @test eltype(wrapped) == eltype(iter)
        @test length(wrapped) == length(iter)

        check = true
        vals = collect(iter)
        wals = collect(wrapped)
        for (v, w) in zip(vals, wals)
            check = check && (v == w)
        end
        @test check
    end

    # TODO
    # SymmetryWrapper, SuperfluidDensity wrapper
end

# # L = 3 hardcoded Honeycomb
# struct HoneycombTestLattice <: AbstractLattice
#     positions::Vector{Vector{Float64}}
# end

# function HoneycombTestLattice()
#     # Taken from LatticePhysics (L = 3)
#     pos = [[0.0, 0.0], [0.5773502691896258, 0.0], [0.8660254037844386, 0.5], [1.4433756729740645, 0.5], [1.7320508075688772, 1.0], [2.309401076758503, 1.0], [0.8660254037844386, -0.5], [1.4433756729740645, -0.5], [1.7320508075688772, 0.0], [2.309401076758503, 0.0], [2.598076211353316, 0.5], [3.1754264805429417, 0.5], [1.7320508075688772, -1.0], [2.309401076758503, -1.0], [2.598076211353316, -0.5], [3.1754264805429413, -0.5], [3.4641016151377544, 0.0], [4.04145188432738, 0.0]]
#     HoneycombTestLattice(pos)
# end
# Base.length(l::HoneycombTestLattice) = length(l.positions)
# MonteCarlo.positions(l::HoneycombTestLattice) = l.positions
# function MonteCarlo.DistanceMask(lattice::HoneycombTestLattice)
#     L = 3
#     wrap = MonteCarlo.generate_combinations(
#         [L * [0.8660254037844386, -0.5], L * [0.8660254037844386, 0.5]]
#     )
#     MonteCarlo.VerboseDistanceMask(lattice, wrap)
# end
# function MonteCarlo.directions(mask::MonteCarlo.VerboseDistanceMask, lattice::HoneycombTestLattice)
#     pos = MonteCarlo.positions(lattice)
#     dirs = [pos[trg] - pos[src] for (src, trg) in first.(mask.targets)]
#     # marked = Set{Int64}()
#     # dirs = Vector{eltype(pos)}(undef, maximum(first(x) for x in mask.targets))
#     # for src in 1:size(mask.targets, 1)
#     #     for (idx, trg) in mask.targets[src, :]
#     #         if !(idx in marked)
#     #             push!(marked, idx)
#     #             dirs[idx] = pos[trg] - pos[src]
#     #         end
#     #     end
#     # end
#     L = 3
#     wrap = MonteCarlo.generate_combinations(
#         [L * [0.8660254037844386, -0.5], L * [0.8660254037844386, 0.5]]
#     )
#     map(dirs) do _d
#         d = round.(_d .+ wrap[1], digits=6)
#         for v in wrap[2:end]
#             new_d = round.(_d .+ v, digits=6)
#             if norm(new_d) < norm(d)
#                 d .= new_d
#             end
#         end
#         d
#     end
# end

# @testset "Mask with basis" begin
#     l = HoneycombTestLattice()
#     L = 3
#     positions = MonteCarlo.positions(l)
#     mask = DistanceMask(l)
#     @test mask isa MonteCarlo.VerboseDistanceMask
#     @test length(mask) == 27
#     @test allunique(MonteCarlo.getorder(mask))
#     # for i in 1:length(l)
#     #     @test allunique(mask[i, :])
#     #     @test allunique(mask[:, i])
#     # end
#     wrap = MonteCarlo.generate_combinations(
#         [L * [0.8660254037844386, -0.5], L * [0.8660254037844386, 0.5]]
#     )

#     # in the same direction
#     dirs = MonteCarlo.directions(mask, l)
#     for dir_idx in 1:length(mask)
#         _dirs = map(MonteCarlo.getdirorder(mask, dir_idx)) do (src, trg)
#             d = round.(positions[trg] .- positions[src] .+ wrap[1], digits=6)
#             for v in wrap[2:end]
#                 new_d = round.(positions[trg] .- positions[src] .+ v, digits=6)
#                 if norm(new_d) < norm(d)
#                     d .= new_d
#                 end
#             end
#             d
#         end
#         @test all(==(dirs[dir_idx]), _dirs)
#     end
#     # for src in 1:length(l)
#     #     for (idx, trg) in MonteCarlo.getorder(mask, src)
#     #         d = round.(positions[trg] .- positions[src] .+ wrap[1], digits=6)
#     #         for v in wrap[2:end]
#     #             new_d = round.(positions[trg] .- positions[src] .+ v, digits=6)
#     #             if norm(new_d) < norm(d)
#     #                 d .= new_d
#     #             end
#     #         end
#     #         @test dirs[idx] ≈ d
#     #     end
#     # end
# end

# @testset "2D Honeycomb" begin
#     for L in (2, 3)
#         l = MonteCarlo.HoneycombLattice(L)
#         @test length(l) == (2L)^2
#
#         bonds = collect(MonteCarlo.neighbors(l, Val(true)))
#         @test length(bonds) == 3 * (2L)^2
#         @test allunique(bonds)
#         for i in 0:length(l)-1
#             # same source
#             @test all(bonds[3i+1][1] == x[1] for x in bonds[3i .+ (2:3)])
#             # different target
#             @test allunique(x[2] for x in bonds[3i .+ (1:3)])
#         end
#
#         reduced_bonds = collect(MonteCarlo.neighbors(l, Val(false)))
#         @test length(reduced_bonds) == round(Int64, 1.5 * (2L)^2)
#         # If directed is false, only one of (i, j) and (j, i)
#         # should be kept.
#         mirrored_bonds = [[trg, src] for (src, trg) in reduced_bonds]
#         all_bonds = vcat(reduced_bonds, mirrored_bonds)
#         @test sort(all_bonds) == sort(bonds)
#
#         # For L < 3 NNN wrap around the lattice
#         L < 3 && continue
#         l = MonteCarlo.HoneycombLattice(L, true)
#         @test size(l.NNNs) == (6, length(l))
#         for i in 1:length(l)
#             @test allunique(l.NNNs[:, i])
#         end
#
#         positions = MonteCarlo.positions(l)
#         mask = MonteCarlo.DistanceMask(l)
#         @test mask isa MonteCarlo.VerboseDistanceMask
#         for i in 1:length(l)
#             @test allunique(mask[i, :])
#             @test allunique(mask[:, i])
#         end
#         wrap = MonteCarlo.generate_combinations(
#             [L * [0.8660254037844386, -0.5], L * [0.8660254037844386, 0.5]]
#         )
#
#         # in the same direction
#         dirs = MonteCarlo.directions(mask, l)
#         for src in 1:length(l)
#             for (idx, trg) in MonteCarlo.getorder(mask, src)
#                 d = round.(positions[trg] .- positions[src] .+ wrap[1], digits=6)
#                 for v in wrap[2:end]
#                     new_d = round.(positions[trg] .- positions[src] .+ v, digits=6)
#                     if norm(new_d) < norm(d)
#                         d .= new_d
#                     end
#                 end
#                 @test dirs[idx] ≈ d
#             end
#         end
#     end
# end

# # TODO: maybe if LatticePhysics gets registered one day?
# @testset "LatPhys" begin
#     for uc in [
#             getUnitcellSquare(), getUnitcellTriangular(),
#             getUnitcellKagome(), getUnitcellHoneycomb(),
#             getUnitcellSC(), getUnitcellFCC(),
#             getUnitcellDiamond(),
#             getUnitcell_9_3_a() # <- This one takes a bit longer
#         ]
#         l = getLatticePeriodic(uc, 4)
#         @info l
#         lattice = MonteCarlo.LatPhysLattice(l)
#
#         @test length(lattice) == numSites(l)
#         @test ndims(lattice) == ndims(l)
#
#         # ...
#
#         mask = MonteCarlo.DistanceMask(lattice)
#         for i in 1:length(lattice)
#             @test allunique(mask[i, :])
#             @test allunique(mask[:, i])
#         end
#
#         positions = point.(sites(l))
#         wrap = MonteCarlo.generate_combinations(latticeVectors(l))
#
#         # in the same direction
#         dirs = MonteCarlo.directions(mask, lattice)
#         for src in 1:length(lattice)
#             for (idx, trg) in MonteCarlo.getorder(mask, src)
#                 d = round.(positions[trg] .- positions[src] .+ wrap[1], digits=6)
#                 for v in wrap[2:end]
#                     new_d = round.(positions[trg] .- positions[src] .+ v, digits=6)
#                     if norm(new_d) < norm(d)
#                         d .= new_d
#                     end
#                 end
#                 @test dirs[idx] ≈ d
#             end
#         end
#     end
# end


text:
using CUDAdrv, CUDAnative, CuArrays

using Test

function vadd(gpu, a, b, c)
    i = threadIdx().x + blockDim().x * ((blockIdx().x-1) + (gpu-1) * gridDim().x)
    c[i] = a[i] + b[i]
    return
end

gpus = Int(length(devices()))

dims = (gpus,3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)

# FIXME: CuArray doesn't tie in with unified memory yet
buf_a = Mem.alloc(sizeof(a), true)
Mem.upload!(buf_a, a)
d_a = CuArray{Float32,3}(buf_a, dims)
buf_b = Mem.alloc(sizeof(a), true)
Mem.upload!(buf_b, b)
d_b = CuArray{Float32,3}(buf_b, dims)
buf_c = Mem.alloc(sizeof(a), true)
d_c = CuArray{Float32,3}(buf_c, dims)

len = prod(dims)
blocks = gpus
threads = len ÷ blocks

for (gpu,dev) in enumerate(devices())
    @debug "Allocating slice $gpu on device $(name(dev))"
    device!(dev)
    @cuda blocks=blocks÷gpus threads=threads vadd(gpu, d_a, d_b, d_c)
end

@debug "Synchronizing devices"
for dev in devices()
    # NOTE: normally you'd use events and wait for them
    device!(dev)
    synchronize()
end

c = Array(d_c)
@test a+b ≈ c


text:
export reluActivation, reluActivation!

"""
 relu activation A = relu(Y)

 Input:

   Y - array of features

 Optional Input:

   doDerivative - flag for computing derivative, set via varargin
                  Ex: reluActivation(Y,true);

 Output:

  A  - activation
  dA - derivatives
"""
function reluActivation(Y::Array{T},doDerivative::Bool=false) where {T}

    A = max.(Y,zero(T));

    if doDerivative
        dA = sign.(A);
    else
        dA = zeros(T,0)
    end

    return A,dA
end

function reluActivation!(A::Array{T},dA,doDerivative::Bool=false) where {T}
    A .= max.(A,zero(T));
    if doDerivative
        if isempty(dA)
            dA = sign.(A);
        else
            dA .= sign.(A);
        end
    end
    return A,dA
end


text:
export coordarray,
    rmsd,
    displacements,
    distance

"""
Get the atomic coordinates of a `StructuralElementOrList` as a 2D `Array` with
each column corresponding to one atom.
"""
function coordarray(element::StructuralElementOrList, selector_functions::Function...)
    atoms = collectatoms(element, selector_functions...)
    coords = zeros(3, length(atoms))
    for j in eachindex(atoms)
        coords[1,j] = x(atoms[j])
        coords[2,j] = y(atoms[j])
        coords[3,j] = z(atoms[j])
    end
    return coords
end

coordarray(coords::Array{Float64}, selector_functions::Function...) = coords

"""
Get the root-mean-square deviation (RMSD) between two `StructuralElementOrList`s
or two coordinate `Array`s of the same size. Assumes they are already
superimposed.
"""
function rmsd(coords_one::Array{Float64}, coords_two::Array{Float64})
    @assert size(coords_one) == size(coords_two) "Sizes of coordinate arrays are different - cannot calculate RMSD"
    difference = coords_one - coords_two
    return sqrt(sum(difference .* difference) / size(coords_one, 2))
end

rmsd(element_one::StructuralElementOrList, element_two::StructuralElementOrList, selector_functions::Function...) = rmsd(coordarray(element_one, selector_functions...), coordarray(element_two, selector_functions...))

"""
Get the displacements between atomic coordinates from two
`StructuralElementOrList`s or two coordinate `Array`s of the same size. Assumes
they are already superimposed.
"""
function displacements(coords_one::Array{Float64}, coords_two::Array{Float64})
    @assert size(coords_one) == size(coords_two) "Sizes of coordinate arrays are different - cannot calculate displacements"
    difference = coords_one - coords_two
    return sqrt(sum(difference .* difference, 1))[:]
end

displacements(element_one::StructuralElementOrList, element_two::StructuralElementOrList, selector_functions::Function...) = displacements(coordarray(element_one, selector_functions...), coordarray(element_two, selector_functions...))

"Get the minimum distance between two `StructuralElementOrList`s."
function distance(element_one::StructuralElementOrList, element_two::StructuralElementOrList, selector_functions::Function...)
    coords_one = coordarray(element_one, selector_functions...)
    coords_two = coordarray(element_two, selector_functions...)
    count_one = size(coords_one, 2)
    count_two = size(coords_two, 2)
    sq_dists = zeros(count_one, count_two)
    for i in 1:count_one, j in 1:count_two
        @inbounds sq_dists[i,j] = (coords_one[1,i] - coords_two[1,j]) ^ 2 + (coords_one[2,i] - coords_two[2,j]) ^ 2 + (coords_one[3,i] - coords_two[3,j]) ^ 2
    end
    return sqrt(minimum(sq_dists))
end

distance(atom_one::AbstractAtom, atom_two::AbstractAtom) = sqrt((x(atom_one) - x(atom_two)) ^ 2 + (y(atom_one) - y(atom_two)) ^ 2 + (z(atom_one) - z(atom_two)) ^ 2)


text:
# Default nlsolve behavior, should move to FiniteDiff.jl

Base.@pure determine_chunksize(u,alg::DiffEqBase.DEAlgorithm) = determine_chunksize(u,get_chunksize(alg))
Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

macro swap!(x,y)
  quote
    local tmp = $(esc(x))
    $(esc(x)) = $(esc(y))
    $(esc(y)) = tmp
  end
end

macro cache(expr)
  name = expr.args[2].args[1].args[1]
  fields = [x for x in expr.args[3].args if typeof(x)!=LineNumberNode]
  cache_vars = Expr[]
  jac_vars = Pair{Symbol,Expr}[]
  for x in fields
    if x.args[2] == :uType || x.args[2] == :rateType ||
       x.args[2] == :kType || x.args[2] == :uNoUnitsType
      push!(cache_vars,:(c.$(x.args[1])))
    elseif x.args[2] == :DiffCacheType
      push!(cache_vars,:(c.$(x.args[1]).du))
      push!(cache_vars,:(c.$(x.args[1]).dual_du))
    end
  end
  quote
    $expr
    $(esc(:full_cache))(c::$name) = tuple($(cache_vars...))
  end
end

# Nest one layer of value in order to get rid of possible Dual{Complex} or Complex{Dual} issues
# value should recurse for anything else.
function constvalue(::Type{T}) where T
  _T = DiffEqBase.value(T)
  return _T <: Complex ? DiffEqBase.value(real(_T)) : DiffEqBase.value(_T)
end
function constvalue(x)
  _x = DiffEqBase.value(x)
  return _x isa Complex ? DiffEqBase.value(real(_x)) : DiffEqBase.value(_x)
end

function diffdir(integrator::DiffEqBase.DEIntegrator)
  difference = maximum(abs, integrator.uprev)*sqrt(eps(typeof(integrator.t)))
  dir = integrator.tdir > zero(integrator.tdir) ?
          integrator.t > integrator.sol.prob.tspan[2] - difference ? -true :  true :
          integrator.t < integrator.sol.prob.tspan[2] + difference ?  true : -true
end