text:
"""
    refine_t!(t::Triangulation, split_indices = [], k)

Refine a triangulation, also updating the images of the simplices
in the triangulation.

`t.points::Array{Float64, 2} is a (n_triangulation_vertices x embeddingdim) sized
    array where each row is a vertex of the triangulation.
`t.simplex_inds::Array{Int, 2}` is an array of size
    (n_trinagulation_simplices x (embeddingdim+1)). Each row of simplex_indices contains
    the indices of the vertices (rows of the vertices array) furnishing the corresponding
    simplex.
`t.impoints::Array{Float64, 2} is a (n_triangulation_vertices x embeddingdim) sized
        array where each row is a vertex of the triangulation.
`split_indices::Array{Int, 1}` are the row numbers of simplex_indices indicating which
    simplices should be split.
"""
function refine_t!(t::Triangulation, split_indices::Vector{Int}, k::Int)

    # The number of simplices to split
    n_split_simplices = length(split_indices)
    n_points_beforesplit = size(t.points, 1)
    # Indices of the simplices that are not split
    untouched_inds = complementary(split_indices, size(t.simplex_inds, 1))

    # The dimension of the space
    E = size(t.points, 2)

    # Rules for forming the strictly new vertices of the subtriangulation
    rules, subtriangulation = simplicial_subdivision(k, E)

    # How many new vertices are created each split?
    nverts_persplit = size(rules, 1)

    # We need an array that can accomodate all of them. Each row in this
    # array will be a new vertex. Stacks of nverts_persplit * E arrays.
    # We have as many as we have simplices to split.
    new_vertices = zeros(nverts_persplit * n_split_simplices, E)
    new_imagevertices = zeros(nverts_persplit * n_split_simplices, E)

    # Fill the array by looping over all simplices that we want to split
    for i = 1:n_split_simplices
        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = nverts_persplit * (i -1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Get the vertices of the simplex currently being splitted. Each of the
        # nverts_persplit new vertices will be a linear combination
        # of these vertices. Each row is a vertex.
        vertices = t.points[t.simplex_inds[simplex_idx, :], :]
        imagevertices = t.impoints[t.simplex_inds[simplex_idx, :], :]

        # Generate the strictly new vertices for each sub
        for j = 1:nverts_persplit
            # The index for a particular new vertex of the ith new simplex
            ind_newvertex = ind + j
            # Go the jth new subsimplex of the ith splitted simplex. The entries of this
            # vector
            jth_subsimplex = rules[j, :]

            # Pick the corresponding original vertices with indices contained in rules[j, :]
            original_vertices = vertices[rules[j, :], :]
            original_vertices_image = imagevertices[rules[j, :], :]

            new_vertices[ind_newvertex, :] = sum(original_vertices, dims=1) ./ k
            new_imagevertices[ind_newvertex, :] = sum(original_vertices_image, dims=1) ./ k

        end
    end

    # Find the unique new vertices
    unique_newverts = unique(new_vertices, 1)
    unique_newimverts = Array{Float64}(size(unique_newverts, 1), E)

    Ind = Vector{Int}(size(new_vertices, 1))
    for i = 1:size(unique_newverts, 1)

        count = 0

        for j = 1:size(new_vertices, 1) # run over the possibly repeated vertices
            if unique_newverts[i, :] == new_vertices[j, :]
                Ind[j] = i

                count = count + 1
                if count == 1
                    unique_newimverts[i, :] = new_imagevertices[j, :]
                end
            end
        end
    end
    #@show unique(new_imagevertices, 1) - unique_newimverts
    #@show size(unique(new_imagevertices)), size(unique_newimverts)
    #assert(all(unique(new_imagevertices, 1) - unique_newimverts .≈ 0.0))

    # Start counting indices from the new vertices.
    Ind = Ind + n_points_beforesplit
    nsimplices_persplit = size(subtriangulation, 1)

    # The subsimplices formed by the splitting. Each row contains E + 1 indices referencing
    # the vertices furnishing that particular subsimplex (now found in the updated
    # t.points array).
    subtriangs_inds = Array{Int64}(nsimplices_persplit * n_split_simplices, E + 1)

    # For each simplex that we need to split,
    for i = 1:n_split_simplices
        # The beginning of the stack we need to fill.
        index = nsimplices_persplit * (i - 1)

        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = nverts_persplit * (i - 1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Pick the indices of the original vertices. Should be a column vector.
        inds_original_vertices = t.simplex_inds[simplex_idx, :]

        # Indices of the new vertices. Should be a column vector
        inds_newverts = Ind[(ind + 1):(ind + nverts_persplit)]
        inds_allverts = vcat(inds_original_vertices, inds_newverts)

        # Go through each of the subsimplices formed by splitting that simplex
        for j = 1:nsimplices_persplit
            subtriangs_inds[index + j, :] = inds_allverts[subtriangulation[j, :]]
        end
    end

    # Append new values
    t.points = vcat(t.points, unique_newverts)
    t.impoints = vcat(t.impoints, unique_newimverts)
    t.simplex_inds = vcat(t.simplex_inds[untouched_inds , :], subtriangs_inds)

    # Calculate centroids, radii and volumes for the new simplices
    t.centroids, t.radii = centroids_radii2(t.points, t.simplex_inds)
    t.centroids_im, t.radii_im = centroids_radii2(t.impoints, t.simplex_inds)

    @assert length(t.radii) == size(t.simplex_inds, 1)
    @assert size(t.centroids, 1) == size(t.simplex_inds, 1)

    t.volumes = simplex_volumes(t.points, t.simplex_inds)
    t.volumes_im = simplex_volumes(t.impoints, t.simplex_inds)

end

"""
    refine_t!(t::Triangulation, split_indices = [], k)

Refine a triangulation, also updating the images of the simplices
in the triangulation.

`t.points::Array{Float64, 2} is a (n_triangulation_vertices x embeddingdim) sized
    array where each row is a vertex of the triangulation.
`t.simplex_inds::Array{Int, 2}` is an array of size
    (n_trinagulation_simplices x (embeddingdim+1)). Each row of simplex_indices contains
    the indices of the vertices (rows of the vertices array) furnishing the corresponding
    simplex.
`t.impoints::Array{Float64, 2} is a (n_triangulation_vertices x embeddingdim) sized
        array where each row is a vertex of the triangulation.
`split_indices::Array{Int, 1}` are the row numbers of simplex_indices indicating which
    simplices should be split.
"""
function refine_t(t::Triangulation, split_indices::Vector{Int}, k::Int)

    # The number of simplices to split
    n_split_simplices = length(split_indices)
    n_points_beforesplit = size(t.points, 1)
    # Indices of the simplices that are not split
    untouched_inds = complementary(split_indices, size(t.simplex_inds, 1))

    # The dimension of the space
    E = size(t.points, 2)

    # Rules for forming the strictly new vertices of the subtriangulation
    rules, subtriangulation = simplicial_subdivision(k, E)

    # How many new vertices are created each split?
    nverts_persplit = size(rules, 1)

    # We need an array that can accomodate all of them. Each row in this
    # array will be a new vertex. Stacks of nverts_persplit * E arrays.
    # We have as many as we have simplices to split.
    new_vertices = zeros(nverts_persplit * n_split_simplices, E)
    new_imagevertices = zeros(nverts_persplit * n_split_simplices, E)

    # Fill the array by looping over all simplices that we want to split
    for i = 1:n_split_simplices
        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = nverts_persplit * (i - 1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Get the vertices of the simplex currently being splitted. Each of the
        # nverts_persplit new vertices will be a linear combination
        # of these vertices. Each row is a vertex.
        vertices = t.points[t.simplex_inds[simplex_idx, :], :]
        imagevertices = t.impoints[t.simplex_inds[simplex_idx, :], :]

        # Generate the strictly new vertices for each sub
        for j = 1:nverts_persplit
            # The index for a particular new vertex of the ith new simplex
            ind_newvertex = ind + j
            # Go the jth new subsimplex of the ith splitted simplex. The entries of this
            # vector
            jth_subsimplex = rules[j, :]

            # Pick the corresponding original vertices with indices contained in rules[j, :]
            original_vertices = vertices[rules[j, :], :]
            original_vertices_image = imagevertices[rules[j, :], :]

            new_vertices[ind_newvertex, :] = sum(original_vertices, dims=1) ./ k
            new_imagevertices[ind_newvertex, :] = sum(original_vertices_image, dims=1) ./ k

        end
    end

    # Find the unique new vertices
    unique_newverts = unique(new_vertices, 1)
    unique_newimverts = Array{Float64}(size(unique_newverts, 1), E)

    Ind = Vector{Int}(size(new_vertices, 1))
    for i = 1:size(unique_newverts, 1)

        count = 0

        for j = 1:size(new_vertices, 1) # run over the possibly repeated vertices
            if unique_newverts[i, :] == new_vertices[j, :]
                Ind[j] = i

                count = count + 1
                if count == 1
                    unique_newimverts[i, :] = new_imagevertices[j, :]
                end
            end
        end
    end
    #@show unique(new_imagevertices, 1) - unique_newimverts
    #@show size(unique(new_imagevertices)), size(unique_newimverts)
    #assert(all(unique(new_imagevertices, 1) - unique_newimverts .≈ 0.0))

    # Start counting indices from the new vertices.
    Ind = Ind + n_points_beforesplit
    nsimplices_persplit = size(subtriangulation, 1)

    # The subsimplices formed by the splitting. Each row contains E + 1 indices referencing
    # the vertices furnishing that particular subsimplex (now found in the updated
    # t.points array).
    subtriangs_inds = Array{Int64}(nsimplices_persplit * n_split_simplices, E + 1)

    # For each simplex that we need to split,
    for i = 1:n_split_simplices
        # Identify the beginning index of the stack we need to fill.
        index = nsimplices_persplit * (i - 1)

        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_vertices
        ind = nverts_persplit * (i - 1)

        # Index of the simplex we need to split
        simplex_idx = split_indices[i]

        # Pick the indices of the original vertices. Should be a column vector.
        inds_original_vertices = t.simplex_inds[simplex_idx, :]

        # Indices of the new vertices. Should be a column vector
        inds_newverts = Ind[(ind + 1):(ind + nverts_persplit)]
        inds_allverts = vcat(inds_original_vertices, inds_newverts)

        # Go through each of the subsimplices formed by splitting that simplex
        for j = 1:nsimplices_persplit
            subtriangs_inds[index + j, :] = inds_allverts[subtriangulation[j, :]]
        end
    end

    # Append new values
    t.points = vcat(t.points, unique_newverts)
    t.impoints = vcat(t.impoints, unique_newimverts)
    t.simplex_inds = vcat(t.simplex_inds[untouched_inds , :], subtriangs_inds)

    # Calculate centroids, radii and volumes for the new simplices
    t.centroids, t.radii = centroids_radii2(t.points, t.simplex_inds)
    t.centroids_im, t.radii_im = centroids_radii2(t.impoints, t.simplex_inds)

    @assert length(t.radii) == size(t.simplex_inds, 1)
    @assert size(t.centroids, 1) == size(t.simplex_inds, 1)

    t.volumes = simplex_volumes(t.points, t.simplex_inds)
    t.volumes_im = simplex_volumes(t.impoints, t.simplex_inds)

end


text:
################################################################################
# InteratomicPotentials API implementations for trainable potentials
################################################################################

function get_parameters(tp::TrainablePotential{P}) where {names,P<:NamedTuple{names}}
    P((getproperty(tp, name) for name ∈ names))
end
function set_parameters(tp::TP, p::P) where {P,TP<:TrainablePotential{P}}
    TP(; (name => getproperty(tp, name) for name ∈ propertynames(tp))..., p...)
end

function serialize_parameters(tp::TrainablePotential)
    # default beharior assumes parameters are all scalar, should be overloaded if this in not the case
    collect(get_parameters(tp))
end
function deserialize_parameters(tp::TrainablePotential{P}, p::AbstractVector) where {P}
    set_parameters(tp, P(p))
end

function get_hyperparameters(tp::TrainablePotential{P,HP}) where {P,names,HP<:NamedTuple{names}}
    HP((getproperty(tp, name) for name ∈ names))
end
function set_hyperparameters(tp::TP, p::HP) where {P,HP,TP<:TrainablePotential{P,HP}}
    TP(; (name => getproperty(tp, name) for name ∈ propertynames(tp))..., p...)
end

function serialize_hyperparameters(tp::TrainablePotential)
    # default beharior assumes hyperparameters are all scalar, should be overloaded if this in not the case
    collect(get_hyperparameters(tp))
end
function deserialize_hyperparameters(tp::TrainablePotential{P,HP}, p::AbstractVector) where {P,HP}
    set_hyperparameters(tp, HP(p))
end

ext: jl


text:
using Documenter

push!(LOAD_PATH, "..\\src\\")

using DelaunayMeshes

makedocs(
  format = :html,
  sitename = "DelaunayMeshes",
  pages = ["User Guide" => "UserGuide.md", "API" => "index.md"]
)

deploydocs(
    repo = "github.com/matthias-vigelius/DelaunayMeshes.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    julia = "release"
)


text:
function fastnorm(u)
  # dest[1] = ...
  n = length(u)
  T = eltype(u)
  s = zero(T)
  @fastmath @inbounds @simd for i in 1:n
      s += u[i]^2
  end
  @fastmath @inbounds return sqrt(s)
end

"""
    $TYPEDSIGNATURES

Calculate the Kernel Embedding MMD 'distance' between sample points (or kernel density estimates).

Notes
- `bw::Vector=[0.001;]` controls the mmd kernel bandwidths.

Related

`KDE.kld`
"""
function mmd( p1::AbstractMatrix{<:Real},
              p2::AbstractMatrix{<:Real},
              varType::Union{InstanceType{InferenceVariable},InstanceType{FunctorInferenceType}};
              bw::AbstractVector{<:Real}=[0.001;] )
  #
  manis = convert(AMP.Manifold, varType)
  mmd(p1, p2, manis, bw=bw)
end

function mmd( p1::BallTreeDensity,
              p2::BallTreeDensity,
              nodeType::Union{InstanceType{InferenceVariable},InstanceType{FunctorInferenceType}};
              bw::AbstractVector{<:Real}=[0.001;])
  #
  mmd(getPoints(p1), getPoints(p2), nodeType, bw=bw)
end

# moved to CalcFactor.jl

function sampleFactor(ccwl::CommonConvWrapper,
                      N::Int  )
  #
  cf = CalcFactor( ccwl.usrfnc!, _getFMdThread(ccwl), 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)
  sampleFactor(cf, N)
end

# part of consolidation, see #927
function sampleFactor!( ccwl::CommonConvWrapper,
                        N::Int,
                        fmd::FactorMetadata=_getFMdThread(ccwl),
                        vnd=nothing )
  #
  # depr warning added before IIF v0.20
  vnd !== nothing ? @warn("sampleFactor! no longer accepts vnd::Vector as meaningful input.") : nothing

  # build a CalcFactor object and get fresh samples.
  cf = CalcFactor( ccwl.usrfnc!, fmd, 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)
  # TODO make this an in-place operation as far possible
  ccwl.measurement = sampleFactor(cf, N)

  nothing
end

function sampleFactor(dfg::AbstractDFG,
                      sym::Symbol,
                      N::Int=1 )
  #
  # fct = getFactor(dfg, sym)
  return sampleFactor(_getCCW(dfg, sym), N)
  # usrfnc = getFactorType(fct)
  # variables = getVariable.(dfg, getVariableOrder(fct))

  # FIXME why not use ccwl.cpt[thrid].fmd rather than building a new one
  # fmd = _getFMdThread() # what about shared mem concurrancy?
  # FIXME fix/avoid getSample issue in testMultihypoFMD.jl: ummm, can we sample without knowing the hypotheses?
  # also deconv?
  # not really, because that would imply stochastic dependency on association before noise process??
  # fmd = FactorMetadata(variables, getLabel.(variables), Vector{Matrix{Float64}}(), :null, nothing)
  # sampleFactor(usrfnc, N, fmd )
end

"""
    $(SIGNATURES)

Update cliq `cliqID` in Bayes (Juction) tree `bt` according to contents of `urt` -- intended use is to update main clique after a upward belief propagation computation has been completed per clique.
"""
function updateFGBT!( fg::AbstractDFG,
                      cliq::TreeClique,
                      IDvals::Dict{Symbol, TreeBelief};
                      dbg::Bool=false,
                      fillcolor::String="",
                      logger=ConsoleLogger()  )
  #
  # if dbg
  #   # TODO find better location for the debug information (this is old code)
  #   cliq.attributes["debug"] = deepcopy(urt.dbgUp)
  # end
  if fillcolor != ""
    setCliqueDrawColor!(cliq, fillcolor)
  end
  for (id,dat) in IDvals
    with_logger(logger) do
      @info "updateFGBT! up -- update $id, inferdim=$(dat.inferdim)"
    end
    updvert = DFG.getVariable(fg, id)
    setValKDE!(updvert, deepcopy(dat), true) ## TODO -- not sure if deepcopy is required
  end
  with_logger(logger) do
    @info "updateFGBT! up -- updated $(getLabel(cliq))"
  end
  nothing
end

#


text:
#%%
using InteractiveUtils

"""
Julia 100-knock
"""

#%%
@show VERSION
versioninfo()
#%%
Z = zeros(10)
Z[5] = 1
@show Z
#%%
Z = [10:100]
#%%
Z = reshape(collect(1:9), (3, 3))
#%%
x=[1, 2, 0, 0, 4, 0]
nonzeroidx = collect(1:length(x))[iszero.(x)]
@show nonzeroidx

#%%

function gomagoma(n)
    for i in 1:n
        @show i
        j = i+1
        k = i-1
    end
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using MutableArithmetics
const MA = MutableArithmetics

function _resize_zeros!(v::Vector, new_len)
    old_len = length(v)
    if old_len < new_len
        resize!(v, new_len)
        for i in (old_len + 1):new_len
            v[i] = zero(eltype(v))
        end
    end
end

"""
    add_conv(out::Vector{T}, E::Vector{T}, k::Vector{T})
Returns the vector `out + fastconv(E, k)`. Note that only
`MA.mutable_buffered_operate!` is implemented.
"""
function add_conv end

# The buffer we need is the buffer needed by the `MA.add_mul` operation.
# For instance, `BigInt`s need a `BigInt` buffer to store `E[x] * k[i]` before
# adding it to `out[j]`.
function MA.buffer_for(::typeof(add_conv), ::Type{Vector{T}}, ::Type{Vector{T}}, ::Type{Vector{T}}) where {T}
    return MA.buffer_for(MA.add_mul, T, T, T)
end
function MA.mutable_buffered_operate!(buffer, ::typeof(add_conv), out::Vector{T}, E::Vector{T}, k::Vector{T}) where {T}
    for x in eachindex(E)
        for i in eachindex(k)
            j = x + i - 1
            out[j] = MA.buffered_operate!(buffer, MA.add_mul, out[j], E[x], k[i])
        end
    end
    return out
end

"""
    @register_mutable_arithmetic
Register polynomial type (with vector based backend) to work with MutableArithmetics
"""
macro register_mutable_arithmetic(name)
    poly = esc(name)
    quote
        MA.mutability(::Type{<:$poly}) = MA.IsMutable()

        function MA.promote_operation(::Union{typeof(+), typeof(*)},
                                      ::Type{$poly{S,X}}, ::Type{$poly{T,X}}) where {S,T,X}
            R = promote_type(S,T)
            return $poly{R,X}
        end

        function MA.buffer_for(::typeof(MA.add_mul),
                               ::Type{<:$poly{T,X}},
                               ::Type{<:$poly{T,X}}, ::Type{<:$poly{T,X}}) where {T,X}
            V = Vector{T}
            return MA.buffer_for(add_conv, V, V, V)
        end

        function MA.mutable_buffered_operate!(buffer, ::typeof(MA.add_mul),
                                              p::$poly, q::$poly, r::$poly)
            ps, qs, rs = coeffs(p), coeffs(q), coeffs(r)
            _resize_zeros!(ps, length(qs) + length(rs) - 1)
            MA.mutable_buffered_operate!(buffer, add_conv, ps, qs, rs)
            return p
        end
    end
end


text:
immutable ThreeBoxSITB <: Cut{3}; groups; end
export ThreeBoxSITB

function getcut(prob::Problem, vars::Variables, ::ThreeBoxSITB, i, j, k)
    cˣ, cʸ = vars.cˣ, vars.cʸ
    ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
    zˣ, zʸ = vars.zˣ, vars.zʸ

    lbˣ, lbʸ = prob.wlb, prob.hlb
    Lˣ, Lʸ = prob.W, prob.H
    @LinearConstraints(begin
        cˣ[i] ≥ 0.5ℓˣ[i] + lbˣ[k]*zˣ[k,i] + lbˣ[j]*(zˣ[j,k] + zˣ[k,i] - 1)
        cʸ[i] ≥ 0.5ℓʸ[i] + lbʸ[k]*zʸ[k,i] + lbʸ[j]*(zʸ[j,k] + zʸ[k,i] - 1)
        cˣ[i] + 0.5ℓˣ[i] + lbˣ[k]*zˣ[i,k] + lbˣ[j]*(zˣ[i,k] + zˣ[k,j] - 1) ≤ Lˣ
        cʸ[i] + 0.5ℓʸ[i] + lbʸ[k]*zʸ[i,k] + lbʸ[j]*(zʸ[i,k] + zʸ[k,j] - 1) ≤ Lʸ
    end)
end

immutable ThreeBoxNonoverlap <: Cut{3}; groups; end
export ThreeBoxNonoverlap

function getcut(prob::Problem, vars::Variables, ::ThreeBoxNonoverlap, i, j, k)
    cˣ, cʸ = vars.cˣ, vars.cʸ
    ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
    zˣ, zʸ = vars.zˣ, vars.zʸ

    lbˣ, lbʸ = prob.wlb, prob.hlb
    Lˣ, Lʸ = prob.W, prob.H
    @LinearConstraints(begin
        cˣ[i] - cˣ[j] + 0.5*(ℓˣ[i]+ℓˣ[j]) ≤ Lˣ*(1-zˣ[i,j]) - lbˣ[k]*(zˣ[i,k]+zˣ[k,j]-1)
        cʸ[i] - cʸ[j] + 0.5*(ℓʸ[i]+ℓʸ[j]) ≤ Lʸ*(1-zʸ[i,j]) - lbʸ[k]*(zʸ[i,k]+zʸ[k,j]-1)
    end)
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()
config_with_basis = MOIT.TestConfig(basis = true)

@testset "Split" begin
    T = Float64
    bridged_mock = MOIB.Constraint.SplitInterval{T}(mock)
    MOIT.basic_constraint_tests(bridged_mock, config, include=[
        (MOI.SingleVariable, MOI.Interval{T}),
        (MOI.ScalarAffineFunction{T}, MOI.Interval{T}),
        (MOI.ScalarQuadraticFunction{T}, MOI.Interval{T}),
        (MOI.SingleVariable, MOI.EqualTo{T}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        (MOI.ScalarQuadraticFunction{T}, MOI.EqualTo{T}),
        (MOI.VectorOfVariables, MOI.Zeros),
        (MOI.VectorAffineFunction{T}, MOI.Zeros),
        (MOI.VectorQuadraticFunction{T}, MOI.Zeros)
    ])
end

@testset "Interval" begin
    bridged_mock = MOIB.Constraint.SplitInterval{Float64}(mock)
    MOIU.set_mock_optimize!(mock,
         (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [5.0, 5.0], con_basis =
             [(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [MOI.NONBASIC],
              (MOI.SingleVariable, MOI.GreaterThan{Float64})                => [MOI.BASIC, MOI.BASIC]],
              (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [-1]),
         (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [2.5, 2.5], con_basis =
             [(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.NONBASIC],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [MOI.BASIC],
              (MOI.SingleVariable, MOI.GreaterThan{Float64})                => [MOI.BASIC, MOI.BASIC]],
              (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [1],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [0]),
         (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0],
                  (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [1],
                  (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [0],
              con_basis =
             [(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.NONBASIC],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [MOI.BASIC],
              (MOI.SingleVariable, MOI.GreaterThan{Float64})                => [MOI.BASIC, MOI.BASIC]]),
         (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [6.0, 6.0], con_basis =
             [(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [MOI.NONBASIC],
              (MOI.SingleVariable, MOI.GreaterThan{Float64})                => [MOI.BASIC, MOI.BASIC]]))
    MOIT.linear10test(bridged_mock, config_with_basis)
    MOIU.set_mock_optimize!(mock,
         (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.0, 0.0],
               (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0],
               (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [-1],
              con_basis =
             [(MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [MOI.BASIC],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [MOI.BASIC],
              (MOI.SingleVariable, MOI.GreaterThan{Float64})                => [MOI.NONBASIC, MOI.NONBASIC]],
              (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0],
              (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [0]))
    MOIT.linear10btest(bridged_mock, config_with_basis)

    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}()))
    vis = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    newf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], vis), 0.0)
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, newf)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ newf

    MOI.modify(bridged_mock, ci, MOI.ScalarCoefficientChange(vis[2], 1.0))
    modified_f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(2), vis), 0.0)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ modified_f

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, 2.0)
        @test MOI.get(bridged_mock, attr, ci) == 2.0
        bridge = MOIB.bridge(bridged_mock, ci)
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(bridged_mock, attr, bridge.lower) == 2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == 2.0
        else
            @test MOI.get(bridged_mock, attr, bridge.lower) == 2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == 0.0
        end
        MOI.set(bridged_mock, attr, ci, -2.0)
        @test MOI.get(bridged_mock, attr, ci) == -2.0
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(bridged_mock, attr, bridge.lower) == -2.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == -2.0
        else
            @test MOI.get(bridged_mock, attr, bridge.lower) == 0.0
            @test MOI.get(bridged_mock, attr, bridge.upper) == -2.0
        end
    end

    test_delete_bridge(bridged_mock, ci, 2, ((MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
                                            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64},    0)))
end
@testset "EqualTo{$T}" for T in [Float64, Int]
    bridged_mock = MOIB.Constraint.SplitInterval{T}(mock)
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [one(T), one(T)],
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}) => zeros(T, 2),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T})    => zeros(T, 1)))
    MOIT.linear13test(bridged_mock, MOIT.TestConfig{T}())

    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()))

    test_delete_bridge(bridged_mock, ci, 2, (
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}, 1),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T},    0)))
end
@testset "Zeros" begin
    bridged_mock = MOIB.Constraint.SplitInterval{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
        [1.0, 0.0, 2.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.0, 0.0]],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) => [[-3, -1]])
    MOIT.lin1vtest(bridged_mock, config)

    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}()))

    test_delete_bridge(bridged_mock, ci, 3, (
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),
        (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives, 0)))
end


text:
function f(N)
x = 0.
    for ii = 0:N-1
        x = 0.5*x + mod(ii,10)
    end
return x
end

println("simple_iter:   N: $N")
o=@benchmark f(N)
println(o)

max_stars_repo_name: scienceopen/numba-examples

max_issues_repo_name: scienceopen/numba-examples

max_forks_repo_name: scienceopen/numba-examples


text:
# Functions to get and check values

"""
    has_focal_mechanism(event) -> ::Bool

Return `true` if `event` contains one or more focal mechanisms defined.
"""
has_focal_mechanism(event::Event) = !isempty(event.focal_mechanism)

"""
    has_magnitude(event) --> ::Bool

Return `true` if `event` has one or more magnitudes defined.
"""
has_magnitude(event::Event) = !isempty(event.magnitude)

"""
    has_origin(event) -> ::Bool

Return `true` if `event` has one or more origins defined.
"""
has_origin(event::Event) = !isempty(event.origin)

"""
    preferred_focal_mechanism(event; verbose=false) -> focal_mechanism

Return the preferred focal mechanism for an `event`.  This may be defined if
there is more than one focal mechanism given for an event, and the
`preferred_focal_mechanism_id` field is set.  If there is only one focal mechanism
for this `event`, then that is returned.  If there is no focal mechanism
associated with this event which matches the stated `preferred_focal_mechanism_id`,
then the first focal mechanism is returned, and a warning given is `verbose` is `true`.
"""
function preferred_focal_mechanism(e::Event; verbose=false)
    has_focal_mechanism(e) ||
        throw(ArgumentError("event contains no focal mechanisms"))
    length(e.focal_mechanism) == 1 && return first(e.focal_mechanism)
    preferred_id = e.preferred_focal_mechanism_id
    ind = findfirst(x -> x.public_id === preferred_id, e.focal_mechanism)
    focal_mechanism = if ind === nothing
        verbose &&
            @warn("no focal mechanism with preferred id; returning the first focal mechanism")
        first(e.focal_mechanism)
    else
        e.focal_mechanism[ind]
    end
    focal_mechanism
end

"""
    preferred_focal_mechanisms(quakeml::EventParameters; verbose=false) -> focal_mechanisms

Return the preferred focal mechanisms for the events contained within `quakeml`.
See [`preferred_magnitude`](@ref).
"""
preferred_focal_mechanisms(ep::EventParameters; verbose=false) =
    preferred_focal_mechanism.(ep.event; verbose=verbose)

"""
    preferred_magnitude(event; verbose=false) -> magnitude

Return the preferred magnitude for an `event`.  This may be defined if
there is more than one magnitude given for an event, and the
`preferred_magnitude_id` field is set.  If there is only one magnitude
for this `event`, then that is returned.  If there is no magnitude
associated with this event which matches the stated `preferred_magnitude_id`,
then the first magnitude is returned, and a warning given is `verbose` is `true`.
"""
function preferred_magnitude(e::Event; verbose=false)
    has_magnitude(e) || throw(ArgumentError("event contains no magnitudes"))
    length(e.magnitude) == 1 && return first(e.magnitude)
    preferred_id = e.preferred_magnitude_id
    ind = findfirst(x -> x.public_id === preferred_id, e.magnitude)
    magnitude = if ind === nothing
        verbose && @warn("no magnitude with preferred id; returning the first magnitude")
        first(e.magnitude)
    else
        e.magnitude[ind]
    end
    magnitude
end

"""
    preferred_magnitueds(quakeml::EventParameters; verbose=false) -> magnitudes

Return the preferred magnitudes for the events contained within `quakeml`.
See [`preferred_magnitude`](@ref).
"""
preferred_magnitudes(ep::EventParameters; verbose=false) =
    preferred_magnitude.(ep.event; verbose=verbose)

"""
    preferred_origin(event; verbose=false) -> origin

Return the preferred origin for an `event`.  This may be defined if
there is more than one origin given for an event, and the `preferred_origin_id`
field is set.  If there is only one origin for this `event`, then that is
returned.  If there is no origin associated with this event which matches
the stated `preferred_origin_id`, then the first origin is returned
and a warning is given when `verbose=true`
"""
function preferred_origin(e::Event; verbose=false)
    has_origin(e) || throw(ArgumentError("event contains no origins"))
    length(e.origin) == 1 && return first(e.origin)
    preferred_id = e.preferred_origin_id
    ind = findfirst(x -> x.public_id === preferred_id, e.origin)
    origin = if ind === nothing
        verbose && @warn("no origin with preferred id; returning the first origin")
        first(e.origin)
    else
        e.origin[ind]
    end
    origin
end

"""
    preferred_origins(quakeml::EventParameters; verbose=false) -> origins

Return the preferred origins for the events contained within `quakeml`.
See [`preferred_origin`](@ref).
"""
preferred_origins(ep::EventParameters; verbose=false) =
    preferred_origin.(ep.event; verbose=verbose)

max_stars_repo_name: anowacki/QuakeML


text:
"""
    MulStyle

indicates that an `Applied` object should be materialized
via `Mul`.
"""
struct MulStyle <: AbstractArrayApplyStyle end

Mul(M::Applied) = Mul(M.args...)
arguments(M::Mul) = (M.A, M.B)

similar(M::Applied{MulStyle}, ::Type{T}) where T = similar(Mul(M), T)
copy(M::Applied{MulStyle}) = mul(arguments(M)...)
@inline copyto!(dest::AbstractArray, M::Applied{MulStyle}) = copyto!(dest, Mul(M))

const MulArray{T, N, Args} = ApplyArray{T, N, typeof(*), Args}

const MulVector{T, Args} = MulArray{T, 1, Args}
const MulMatrix{T, Args} = MulArray{T, 2, Args}

check_applied_axes(A::Applied{<:Any,typeof(*)}) = check_mul_axes(A.args...)

size(M::Applied{<:Any,typeof(*)}, p::Int) = size(M)[p]
axes(M::Applied{<:Any,typeof(*)}, p::Int) = axes(M)[p]
ndims(M::Applied{<:Any,typeof(*)}) = ndims(last(M.args))

_mul_ndims(::Type{Tuple{A}}) where A = ndims(A)
_mul_ndims(::Type{Tuple{A,B}}) where {A,B} = ndims(B)
ndims(::Type{<:Applied{<:Any,typeof(*),Args}}) where Args = _mul_ndims(Args)

length(M::Applied{<:Any,typeof(*)}) = prod(size(M))
size(M::Applied{<:Any,typeof(*)}) = length.(axes(M))

@inline _eltypes() = tuple()
@inline _eltypes(A, B...) = tuple(eltype(A), _eltypes(B...)...)

@inline eltype(M::Applied{<:Any,typeof(*)}) = _mul_eltype(_eltypes(M.args...)...)

@inline mulaxes1(::Tuple{}) = ()
@inline mulaxes1(::Tuple{}, B, C...) = mulaxes1(B, C...)
@inline mulaxes1(A::Tuple, C...) = first(A)
@inline mulaxes2(::Tuple{}) = ()
@inline mulaxes2(::Tuple{}, B, C...) = mulaxes2(B, C...)
@inline mulaxes2(A::Tuple{<:Any}, C...) = ()
@inline mulaxes2(A::Tuple{<:Any,<:Any}, C...) = last(A)

@inline _combine_axes(::Tuple{}, ::Tuple{}) = ()
@inline _combine_axes(a, ::Tuple{}) = (a,)
@inline _combine_axes(a, b) = (a,b)
@inline mulaxes(ax...) = _combine_axes(mulaxes1(ax...), mulaxes2(reverse(ax)...))

@inline axes(M::Applied{<:Any,typeof(*)}) = mulaxes(map(axes,M.args)...)
@inline axes(M::Applied{<:Any, typeof(*), Tuple{}}) = ()

###
# show
###

function _applyarray_summary(io::IO, ::typeof(*), args)
    print(io, "(")
    summary(io, first(args))
    print(io, ")")
    for a in tail(args)
        print(io, " * (")
        summary(io, a)
        print(io, ")")
    end
end

####
# Matrix * Array
####

combine_mul_styles(a) = a
combine_mul_styles(a, b) = error("Overload for $a and $b")
combine_mul_styles(::T, ::T) where T = T()
combine_mul_styles(::MulStyle, ::MulStyle) = DefaultArrayApplyStyle()
combine_mul_styles(::MulStyle, ::DefaultArrayApplyStyle) = DefaultArrayApplyStyle()
combine_mul_styles(::DefaultArrayApplyStyle, ::MulStyle) = DefaultArrayApplyStyle()
combine_mul_styles(::DefaultArrayApplyStyle, ::DefaultApplyStyle) = DefaultApplyStyle()
combine_mul_styles(::DefaultApplyStyle, ::DefaultArrayApplyStyle) = DefaultApplyStyle()
combine_mul_styles(a, b, c...) = combine_mul_styles(combine_mul_styles(a, b), c...)
# We need to combine all branches to determine whether it can be  simplified
ApplyStyle(::typeof(*), a) = DefaultApplyStyle()
ApplyStyle(::typeof(*), a::AbstractArray) = DefaultArrayApplyStyle()
_mul_ApplyStyle(a, b...) = combine_mul_styles(ApplyStyle(*, a, most(b)...),  ApplyStyle(*, b...))
ApplyStyle(::typeof(*), a, b...) = _mul_ApplyStyle(a, b...)
ApplyStyle(::typeof(*), a::Type{<:AbstractArray}, b::Type{<:AbstractArray}...) = _mul_ApplyStyle(a, b...)
ApplyStyle(::typeof(*), ::Type{<:AbstractArray}) = DefaultArrayApplyStyle()
ApplyStyle(::typeof(*), ::Type{<:Number}, ::Type{<:AbstractArray}) = DefaultArrayApplyStyle()
ApplyStyle(::typeof(*), ::Type{<:AbstractArray}, ::Type{<:Number}) = DefaultArrayApplyStyle()
ApplyStyle(::typeof(*), ::Type{<:AbstractArray}, ::Type{<:AbstractArray}) = MulStyle()

# arguments for something that is a *
@inline _mul_arguments(::ApplyLayout{typeof(*)}, A) = arguments(A)
@inline _mul_arguments(::DualLayout{ApplyLayout{typeof(*)}}, A) = arguments(A)
@inline _mul_arguments(_, A) = (A,)
@inline _mul_arguments(A) = _mul_arguments(MemoryLayout(A), A)

@inline __flatten(A::Tuple{<:Any}, B::Tuple) = (A..., _flatten(B...)...)
@inline __flatten(A::Tuple, B::Tuple) = _flatten(A..., B...)
@inline __flatten(A::Tuple{<:Any}, ::Tuple{}) = A
@inline __flatten(A::Tuple, B::Tuple{}) = _flatten(A...)

@inline _flatten() = ()
@inline _flatten(A, B...) = __flatten(_mul_arguments(A), B)
@inline _flatten(A) = __flatten(_mul_arguments(A), ())
@inline flatten(A) = _mul(_flatten(_mul_arguments(A)...)...)

####
# MulArray
#####

_mul(A) = A
_mul(A,B,C...) = lazymaterialize(*,A,B,C...)

_mul_colsupport(j, Z) = colsupport(Z,j)
_mul_colsupport(j, Z::AbstractArray) = colsupport(Z,j)
_mul_colsupport(j, Z, Y...) = axes(Z,1) # default is return all

_mul_colsupport(j, Z::AbstractArray, Y...) = _mul_colsupport(colsupport(Z,j), Y...)

colsupport(B::Applied{<:Any,typeof(*)}, j) = _mul_colsupport(j, reverse(B.args)...)
colsupport(B::MulArray, j) = _mul_colsupport(j, reverse(B.args)...)

_mul_rowsupport(j, A) = rowsupport(A,j)
_mul_rowsupport(j, A::AbstractArray) = rowsupport(A,j)
_mul_rowsupport(j, A, B...) = axes(A,1) # default is return all
_mul_rowsupport(j, A::AbstractArray, B...) = _mul_rowsupport(rowsupport(A,j), B...)

rowsupport(B::Applied{<:Any,typeof(*)}, j) = _mul_rowsupport(j, B.args...)
rowsupport(B::MulArray, j) = _mul_rowsupport(j, B.args...)

function getindex(M::Applied{<:Any,typeof(*)}, k...)
    A,Bs = first(M.args), tail(M.args)
    B = _mul(Bs...)
    Mul(A, B)[k...]
end

_flatten(A::MulArray, B...) = _flatten(Applied(A), B...)
flatten(A::MulArray) = ApplyArray(flatten(Applied(A)))

adjoint(A::MulMatrix) = ApplyArray(*, reverse(map(adjoint,A.args))...)
transpose(A::MulMatrix) = ApplyArray(*, reverse(map(transpose,A.args))...)

###
# sub materialize
###

# determine rows/cols of multiplication
_mul_args_rowsupport(a,kr) = rowsupport(a,kr)
_mul_args_colsupport(a,kr) = colsupport(a,kr)
__mul_args_rows(kr, a) = (kr,)
__mul_args_rows(kr, a, b...) =
    (kr, __mul_args_rows(_mul_args_rowsupport(a,kr), b...)...)
_mul_args_rows(kr, a, b...) = __mul_args_rows(_mul_args_rowsupport(a,kr), b...)
__mul_args_cols(jr, z) = (jr,)
__mul_args_cols(jr, z, y...) =
    (__mul_args_cols(_mul_args_colsupport(z,jr), y...)..., jr)
_mul_args_cols(jr, z, y...) = __mul_args_cols(_mul_args_colsupport(z,jr), y...)

sublayout(::ApplyLayout{typeof(*)}, _...) = ApplyLayout{typeof(*)}()
# matrix-indexing loses the multiplication structure as we don't support tensor multiplication
sublayout(::ApplyLayout{typeof(*)}, ::Type{<:Tuple{AbstractMatrix}}) = UnknownLayout()
sublayout(::ApplyLayout{typeof(*)}, ::Type{<:Tuple{AbstractVector{<:CartesianIndex}}}) = UnknownLayout()

call(::ApplyLayout{typeof(*)}, V::SubArray) = *

function _mat_mul_arguments(args, (kr,jr)::Tuple{Any,Any})
    kjr = intersect.(_mul_args_rows(kr, args...), _mul_args_cols(jr, reverse(args)...))
    map(view, args, (kr, kjr...), (kjr..., jr))
end

_vec_mul_view(a...) = view(a...)
_vec_mul_view(a::AbstractVector, kr, ::Colon) = view(a, kr)

# this is a vector view of a MulVector
function _vec_mul_arguments(args, (kr,))
    kjr = intersect.(_mul_args_rows(kr, args...), _mul_args_cols(Base.OneTo(1), reverse(args)...))
    _vec_mul_view.(args, (kr, kjr...), (kjr..., :))
end

# this is a vector view of a MulMatrix
_vec_mul_arguments(args, (kr,jr)::Tuple{AbstractVector,Number}) =
    _mat_mul_arguments(args, (kr,jr))

# this is a row-vector view
_vec_mul_arguments(args, (kr,jr)::Tuple{Number,AbstractVector}) =
    _vec_mul_arguments(reverse(map(permutedims, args)), (jr,kr))

_mat_mul_arguments(V) = _mat_mul_arguments(arguments(parent(V)), parentindices(V))
_vec_mul_arguments(V) = _vec_mul_arguments(arguments(parent(V)), parentindices(V))

arguments(::ApplyLayout{typeof(*)}, V::SubArray{<:Any,2}) = _mat_mul_arguments(V)
arguments(::ApplyLayout{typeof(*)}, V::SubArray{<:Any,1}) = _vec_mul_arguments(V)

@inline sub_materialize(lay::ApplyLayout{typeof(*)}, V) = apply(*, map(sub_materialize, arguments(lay, V))...)

##
# adjoint Mul
##

adjointlayout(::Type, ::ApplyLayout{typeof(*)}) = ApplyLayout{typeof(*)}()
transposelayout(::ApplyLayout{typeof(*)}) = ApplyLayout{typeof(*)}()

call(::ApplyLayout{typeof(*)}, V::Adjoint) = *
call(::ApplyLayout{typeof(*)}, V::Transpose) = *

arguments(::ApplyLayout{typeof(*)}, V::Adjoint) = reverse(adjoint.(arguments(V')))
arguments(::ApplyLayout{typeof(*)}, V::Transpose) = reverse(transpose.(arguments(V')))

permutedims(A::ApplyArray{<:Any,2,typeof(*)}) = ApplyArray(*, reverse(map(permutedims, A.args))...)

##
# * specialcase
##

for op in (:*, :\)
    @eval broadcasted(::LazyArrayStyle{N}, ::typeof($op), a::Number, b::ApplyArray{<:Number,N,typeof(*)}) where N =
        ApplyArray(*, broadcast($op,a,first(b.args)), tail(b.args)...)
end

broadcasted(::LazyArrayStyle{N}, ::typeof(/), b::ApplyArray{<:Number,N,typeof(*)}, a::Number) where N =
        ApplyArray(*, most(b.args)..., broadcast(/,last(b.args),a))

for Typ in (:Lmul, :Rmul)
    @eval $Typ(M::Applied{<:Any,typeof(*)}) = $Typ(M.args...)
end

##
# L/Rmul
##

@inline ApplyArray(M::Lmul) = ApplyArray(*, M.A, M.B)
@inline ApplyArray(M::Rmul) = ApplyArray(*, M.A, M.B)
@inline ApplyArray(M::Mul) = ApplyArray(*, M.A, M.B)
@inline ApplyArray(M::Mul{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}) = ApplyArray(*, arguments(M.A)..., arguments(M.B)...)
@inline ApplyArray(M::Mul{ApplyLayout{typeof(*)}}) = ApplyArray(*, arguments(M.A)..., M.B)
@inline ApplyArray(M::Mul{<:Any,ApplyLayout{typeof(*)}}) = ApplyArray(*, M.A, arguments(M.B)...)

# Support QuasiArrays

lazymaterialize(::typeof(*), a::AbstractArray) = a
lazymaterialize(F::Function, args::AbstractArray...) = copy(ApplyArray(F, args...))
lazymaterialize(M::Mul) = lazymaterialize(*, M.A, M.B)

###
# Simplify
# Here we implement a simple routine for simplifying multiplication by expanding what can be expanded

#
# The method is given *(a, b, ..., y, z) for see if *(a, b, ..., y) can be simplified.
# If so, simplify and start over. If not, see if *(b, ..., y, z) can  be simplified. If so,
# simplify and start over. Of not, return a lazy version.
# In ContinuumArrays we use this for simplifying differential operators
###

@inline _or(::Val{true}, ::Val{true}) = Val(true)
@inline _or(::Val{true}, ::Val{false}) = Val(true)
@inline _or(::Val{false}, ::Val{true}) = Val(true)
@inline _or(::Val{false}, ::Val{false}) = Val(false)

@inline _not(::Val{true}) = Val(false)
@inline _not(::Val{false}) = Val(true)

@inline _islazy(::AbstractLazyLayout) = Val(true)
@inline _islazy(_) = Val(false)
@inline islazy(A) = _islazy(MemoryLayout(A))
@inline simplifiable(M::Mul) = _not(_or(islazy(M.A), islazy(M.B)))

@inline simplifiable(::typeof(*), a) = Val(false)
@inline simplifiable(::typeof(*), a, b) = simplifiable(Mul(a,b))
@inline simplifiable(::typeof(*), a...) = _most_simplifiable(*, simplifiable(*, most(a)...), a)
@inline _most_simplifiable(::typeof(*), ::Val{true}, a) = Val(true)
@inline _most_simplifiable(::typeof(*), ::Val{false}, a) = simplifiable(*, tail(a)...)

# Flatten first
@inline simplify(::typeof(*), args...) = _simplify(*, _flatten(args...)...)
@inline _simplify(::typeof(*), a, b) = _twoarg_simplify(*, simplifiable(*, a, b), a, b)
@inline _twoarg_simplify(::typeof(*), ::Val{false}, a, b) = lazymaterialize(*, a, b)
@inline _twoarg_simplify(::typeof(*), ::Val{true}, a, b) = mul(a,b)
@inline _simplify(::typeof(*), args...) = _most_simplify(simplifiable(*, most(args)...), args)
@inline _most_simplify(::Val{true}, args) = *(_mul_arguments(_simplify(*, most(args)...))..., last(args))
@inline _most_simplify(::Val{false}, args) = _tail_simplify(simplifiable(*, tail(args)...), args)
@inline _tail_simplify(::Val{true}, args) = *(first(args), _mul_arguments(_simplify(*, tail(args)...))...)
@inline _tail_simplify(::Val{false}, args) = lazymaterialize(*, args...)

simplify(M::Mul) = simplify(*, M.A, M.B)
simplify(M::Applied{<:Any,typeof(*)}) = simplify(*, arguments(M)...)

@inline copy(M::Mul{<:AbstractLazyLayout,<:AbstractLazyLayout}) = simplify(M)
@inline copy(M::Mul{<:AbstractLazyLayout}) = simplify(M)
@inline copy(M::Mul{<:Any,<:AbstractLazyLayout}) = simplify(M)
@inline copy(M::Mul{<:AbstractQLayout,<:AbstractLazyLayout}) = simplify(M)
@inline copy(M::Mul{<:AbstractLazyLayout,<:AbstractQLayout}) = simplify(M)

@inline simplifiable(M::Mul{<:DualLayout,<:AbstractLazyLayout,<:AbstractMatrix,<:AbstractVector}) = Val(true)
@inline copy(M::Mul{<:DualLayout,<:AbstractLazyLayout,<:AbstractMatrix,<:AbstractVector}) = copy(Dot(M))

applylayout(::Type{typeof(*)}, ::DualLayout{Lay}, args...) where Lay = DualLayout{typeof(applylayout(typeof(*), Lay(), args...))}()

#TODO: Why not all DiagonalLayout?
@inline simplifiable(M::Mul{<:DiagonalLayout{<:AbstractFillLayout}}) = Val(true)
@inline simplifiable(M::Mul{<:DiagonalLayout{<:AbstractFillLayout},<:DiagonalLayout{<:AbstractFillLayout}}) = Val(true)
@inline simplifiable(M::Mul{<:Any,<:DiagonalLayout{<:AbstractFillLayout}}) = Val(true)
@inline copy(M::Mul{<:DiagonalLayout{<:AbstractFillLayout},<:AbstractLazyLayout}) = copy(mulreduce(M))
@inline copy(M::Mul{<:AbstractLazyLayout,<:DiagonalLayout{<:AbstractFillLayout}}) = copy(mulreduce(M))
@inline copy(M::Mul{BroadcastLayout{typeof(*)},<:DiagonalLayout{<:AbstractFillLayout}}) = copy(mulreduce(M))

@inline simplifiable(M::Mul{<:Union{ZerosLayout,DualLayout{ZerosLayout}},<:Union{ZerosLayout,DualLayout{ZerosLayout}}}) = Val(true)
@inline simplifiable(M::Mul{<:Union{ZerosLayout,DualLayout{ZerosLayout}}}) = Val(true)
@inline simplifiable(M::Mul{<:Any,<:Union{ZerosLayout,DualLayout{ZerosLayout}}}) = Val(true)
@inline simplifiable(M::Mul{<:Union{ZerosLayout,DualLayout{ZerosLayout}},<:DiagonalLayout{<:AbstractFillLayout}}) = Val(true)
@inline simplifiable(M::Mul{<:DiagonalLayout{<:AbstractFillLayout},<:Union{ZerosLayout,DualLayout{ZerosLayout}}}) = Val(true)

# inv

function _inv(Lay::ApplyLayout{typeof(*)}, _, A)
    args = arguments(Lay, A)
    map(checksquare,args)
    *(reverse(map(inv, arguments(Lay, A)))...)
end

##
# getindex
##
_reverse_mul_vec(z) = z
_reverse_mul_vec(z, y, w...) = _reverse_mul_vec(y*z, w...)
function getindex(M::ApplyMatrix{<:Any,typeof(*)}, ::Colon, j::Integer)
    rargs = reverse(M.args)
    _reverse_mul_vec(first(rargs)[:,j], tail(rargs)...)
end