text:
# sampling functions

"""
    sample_trajectory!(x, sampler; options=MDOptions())

In place applciation of the `sampler` to `x`.  Number of iterations are
set using the `options` argument.

### Fields

* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
* `constraints` - Constraints on the trajectory
### Optional Fields
* `options`   - Sampling options, including number of iteration
"""
function sample_trajectory!(x::Tx, sampler::S, constraints::C; options = MDOptions()) where {Tx,S<:AbstractSampler,C<:AbstractConstraints}

    state = InitState!(x, sampler)
    for i = 1:options.n_iters
        if (mod(i - 1, constraints.n_before) == 0)
            constraints.before_update!(state)
        end
        UpdateState!(state, sampler)
        if (mod(i, constraints.n_after) == 0)
            constraints.after_update!(state)
        end
        # constraint!(state, i)
    end
    x
end

"""
    sample_trajectory(x₀, sampler; options=MDOptions())

Run the `sampler` starting at `x₀`.  Number of iterations and interval between
saves are set using the `options` argument.  For Metropolis samplers, the
running acceptance rates are also resturned.

### Fields
* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
* `constraints` - Constraints on the trajectory
### Optional Fields
* `options`   - Sampling options, including number of iteration
"""
function sample_trajectory(x₀::Tx, sampler::S, constraints::C; options = MDOptions()) where {Tx,S<:MetropolisSampler,C<:AbstractConstraints}

    n_accept = Int(0)

    state = InitState(x₀, sampler)

    # allocate memory for samples
    samples = Tx[similar(x₀) for i = 1:options.n_save]
    acceptance_rates = zeros(options.n_save)
    save_index = 1
    for i = 1:options.n_iters
        if (mod(i - 1, constraints.n_before) == 0)
            constraints.before_update!(state)
        end
        UpdateState!(state, sampler)
        if (mod(i, constraints.n_after) == 0)
            constraints.after_update!(state)
        end
        n_accept += state.accept
        if (mod(i, options.n_save_iters) == 0)
            @. samples[save_index] = deepcopy(state.x)
            acceptance_rates[save_index] = n_accept / i
            save_index += 1
        end
    end
    return samples, acceptance_rates
end

function sample_trajectory(x₀::Tx, sampler::S, constraints::C; options = MDOptions()) where {Tx,S<:NonMetropolisSampler,C<:AbstractConstraints}

    state = InitState(x₀, sampler)
    samples = Tx[similar(x₀) for i = 1:options.n_save]
    save_index = 1
    for i = 1:options.n_iters
        if (mod(i - 1, constraints.n_before) == 0)
            constraints.before_update!(state)
        end
        UpdateState!(state, sampler)
        if (mod(i, constraints.n_after) == 0)
            constraints.after_update!(state)
        end
        if (mod(i, options.n_save_iters) == 0)
            @. samples[save_index] = deepcopy(state.x)
            save_index += 1
        end
    end
    return samples
end

"""
    sample_observables(x₀, sampler, observables; options=MDOptions())

Run the `sampler` starting at `x₀`, evaluating the trajectory on a tuple of
`observables` scalar functions.  Number of iterations and interval between saves
are set using the `options` argument.  Only the computed observables are returned.

### Fields
* `x`         - Starting position for sampler, modified in place
* `sampler`   - Desired sampler
* `constraints` - Constraints on the trajectory
* `observables` - Observables on which to evaluate the trajectory
### Optional Fields
* `TO`- Observable data type, if needed, should be entered as the first argument
* `options`   - Sampling options, including number of iteration
"""
function sample_observables(x₀::Tx, sampler::S, constraints::C, observables::Tuple{Vararg{<:Function,NO}};
    options = MDOptions()) where {Tx,S<:AbstractSampler,NO,C<:AbstractConstraints}

    state = InitState(x₀, sampler)
    observable_samples = zeros(NO, options.n_save)
    save_index = 1
    for i = 1:options.n_iters
        if (mod(i - 1, constraints.n_before) == 0)
            constraints.before_update!(state)
        end
        UpdateState!(state, sampler)
        if (mod(i, constraints.n_after) == 0)
            constraints.after_update!(state)
        end
        if (mod(i, options.n_save_iters) == 0)
            ntuple(k -> observable_samples[k, save_index] = (observables[k])(state.x), NO)
            save_index += 1
        end
    end
    return observable_samples

end

function sample_observables(TO::Type, x₀::Tx, sampler::S, constraints::C, observables::Tuple{Vararg{<:Function,NO}};
    options = MDOptions()) where {Tx,S<:AbstractSampler,NO,C<:AbstractConstraints}

    state = InitState(x₀, sampler)
    observable_samples = zeros(TO, NO, options.n_save)
    save_index = 1
    for i = 1:options.n_iters
        if (mod(i - 1, constraints.n_before) == 0)
            constraints.before_update!(state)
        end
        UpdateState!(state, sampler)
        if (mod(i, constraints.n_after) == 0)
            constraints.after_update!(state)
        end
        if (mod(i, options.n_save_iters) == 0)
            ntuple(k -> observable_samples[k, save_index] = (observables[k])(state.x), NO)
            save_index += 1
        end
    end
    return observable_samples
end


text:
# 1D arrays stored as Vectors

# Uninitialized Arrays

arr = [];

# Typed Arrays

arr = Int64[];

# Initialized Arrays

arr = [1, 2, 3, 4];
arr = ['a', 1, 2, 3.4];

# Like strings, indexing starts at 1
# - Trying to access invalid indexes will throw BoundsError

# arr[0];
arr[1];
arr[4];
# arr[5];

# Array concatenation

# Commas
arr = [1:2, 3:4]; # creates array of 2 elements, each a UnitRange
println(arr);

# Spaces or tabs
arr = [1:2 3:4]; # creates 2D array, each space creates new column
println(arr);

println(arr[1][1]);

# Semicolons or newlines
arr = [1:2; 3:4]; # creates array of 4 elements
println(arr);
arr = [1:2
        3:4];
println(arr);

# using both
arr = [1 2; 3 4];
println(arr);

# Matrices

mat = [1 2 3; 4 5 6];

# Matrix access by mat[row, col] notation
mat[2, 3];

# Multidimensional arrays

muldimarr = zeros(2, 2, 3);
println(ndims(muldimarr));

# Array slicing

arr = [1, 2, 3, 4, 5, 6, 7];
sliceArr = arr[3:5];
println(sliceArr);

mat = [1 2 3; 4 5 6];
sliceMat = mat[1:2, 2:3];
println(sliceMat);

# Copies and using arrays

a = [1, 2, 3];
# Pointer
b = a;

b[2] = 4;
println(a);
println(b);

# Copy
c = copy(a);
c[2] = 2;
println(a);
println(c);

ext: jl

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
abstract type ElectricLoad <: StaticInjection end
abstract type StaticLoad <: ElectricLoad end
abstract type ControllableLoad <: ElectricLoad end


text:
export StraightLineHomotopy

"""
    StraightLineHomotopy(G::System, F::System; gamma = 1.0)
    StraightLineHomotopy(G::AbstractSystem, F::AbstractSystem; gamma = 1.0)

Constructs the straight line homotopy ``H(x, t) = γ t G(x) + (1-t) F(x)`` where
``γ`` is `gamma`.
"""
struct StraightLineHomotopy{S<:AbstractSystem,T<:AbstractSystem} <: AbstractHomotopy
    start::S
    target::T
    γ::ComplexF64

    u::Vector{ComplexF64}
    ū::Vector{ComplexDF64}
    v̄::Vector{ComplexDF64}
    U::Matrix{ComplexF64}

    dv_start::TaylorVector{4,ComplexF64}
    dv_target::TaylorVector{4,ComplexF64}
end

function StraightLineHomotopy(
    start::System,
    target::System;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
)
    StraightLineHomotopy(
        fixed(start; compile = compile),
        fixed(target; compile = compile);
        kwargs...,
    )
end
function StraightLineHomotopy(
    start::AbstractSystem,
    target::AbstractSystem;
    γ = 1.0,
    gamma = γ,
)
    size(start) == size(target) || throw(
        ArgumentError(
            "Start and target do not have the same size, got $(size(start)) and $(size(target))",
        ),
    )

    m, n = size(start)
    u = zeros(ComplexF64, m)
    ū = zeros(ComplexDF64, m)
    v̄ = zeros(ComplexDF64, m)
    U = zeros(ComplexF64, m, n)

    dv_start = TaylorVector{4}(ComplexF64, m)
    dv_target = TaylorVector{4}(ComplexF64, m)

    StraightLineHomotopy(
        start,
        target,
        ComplexF64(gamma),
        u,
        ū,
        v̄,
        U,
        dv_start,
        dv_target,
    )
end

Base.size(H::StraightLineHomotopy) = size(H.start)

function Base.show(io::IO, mime::MIME"text/plain", H::StraightLineHomotopy)
    println(io, typeof(H), ":")
    println(io, "γ: ", H.γ)
    println(io, "\nG: ")
    show(io, mime, H.start)
    println(io, "\n\nF:")
    show(io, mime, H.target)
end

function ModelKit.evaluate!(
    u,
    H::StraightLineHomotopy,
    x::Vector{ComplexDF64},
    t,
    p = nothing,
) where {T}
    evaluate!(H.v̄, H.start, x)
    evaluate!(H.ū, H.target, x)
    ts, tt = H.γ .* t, 1 - t
    for i = 1:size(H, 1)
        @inbounds u[i] = ts * H.v̄[i] + tt * H.ū[i]
    end
end

function ModelKit.evaluate!(u, H::StraightLineHomotopy, x, t, p = nothing)
    evaluate!(u, H.start, x)
    evaluate!(H.u, H.target, x)
    ts, tt = H.γ .* t, 1 - t
    for i = 1:size(H, 1)
        @inbounds u[i] = ts * u[i] + tt * H.u[i]
    end
    u
end

function ModelKit.evaluate_and_jacobian!(
    u,
    U,
    H::StraightLineHomotopy,
    x::AbstractVector{T},
    t,
    p = nothing,
) where {T}
    evaluate_and_jacobian!(u, U, H.start, x)
    evaluate_and_jacobian!(H.u, H.U, H.target, x)
    ts, tt = H.γ .* t, 1 - t
    for i = 1:size(H, 1)
        @inbounds u[i] = ts * u[i] + tt * H.u[i]
    end
    for j = 1:size(H, 2), i = 1:size(H, 1)
        @inbounds U[i, j] = ts * U[i, j] + tt * H.U[i, j]
    end
    nothing
end
function ModelKit.taylor!(u, ::Val{1}, H::StraightLineHomotopy, x, t)
    evaluate!(u, H.start, x)
    evaluate!(H.u, H.target, x)
    for i = 1:size(H, 1)
        @inbounds u[i] = H.γ * u[i] - H.u[i]
    end
    u
end
function ModelKit.taylor!(
    u,
    v::Val{K},
    H::StraightLineHomotopy,
    tx::TaylorVector,
    t,
) where {K}
    taylor!(H.dv_start, v, H.start, tx)
    taylor!(H.dv_target, v, H.target, tx)

    for i = 1:size(H, 1)
        start = H.γ * (H.dv_start[i, K] + t * H.dv_start[i, K+1])
        target = (1 - t) * H.dv_target[i, K+1] - H.dv_target[i, K]
        u[i] = start + target
    end
    u
end


text:
f0 = Complex(0.0, 0.0)
f1 = Complex(1.0, 0.0)
f1i = Complex(0.0, 1.0)
f11 = Complex(1.0, 1.0)
f2h = Complex(2.0, 0.5)

d0 = Complex(Double64(0), Double64(0))
d1 = Complex(Double64(1), Double64(0))
d1i = Complex(Double64(0), Double64(1))
d11 = Complex(Double64(1), Double64(1))
d2h = Complex(Double64(2), Double64(0.5))

@testset "Complex" begin
  @test isapprox(log(f11), Complex{Float64}(log(d11)))
  @test isapprox(exp(f11), Complex{Float64}(exp(d11)))

  @test isapprox(sin(f11), Complex{Float64}(sin(d11)))
  @test isapprox(cos(f11), Complex{Float64}(cos(d11)))
  @test isapprox(tan(f11), Complex{Float64}(tan(d11)))
  @test isapprox(csc(f11), Complex{Float64}(csc(d11)))
  @test isapprox(sec(f11), Complex{Float64}(sec(d11)))
  @test isapprox(cot(f11), Complex{Float64}(cot(d11)))

  @test isapprox(sinh(f11), Complex{Float64}(sinh(d11)))
  @test isapprox(cosh(f11), Complex{Float64}(cosh(d11)))
  @test isapprox(tanh(f11), Complex{Float64}(tanh(d11)))
  @test isapprox(csch(f11), Complex{Float64}(csch(d11)))
  @test isapprox(sech(f11), Complex{Float64}(sech(d11)))
  @test isapprox(coth(f11), Complex{Float64}(coth(d11)))

  @test isapprox(asin(f11), Complex{Float64}(asin(d11)))
  @test isapprox(acos(f11), Complex{Float64}(acos(d11)))
  @test isapprox(atan(f2h), Complex{Float64}(atan(d2h)))
  @test isapprox(acsc(f11), Complex{Float64}(acsc(d11)))
  @test isapprox(asec(f11), Complex{Float64}(asec(d11)))
  @test isapprox(acot(f2h), Complex{Float64}(acot(d2h)))

  @test isapprox(asinh(f11), Complex{Float64}(asinh(d11)))
  @test isapprox(acosh(f11), Complex{Float64}(acosh(d11)))
  @test isapprox(atanh(f11), Complex{Float64}(atanh(d11)))
  @test isapprox(acsch(f11), Complex{Float64}(acsch(d11)))
  @test isapprox(asech(f11), Complex{Float64}(asech(d11)))
  @test isapprox(acoth(f11), Complex{Float64}(acoth(d11)))
end

@testset "Complex Accuracy" begin
  x = ComplexDF64(sqrt(df64"2"), cbrt(df64"3"))
  y = Complex{BigFloat}(x)

  @test abs(sin(y) - sin(x)) < 1.0e-31
  @test abs(cos(y) - cos(x)) < 1.0e-31
  @test abs(tan(y) - tan(x)) < 1.0e-31
  @test abs(csc(y) - csc(x)) < 1.0e-31
  @test abs(sec(y) - sec(x)) < 1.0e-31
  @test abs(tan(y) - tan(x)) < 1.0e-31

  @test abs(sinh(y) - sinh(x)) < 1.0e-31
  @test abs(cosh(y) - cosh(x)) < 1.0e-31
  @test abs(tanh(y) - tanh(x)) < 1.0e-31
  @test abs(csch(y) - csch(x)) < 1.0e-31
  @test abs(sech(y) - sech(x)) < 1.0e-31
  @test abs(tanh(y) - tanh(x)) < 1.0e-31

  @test abs(asin(y) - asin(x)) < 1.0e-31
  @test abs(acos(y) - acos(x)) < 1.0e-31
  @test abs(atan(y) - atan(x)) < 1.0e-31
  @test abs(acsc(y) - acsc(x)) < 1.0e-31
  @test abs(asec(y) - asec(x)) < 1.0e-31
  @test abs(atan(y) - atan(x)) < 1.0e-31

  @test abs(asinh(y) - asinh(x)) < 1.0e-31
  @test abs(acosh(y) - acosh(x)) < 1.0e-31
  @test abs(atanh(y) - atanh(x)) < 1.0e-31
  @test abs(acsch(y) - acsch(x)) < 1.0e-31
  @test abs(asech(y) - asech(x)) < 1.0e-31
  @test abs(atanh(y) - atanh(x)) < 1.0e-31

end

@testset "Complex exp(huge)" begin
    x = Double64(-1.168869267851442e79) + Double64(4.140780141319954e87)*im
    @test ComplexDF64(exp(x)) == exp(ComplexDF64(x))
end


text:
module TestZygote

if VERSION ≥ v"1.6-"
    try
        using Zygote
    catch
        @info "Failed to import Zygote. Trying again with `@stdlib`..."
        push!(LOAD_PATH, "@stdlib")
    end
end

include("preamble.jl")
using Zygote

include("../examples/zygote.jl")
@test g isa Matrix
@test size(g) == (d, d)

end  # module


text:
without_basesize(; basesize = nothing, kw...) = kw

function ThreadsX.reduce(op, itr; kw...)
    xf, reducible = extract_transducer(itr)
    result = foldxt(
        op,
        xf,
        reducible;
        init = Init(op),
        basesize = default_basesize(reducible),
        kw...,
    )
    result === Init(op) && return reduce_empty(op, eltype(itr))
    return result
end

function ThreadsX.mapreduce(f, op, itr; kw...)
    xf, reducible = extract_transducer(itr)
    result = foldxt(
        op,
        opcompose(xf, Map(f)),
        reducible;
        init = Init(op),
        basesize = default_basesize(reducible),
        kw...,
    )
    result === Init(op) && return mapreduce_empty(f, op, eltype(itr))
    return result
end

function ThreadsX.mapreduce(f, op, itr, itrs...; kw...)
    if isempty(itr)
        # `Base` just does `reduce(op, map(f, ...))`:
        return mapreduce(f, op, itr, itrs...; without_basesize(; kw...)...)
    end
    return foldxt(
        op,
        MapSplat(f),
        zip(itr, itrs...);
        basesize = default_basesize(itr),
        kw...,
    )
end

# Maybe refactor the function based on mapreduce into AbstractReducers.jl?

ThreadsX.sum(itr; kw...) = ThreadsX.sum(identity, itr; kw...)
ThreadsX.sum(f, itr; kw...) = ThreadsX.mapreduce(f, add_sum, itr; simd = Val(true), kw...)

ThreadsX.prod(itr; kw...) = ThreadsX.prod(identity, itr; kw...)
ThreadsX.prod(f, itr; kw...) = ThreadsX.mapreduce(f, mul_prod, itr; simd = Val(true), kw...)

ThreadsX.count(itr; kw...) = ThreadsX.count(identity, itr; kw...)
ThreadsX.count(f, itr; kw...) =
    ThreadsX.sum(x -> Int(f(x)::Bool), itr; init = 0, simd = Val(true), kw...)

ThreadsX.maximum(itr; kw...) = ThreadsX.maximum(identity, itr; kw...)
ThreadsX.maximum(f, itr; kw...) = ThreadsX.mapreduce(f, max, itr; simd = Val(true), kw...)

ThreadsX.minimum(itr; kw...) = ThreadsX.minimum(identity, itr; kw...)
ThreadsX.minimum(f, itr; kw...) = ThreadsX.mapreduce(f, min, itr; simd = Val(true), kw...)

asbool(f) = x -> f(x)::Bool

# TODO: `any` and `all` should be done with "unordered" version
ThreadsX.any(itr; kw...) = ThreadsX.any(identity, itr; kw...)
ThreadsX.any(f, itr; kw...) = foldxt(
    right,  # no need to use `|`
    opcompose(Map(asbool(f)), ReduceIf(identity)),
    simd = Val(true),
    basesize = default_basesize(ThreadsX.any, f, itr),
    itr;
    kw...,
    init = false,
)

ThreadsX.all(itr; kw...) = ThreadsX.all(identity, itr; kw...)
ThreadsX.all(f, itr; kw...) = foldxt(
    right,  # no need to use `&`
    opcompose(Map(asbool(f)), ReduceIf(!)),
    itr;
    simd = Val(true),
    basesize = default_basesize(ThreadsX.all, f, itr),
    kw...,
    init = true,
)

ThreadsX.findfirst(itr; kw...) = ThreadsX.findfirst(identity, itr; kw...)
ThreadsX.findfirst(f, array::AbstractArray; kw...) = foldxt(
    right,
    ReduceIf(i -> f(@inbounds array[i])),
    keys(array);
    init = nothing,
    simd = Val(true),
    basesize = default_basesize(ThreadsX.findfirst, f, array),
    kw...,
)

ThreadsX.findlast(itr; kw...) = ThreadsX.findlast(identity, itr; kw...)
function ThreadsX.findlast(f, array::AbstractArray; kw...)
    idx = keys(array)
    return foldxt(
        right,
        opcompose(Map(i -> (@inbounds idx[i])), ReduceIf(i -> f(@inbounds array[i]))),
        lastindex(idx):-1:firstindex(idx);
        init = nothing,
        simd = Val(true),
        basesize = default_basesize(ThreadsX.findlast, f, array),
        kw...,
    )
end

ThreadsX.findall(itr; kw...) = ThreadsX.findall(identity, itr; kw...)
function ThreadsX.findall(f, array::AbstractArray; kw...)
    idxs = tcollect(
        Filter(i -> f(@inbounds array[i])),
        keys(array);
        simd = Val(true),
        kw...,
    )
    isempty(idxs) && return keytype(array)[]
    return idxs
end

_minmax((min0, max0), (min1, max1)) = (min(min0, min1), max(max0, max1))

ThreadsX.extrema(itr; kw...) = ThreadsX.extrema(identity, itr; kw...)
ThreadsX.extrema(f, itr; kw...) = foldxt(
    asmonoid(_minmax),
    Map(x -> (y = f(x); (y, y))),
    itr;
    simd = Val(true),
    basesize = default_basesize(ThreadsX.extrema, f, itr),
    kw...,
)

ThreadsX.issorted(
    itr;
    lt = isless,
    by = identity,
    rev::Union{Bool,Nothing} = nothing,
    order::Base.Ordering = Base.Forward,
    kw...,
) = ThreadsX.issorted(itr, Base.ord(lt, by, rev, order); kw...)

ThreadsX.issorted(itr, order::Base.Ordering; kw...) =
    ThreadsX.all(
        zip(
            view(itr, firstindex(itr):lastindex(itr)-1),
            view(itr, firstindex(itr)+1:lastindex(itr)),
        );
        kw...,
    ) do (x1, x2)
        !Base.lt(order, x2, x1)
    end

struct PushUnique{F} <: Function
    f::F
end
PushUnique(::Type{T}) where {T} = PushUnique{Type{T}}(T)

function (f!::PushUnique)((ys, seen), x)
    fx = f!.f(x)
    return fx in seen ? (ys, seen) : (push!!(ys, x), push!!(seen, fx))
end

# TODO: do this with public API of Transducers or make it public
function Transducers.combine(f!::PushUnique, (ys1, seen1), (ys2, seen2))
    seen3 = setdiff!(seen2, seen1)
    isempty(seen3) && return (ys1, seen1)
    return (append!!(Filter(x -> f!.f(x) in seen3), ys1, ys2), union!!(seen1, seen2))
end
# * Add an option to avoid re-compute `f(x)` in combine?
# * Iterate over `seen3` if `length(seen3) << length(ys2)`?

# Manually create a singleton callable since closure captures the
# types as `DataType` and causes type-instability:
struct InitUnique{X,Y} <: Function end
@inline function (::InitUnique{X,Y})() where {X,Y}
    if isbitstype(Y) || Base.isbitsunion(Y)
        return (X[], Set{Y}())
    else
        return (X[], Empty(Set))
    end
end

ThreadsX.unique(itr; kw...) = ThreadsX.unique(identity, itr; kw...)
function ThreadsX.unique(f::F, itr::AbstractVector{X}; kw...) where {F,X}
    # Using inference as an optimization. The result of this inference
    # does not affect the result:
    Y = Core.Compiler.return_type(f, Tuple{X})
    ys, = foldxt(
        PushUnique(f),
        Map(identity),
        itr;
        basesize = default_basesize(itr),
        kw...,
        init = OnInit(InitUnique{X,Y}()),
    )
    return ys::Vector{X}
end

ThreadsX.Set(itr; kw...) = tcopy(Set, itr; kw...)


text:
using Graphics

export
          # Structural:
          Font,
          # General utility:
          Point, Square,
          # Drawable Shapes:
          NBox, Circle, Arc, TextLine,
          # Constructors:
          Border, BoxOutline
# ======================================================================================
# ======================================================================================
abstract type Geo end

#=mutable struct Point
    x::Float64
    y::Float64
    Point(x,y) = new(x,y)
end=#

#-------------------------------------------------------------------------------
# currently not being used..
#-------------------------------------------------------------------------------
macro import_fields(t)
  tt = eval(t)
  fields = fieldnames(tt)
  ex = :()
  for i = 1:length(fields)
    ft = fieldtype(tt, fields[i])
    if i==1
      ex = :($(fields[i])::$(ft))
    else
      ex = :($ex ; $(fields[i]) :: $(ft))
    end
  end
  return ex
end
#-------------------------------------------------------------------------------
mutable struct BoxOutline <: Geo
    left::Float64
    top::Float64
    right::Float64
    bottom::Float64
    width::Float64
    height::Float64
    BoxOutline(left,top,right,bottom,width,height) = new(left,top,right,bottom,width,height)
    BoxOutline() = new(0,0,0,0,0,0)
end
mutable struct Square <: Geo
    left::Float64
    top::Float64
    width::Float64
    height::Float64
end

#-------------------------------------------------------------------------------
abstract type Draw <: Geo end
mutable struct Border <: Geo
    left::Float64
    top::Float64
    right::Float64
    bottom::Float64
    width::Float64
    height::Float64

    style::Any # maybe AbstractString
    color::Vector{Float32} # this may be an array of arrays in the case that each side has a different color
    radius::Union{Nothing,Array} # Nullable{Array}
    Border(left,top, right,bottom, width,height, style,color,radius) = new(left,top, right,bottom, width,height, style,color,radius)
    Border() = new(0,0,0,0,0,0, 0,[],[0,0,0,0])
end

import Base.get
get(border::Union{Nothing,Border}, default::Border) = border == nothing ? default : border
get(marg_or_pad::Union{Nothing,BoxOutline}, default::BoxOutline) = marg_or_pad == nothing ? default : marg_or_pad
# get(padding::) = padding == nothing ? BoxOutline(0,0,0,0,0,0) : padding
# ==============================================================================  <: Shape

mutable struct BasicShape <: Draw
    flags::BitArray{1}
    color::Vector{Any}
    gradient::Vector{Any}
    opacity::Any
    padding::Union{Nothing,BoxOutline}
    border::Union{Nothing,Border}
    margin::Union{Nothing,BoxOutline}
    offset::Union{Nothing,Point}
    left::Float64
    top::Float64
    width::Float64
    height::Float64
    BasicShape() = new(falses(64), [], [], 1, nothing, nothing, nothing, nothing,0,0,0,0)
end
function Shape()
    return (falses(64), [], [], 1, nothing, nothing, nothing, nothing,0,0,0,0)
end
# ==============================================================================
# ==============================================================================
# ==============================================================================
# ==============================================================================
mutable struct NBox <: Draw
    flags::BitArray{1}
    color::Vector{Any}
    properties::Vector{Any}
    gradient::Vector{Any}
    opacity::Any
    padding::Union{Nothing,BoxOutline}
    border::Union{Nothing,Border}
    margin::Union{Nothing,BoxOutline}
    offset::Union{Nothing,Point}
    left::Float64
    top::Float64
    width::Float64
    height::Float64
    NBox() = new( falses(64), [], [], [], 1, nothing, nothing, nothing, nothing,0,0,0,0 )
end

mutable struct Circle <: Draw
    flags::BitArray{1}
    color::Vector{Any}
    gradient::Vector{Any}
    opacity::Any
    padding::Union{Nothing,BoxOutline}
    border::Union{Nothing,Border}
    margin::Union{Nothing,BoxOutline}
    offset::Union{Nothing,Point}
    left::Float64
    top::Float64
    width::Float64
    height::Float64
        radius::Float64
    Circle() = new(falses(64), [], [], 1, nothing, nothing, nothing, nothing,0,0,0,0,0,)
end
# Box, RoundBox, Arc, Circle, Line, Curve, TextElement, Ellipse
# Polygon, Polyline, Path
mutable struct Arc <: Draw
    flags::BitArray{1}
    color::Vector{Any}
    gradient::Vector{Any}
    opacity::Any
    padding::Union{Nothing,BoxOutline}
    border::Union{Nothing,Border}
    margin::Union{Nothing,BoxOutline}
    offset::Union{Nothing,Point}
    left::Float64
    top::Float64
    width::Float64
    height::Float64
        radius::Float64
        origin::Point
    startAngle::Float64
    stopAngle::Float64
    Arc() = new(0, Point(0,0), 0,0)
end

#=---------------------------------=#
mutable struct Font #<: Draw
    flags::BitArray{1}
    color::Vector{Any}
    gradient::Vector{Any}
    opacity::Any
    padding::Union{Nothing,BoxOutline}
    border::Union{Nothing,Border}
    margin::Union{Nothing,BoxOutline}
    offset::Union{Nothing,Point}
    size::Float64
    lineHeight::Float16
    family::String
    fill::Vector{Float32}
    lineWidth::Float64
    Font() = new( falses(64),
        [], [], 1, nothing, nothing, nothing, nothing,
        12, 1.4,  "Sans", [], 1 )
end
#=---------------------------------=#
# This is inserted as if a normal Element node and yet is also treated as a shape.
mutable struct TextLine <: Draw
    flags::BitArray{1}
    parent::Any # Because Element()s have parents
    text::String #
    left::Float64
    top::Float64
    width::Float64
    height::Float64
    TextLine(parent, text) = new(falses(64), parent, text, 0, 0, 0, 0)
    TextLine(parent, text, left, top, width, height) = new(falses(64), parent, text, left, top, width, height)
    TextLine(parent, text, left, top) = new(falses(64), parent, text, left, top,0,0)
end

#==============================================================================#

#==============================================================================#


text:
module VisualizationTests

using Test

@testset "Vtk" begin include("VtkTests.jl") end

end # module


text:
# demo for double loop
for i = 1:2, j = 3:4
    println((i, j))
end

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
using Plots

export ellipsoid_im

"""
`phantom = ellipsoid_im(ig, params;
    oversample=1, checkfov=false, how=:slow, showmem=false, hu_scale=1, return_params=false)`

generate ellipsoid phantom image from parameters:
    `[x_center, y_center, z_center, x_radius, y_radius, z_radius,
        xy_angle_degrees, z_angle_degrees, density]`

in
* `ig`            image_geom()
* `params`        `[9 N]` ellipsoid parameters.
                            if empty, use 3d shepp-logan

option
* `oversample::Integer`       oversampling factor (default:1)
* `checkfov::Bool`            warn if any ellipsoid is out of fov
* `how::Symbol`              `:fast` does it fast -- to do, only works slow
                                `:lowmem` uses less memory than :fast but a bit slower
                                `:slow` default
* `showmem::Bool`
* `hu_scale::Real`            use 1000 to scale shepp-logan to HU (default:1)
* `return_params::Bool`       if true, return both phantom and params

out
* `phantom`         `[nx ny nz]` image
* `params`          `[9 N]` ellipsoid parameters (only return if return_params=true)
"""
function ellipsoid_im(ig::MIRT_image_geom,
    params::AbstractArray{<:Real,2};
    oversample::Integer=1,
    checkfov::Bool=false,
    how::Symbol=:slow,
    showmem::Bool=false,
    hu_scale::Real=1,
    return_params::Bool=false)

    size(params,2) != 9 && throw("bad cuboid parameter vector size")
    params[:,9] .*= hu_scale
    args = (ig.nx, ig.ny, ig.nz, params, ig.dx, ig.dy, ig.dz,
            ig.offset_x, ig.offset_y, ig.offset_z, oversample, showmem)
    phantom = zeros(Float32, ig.nx, ig.ny, ig.nz)

    if checkfov
    	if !ellipsoid_im_check_fov(args...)
    		throw("ellipsoid exceeds FOV")
    	end
    end

    if how == :lowmem
        phantom += ellipsoid_im_lowmem(args...)
    elseif how == :fast
        phantom += ellipsoid_im_fast(args...)
    elseif how == :slow
        phantom += ellipsoid_im_slow(args...)
    else
        throw("bad how $how")
    end

    if return_params
        return (phantom, params)
    end
    return phantom
end

"""
`ellipsoid_im_slow()`

brute force fine grid - can use lots of memory
"""
function ellipsoid_im_slow(nx, ny, nz, params, dx, dy, dz,
    offset_x, offset_y, offset_z, over, showmem)

    phantom = zeros(Float32, nx*over, ny*over, nz*over)

    wx = (nx*over - 1)/2 + offset_x*over
    wy = (ny*over - 1)/2 + offset_y*over
    wz = (nz*over - 1)/2 + offset_z*over
    x = ((0:(nx * over - 1)) .- wx) * dx / over
    y = ((0:(ny * over - 1)) .- wy) * dy / over
    z = ((0:(nz * over - 1)) .- wz) * dz / over
    xmax = maximum(x)
    xmin = minimum(x)
    ymax = maximum(y)
    ymin = minimum(y)
    zmax = maximum(z)
    zmin = minimum(z)
    (xx, yy, zz) = ndgrid(x, y, z)

    #ticker reset
    np = size(params)[1]
    for ip in 1:np
        #ticker(mfilename, ip, np)

        par = params[ip, :]
        cx = par[1]
        cy = par[2]
        cz = par[3]
        rx = par[4]
        ry = par[5]
        rz = par[6]

        azim = par[7] * (pi/180)
        polar = par[8] * (pi/180)
        (xr, yr, zr) = rot3(xx .- cx, yy .- cy, zz .- cz, azim, polar)
        tmp = ((xr / rx).^2 + (yr / ry).^2 + (zr / rz).^2) .<=1
        value = Float32(par[9])
        phantom += value * tmp
    end
    return phantom
end

#=
"""
`ellipsoid_im_fast()`

currently not working
only slow option works
"""
function ellipsoid_im_fast(nx, ny, nz, params, dx, dy, dz,
    offset_x, offset_y, offset_z, over, showmem)

    phantom = zeros(Float32, nx, ny, nz)

    wx = (nx-1)/2 + offset_x
    wy = (ny-1)/2 + offset_y
    wz = (nz-1)/2 + offset_z
    x = ((0:nx-1) .- wx) * dx
    y = ((0:ny-1) .- wy) * dy
    z = ((0:nz-1) .- wz) * dz
    xmax = maximum(x)
    xmin = minimum(x)
    ymax = maximum(y)
    ymin = minimum(y)
    zmax = maximum(z)
    zmin = minimum(z)
    (xx, yy, zz) = ndgrid(x, y, z)

    if over > 1
        tmp = ((1:over) - (over+1)/2) / over
        (xf, yf, zf) = ndgrid(tmp*dx, tmp*dy, tmp*dz)
        xf = xf[:]'
        yf = yf[:]'
        zf = zf[:]'
    end

    hx = abs(dx) / 2
    hy = abs(dy) / 2
    hz = abs(dz) / 2

    #ticker reset
    np = size(params)[1]
    for ip in 1:np
        #ticker(mfilename, ip, np)
        par = params[ip,:]
        cx = par[1]
        cy = par[2]
        cz = par[3]
        rx = par[4]
        ry = par[5]
        rz = par[6]

        azim = par[7] * (pi/180)
        polar = par[8] * (pi/180)

        xs = xx .- cx
        ys = yy .- cy
        zs = zz .- cz

        (xr, yr, zr) = rot3(xs, ys, zs, azim, polar)
        if over == 1
            vi = ((xr / rx).^2 + (yr / ry).^2 + (zr / rz).^2) .<= 1
            value = Float32(par[9])
            phantom += value * vi
        end

        #= check all 8 corners of the cube that bounds
        all possible 3D rotations of the voxel =#
        hh = sqrt(hx^2 + hy^2 + hz^2)
        vi = true # voxels entirely inside the ellipsoid
        #vo = true # voxels entirely outside the ellipsoid

        for ix in -1:1, iy in -1:1, iz in -1:1
            xo = xr .+ ix*hh
            yo = yr .+ iy*hh
            zo = zr .+ iz*hh
            vi = vi .& (((xo/rx).^2 + (yo/ry).^2 + (zo/rz).^2) .< 1)
        end

        vo = false #to do. for now, "outside" test is failing
        if ip == 1
            throw("todo:must debug this")
        end
        if any(vi[:] & vo[:])
            throw("bug")
        end
        #=
        if 0
        	% coordinates of "outer" corner of each voxel, relative to ellipsoid center
        	sign_p = @(x) (x >= 0) * sqrt(3); % modified sign()
        	xo = xr + sign_p(xr) * hx;
        	yo = yr + sign_p(yr) * hy;
        	zo = zr + sign_p(zr) * hz;

        	% voxels that are entirely inside the ellipsoid
        	vi = (xo / rx).^2 + (yo / ry).^2 + (zo / rz).^2 <= 1;
        end

        if 0
        	% coordinates of "inner" corner of each pixel, relative to ellipse center
        	sign_n = @(x) (x > 0) * sqrt(3); % modified sign()
        	xi = xr - sign_n(xs) * hx;
        	yi = yr - sign_n(ys) * hy;
        	zi = zr - sign_n(zs) * hz;

        	% voxels that are entirely outside the ellipsoid
        	vo = (max(abs(xi),0) / rx).^2 + (max(abs(yi),0) / ry).^2 ...
        		+ (max(abs(zi),0) / rz).^2 >= 1;
        end
        =#

        edge = !vi & !vo
        if showmem
            println("edge fraction ",
                sum(edge[:]) / prod(size(edge)), "=",
                sum(edge[:]), "/", prod(size(edge)))
        end

        x = xx[edge] .- cx
        y = yy[edge] .- cy
        z = zz[edge] .- cz
        x = x .+ xf'
    	y = y .+ yf'
    	z = z .+ zf'
        (xr, yr, zr) = rot3(x, y, z, azim, polar)
        in = ((xr/rx).^2 + (yr/ry).^2 + (zr/rz).^2) .<= 1
        tmp = mean(in, 2)
        gray = Float32(vi)
        gray[edge] = tmp
        #=
        if 0 % todo: help debug
            ee = (gray > 0) & (gray < 1);
            im(ee)
            im(vi)
	    prompt
		    %minmax(tmp)
		    %clf, im pl 1 2
		    %im(tmp)
	    end=#

        value = Float32(par[9])
        phantom += value*gray
    end
    return phantom
end

"""
`ellipsoid_im_lowmem()`

this version does 'one slice at a time' to reduce memory
"""
function ellipsoid_im_lowmem(nx, ny, nz, params, dx, dy, dz,
    offset_x, offset_y, offset_z, over, showmem)

    phantom = zeros(Float32, nx, ny, nz)
    for iz in 1:nz
        offset_z_new = (nz-1)/2 + offset_z - (iz-1)
        phantom[:,:,iz] = ellipsoid_im_fast(nx, ny, 1, params, dx, dy, dz,
                            offset_x, offset_y, offset_z_new, over, showmem && iz == 1)
    end
    return phantom
end
=#

"""
`ellipsoid_im_check_fov()`
"""
function ellipsoid_im_check_fov(nx, ny, nz, params,
        dx, dy, dz, offset_x, offset_y, offset_z, over, showmem)
    wx = (nx - 1)/2 + offset_x
    wy = (ny - 1)/2 + offset_y
    wz = (nz - 1)/2 + offset_z
    xx = ((0:nx - 1) .- wx) * dx
    yy = ((0:ny - 1) .- wy) * dy
    zz = ((0:nz - 1) .- wz) * dz

    xmax = maximum(xx)
    xmin = minimum(xx)
    ymax = maximum(yy)
    ymin = minimum(yy)
    zmax = maximum(zz)
    zmin = minimum(zz)

    ok = true

    for ip in 1:size(params)[1]
        par = params[ip, :]
        cx = par[1]
        cy = par[2]
        cz = par[3]
        rx = par[4]
        ry = par[5]
        rz = par[6]

        if (cx + rx > xmax) || (cx - rx < xmin)
            throw("fov: x range $xmin $xmax, cx=$cx, rx=$rx")
            ok = false
        end
        if (cy + ry > ymax) || (cy - ry < ymin)
            throw("fov: y range $ymin $ymax, cy=$cy, ry=$ry")
            ok = false
        end

        if (cz + rz > zmax) || (cz - rz < zmin)
            throw("fov: z range $zmin $zmax, cz=$cz, rz=$rz")
		    ok = false
	    end
    end
    return ok
end

"""
`phantom = ellipsoid_im(nx, dx, params; args...)`

specifying voxel size of `dx` and ellipsoid `params`
"""
function ellipsoid_im(nx::Integer, dx::Real, params;args...)
    ig = image_geom(nx=nx, dx=1)
    return ellipsoid_im(ig, params; args...)
end

"""
`phantom = ellipsoid_im(nx::Integer, params; args...)`

voxel size of `dx=1` and ellipsoid `params`
"""
function ellipsoid_im(nx::Integer, params; args...)
    return ellipsoid_im(nx, 1., params; args...)
end

"""
`phantom = ellipsoid_im(nx::Integer, ny::Integer=nx, dx::Real=1, args...)`

image of size `nx` by `ny` (default `nx`) with specified `dx` (default 1)
defaults to `:zhu`
"""
function ellipsoid_im(nx::Integer, ny::Integer=nx, dx::Real=1, args...)
    if image_geom(nx=nx, ny=ny, dx=dx, args...)
        return ellipsoid_im(ig, :zhu; args...)
    end
end

"""
`phantom = ellipsoid_im(ig, ptype; args...)`

`ptype = :zhu | :kak | :e3d | :spheroid`
"""
function ellipsoid_im(ig::MIRT_image_geom, ptype::Symbol; args...)
    xfov = ig.fovs[1]
    yfov = ig.fovs[2]
    zfov = ig.fovs[3]

    if ptype == :zhu
        params = shepp_logan_3d_parameters(xfov, yfov, zfov, :zhu)
    elseif ptype == :kak
        params = shepp_logan_3d_parameters(xfov, yfov, zfov, :kak)
    elseif ptype == :e3d
        params = shepp_logan_3d_parameters(xfov, yfov, zfov, :e3d)
    elseif ptype == :spheroid
        params = spheroid_params(xfov, yfov, zfov, ig.dx, ig.dy, ig.dz)
    else
        throw("bad phantom symbol $ptype")
    end
    return ellipsoid_im(ig, params; args...)
end

"""
`phantom = ellipsoid_im(ig; args...)`

`:zhu` (default) for given image geometry `ig`
"""
function ellipsoid_im(ig::MIRT_image_geom; args...)
    return ellipsoid_im(ig, :zhu; args...)
end

"""
`rot3`
"""
function rot3(x, y, z, azim, polar)
    if polar != 0
        throw("z (polar) rotation not done")
    end
    xr = cos(azim) * x + sin(azim) * y
    yr = -sin(azim) * x + cos(azim) * y
    zr = z
    return (xr, yr, zr)
end

"""
`shepp_logan_3d_parameters()`

most of these values are unitless 'fractions of field of view'
"""
function shepp_logan_3d_parameters(xfov, yfov, zfov, ptype)
    # parameters from Kak and Slaney, typos ?
    ekak = Float32[
        0       0       0       0.69    0.92    0.9     0   0   2.0
        0       0       0       0.6624  0.874   0.88    0   0   -0.98
        -0.22   0       -0.25   0.41    0.16    0.21    0   0   -0.98
        0.22    0       -0.25   0.31    0.11    0.22    72  0  -0.02
        0       0.1     -0.25   0.046   0.046   0.046   0   0   0.02
        0       0.1     -0.25   0.046   0.046   0.046   0   0   0.02
        -0.8    -0.65   -0.25   0.046   0.023   0.02    0   0   0.01
        0.06    -0.065  -0.25   0.046   0.023   0.02    90  0  0.01
        0.06    -0.105  0.625   0.56    0.04    0.1     90  0  0.02
        0       0.1     -0.625  0.056   0.056   0.1     0   0   -0.02
    ]

    #parameters from leizhu@standford.edu, says kak&slaney are incorrect
    ezhu = Float32[
    0      0         0       0.69        0.92        0.9     0      0       2.0
	0      -0.0184   0       0.6624	     0.874       0.88	 0      0	    -0.98
	-0.22  0	     -0.25   0.41	     0.16	     0.21	 -72    0	    -0.02
	0.22   0	     -0.25	 0.31	     0.11	     0.22	 72     0	     -0.02
	0	   0.35	     -0.25	 0.21	     0.25	     0.35	 0      0	     0.01
	0	   0.1	     -0.25	 0.046	     0.046	     0.046	 0      0	     0.01
	-0.08  -0.605	 -0.25	 0.046	     0.023	     0.02	 0      0	     0.01
	0	   -0.1	     -0.25	 0.046	     0.046	     0.046	 0      0	     0.01
	0	   -0.605	 -0.25   0.023	     0.023     	 0.023	 0      0	     0.01
	0.06   -0.605	 -0.25	 0.046	     0.023	     0.02	 -90    0     	 0.01
	0.06   -0.105	 0.0625	 0.056	     0.04	     0.1	 -90    0      	 0.02
	0	   0.1	     0.625	 0.056	     0.056	     0.1	 0      0	     -0.02
    ]

    if ptype == :zhu
        params = ezhu
    elseif ptype == :kak
        params = ekak
    #elseif ptype == :shepp_logan_e3d || ptype == :e3d
        #params = e3d[:, [5:7 2:4 8 1]]
    else
        throw("unknown parameter type $ptype")
    end

    params[:,[1,4]] .*= xfov/2
    params[:,[2,5]] .*= yfov/2
    params[:,[3,6]] .*= zfov/2
    return params
end

"""
`spheroid_params()`
"""
function spheroid_params(xfov, yfov, zfov, dx, dy, dz)
    #xfov = nx * dx, number * size
    xradius = (xfov/2) - dx
    yradius = (yfov/2) - dy
    zradius = (zfov/2) - dz
    params =[
    0 0 0 xradius yradius zradius 0 0 1
    ]
    return params
end

"""
`ellipsoid_im()`

show docstrings
"""
function ellipsoid_im()
    @doc ellipsoid_im
end

"""
`ellipsoid_im_show()`
"""
function ellipsoid_im_show()
    ig = image_geom(nx=512, nz=64, dz=0.625, fov=500)
    ig = ig.down(8)

    x = ellipsoid_im(ig, :zhu; hu_scale=1000, how=:slow)
    p3 = jim(x, title="zhu") #clim=[900,1100]

    spheroid = ellipsoid_im(ig, :spheroid; how=:slow)
    p2 = jim(spheroid, title="spheroid")

    #x2 = ellipsoid_im(ig, :zhu; how=:lowmem)
    #p2 = jim(x2, title="zhu, lowmem")

    plot(p2, p3)
end

"""
`ellipsoid_im_test()`
"""
function ellipsoid_im_test()
    ig = image_geom(nx=512, nz=64*2, dz=0.625, fov=500)
    ig = ig.down(8)

    #test different ellipses
    e1 = ellipsoid_im(ig)
    e2 = ellipsoid_im(ig, :zhu)
    e3 = ellipsoid_im(ig, :kak)
    e4 = ellipsoid_im(ig; checkfov=true)
    e5 = ellipsoid_im(ig; showmem=true)
    e6 = ellipsoid_im(ig; return_params=true)

    #ell1 = ellipsoid_im(ig, :zhu; how=:fast) # fast doesn't work
    #ell2 = ellipsoid_im(ig, :zhu; how=:lowmem) # lowmem calls fast - doesn't work
end

"""
`ellipsoid_im(:test)`

`ellipsoid_im(:show)`

run tests
"""
function ellipsoid_im(test::Symbol)
	if test == :show
		return ellipsoid_im_show()
	end
	test != :test && throw(ArgumentError("test $test"))
	ellipsoid_im()
    ellipsoid_im(:show)
	ellipsoid_im_test()
	true
end

#ellipsoid_im(:show)
#ellipsoid_im(:test)


text:
"""
"""
struct RestrictedDiscreteModel{Dc,Dp} <: DiscreteModel{Dc,Dc}
  model::DiscreteModelPortion{Dc,Dp}
end

function RestrictedDiscreteModel(model::DiscreteModel, cell_to_oldcell::Vector{Int})
  _model = DiscreteModelPortion(model,cell_to_oldcell)
  RestrictedDiscreteModel(_model)
end

function RestrictedDiscreteModel(model::DiscreteModel, cell_to_mask::Vector{Bool})
  _model = DiscreteModelPortion(model,cell_to_mask)
  RestrictedDiscreteModel(_model)
end

get_grid(model::RestrictedDiscreteModel) = get_grid(model.model)

get_grid_topology(model::RestrictedDiscreteModel) = get_grid_topology(model.model)

get_face_labeling(model::RestrictedDiscreteModel) = get_face_labeling(model.model)

get_face_to_oldface(model::RestrictedDiscreteModel,d::Integer) = get_face_to_oldface(model.model,d)

get_cell_to_oldcell(model::RestrictedDiscreteModel) = get_cell_to_oldcell(model.model)

get_oldmodel(model::RestrictedDiscreteModel) = get_oldmodel(model.model)

function RestrictedTriangulation(model::RestrictedDiscreteModel)
  oldtrian = Triangulation(get_oldmodel(model))
  cell_to_oldcell = get_cell_to_oldcell(model)
  RestrictedTriangulation(oldtrian,cell_to_oldcell)
end

function Triangulation(model::RestrictedDiscreteModel)
  RestrictedTriangulation(model)
end

function get_triangulation(model::RestrictedDiscreteModel)
  RestrictedTriangulation(model)
end

function Triangulation(::Type{ReferenceFE{d}},model::RestrictedDiscreteModel) where d
  @notimplemented
end

function Triangulation(model::RestrictedDiscreteModel,cell_to_oldcell::AbstractVector{<:Integer})
  @notimplemented
end

function Triangulation(model::RestrictedDiscreteModel,cell_to_mask::AbstractVector{Bool})
  @notimplemented
end

function BoundaryTriangulation(model::RestrictedDiscreteModel,face_to_mask::Vector{Bool},icell_around::Integer)
  d = num_cell_dims(model)-1
  face_to_oldface = get_face_to_oldface(model,d)
  oldmodel = get_oldmodel(model)
  num_oldfaces = num_faces(oldmodel,d)
  oldface_to_mask = fill(false,num_oldfaces)
  oldface_to_mask[face_to_oldface] .= face_to_mask
  BoundaryTriangulation(oldmodel,oldface_to_mask)
end

function InterfaceTriangulation(model_in::RestrictedDiscreteModel,model_out::RestrictedDiscreteModel)
  cells_in = get_cell_to_oldcell(model_in)
  cells_out = get_cell_to_oldcell(model_out)
  model = get_oldmodel(model_in)
  InterfaceTriangulation(model,cells_in,cells_out)
end


text:
original_decks = open("cards.txt") do f
    map(m -> parse.(Int, split(strip(m.captures[1]), "\n")), eachmatch(r"Player .:\n((?:\d+\n?)+)", read(f, String)))
end

function play_combat(decks::Array{Array{Int64,1},1})::Int
    while minimum(map(d -> length(d), decks)) > 0
        topcards = map(d -> popfirst!(d), decks)
        winner = argmax(topcards)
        decks[winner] = append!(decks[winner], winner == 1 ? topcards : reverse(topcards))
    end

    return argmax(map(d -> length(d), decks))
end

function return_first(x::Array{Int64,1}, y::Int64)::Array{Int64,1}
    return deepcopy(x[1:y])
end

function play_recursive_combat(decks::Array{Array{Int64,1},1})::Int
    # println("Starting a game")

    history = Array{Array{Array{Int64,1},1},1}()
    while minimum(map(d -> length(d), decks)) > 0
        # for i in 1:2
        #     println("Player $(i)'s deck: $(join(decks[i], ", "))")
        # end
        if decks in history
            # println("We've already seen these decks, player 1 wins the game")
            return 1
        end

        push!(history, deepcopy(decks))
        topcards = map(d -> popfirst!(d), decks)
        # for i in 1:2
        #     println("Player $(i) plays: $(topcards[i])")
        # end

        if all(map(d -> length(d), decks) .>= topcards)
            new_decks = map(return_first, decks, topcards)
            winner = play_recursive_combat(new_decks)
        else
            winner = argmax(topcards)
            # println("Player $(winner) wins the round")
        end

        decks[winner] = append!(decks[winner], winner == 1 ? topcards : reverse(topcards))
    end

    # println("Player $(argmax(map(d -> length(d), decks))) wins the game")
    return argmax(map(d -> length(d), decks))

end

decks = deepcopy(original_decks)
winner = play_combat(decks)
println(sum(decks[winner] .* collect(length(decks[winner]):-1:1)))

decks = deepcopy(original_decks)
winner = play_recursive_combat(decks)
println(sum(decks[winner] .* collect(length(decks[winner]):-1:1)))


text:
type SepShot3C
    ux
    uy
    uz
    up
    us
    vp
    vs
    fmin
    fmax
    decomp
    H
    verbose
end

function Sep1Shot(shot)

    run( `wavesep H=$(shot.H)
        ux=$(shot.ux) uy=$(shot.uy) uz=$(shot.uz)
	up=$(shot.up) us=$(shot.us)
	vp=$(shot.vp) vs=$(shot.vs)
	fmin=$(shot.fmin) fmax=$(shot.fmax)
	verbose=$(shot.verbose)`)
    return(1)

end

function WaveSep(ux, uy, uz, up, us;
                 H=true, decomp=true, vp=2000, vs=1000, fmin=0, fmax=80,
                 verbose=true, isx=0, isy=0)

    # H :  flag for exact Helmholtz decomp op (true), or adjoint of inverse
    # Helmholtz decomp op (false)
    # decomp : flag for decomposition from data components to potentials (true),
    # or recomposition from potentials to data components (false)
    # vp :  the P-wave velocity
    # vs : the S-wave velocity
    # fmin : min frequency to process (Hz)
    # fmax : max frequency to process (Hz)
    # verbose : flag for error / debugging messages
    # isx : array of source bin numbers in the x direction
    # isy : array of source bin numbers in the y direction

    nshot = length(isx)

    shot_list = Array(SepShot3C,nshot)
    for ishot = 1 : nshot
	shot_list[ishot] = SepShot3C(0,0,0,0,0,0,0,0,0,0,0,0)
	shot_list[ishot].ux = join([ux "_shot_" int(floor(isx[ishot])) "_" int(floor(isy[ishot]))])
	shot_list[ishot].uy = join([uy "_shot_" int(floor(isx[ishot])) "_" int(floor(isy[ishot]))])
	shot_list[ishot].uz = join([uz "_shot_" int(floor(isx[ishot])) "_" int(floor(isy[ishot]))])
	shot_list[ishot].up = join([up "_shot_" int(floor(isx[ishot])) "_" int(floor(isy[ishot]))])
	shot_list[ishot].us = join([us "_shot_" int(floor(isx[ishot])) "_" int(floor(isy[ishot]))])
	shot_list[ishot].vp = vp
	shot_list[ishot].vs = vs
	shot_list[ishot].fmin = fmin
	shot_list[ishot].fmax = fmax
	shot_list[ishot].decomp = decomp == true ? "y" : "n"
	shot_list[ishot].H = H == true ? "y" : "n"
	shot_list[ishot].verbose = verbose == true ? "y" : "n"
    end
    if (decomp == "y")
	for ishot = 1 : nshot
	    SeisWindow(ux, shot_list[ishot].ux, key=["isx","isy"],
                       minval=[isx[ishot],isy[ishot]],
                       maxval=[isx[ishot],isy[ishot]])
	    SeisWindow(uy, shot_list[ishot].uy, key=["isx","isy"],
                       minval=[isx[ishot],isy[ishot]],
                       maxval=[isx[ishot],isy[ishot]])
	    SeisWindow(uz, shot_list[ishot].uz, key=["isx","isy"],
                       minval=[isx[ishot],isy[ishot]],
                       maxval=[isx[ishot],isy[ishot]])
	end
	a = pmap(Sep1Shot,shot_list)
	j = 1
	for ishot = 1 : nshot
	    up_shot,h_shot,e = SeisRead(shot_list[ishot].up)
	    for itrace = 1 : size(up_shot,2)
		h_shot[ishot].tracenum = j + itrace
	    end
	    SeisWrite(up,up_shot,h_shot,itrace=j)
	    us_shot,h_shot,e = SeisRead(shot_list[ishot].us)
	    for itrace = 1 : size(us_shot,2)
		h_shot[ishot].tracenum = j + itrace
	    end
	    SeisWrite(us,us_shot,h_shot,itrace=j)
	    j += size(up_shot,2)
	    SeisRemove(shot_list[ishot].up)
	    SeisRemove(shot_list[ishot].us)
	    SeisRemove(shot_list[ishot].ux)
	    SeisRemove(shot_list[ishot].uy)
	    SeisRemove(shot_list[ishot].uz)

	end
    else
	for ishot = 1 : nshot
	    SeisWindow(up, shot_list[ishot].up, key=["isx","isy"],
                       minval=[isx[ishot],isy[ishot]],
                       maxval=[isx[ishot],isy[ishot]])
	    SeisWindow(us, shot_list[ishot].us, key=["isx","isy"],
                       minval=[isx[ishot],isy[ishot]],
                       maxval=[isx[ishot],isy[ishot]])
	end
	a = pmap(Sep1Shot,shot_list)
	j = 1
	for ishot = 1 : nshot
	    ux_shot,h_shot,e = SeisRead(shot_list[ishot].ux)
	    SeisWrite(ux,ux_shot,h_shot,itrace=j)
	    uy_shot,h_shot,e = SeisRead(shot_list[ishot].uy)
	    SeisWrite(uy,uy_shot,h_shot,itrace=j)
	    uz_shot,h_shot,e = SeisRead(shot_list[ishot].uz)
	    SeisWrite(uz,uz_shot,h_shot,itrace=j)
	    j += size(ux_shot,2)
	    SeisRemove(shot_list[ishot].up)
	    SeisRemove(shot_list[ishot].us)
	    SeisRemove(shot_list[ishot].ux)
	    SeisRemove(shot_list[ishot].uy)
	    SeisRemove(shot_list[ishot].uz)
	end
    end
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using Weave

path = "C:\\Users\\Jasmine\\.julia\\dev\\DCDC\\doc"
weave(joinpath(path,"BandwidthSelection.jmd"),doctype="pandoc",out_path=path)


text:
using PaddedMatrices: AbstractMutableFixedSizeMatrix
###
### We are solving for A in
### B = A * L
### or
### B = A * U
### where L is lower triangular and U is upper triangular.
### That is, we're solving the system of triangular equations
### A = B / L
### or
### A = B / U
### Our strategy for doing so involves solving for blocks
### from A at a time, iterating through earlier/later
### columns to update when solving  U/K
### A[:,k] = ( B[:,k] - ∑ⱼ₌ₖ₊₁ᴷ A[:,j] * L[j,k] ) / L[k,k]
### or
### A[:,k] = ( B[:,k] - ∑ⱼ₌₁ᵏ⁻¹ A[:,j] * U[j,k] ) / U[k,k]
### eg a 2vec x 3 block from each will allow
###
Base.@kwdef mutable struct A_rdiv_U_config{T}
    R::Union{Symbol,Int} = -1
    C::Int = -1
    K::Union{Symbol,Int} = -1
    Astride::Union{Symbol,Int} = -1
    Bstride::Union{Symbol,Int} = -1
    Ustride::Int = -1
    isL′::Bool = false
    invdiagptr::Bool = false
    Bsym::Symbol = :ptrB
    Asym::Symbol = :ptrA
    Utrisym::Symbol = :ptrUtri
    Udiagsym::Symbol = :ptrUdiag
    maskload::Bool = true
    loadB::Bool = true
    storeA::Bool = true
    reduce_sym::Union{Symbol,Nothing} = nothing
    maskexpr::Symbol = :__mask__
end

# Ustride
# if ilL′

# K is number of earlier iterations
# divides B by U, stores in A
function A_rdiv_U_kernel_quote(
    config::A_rdiv_U_config{T}
) where {T}
    @unpack R, C, K, Astride, Bstride, Ustride, isL′, invdiagptr, Bsym, Asym, Utrisym, Udiagsym, maskload, loadB, storeA, reduce_sym, maskexpr = config
    # K is either a symbol or integer, indicating number of preceding columns
    size_T = sizeof(T)
    if R isa Integer
        W, Wshift = VectorizationBase.pick_vector_width_shift(R, T)
        Wm1 = W - 1
        Riter = R >>> Wshift
        Rrem = R & Wm1
        mask = VectorizationBase.mask(T, Rrem)
    else # We assume this is meant to handle a single vector remainder
        W, Wshift = VectorizationBase.pick_vector_width_shift(T)
        Wm1 = W - 1
        Riter = 0
        Rrem = 1
        mask = maskexpr# :(VectorizationBase.mask($T, $R & $Wm1))
    end
    q = quote end
    if loadB
        if K isa Symbol && Bstride isa Symbol
            push!(q.args, :( BsymK = $Bsym + $size_T * $Bstride * $K ))
        elseif K isa Symbol
            push!(q.args, :( BsymK = $Bsym + $(size_T*Bstride)*$K ))
        elseif Bstride isa Symbol
            push!(q.args, :( BsymK = $Bsym + $(size_T*K)*$Bstride ))
        else
            push!(q.args, :( BsymK = $Bsym + $(size_T*Bstride*K) ))
        end
        for c ∈ 0:C-1
            if Bstride isa Symbol
                for r ∈ 0:Riter-1
                    push!(q.args, :($(Symbol(:A_,r,:_,c)) = vload(Vec{$W,$T}, BsymK + $(size_T * (r*W)) + $(c*size_T)*$Bstride)) )
                end
                if Rrem > 0
                    # Only need to mask if we're on last column
                    if maskload && c == C-1
                        push!(q.args, :($(Symbol(:A_,Riter,:_,c)) = vload(Vec{$W,$T}, BsymK + $(size_T * (Riter*W)) + $(c*size_T)*$Bstride, $mask )))
                    else
                        push!(q.args, :($(Symbol(:A_,Riter,:_,c)) = vload(Vec{$W,$T}, BsymK + $(size_T * (Riter*W)) + $(c*size_T)*$Bstride )))
                    end
                end
            else
                for r ∈ 0:Riter-1
                    push!(q.args, :($(Symbol(:A_,r,:_,c)) = vload(Vec{$W,$T}, BsymK + $(size_T * (r*W + c*Bstride)) )))
                end
                if Rrem > 0
                    # Only need to mask if we're on last column
                    if maskload && c == C-1
                        push!(q.args, :($(Symbol(:A_,Riter,:_,c)) = vload(Vec{$W,$T}, BsymK + $(size_T * (Riter*W + c*Bstride)), $mask ) ))
                    else
                        push!(q.args, :($(Symbol(:A_,Riter,:_,c)) = vload(Vec{$W,$T}, BsymK + $(size_T * (Riter*W + c*Bstride)) ) ))
                    end
                end
            end
        end
    end
    if !isL′
        if Ustride isa Integer
            coff = Ustride
            for c ∈ 0:C-1
                push!(q.args, Expr(:(=), Symbol(:Uoffset_,c), coff*size_T))
                coff += c
            end
        else
            push!(q.args, Expr(:(=), :Uoffset_0, 0))
            for c ∈ 0:C-2
                push!(q.args, Expr(:(=), Symbol(:Uoffset_,c+1), :($(Symbol(:Uoffset_,c)) + $size_T*($c+K))))
            end
        end
    end
    # We don't need to mask these loads.
    Riterl = Rrem > 0 ? Riter : Riter-1
    # Updating based on all earlier columns
    if K isa Symbol || K > 0
        loopbody = quote end
        for r ∈ 0:Riterl
            push!(loopbody.args, :($(Symbol(:A_,r,:_j)) = vload(Vec{$W,$T}, $Asym + (j*$Astride + $(r*W))*$size_T)))
        end
        for c ∈ 0:C-1
            if isL′ # U is actually a transposed lower triangular matrix
                # c increase corresponds row increase
                push!(loopbody.args, :($(Symbol(:U_j_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Utrisym + $((c-1)*size_T)))))
            else
                push!(loopbody.args, :($(Symbol(:U_j_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Utrisym + $(Symbol(:Uoffset_,c)) + (j*$size_T)))))
            end
            for r ∈ 0:Riterl
                push!(loopbody.args, :($(Symbol(:A_,r,:_,c)) = vfnmadd($(Symbol(:A_,r,:_j)), $(Symbol(:U_j_,c)), $(Symbol(:A_,r,:_,c)))))
            end
        end
        if isL′
            if Ustride isa Symbol
                push!(loopbody.args, :($Utrisym += $size_T * ($Ustride - j - 2)))
            else
                push!(loopbody.args, :($Utrisym += $size_T * ($(Ustride-2) - j)))
            end
        end
        mainloop = quote
            @inbounds for j ∈ 0:$(K isa Symbol ? Expr(:call, :(-), K, 1) : K-1)
                $loopbody
            end
        end
        push!(q.args, mainloop)
    end
    # Update using just-calculated block
    scaleop = invdiagptr ? :vmul : :vfdiv
    if isL′
        trioff = 0#size_T
        #push!(q.args, :(ptrUtr -= $size_T))
    end
    for j ∈ 0:max(C-2,0)
        if j == 0
            vUjj = Symbol(:vU_,j,:_,j)
            push!(q.args, Expr(:(=), vUjj, :(vbroadcast(Vec{$W,$T}, VectorizationBase.load($Udiagsym + $(size_T*j))))))
            for r ∈ 0:Riterl
                Arj = Symbol(:A_,r,:_0)
                push!(q.args, Expr(:(=), Arj, :($scaleop($Arj, $vUjj))))
            end
            # if storeA
            #     for r ∈ 0:Riter-1
            #         push!(q.args, :( vstore!( AsymK + $(size_T * (r*W)), $(Symbol(:A_,r,:_0)) ) ))
            #     end
            #     if Rrem > 0
            #         push!(q.args, :( vstore!( AsymK + $(size_T * (Riter*W)), $(Symbol(:A_,Riter,:_0)), $mask ) ))
            #     end
            # end
            if reduce_sym isa Symbol
                for r ∈ 0:Riter-1
                    Arj = Symbol(:A_,r,:_0)
                    reduce_r = Symbol(reduce_sym, :_, r % 4)
                    push!(q.args, Expr(:(=), reduce_r, :(vmuladd($Arj, $Arj, $reduce_r))))
                end
                if Rrem > 0
                    Arj = Symbol(:A_,Riter,:_0)
                    reduce_r = Symbol(reduce_sym, :_, Riter % 4)
                    push!(q.args, Expr(:(=), reduce_r, :(vifelse($mask, vadd(vmul($Arj,$Arj),$reduce_r),$reduce_r))))
                end
            end
        end
        if isL′ && j > 0
            if K isa Symbol
                if Ustride isa Symbol
                    push!(q.args, :($Utrisym += $size_T*($Ustride - $K - $j)))
                else
                    push!(q.args, :($Utrisym += $size_T*($(Ustride-j) - $K)))
                end
            else
                trioff += size_T*(Ustride - K - j)
            end
        end
        for c ∈ j+1:C-1
            if isL′ # U is actually a transposed lower triangular matrix
                # c increase corresponds row increase
                #push!(q.args, :($(Symbol(:U_,j,:_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Utrisym + $((c-1-j)*size_T + trioff)))))
                push!(q.args, :($(Symbol(:U_,j,:_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Utrisym + $((c-1-j)*size_T + trioff)))))
            else
                push!(q.args, :($(Symbol(:U_,j,:_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Utrisym + $(Symbol(:Uoffset_,c)) + $(K isa Symbol ? :($size_T*($K+$j)) : size_T*(K+j) )))))
            end
            for r ∈ 0:Riterl
                push!(q.args, :($(Symbol(:A_,r,:_,c)) = vfnmadd($(Symbol(:A_,r,:_,j)), $(Symbol(:U_,j,:_,c)), $(Symbol(:A_,r,:_,c)))))
            end
            if c == j+1
                vUjj = Symbol(:vU_,c,:_,c)
                push!(q.args, Expr(:(=), vUjj, :(vbroadcast(Vec{$W,$T}, VectorizationBase.load($Udiagsym + $(size_T*c))))))
                for r ∈ 0:Riterl
                    push!(q.args, :($(Symbol(:A_,r,:_,c)) = $scaleop($(Symbol(:A_,r,:_,c)), $vUjj)))
                end
                # if storeA
                #     for r ∈ 0:Riter-1
                #         push!(q.args, :( vstore!( AsymK + $(size_T * (r*W + c*Astride)), $(Symbol(:A_,r,:_,c)) ) ))
                #     end
                #     if Rrem > 0
                #         push!(q.args, :( vstore!( AsymK + $(size_T * (Riter*W + c*Astride)), $(Symbol(:A_,Riter,:_,c)), $mask ) ))
                #     end
                # end
                if reduce_sym isa Symbol
                    for r ∈ 0:Riter-1
                        Arj = Symbol(:A_,r,:_,c)
                        reduce_r = Symbol(reduce_sym, :_, r % 4)
                        push!(q.args, Expr(:(=), reduce_r, :(vmuladd($Arj, $Arj, $reduce_r))))
                    end
                    if Rrem > 0
                        Arj = Symbol(:A_,Riter,:_,c)
                        reduce_r = Symbol(reduce_sym, :_, Riter % 4)
                        push!(q.args, Expr(:(=), reduce_r, :(vifelse($mask, vmuladd($Arj,$Arj,$reduce_r), $reduce_r))))
                    end
                end
            end
        end
        #if j == 0 && isL′
        #    trioff -= size_T
        #end
    end
    ### Store in A
    ### Seems faster when I place all the stores together at the end???
    if storeA
        # if K isa Symbol
        push!(q.args, :(AsymK = $Asym + $size_T*$Astride*$K))
        # else
            # push!(q.args, :(AsymK = $Asym + $(size_T*Astride*K)))
        # end
    # end
    # if storeA
        for c ∈ 0:C-1
            for r ∈ 0:Riter-1
                push!(q.args, :( vstore!( AsymK + $size_T * ($(r*W) + $c*$Astride), $(Symbol(:A_,r,:_,c)) ) ))
            end
            if Rrem > 0
                push!(q.args, :( vstore!( AsymK + $size_T * ($(Riter*W) + $c*$Astride), $(Symbol(:A_,Riter,:_,c)), $mask ) ))
            end
        end
    end
    q
end

Base.@kwdef mutable struct A_rdiv_L_config{T}
    R::Union{Symbol,Int} = -1
    C::Int = -1
    Kmax::Union{Symbol,Int} = -1
    Astride::Union{Symbol,Int} = -1
    Bstride::Union{Symbol,Int} = -1
    isU′::Bool = false
    invdiagptr::Bool = false
    Bsym::Symbol = :ptrB
    Asym::Symbol = :ptrA
    Ltrisym::Symbol = :ptrLtri
    Ldiagsym::Symbol = :ptrLdiag
    maskload::Bool = true
    loadB::Bool = true
    storeA::Bool = true
    calc_product::Int = 0
    maskexpr::Symbol = :__mask__
    increment_store::Bool = false
    v∂Lstride::Int = VectorizationBase.pick_vector_width(T)
end

# ptr and Lstride is so that the first of the reverse iterations has offset 0
# K is the amount of preceding columns.
function A_rdiv_L_kernel_quote(
    config::A_rdiv_L_config{T}
) where {T}
    @unpack R, C, Kmax, Astride, Bstride, isU′, invdiagptr, Bsym, Asym, Ltrisym, Ldiagsym, maskload, loadB, storeA, calc_product, maskexpr, increment_store, v∂Lstride = config
    # K is either a symbol or integer, indicating number of preceding columns
    size_T = sizeof(T)
    if isU′
        K = Kmax
    end
    if R isa Integer
        W, Wshift = VectorizationBase.pick_vector_width_shift(R, T)
        Wm1 = W - 1
        Riter = R >>> Wshift
        Rrem = R & Wm1
        mask = VectorizationBase.mask(T, Rrem)
    else # We assume this is meant to handle a single vector remainder
        W, Wshift = VectorizationBase.pick_vector_width_shift(T)
        Wm1 = W - 1
        Riter = 0
        Rrem = 1
        mask = maskexpr #:(VectorizationBase.mask($T, $R & $Wm1))
    end
    V = Vec{W,T}
    q = quote end
    if loadB
        for c ∈ 0:C-1
            for r ∈ 0:Riter-1
                set_asym = Symbol(:A_,r,:_,c)
                push!(q.args, :($set_asym = vload(Vec{$W,$T}, $Bsym + $size_T * ($(r*W) + $c*$Bstride)) ) )
            end
            if Rrem > 0
                # Only need to mask if we're on last column
                set_asym = Symbol(:A_,Riter,:_,c)
                if maskload && c == C-1
                    push!(q.args, :($set_asym = vload(Vec{$W,$T}, $Bsym + $size_T * ($(Riter*W) + $c*$Bstride), $mask) ))
                else
                    push!(q.args, :($set_asym = vload(Vec{$W,$T}, $Bsym + $size_T * ($(Riter*W) + $c*$Bstride)) ))
                end
            end
        end
    end
    # We don't need to mask these loads.
    Riterl = Rrem > 0 ? Riter : Riter-1
    if !isU′
        Ltrisym2 = Ltrisym
        coff = 0
        for c ∈ 0:C-1
            if Kmax isa Symbol
                #push!(q.args, Expr(:(=), Symbol(:Loffset_,c), :($(2c*size_T)*$Kmax + $((coff-2c*C)*size_T))))
                #coff -=  c + 2
                if c == 0
                    push!(q.args, Expr(:(=), Symbol(:Loffset_,c), 0  ))
                else
                    # push!(q.args, Expr(:(=), Symbol(:Loffset_,c), :($Kmax*$(c*size_T) - $(coff*size_T))  ))
                    push!(q.args, Expr(:(=), Symbol(:Loffset_,c), :( ($Kmax*$(c*size_T) +  $((coff-C*c)*size_T)))))
                end
            else
                push!(q.args, Expr(:(=), Symbol(:Loffset_,c), size_T * ((Kmax-C)*c + coff)))
            end
            coff += C - c - 2
        end
    end
    # Updating based on all earlier columns
    if Kmax isa Symbol || (C < Kmax)
        if isU′
            Ltrisym2 = Symbol(:Ltrisym, 2)
            CtriL = (C*(C+1)) >>> 1
            if K isa Symbol
#                push!(q.args, :(KpC = $K+$C))
                #                push!(q.args, Expr(:(=), Ltrisym2, :($Ltrisym + ((KpC*(KpC+1))>>>1) - (($K*($K+1))>>>1))))
                push!(q.args, Expr(:(=), Ltrisym2, :($CtriL + $K*$C)))
            else
#                KpC = K + C
                #                push!(q.args, Expr(:(=), Ltrisym2, :($Ltrisym + $(((KpC*(KpC+1))>>>1)-((K*(K+1))>>>1)))))
                push!(q.args, Expr(:(=), Ltrisym2, :($(CtriL + K*C))))
            end
        else
            push!(q.args, Expr(:(=), :ptrAloop, :($Asym + $(C*size_T)*$Astride)))
            push!(q.args, Expr(:(=), :Ltripointer2, :($Ltrisym2 + $((C-1)*size_T))))
        end
        loopbody = quote end
        for r ∈ 0:Riterl
            push!(loopbody.args, :($(Symbol(:A_,r,:_j)) = vload(Vec{$W,$T}, ptrAloop + j*$Astride*$size_T + $(r*W*size_T))))
        end
        if isU′
            push!(loopbody.args, :($Ltrisym += j))
        end
        for c ∈ 0:C-1
            if isU′ # U is actually a transposed lower triangular matrix
                # c increase corresponds row increase
                push!(loopbody.args, :($(Symbol(:L_j_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Ltrisym2 + $(c*size_T)))))
            else
                push!(loopbody.args, :($(Symbol(:L_j_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load(Ltripointer2 + $(Symbol(:Loffset_,c)) + j*$size_T))))
            end
            for r ∈ 0:Riterl
                push!(loopbody.args, :($(Symbol(:A_,r,:_,c)) = vfnmadd($(Symbol(:A_,r,:_j)), $(Symbol(:L_j_,c)), $(Symbol(:A_,r,:_,c)))))
            end
        end
        if isU′
            jloopstart = K isa Symbol ? :($K - $C) : K - C
            jloopend = Kmax isa Symbol ? :($Kmax - 1) : Kmax - 1
            # @show jloopend, isU′
        else
            jloopstart = 0#C-1
            jloopend = Kmax isa Symbol ? :($Kmax - $(C+1)) : Kmax - C - 1
            # @show jloopend, isU′
        end
        mainloop = quote
            @inbounds for j ∈ $jloopstart:$jloopend
                $loopbody
            end
        end
        push!(q.args, mainloop)
    end
    # Update using just-calculated block
    finish_block = quote end
    scaleop = invdiagptr ? :vmul : :vfdiv
    if isU′
        Uoffset = ((K*(K-1))>>>1) + (C*K) + (((C-2)*(C-1))>>>1)
    end
    for nj ∈ 1:C
        j = C - nj
        vLjj = Symbol(:vL_,j,:_,j)
        push!(q.args, Expr(:(=), vLjj, :(vbroadcast(Vec{$W,$T}, VectorizationBase.load($Ldiagsym + $(size_T*j))))))
        for r ∈ 0:Riterl
            push!(q.args, :($(Symbol(:A_,r,:_,j)) = $scaleop($(Symbol(:A_,r,:_,j)), $vLjj)))
        end
        for c ∈ 0:j-1
            if isU′ # U is actually a transposed lower triangular matrix
                # c increase corresponds row increase
                push!(q.args, :($(Symbol(:L_,j,:_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Ltrisym + $((c+Uoffset)*size_T)))))
            else
                push!(q.args, :($(Symbol(:L_,j,:_,c)) = vbroadcast(Vec{$W,$T}, VectorizationBase.load($Ltrisym + $(Symbol(:Loffset_,c)) + $((j-1)*size_T)))))
            end
            for r ∈ 0:Riterl
                push!(q.args, :($(Symbol(:A_,r,:_,c)) = vfnmadd($(Symbol(:A_,r,:_,j)), $(Symbol(:L_,j,:_,c)), $(Symbol(:A_,r,:_,c)))))
            end
        end
        if isU′
            Uoffset -= K + j - 1
#            push!(q.args, :($Ltrisym += $(size_T*(j+1))))
        end
    end
    push!(q.args, finish_block)
    if calc_product > 0 # calculate lower tirangle of A' * B
        # calc_product is an indicator
        # 0 indicates don't calc product
        # otherwise, it indicates number of columns to left of kernel
        # outer loop iterates over column of v∂L; we update these columns, iterating over rows of A and B inside to maximize out of order execution
        for c ∈ 0:C-1
            b2c = binomial2(c+1)
            # We load elements from v∂L of this column
            # diagonal
            diagv∂L = Symbol(:v∂L_, c, :_, c)
            push!(q.args, Expr(:(=), diagv∂L, :(vload($V, ptrv∂Ldiag + $(c*size_T*v∂Lstride)))))
            # now, the off diagonal
            for cc ∈ c:C-2
                subdiagv∂L = Symbol(:v∂L_, cc+1, :_, c)
                push!(q.args, Expr(:(=), subdiagv∂L, :(vload($V, ptrv∂Ltri + $(size_T*v∂Lstride)*($Kmax*$c + $(cc-c - b2c))))))
            end
            # Now, we load vectors from B one at a time, and multiply
            for r ∈ 0:Riterl
                mask_this_iter = maskload && r == Riterl && Rrem > 0
                # load
                set_bsym = Symbol(:B_,r,:_,c)
                # if mask_this_iter
                #     push!(q.args, :($set_bsym = vload($V, $Bsym + $size_T * ($(Riterl*W) + $c*$Bstride), $mask) ))
                # else
                    push!(q.args, :($set_bsym = vload($V, $Bsym + $size_T * ($(r*W) + $c*$Bstride)) ) )
                # end
                # multuplication; diag first
                if mask_this_iter
                    push!(q.args, Expr(:(=), diagv∂L, :(vifelse($mask,vmuladd($(Symbol(:A_,r,:_,c)),$set_bsym,$diagv∂L),$diagv∂L))))
                    for cc ∈ c+1:C-1
                        subdiagv∂L = Symbol(:v∂L_, cc, :_, c)
                        push!(q.args, Expr(:(=), subdiagv∂L, :(vifelse($mask,vmuladd($(Symbol(:A_,r,:_,cc)),$set_bsym,$subdiagv∂L),$subdiagv∂L))))
                    end
                else
                    push!(q.args, Expr(:(=), diagv∂L, :(vmuladd($(Symbol(:A_,r,:_,c)),$set_bsym,$diagv∂L))))
                    for cc ∈ c+1:C-1
                        subdiagv∂L = Symbol(:v∂L_, cc, :_, c)
                        push!(q.args, Expr(:(=), subdiagv∂L, :(vmuladd($(Symbol(:A_,r,:_,cc)),$set_bsym,$subdiagv∂L))))
                    end
                end
            end
            # Now, store the results
            push!(q.args, :(vstore!(ptrv∂Ldiag + $(c*size_T*v∂Lstride), $diagv∂L)))
            # now, store the off diagonal
            for cc ∈ c:C-2
                subdiagv∂L = Symbol(:v∂L_, cc+1, :_, c)
                push!(q.args, :(vstore!(ptrv∂Ltri + $(size_T*v∂Lstride)*($Kmax*$c + $(cc-c - b2c)), $subdiagv∂L)))
            end
        end
    end
    # take a break from calc_product while that part of A is still close in memory, in case A === B
    # Store in A
    if storeA
        if increment_store
            for c ∈ 0:C-1
                for r ∈ 0:Riter-1
                    tempptrA = gensym(Asym)
                    astore_q = quote
                        $tempptrA = $Asym + $size_T * ($(r*W) + $c*$Astride)
                        vstore!( $tempptrA, vadd(vload($V, $tempptrA), $(Symbol(:A_,r,:_,c)) ) )
                    end
                    push!(q.args, astore_q )
                end
                if Rrem > 0
                    tempptrA = gensym(Asym)
                    astore_q = quote
                        $tempptrA = $Asym + $size_T * ($(Riter*W) + $c*$Astride)
                        vstore!( $tempptrA, vadd(vload($V, $tempptrA, $mask), $(Symbol(:A_,Riter,:_,c)), $mask ) )
                    end
                    push!(q.args, astore_q)
                end
            end
        else
            for c ∈ 0:C-1
                for r ∈ 0:Riter-1
                    push!(q.args, :( vstore!( $Asym + $size_T * ($(r*W) + $c*$Astride), $(Symbol(:A_,r,:_,c)) ) ))
                end
                if Rrem > 0
                    push!(q.args, :( vstore!( $Asym + $size_T * ($(Riter*W) + $c*$Astride), $(Symbol(:A_,Riter,:_,c)), $mask ) ))
                end
            end
        end
    end
    if calc_product > 0 && (Kmax isa Symbol || (Kmax < calc_product))
        loopbody = quote end
        # load elements from v∂L
        for c ∈ 0:C-1
            subdiagv∂L = Symbol(:v∂L_, c, :_x)
            push!(loopbody.args, Expr(:(=), subdiagv∂L, :(vload($V, ptrv∂Ltri + $(v∂Lstride*size_T)*($c - decrement)))))
        end
        # Now, we load vectors from B one at a time, and multiply
        for r ∈ 0:Riterl
            # load
            set_bsym = Symbol(:B_,r,:_x)
            mask_this_iter = maskload && r == Riterl && Rrem > 0
            if mask_this_iter
                push!(loopbody.args, :($set_bsym = vload($V, $Bsym + $size_T * ($(Riter*W) - colleft*$Bstride), $mask) ))
            else
                push!(loopbody.args, :($set_bsym = vload($V, $Bsym + $size_T * ($(r*W) - colleft*$Bstride)) ) )
            end
            for c ∈ 0:C-1
                subdiagv∂L = Symbol(:v∂L_, c, :_x)
                if mask_this_iter
                    push!(loopbody.args, Expr(:(=), subdiagv∂L, :(vifelse($mask,vmuladd($(Symbol(:A_,r,:_,c)),$set_bsym,$subdiagv∂L),$subdiagv∂L) )))
                else
                    push!(loopbody.args, Expr(:(=), subdiagv∂L, :(vmuladd($(Symbol(:A_,r,:_,c)),$set_bsym,$subdiagv∂L))))
                end
            end
        end
        # now, store the sub-diagonal
        for c ∈ 0:C-1
            subdiagv∂L = Symbol(:v∂L_, c, :_x)
            push!(loopbody.args, :(vstore!(ptrv∂Ltri + $(v∂Lstride*size_T)*($c - decrement),$subdiagv∂L)))
        end

        calc_prod_quote = quote
            decrement = $Kmax
            numtrianglerows = $Kmax
            colleft = 0
            # do one column at a time
            while numtrianglerows < $calc_product
                colleft += 1
                $loopbody
                decrement += numtrianglerows
                numtrianglerows += 1
            end
        end
        push!(q.args, calc_prod_quote)
    end
    q
end

function div_triangle_loads_crfirst(ncol, aloads, colblock, ncolblock, colrem)
    necessary_loads = aloads*ncol + binomial2(ncol+1)
    # if we load the remainder first, it has to be reloaded on all future ncolblocks
    repeat_aloads = (colrem*ncolblock + binomial2(ncolblock)*colblock) * aloads
    necessary_loads + repeat_aloads
end
function div_triangle_loads_crfirst(ncol, aloads, colblock)#, ncolblock = div(ncol, colblock), ncolrem = rem(ncol, colblock))
#    if ncol == typemax(typeof(ncol))
#        ncol = colblock
#        ncolblock, colrem = 1, 0
#    else
        ncolblock, colrem = divrem(ncol, colblock)
#    end
    ncolblock, colrem = divrem(ncol, colblock)
    div_triangle_loads_crfirst(ncol, aloads, colblock, ncolblock, colrem)
end
function div_triangle_loads_crlast(ncol, aloads, colblock, ncolblock, colrem)
    necessary_loads = aloads*ncol + binomial2(ncol+1)
    repeat_aloads = colblock*binomial2(ncolblock)*aloads
    # if we load the remainder last, it has to load all (colblock*ncolblock) completed columns
    if colrem > 0
        repeat_aloads += aloads*colblock*ncolblock
    end
    necessary_loads + repeat_aloads
end
function div_triangle_loads_crlast(ncol, aloads, colblock)
    ncolblock, colrem = divrem(ncol, colblock)
    div_triangle_loads_crlast(ncol, aloads, colblock, ncolblock, colrem)
end

# returns:
# Number of loads, whether there is a remainder
function div_ul_loads(ncol, aloads, colblock, ncolblock, colrem)
    tfirst = div_triangle_loads_crfirst(ncol, aloads, colblock, ncolblock, colrem)
    twotfirst = tfirst + tfirst
    if colrem == 0
        twotfirst - (aloads * colblock), true, false
    else
        tlast = div_triangle_loads_crlast(ncol, aloads, colblock, ncolblock, colrem)
        alt_strat = tfirst + tlast - (aloads * colblock)
        primary_strat = twotfirst - (aloads * colrem)
        min(primary_strat, alt_strat), false, alt_strat < primary_strat
    end
end
function div_ul_loads(ncol, aloads, colblock)
    ncolblock, colrem = divrem(ncol, colblock)
    div_ul_loads(ncol, aloads, colblock, ncolblock, colrem)
end

not_max(x::T) where {T} = x != typemax(T)
function div_triangle_blocking_structure(
    rows = typemax(UInt), cols = typemax(UInt), ::Type{T} = Float64; reg_count = VectorizationBase.REGISTER_COUNT, reg_size = VectorizationBase.REGISTER_SIZE, verbose = false
) where {T}
    if cols == typemax(typeof(cols))
        W, rc, cc, mi = PaddedMatrices.pick_kernel_size(T, rows, cols,  W = reg_size ÷ sizeof(T), NREG = reg_count)
        return W, rc, cc
    end
    W = div(reg_size, sizeof(T))
    if not_max(rows)
        while 2rows < W
            W >>>= 1
        end
    end
    max_aloads = (reg_count >>> 1) - 1
    ratios = Vector{Float64}(undef, max_aloads)
    row_counts = Vector{Int}(undef, max_aloads)
    col_counts = Vector{Int}(undef, max_aloads)
    for aloads in 1:max_aloads
        ncol = min(div(reg_count - aloads - 1, aloads), cols)
        loads = div_triangle_loads_crfirst(cols, aloads, ncol)
        nrow = aloads * W
        ratio = aloads / loads
        if not_max(rows)
            complete_rows, row_rem = divrem(rows, nrow)
            if row_rem > 0
                rem_aloads = cld(row_rem, W)
                ratio = ( ratio * complete_rows * nrow + rem_aloads * row_rem / div_triangle_loads_crfirst(cols, rem_aloads, ncol) ) / rows
            end
        end
        verbose && @show nrow, ncol, ratio
        if cols % ncol == 0 # bonus for there being no column remainder.
            ratio *= 1.05
        end
        ratios[aloads] = ratio
        row_counts[aloads] = nrow
        col_counts[aloads] = ncol
        if nrow >= rows
            ratios[aloads+1:end] .= 0
            break
        end
    end
    bestratio, bestind = findmax(ratios)
    W, row_counts[bestind], col_counts[bestind]
end

function div_ul_blocking_structure(rows = typemax(UInt), cols = typemax(UInt), ::Type{T} = Float64; reg_count = VectorizationBase.REGISTER_COUNT, reg_size = VectorizationBase.REGISTER_SIZE) where {T}
    W = div(reg_size, sizeof(T))
    if not_max(rows)
        while 2rows < W
            W >>>= 1
        end
    end
    max_aloads = (reg_count >>> 1) - 1
    ratios = Vector{Float64}(undef, max_aloads)
    row_counts = Vector{Int}(undef, max_aloads)
    col_counts = Vector{Int}(undef, max_aloads)
    # is the remainder 0, are we using the alt strat (doing the remainder first then last) instead of remainder first for both divisions
    strategy = Vector{Tuple{Bool,Bool}}(undef, max_aloads)
    if not_max(cols)
        for aloads in 1:max_aloads
            ncol = min(div(reg_count - aloads - 1, aloads), cols)
            loads, b1, b2 = div_ul_loads(cols, aloads, ncol)
            strategy[aloads] = (b1,b2)
            ratio = aloads / loads
            nrow = aloads * W
            if not_max(rows)
                complete_rows, row_rem = divrem(rows, nrow)
                if row_rem > 0
                    rem_aloads = cld(row_rem, W)
                    rem_loads, b1, b2 = div_ul_loads(cols, rem_aloads, ncol)
                    ratio = ( ratio * complete_rows * nrow + row_rem * rem_aloads / rem_loads ) / rows
                end
            end
            # @show aloads, ncol, ratio
            if cols % ncol == 0
                ratio *= 1.05
            end
            ratios[aloads] = ratio
            row_counts[aloads] = nrow
            col_counts[aloads] = ncol
            if nrow >= rows
                ratios[aloads+1:end] .= 0
                break
            end
        end
    else
        for aloads in 1:max_aloads
            ncol = div(reg_count - aloads - 1, aloads)
            # ncol = min(cols, ncol) # but cols is max, so of course ncol < col?
            nrow = aloads * W
            row_counts[aloads] = nrow
            col_counts[aloads] = ncol
            loads = aloads + ncol
            ratio = aloads * ncol / loads
            if not_max(rows)
                complete_rows, row_rem = divrem(rows, nrow)
                if row_rem > 0
                    rem_aloads = cld(row_rem, W)
                    ratio = ( ratio * complete_rows * nrow + row_rem * (rem_aloads * ncol) / (rem_aloads + ncol) ) / rows
                end
            end
            ratios[aloads] = ratio
            if nrow >= rows
                ratios[aloads+1:end] .= 0
                break
            end
        end
    end
    bestratio, bestind = findmax(ratios)
    W, row_counts[bestind], col_counts[bestind], strategy[bestind]
end

function A_div_U_rowiter(
    Mk, Nk, col_rem, ::Type{T}, CP, AP, n_col_reps
) where {T}
    size_T = sizeof(T)
    config = A_rdiv_U_config{T}()
    config.Mk = Mk
    config.C = col_rem
    config.K = 0
    config.Astride = CP
    config.Bstride = AP
    config.Ustride = 0
    config.isL′ = false
    config.invdiagptr = true
    if col_rem > 0
        row_iter = A_rdiv_U_kernel_quote( config )
        push!(row_iter.args, :(ptrUtri += $(binomial2(col_rem)*size_T)))
        push!(row_iter.args, :(ptrUdiag += $(col_rem*size_T)))
        base_K = col_rem
    else
        row_iter = quote end
        base_K = 0
    end
    if n_col_reps > 1
        config.C = Nk
        config.K = :K
        config.Ustride = :K
        iterquote = A_rdiv_U_kernel_quote( config )
        row_iter_loop = quote
            K = $base_K
            for crep ∈ 0:$(n_col_reps-1)
                $iterquote
                ptrUtri += $(size_T*binomial2(Nk)) + $(size_T*Nk)*K
                ptrUdiag += $(size_T*Nk)
                K += $Nk
            end
        end
        push!(row_iter.args, row_iter_loop)
    elseif n_col_reps == 1
        config.C = Nk
        config.K = 0
        config.Ustride = 0
        row_iter_single = A_rdiv_U_kernel_quote( config )
        push!(row_iter.args, row_iter_single)
    end
    row_iter
end
function A_rdiv_U_quote(
    M, N, T, CP, AP, NBL
)
    # W = vector width
    # Mk = kernel rows
    # Nk = kernel colums
    W, Mk, Nk = div_triangle_blocking_structure(M, N, T)
    Wm1 = W - 1
    n_row_reps, row_rem = divrem(M, Mk)
    total_row_iterations = n_row_reps + (row_rem > 0)
    n_col_reps, col_rem = divrem(N, Nk)
    total_col_iterations = n_col_reps + (col_rem > 0)
    Nl = ( N + Wm1 ) & ~Wm1
    Nl = Nl > NBL ? N : Nl # Don't segfault
    size_T = sizeof(T)
    q = quote
        invdiag = FixedSizeVector{$N,$T,$Nl,$Nl}(undef)
        LoopVectorization.@vvectorize $T for n ∈ 1:$Nl
            invdiag[n] = one($T) / B[n]
        end
        ptrB = pointer(A)
        ptrA = pointer(C)
    end
    Mk1 = n_row_reps == 0 ? row_rem : Mk
    row_iter = A_div_U_rowiter(
        Mk1, Nk, col_rem, T, CP, AP, n_col_reps
    )
    if n_row_reps > 1
        row_loops = quote
            for rrep ∈ 1:$n_row_reps
                ptrUdiag = pointer(invdiag)
                ptrUtri = pointer(B) + $(size_T * N)
                $row_iter
                ptrB += $(size_T*Mk)
                ptrA += $(size_T*Mk)
            end
        end
        push!(q.args, row_loops)
    else
        push!(q.args, :(ptrUdiag = pointer(invdiag)))
        push!(q.args, :(ptrUtri = pointer(B) + $(size_T * N)))
        push!(q.args, row_iter)
        row_rem > 0 && push!(q.args, :(ptrB += $(size_T*Mk); ptrA += $(Size_T*Mk)))
    end
    if row_rem > 0 && n_row_reps > 0
        push!(q.args, :(ptrUdiag = pointer(invdiag)))
        push!(q.args, :(ptrUtri = pointer(B) + $(size_T * N)))
        push!(q.args, A_div_U_rowiter( row_rem, Nk, col_rem, T, CP, AP, n_col_reps ))
    end
    push!(q.args, :C)
    q
end
@generated function A_rdiv_B!(
    C::AbstractMutableFixedSizeMatrix{M,N,T,CP},
    A::AbstractMutableFixedSizeMatrix{M,N,T,AP},
    B::AbstractUpperTriangularMatrix{N,T,NBL}
#) where {M,N,T,NBL,CP,AP}
) where {M,N,T,CP,AP,NBL}
    A_rdiv_U_quote(M, N, T, CP, AP, NBL)
end

function A_div_L′_rowiter(
    Mk, Nk, col_rem, ::Type{T}, CP, AP, n_col_reps
) where {T}
    N = Nk * n_col_reps + col_rem
    size_T = sizeof(T)
    config = A_rdiv_U_config{T}()
    config.R = Mk
    config.C = col_rem
    config.K = 0
    config.Astride = CP
    config.Bstride = AP
    config.Ustride = N
    config.isL′ = true
    config.invdiagptr = true
    if col_rem > 0
        row_iter = A_rdiv_U_kernel_quote( config )
        push!(row_iter.args, :(ptrUdiag += $(col_rem*size_T)))
        base_K = col_rem
        KmZ = false
    else
        row_iter = quote end
        base_K = 0
        KmZ = true
    end
    if n_col_reps > 1
        config.C = Nk
        config.K = :K
        iterquote = A_rdiv_U_kernel_quote( config )
        row_iter_loop = quote
            K = $base_K
            for crep ∈ 0:$(n_col_reps-1)
                ptrUtri = ptrUtribase + K*$size_T
                $iterquote
                ptrUdiag += $(size_T*Nk)
                K += $Nk
            end
        end
        push!(row_iter.args, row_iter_loop)
    elseif n_col_reps == 1
        config.C = Nk
        config.K = col_rem
        row_iter_single = A_rdiv_U_kernel_quote( config )
        push!(row_iter.args, row_iter_single)
    end
    row_iter
end
function A_rdiv_L′_quote(
    M, N, T, CP, AP, NBL
)
    # W = vector width
    # Mk = kernel rows
    # Nk = kernel colums
    W, Mk, Nk = div_triangle_blocking_structure(M, N, T)
    Wm1 = W - 1
    n_row_reps, row_rem = divrem(M, Mk)
    total_row_iterations = n_row_reps + (row_rem > 0)
    n_col_reps, col_rem = divrem(N, Nk)
    total_col_iterations = n_col_reps + (col_rem > 0)
    Nl = ( N + W - 1 ) & ~Wm1
    Nl = Nl > NBL ? N : Nl # Don't segfault
    size_T = sizeof(T)
    q = quote
        B = Badj.parent
        invdiag = FixedSizeVector{$N,$T,$Nl,$Nl}(undef)
        LoopVectorization.@vvectorize $T for n ∈ 1:$Nl
            invdiag[n] = one($T) / B[n]
        end
        ptrB = pointer(A)
        ptrA = pointer(C)
        ptrUtribase = pointer(B) + $(N*size_T)
    end
    Mk1 = n_row_reps == 0 ? row_rem : Mk
    row_iter = A_div_L′_rowiter(
        Mk1, Nk, col_rem, T, CP, AP, n_col_reps
    )
    if n_row_reps > 1
        row_loops = quote
            for rrep ∈ 1:$n_row_reps
                ptrUdiag = pointer(invdiag)
                ptrUtri = ptrUtribase#pointer(B) + $(size_T * N)
                $row_iter
                ptrB += $(size_T*Mk)
                ptrA += $(size_T*Mk)
            end
        end
        push!(q.args, row_loops)
    else
        push!(q.args, :(ptrUdiag = pointer(invdiag)))
        push!(q.args, :(ptrUtri = ptrUtribase))
        push!(q.args, row_iter)
        push!(q.args, :(ptrB += $(size_T*Mk); ptrA += $(size_T*Mk) ))
    end
    if row_rem > 0 && n_row_reps > 0 # then n_row_reps == 1 and row_rem > 0
        push!(q.args, :(ptrUdiag = pointer(invdiag)))
        push!(q.args, :(ptrUtri = ptrUtribase))
        push!(q.args, A_div_L′_rowiter( row_rem, Nk, col_rem, T, CP, AP, n_col_reps ))
    end
    push!(q.args, :C)
    q
end
@generated function A_rdiv_B!(
    C::AbstractMutableFixedSizeMatrix{M,N,T,CP},
    A::AbstractMutableFixedSizeMatrix{M,N,T,AP},
    Badj::Adjoint{T,<:AbstractLowerTriangularMatrix{N,T,NBL}}
# ) where {M,N,T,CP,AP,NBL}
) where {M,N,T,NBL,CP,AP}
    A_rdiv_L′_quote(
        M, N, T, CP, AP, NBL
    )
end

function A_div_L_rowiter(
    Mk, Nk, col_rem, ::Type{T}, CP, AP, n_col_reps
) where {T}
    N = Nk * n_col_reps + col_rem
    size_T = sizeof(T)
    config = A_rdiv_L_config{T}()
    config.R = Mk
    config.C = col_rem
    config.Kmax = col_rem
    config.Astride = CP
    config.Bstride = AP
    config.isU′ = false
    config.invdiagptr = true
    if col_rem > 0
        row_iter = A_rdiv_L_kernel_quote( config )
        #pushfirst!(row_iter.args, :(ptrUtri = ptrUtribase))
        fullcols = Nk * n_col_reps
        # handle following in A_rdiv_L_quote
        # pushfirst!(row_iter.args, :(ptrA -= $(CP*col_rem*size_T)))
        # pushfirst!(row_iter.args, :(ptrB -= $(AP*col_rem*size_T)))
        push!(row_iter.args, :(ptrLdiag -= $(col_rem*size_T)))
        push!(row_iter.args, :(ptrLtri -= $((binomial2(Nk) + Nk*col_rem)*size_T)))
        base_K = col_rem
        KmZ = false
    else
        row_iter = quote end
        base_K = 0
        KmZ = true
    end
    if n_col_reps > 1
        config.C = Nk
        config.Kmax = :K
        iterquote = A_rdiv_L_kernel_quote( config )
        row_iter_loop = quote
            K = $col_rem
            for crep ∈ 0:$(n_col_reps-1)
                K += $Nk
                $iterquote
                ptrLdiag -= $(size_T*Nk)
                #K += $Nk
                ptrA -= $(Nk*CP*size_T)
                ptrB -= $(Nk*AP*size_T)
                ptrLtri -= $size_T*($Nk*K + $(binomial2(Nk)))  # = ptrLtribase + K*$size_T
            end
        end
        push!(row_iter.args, row_iter_loop)
    elseif n_col_reps == 1
        config.C = Nk
        config.Kmax = N
        row_iter_single = A_rdiv_L_kernel_quote( config )
        push!(row_iter.args, row_iter_single)
    end
    row_iter
end
function A_rdiv_L_quote(
    M, N, T, CP, AP, NBL
)
    # W = vector width
    # Mk = kernel rows
    # Nk = kernel colums
    W, Mk, Nk = div_triangle_blocking_structure(M, N, T)
    Wm1 = W - 1
    n_row_reps, row_rem = divrem(M, Mk)
    total_row_iterations = n_row_reps + (row_rem > 0)
    n_col_reps, col_rem = divrem(N, Nk)
    total_col_iterations = n_col_reps + (col_rem > 0)
    Nl = ( N + W - 1 ) & ~Wm1
    Nl = Nl > NBL ? N : Nl # Don't segfault
    size_T = sizeof(T)
    startoffset = (total_col_iterations-1) * Nk
    q = quote
        invdiag = FixedSizeVector{$N,$T,$Nl,$Nl}(undef)
        LoopVectorization.@vvectorize $T for n ∈ 1:$Nl
            invdiag[n] = one($T) / B[n]
        end
        ptrB_base = pointer(A) + $(size_T*AP*startoffset)
        ptrA_base = pointer(C) + $(size_T*CP*startoffset)
        ptrLtribase = pointer(B) + $(size_T * (N + binomial2(startoffset) + startoffset * (N - startoffset))) # diag + triangle + subtriangle
        ptrLdiagbase = pointer(invdiag) + $(size_T * startoffset)
    end
    Mk1 = n_row_reps == 0 ? row_rem : Mk
    row_iter = A_div_L_rowiter(
        Mk1, Nk, col_rem, T, CP, AP, n_col_reps
    )
    if n_row_reps > 1
        row_loops = quote
            for rrep ∈ 1:$n_row_reps
                ptrLdiag = ptrLdiagbase
                ptrLtri = ptrLtribase
                ptrA = ptrA_base
                ptrB = ptrB_base
                $row_iter
                ptrB_base += $(size_T*Mk)
                ptrA_base += $(size_T*Mk)
            end
        end
        push!(q.args, row_loops)
    else
        push!(q.args, :(ptrLdiag = ptrLdiagbase))
        push!(q.args, :(ptrLtri = ptrLtribase))
        push!(q.args, :(ptrA = ptrA_base))
        push!(q.args, :(ptrB = ptrB_base))
        push!(q.args, row_iter)
    end
    if row_rem > 0 && n_row_reps > 0
        push!(q.args, :(ptrLdiag = ptrLdiagbase))
        push!(q.args, :(ptrLtri = ptrLtribase))
        push!(q.args, :(ptrA = ptrA_base + $(size_T*Mk)))
        push!(q.args, :(ptrB = ptrB_base + $(size_T*Mk)))
        push!(q.args, A_div_L_rowiter( row_rem, Nk, col_rem, T, CP, AP, n_col_reps ))
    end
    push!(q.args, :C)
    q
end

@generated function A_rdiv_B!(
    C::AbstractMutableFixedSizeMatrix{M,N,T,CP},
    A::AbstractMutableFixedSizeMatrix{M,N,T,AP},
    B::AbstractLowerTriangularMatrix{N,T,NBL}
# ) where {M,N,T,CP,AP,NBL}
) where {M,N,T,NBL,CP,AP}
    A_rdiv_L_quote(
        M, N, T, CP, AP, NBL
    )
end
function A_rdiv_B_expr(
    C::AbstractMutableFixedSizeMatrix{M,N,T,CP},
    A::AbstractMutableFixedSizeMatrix{M,N,T,AP},
    B::AbstractLowerTriangularMatrix{N,T,NBL}
) where {M,N,T,CP,AP,NBL}
#) where {M,N,T,NBL,CP,AP}
    PaddedMatrices.prettify(A_rdiv_L_quote(
        M, N, T, CP, AP, NBL
    ))
end

function A_rdiv_B!(
    C::AbstractMutableFixedSizeMatrix{M,N,T},
    A::AbstractMutableFixedSizeMatrix{M,N,T},
    B::Adjoint{T,<:AbstractUpperTriangularMatrix{N,T}}
) where {M,N,T,CP,AP,NBL}
#) where {M,N,T,NBL,CP,AP}

end

# function A_rdiv_B!(
#     C::AbstractMutableFixedSizeMatrix{M,N,T},
#     A::AbstractMutableFixedSizeMatrix{M,N,T},
#     B::AbstractSymmetricMatrixU{N,T}
# ) where {M,N,T,CP,AP,NBL}
# #) where {M,N,T,NBL,CP,AP}
# end
# function A_rdiv_B!(
#     C::AbstractMutableFixedSizeMatrix{M,N,T},
#     A::AbstractMutableFixedSizeMatrix{M,N,T},
#     B::AbstractSymmetricMatrixL{N,T}
# ) where {M,N,T,CP,AP,NBL}
# #) where {M,N,T,NBL,CP,AP}
# end


text:
"""
# Algorithms with numbers

- Basic arithmetic
    - Addition
    - Bases and logs
    - Multiplication and division
- Modular arithmetic
    - Substitution rule
    - Modular addition and multiplication
    - Two's complement
    - Modular exponentiation
    - Euclid's algorithm for greatest common divisor
    - Euclid's rule
    - An extension of Euclid's algorithm
    - Modular division
    - Is your social security number a prime?
- Primality testing
    - Fermat's little theorem
    - Hey, that was group theory!
    - Generating random primes
    - Lagrange's prime number theorem
    - Carmichael numbers
    - Randomized algorithms: a virtual chapter
- Cryptography
    - An application of number theory?
    - Private-key schemes: one-time pad and AES
    - RSA
- Universal Hashing
    - Hash tables
    - Families of hash functions
"""
module Chapter1

export
    multiply,
    divide,
    modexp,
    euclid,
    extendedeuclid,
    primality,
    primality2

"""
    multiply(x, y)

Compute the product of two integers `x` and `y`.

# Examples
```jldoctest
julia> multiply(7, 8)
56
```
"""
function multiply(x::Integer, y::Integer)
    if y == 0
        return 0
    end

    z = multiply(x, y ÷ 2)
    if iseven(y)
        return 2z
    else
        return x + 2z
    end
end

"""
    divide(x, y)

Compute the quotient and remainder of `x ÷ y`.

# Examples
```jldoctest
julia> divide(5, 3)
(1, 2)
```
"""
function divide(x::Integer, y::Integer)
    if x == 0
        return (q,r) = (0,0)
    end

    (q,r) = divide(x ÷ 2, y)
    q = 2q
    r = 2r
    if isodd(x)
        r = r + 1
    end
    if r >= y
        r = r - y
        q = q + 1
    end
    return (q,r)
end

"""
    modexp(x, y, N)

Compute the modular exponentation `x^y mod N`.

# Examples
```jldoctest
julia> modexp(2, 3, 5)
3
```
"""
function modexp(x::Integer, y::Integer, N::Integer)
    if y == 0
        return 1
    end

    z = modexp(x, y ÷ 2, N)
    if iseven(y)
        return mod(z^2, N)
    else
        return mod(x * z^2, N)
    end
end

"""
    euclid(a, b)

Using Euclid's algorithm, compute the greatest common divisor
between two integers `a` and `b` where `a >= b >= 0`.

# Examples
```jldoctest
julia> euclid(25, 20)
5
```
"""
function euclid(a::Integer, b::Integer)
    @assert a >= b >= 0
    if b == 0
        return a
    end

    return euclid(b, mod(a, b))
end

"""
    extendedeuclid(a, b)

A simple extension of Euclid's algorithm.

- Input: Two positive integers `a` and `b` with `a >= b >= 0`.
- Ouput: Integers `x`, `y`, `d` such that `d = gcd(a,b)` and `ax + by = d`.

# Examples
```jldoctest
julia> extendedeuclid(9, 6)
(1, -1, 3)
```
"""
function extendedeuclid(a::Integer, b::Integer)
    @assert a >= b && b >= 0
    if b == 0
        return (1, 0, a)
    end

    (x′, y′, d) = extendedeuclid(b, mod(a, b))
    return (y′, x′ - a ÷ b * y′, d)
end

"""
    isprime(N)

Determine whether an integer `N` is prime using Fermat's little theorem.

# Examples
```jldoctest
julia> primality(5)
true
```
"""
function primality(N::Integer)
    if N < 2
        return false
    end

    a = rand(1:(N - 1))
    return powermod(a, N - 1, N) == 1
end

"""
    primality2(N)

Determine whether an integer `N` is prime with low error probability.
"""
function primality2(N::Integer)
    if N < 2
        return false
    end

    a = rand(1:(N - 1), N - 1)
    if all([powermod(a[i], N - 1, N) == 1 for i ∈ 1:(N - 1)])
        return true
    else
        return false
    end
end

end # module


text:
using BenchmarkTools, JLD2, Arpack, TSVD # run PROPACK, tsvd (remove k=) then change to svds and do ARPACK
import PROPACK.tsvd as tsvdpro

H1 = 0:4500
H2 = 0:4500

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
puls = load("puls_57600.jld2", "puls")
Puls_L = size(puls,1)
Hank = Array{Float64}(undef,length(H1)*Puls_L,length(H2))

#---------- Cell Definition -----------------#
# Cell = Construct("LG M50") #Alternative "Doyle_94"
# Init_SOC = 0.75

# #Arrhenius
# Cell.Const.T = 298.15
# Arr_Factor = (1/Cell.Const.T_ref-1/Cell.Const.T)/R

# #Set Cell Constants
# Cell.Const.SOC = Init_SOC
# Cell.RA.Tlen = 57600
# Cell.Const.κ = Cell.Const.κf(Cell.Const.ce0)*exp(Cell.Const.Ea_κ*Arr_Factor)
# Cell.RA.Nfft = Cell.RA.Nfft!(Cell.RA.Fs, Cell.RA.Tlen)
# Cell.RA.f = Cell.RA.f!(Cell.RA.Nfft)
# Cell.RA.s = Cell.RA.s!(Cell.RA.Fs,Cell.RA.Nfft,Cell.RA.f)
# Cell.Neg.β = Cell.Neg.β!(Cell.RA.s)
# Cell.Pos.β = Cell.Pos.β!(Cell.RA.s)

function SVD!(Hank,puls,Hlen1,Hlen2)
    A = T = P = tuple()
    Puls_L = size(puls,1)
    for lp1 in 1:length(Hlen2), lp2 in 1:length(Hlen1)
        Hank[Puls_L*(lp2-1)+1:Puls_L*lp2,lp1] .= @view puls[:,Hlen2[lp1]+Hlen1[lp2]+1]
    end

    A = @benchmark svds(Hank; nsv=6)[1]
    T = @benchmark tsvd(Hank, 6)
    P = @benchmark tsvdpro(Hank, k=6)
    return A,T,P
end

# @inline function DRA_loop(Cell)
#     x = Time = tuple()
#     for i in 1500:500:5000
#         Cell.RA.H1 = 0:i
#         Cell.RA.H2 = 0:i
#         x = @benchmark DRA(Cell,Cell.RA.s,Cell.RA.f)
#         Time = flatten_(Time,x)
#     end
#     return Time
# end

A,T,P = SVD!(Hank,puls,H1,H2)

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
#!/usr/bin/env julia

function process_group(answers::Vector{String})::Int
  h = zeros(Int, 'z' - 'a' + 1)

  for answer in answers
    for c in answer
      h[1 + c - 'a'] += 1
    end
  end

  # task one
  # count(a -> a > 0, h)
  ppl_cnt = length(answers)
  count(a -> a == ppl_cnt, h)
end

function main()
  group = Vector{String}()
  res = 0
  for line in eachline("./input")
    if length(line) == 0
      res += process_group(group)
      group = Vector{String}()
      continue
    end

    push!(group, line)
  end

  res += process_group(group)
  println(res)
end

main()