text:
# Returns the length (width) of the piece of smallest length (width) that fits
# the plate defined by LxW. Errors if there is no piece that fits the plate.
function _min_l_of_fit_piece(
	sllw :: SortedLinkedLW{D, S}, L :: S, W :: S
) :: S where {D, S}
	for (i, l) in enumerate(sllw.sl)
		l > L && break
		sllw.w[sllw.sli2pii[i]] <= W && return l
	end
	@error "No piece fits the plate"
end
function _min_w_of_fit_piece(
	sllw :: SortedLinkedLW{D, S}, L :: S, W :: S
) :: S where {D, S}
	for (i, w) in enumerate(sllw.sw)
		w > W && break
		sllw.l[sllw.swi2pii[i]] <= L && return w
	end
	@error "No piece fits the plate"
end

# Returns a vector with length equal to highest_plate, and in which position
# `p` gives the range of `cuts` that have `p` as their first field (i.e., the
# parent plate). This range is shifted by offset.
# NOTE: `highest_plate >= maximum(getindex.(cuts, 1))` or undefined behavior.
# NOTE: if some plate `p` (1 < p <= highest_plate) is not the first field of
# any element in `cuts`, then the range of `p` in the returned vector will be
# an empty range and, consequently, the returned array will probably fail a
# `issorted` check, because no empty range is considered greater than a
# non-empty range;
function _build_pp_ranges(
	cuts :: AbstractVector{NTuple{3, P}},
	highest_plate :: P,
	offset :: P = zero(P)
) :: Vector{UnitRange{P}} where {P}
	# If `cuts` is empty, return a vector of empty ranges of the expected size.
	empty_r = one(P):zero(P)
	isempty(cuts) && return fill(empty_r, highest_plate)
	# If there is a cut, there should be at least one plate.
	@assert highest_plate >= one(P)
	ranges = Vector{UnitRange{P}}(undef, highest_plate)
	curr_pp = first(cuts)[PARENT]
	# The original plate (index one) is always present.
	@assert isone(curr_pp)
	first_idx_curr_pp = one(P)
	# The first field of `cuts` elements (i.e., PARENT) never decreases, and when
	# it increases, the range in which they kept the same value is saved.
	for (i, cut) in pairs(cuts)
		pp = cut[PARENT]
		if pp < curr_pp
			@error "The parameter `cuts` must be sorted by its first field."
		elseif pp > curr_pp
			ranges[curr_pp] = (offset + first_idx_curr_pp):(offset + i - 1)
			# If pp - curr_pp > 1 (i.e., pp values were skipped) save empty
			# ranges for such values.
			for j = (curr_pp+1):(pp-1); ranges[j] = empty_r; end
			# Update our sentinels.
			curr_pp = pp
			first_idx_curr_pp = i
		end
		# If `curr_pp` is equal to `highest_plate` we assume `cuts` is
		# non-decreasing and break early. The `ranges[highest_plate]` is
		# saved after the loop.
		curr_pp == highest_plate && break
	end
	# The last `curr_pp` has no change of value to trigger its writting.
	ranges[curr_pp] = (offset + first_idx_curr_pp):(offset + length(cuts))
	return ranges
end

function _get_LW(
	bp :: ByproductPPG2KP{D, S, P},
	plate_idx :: P
) :: Tuple{S, S} where {D, S, P}
	return bp.pli_lwb[plate_idx][LENGTH], bp.pli_lwb[plate_idx][WIDTH]
end

# Find a cut in `bp.cuts[cuts_by_pp[pp]]` with a FIRST_CHILD of size `cut_size`
# in the dimension `dim`, `push!` it to `cuts_done`, and return its
# `FIRST_CHILD` and `SECOND_CHILD` (a convenience, as we could use the last
# element of `cuts_done` to index `bp.cuts` and get them).
function _make_cut!(
	cuts_done :: Vector{P},
	pp :: P,
	bp :: ByproductPPG2KP{D, S, P},
	cuts_by_pp :: Vector{UnitRange{P}},
	cut_size :: S,
	dim :: Dimension
) where {D, S, P}
	idx_in_view = findfirst(
		c -> bp.pli_lwb[c[FIRST_CHILD]][dim] == cut_size,
		view(bp.cuts, cuts_by_pp[pp])
	) :: Int
	cut_idx = cuts_by_pp[pp].start + convert(P, idx_in_view) - one(P)
	push!(cuts_done, cut_idx)
	@assert bp.cuts[cut_idx][PARENT] == pp
	return bp.cuts[cut_idx][FIRST_CHILD], bp.cuts[cut_idx][SECOND_CHILD]
end

# Takes a `plate_idx`, and then successively chip the plate removing trims of
# `trim_size` (of dimension `dim`) while the plate is not smaller than
# `size_limit`. Returns a plate with `dim` smaller than `size_limit +
# trim_size`. Save the `trim_cuts` to `cuts_done`. The number of trims is
# not computed beforehand because of `--PPG2KP-round2disc` (i.e., a trim
# can reduce the dimension by more than `trim_size`).
function _successive_trims!(
	cuts_done :: Vector{P},
	plate_idx :: P,
	bp :: ByproductPPG2KP{D, S, P},
	cuts_by_pp :: Vector{UnitRange{P}},
	size_limit :: S,
	trim_size :: S,
	dim :: Dimension
) :: P where {D, S, P}
	pp = plate_idx
	safe_limit = size_limit + trim_size
	while bp.pli_lwb[pp][dim] >= safe_limit
		# We throw away the info about which was the first child (i.e., the trim
		# plate) and update the parent plate to be the second child.
		_, pp = _make_cut!(cuts_done, pp, bp, cuts_by_pp, trim_size, dim)
		@assert !iszero(pp)
	end
	@assert bp.pli_lwb[pp][dim] >= size_limit
	return pp
end

# Trim the dimension `dim` of plate `pp` to size `s` (if `bm == FURINI`) or
# close enough (if `bm == BECKER`, and by 'close enough' we mean that no
# piece would fit in the trim of the 'close enough' plate to exactly size `s`).
# The cuts needed are added to `cuts_done`. For a plate index `pli` the
# `cuts_by_pp[pli]` has the range of cuts (that divide `dim`) with `pli` as
# the PARENT plate of the cut.
function _safe_trim_dim!(
	cuts_done :: Vector{P},
	pp :: P,
	s :: S,
	dim :: Dimension,
	cuts_by_pp :: Vector{UnitRange{P}},
	sllw :: SortedLinkedLW{D, S},
	bp :: ByproductPPG2KP{D, S, P},
	bm :: BaseModel
) where {D, S, P}
	# rS: size of `dim` in plate `pp`
	rS = bp.pli_lwb[pp][dim]
	@assert rS >= s
	if s != rS # Check if trimming is actually necessary.
		# If the cut is in the first half of the plate it is guaranteed to exist.
		if s <= div(rS, 2)
			# In the case of a 'clean' trim cut, the trim is the SECOND_CHILD that
			# is thrown away, and the PARENT is updated to the FIRST_CHILD.
			pp, _ = _make_cut!(cuts_done, pp, bp, cuts_by_pp, s, dim)
		else # If a trim cut may be needed.
			L, W = _get_LW(bp, pp)
			if dim == LENGTH
				min_s = _min_l_of_fit_piece(sllw, L, W)
			else
				min_s = _min_w_of_fit_piece(sllw, L, W)
			end
			# If the difference between the plate and the piece allows for copies of
			# the smallest-dim-size piece-that-fits, reduce the difference by
			# successive trim cuts.
			if rS - s >= min_s
				pp = _successive_trims!(
					cuts_done, pp, bp, cuts_by_pp, s, min_s, dim
				)
				rS = bp.pli_lwb[pp][dim]
				# If bm == BECKER and the other dim has a difference smaller than
				# the smallest-other-dim-size piece-that-fits, then it is guaranteed
				# that the extraction will be in bp.np.
			end
			# Otherwise, if bm == FURINI, the match of piece and plate need
			# to be exact, so we may need a final trim cut. It was done here,
			# after the successive trim cuts, because here it is guaranteed to
			# exist. If it was done before, it could be the cut did not exist
			# because the symmetry removal used by Furini2016 (i.e, if there was a
			# cut of length exactly L - l, in the first half of the plate, then a
			# cut of length l, in the second half of the plate, would not exist).
			if bm == FURINI && s != rS
				pp, _ = _make_cut!(cuts_done, pp, bp, cuts_by_pp, s, dim)
			end
		end
	end
	return pp
end

function _sell_piece!(
	pieces_sold :: Vector{P},
	cuts_done :: Vector{P},
	plate_idx :: P,
	piece_idx :: D,
	sllw :: SortedLinkedLW{D, S},
	bp :: ByproductPPG2KP{D, S, P},
	hcuts_by_pp :: Vector{UnitRange{P}},
	vcuts_by_pp :: Vector{UnitRange{P}},
	bm :: BaseModel
) :: Nothing where {D, S, P}
	# The piece-to-be-sold dimensions.
	piece_l, piece_w = bp.l[piece_idx], bp.w[piece_idx]
	# pp: the current parent plate that is updated as we trim the given plate
	pp = plate_idx
	pp = _safe_trim_dim!(
		cuts_done, pp, piece_l, LENGTH, hcuts_by_pp, sllw, bp, bm
	)
	pp = _safe_trim_dim!(
		cuts_done, pp, piece_w, WIDTH, vcuts_by_pp, sllw, bp, bm
	)
	# Now a extraction from the current plate to piece_idx MUST exist.
	np_idx = findfirst(==((pp, piece_idx)), bp.np) :: Int
	push!(pieces_sold, convert(P, np_idx))
	return
end

"""
    cuts_and_extractions_from_2_staged_solution(pattern, bp, bm)

Given a 2-staged `pattern` and a PPG2KP model `bp` (byproduct), returns the
cuts (`bp.cuts` indexes) and extractions (`bp.np` indexes) that make up a
valid solution representing such `pattern`. Note: this is intended to be
called over non-priced PPG2KP models; if the model is priced the variables
assumed to exist in this code may not exit anymore (especially if Furini's
multistep pricing is used).

The `pattern` parameter has many restrictions:

1) It is 2-staged, what is enforced by being a vector of vectors.
2) Each inner vector is a vertical stripe with the width of the first piece
   inside such inner vector. No horizontal stripes allowed. No piece inside
   an inner vector has width larger than the first piece of the inner vector.
3) The first inner vector has the largest stripe width.
4) No inner vectors should be empty. Every number is a valid piece type index.
"""
function cuts_and_extractions_from_2_staged_solution(
	pattern :: Vector{Vector{D}},
	bp :: ByproductPPG2KP{D, S, P},
	bm :: BaseModel
) :: Tuple{Vector{P}, Vector{P}} where {D, S, P}
	cuts_done, pieces_sold = P[], P[]
	isempty(pattern) && return cuts_done, pieces_sold
	@assert all((!isempty).(pattern)) # No stripe should be empty.
	# @assert explanation: the first stripe is the one of largest width.
	@assert all(fp -> bp.w[fp] <= bp.w[pattern[1][1]], getindex.(pattern, 1))
	# @assert explanation: each of the inner vectors has no piece with
	# width greater than their width (i.e., the width of the first element).
	@assert all(a -> all(p -> bp.w[p] <= bp.w[a[1]], a), pattern)
	sllw = SortedLinkedLW(D, bp.l, bp.w)
	highest_plate = length(bp.pli_lwb)
	# Compute cut ranges that share the same orientation (vertical or horizontal)
	# and the same parent plate (pp). This is O(n) and greatly reduces the effort
	# spent finding the desired cuts in the rest of the algorithm.
	hcuts_by_pp = _build_pp_ranges(
		(@view bp.cuts[one(P):(bp.first_vertical_cut_idx - one(P))]), highest_plate
	)
	num_cuts = length(bp.cuts)
	vcuts_by_pp = _build_pp_ranges(
		(@view bp.cuts[bp.first_vertical_cut_idx:num_cuts]), highest_plate,
		bp.first_vertical_cut_idx - one(P)
	)
	rp = one(P) # The initial `r`emaining `p`late is the original plate.
	num_stripes = length(pattern)
	# stripe_plates: store the plate for each stripe, used after but not updated.
	stripe_plates = Vector{P}(undef, num_stripes)
	# The order we make the cuts is not specially relevant, except that we only
	# have the guarantee that the cuts really exist if they are in the first
	# half of the plate, and this happens naturally if we left the largest
	# width stripe for last.
	for vstripe_idx in 2:num_stripes
		# Cut a vertical (i.e., keep the length divides the width) stripe.
		stripe_plates[vstripe_idx], rp = _make_cut!(
			cuts_done, rp, bp, vcuts_by_pp, bp.w[pattern[vstripe_idx][1]], WIDTH
		)
	end
	# If it exists a clean cut for the last stripe do it.
	largest_w = bp.w[pattern[1][1]]
	if largest_w <= (bp.pli_lwb[rp][WIDTH]+1)/2
		stripe_plates[1], _ = _make_cut!(
			cuts_done, rp, bp, vcuts_by_pp, largest_w, WIDTH
		)
	else # Otherwise, just use the remaining plate as last stripe.
		stripe_plates[1] = rp
	end
	# Now all the vertical stripes are cut. The stripes need to be cut
	# horizontally creating a new plate for each piece inside the stripe,
	# and then these plates may need a vertical trim (because they are smaller
	# than the stripe width) before they are finally sold as pieces.
	rp = zero(P) # From there on, only the `stripe_rp` is relevant.
	for vstripe_idx in 1:num_stripes # the outer order is not relevant
		stripe_rp = stripe_plates[vstripe_idx]
		vstripe = pattern[vstripe_idx]
		_, idx_largest_l_piece_in_vstripe = findmax(bp.l[vstripe])
		# The inner order is not relevant, except the piece with the largest
		# length should be specially treated after. For the same reason we
		# treat the stripe of largest width specially above.
		for (vstripe_idx, piece_idx) in pairs(vstripe)
			vstripe_idx == idx_largest_l_piece_in_vstripe && continue
			plate_idx, stripe_rp = _make_cut!(
				cuts_done, stripe_rp, bp, hcuts_by_pp, bp.l[piece_idx], LENGTH
			)
			_sell_piece!(
				pieces_sold, cuts_done, plate_idx, piece_idx, sllw, bp,
				hcuts_by_pp, vcuts_by_pp, bm
			)
		end
		# The largest length piece of the largest width stripe can now just be
		# sold. It is cut after the loop because then there is guarantee that all
		# `_make_cut!` calls inside the loop happen in the first half of the
		# current stripe_rp and, therefore, are safe to make. For this last
		# piece a `_make_cut!` call is not necessary and `_sell_piece!` will
		# take care of the edge cases.
		piece_idx = vstripe[idx_largest_l_piece_in_vstripe]
		_sell_piece!(
			pieces_sold, cuts_done, stripe_rp, piece_idx, sllw, bp,
			hcuts_by_pp, vcuts_by_pp, bm
		)
	end

	return cuts_done, pieces_sold
end
export cuts_and_extractions_from_2_staged_solution

"""
    save_mip_start(model)

TODO: document. Returns the four last arguments of raw_mip_start!
"""
@timeit TIMER function save_mip_start(model)
	nz_pe_idxs, nz_pe_vals = gather_nonzero(model[:picuts], Int)
	nz_cm_idxs, nz_cm_vals = gather_nonzero(model[:cuts_made], Int)
	return (nz_pe_idxs, nz_pe_vals, nz_cm_idxs, nz_cm_vals)
end

"""
    unset_mip_start!(model, nz_pe_idxs, nz_cm_idxs)

TODO: document. Just take the first and third returned values from
`save_mip_start!` (or the first and third parameters of `raw_mip_start!`
because the specific values of the old start are irrelevant).
"""
@timeit TIMER function unset_mip_start!(model, nz_pe_idxs, nz_cm_idxs)
	pe_nothings = Iterators.repeated(nothing, length(nz_pe_idxs))
	cm_nothings = Iterators.repeated(nothing, length(nz_cm_idxs))
	raw_mip_start!(model, nz_pe_idxs, pe_nothings, nz_cm_idxs, cm_nothings)
	return
end

"""
    raw_mip_start!(model, nz_pe_idxs, nz_pe_vals, nz_cm_idxs, nz_cm_vals)

TODO: document.
"""
@timeit TIMER function raw_mip_start!(
	model, nz_pe_idxs, nz_pe_vals, nz_cm_idxs, nz_cm_vals
)
	@assert length(nz_pe_idxs) == length(nz_pe_vals)
	@assert length(nz_cm_idxs) == length(nz_cm_vals)
	pe = model[:picuts]
	cm = model[:cuts_made]
	# The MOI.set is used instead of JuMP.set_start_value because the last do
	# not accept nothing as a value (to unset the MIP start).
	MOI.set.(
		backend(model), MOI.VariablePrimalStart(), JuMP.index.(pe[nz_pe_idxs]),
		((eltype(nz_pe_vals) <: Number) ? Float64.(nz_pe_vals) : nz_pe_vals)
	)
	MOI.set.(
		backend(model), MOI.VariablePrimalStart(), JuMP.index.(cm[nz_cm_idxs]),
		((eltype(nz_cm_vals) <: Number) ? Float64.(nz_cm_vals) : nz_cm_vals)
	)
	#set_start_value.(pe[nz_pe_idxs], nz_pe_vals)
	#set_start_value.(cm[nz_cm_idxs], nz_cm_vals)
	return
end
export raw_mip_start!

"""
    mip_start_by_heuristic!(model, seed, p, bp, bm)

TODO: document. Returns two tuples, the first is the result of internal
`fast_iterated_greedy` the second tuple contain the last four arguments
of the raw_mip_start! used to MIP-start the model (this way
unset_mip_start! can be easily called over the model).
"""
function mip_start_by_heuristic!(
	model, bp :: ByproductPPG2KP{D, S, P}, p, seed, bm :: BaseModel
) where {D, S, P}
	rng = Xoroshiro128Plus(seed)
	bkv, selected, shelves = fast_iterated_greedy(
		bp.d, p, bp.l, bp.w, bp.L, bp.W, rng
	)
	raw_ws = @timeit TIMER "create_raw_ws_from_heuristic" begin
		cuts, extractions = cuts_and_extractions_from_2_staged_solution(
			shelves, bp, bm
		)
		qt_cuts = unify!(D, cuts)
		qt_extractions = unify!(D, extractions)
		(extractions, qt_extractions, cuts, qt_cuts)
	end
	raw_mip_start!(model, raw_ws...)
	return (bkv, selected, shelves), raw_ws
end
export mip_start_by_heuristic!


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Text cell with ANSI escape sequences.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export AnsiTextCell

"""
    AnsiTextCell

A text cell that supports rendering ANSI escape sequences without interfering
with the table layout.

# Fields

**Public**

- `string::String`: The string with the cell text that can contain ANSI escape
    sequences.

**Private**

- `_rendered_lines::Union{Nothing, Vector{String}}`: The lines with the rendered
    strings.
- `_stripped_lines::Union{Nothing, Vector{String}}`: The lines with the
    printable text.
- `_crops::Dict{Int, Int}`: Dictionary with the number of characters that must
    be cropped.
- `_left_pads::Dict{Int, Int}`: Left padding to be applied to each line.
- `_right_pads::Dict{Int, Int}`: Right padding to be applied to each line.
- `_suffixes::Dict{Int, String}`: Suffixed to be applied to each line.
"""
mutable struct AnsiTextCell <: CustomTextCell
    string::String

    # Private
    _rendered_lines::Union{Nothing, Vector{String}}
    _stripped_lines::Union{Nothing, Vector{String}}
    _crops::Union{Nothing, Vector{Int}}
    _left_pads::Union{Nothing, Vector{Int}}
    _right_pads::Union{Nothing, Vector{Int}}
    _suffixes::Union{Nothing, Vector{String}}
end

"""
    AnsiTextCell(string::AbstractString)

Create an [`AnsiTextCell`](@ref) using `string`.

    AnsiTextCell(renderfn[; context])

Create an [`AnsiTextCell`](@ref) using a render function.

`renderfn` is a function with the following signature:

    renderfn(io)::String

that renders a string that can contain ANSI sequences into `io`.

`context` is a tuple of context arguments passed to an `IOContext` that
`renderfn` receives. See [`IOContext`](@ref) for details on what arguments are
available.

Useful for supporting packages that have rich terminal outputs.

## Examples

Below are examples for wrappers around `AnsiTextCell` to print rich data into
tables that make use of packages with rich terminal output.

### Crayons.jl

Apply custom decoration to text inside a cell.

```julia
using Crayons, PrettyTables

b = crayon"blue bold"
y = crayon"yellow bold"
g = crayon"green bold"

pretty_table([AnsiTextCell("\$(g)This \$(y)is \$(b)awesome!") for _ in 1:5, _ in 1:5])
```

### ImageInTerminal.jl

Show images inside a table.

```julia
using ImageInTerminal, PrettyTables

function ImageCell(img, size)
    return AnsiTextCell(
        io -> ImageInTerminal.imshow(io, img),
        context = (:displaysize => size,),)
end

using TestImages
img = testimage("lighthouse")
pretty_table([ImageCell(img, (20, 20)) ImageCell(img, (40, 40))])
```

### UnicodePlots.jl

Show a variety of plots in a table.

```julia
using UnicodePlots, PrettyTables

function UnicodePlotCell(p)
    return AnsiTextCell(
        io -> show(io, p),
        context = (:color => true,)
    )
end

pretty_table([
    UnicodePlotCell(barplot(Dict("x" => 10, "y" => 20)))
    UnicodePlotCell(boxplot([1, 3, 3, 4, 6, 10]))
])
```

### CommonMark.jl

Use rich Markdown inside tables.

```julia
using CommonMark, PrettyTables

function MarkdownCell(md)
    return AnsiTextCell(
        renderfn = io -> display(TextDisplay(io), md),
        context = (:color => true,)
    )
end

pretty_table([MarkdownCell(cm"**Hi**") MarkdownCell(cm"> quote")])
```
"""
function AnsiTextCell(string::AbstractString)
    return AnsiTextCell(
        string,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing
    )
end

function AnsiTextCell(renderfn::Function; context::Tuple = ())
    # Render the text and create the text cell.
    io = IOBuffer()
    renderfn(IOContext(io, context...))
    rendered = String(take!(io))
    return AnsiTextCell(rendered)
end

################################################################################
#                                     API
################################################################################

function get_printable_cell_line(cell::AnsiTextCell, l::Int)
    if l > length(cell._stripped_lines)
        return ""
    else
        lpad = cell._left_pads[l]
        rpad = cell._right_pads[l]
        line = " "^lpad * cell._stripped_lines[l] * " "^rpad

        # Compute the total size of the string.
        line_width = lpad + textwidth(cell._stripped_lines[l]) + rpad

        # Return the cropped string.
        return _crop_str(line, line_width - cell._crops[l], line_width)
    end
end

function get_rendered_line(cell::AnsiTextCell, l::Int)
    if l > length(cell._stripped_lines)
        return ""
    else
        lpad = cell._left_pads[l]
        rpad = cell._right_pads[l]
        suffix = cell._suffixes[l]

        # Compute the total size of the string.
        line_width = lpad + textwidth(cell._stripped_lines[l]) + rpad

        # Create the rendered line.
        line = " "^lpad * cell._rendered_lines[l] * " "^rpad

        # Crop the rendered line.
        cropped_line = _crop_str(line, line_width - cell._crops[l], line_width)

        # We must reset everything after rendering the cell to avoid messing
        # with the decoration of the table.
        return cropped_line * "\e[0m" * suffix
    end
end

function append_suffix_to_line!(cell::AnsiTextCell, l::Int, suffix::String)
    cell._suffixes[l] *= suffix
    return nothing
end

function apply_line_padding!(cell::AnsiTextCell, l::Int, left_pad::Int, right_pad::Int)
    cell._left_pads[l] += left_pad
    cell._right_pads[l] += right_pad
    return nothing
end

function crop_line!(cell::AnsiTextCell, l::Int, num::Int)
    if num ≥ 0
        cell._crops[l] = num
    end

    return nothing
end

function parse_cell_text(cell::AnsiTextCell; kwargs...)
    # Regex to remove escape sequences from the string.
    r_ansi_escape = r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])"

    _rendered_lines = map(String, split(cell.string, '\n'))
    _stripped_lines = map(l -> replace(l, r_ansi_escape => ""), _rendered_lines)

    num_lines = length(_rendered_lines)

    cell._rendered_lines = _rendered_lines
    cell._stripped_lines = _stripped_lines
    cell._crops          = zeros(Int, num_lines)
    cell._left_pads      = zeros(Int, num_lines)
    cell._right_pads     = zeros(Int, num_lines)
    cell._suffixes       = fill("", num_lines)

    return cell._stripped_lines
end

function reset!(cell::AnsiTextCell)
    cell._rendered_lines = nothing
    cell._stripped_lines = nothing
    cell._crops          = nothing
    cell._left_pads      = nothing
    cell._right_pads     = nothing
    cell._suffixes       = nothing

    return nothing
end

ext: jl


text:
module mocapviz

using Quaternions
using MeshCat
using CoordinateTransformations: LinearMap, AffineMap, Translation
using Rotations: rotation_between
using GeometryBasics: Cylinder, HyperSphere, Point
using LinearAlgebra: norm, diagm
using DSP   # conv
using Colors   # RGB, RGBA, colorant_str
using Formatting, ArgCheck

import MeshCat: GenericMaterial

Base.convert(::Type{Tuple}, c::RGB) = (c.r, c.g, c.b)
transparency(c::RGB, α) = RGBA(convert(Tuple, c)..., α)
darkencol(x::RGB) = whitebalance(x, colorant"white", colorant"grey55")
darkencol(x::RGBA) = let u=convert(RGB, x); y=darkencol(u); y=RGBA(y.r, y.g, y.b, x.alpha); y; end
darkencol(x::GenericMaterial) = GenericMaterial(_type=x._type, color=darkencol(x.color))   # includes MeshPhongMaterial

const blackmesh = MeshPhongMaterial(color=RGBA(0.1, 0.1, 0.1, 0.4))
const greymesh = MeshPhongMaterial(color=RGBA(0.3, 0.3, 0.3, 0.7))
const yellowmesh = MeshPhongMaterial(color=RGBA(204/255, 204/255, 0., 0.7))
const redmesh = MeshPhongMaterial(color=RGBA(224/255, 131/255, 94/255, 0.7))

const bluemesh = MeshPhongMaterial(color=transparency(colorant"#4C9DFF", 0.8));
const greenmesh = MeshPhongMaterial(color=transparency(colorant"#4CE54C", 0.8));
const deepredmesh = MeshPhongMaterial(color=transparency(colorant"#E56666", 0.8));

# note that the y-axis (vertical) in three.js is the *last*/*3rd* dimension.
project_ground!(x::AbstractVector) = begin; (length(x)==0 && return); @assert length(x) == 3; x[3] = 0; end;
project_ground(x::AbstractVector) = begin; x = copy(x); project_ground!(x); x; end
project_ground!(X::AbstractMatrix) = begin; (length(X)==0 && return); @assert size(X,2) == 3; X[:,3] .= 0; end;
project_ground(X::AbstractMatrix) = begin; X = copy(X); project_ground!(X); X; end

# initialise lines and dots
mutable struct LineCollection{T, R}
    data::Array{T}
    vnames::Array{String}
    attributes::Array{R}
end
mutable struct JointCollection{T, R}
    data::Array{T}
    vnames::Array{String}
    attributes::Array{R}
end

Base.length(x::LineCollection) = length(x.data)
Base.length(x::JointCollection) = length(x.data)

function update_pos!(lines::LineCollection, origin::Matrix, extremity::Matrix)
    @argcheck length(lines) == size(origin, 1) == size(extremity, 1)
    lines.data = [Cylinder(Point{3}(origin[i,:]), Point{3}(extremity[i,:]), l.r)
                  for (i,l) in enumerate(lines.data)]
end

function update_pos!(lines::LineCollection, data::Matrix, parents::Vector)
    parents = parents[1:length(lines)]
    parent = data[parents,:]
    data = data[2:end, :]  # remove root (has no line)
    update_pos!(lines, parent, data)
end

function update_pos!(joints::JointCollection, data::Matrix)
    @argcheck length(joints) == size(data, 1)
    joints.data = [HyperSphere(Point{3}(data[i,:]), l.r) for (i,l) in enumerate(joints.data)]
end

# 3D plotting utilities: map line from x₁ → x₂ into affine transformation of e₁.
function line_into_transform(x1, x2)
    scaling, position, rotation = line_into_transforms_indiv(x1, x2)
    return AffineMap(rotation, position) ∘ LinearMap(diagm(0=>scaling))
end

function line_into_transforms_indiv(x1, x2)
    scaling = [norm(x2-x1), 1., 1]
    position = x1
    rotation = rotation_between([1., 0, 0], x2 - x1)
    return scaling, position, rotation
end

line_into_transform(x::Cylinder) = line_into_transform(x.origin, x.extremity)
line_into_transforms_indiv(x::Cylinder) = line_into_transforms_indiv(x.origin, x.extremity)

function setobj_collection!(v::AbstractVisualizer, args...)
    for arg in args
        setobj_collection!(v, arg)
    end
end

function setobj_collection!(v::AbstractVisualizer, coll::Union{LineCollection,JointCollection})
    for i in 1:length(coll)
        setobject!(v[coll.vnames[i]], coll.data[i], coll.attributes[i])
    end
end

function settransform_collection!(v::AbstractVisualizer, args...; anim::Bool=false)
    for arg in args
        settransform_collection!(v, arg; anim=anim)
    end
end

function settransform_collection!(v::AbstractVisualizer, L::LineCollection; anim::Bool=false)
    if false #anim
        # for i in 1:length(L)
        #     scaling, position, rotation = line_into_transforms_indiv(L.data[i])
        #     anim_settransform!(v[L.vnames[i]], scaling, position, rotation)
        # end
    else
        for i in 1:length(L)
            tform = line_into_transform(L.data[i])
            settransform!(v[L.vnames[i]], tform)
        end
    end
end

function settransform_collection!(v::AbstractVisualizer, J::JointCollection; anim::Bool=false)
    for i in 1:length(J)
        settransform!(v[J.vnames[i]], Translation(J.data[i].center))
    end
end

# # HACK: alternative method for rotation∘scale because MeshCat *animation* lowering is broken
# # -> This was used for a previous version of MeshCat. I don't think it's needed anymore, HOWEVER, the angles
# # between joints is a little off. I don't have time to chase this down.
# function anim_settransform!(vis::MeshCat.AbstractVisualizer, scaling, position, rotation)
#     clip = MeshCat.getclip!(vis)
#     MeshCat._setprop!(clip, vis.current_frame, "scale", "vector3", scaling)
#     MeshCat._setprop!(clip, vis.current_frame, "position", "vector3", position)
#     MeshCat._setprop!(clip, vis.current_frame, "quaternion", "quaternion", MeshCat.js_quaternion(rotation))
# end

function create_animation(data::Vector, names::Union{String, Array{String}}="dataset";
    vis=nothing, parents=[1,2,3,4,1,6,7,8,1,10,11,12,12,14,15,16,12,18,19,20],
    jointmesh::Union{AbstractMaterial, Vector{T} where T <: AbstractMaterial}=greymesh,
    linemesh::Union{AbstractMaterial, Vector{T} where T <: AbstractMaterial}=yellowmesh, scale=0.1,
    camera::Symbol=:front, path::Union{Nothing, AbstractMatrix{T} where T <: Real, Vector}=nothing)

    # dimensions
    Ts = [size(d,1) for d in data]
    ls = [size(d,2) for d in data]
    n  = length(ls)
    maxT = maximum(Ts)
    ispath = !(path === nothing)

    # convenience expansions if args are singular
    (names isa String) && (names = [names * string(i) for i in 1:n])
    (jointmesh isa AbstractMaterial) && (jointmesh = repeat([jointmesh], n))
    (linemesh isa AbstractMaterial) && (linemesh = repeat([linemesh], n))
    ispath && (path isa AbstractMatrix) && (path = [path])
    npaths = ispath ? length(path) : 0
    (npaths >= 1) && (pathmesh = map(darkencol, linemesh));
    (npaths == 1) && (pathmesh = [blackmesh]);
    # assertions
    @argcheck length(data) == length(names)
    @assert (all([ndims(d) for d in data] .== 3) && all([size(d,3) for d in data] .==3)) "Need NxJx3 matrices."
    @argcheck camera in [:front, :back, :static]
    ispath && @argcheck size(path[1]) == (Ts[1], 12)
    # @assert all(ls .<= length(parents)-1) "more joints given (dim 1) than are assigned parents."
    # any(ls .> length(parents)) && @warn format("({:d}/{:d}) datasets have fewer joints than specified in parents.",
    #     sum(ls .> length(parents)), n)

    if vis === nothing
        vis = Visualizer()
        open(vis)
    else
        @assert (vis isa Visualizer) "Don't know what to do with vis::$(typeof(vis)). Expecting Visualizer."
    end

    gen_numstr(prefix) = num -> prefix * format("{:02d}", num);
    dotstr, linestr, pathstr = gen_numstr("dot_"), gen_numstr("line_"), gen_numstr("path_")

    # Initialise objects for each dataset
    objs = map(1:n) do i
        if i <= npaths
            obj_path = LineCollection(repeat([Cylinder(zero(Point{3, Float64}), Point{3}([1.,0,0]), 0.02)], 10), map(pathstr,  1:10), repeat([pathmesh[i]], 10))
            setobj_collection!(vis[names[i]], obj_path)
        else
            obj_path = nothing
        end

        nj, nl = ls[i], ls[i]-1
        obj_lc = LineCollection(repeat([Cylinder(zero(Point{3, Float64}), Point{3}([1.,0,0]), 0.03)], nl), map(linestr, 1:nl), repeat([linemesh[i]], nl))
        obj_jc = JointCollection(repeat([HyperSphere(Point{3}(0.), 0.04)], nj), map(dotstr, 1:nj), repeat([jointmesh[i]], nj))
        setobj_collection!(vis[names[i]], obj_lc, obj_jc)  # note order: lines at back, joints at front
        (obj_lc, obj_jc, obj_path)
    end

    #=========== BEGIN ANIMATION CONSTRUCTION ===============================#
    anim = Animation()

    # camera trajectory (smoothed (MA-4) over skeleton traj)
    cam_ix = 1
    pos = hcat([data[cam_ix][tt,1,1]*scale for tt in 1:maxT], [data[cam_ix][tt,1,3]*scale for tt in 1:maxT])
    pos = hcat(conv(pos[:,1], Windows.rect(4)/4), conv(pos[:,2], Windows.rect(4)/4))[1:maxT,:]
    cc = [scale, scale] * (camera == :front ? -1 : +1)

    reord = [1,3,2] # need to permute z/y axis for three.js setup
    for tt in 1:maxT
        atframe(anim, tt-1) do
            for i in 1:n
                if Ts[i] >= tt
                    lines, joints, pathlines = objs[i]
                    update_pos!(lines, data[i][tt,:,reord]*scale, parents)   # Note that update_pos only updates the line/joint collection position.
                    update_pos!(joints, data[i][tt,:,reord]*scale);          # Once these objects are correct, `settransform_collection` pushes => animation.
                    settransform_collection!(vis[names[i]], lines, joints; anim=true)  # note order: lines at back, joints at front

                    if i <= npaths
                        trail = randn(5,3)*0.001   # case i=1: all segments need positive length (Meshcat reqmt.)
                        trail_ixs = max(tt-41, 1):10:(tt-1)
                        if false
                            trail[(6-length(trail_ixs)):5, :] = project_ground(data[i][trail_ixs, 1, reord])
                        else
                            trail[(6-length(trail_ixs)):5, :] = hcat(path[i][trail_ixs, 1], path[i][trail_ixs, 7], zeros(length(trail_ixs)))
                        end
                        path_matrix = vcat(trail, hcat(reshape(path[i][tt, :], 6, 2), zeros(6)))*scale
                        update_pos!(pathlines, path_matrix[1:10,:], path_matrix[2:11,:])
                        settransform_collection!(vis[names[i]], pathlines; anim=true)
                    end
                end
            end
            # one-per-frame updates
            if camera != :static
                settransform!(vis["/Cameras/default"], Translation(pos[tt,1] + cc[1], pos[tt,2] + cc[2], 1))
            end
        end
    end

    setanimation!(vis, anim)
    return vis
end

function plot_2d_skeleton(joints::AbstractMatrix)
    @argcheck size(joints) == (21, 3)
    parents=[1,2,3,4,1,6,7,8,1,10,11,12,12,14,15,16,12,18,19,20]
    scatter(joints[:,1], joints[:,2]);
    [plot([[joints[p,j], joints[e+1,j]] for j in 1:2]..., c="k") for (e,p) in enumerate(parents)]
    [text(joints[i,1]+0.2, joints[i,2]-0.7, i) for i in 1:21]
    gca().set_aspect("equal")
    gca().axis("off");
end

function plot_2d_skeleton(root_y::AbstractFloat, joints::AbstractMatrix)
    @argcheck size(joints) == (20, 3)
    joints = vcat([0 root_y 0], joints);
    plot_2d_skeleton(joints)
end

end

max_stars_repo_name: ornithos/pytorch-mtds-mocap

max_issues_repo_name: ornithos/pytorch-mtds-mocap

max_forks_repo_name: ornithos/pytorch-mtds-mocap


text:
abstract type Wavelet end

kernelmatrix(w, n::Int, m::Int = n) = map(i -> w(abs(i[1]-i[2])), CartesianIndices((n,m)))

function kernelmatrix(v::Vector, n::Int)
    m = []
    idx = ceil(Int, length(v)/2)
    for i in 1:n
        if i <= idx
            push!(m, vcat(v[idx-i+1:end], zeros(n-length(v[idx-i+1:end]))))
        elseif i <= n - idx
            vs = vcat(zeros(i - idx), v)
            push!(m, vcat(vs, zeros(n - length(vs))))
        else
            push!(m, vcat(zeros(i - idx), v[1:n-i+idx]))
        end
    end

    Matrix(hcat(m...)')
end

struct Ricker <: Wavelet
  f::Real
  dt::Real
end

function (w::Ricker)(h::Int)
  t = w.dt * h

  x = (π * w.f * t) .^ 2
  A = (1 .- 2x) .* exp.(-x)
end

struct Ormsby <: Wavelet
  f::NTuple{4,Real}
  dt::Real
end

function (w::Ormsby)(h::Int)
  t = w.dt * h
  f = collect(w.f)

  z = π * f
  s = @. sinc(z * t)^2
  A = (s[4] * z[4]^2 - s[3] * z[3]^2)/(z[4]-z[3]) - (s[2] * z[2]^2 - s[1] * z[1]^2)/(z[2]-z[1])
end


text:
include("yeast_dyn.jl")
include("shuffle.jl")
include("coop.jl")
srand(1)

## listt = getpointattr(interc,selfact)
println("\n yeast's avg coop = ",cooppointattr(interc,selfact)[1])
## println("\n yeast's avg coop = ",coopbasinweighted(interc,selfact)[1]/length(selfact))

## --- a (wrong?) analytical model for mean cooperativity
## --- gamma = 2.5 is the best fit to the data for n=100
## --- with basin weighted averaging - coopbasinweighted()
function coop_anal(p,a)
return p^a + (1-p)^a
end

# --- calculate average cooperativity over random networks
tstart = time()
n = 11*3
selfact = rand(0:1,n)
n = length(selfact)
ne = 29*3                                # number of edges
ntotal = 100;                         # ensemble size
output = Array(Float64,0,4)
netdict = Dict();
for (np = 15*3:1:15*3+1)                      # scan np = 0 .. ne
    listc=Float64[]                    # list of cooperativity values
    lista=Float64[]                    # list of active node fractions
    notconn=0;
    nofp=0;
    samenw=0;
    interc = randomnkp(selfact,ne,np) # generate one sample with given (nkp)
    netdict = Dict();
    i=0
    while (i<ntotal)
#        (intersh,selfsh) = (switchnet(interc),selfact) # preserve in/out degrees
        (intersh,selfsh) = shufflenet(interc,selfact) # shuffle
        if (isconnected(symadj(intersh))) # check if connected
#            netkey=join(sort(real(eigvals(intersh+diagm(selfsh)))))
            netkey=join(reshape(intersh+diagm(selfsh),n*n,1))
#            if (haskey(netdict,netkey)) samenw += 1  # generated before
#            else
            netdict[netkey] = 1      # record this network
            (ncop,actratio) = coopattr(intersh,selfsh) # get cooperativity
##            (ncop,actratio) = cooppointattr(intersh,selfsh) # get cooperativity
##                (ncop,actratio) = coopbasinweighted(intersh,selfsh) # get cooperativity
            if (ncop >= 0.99)  # will return -1 if no nontrivial fp
                push!(listc,ncop)
                println(ncop)
                push!(lista,actratio)
            else
                nofp += 1
            end
#                        end
        else
            notconn += 1
        end
    end
    println("\n", np/ne,"\t",mean(listc),"\t",std(listc),"\t",samenw,"\t",mean(lista))
    output = [output;[np/ne,mean(listc),std(listc),mean(lista)]']
end
tend = time()
println("time:",tend-tstart)
plot(output[:,1],output[:,2])

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
export Init
"""
    Init(data ; fast=false)

Generate an initial data to be used in the function `Problem`.

`data` should contain either

- a function `η` and a function `v` (in this order)
- a Namedtuple with a function `η` and a function `v`
- a mesh and two vectors representing η(mesh.x) and v(mesh.x) (in this order)
- a mesh and a Namedtuple with a vector `η` and a vector `v` as above
- an array of collocation points and two vectors representing `η(x)` and `v(x)` (in this order)
- a mesh and a Namedtuple with a vector `η` and a vector `v` as above

In the last four cases, an optional keyword argument `fast` can be set to `true`,
(default is `false`), in which case the algorithm is faster and uses less allocations, but is less precise.

In the last two cases, the collocation points must be regularly spaced, otherwise an `ErrorException` is raised.

"""
struct Init <: InitialData

    η
    v

    function Init(η , v)
        new( x->η(x) , x->v(x) )
    end

    function Init(p :: NamedTuple)
        new( x->p.η(x) , x->p.v(x) )
    end

    function Init(mesh :: Mesh, η0 , v0 ; fast = false)
        new( x->interpolate(mesh,η0,x;fast=fast) , x->interpolate(mesh,v0,x;fast=fast) )
    end

    function Init(mesh :: Mesh, p :: NamedTuple ; fast = false)
        new( x->interpolate(mesh,p.η,x;fast=fast) , x->interpolate(mesh,p.v,x;fast=fast) )
    end

    function Init(x :: Array, η0 , v0 ; fast = false)
        y=x[2:end]-x[1:end-1];y.-=y[1];
        if maximum(abs.(y))>8*eps(maximum(x))
            @error("Collocation points must be equally spaced.")
        end
        mesh=Mesh(x)
        Init(mesh, η0 , v0 ; fast = fast)
    end

    function Init(x :: Array, p :: NamedTuple ; fast = false)
        y=x[2:end]-x[1:end-1];y.-=y[1];
        if maximum(abs.(y))>8*eps(maximum(x))
            @error("Collocation points must be equally spaced.")
        end
        mesh=Mesh(x)
        Init(mesh, p ; fast = fast)
    end

end


text:
# ---
# title: 814. Binary Tree Pruning
# id: problem814
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/binary-tree-pruning/description/>
# hidden: true
# ---
#
# We are given the head node `root` of a binary tree, where additionally every
# node's value is either a 0 or a 1.
#
# Return the same tree where every subtree (of the given tree) not containing a
# 1 has been removed.
#
# (Recall that the subtree of a node X is X, plus every node that is a
# descendant of X.)
#
#
#
#     **Example 1:**
#     Input: [1,null,0,0,1]
#     Output: [1,null,0,null,1]
#
#     Explanation:
#     Only the red nodes satisfy the property "every subtree not containing a 1".
#     The diagram on the right represents the answer.
#
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/04/06/1028_2.png)
#
#
#
#     **Example 2:**
#     Input: [1,0,1,0,0,0,1]
#     Output: [1,null,1,null,1]
#
#
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/04/06/1028_1.png)
#
#
#
#     **Example 3:**
#     Input: [1,1,0,1,1,0,1,0]
#     Output: [1,1,0,1,1,null,1]
#
#
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/04/05/1028.png)
#
#
# **Note:**
#
#   * The binary tree will have at most `200 nodes`.
#   * The value of each node will only be `0` or `1`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Documenter, ModularForms

makedocs(sitename="ModularForms.jl: computation with modular forms in Julia")

deploydocs(
    repo = "github.com/aghitza/ModularForms.jl.git",
)


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

module BitArrayTests

using Base.Test
using Base: findprevnot, findnextnot

tc{N}(r1::NTuple{N,Any}, r2::NTuple{N,Any}) = all(x->tc(x...), [zip(r1,r2)...])
tc{N}(r1::BitArray{N}, r2::Union{BitArray{N},Array{Bool,N}}) = true
tc(r1::RowVector{Bool,BitVector}, r2::Union{RowVector{Bool,BitVector},RowVector{Bool,Vector{Bool}}}) = true
tc{T}(r1::T, r2::T) = true
tc(r1,r2) = false

bitcheck(b::BitArray) = Base._check_bitarray_consistency(b)
bitcheck(x) = true

function check_bitop_call(ret_type, func, args...)
    r1 = func(args...)
    r2 = func(map(x->(isa(x, BitArray) ? Array(x) : x), args)...)
    ret_type ≢ nothing && !isa(r1, ret_type) && @show ret_type, r1
    ret_type ≢ nothing && @test isa(r1, ret_type)
    @test tc(r1, r2)
    @test isequal(r1, ret_type ≡ nothing ? r2 : convert(ret_type, r2))
    @test bitcheck(r1)
end
macro check_bit_operation(ex, ret_type)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop_call, esc(ret_type), map(esc, ex.args)...)
end
macro check_bit_operation(ex)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop_call, nothing, map(esc, ex.args)...)
end

let t0 = time()
    global timesofar
    function timesofar(str)
        return # no-op, comment to see timings
        t1 = time()
        println(str, ": ", t1-t0, " seconds")
        t0 = t1
    end
end

@testset "empty bitvector" begin
    @test BitVector() == BitVector(0)
end

# vectors size
v1 = 260
# matrices size
n1, n2 = 17, 20
# arrays size
s1, s2, s3, s4 = 5, 8, 3, 7

allsizes = [((), BitArray{0}), ((v1,), BitVector),
            ((n1,n2), BitMatrix), ((s1,s2,s3,s4), BitArray{4})]

@testset "trues and falses for size $sz" for (sz,T) in allsizes
    a = falses(sz...)
    @test a == falses(sz)
    @test !any(a)
    @test sz == size(a)
    b = trues(sz...)
    @test b == trues(sz)
    @test all(b)
    @test sz == size(b)
    c = trues(a)
    @test all(c)
    @test !any(a)
    @test sz == size(c)
    d = falses(b)
    @test !any(d)
    @test all(b)
    @test sz == size(d)
end

@testset "Conversions for size $sz" for (sz, T) in allsizes
    b1 = rand!(falses(sz...))
    @test isequal(BitArray(Array(b1)), b1)
    @test isequal(convert(Array{Float64,ndims(b1)}, b1),
                  convert(Array{Float64,ndims(b1)}, Array(b1)))
    @test isequal(convert(AbstractArray{Float64,ndims(b1)}, b1),
                  convert(AbstractArray{Float64,ndims(b1)}, Array(b1)))

    i1 = rand!(zeros(Bool, sz...), false:true)
    @test isequal(Array(BitArray(i1)), i1)
end

timesofar("conversions")

@testset "utility functions" begin
    b1 = bitrand(v1)
    @test isequal(fill!(b1, true), trues(size(b1)))
    @test isequal(fill!(b1, false), falses(size(b1)))

    for (sz,T) in allsizes
        @test isequal(Array(trues(sz...)), ones(Bool, sz...))
        @test isequal(Array(falses(sz...)), zeros(Bool, sz...))

        b1 = rand!(falses(sz...))
        @test isa(b1, T)

        @check_bit_operation length(b1) Int
        @check_bit_operation ndims(b1)  Int
        @check_bit_operation size(b1)   Tuple{Vararg{Int}}

        b2 = similar(b1)
        u1 = Array(b1)
        @check_bit_operation copy!(b2, b1) T
        @check_bit_operation copy!(b2, u1) T
    end

    @testset "copy!" begin
        for n in [1; 1023:1025]
            b1 = falses(n)
            for m in [1; 10; 1023:1025]
                u1 = ones(Bool, m)
                for fu! in [u->fill!(u, true), u->rand!(u)]
                    fu!(u1)
                    c1 = convert(Vector{Int}, u1)
                    for i1 in [1; 10; 53:65; 1013:1015; 1020:1025], i2 in [1; 3; 10; 511:513], l in [1; 5; 10; 511:513; 1023:1025]
                        for fb! in [b->fill!(b, false), b->rand!(b)]
                            fb!(b1)
                            if i1 < 1 || i1 > n || (i2 + l - 1 > m) || (i1 + l - 1 > n)
                                @test_throws BoundsError copy!(b1, i1, u1, i2, l)
                            else
                                @check_bit_operation copy!(b1, i1, u1, i2, l) BitArray
                                @check_bit_operation copy!(b1, i1, c1, i2, l) BitArray
                            end
                        end
                    end
                end
            end
        end
    end

    @test_throws BoundsError size(trues(5), 0)

    @testset "reshape and resize!" begin
        b1 = bitrand(n1, n2)
        @check_bit_operation reshape(b1, (n2,n1)) BitMatrix
        @test_throws DimensionMismatch reshape(b1, (1,n1))

        b1 = bitrand(s1, s2, s3, s4)
        @check_bit_operation reshape(b1, (s3,s1,s2,s4)) BitArray{4}
        @test_throws DimensionMismatch reshape(b1, (1,n1))

        b1 = bitrand(v1)
        @test_throws BoundsError resize!(b1, -1)
        @check_bit_operation resize!(b1, v1 ÷ 2) BitVector
        gr(b) = (resize!(b, v1)[(v1÷2):end] = 1; b)
        @check_bit_operation gr(b1) BitVector
    end

    @testset "sizeof (issue #7515)" begin
        @test sizeof(BitArray(64)) == 8
        @test sizeof(BitArray(65)) == 16
    end
end

timesofar("utils")

@testset "Constructors" begin
    @testset "non-Int dims constructors" begin
        b1 = BitArray(Int32(v1))
        b2 = BitArray(Int64(v1))
        @test size(b1) == size(b2)

        for c in [trues, falses]
            b1 = c(Int32(v1))
            b2 = c(Int64(v1))
            @test b1 == b2
        end
    end

    @testset "constructors from iterables" begin
        for g in ((x%7==3 for x = 1:v1),
                  (x%7==3 for x = 1:v1 if x>5),
                  ((x+y)%5==2 for x = 1:n1, y = 1:n2),
                  ((x+y+z+t)%5==2 for x = 1:s2, y = 1:s2, z = 1:s3, t = 1:s4),
                  ((x+y)%5==2 for x = 1:n1 for y = 1:n2))
            @test BitArray(g) == BitArray(collect(g))
        end
    end

    @testset "one" begin
        @test Array(one(BitMatrix(2,2))) == eye(2,2)
        @test_throws DimensionMismatch one(BitMatrix(2,3))
    end
end

timesofar("constructors")

@testset "Indexing" begin
    @testset "0d for size $sz" for (sz,T) in allsizes
        b1 = rand!(falses(sz...))
        @check_bit_operation getindex(b1)         Bool
        @check_bit_operation setindex!(b1, true)  T
        @check_bit_operation setindex!(b1, false) T
    end

    @testset "linear for size $sz" for (sz,T) in allsizes[2:end]
        l = *(sz...)
        b1 = rand!(falses(sz...))
        for j = 1:l
            @check_bit_operation getindex(b1, j) Bool
        end
        for j in [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, l-1, l]
            @check_bit_operation getindex(b1, 1:j)   BitVector
            @check_bit_operation getindex(b1, j+1:l) BitVector
        end
        for j in [1, 63, 64, 65, 127, 128, 129, div(l,2)]
            m1 = j:(l-j)
            @check_bit_operation getindex(b1, m1) BitVector
        end

        t1 = find(bitrand(l))
        @check_bit_operation getindex(b1, t1)        BitVector

        for j = 1:l
            x = rand(Bool)
            @check_bit_operation setindex!(b1, x, j) T
        end

        y = rand(0.0:1.0)
        @check_bit_operation setindex!(b1, y, 100) T

        for j in [1, 63, 64, 65, 127, 128, 129, 191, 192, 193, l-1]
            x = rand(Bool)
            @check_bit_operation setindex!(b1, x, 1:j) T
            b2 = bitrand(j)
            for bb in (b2, view(b2, 1:j), view(Array{Any}(b2), :))
                @check_bit_operation setindex!(b1, bb, 1:j) T
            end
            x = rand(Bool)
            @check_bit_operation setindex!(b1, x, j+1:l) T
            b2 = bitrand(l-j)
            @check_bit_operation setindex!(b1, b2, j+1:l) T
        end
        for j in [1, 63, 64, 65, 127, 128, 129, div(l,2)]
            m1 = j:(l-j)
            x = rand(Bool)
            @check_bit_operation setindex!(b1, x, m1) T
            b2 = bitrand(length(m1))
            @check_bit_operation setindex!(b1, b2, m1) T
        end
        x = rand(Bool)
        @check_bit_operation setindex!(b1, x, 1:100) T
        b2 = bitrand(100)
        @check_bit_operation setindex!(b1, b2, 1:100) T

        y = rand(0.0:1.0)
        @check_bit_operation setindex!(b1, y, 1:100) T

        t1 = find(bitrand(l))
        x = rand(Bool)
        @check_bit_operation setindex!(b1, x, t1) T
        b2 = bitrand(length(t1))
        @check_bit_operation setindex!(b1, b2, t1) T

        y = rand(0.0:1.0)
        @check_bit_operation setindex!(b1, y, t1) T
    end

    @testset "multidimensional" begin
        rand_m1m2() = rand(1:n1), rand(1:n2)
        b1 = bitrand(n1, n2)
        m1, m2 = rand_m1m2()
        b2 = bitrand(m1, m2)
        @check_bit_operation copy!(b1, b2) BitMatrix

        function gen_getindex_data(c)
            m1, m2 = rand_m1m2()
            put!(c, (m1, m2, Bool))
            m1, m2 = rand_m1m2()
            put!(c, (m1, 1:m2, BitVector))
            put!(c, (m1, :, BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (m1, randperm(m2), BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (1:m1, m2, BitVector))
            put!(c, (:, m2, BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (1:m1, 1:m2, BitMatrix))
            put!(c, (:, :, BitMatrix))
            m1, m2 = rand_m1m2()
            put!(c, (1:m1, randperm(m2), BitMatrix))
            put!(c, (:, randperm(m2), BitMatrix))
            m1, m2 = rand_m1m2()
            put!(c, (randperm(m1), m2, BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (randperm(m1), 1:m2, BitMatrix))
            put!(c, (randperm(m1), :, BitMatrix))
            m1, m2 = rand_m1m2()
            put!(c, (randperm(m1), randperm(m2), BitMatrix))
        end

        for (k1, k2, T) in Channel(gen_getindex_data)
            # println(typeof(k1), " ", typeof(k2), " ", T) # uncomment to debug
            @check_bit_operation getindex(b1, k1, k2) T
        end

        m1, m2 = rand_m1m2()
        @check_bit_operation getindex(b1, 1:m1, m2, 1) BitVector
        @check_bit_operation getindex(b1, :, randperm(m2), 1) BitMatrix

        b1 = bitrand(s1, s2, s3, s4)
        function gen_getindex_data4(c)
            m1, m2, m3, m4 = (:, :, :, :)
            put!(c, (m1, m2, m3, m4, BitArray{4}))

            m1, m2, m3, m4 = (2, :, :, :)
            put!(c, (m1, m2, m3, m4, BitArray{3}))
            m1, m2, m3, m4 = (:, :, 2, :)
            put!(c, (m1, m2, m3, m4, BitArray{3}))
            m1, m2, m3, m4 = (:, :, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{3}))

            m1, m2, m3, m4 = (2, :, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (:, 2, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (2, :, 2, :)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (2, 2, :, :)
            put!(c, (m1, m2, m3, m4, BitArray{2}))

            m1, m2, m3, m4 = (:, 2, 2, 2)
            put!(c, (m1, m2, m3, m4, BitArray{1}))
            m1, m2, m3, m4 = (2, 2, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{1}))

            m1, m2, m3, m4 = (4, 3:7, 2:2, 2)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (1:2, 5, 1, 2:7)
            put!(c, (m1, m2, m3, m4, BitArray{2}))

            m1, m2, m3, m4 = (2:3, 2:7, 1:2, 4:6)
            put!(c, (m1, m2, m3, m4, BitArray{4}))
        end

        for (k1, k2, k3, k4, T) in Channel(gen_getindex_data4)
            #println(typeof(k1), " ", typeof(k2), " ", typeof(k3), " ", typeof(k4), " ", T) # uncomment to debug
            @check_bit_operation getindex(b1, k1, k2, k3, k4) T
        end

        b1 = bitrand(n1, n2)
        function gen_setindex_data(c)
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), m1, m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), m1, 1:m2))
            put!(c, (rand(Bool), m1, :))
            put!(c, (bitrand(m2), m1, 1:m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), m1, randperm(m2)))
            put!(c, (bitrand(m2), m1, randperm(m2)))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), 1:m1, m2))
            put!(c, (rand(Bool), :, m2))
            put!(c, (bitrand(m1), 1:m1, m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), 1:m1, 1:m2))
            put!(c, (rand(Bool), :, :))
            put!(c, (bitrand(m1, m2), 1:m1, 1:m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), 1:m1, randperm(m2)))
            put!(c, (rand(Bool), :, randperm(m2)))
            put!(c, (bitrand(m1, m2), 1:m1, randperm(m2)))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), randperm(m1), m2))
            put!(c, (bitrand(m1), randperm(m1), m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), randperm(m1), 1:m2))
            put!(c, (rand(Bool), randperm(m1), :))
            put!(c, (bitrand(m1,m2), randperm(m1), 1:m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), randperm(m1), randperm(m2)))
            put!(c, (bitrand(m1,m2), randperm(m1), randperm(m2)))
        end

        for (b2, k1, k2) in Channel(gen_setindex_data)
            # println(typeof(b2), " ", typeof(k1), " ", typeof(k2)) # uncomment to debug
            for bb in ((b2 isa AbstractArray) ? (b2, view(b2, :), view(Array{Any}(b2), :)) : (b2,))
                @check_bit_operation setindex!(b1, bb, k1, k2) BitMatrix
            end
        end

        m1, m2 = rand_m1m2()
        b2 = bitrand(1, 1, m2)
        @check_bit_operation setindex!(b1, b2, m1, 1:m2) BitMatrix
        x = rand(Bool)
        b2 = bitrand(1, m2, 1)
        @check_bit_operation setindex!(b1, x, m1, 1:m2, 1)  BitMatrix
        @check_bit_operation setindex!(b1, b2, m1, 1:m2, 1) BitMatrix

        b1 = bitrand(s1, s2, s3, s4)
        function gen_setindex_data4(c)
            m1, m2, m3, m4 = (:, :, :, :)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s1, s2, s3, s4), m1, m2, m3, m4))

            m1, m2, m3, m4 = (2, :, :, :)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s2, s3, s4), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, :, 2, :)
            put!(c, (bitrand(s1, s2, s4), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, :, :, 2)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s1, s2, s3), m1, m2, m3, m4))

            m1, m2, m3, m4 = (2, :, :, 2)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s2, s3), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, 2, :, 2)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s1, s3), m1, m2, m3, m4))
            m1, m2, m3, m4 = (2, :, 2, :)
            put!(c, (bitrand(s2, s4), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, 2, 2, :)
            put!(c, (bitrand(s1, s4), m1, m2, m3, m4))

            m1, m2, m3, m4 = (:, 2, 2, 2)
            put!(c, (bitrand(s1), m1, m2, m3, m4))
            m1, m2, m3, m4 = (2, 2, :, 2)
            put!(c, (bitrand(s3), m1, m2, m3, m4))

            m1, m2, m3, m4 = (4, 3:7, 2:2, 2)
            put!(c, (bitrand(5, 1), m1, m2, m3, m4))
            m1, m2, m3, m4 = (1:2, 5, 1, 2:7)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(2, 6), m1, m2, m3, m4))

            m1, m2, m3, m4 = (2:3, 2:7, 1:2, 4:6)
            put!(c, (bitrand(2, 6, 2, 3), m1, m2, m3, m4))
        end

        for (b2, k1, k2, k3, k4) in Channel(gen_setindex_data4)
            # println(typeof(b2), " ", typeof(k1), " ", typeof(k2), " ", typeof(k3), " ", typeof(k4)) # uncomment to debug
            @check_bit_operation setindex!(b1, b2, k1, k2, k3, k4) BitArray{4}
        end

        for p1 = [rand(1:v1) 1 63 64 65 191 192 193]
            for p2 = [rand(1:v1) 1 63 64 65 191 192 193]
                for n = 0 : min(v1 - p1 + 1, v1 - p2 + 1)
                    b1 = bitrand(v1)
                    b2 = bitrand(v1)
                    @check_bit_operation copy!(b1, p1, b2, p2, n) BitVector
                end
            end
        end
    end

    @testset "logical indexing" begin
        b1 = bitrand(n1, n2)

        t1 = bitrand(n1, n2)
        @test isequal(Array(b1[t1]), Array(b1)[t1])
        @test isequal(Array(b1[t1]), Array(b1)[Array(t1)])

        t1 = bitrand(n1)
        t2 = bitrand(n2)
        @test isequal(Array(b1[t1, t2]), Array(b1)[t1, t2])
        @test isequal(Array(b1[t1, t2]), Array(b1)[Array(t1), Array(t2)])

        b1 = bitrand(n1, n2)
        t1 = bitrand(n1, n2)
        @check_bit_operation setindex!(b1, true, t1) BitMatrix

        t1 = bitrand(n1, n2)
        b2 = bitrand(countnz(t1))
        @check_bit_operation setindex!(b1, b2, t1) BitMatrix

        m1 = rand(1:n1)
        m2 = rand(1:n2)
        t1 = bitrand(n1)
        b2 = bitrand(countnz(t1), m2)
        k2 = randperm(m2)
        @check_bit_operation setindex!(b1, b2, t1, 1:m2)       BitMatrix
        @check_bit_operation setindex!(b1, b2, t1, n2-m2+1:n2) BitMatrix
        @check_bit_operation setindex!(b1, b2, t1, k2)         BitMatrix

        t2 = bitrand(n2)
        b2 = bitrand(m1, countnz(t2))
        k1 = randperm(m1)
        @check_bit_operation setindex!(b1, b2, 1:m1, t2)       BitMatrix
        @check_bit_operation setindex!(b1, b2, n1-m1+1:n1, t2) BitMatrix
        @check_bit_operation setindex!(b1, b2, k1, t2)         BitMatrix
    end
end

timesofar("indexing")

@testset "Deque Functionality" begin
    b1 = BitArray(0)
    i1 = Bool[]
    for m = 1:v1
        x = rand(Bool)
        push!(b1, x)
        push!(i1, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    for m1 = 0:v1, m2 = [0, 1, 63, 64, 65, 127, 128, 129]
        b1 = bitrand(m1)
        b2 = bitrand(m2)
        i1 = Array(b1)
        i2 = Array(b2)
        @test isequal(Array(append!(b1, b2)), append!(i1, i2))
        @test isequal(Array(append!(b1, i2)), append!(i1, b2))
        @test bitcheck(b1)
    end

    for m1 = 0:v1, m2 = [0, 1, 63, 64, 65, 127, 128, 129]
        b1 = bitrand(m1)
        b2 = bitrand(m2)
        i1 = Array(b1)
        i2 = Array(b2)
        @test isequal(Array(prepend!(b1, b2)), prepend!(i1, i2))
        @test isequal(Array(prepend!(b1, i2)), prepend!(i1, b2))
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = 1:v1
        jb = pop!(b1)
        ji = pop!(i1)
        @test jb == ji
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    @test length(b1) == 0

    b1 = BitArray(0)
    i1 = Bool[]
    for m = 1:v1
        x = rand(Bool)
        unshift!(b1, x)
        unshift!(i1, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = 1:v1
        jb = shift!(b1)
        ji = shift!(i1)
        @test jb == ji
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end
    @test length(b1) == 0

    b1 = BitArray(0)
    @test_throws BoundsError insert!(b1, 2, false)
    @test_throws BoundsError insert!(b1, 0, false)
    i1 = Array(b1)
    for m = 1:v1
        j = rand(1:m)
        x = rand(Bool)
        @test insert!(b1, j, x) === b1
        insert!(i1, j, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
        x = rand(0:1)
        @test insert!(b1, j, x) === b1
        insert!(i1, j, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = v1:-1:1
        j = rand(1:m)
        b = splice!(b1, j)
        i = splice!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
        @test b == i
        @test bitcheck(b)
    end
    @test length(b1) == 0

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = v1:-1:1
        j = rand(1:m)
        deleteat!(b1, j)
        deleteat!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end
    @test length(b1) == 0
    b1 = bitrand(v1)
    @test_throws ArgumentError deleteat!(b1, [1, 1, 2])
    @test_throws BoundsError deleteat!(b1, [1, length(b1)+1])

    b1 = bitrand(v1)
    i1 = Array(b1)
    for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
        b = splice!(b1, j)
        i = splice!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
        deleteat!(b1, j)
        deleteat!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:v1, m2 = m1:v1
        b2 = copy(b1)
        i2 = copy(i1)
        b = splice!(b2, m1:m2)
        i = splice!(i2, m1:m2)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:v1, m2 = m1:v1
        b2 = copy(b1)
        i2 = copy(i1)
        deleteat!(b2, m1:m2)
        deleteat!(i2, m1:m2)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = [1:7:v1; v1+1], m2 = [(m1-1):5:(v1-1); v1], v2 = [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, rand(1:v1)]
        b2 = copy(b1)
        i2 = copy(i1)
        b3 = bitrand(v2)
        i3 = Array(b3)
        b = splice!(b2, m1:m2, b3)
        i = splice!(i2, m1:m2, i3)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)

        b2 = copy(b1)
        i2 = copy(i1)
        i3 = map(Int, bitrand(v2))
        b = splice!(b2, m1:m2, i3)
        i = splice!(i2, m1:m2, i3)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
        b2 = copy(b1)
        i2 = copy(i1)
        i3 = Dict(j => rand(0:1) for j = 1:v2)
        b = splice!(b2, m1:m2, values(i3))
        i = splice!(i2, m1:m2, values(i3))
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:v1, v2 = [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, rand(1:v1)]
        b2 = copy(b1)
        i2 = copy(i1)
        b3 = bitrand(v2)
        i3 = Array(b3)
        b = splice!(b2, m1, b3)
        i = splice!(i2, m1, i3)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:(v1-1), m2 = (m1+1):v1
        locs = bitrand(m2-m1+1)
        m = [m1:m2...][locs]
        b2 = copy(b1)
        i2 = copy(i1)
        deleteat!(b2, m)
        deleteat!(i2, m)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    empty!(b1)
    empty!(i1)
    @test isequal(Array(b1), i1)
    @test bitcheck(b1)
end

timesofar("dequeue")

@testset "Unary operators" begin
    b1 = bitrand(n1, n2)
    @check_bit_operation broadcast(~, b1)  BitMatrix
    @check_bit_operation broadcast(!, b1)  BitMatrix
    @check_bit_operation (-)(b1)  Matrix{Int}
    @check_bit_operation broadcast(sign, b1) BitMatrix
    @check_bit_operation real(b1) BitMatrix
    @check_bit_operation imag(b1) BitMatrix
    @check_bit_operation conj(b1) BitMatrix

    b0 = falses(0)
    @check_bit_operation broadcast(~, b0)  BitVector
    @check_bit_operation broadcast(!, b0)  BitVector
    @check_bit_operation (-)(b0)  Vector{Int}
    @check_bit_operation broadcast(sign, b0) BitVector

    @testset "flipbits!" begin
        b1 = bitrand(n1, n2)
        i1 = Array(b1)
        @test flipbits!(b1) == .~i1
        @test bitcheck(b1)
    end
end

timesofar("unary arithmetic")

end # module

max_stars_repo_name: jamesonquinn/julia

max_issues_repo_name: jamesonquinn/julia

max_forks_repo_name: jamesonquinn/julia


text:
using DataStructures: OrderedDict
using DataFrames: DataFrame
using Dates: AbstractTime

struct Simulation
    base::Union{String,Symbol,Nothing}
    index::OrderedDict{Symbol,Any}
    target::OrderedDict{Symbol,Any}
    mapping::OrderedDict{Symbol,Any}
    meta::OrderedDict{Symbol,Any}
    result::DataFrame
end

simulation(s::System; config=(), base=nothing, index=nothing, target=nothing, meta=nothing) = begin
    sb = s[base]
    I = parseindex(index, sb)
    T = parsetarget(target, sb)
    #HACK: ignore unavailable properties (i.e. handle default :time in target)
    I = filtersimulationdict(I, sb)
    T = filtersimulationdict(T, sb)
    IT = merge(I, T)
    M = parsemeta(meta, s.context.config)
    Simulation(base, I, T, IT, M, DataFrame())
end

parsesimulationkey(p::Pair, s) = [p]
parsesimulationkey(a::Symbol, s) = [a => a]
parsesimulationkey(a::String, s) = [Symbol(a) => a]
parsesimulationkey(a::String, s::System) = begin
    A = split(a, '.')
    # support wildcard names (i.e. "s.*" expands to ["s.a", "s.b", ...])
    if A[end] == "*"
        a0 = join(A[1:end-1], '.')
        ss = s[a0]
        p(n) = let k = join(filter(!isempty, [a0, string(n)]), '.')
            Symbol(k) => k
        end
        [p(n) for n in fieldnamesunique(ss)]
    else
        [Symbol(a) => a]
    end
end
parsesimulationkey(a::Vector, s) = parsesimulationkey.(a, Ref(s)) |> Iterators.flatten |> collect

parsesimulation(a::Vector, s) = OrderedDict(parsesimulationkey.(a, Ref(s)) |> Iterators.flatten)
parsesimulation(a::Tuple, s) = parsesimulation(collect(a), s)
parsesimulation(::Tuple{}, s) = parsesimulation([], s)
parsesimulation(a, s) = parsesimulation([a], s)

parseindex(::Nothing, s) = parsesimulation(:time => "context.clock.time", s)
parseindex(I, s) = parsesimulation(I, s)

parsetarget(::Nothing, s) = parsesimulation("*", s)
parsetarget(T, s) = parsesimulation(T, s)

filtersimulationdict(m::OrderedDict, s::System) = filter(m) do (k, v); hasproperty(s, v) end

extract(s::System, m::Simulation) = extract(s[m.base], m.mapping)
extract(s::System, m::OrderedDict{Symbol,Any}) = begin
    K = keys(m)
    V = (value(s[k]) for k in values(m))
    #HACK: Any -- prevent type promotion with NoUnits
    d = OrderedDict{Symbol,Any}(zip(K, V))
    filter!(p -> extractable(s, p), d)
    [d]
end
extract(b::Bundle{S}, m::OrderedDict{Symbol,Any}) where {S<:System} = begin
    [extract(s, m) for s in collect(b)]
end
extractable(s::System, p) = begin
    # only pick up variables of simple types by default
    p[2] isa Union{Number,Symbol,AbstractString,AbstractTime}
end

parsemetadata(p::Pair, c) = p
parsemetadata(a::Symbol, c) = c[a]
parsemeta(a::Vector, c) = merge(parsemeta(nothing, c), OrderedDict.(parsemetadata.(a, Ref(c)))...)
parsemeta(a::Tuple, c) = parsemeta(collect(a), c)
parsemeta(a, c) = parsemeta([a], c)
parsemeta(::Nothing, c) = OrderedDict()

update!(m::Simulation, s::System, snatch!) = begin
    D = extract(s, m)
    snatch!(D, s)
    !isempty(D) && append!(m.result, D)
end

format!(m::Simulation; nounit=false, long=false) = begin
    r = m.result
    for (k, v) in m.meta
        r[!, k] .= v
    end
    if nounit
        r = deunitfy.(r)
    end
    if long
        i = collect(keys(m.index))
        t = setdiff(propertynames(r), i)
        r = DataFrames.stack(r, t, i)
    end
    r
end

using ProgressMeter: Progress, ProgressUnknown, ProgressMeter
const barglyphs = ProgressMeter.BarGlyphs("[=> ]")
progress!(s::System, M::Vector{Simulation}; stop=nothing, snap=nothing, snatch=nothing, callback=nothing, verbose=true, kwargs...) = begin
    probe(a::Union{Symbol,String}) = s -> s[a]'
    probe(a::Function) = s -> a(s)
    probe(a) = s -> a

    stopprobe(::Nothing) = probe(0)
    stopprobe(a) = probe(a)

    snapprobe(::Nothing) = probe(true)
    snapprobe(a::Quantity) = s -> let c = s.context.clock; (c.time' - c.init') % a |> iszero end
    snapprobe(a) = probe(a)

    stop = stopprobe(stop)
    snap = snapprobe(snap)
    snatch = isnothing(snatch) ? (D, s) -> nothing : snatch
    callback = isnothing(callback) ? (s, m) -> nothing : callback

    count(v::Number) = v
    count(v::Quantity) = ceil(Int, v / s.context.clock.step')
    count(v::Bool) = nothing
    count(v) = error("unrecognized stop condition: $v")
    n = count(stop(s))

    dt = verbose ? 1 : Inf
    showspeed = true
    if n isa Number
        p = Progress(n; dt, barglyphs, showspeed)
        check = s -> p.counter < p.n
    else
        p = ProgressUnknown(; dt, desc="Iterations:", showspeed)
        check = s -> !stop(s)
    end

    snap(s) && update!.(M, s, snatch)
    while check(s)
        update!(s)
        snap(s) && for m in M
            update!(m, s, snatch)
            callback(s, m)
        end
        ProgressMeter.next!(p)
    end
    ProgressMeter.finish!(p)
    format!.(M; kwargs...)
end

"""
    simulate!([f,] s[, layout]; <keyword arguments>) -> DataFrame

Run simulations with an existing instance of system `s`. The instance is altered by internal updates for running simulations.

See also: [`simulate`](@ref)
"""
simulate!(s::System; base=nothing, index=nothing, target=nothing, meta=nothing, kwargs...) = begin
    simulate!(s, [(; base, index, target, meta)]; kwargs...) |> only
end
simulate!(s::System, layout::Vector; kwargs...) = begin
    M = [simulation(s; l...) for l in layout]
    progress!(s, M; kwargs...)
end
simulate!(f::Function, s::System, args...; kwargs...) = simulate!(s, args...; snatch=f, kwargs...)

"""
    simulate([f,] S[, layout, [configs]]; <keyword arguments>) -> DataFrame

Run simulations by making instance of system `S` with given configuration to generate an output in the form of DataFrame. `layout` contains a list of variables to be saved in the output. A layout of single simulation can be specified in the layout arguments placed as keyword arguments. `configs` contains a list of configurations for each run of simulation. Total number of simulation runs equals to the size of `configs`. For a single configuration, `config` keyword argument may be preferred. Optional callback function `f` allows do-block syntax to specify `snatch` argument for finer control of output format.

See also: [`instance`](@ref), [`@config`](@ref)

# Arguments
- `S::Type{<:System}`: type of system to be simulated.
- `layout::Vector`: list of output layout definition in a named tuple `(; base, index, target, meta)`.
- `configs::Vector`: list of configurations for defining multiple runs of simluations.

# Keyword Arguments
## Layout
- `base=nothing`: base system where `index` and `target` are populated; default falls back to the instance of `S`.
- `index=nothing`: variables to construct index columns of the output; default falls back to `context.clock.time`.
- `target=nothing`: variables to construct non-index columns of the output; default includes most variables in the root instance.
- `meta=nothing`: name of systems in the configuration to be included in the output as metadata.

## Configuration
- `config=()`: a single configuration for the system (can't be used with `configs`).
- `configs=[]`: multiple configurations for the system (can't be used with `config` nor `configs` argument).
- `seed=nothing`: random seed for resetting each simulation run.

## Progress
- `stop=nothing`: condition checked before calling updates for the instance; default stops with no update.
- `snap=nothing`: condition checked to decide if a snapshot of current update is saved in the output; default snaps all updates.
- `snatch=nothing`: callback for modifying intermediate output; list of DataFrame `D` collected from current update and the instance of system `s` are provided.
- `verbose=true`: shows a progress bar.

## Format
- `nounit=false`: remove units from the output.
- `long=false`: convert output table from wide to long format.

# Examples
```julia-repl
julia> @system S(Controller) begin
           a => 1 ~ preserve(parameter)
           b(a) ~ accumulate
       end;

julia> simulate(S; stop=1)
2×3 DataFrame
 Row │ time       a        b
     │ Quantity…  Float64  Float64
─────┼─────────────────────────────
   1 │    0.0 hr      1.0      0.0
   2 │    1.0 hr      1.0      1.0
```
"""
simulate(; system, kw...) = simulate(system; kw...)
simulate(S::Type{<:System}; base=nothing, index=nothing, target=nothing, meta=nothing, kwargs...) = begin
    simulate(S, [(; base, index, target, meta)]; kwargs...) |> only
end
simulate(S::Type{<:System}, layout::Vector; config=(), configs=[], options=(), seed=nothing, kwargs...) = begin
    if isempty(configs)
        s = instance(S; config, options, seed)
        simulate!(s, layout; kwargs...)
    elseif isempty(config)
        simulate(S, layout, configs; options, seed, kwargs...)
    else
        @error "redundant configurations" config configs
    end
end
simulate(S::Type{<:System}, layout::Vector, configs::Vector; verbose=true, kwargs...) = begin
    n = length(configs)
    R = Vector(undef, n)
    dt = verbose ? 1 : Inf
    p = Progress(n; dt, barglyphs)
    Threads.@threads for i in 1:n
        R[i] = simulate(S, layout; config=configs[i], verbose=false, kwargs...)
        ProgressMeter.next!(p)
    end
    ProgressMeter.finish!(p)
    [vcat(r...) for r in eachrow(hcat(R...))]
end
simulate(f::Function, S::Type{<:System}, args...; kwargs...) = simulate(S, args...; snatch=f, kwargs...)

export simulate, simulate!


text: abstract type AbstractArbMatrix{P, T}  <: AbstractMatrix{T} end


text:
# # Polyhedral QP layer

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/polyhedral_project.jl)

# We use DiffOpt to define a custom network layer which, given an input matrix `y`,
# computes its projection onto a polytope defined by a fixed number of inequalities:
# `a_i^T x ≥ b_i`.
# A neural network is created using Flux.jl and trained on the MNIST dataset,
# integrating this quadratic optimization layer.
#
# The QP is solved in the forward pass, and its DiffOpt derivative is used in the backward pass expressed with `ChainRulesCore.rrule`.

# This example is similar to the custom ReLU layer, except that the layer is parameterized
# by the hyperplanes `(w,b)` and not a simple stateless function.
# This also means that `ChainRulesCore.rrule` must return the derivatives of the output with respect to the
# layer parameters to allow for backpropagation.

using JuMP
import DiffOpt
import Ipopt
import ChainRulesCore
import Flux
import MLDatasets
import Statistics
using Base.Iterators: repeated
using LinearAlgebra
using Random

Random.seed!(42)

# ## The Polytope representation and its derivative

struct Polytope{N}
    w::NTuple{N, Vector{Float64}}
    b::Vector{Float64}
end

Polytope(w::NTuple{N}) where {N} = Polytope{N}(w, randn(N))

# We define a "call" operation on the polytope, making it a so-called functor.
# Calling the polytope with a matrix `y` operates an Euclidean projection of this matrix onto the polytope.
function (polytope::Polytope{N})(
    y::AbstractMatrix;
    model = direct_model(DiffOpt.diff_optimizer(Ipopt.Optimizer))
) where {N}
    layer_size, batch_size = size(y)
    empty!(model)
    set_silent(model)
    @variable(model, x[1:layer_size, 1:batch_size])
    @constraint(model,
        greater_than_cons[idx in 1:N, sample in 1:batch_size],
        dot(polytope.w[idx], x[:, sample]) ≥ polytope.b[idx]
    )
    @objective(model, Min, dot(x - y, x - y))
    optimize!(model)
    return JuMP.value.(x)
end

# The `@functor` macro from Flux implements auxiliary functions for collecting the parameters of
# our custom layer and operating backpropagation.
Flux.@functor Polytope

# Define the reverse differentiation rule, for the function we defined above.
# Flux uses ChainRules primitives to implement reverse-mode differentiation of the whole network.
# To learn the current layer (the polytope the layer contains),
# the gradient is computed with respect to the `Polytope` fields in a ChainRulesCore.Tangent type
# which is used to represent derivatives with respect to structs.
# For more details about backpropagation, visit [Introduction, ChainRulesCore.jl](https://juliadiff.org/ChainRulesCore.jl/dev/).

function ChainRulesCore.rrule(polytope::Polytope{N}, y::AbstractMatrix) where {N}
    model = direct_model(DiffOpt.diff_optimizer(Ipopt.Optimizer))
    xv = polytope(y; model = model)
    function pullback_matrix_projection(dl_dx)
        layer_size, batch_size = size(dl_dx)
        dl_dx = ChainRulesCore.unthunk(dl_dx)
        ##  `dl_dy` is the derivative of `l` wrt `y`
        x = model[:x]
        ## grad wrt input parameters
        dl_dy = zeros(size(dl_dx))
        ## grad wrt layer parameters
        dl_dw = zero.(polytope.w)
        dl_db = zero(polytope.b)
        ## set sensitivities
        MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, dl_dx)
        ## compute grad
        DiffOpt.reverse_differentiate!(model)
        ## compute gradient wrt objective function parameter y
        obj_expr = MOI.get(model, DiffOpt.ReverseObjectiveFunction())
        dl_dy .= -2 * JuMP.coefficient.(obj_expr, x)
        greater_than_cons = model[:greater_than_cons]
        for idx in 1:N, sample in 1:batch_size
            cons_expr = MOI.get(model,
                DiffOpt.ReverseConstraintFunction(),
                greater_than_cons[idx, sample])
            dl_db[idx] -= JuMP.constant(cons_expr)/batch_size
            dl_dw[idx] .+= JuMP.coefficient.(cons_expr, x[:,sample])/batch_size
        end
        dself = ChainRulesCore.Tangent{Polytope{N}}(; w = dl_dw, b = dl_db)
        return (dself, dl_dy)
    end
    return xv, pullback_matrix_projection
end

# ## Define the Network

layer_size = 20
m = Flux.Chain(
    Flux.Dense(784, layer_size), # 784 being image linear dimension (28 x 28)
    Polytope((randn(layer_size), randn(layer_size), randn(layer_size))),
    Flux.Dense(layer_size, 10), # 10 being the number of outcomes (0 to 9)
    Flux.softmax,
)

# ## Prepare data

M = 500 # batch size
## Preprocessing train data
imgs = MLDatasets.MNIST.traintensor(1:M)
labels = MLDatasets.MNIST.trainlabels(1:M);
train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), M)) # stack images
train_Y = Flux.onehotbatch(labels, 0:9);
## Preprocessing test data
test_imgs = MLDatasets.MNIST.testtensor(1:M)
test_labels = MLDatasets.MNIST.testlabels(1:M)
test_X = float.(reshape(test_imgs, size(test_imgs, 1) * size(test_imgs, 2), M))
test_Y = Flux.onehotbatch(test_labels, 0:9);

# Define input data
# The original data is repeated `epochs` times because `Flux.train!` only
# loops through the data set once

epochs = 50
dataset = repeated((train_X, train_Y), epochs);

# ## Network training

# training loss function, Flux optimizer
custom_loss(x, y) = Flux.crossentropy(m(x), y)
opt = Flux.ADAM()
evalcb = () -> @show(custom_loss(train_X, train_Y))

# Train to optimize network parameters

@time Flux.train!(custom_loss, Flux.params(m), dataset, opt, cb = Flux.throttle(evalcb, 5));

# Although our custom implementation takes time, it is able to reach similar
# accuracy as the usual ReLU function implementation.

# ## Accuracy results

# Average of correct guesses
accuracy(x, y) = Statistics.mean(Flux.onecold(m(x)) .== Flux.onecold(y));

# Training accuracy

accuracy(train_X, train_Y)

# Test accuracy

accuracy(test_X, test_Y)

# Note that the accuracy is low due to simplified training.
# It is possible to increase the number of samples `N`,
# the number of epochs `epoch` and the connectivity `inner`.


text:
if normpath(@__DIR__) ∉ LOAD_PATH
    pushfirst!(LOAD_PATH, normpath(@__DIR__, "../.."))
    pushfirst!(LOAD_PATH, normpath(@__DIR__))
end

using DECAES
using Statistics
using LaTeXStrings
using CairoMakie
set_theme!(theme_minimal(); resolution = (800,600), font = "CMU Serif")

const AbstractTensor{N} = AbstractArray{T,N} where {T}

gridsize(n) = (imax = floor(Int, sqrt(n)); jmax = ceil(Int, n / imax); (imax, jmax))

function global_reduction(reducer, img::AbstractTensor{4}, mask = (img[:,:,:,1] .> 0))
    [reducer(img[mask,echo]) for echo in axes(img, 4)]
end

function plot_echoes(img::AbstractTensor{4}; echoes)
    @assert size(img, 3) == 1
    fig = Figure()
    nplots = length(echoes)
    limits = extrema(vec(img))
    for (i,I) in enumerate(CartesianIndices(gridsize(nplots)))
        i > nplots && break
        ax = Axis(fig[I[1], I[2]]; aspect = AxisAspect(1))
        hidedecorations!(ax)
        heatmap!(img[:,:,1,echoes[i]]; limits = limits)
        text!("$(echoes[i])"; color = :white, position = (1.5, 1.5), align = (:left, :baseline))
    end
    return fig
end

function plot_band!(ax, ydata::AbstractTensor{4}, xdata = 1:size(ydata, 4); meankwargs, bandkwargs)
    μ = global_reduction(mean, ydata)
    σ = global_reduction(std, ydata)
    (meankwargs !== nothing) && scatterlines!(ax, xdata, μ; meankwargs...)
    (bandkwargs !== nothing) && band!(ax, xdata, μ .- σ, μ .+ σ; bandkwargs...)
end

function plot_global_signal(img::AbstractTensor{4};
        meankwargs = (label = L"\mu", color = :red, markersize = 5, linewidth = 1),
        bandkwargs = (label = L"\mu \pm \sigma", color = (:red, 0.5)),
    )
    fig = Figure()
    ax = Axis(fig[1,1])
    plot_band!(ax, img; meankwargs = meankwargs, bandkwargs = bandkwargs)
    axislegend(ax)
    return fig
end

function compare_refcon_correction(
        img::AbstractTensor{4},
        mask::AbstractTensor{3},
        t2maps_norefcon::AbstractDict,
        t2dist_norefcon::AbstractTensor{4},
        t2maps_refcon::AbstractDict,
        t2dist_refcon::AbstractTensor{4},
    )
    DECAES.@unpack t2times, echotimes = t2maps_norefcon
    t2times, echotimes = 1000 .* t2times, 1000 .* echotimes
    globalmean(x) = vec(mean(x[mask, :]; dims = 1))

    fig = Figure()
    ax = Axis(fig[1,1]; ylabel = "mean signal [a.u.]", xlabel = "echo time [ms]")
    lines!(ax, echotimes, globalmean(img); label = "data", markersize = 5, linewidth = 3, color = :red)
    lines!(ax, echotimes, globalmean(t2maps_norefcon["decaycurve"]); label = "no correction", markersize = 5, linewidth = 3, color = :blue)
    lines!(ax, echotimes, globalmean(t2maps_refcon["decaycurve"]); label = "w/ correction", markersize = 5, linewidth = 3, color = :green)
    hideydecorations!(ax, label = false)
    axislegend(ax)

    ax = Axis(fig[1,2]; ylabel = "mean difference [a.u.]", xlabel = "echo time [ms]")
    hlines!(ax, [0.0]; label = "zero", linewidth = 3, linestyle = :dot, color = :red)
    lines!(ax, echotimes, globalmean(img .- t2maps_norefcon["decaycurve"]); label = "no correction", markersize = 5, linewidth = 3, color = :blue)
    lines!(ax, echotimes, globalmean(img .- t2maps_refcon["decaycurve"]); label = "w/ correction", markersize = 5, linewidth = 3, color = :green)
    hideydecorations!(ax, label = false)
    axislegend(ax)

    ax = Axis(fig[2,1]; xlabel = "fit error [a.u.]", ylabel = "relative count")
    hist!(ax, t2maps_norefcon["resnorm"][mask]; label = "no correction", color = (:blue, 0.5), scale_to = 0.9, offset = 1, bins = 50)
    hist!(ax, t2maps_refcon["resnorm"][mask]; label = "w/ correction", color = (:green, 0.5), scale_to = 0.9, offset = 0, bins = 50)
    hideydecorations!(ax, label = false)
    axislegend(ax)

    ax = Axis(fig[2,2]; xlabel = "T2 times [ms]", ylabel = "mean T2 dist. [a.u.]", xscale = log10)
    band!(ax, t2times, zero(t2times), globalmean(t2dist_norefcon); label = "no correction", linewidth = 3, color = (:blue, 0.5))
    band!(ax, t2times, zero(t2times), globalmean(t2dist_refcon); label = "w/ correction", linewidth = 3, color = (:green, 0.5))
    scatter!(ax, t2times, globalmean(t2dist_norefcon); markersize = 5, color = :blue)
    scatter!(ax, t2times, globalmean(t2dist_refcon); markersize = 5, color = :green)
    hideydecorations!(ax, label = false)
    axislegend(ax)

    return fig
end


text:
module JuliaSets
using PyPlot

export visJuliaSet, generatePlot

function visJuliaSet{T <: Real}(R::Function, x::Array{T, 1}, y::Array{T, 1},
	n_iter::Int64=1000, escape_tol::Float64=1e6)

	A = zeros(length(x), length(y))
	for i in 1:length(x)
  		for j in 1:length(y)
    		z = x[i] + im * y[j]
    		for k in 1:n_iter
      			z = R(z)
      			if abs(z) > escape_tol # diverged
        			A[i,j] = k
        			break
      			end
    		end
    		if A[i,j] == 0 # didn't diverge
      			A[i,j] = escape_tol + 1
    		end
  		end
	end
	A
end

function generatePlot(A::Array{Float64,2})
	matshow(A')
end

end # module

max_stars_repo_name: Jgmorton/JuliaSets

max_issues_repo_name: Jgmorton/JuliaSets

max_forks_repo_name: Jgmorton/JuliaSets


text:
# using Base.Test
#include("../src/RunTests.jl")

using PhysicalConstant, Measurements, Unitful
if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

using PhysicalConstant.CODATA2019

@testset begin
    @test ustrip(big(h)) == big"6.626070040e-34"
    @test setprecision(BigFloat, 768) do; precision(ustrip(big(c))) end == 768
    @test measurement(h) === measurement(h)
    @test iszero(measurement(α) - measurement(α))
    @test isone(measurement(BigFloat, atm) / measurement(BigFloat, atm))
    @test iszero(measurement(BigFloat, ħ) - (measurement(BigFloat, h) / 2big(pi)))
    @test isone(measurement(BigFloat, ħ) / (measurement(BigFloat, h) / 2big(pi)))
end


text:
module FieldArraysTests

using Test
using Gridap.Arrays
using Gridap.Fields
using Gridap.Fields: MockField, MockBasis
using Gridap.Fields: Valued, FieldOpKernel
using FillArrays
using Gridap.TensorValues

np = 4
p = Point(1,2)
x = fill(p,np)

v = 3.0
d = 2
f = MockField{d}(v)
fx = fill(v,np)
∇fx = fill(VectorValue(v,0.0),np)

l = 10
af = Fill(f,l)
ax = fill(x,l)
afx = fill(fx,l)
a∇fx = fill(∇fx,l)
test_array_of_fields(af,ax,afx,grad=a∇fx)

ag = apply_to_field_array(FieldOpKernel(+),af,af)
gx = fill(v+v,np)
∇gx = fill(VectorValue(v+v,0.0),np)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

struct FieldPlaceHolder <: Field end

ag = apply_to_field_array(FieldPlaceHolder,FieldOpKernel(+),af,af)
test_array(evaluate_field_array(ag,ax),agx)

c = 0.5
ac = fill(c,l)
ag = apply_to_field_array(FieldOpKernel(*),ac,af)
gx = fill(c*v,np)
∇gx = fill(VectorValue(c*v,0.0),np)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

w = 2.0
aw = fill(w,l)
ag = apply_to_field_array(FieldOpKernel(+),af,aw)
gx = fill(v+w,np)
∇gx = fill(VectorValue(v,0.0),np)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

ag = apply_to_field_array(FieldPlaceHolder,FieldOpKernel(+),af,aw)
test_array(evaluate_field_array(ag,ax),agx)

l = 10
af = Fill(f,l)
ax = Fill(x,l)
ag = apply_to_field_array(FieldOpKernel(+),af,af)
r1 = evaluate(ag,ax)
@test isa(r1,Fill)

np = 4
p = Point(1,2)
x = fill(p,np)
v = 2.0
d = 2
ndof = 8
wi = 3.0
w = fill(wi,ndof)
r = fill(v+wi,np,ndof)
f = MockBasis{d}(v,ndof)
∇fx = evaluate(∇(f),x)
af = Fill(f,l)
ax = fill(x,l)
aw = fill(w,l)
ag = apply_to_field_array(FieldOpKernel(+),af,aw)
agx = fill(r,l)
a∇gx = fill(∇fx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

v = 2.0
d = 2
wi = 3.0
w = fill(wi,ndof)
r = fill(v+wi,np,ndof)
f = MockField{d}(v)
∇r = fill(VectorValue(v,0.0),np,ndof)
af = Fill(f,l)
ax = fill(x,l)
aw = fill(w,l)
ag = apply_to_field_array(FieldOpKernel(+),af,aw)
agx = fill(r,l)
a∇gx = fill(∇r,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

# lazy_append

np = 4
p = Point(1,2)
x = fill(p,np)

v = 3.0
d = 2
f = MockField{d}(v)
fx = evaluate(f,x)
∇fx = evaluate(∇(f),x)

l = 10
af = Fill(f,l)
ax = fill(x,l)
afx = fill(fx,l)
a∇fx = fill(∇fx,l)

np = 5
p = Point(4,3)
x = fill(p,np)

v = 5.0
d = 2
f = MockField{d}(v)
fx = evaluate(f,x)
∇fx = evaluate(∇(f),x)

l = 15
bf = Fill(f,l)
bx = fill(x,l)
bfx = fill(fx,l)
b∇fx = fill(∇fx,l)

cf = lazy_append(af,bf)
cx = lazy_append(ax,bx)

cfx = evaluate(cf,cx)
∇cfx = evaluate(∇(cf),cx)
rfx = vcat(afx,bfx)
∇rfx = vcat(a∇fx,b∇fx)
test_array_of_fields(cf,cx,rfx,grad=∇rfx)
@test isa(cfx, AppendedArray)
@test isa(∇cfx, AppendedArray)

end # module


text:
"""
FHIR version R4.

## Summary
abstract type Authentication <: Any
"""
struct R4 <: FHIRVersion
end

module R4Types

import Dates
import JSON3
import SaferIntegers
import StructTypes
import TimeZones

include("autogenerated/autogenerated-R4.jl")

end # end module R4Types


text: println("testing suite is hooked up")


text:
# Write a for loop that counts backward using Julia’s range notation.
for i = 10:-1:1
  println(i)
end

#Write an iteration over a multidimensional array like [1 2 3; 4 5 6; 7 8 9]. In what order does it get printed out?
for i = [1 2 3; 4 5 6; 7 8 9]
  println(i)
end

# Write a factorial function as a parallel for loop.
function factorial(n)
  @parallel (*) for i = 1:n
    i
  end
end

# Add a method for concat that can concatenate an integer with a matrix. concat(5,[1 2 ; 3 4]) should produce [ 5 5 1 2 ; 5 5 3 4].
function concat(a :: Int64, b :: Array)
  hcat(fill(a, (size(b)[1], 2)), b)
end

# You can extend built-in functions with new methods too. Add a new method for + to make "jul" + "ia" work.
function +(a :: ASCIIString, b :: ASCIIString)
  "$a$b"
end

max_stars_repo_name: hemslo/seven-in-seven

max_issues_repo_name: hemslo/seven-in-seven

max_forks_repo_name: hemslo/seven-in-seven


text:
include("imports.jl")

@testset ExtendedTestSet "`NormalizeRow`" begin
    cols = [:col1, :col2, :col3]
    item = TabularItem(NamedTuple(zip(cols, [1, "a", 10])), cols)
    cols_to_normalize = [:col1, :col3]
    col1_mean, col1_std = 10, 100
    col3_mean, col3_std = 100, 10
    normdict = Dict(:col1 => (col1_mean, col1_std), :col3 => (col3_mean, col3_std))

    tfm = NormalizeRow(normdict, cols_to_normalize)
    testapply(tfm, item)
    titem = apply(tfm, item)
    @test titem.data[:col1] == (item.data[:col1] - col1_mean)/col1_std
    @test titem.data[:col3] == (item.data[:col3] - col3_mean)/col3_std
end

@testset ExtendedTestSet "`FillMissing`" begin
    cols = [:col1, :col2, :col3]
    item = TabularItem(NamedTuple(zip(cols, [1, missing, missing])), cols)
    cols_to_fill = [:col1, :col3]
    col1_fmval = 1000.
    col2_fmval = "d"
    col3_fmval = 1000.
    fmdict = Dict(:col1 => col1_fmval, :col2 => col2_fmval, :col3 => col3_fmval)

    tfm1 = FillMissing(fmdict, cols_to_fill)
    @test_nowarn apply(tfm1, item)
    titem = apply(tfm1, item)
    @test titem.data[:col1] == coalesce(item.data[:col1], col1_fmval)
    @test titem.data[:col3] == coalesce(item.data[:col3], col3_fmval)
    @test ismissing(titem.data[:col2])

    push!(cols_to_fill, :col2)
    tfm2 = FillMissing(fmdict, cols_to_fill)
    testapply(tfm2, item)
    titem2 = apply(tfm2, item)
    @test titem2.data[:col2] == coalesce(item.data[:col2], "d")
end

@testset ExtendedTestSet "`Categorify`" begin
    cols = [:col1, :col2, :col3, :col4]
    item = TabularItem(NamedTuple(zip(cols, [1, "a", "A", missing])), cols)
    cols_to_categorify = [:col2, :col3, :col4]

    categorydict = Dict(:col2 => ["a", "b", "c"], :col3 => ["C", "B", "A"], :col4 => [missing, 10, 20])
    tfm = Categorify(categorydict, cols_to_categorify)
    @test !any(ismissing, tfm.dict[:col4])
    @test_nowarn apply(tfm, item)
    testapply(tfm, item)
    titem = apply(tfm, item)
    @test titem.data[:col2] == 2
    @test titem.data[:col3] == 4
    @test titem.data[:col4] == 1
end


text:
Base.IndexStyle(::Type{<:SimpleSDMLayer}) = IndexCartesian()

"""
    Base.CartesianIndices(layer::T) where {T <: SimpleSDMLayer}

This function is used to access the positions in a layer as cartesian indices.
Internally, although it may be more convenient to access positions by their
coordinates, the value of the raster are extracted using their cartesian
coordinates.
"""
function Base.CartesianIndices(layer::T) where {T <: SimpleSDMLayer}
    return CartesianIndices(layer.grid)
end

"""
    Base.getindex(layer::T, i::CartesianIndex{2}) where {T <: SimpleSDMLayer}

Returns the value stored at a given cartesian index.
"""
function Base.getindex(layer::T, i::CartesianIndex{2}) where {T <: SimpleSDMLayer}
    return layer.grid[i]
end

"""
    Base.getindex(layer::T, i::Integer) where {T <: SimpleSDMLayer}

Returns the value stored at a linear index (a good candidate for deprecation as
we have a better iteration interface now).
"""
function Base.getindex(layer::T, i::Integer) where {T <: SimpleSDMLayer}
    return layer.grid[CartesianIndices(layer)[i]]
end

"""
    Base.getindex(layer::T, i::Integer, j::Integer) where {T <: SimpleSDMLayer}

Standard abstract array accession, where the dimensions follow the dimensions of
the underlying grid.
"""
function Base.getindex(layer::T, i::Integer, j::Integer) where {T <: SimpleSDMLayer}
    return layer.grid[CartesianIndex(i,j)]
end

"""
    Base.getindex(layer::T, i::Array{CartesianIndex{2}}) where {T <: SimpleSDMLayer}

Returns an array of values based on an array of cartesian indices - this can be
a vector or a matrix, and the elements can be in any order. This will *not*
return a raster.
"""
function Base.getindex(layer::T, i::Array{CartesianIndex{2}}) where {T <: SimpleSDMLayer}
    return layer.grid[i]
end

"""
    Base.getindex(layer::T, longitude::AbstractFloat, latitude::AbstractFloat) where {T <: SimpleSDMLayer}

Returns a value by longitude and latitude.
"""
function Base.getindex(layer::T, longitude::AbstractFloat, latitude::AbstractFloat) where {T <: SimpleSDMLayer}
    return layer[Point(longitude,latitude)]
end

"""
    Base.getindex(layer::T, c::Point) where {T <: SimpleSDMLayer}

Access a value through a `Point` (from `GeometryBasics`), which has the
longitude first and the latitude last -- this follows the GeoJSON convention.
"""
function Base.getindex(layer::T, c::Point) where {T <: SimpleSDMLayer}
    return layer[_point_to_cartesian(layer, c)]
end

"""
    Base.getindex(layer::T, c::Array{<:Point}) where {T <: SimpleSDMLayer}

Access a value through an array of `Point` (from `GeometryBasics`), which has
the longitude first and the latitude last. The array can be in any order, so
this method will not return a raster.
"""
function Base.getindex(layer::T, c::Array{<:Point}) where {T <: SimpleSDMLayer}
    return layer[[_point_to_cartesian(layer, i) for i in c]]
end

"""
    Base.getindex(::T, ::Nothing) where {T <: SimpleSDMLayer}

If the user requests a point that is out of bounds, its cartesian coordinate
will be matched to `nothing`, and then we return `nothing`.
"""
function Base.getindex(::T, ::Nothing) where {T <: SimpleSDMLayer}
    return nothing
end

function Base.setindex!(layer::T, v, i::CartesianIndex{2}) where {T <: SimpleSDMResponse}
    return setindex!(layer.grid, v, i)
end

function Base.setindex!(layer::T, v, i::Array{CartesianIndex{2}}) where {T <: SimpleSDMResponse}
    return setindex!(layer.grid, v, i)
end

function Base.setindex!(layer::T, v, point::Point) where {T <: SimpleSDMResponse}
    return setindex!(layer, v, _point_to_cartesian(layer, point))
end

function Base.setindex!(layer::T, v, points::Array{<:Point}) where {T <: SimpleSDMResponse}
    i = broadcast(point -> _point_to_cartesian(layer, point), points)
    return setindex!(layer, v, i)
end

function Base.setindex!(layer::T, v, longitude::Float64, latitude::Float64) where {T <: SimpleSDMResponse}
    return setindex!(layer, v, Point(longitude, latitude))
end


text:
#############
# Exp & Log #
#############

struct Exp{N} <: Bijector{N} end
struct Log{N} <: Bijector{N} end
up1(::Exp{N}) where {N} = Exp{N + 1}()
up1(::Log{N}) where {N} = Log{N + 1}()

Exp() = Exp{0}()
Log() = Log{0}()

(b::Exp{0})(y::Real) = exp(y)
(b::Log{0})(x::Real) = log(x)

(b::Exp{0})(y::AbstractArray{<:Real}) = exp.(y)
(b::Log{0})(x::AbstractArray{<:Real}) = log.(x)

(b::Exp{1})(y::AbstractVector{<:Real}) = exp.(y)
(b::Exp{1})(y::AbstractMatrix{<:Real}) = exp.(y)
(b::Log{1})(x::AbstractVector{<:Real}) = log.(x)
(b::Log{1})(x::AbstractMatrix{<:Real}) = log.(x)

(b::Exp{2})(y::AbstractMatrix{<:Real}) = exp.(y)
(b::Log{2})(x::AbstractMatrix{<:Real}) = log.(x)

(b::Exp{2})(y::AbstractArray{<:AbstractMatrix{<:Real}}) = map(b, y)
(b::Log{2})(x::AbstractArray{<:AbstractMatrix{<:Real}}) = map(b, x)

inv(b::Exp{N}) where {N} = Log{N}()
inv(b::Log{N}) where {N} = Exp{N}()

logabsdetjac(b::Exp{0}, x::Real) = x
logabsdetjac(b::Exp{0}, x::AbstractVector) = x
logabsdetjac(b::Exp{1}, x::AbstractVector) = sum(x)
logabsdetjac(b::Exp{1}, x::AbstractMatrix) = vec(sum(x; dims = 1))
logabsdetjac(b::Exp{2}, x::AbstractMatrix) = sum(x)
logabsdetjac(b::Exp{2}, x::AbstractArray{<:AbstractMatrix}) = map(x) do x
    logabsdetjac(b, x)
end

logabsdetjac(b::Log{0}, x::Real) = -log(x)
logabsdetjac(b::Log{0}, x::AbstractVector) = .-log.(x)
logabsdetjac(b::Log{1}, x::AbstractVector) = - sum(log, x)
logabsdetjac(b::Log{1}, x::AbstractMatrix) = - vec(sum(log, x; dims = 1))
logabsdetjac(b::Log{2}, x::AbstractMatrix) = - sum(log, x)
logabsdetjac(b::Log{2}, x::AbstractArray{<:AbstractMatrix}) = map(x) do x
    logabsdetjac(b, x)
end