text:
module NAFEMS_examples
using FinEtools
using FinEtools
using FinEtools.AlgoDeforLinearModule: ssit

function NAFEMS_FV32_algo()
    println("""
    FV32: Cantilevered tapered membrane
    This is a test recommended by the National Agency for Finite Element Methods and
    Standards (U.K.): Test FV32 from NAFEMS publication TNSB, Rev. 3, “The Standard
    NAFEMS Benchmarks,” October 1990.

    Reference solution: 44.623	130.03	162.70	246.05	379.90	391.44 for the first
    six modes.
    """)

    t0 = time()

    E = 200*phun("GPA");
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    L = 10*phun("M");
    W0 = 5*phun("M");
    WL = 1*phun("M");
    H = 0.05*phun("M");
    nL, nW, nH = 8, 4, 1;# How many element edges per side?
    neigvs=20                   # how many eigenvalues
    Reffs = [44.623	130.03	162.70	246.05	379.90	391.44]

    fens,fes =H20block(1.0, 2.0, 1.0, nL, nW, nH)
    for i = 1:count(fens)
        xi, eta, theta = fens.xyz[i,:];
        eta = eta - 1.0
        fens.xyz[i,:] = [xi*L eta*(1.0 - 0.8*xi)*W0/2 theta*H/2];
    end
    # File =  "mesh.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastIso(MR, rho, E, nu, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,2)),
    material), "femm_mass"=>FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,3)),
    material))

    nl1 = selectnode(fens; plane=[1.0 0.0 0.0 0.0], thickness=H/1.0e4)
    ebc1 = FDataDict("node_list"=>nl1, "component"=>1, "displacement"=>0.0)
    ebc2 = FDataDict("node_list"=>nl1, "component"=>2, "displacement"=>0.0)
    ebc3 = FDataDict("node_list"=>nl1, "component"=>3, "displacement"=>0.0)

    nl4 = selectnode(fens; plane=[0.0 0.0 1.0 0.0], thickness=H/1.0e4)
    ebc4 = FDataDict("node_list"=>nl4, "component"=>3, "displacement"=>0.0)

    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region1], "essential_bcs"=>[ebc1 ebc2 ebc3 ebc4],
    "neigvs"=>neigvs)

    # Solve
    modeldata = FinEtools.AlgoDeforLinearModule.modal(modeldata)

    fs = modeldata["omega"]/(2*pi)
    println("Eigenvalues: $fs [Hz]")
    println("Percentage frequency errors: $((vec(fs[1:6]) - vec(Reffs))./vec(Reffs)*100)")

    modeldata["postprocessing"] = FDataDict("file"=>"FV32-modes", "mode"=>1:10)
    modeldata=FinEtools.AlgoDeforLinearModule.exportmode(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["file"]*"1.vtk")`)

    true

end # NAFEMS_FV32_algo

function NAFEMS_TEST13H_vib()
    # Harmonic forced vibration problem is solved for a homogeneous square plate,
    # simply-supported on the circumference.
    # This is the TEST 13H from the Abaqus v 6.12 Benchmarks manual.
    # The test is recommended by the National Agency for Finite Element Methods and Standards (U.K.):
    # Test 13 from NAFEMS “Selected Benchmarks for Forced Vibration,” R0016, March 1993.
    #
    #
    # The plate is discretized with hexahedral solid elements. The simple support
    # condition is approximated by distributed rollers on the boundary.
    # Because only the out of plane displacements are prevented, the structure
    # has three rigid body modes in the plane of the plate.
    #
    #
    # The nonzero benchmark frequencies are (in hertz): 2.377, 5.961, 5.961,
    # 9.483, 12.133, 12.133, 15.468, 15.468 [Hz].

    println("""
    Homogeneous square plate, simply-supported on the circumference from
    the test 13 from NAFEMS “Selected Benchmarks for Forced Vibration,” R0016, March 1993.
    The nonzero benchmark frequencies are (in hertz): 2.377, 5.961, 5.961,
    9.483, 12.133, 12.133, 15.468, 15.468 [Hz].
    """)

    # t0 = time()

    E = 200*phun("GPa");# Young's modulus
    nu = 0.3;# Poisson ratio
    rho = 8000*phun("KG*M^-3");# mass density
    L = 10.0*phun("M"); # side of the square plate
    t = 0.05*phun("M"); # thickness of the square plate
    nL = 8; nt = 4;
    tolerance = t/nt/100;
    neigvs = 11;
    OmegaShift = (2*pi*0.5) ^ 2; # to resolve rigid body modes

    MR = DeforModelRed3D
    fens,fes  = H8block(L, L, t, nL, nL, nt)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
    nl = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    nl = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    nl = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    nl = selectnode(fens, box=[-Inf Inf L L -Inf Inf], inflate=tolerance)
    setebc!(u, nl, true, 3)
    applyebc!(u)
    numberdofs!(u)
    println("nfreedofs = $(u.nfreedofs)")

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearMSH8(MR, IntegData(fes, GaussRule(3,2)), material)
    femm = associategeometry!(femm, geom)
    K = stiffness(femm, geom, u)
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(3,3)), material)
    M = mass(femm, geom, u)

    if true
        t0 = time()
        d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
        d = d - OmegaShift;
        fs = real(sqrt.(complex(d)))/(2*pi)
        println("Reference Eigenvalues: $fs [Hz]")
        println("eigs solution ($(time() - t0) sec)")
    end
    true

end # NAFEMS_TEST13H_vib

function allrun()
    println("#####################################################")
    println("# NAFEMS_FV32_algo ")
    NAFEMS_FV32_algo()
    println("#####################################################")
    println("# NAFEMS_TEST13H_vib ")
    NAFEMS_TEST13H_vib()
    return true
end # function allrun

end # module NAFEMS_examples


text:
#
# Aggregate combinators lifted from Julia functions.
#

function AggregateQuery(fn, argtype::Type, restype::Type, haszero::Bool, q::Query)
    @assert method_exists(fn, (argtype,))
    @assert fits(domain(output(q)), Domain(argtype))
    argq = RecordQuery(q)
    sig = AggregateSig(fn, argtype, restype, haszero)
    ity = Input(domain(output(argq)))
    oty = Output(restype) |> setoptional(!haszero && isoptional(output(q)))
    return argq >> Query(sig, ity, oty)
end

AggregateQuery(fn, argtype::Type, restype::Type, q::Query) =
    AggregateQuery(fn, argtype, restype, true, q)

AggregateQuery(fn, argtype::Type, haszero::Bool, q::Query) =
    let restype = Union{Base.return_types(fn, (Vector{argtype},))...}
        AggregateQuery(fn, argtype, restype, haszero, q)
    end

AggregateQuery(fn, argtype::Type, q::Query) =
    AggregateQuery(fn, argtype, true, q)

AggregateQuery(fn, haszero::Bool, q::Query) =
    let argtype = datatype(domain(output(q)))
        AggregateQuery(fn, argtype, haszero, q)
    end

AggregateQuery(fn, q::Query) =
    AggregateQuery(fn, true, q)

immutable AggregateSig <: AbstractPrimitive
    fn::Function
    argtype::Type
    restype::Type
    haszero::Bool
end

ev(sig::AggregateSig, dv::DataVector) =
    let col = column(dv, 1)
        sig.haszero || !isoptional(col) ?
            plain_aggregate_impl(sig.fn, sig.restype, length(dv), col) :
            aggregate_impl(sig.fn, sig.restype, length(dv), col)
    end

function plain_aggregate_impl{T}(fn::Function, otype::Type{T}, len::Int, arg::Column)
    vals = Vector{T}(len)
    cr = cursor(arg)
    while !done(arg, cr)
        next!(arg, cr)
        vals[cr.pos] = fn(cr)
    end
    return PlainColumn(vals)
end

function aggregate_impl{T}(fn::Function, otype::Type{T}, len::Int, arg::Column)
    argoffs = offsets(arg)
    size = 0
    @inbounds for k = 1:len
        if argoffs[k] < argoffs[k+1]
            size += 1
        end
    end
    if size == len
        col = plain_aggregate_impl(fn, otype, len, arg)
        return OptionalColumn(col.offs, col.vals)
    end
    offs = Vector{Int}(len+1)
    offs[1] = 1
    vals = Vector{T}(size)
    cr = cursor(arg)
    n = 1
    while !done(arg, cr)
        next!(arg, cr)
        if !isempty(cr)
            vals[n] = fn(cr)
            n += 1
        end
        offs[cr.pos+1] = n
    end
    return Column(offs, vals)
end

ext: jl

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# logarithmic barrier function

export LogBarrier

"""
**Logarithmic barrier**

    LogBarrier(a=1, b=0, μ=1)

Returns the function
```math
f(x) = -μ⋅∑_i\\log(a⋅x_i+b),
```
for a nonnegative parameter `μ`.
"""
struct LogBarrier{R, S, T} <: ProximableFunction
    a::R
    b::S
    mu::T
    function LogBarrier{R, S, T}(a::R, b::S, mu::T) where {R, S, T}
        if mu <= 0
            error("parameter mu must be positive")
        else
            new(a, b, mu)
        end
    end
end

is_separable(f::LogBarrier) = true
is_convex(f::LogBarrier) = true

LogBarrier(a::R=1, b::S=0, mu::T=1) where {R, S, T} = LogBarrier{R, S, T}(a, b, mu)

function (f::LogBarrier)(x::AbstractArray{T}) where T <: Real
    sumf = T(0)
    v = T(0)
    for i in eachindex(x)
        v = f.a * x[i] + f.b
        if v <= T(0)
            return +Inf
        end
        sumf += log(v)
    end
    return -f.mu * sumf
end

function prox!(y::AbstractArray{T}, f::LogBarrier, x::AbstractArray{T}, gamma::Real=1) where T <: Real
    par = 4 * gamma * f.mu * f.a * f.a
    sumf = T(0)
    z = T(0)
    v = T(0)
    for i in eachindex(x)
        z = f.a * x[i] + f.b
        v = (z + sqrt(z * z + par)) / 2
        y[i] = (v - f.b) / f.a
        sumf += log(v)
    end
    return -f.mu * sumf
end

function prox!(y::AbstractArray{T}, f::LogBarrier, x::AbstractArray{T}, gamma::AbstractArray) where T <: Real
    par = 4 * f.mu * f.a * f.a
    sumf = T(0)
    z = T(0)
    v = T(0)
    for i in eachindex(x)
        par_i = gamma[i] * par
        z = f.a * x[i] + f.b
        v = (z + sqrt(z * z + par_i)) / 2
        y[i] = (v - f.b) / f.a
        sumf += log(v)
    end
    return -f.mu * sumf
end

function gradient!(y::AbstractArray{T}, f::LogBarrier, x::AbstractArray{T}) where T <: Real
    sum = T(0)
    for i in eachindex(x)
        logarg = f.a * x[i] + f.b
        y[i] = -f.mu * f.a / logarg
        sum += log(logarg)
    end
    sum *= -f.mu
    return sum
end

fun_name(f::LogBarrier) = "logarithmic barrier"
fun_dom(f::LogBarrier) = "AbstractArray{Real}"
fun_expr(f::LogBarrier) = "x ↦ -μ * sum( log(a*x_i+b), i=1,...,n )"
fun_params(f::LogBarrier) = "a = $(f.a), b = $(f.b), μ = $(f.mu)"

function prox_naive(f::LogBarrier, x::AbstractArray{T,1}, gamma::Union{Real, AbstractArray}=1) where T <: Real
    asqr = f.a * f.a
    z = f.a * x .+ f.b
    y = ((z .+ sqrt.(z .* z .+ 4 * gamma * f.mu * asqr)) / 2 .- f.b) / f.a
    fy = -f.mu * sum(log.(f.a .* y .+ f.b))
    return y, fy
end


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

"""
    Base.Broadcast

Module containing the broadcasting implementation.
"""
module Broadcast

using .Base.Cartesian
using .Base: Indices, OneTo, tail, to_shape, isoperator, promote_typejoin,
             _msk_end, unsafe_bitgetindex, bitcache_chunks, bitcache_size, dumpbitcache, unalias
import .Base: copy, copyto!, axes
export broadcast, broadcast!, BroadcastStyle, broadcast_axes, broadcastable, dotview, @__dot__, broadcast_preserving_zero_d

## Computing the result's axes: deprecated name
const broadcast_axes = axes

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

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

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

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

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

    Base.BroadcastStyle(::Type{<:MyContainer}) = Broadcast.Style{MyContainer}()
"""
struct Style{T} <: BroadcastStyle end

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    BroadcastStyle(::Primary, ::Secondary) = Primary()

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

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

### Lazy-wrapper for broadcasting

# `Broadcasted` wrap the arguments to `broadcast(f, args...)`. A statement like
#    y = x .* (x .+ 1)
# will result in code that is essentially
#    y = copy(Broadcasted(*, x, Broadcasted(+, x, 1)))
# `broadcast!` results in `copyto!(dest, Broadcasted(...))`.

# The use of `Nothing` in place of a `BroadcastStyle` has a different
# application, in the fallback method
#    copyto!(dest, bc::Broadcasted) = copyto!(dest, convert(Broadcasted{Nothing}, bc))
# This allows methods
#    copyto!(dest::DestType,  bc::Broadcasted{Nothing})
# that specialize on `DestType` to be easily disambiguated from
# methods that instead specialize on `BroadcastStyle`,
#    copyto!(dest::AbstractArray, bc::Broadcasted{MyStyle})

struct Broadcasted{Style<:Union{Nothing,BroadcastStyle}, Axes, F, Args<:Tuple} <: Base.AbstractBroadcasted
    f::F
    args::Args
    axes::Axes          # the axes of the resulting object (may be bigger than implied by `args` if this is nested inside a larger `Broadcasted`)
end

Broadcasted(f::F, args::Args, axes=nothing) where {F, Args<:Tuple} =
    Broadcasted{typeof(combine_styles(args...))}(f, args, axes)
function Broadcasted{Style}(f::F, args::Args, axes=nothing) where {Style, F, Args<:Tuple}
    # using Core.Typeof rather than F preserves inferrability when f is a type
    Broadcasted{Style, typeof(axes), Core.Typeof(f), Args}(f, args, axes)
end

Base.convert(::Type{Broadcasted{NewStyle}}, bc::Broadcasted{Style,Axes,F,Args}) where {NewStyle,Style,Axes,F,Args} =
    Broadcasted{NewStyle,Axes,F,Args}(bc.f, bc.args, bc.axes)

function Base.show(io::IO, bc::Broadcasted{Style}) where {Style}
    print(io, Broadcasted)
    # Only show the style parameter if we have a set of axes — representing an instantiated
    # "outermost" Broadcasted. The styles of nested Broadcasteds represent an intermediate
    # computation that is not relevant for dispatch, confusing, and just extra line noise.
    bc.axes isa Tuple && print(io, '{', Style, '}')
    print(io, '(', bc.f, ", ", bc.args, ')')
    nothing
end

## Allocating the output container
Base.similar(bc::Broadcasted, ::Type{T}) where {T} = similar(bc, T, axes(bc))
Base.similar(::Broadcasted{DefaultArrayStyle{N}}, ::Type{ElType}, dims) where {N,ElType} =
    similar(Array{ElType}, dims)
Base.similar(::Broadcasted{DefaultArrayStyle{N}}, ::Type{Bool}, dims) where N =
    similar(BitArray, dims)
# In cases of conflict we fall back on Array
Base.similar(::Broadcasted{ArrayConflict}, ::Type{ElType}, dims) where ElType =
    similar(Array{ElType}, dims)
Base.similar(::Broadcasted{ArrayConflict}, ::Type{Bool}, dims) =
    similar(BitArray, dims)

@inline Base.axes(bc::Broadcasted) = _axes(bc, bc.axes)
_axes(::Broadcasted, axes::Tuple) = axes
@inline _axes(bc::Broadcasted, ::Nothing)  = combine_axes(bc.args...)
_axes(bc::Broadcasted{<:AbstractArrayStyle{0}}, ::Nothing) = ()

@inline Base.axes(bc::Broadcasted{<:Any, <:NTuple{N}}, d::Integer) where N =
    d <= N ? axes(bc)[d] : OneTo(1)

BroadcastStyle(::Type{<:Broadcasted{Style}}) where {Style} = Style()
BroadcastStyle(::Type{<:Broadcasted{S}}) where {S<:Union{Nothing,Unknown}} =
    throw(ArgumentError("Broadcasted{Unknown} wrappers do not have a style assigned"))

argtype(::Type{Broadcasted{Style,Axes,F,Args}}) where {Style,Axes,F,Args} = Args
argtype(bc::Broadcasted) = argtype(typeof(bc))

@inline Base.eachindex(bc::Broadcasted) = _eachindex(axes(bc))
_eachindex(t::Tuple{Any}) = t[1]
_eachindex(t::Tuple) = CartesianIndices(t)

Base.IndexStyle(bc::Broadcasted) = IndexStyle(typeof(bc))
Base.IndexStyle(::Type{<:Broadcasted{<:Any,<:Tuple{Any}}}) = IndexLinear()
Base.IndexStyle(::Type{<:Broadcasted{<:Any}}) = IndexCartesian()

Base.LinearIndices(bc::Broadcasted{<:Any,<:Tuple{Any}}) = axes(bc)[1]

Base.ndims(::Broadcasted{<:Any,<:NTuple{N,Any}}) where {N} = N
Base.ndims(::Type{<:Broadcasted{<:Any,<:NTuple{N,Any}}}) where {N} = N

Base.size(bc::Broadcasted) = map(length, axes(bc))
Base.length(bc::Broadcasted) = prod(size(bc))

function Base.iterate(bc::Broadcasted)
    iter = eachindex(bc)
    iterate(bc, (iter,))
end
Base.@propagate_inbounds function Base.iterate(bc::Broadcasted, s)
    y = iterate(s...)
    y === nothing && return nothing
    i, newstate = y
    return (bc[i], (s[1], newstate))
end

Base.IteratorSize(::Type{<:Broadcasted{<:Any,<:NTuple{N,Base.OneTo}}}) where {N} = Base.HasShape{N}()
Base.IteratorEltype(::Type{<:Broadcasted}) = Base.EltypeUnknown()

## Instantiation fills in the "missing" fields in Broadcasted.
instantiate(x) = x

"""
    Broadcast.instantiate(bc::Broadcasted)

Construct and check the axes for the lazy Broadcasted object `bc`.

Custom [`BroadcastStyle`](@ref)s may override this default in cases where it is fast and easy
to compute and verify the resulting `axes` on-demand, leaving the `axis` field
of the `Broadcasted` object empty (populated with [`nothing`](@ref)).
"""
@inline function instantiate(bc::Broadcasted{Style}) where {Style}
    if bc.axes isa Nothing # Not done via dispatch to make it easier to extend instantiate(::Broadcasted{Style})
        axes = combine_axes(bc.args...)
    else
        axes = bc.axes
        check_broadcast_axes(axes, bc.args...)
    end
    return Broadcasted{Style}(bc.f, bc.args, axes)
end
instantiate(bc::Broadcasted{<:AbstractArrayStyle{0}}) = bc
# Tuples don't need axes, but when they have axes (for .= assignment), we need to check them (#33020)
instantiate(bc::Broadcasted{Style{Tuple}, Nothing}) = bc
function instantiate(bc::Broadcasted{Style{Tuple}})
    check_broadcast_axes(bc.axes, bc.args...)
    return bc
end
## Flattening

"""
    bcf = flatten(bc)

Create a "flat" representation of a lazy-broadcast operation.
From
   f.(a, g.(b, c), d)
we produce the equivalent of
   h.(a, b, c, d)
where
   h(w, x, y, z) = f(w, g(x, y), z)
In terms of its internal representation,
   Broadcasted(f, a, Broadcasted(g, b, c), d)
becomes
   Broadcasted(h, a, b, c, d)

This is an optional operation that may make custom implementation of broadcasting easier in
some cases.
"""
function flatten(bc::Broadcasted{Style}) where {Style}
    isflat(bc) && return bc
    # concatenate the nested arguments into {a, b, c, d}
    args = cat_nested(bc)
    # build a function `makeargs` that takes a "flat" argument list and
    # and creates the appropriate input arguments for `f`, e.g.,
    #          makeargs = (w, x, y, z) -> (w, g(x, y), z)
    #
    # `makeargs` is built recursively and looks a bit like this:
    #     makeargs(w, x, y, z) = (w, makeargs1(x, y, z)...)
    #                          = (w, g(x, y), makeargs2(z)...)
    #                          = (w, g(x, y), z)
    let makeargs = make_makeargs(()->(), bc.args), f = bc.f
        newf = @inline function(args::Vararg{Any,N}) where N
            f(makeargs(args...)...)
        end
        return Broadcasted{Style}(newf, args, bc.axes)
    end
end

const NestedTuple = Tuple{<:Broadcasted,Vararg{Any}}
isflat(bc::Broadcasted) = _isflat(bc.args)
_isflat(args::NestedTuple) = false
_isflat(args::Tuple) = _isflat(tail(args))
_isflat(args::Tuple{}) = true

cat_nested(t::Broadcasted, rest...) = (cat_nested(t.args...)..., cat_nested(rest...)...)
cat_nested(t::Any, rest...) = (t, cat_nested(rest...)...)
cat_nested() = ()

"""
    make_makeargs(makeargs_tail::Function, t::Tuple) -> Function

Each element of `t` is one (consecutive) node in a broadcast tree.
Ignoring `makeargs_tail` for the moment, the job of `make_makeargs` is
to return a function that takes in flattened argument list and returns a
tuple (each entry corresponding to an entry in `t`, having evaluated
the corresponding element in the broadcast tree). As an additional
complication, the passed in tuple may be longer than the number of leaves
in the subtree described by `t`. The `makeargs_tail` function should
be called on such additional arguments (but not the arguments consumed
by `t`).
"""
@inline make_makeargs(makeargs_tail, t::Tuple{}) = makeargs_tail
@inline function make_makeargs(makeargs_tail, t::Tuple)
    makeargs = make_makeargs(makeargs_tail, tail(t))
    (head, tail...)->(head, makeargs(tail...)...)
end
function make_makeargs(makeargs_tail, t::Tuple{<:Broadcasted, Vararg{Any}})
    bc = t[1]
    # c.f. the same expression in the function on leaf nodes above. Here
    # we recurse into siblings in the broadcast tree.
    let makeargs_tail = make_makeargs(makeargs_tail, tail(t)),
            # Here we recurse into children. It would be valid to pass in makeargs_tail
            # here, and not use it below. However, in that case, our recursion is no
            # longer purely structural because we're building up one argument (the closure)
            # while destructuing another.
            makeargs_head = make_makeargs((args...)->args, bc.args),
            f = bc.f
        # Create two functions, one that splits of the first length(bc.args)
        # elements from the tuple and one that yields the remaining arguments.
        # N.B. We can't call headargs on `args...` directly because
        # args is flattened (i.e. our children have not been evaluated
        # yet).
        headargs, tailargs = make_headargs(bc.args), make_tailargs(bc.args)
        return @inline function(args::Vararg{Any,N}) where N
            args1 = makeargs_head(args...)
            a, b = headargs(args1...), makeargs_tail(tailargs(args1...)...)
            (f(a...), b...)
        end
    end
end

@inline function make_headargs(t::Tuple)
    let headargs = make_headargs(tail(t))
        return @inline function(head, tail::Vararg{Any,N}) where N
            (head, headargs(tail...)...)
        end
    end
end
@inline function make_headargs(::Tuple{})
    return @inline function(tail::Vararg{Any,N}) where N
        ()
    end
end

@inline function make_tailargs(t::Tuple)
    let tailargs = make_tailargs(tail(t))
        return @inline function(head, tail::Vararg{Any,N}) where N
            tailargs(tail...)
        end
    end
end
@inline function make_tailargs(::Tuple{})
    return @inline function(tail::Vararg{Any,N}) where N
        tail
    end
end

## Broadcasting utilities ##

## logic for deciding the BroadcastStyle

"""
    combine_styles(cs...) -> BroadcastStyle

Decides which `BroadcastStyle` to use for any number of value arguments.
Uses [`BroadcastStyle`](@ref) to get the style for each argument, and uses
[`result_style`](@ref) to combine styles.

# Examples

```jldoctest
julia> Broadcast.combine_styles([1], [1 2; 3 4])
Base.Broadcast.DefaultArrayStyle{2}()
```
"""
function combine_styles end

combine_styles() = DefaultArrayStyle{0}()
combine_styles(c) = result_style(BroadcastStyle(typeof(c)))
combine_styles(c1, c2) = result_style(combine_styles(c1), combine_styles(c2))
@inline combine_styles(c1, c2, cs...) = result_style(combine_styles(c1), combine_styles(c2, cs...))

"""
    result_style(s1::BroadcastStyle[, s2::BroadcastStyle]) -> BroadcastStyle

Takes one or two `BroadcastStyle`s and combines them using [`BroadcastStyle`](@ref) to
determine a common `BroadcastStyle`.

# Examples

```jldoctest
julia> Broadcast.result_style(Broadcast.DefaultArrayStyle{0}(), Broadcast.DefaultArrayStyle{3}())
Base.Broadcast.DefaultArrayStyle{3}()

julia> Broadcast.result_style(Broadcast.Unknown(), Broadcast.DefaultArrayStyle{1}())
Base.Broadcast.DefaultArrayStyle{1}()
```
"""
function result_style end

result_style(s::BroadcastStyle) = s
result_style(s1::S, s2::S) where S<:BroadcastStyle = S()
# Test both orders so users typically only have to declare one order
result_style(s1, s2) = result_join(s1, s2, BroadcastStyle(s1, s2), BroadcastStyle(s2, s1))

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

# Indices utilities

"""
    combine_axes(As...) -> Tuple

Determine the result axes for broadcasting across all values in `As`.

```jldoctest
julia> Broadcast.combine_axes([1], [1 2; 3 4; 5 6])
(Base.OneTo(3), Base.OneTo(2))

julia> Broadcast.combine_axes(1, 1, 1)
()
```
"""
@inline combine_axes(A, B...) = broadcast_shape(axes(A), combine_axes(B...))
@inline combine_axes(A, B) = broadcast_shape(axes(A), axes(B))
combine_axes(A) = axes(A)

# shape (i.e., tuple-of-indices) inputs
broadcast_shape(shape::Tuple) = shape
broadcast_shape(shape::Tuple, shape1::Tuple, shapes::Tuple...) = broadcast_shape(_bcs(shape, shape1), shapes...)
# _bcs consolidates two shapes into a single output shape
_bcs(::Tuple{}, ::Tuple{}) = ()
_bcs(::Tuple{}, newshape::Tuple) = (newshape[1], _bcs((), tail(newshape))...)
_bcs(shape::Tuple, ::Tuple{}) = (shape[1], _bcs(tail(shape), ())...)
function _bcs(shape::Tuple, newshape::Tuple)
    return (_bcs1(shape[1], newshape[1]), _bcs(tail(shape), tail(newshape))...)
end
# _bcs1 handles the logic for a single dimension
_bcs1(a::Integer, b::Integer) = a == 1 ? b : (b == 1 ? a : (a == b ? a : throw(DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths $a and $b"))))
_bcs1(a::Integer, b) = a == 1 ? b : (first(b) == 1 && last(b) == a ? b : throw(DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths $a and $(length(b))")))
_bcs1(a, b::Integer) = _bcs1(b, a)
_bcs1(a, b) = _bcsm(b, a) ? axistype(b, a) : (_bcsm(a, b) ? axistype(a, b) : throw(DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths $(length(a)) and $(length(b))")))
# _bcsm tests whether the second index is consistent with the first
_bcsm(a, b) = a == b || length(b) == 1
_bcsm(a, b::Number) = b == 1
_bcsm(a::Number, b::Number) = a == b || b == 1
# Ensure inferrability when dealing with axes of different AbstractUnitRange types
# (We may not want to define general promotion rules between, say, OneTo and Slice, but if
#  we get here we know the axes are at least consistent for the purposes of broadcasting)
axistype(a::T, b::T) where T = a
axistype(a::OneTo, b::OneTo) = OneTo{Int}(a)
axistype(a, b) = UnitRange{Int}(a)

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

## Indexing manipulations
"""
    newindex(argument, I)
    newindex(I, keep, default)

Recompute index `I` such that it appropriately constrains broadcasted dimensions to the source.

Two methods are supported, both allowing for `I` to be specified as either a [`CartesianIndex`](@ref) or
an `Int`.

* `newindex(argument, I)` dynamically constrains `I` based upon the axes of `argument`.
* `newindex(I, keep, default)` constrains `I` using the pre-computed tuples `keeps` and `defaults`.
    * `keep` is a tuple of `Bool`s, where `keep[d] == true` means that dimension `d` in `I` should be preserved as is
    * `default` is a tuple of Integers, specifying what index to use in dimension `d` when `keep[d] == false`.
    Any remaining indices in `I` beyond the length of the `keep` tuple are truncated. The `keep` and `default`
    tuples may be created by `newindexer(argument)`.
"""
Base.@propagate_inbounds newindex(arg, I::CartesianIndex) = CartesianIndex(_newindex(axes(arg), I.I))
Base.@propagate_inbounds newindex(arg, I::Integer) = CartesianIndex(_newindex(axes(arg), (I,)))
Base.@propagate_inbounds _newindex(ax::Tuple, I::Tuple) = (ifelse(Base.unsafe_length(ax[1])==1, ax[1][1], I[1]), _newindex(tail(ax), tail(I))...)
Base.@propagate_inbounds _newindex(ax::Tuple{}, I::Tuple) = ()
Base.@propagate_inbounds _newindex(ax::Tuple, I::Tuple{}) = (ax[1][1], _newindex(tail(ax), ())...)
Base.@propagate_inbounds _newindex(ax::Tuple{}, I::Tuple{}) = ()

# If dot-broadcasting were already defined, this would be `ifelse.(keep, I, Idefault)`.
@inline newindex(I::CartesianIndex, keep, Idefault) = CartesianIndex(_newindex(I.I, keep, Idefault))
@inline newindex(i::Integer, keep::Tuple{Bool}, idefault) = ifelse(keep[1], i, idefault[1])
@inline newindex(i::Integer, keep::Tuple{}, idefault) = CartesianIndex(())
@inline _newindex(I, keep, Idefault) =
    (ifelse(keep[1], I[1], Idefault[1]), _newindex(tail(I), tail(keep), tail(Idefault))...)
@inline _newindex(I, keep::Tuple{}, Idefault) = ()  # truncate if keep is shorter than I

# newindexer(A) generates `keep` and `Idefault` (for use by `newindex` above)
# for a particular array `A`; `shapeindexer` does so for its axes.
@inline newindexer(A) = shapeindexer(axes(A))
@inline shapeindexer(ax) = _newindexer(ax)
@inline _newindexer(indsA::Tuple{}) = (), ()
@inline function _newindexer(indsA::Tuple)
    ind1 = indsA[1]
    keep, Idefault = _newindexer(tail(indsA))
    (Base.length(ind1)::Integer != 1, keep...), (first(ind1), Idefault...)
end

@inline function Base.getindex(bc::Broadcasted, I::Union{Integer,CartesianIndex})
    @boundscheck checkbounds(bc, I)
    @inbounds _broadcast_getindex(bc, I)
end
Base.@propagate_inbounds Base.getindex(
    bc::Broadcasted,
    i1::Union{Integer,CartesianIndex},
    i2::Union{Integer,CartesianIndex},
    I::Union{Integer,CartesianIndex}...,
) =
    bc[CartesianIndex((i1, i2, I...))]
Base.@propagate_inbounds Base.getindex(bc::Broadcasted) = bc[CartesianIndex(())]

@inline Base.checkbounds(bc::Broadcasted, I::Union{Integer,CartesianIndex}) =
    Base.checkbounds_indices(Bool, axes(bc), (I,)) || Base.throw_boundserror(bc, (I,))

"""
    _broadcast_getindex(A, I)

Index into `A` with `I`, collapsing broadcasted indices to their singleton indices as appropriate.
"""
Base.@propagate_inbounds _broadcast_getindex(A::Union{Ref,AbstractArray{<:Any,0},Number}, I) = A[] # Scalar-likes can just ignore all indices
Base.@propagate_inbounds _broadcast_getindex(::Ref{Type{T}}, I) where {T} = T
# Tuples are statically known to be singleton or vector-like
Base.@propagate_inbounds _broadcast_getindex(A::Tuple{Any}, I) = A[1]
Base.@propagate_inbounds _broadcast_getindex(A::Tuple, I) = A[I[1]]
# Everything else falls back to dynamically dropping broadcasted indices based upon its axes
Base.@propagate_inbounds _broadcast_getindex(A, I) = A[newindex(A, I)]

# In some cases, it's more efficient to sort out which dimensions should be dropped
# ahead of time (often when the size checks aren't able to be lifted out of the loop).
# The Extruded struct computes that information ahead of time and stores it as a pair
# of tuples to optimize indexing later. This is most commonly needed for `Array` and
# other `AbstractArray` subtypes that wrap `Array` and dynamically ask it for its size.
struct Extruded{T, K, D}
    x::T
    keeps::K    # A tuple of booleans, specifying which indices should be passed normally
    defaults::D # A tuple of integers, specifying the index to use when keeps[i] is false (as defaults[i])
end
@inline axes(b::Extruded) = axes(b.x)
Base.@propagate_inbounds _broadcast_getindex(b::Extruded, i) = b.x[newindex(i, b.keeps, b.defaults)]
extrude(x::AbstractArray) = Extruded(x, newindexer(x)...)
extrude(x) = x

# For Broadcasted
Base.@propagate_inbounds function _broadcast_getindex(bc::Broadcasted{<:Any,<:Any,<:Any,<:Any}, I)
    args = _getindex(bc.args, I)
    return _broadcast_getindex_evalf(bc.f, args...)
end
# Hack around losing Type{T} information in the final args tuple. Julia actually
# knows (in `code_typed`) the _value_ of these types, statically displaying them,
# but inference is currently skipping inferring the type of the types as they are
# transiently placed in a tuple as the argument list is lispily constructed. These
# additional methods recover type stability when a `Type` appears in one of the
# first two arguments of a function.
Base.@propagate_inbounds function _broadcast_getindex(bc::Broadcasted{<:Any,<:Any,<:Any,<:Tuple{Ref{Type{T}},Vararg{Any}}}, I) where {T}
    args = _getindex(tail(bc.args), I)
    return _broadcast_getindex_evalf(bc.f, T, args...)
end
Base.@propagate_inbounds function _broadcast_getindex(bc::Broadcasted{<:Any,<:Any,<:Any,<:Tuple{Any,Ref{Type{T}},Vararg{Any}}}, I) where {T}
    arg1 = _broadcast_getindex(bc.args[1], I)
    args = _getindex(tail(tail(bc.args)), I)
    return _broadcast_getindex_evalf(bc.f, arg1, T, args...)
end
Base.@propagate_inbounds function _broadcast_getindex(bc::Broadcasted{<:Any,<:Any,<:Any,<:Tuple{Ref{Type{T}},Ref{Type{S}},Vararg{Any}}}, I) where {T,S}
    args = _getindex(tail(tail(bc.args)), I)
    return _broadcast_getindex_evalf(bc.f, T, S, args...)
end

# Utilities for _broadcast_getindex
Base.@propagate_inbounds _getindex(args::Tuple, I) = (_broadcast_getindex(args[1], I), _getindex(tail(args), I)...)
Base.@propagate_inbounds _getindex(args::Tuple{Any}, I) = (_broadcast_getindex(args[1], I),)
Base.@propagate_inbounds _getindex(args::Tuple{}, I) = ()

@inline _broadcast_getindex_evalf(f::Tf, args::Vararg{Any,N}) where {Tf,N} = f(args...)  # not propagate_inbounds

"""
    Broadcast.broadcastable(x)

Return either `x` or an object like `x` such that it supports [`axes`](@ref), indexing, and its type supports [`ndims`](@ref).

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

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

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

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

julia> Broadcast.broadcastable("hello") # Strings break convention of matching iteration and act like a scalar instead
Base.RefValue{String}("hello")
```
"""
broadcastable(x::Union{Symbol,AbstractString,Function,UndefInitializer,Nothing,RoundingMode,Missing,Val,Ptr,Regex,Pair}) = Ref(x)
broadcastable(::Type{T}) where {T} = Ref{Type{T}}(T)
broadcastable(x::Union{AbstractArray,Number,Ref,Tuple,Broadcasted}) = x
# Default to collecting iterables — which will error for non-iterables
broadcastable(x) = collect(x)
broadcastable(::Union{AbstractDict, NamedTuple}) = throw(ArgumentError("broadcasting over dictionaries and `NamedTuple`s is reserved"))

## Computation of inferred result type, for empty and concretely inferred cases only
_broadcast_getindex_eltype(bc::Broadcasted) = Base._return_type(bc.f, eltypes(bc.args))
_broadcast_getindex_eltype(A) = eltype(A)  # Tuple, Array, etc.

eltypes(::Tuple{}) = Tuple{}
eltypes(t::Tuple{Any}) = Tuple{_broadcast_getindex_eltype(t[1])}
eltypes(t::Tuple{Any,Any}) = Tuple{_broadcast_getindex_eltype(t[1]), _broadcast_getindex_eltype(t[2])}
eltypes(t::Tuple) = Tuple{_broadcast_getindex_eltype(t[1]), eltypes(tail(t)).types...}

# Inferred eltype of result of broadcast(f, args...)
combine_eltypes(f, args::Tuple) = Base._return_type(f, eltypes(args))

## Broadcasting core

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

Broadcast the function `f` over the arrays, tuples, collections, [`Ref`](@ref)s and/or scalars `As`.

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

The resulting container type is established by the following rules:

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

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

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

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

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

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

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

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

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

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

```
"""
broadcast(f::Tf, As...) where {Tf} = materialize(broadcasted(f, As...))

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

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

Like [`broadcast`](@ref), but store the result of
`broadcast(f, As...)` in the `dest` array.
Note that `dest` is only used to store the result, and does not supply
arguments to `f` unless it is also listed in the `As`,
as in `broadcast!(f, A, A, B)` to perform `A[:] = broadcast(f, A, B)`.

# Examples
```jldoctest
julia> A = [1.0; 0.0]; B = [0.0; 0.0];

julia> broadcast!(+, B, A, (0, -2.0));

julia> B
2-element Vector{Float64}:
  1.0
 -2.0

julia> A
2-element Vector{Float64}:
 1.0
 0.0

julia> broadcast!(+, A, A, (0, -2.0));

julia> A
2-element Vector{Float64}:
  1.0
 -2.0
```
"""
broadcast!(f::Tf, dest, As::Vararg{Any,N}) where {Tf,N} = (materialize!(dest, broadcasted(f, As...)); dest)

"""
    broadcast_preserving_zero_d(f, As...)

Like [`broadcast`](@ref), except in the case of a 0-dimensional result where it returns a 0-dimensional container

Broadcast automatically unwraps zero-dimensional results to be just the element itself,
but in some cases it is necessary to always return a container — even in the 0-dimensional case.
"""
@inline function broadcast_preserving_zero_d(f, As...)
    bc = broadcasted(f, As...)
    r = materialize(bc)
    return length(axes(bc)) == 0 ? fill!(similar(bc, typeof(r)), r) : r
end
@inline broadcast_preserving_zero_d(f) = fill(f())
@inline broadcast_preserving_zero_d(f, as::Number...) = fill(f(as...))

"""
    Broadcast.materialize(bc)

Take a lazy `Broadcasted` object and compute the result
"""
@inline materialize(bc::Broadcasted) = copy(instantiate(bc))
materialize(x) = x

@inline function materialize!(dest, x)
    return materialize!(dest, instantiate(Broadcasted(identity, (x,), axes(dest))))
end

@inline function materialize!(dest, bc::Broadcasted{Style}) where {Style}
    return materialize!(combine_styles(dest, bc), dest, bc)
end
@inline function materialize!(::BroadcastStyle, dest, bc::Broadcasted{Style}) where {Style}
    return copyto!(dest, instantiate(Broadcasted{Style}(bc.f, bc.args, axes(dest))))
end

## general `copy` methods
@inline copy(bc::Broadcasted{<:AbstractArrayStyle{0}}) = bc[CartesianIndex()]
copy(bc::Broadcasted{<:Union{Nothing,Unknown}}) =
    throw(ArgumentError("broadcasting requires an assigned BroadcastStyle"))

const NonleafHandlingStyles = Union{DefaultArrayStyle,ArrayConflict}

@inline function copy(bc::Broadcasted{Style}) where {Style}
    ElType = combine_eltypes(bc.f, bc.args)
    if Base.isconcretetype(ElType)
        # We can trust it and defer to the simpler `copyto!`
        return copyto!(similar(bc, ElType), bc)
    end
    # When ElType is not concrete, use narrowing. Use the first output
    # value to determine the starting output eltype; copyto_nonleaf!
    # will widen `dest` as needed to accommodate later values.
    bc′ = preprocess(nothing, bc)
    iter = eachindex(bc′)
    y = iterate(iter)
    if y === nothing
        # if empty, take the ElType at face value
        return similar(bc′, ElType)
    end
    # Initialize using the first value
    I, state = y
    @inbounds val = bc′[I]
    dest = similar(bc′, typeof(val))
    @inbounds dest[I] = val
    # Now handle the remaining values
    return copyto_nonleaf!(dest, bc′, iter, state, 1)
end

## general `copyto!` methods
# The most general method falls back to a method that replaces Style->Nothing
# This permits specialization on typeof(dest) without introducing ambiguities
@inline copyto!(dest::AbstractArray, bc::Broadcasted) = copyto!(dest, convert(Broadcasted{Nothing}, bc))

# Performance optimization for the common identity scalar case: dest .= val
@inline function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractArrayStyle{0}})
    # Typically, we must independently execute bc for every storage location in `dest`, but:
    # IF we're in the common no-op identity case with no nested args (like `dest .= val`),
    if bc.f === identity && bc.args isa Tuple{Any} && isflat(bc)
        # THEN we can just extract the argument and `fill!` the destination with it
        return fill!(dest, bc.args[1][])
    else
        # Otherwise, fall back to the default implementation like above
        return copyto!(dest, convert(Broadcasted{Nothing}, bc))
    end
end

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

# Preprocessing a `Broadcasted` does two things:
# * unaliases any arguments from `dest`
# * "extrudes" the arguments where it is advantageous to pre-compute the broadcasted indices
@inline preprocess(dest, bc::Broadcasted{Style}) where {Style} = Broadcasted{Style}(bc.f, preprocess_args(dest, bc.args), bc.axes)
preprocess(dest, x) = extrude(broadcast_unalias(dest, x))

@inline preprocess_args(dest, args::Tuple) = (preprocess(dest, args[1]), preprocess_args(dest, tail(args))...)
preprocess_args(dest, args::Tuple{Any}) = (preprocess(dest, args[1]),)
preprocess_args(dest, args::Tuple{}) = ()

# Specialize this method if all you want to do is specialize on typeof(dest)
@inline function copyto!(dest::AbstractArray, bc::Broadcasted{Nothing})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && bc.args isa Tuple{AbstractArray} # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    bc′ = preprocess(dest, bc)
    @simd for I in eachindex(bc′)
        @inbounds dest[I] = bc′[I]
    end
    return dest
end

# Performance optimization: for BitArray outputs, we cache the result
# in a "small" Vector{Bool}, and then copy in chunks into the output
@inline function copyto!(dest::BitArray, bc::Broadcasted{Nothing})
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    ischunkedbroadcast(dest, bc) && return chunkedcopyto!(dest, bc)
    length(dest) < 256 && return invoke(copyto!, Tuple{AbstractArray, Broadcasted{Nothing}}, dest, bc)
    tmp = Vector{Bool}(undef, bitcache_size)
    destc = dest.chunks
    cind = 1
    bc′ = preprocess(dest, bc)
    for P in Iterators.partition(eachindex(bc′), bitcache_size)
        ind = 1
        @simd for I in P
            @inbounds tmp[ind] = bc′[I]
            ind += 1
        end
        @simd for i in ind:bitcache_size
            @inbounds tmp[i] = false
        end
        dumpbitcache(destc, cind, tmp)
        cind += bitcache_chunks
    end
    return dest
end

# For some BitArray operations, we can work at the level of chunks. The trivial
# implementation just walks over the UInt64 chunks in a linear fashion.
# This requires three things:
#   1. The function must be known to work at the level of chunks (or can be converted to do so)
#   2. The only arrays involved must be BitArrays or scalar Bools
#   3. There must not be any broadcasting beyond scalar — all array sizes must match
# We could eventually allow for all broadcasting and other array types, but that
# requires very careful consideration of all the edge effects.
const ChunkableOp = Union{typeof(&), typeof(|), typeof(xor), typeof(~), typeof(identity),
    typeof(!), typeof(*), typeof(==)} # these are convertible to chunkable ops by liftfuncs
const BroadcastedChunkableOp{Style<:Union{Nothing,BroadcastStyle}, Axes, F<:ChunkableOp, Args<:Tuple} = Broadcasted{Style,Axes,F,Args}
ischunkedbroadcast(R, bc::BroadcastedChunkableOp) = ischunkedbroadcast(R, bc.args)
ischunkedbroadcast(R, args) = false
ischunkedbroadcast(R, args::Tuple{<:BitArray,Vararg{Any}}) = size(R) == size(args[1]) && ischunkedbroadcast(R, tail(args))
ischunkedbroadcast(R, args::Tuple{<:Bool,Vararg{Any}}) = ischunkedbroadcast(R, tail(args))
ischunkedbroadcast(R, args::Tuple{<:BroadcastedChunkableOp,Vararg{Any}}) = ischunkedbroadcast(R, args[1]) && ischunkedbroadcast(R, tail(args))
ischunkedbroadcast(R, args::Tuple{}) = true

# Convert compatible functions to chunkable ones. They must also be green-lighted as ChunkableOps
liftfuncs(bc::Broadcasted{Style}) where {Style} = Broadcasted{Style}(bc.f, map(liftfuncs, bc.args), bc.axes)
liftfuncs(bc::Broadcasted{Style,<:Any,typeof(sign)}) where {Style} = Broadcasted{Style}(identity, map(liftfuncs, bc.args), bc.axes)
liftfuncs(bc::Broadcasted{Style,<:Any,typeof(!)}) where {Style} = Broadcasted{Style}(~, map(liftfuncs, bc.args), bc.axes)
liftfuncs(bc::Broadcasted{Style,<:Any,typeof(*)}) where {Style} = Broadcasted{Style}(&, map(liftfuncs, bc.args), bc.axes)
liftfuncs(bc::Broadcasted{Style,<:Any,typeof(==)}) where {Style} = Broadcasted{Style}((~)∘(xor), map(liftfuncs, bc.args), bc.axes)
liftfuncs(x) = x

liftchunks(::Tuple{}) = ()
liftchunks(args::Tuple{<:BitArray,Vararg{Any}}) = (args[1].chunks, liftchunks(tail(args))...)
# Transform scalars to repeated scalars the size of a chunk
liftchunks(args::Tuple{<:Bool,Vararg{Any}}) = (ifelse(args[1], typemax(UInt64), UInt64(0)), liftchunks(tail(args))...)
ithchunk(i) = ()
Base.@propagate_inbounds ithchunk(i, c::Vector{UInt64}, args...) = (c[i], ithchunk(i, args...)...)
Base.@propagate_inbounds ithchunk(i, b::UInt64, args...) = (b, ithchunk(i, args...)...)
@inline function chunkedcopyto!(dest::BitArray, bc::Broadcasted)
    isempty(dest) && return dest
    f = flatten(liftfuncs(bc))
    args = liftchunks(f.args)
    dc = dest.chunks
    @simd for i in eachindex(dc)
        @inbounds dc[i] = f.f(ithchunk(i, args...)...)
    end
    @inbounds dc[end] &= Base._msk_end(dest)
    return dest
end

@noinline throwdm(axdest, axsrc) =
    throw(DimensionMismatch("destination axes $axdest are not compatible with source axes $axsrc"))

function restart_copyto_nonleaf!(newdest, dest, bc, val, I, iter, state, count)
    # Function barrier that makes the copying to newdest type stable
    for II in Iterators.take(iter, count)
        newdest[II] = dest[II]
    end
    newdest[I] = val
    return copyto_nonleaf!(newdest, bc, iter, state, count+1)
end

function copyto_nonleaf!(dest, bc::Broadcasted, iter, state, count)
    T = eltype(dest)
    while true
        y = iterate(iter, state)
        y === nothing && break
        I, state = y
        @inbounds val = bc[I]
        if val isa T || typeof(val) === T
            @inbounds dest[I] = val
        else
            # This element type doesn't fit in dest. Allocate a new dest with wider eltype,
            # copy over old values, and continue
            newdest = Base.similar(dest, promote_typejoin(T, typeof(val)))
            return restart_copyto_nonleaf!(newdest, dest, bc, val, I, iter, state, count)
        end
        count += 1
    end
    return dest
end

## Tuple methods

@inline function copy(bc::Broadcasted{Style{Tuple}})
    dim = axes(bc)
    length(dim) == 1 || throw(DimensionMismatch("tuple only supports one dimension"))
    N = length(dim[1])
    return ntuple(k -> @inbounds(_broadcast_getindex(bc, k)), Val(N))
end

## scalar-range broadcast operations ##
# DefaultArrayStyle and \ are not available at the time of range.jl
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::OrdinalRange) = r
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::StepRangeLen) = r
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::LinRange) = r

broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::OrdinalRange) = range(-first(r), step=-step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepRangeLen) = StepRangeLen(-r.ref, -r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::LinRange) = LinRange(-r.start, -r.stop, length(r))

broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Real, r::AbstractUnitRange) = range(x + first(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::AbstractUnitRange, x::Real) = range(first(r) + x, length=length(r))
# For #18336 we need to prevent promotion of the step type:
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::AbstractRange, x::Number) = range(first(r) + x, step=step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::AbstractRange) = range(x + first(r), step=step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::StepRangeLen{T}, x::Number) where T =
    StepRangeLen{typeof(T(r.ref)+x)}(r.ref + x, r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::StepRangeLen{T}) where T =
    StepRangeLen{typeof(x+T(r.ref))}(x + r.ref, r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::LinRange, x::Number) = LinRange(r.start + x, r.stop + x, length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::LinRange) = LinRange(x + r.start, x + r.stop, length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r1::AbstractRange, r2::AbstractRange) = r1 + r2

broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::AbstractUnitRange, x::Number) = range(first(r)-x, length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::AbstractRange, x::Number) = range(first(r)-x, step=step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::AbstractRange) = range(x-first(r), step=-step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepRangeLen{T}, x::Number) where T =
    StepRangeLen{typeof(T(r.ref)-x)}(r.ref - x, r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::StepRangeLen{T}) where T =
    StepRangeLen{typeof(x-T(r.ref))}(x - r.ref, -r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::LinRange, x::Number) = LinRange(r.start - x, r.stop - x, length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::LinRange) = LinRange(x - r.start, x - r.stop, length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r1::AbstractRange, r2::AbstractRange) = r1 - r2

broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::AbstractRange) = range(x*first(r), step=x*step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::StepRangeLen{T}) where {T} =
    StepRangeLen{typeof(x*T(r.ref))}(x*r.ref, x*r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::LinRange) = LinRange(x * r.start, x * r.stop, r.len)
# separate in case of noncommutative multiplication
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::AbstractRange, x::Number) = range(first(r)*x, step=step(r)*x, length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::StepRangeLen{T}, x::Number) where {T} =
    StepRangeLen{typeof(T(r.ref)*x)}(r.ref*x, r.step*x, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::LinRange, x::Number) = LinRange(r.start * x, r.stop * x, r.len)

broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::AbstractRange, x::Number) = range(first(r)/x, step=step(r)/x, length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::StepRangeLen{T}, x::Number) where {T} =
    StepRangeLen{typeof(T(r.ref)/x)}(r.ref/x, r.step/x, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::LinRange, x::Number) = LinRange(r.start / x, r.stop / x, r.len)

broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::AbstractRange) = range(x\first(r), step=x\step(r), length=length(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::StepRangeLen) = StepRangeLen(x\r.ref, x\r.step, length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::LinRange) = LinRange(x \ r.start, x \ r.stop, r.len)

broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::UnitRange) = big(r.start):big(last(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::StepRange) = big(r.start):big(r.step):big(last(r))
broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::StepRangeLen) = StepRangeLen(big(r.ref), big(r.step), length(r), r.offset)
broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::LinRange) = LinRange(big(r.start), big(r.stop), length(r))

## CartesianIndices
broadcasted(::typeof(+), I::CartesianIndices{N}, j::CartesianIndex{N}) where N =
    CartesianIndices(map((rng, offset)->rng .+ offset, I.indices, Tuple(j)))
broadcasted(::typeof(+), j::CartesianIndex{N}, I::CartesianIndices{N}) where N =
    I .+ j
broadcasted(::typeof(-), I::CartesianIndices{N}, j::CartesianIndex{N}) where N =
    CartesianIndices(map((rng, offset)->rng .- offset, I.indices, Tuple(j)))
function broadcasted(::typeof(-), j::CartesianIndex{N}, I::CartesianIndices{N}) where N
    diffrange(offset, rng) = range(offset-last(rng), length=length(rng))
    Iterators.reverse(CartesianIndices(map(diffrange, Tuple(j), I.indices)))
end

## In specific instances, we can broadcast masked BitArrays whole chunks at a time
# Very intentionally do not support much functionality here: scalar indexing would be O(n)
struct BitMaskedBitArray{N,M}
    parent::BitArray{N}
    mask::BitArray{M}
    BitMaskedBitArray{N,M}(parent, mask) where {N,M} = new(parent, mask)
end
@inline function BitMaskedBitArray(parent::BitArray{N}, mask::BitArray{M}) where {N,M}
    @boundscheck checkbounds(parent, mask)
    BitMaskedBitArray{N,M}(parent, mask)
end
Base.@propagate_inbounds dotview(B::BitArray, i::BitArray) = BitMaskedBitArray(B, i)
Base.show(io::IO, B::BitMaskedBitArray) = foreach(arg->show(io, arg), (typeof(B), (B.parent, B.mask)))
# Override materialize! to prevent the BitMaskedBitArray from escaping to an overrideable method
@inline materialize!(B::BitMaskedBitArray, bc::Broadcasted{<:Any,<:Any,typeof(identity),Tuple{Bool}}) = fill!(B, bc.args[1])
@inline materialize!(B::BitMaskedBitArray, bc::Broadcasted{<:Any}) = materialize!(SubArray(B.parent, to_indices(B.parent, (B.mask,))), bc)
function Base.fill!(B::BitMaskedBitArray, b::Bool)
    Bc = B.parent.chunks
    Ic = B.mask.chunks
    @inbounds if b
        for i = 1:length(Bc)
            Bc[i] |= Ic[i]
        end
    else
        for i = 1:length(Bc)
            Bc[i] &= ~Ic[i]
        end
    end
    return B
end

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

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

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

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

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

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

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

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

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

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

@inline broadcasted_kwsyntax(f, args...; kwargs...) = broadcasted((args...)->f(args...; kwargs...), args...)
@inline function broadcasted(f, args...)
    args′ = map(broadcastable, args)
    broadcasted(combine_styles(args′...), f, args′...)
end
# Due to the current Type{T}/DataType specialization heuristics within Tuples,
# the totally generic varargs broadcasted(f, args...) method above loses Type{T}s in
# mapping broadcastable across the args. These additional methods with explicit
# arguments ensure we preserve Type{T}s in the first or second argument position.
@inline function broadcasted(f, arg1, args...)
    arg1′ = broadcastable(arg1)
    args′ = map(broadcastable, args)
    broadcasted(combine_styles(arg1′, args′...), f, arg1′, args′...)
end
@inline function broadcasted(f, arg1, arg2, args...)
    arg1′ = broadcastable(arg1)
    arg2′ = broadcastable(arg2)
    args′ = map(broadcastable, args)
    broadcasted(combine_styles(arg1′, arg2′, args′...), f, arg1′, arg2′, args′...)
end
@inline broadcasted(::S, f, args...) where S<:BroadcastStyle = Broadcasted{S}(f, args)

end # module


text:
export device, send_to_host, send_to_device

using ElasticArrays
using Flux
using CUDA
using Adapt
using Random

import CUDA: device

send_to_host(x) = send_to_device(Val(:cpu), x)
send_to_device(::Val{:cpu}, x) = x  # cpu(x) is not very efficient! So by default we do nothing here.

send_to_device(::Val{:cpu}, x::CuArray) = adapt(Array, x)
send_to_device(::Val{:gpu}, x) = Flux.fmap(a -> adapt(CuArray{Float32}, a), x)

const KnownArrayVariants = Union{
    SubArray{<:Any,<:Any,<:Union{ReservoirArrayBuffer,CircularArrayBuffer,ElasticArray}},
    Base.ReshapedArray{
        <:Any,
        <:Any,
        <:SubArray{
            <:Any,
            <:Any,
            <:Union{ReservoirArrayBuffer,CircularArrayBuffer,ElasticArray},
        },
    },
    Base.ReshapedArray{
        <:Any,
        <:Any,
        <:Union{ReservoirArrayBuffer,CircularArrayBuffer,ElasticArray},
    },
    SubArray{
        <:Any,
        <:Any,
        <:Base.ReshapedArray{
            <:Any,
            <:Any,
            <:SubArray{
                <:Any,
                <:Any,
                <:Union{ReservoirArrayBuffer,CircularArrayBuffer,ElasticArray},
            },
        },
    },
}

# https://github.com/JuliaReinforcementLearning/ReinforcementLearningCore.jl/issues/130
send_to_device(::Val{:cpu}, x::KnownArrayVariants) = Array(x)
send_to_device(::Val{:gpu}, x::Union{KnownArrayVariants,ElasticArray}) = CuArray(x)

"""
    device(model)

Detect the suitable running device for the `model`.
Return `Val(:cpu)` by default.
"""
device(x) = device(Flux.trainable(x))
device(x::Function) = nothing
device(::CuArray) = Val(:gpu)
device(::Array) = Val(:cpu)
device(x::Tuple{}) = nothing
device(x::NamedTuple{(),Tuple{}}) = nothing
device(x::ElasticArray) = device(x.data)

function device(x::Random.AbstractRNG)
    if x isa CUDA.CURAND.RNG
        Val(:gpu)
    else
        Val(:cpu)
    end
end

function device(x::Union{Tuple,NamedTuple})
    d1 = device(x[1])
    if isnothing(d1)
        device(Base.tail(x))
    else
        d1
    end
end

# recoganize Torch.jl
# device(x::Tensor) = Val(Symbol(:gpu, x.device))


text:
using Documenter, SurrogateModelOptim

makedocs(;
    modules=[SurrogateModelOptim],
    format = Documenter.HTML(   assets = ["assets/favicon.ico"],
                                prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
                "Home" => "index.md",
                "Manual" => [
                        "smoptimize.md",
                        "surrogate_model.md",
                        "options.md"
                ],
        ],
    repo="https://github.com/MrUrq/SurrogateModelOptim.jl/blob/{commit}{path}#L{line}",
    sitename="SurrogateModelOptim.jl",
    #authors="Magnus Urquhart",
    #strict = true,
    #clean = true,
    #checkdocs = :none,
)

deploydocs(;
    repo="github.com/MrUrq/SurrogateModelOptim.jl",
    target="build",
    deps=nothing,
    make=nothing,
)


text:
function jldopen(@nospecialize(f::Function), args...; kws...)
    jld = jldopen(args...; kws...)
    try
        return f(jld)
    finally
        close(jld)
    end
end

function loadtodict!(d::Dict, g::Union{JLDFile,Group}, prefix::String="")
    for k in keys(g)
        v = g[k]
        if v isa Group
            loadtodict!(d, v, prefix*k*"/")
        else
            d[prefix*k] = v
        end
    end
    return d
end

# Name used in JLD2 file to identify objects stored with `save_object`
const SINGLE_OBJECT_NAME = "single_stored_object"

"""
    save_object(filename, x)

Stores an object `x` in a new JLD2 file at `filename`. If a file exists at this
path, it will be overwritten.

Since the JLD2 format requires that all objects have a name, the object will be
stored as `single_sotred_object`. If you want to store more than one object, use
[`@save`](@ref) macro, [`jldopen`](@ref) or the FileIO API.

# Example

To save the string `hello` to the JLD2 file example.jld2:

    hello = "world"
    save_object("example.jld2", hello)
"""
function save_object(filename, x)
  jldopen(filename, "w") do file
    file[SINGLE_OBJECT_NAME] = x
  end
  return
end

"""
    load_object(filename)

Returns the only available object from the JLD2 file `filename` (The stored
object name is inconsequential). If the file contains more than one or no
objects, the function throws an `ArgumentError`.

For loading more than one object, use [`@load`](@ref) macro, [`jldopen`](@ref)
or the FileIO API.

# Example

To load the only object from the JLD2 file example.jld2:

    hello = "world"
    save_object("example.jld2", hello)
    hello_loaded = load_object("example.jld2")
"""
function load_object(filename)
  jldopen(filename, "r") do file
    all_keys = keys(file)
    length(all_keys) == 0 && throw(ArgumentError("File $filename does not contain any object"))
    length(all_keys) > 1 && throw(ArgumentError("File $filename contains more than one object. Use `load` or `@load` instead"))
    file[all_keys[1]] #Uses HDF5 functionality of treating the file like a dict
  end
end

"""
    jldsave(filename, compress=false; kwargs...)

Creates a JLD2 file at `filename` and stores the variables given as keyword arguments.

# Examples

    jldsave("example.jld2"; a=1, b=2, c)

is equivalent to

    jldopen("example.jld2, "w") do f
        f["a"] = 1
        f["b"] = 2
        f["c"] = c
    end

To choose the io type `IOStream` instead of the default `MmapIO` use
`jldsave{IOStream}(fn; kwargs...)`.
"""
function jldsave(filename::AbstractString, compress=false, iotype::T=MmapIO;
                    kwargs...
                    ) where T<:Union{Type{IOStream},Type{MmapIO}}
    jldopen(filename, "w"; compress=compress, iotype=iotype) do f
        wsession = JLDWriteSession()
        for (k,v) in pairs(kwargs)
            write(f, string(k), v, wsession)
        end
    end
end

jldsave(filename::AbstractString, iotype::Union{Type{IOStream},Type{MmapIO}}; kwargs...) =
    jldsave(filename, false, iotype; kwargs...)


text:
using Metaheuristics
using Test

@testset "Constrained" begin

    function test_result(fitness, tol)
        @test ≈(fitness, 0.0, atol=tol)
    end

    function run_methods_cop(problem)
        desired_accuracy = 1e-4

        ff, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)

        # number of function evaluations
        f_calls = 0

        f(x) = begin
            f_calls += 1
            ff(x)
        end

        information = Information(f_optimum = 0.0)
        options = Options(f_tol = desired_accuracy, h_tol=1e-5, seed = 2, debug=false, store_convergence=true)

        methods = [
                   ECA(ε=0.5, options = options, information = information),
                   #DE(options = options, information = information),
                  ]

        for method in methods
            f_calls = 0
            res = optimize(f, bounds, method)
            show(IOBuffer(), res)
            show(IOBuffer(), "text/plain", res.population)
            show(IOBuffer(), "text/html", res.population)
            show(IOBuffer(), res.population[1])
            show(IOBuffer(), method)
            fitness = minimum( res )
            test_result(fitness, desired_accuracy)
            @test f_calls == res.f_calls
            @test fvals(res) == fvals(res.population)
            @test sum(abs.(hval(res.best_sol))) < 1e-3
            @test !any(gval(res.best_sol) .> 0)
            @test Metaheuristics.diff_check(res, method.information, method.options) isa Bool
            @test Metaheuristics.diversity_stop_check(res, method.information, method.options) isa Bool
        end
    end

    for problem in [:constrained1]
        run_methods_cop(problem)
    end
end


text:
include("GetSegParaFromCSV.jl")
#=
ReadCFGPARA.Readin_cfgParas(machine友好名0=machine友好名0 )
rows = ReadCFGPARA.Paras_tofloatArray( metersPerSeg = Float64(6.001) ,timeUnit=0.426)
paraDictObj = Dict()
chn1SegIDRng=[1,8000]
chn2SegIDRng=[9000,9000+8000]
ReadCFGPARA.initParaDict(paraDict=paraDictObj ,rowsObj=rows  ,chn1SegIDRng=chn1SegIDRng , chn2SegIDRng=chn2SegIDRng)
=#
function TestMain(;machine友好名0="GYYS001" )
    ReadCFGPARA.Readin_cfgParas(machine友好名0=machine友好名0 )

    rows = ReadCFGPARA.Paras_tofloatArray( metersPerSeg = Float64(6.001),timeUnit=0.426 )
    # show(ReadCFGPARA.CsvInfoObj.csvInfoDefencePara )
    rowNum,fieldNum   = size( rows )

    for row =1:rowNum
        for fieldidx = 1:fieldNum
            print( rows[row,fieldidx],"," )
        end
        println()
    end
    paraDictObj = Dict()
    chn1SegIDRng=[1,8000]
    chn2SegIDRng=[9000,9000+8000]
    ReadCFGPARA.initParaDict(paraDict=paraDictObj ,rowsObj=rows  ,chn1SegIDRng=chn1SegIDRng , chn2SegIDRng=chn2SegIDRng)
    for key = chn1SegIDRng[1] :chn1SegIDRng[2] # 9000:17000)
        try
            valobj = paraDictObj[key]
            print(key,": ")
            for idx =1: length(valobj )
                print( valobj[idx] ,",  ")
            end
            println()
        catch( e )
            # println(e)
        end
    end

    for key in   chn2SegIDRng[1] :chn2SegIDRng[2]# 9000:17000)
        try
            valobj = paraDictObj[key]
            print(key,": ")
            for idx =1: length(valobj )
                print( valobj[idx] ,",  ")
            end
            println()
        catch( e )
            # println(e)
        end
    end

end

TestMain(machine友好名0="GYYS001" )

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
push!(LOAD_PATH,"../src/")

using Documenter, DTWDT
using DTWDT.DTWDTfunctions

makedocs(
    modules = [DTWDT],
    format = Documenter.HTML(),
    sitename="DTWDT.jl",
    authors="kurama",
    pages = [
        "Home" => "index.md",
        "Functions" => "Functions.md",
        "Sample results" => "Sampleresults.md"

    ],
)

deploydocs(
    repo="github.com/kura-okubo/DTWDT.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#     Wrapper to IOCTL calls.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@inline _ioctl(io::IOStream, args...) = _ioctl(fd(io), args...)

function _ioctl(fd::Integer, request::Integer, arg::T) where T <: Number
    ret = ccall(:ioctl, Cint, (Cint, Culong, T), fd, request, arg)
    ret < 0 && throw(SystemError("Error in IOCTL call", Libc.errno()))
    return ret
end

function _ioctl(fd::Integer, request::Integer, arg::Base.RefValue{T}) where T
    ret = ccall(:ioctl, Cint, (Cint, Culong, Ref{T}), fd, request, arg)
    ret < 0 && throw(SystemError("Error in IOCTL call", Libc.errno()))
    return ret
end

function _ioctl(fd::Integer, request::Integer, arg::AbstractVector{T}) where T
    ret = ccall(:ioctl, Cint, (Cint, Culong, Ref{T}), fd, request, arg)
    ret < 0 && throw(SystemError("Error in IOCTL call", Libc.errno()))
    return ret
end


text:
function ql(A::Matrix)
    Q,R=qr(A[end:-1:1,end:-1:1])
    Q[end:-1:1,end:-1:1],R[end:-1:1,end:-1:1]
end

joukowsky(z) = 0.5*(z+1/z)

function jacobimatrix(a,b,t0,t1,N)
    J = BandedMatrix(Float64,N,N,1,1)

    for i = 1:min(length(a),N)
        J[i,i] = a[i]
    end
    for i=length(a)+1:N
        J[i,i] = t0
    end

    for i = 1:min(length(b),N-1)
        J[i,i+1] = J[i+1,i] = b[i]
    end

    for i=length(b)+1:N-1
        J[i,i+1] = J[i+1,i] = t1
    end

    J
end

jacobimatrix(a,b,N) = jacobimatrix(a,b,0,.5,N)

function jacobioperator(a,b,t0,t1)
    n = max(length(a),length(b)+1)
    a = [a;zeros(n-length(a))]; b = [b;.5+zeros(n-length(b))]
    SymTriPertToeplitz(ToeplitzOperator([t1],[t0,t1]),SymTriOperator(a-t0,b-t1))
end

jacobioperator(a,b) = jacobioperator(a,b,0,.5)

freejacobioperator() = SymTriPertToeplitz([0.],[.5],0.,.5)

## tridiagonal ql

function tridql!(L::Matrix)
    n=size(L,1)

  # Now we do QL for the compact part in the top left
    cc=Array{eltype(L)}(n)
    ss=Array{eltype(L)}(n-1)

    for i = n:-1:2
        nrm=sqrt(L[i-1,i]^2+L[i,i]^2)
        c,s = L[i,i]/nrm, -L[i-1,i]/nrm
        if i > 2
            L[i-1:i,i-2:i] = [c s; -s c]*L[i-1:i,i-2:i]
            L[i-1,i]=0
        else
            L[i-1:i,i-1:i] = [c s; -s c]*L[i-1:i,i-1:i]
            L[i-1,i]=0
        end
        cc[i]=c
        ss[i-1]=s
    end
    cc[1]=sign(L[1,1])
    L[1,1]=abs(L[1,1])
    cc,ss,L
end

function tridql!(J::BandedMatrix)
    n=size(J,1)
    L=BandedMatrix(copy(J.data),J.m,2,0)

  # Now we do QL for the compact part in the top left
    cc=Array{eltype(J)}(n)
    ss=Array{eltype(J)}(n-1)

    for i = n:-1:2
        nrm=sqrt(J[i-1,i]^2+J[i,i]^2)
        c,s = J[i,i]/nrm, -J[i-1,i]/nrm

        for j=max(i-2,1):i
            L[i,j]=-s*J[i-1,j]+c*J[i,j]
            J[i-1,j]=c*J[i-1,j]+s*J[i,j]
        end
        cc[i]=c
        ss[i-1]=s
    end
    cc[1]=sign(J[1,1])
    L[1,1]=abs(J[1,1])
    cc,ss,L
end

#Finds NxN truncation of C such that C'(Q_k(s)) =  (P_k(s)),
# where P_k has Jacobi coeffs a,b and Q_k has Jacobi coeffs c,d
function connectioncoeffsmatrix(a::AbstractVector, b::AbstractVector, c::AbstractVector, d::AbstractVector, N)
  if N>max(length(a),length(b)+1,length(c),length(d)+1)
    a = [a;zeros(N-length(a))]; b = [b;.5+zeros(N-length(b))]
    c = [c;zeros(N-length(c))]; d = [d;.5+zeros(N-length(d))]
  end

  C = zeros(eltype(a),N,N)
  C[1,1] = 1
  C[1,2] = (c[1]-a[1])/b[1]
  C[2,2] = d[1]/b[1]
  for j = 3:N
    C[1,j] = ((c[1]-a[j-1])*C[1,j-1] + d[1]*C[2,j-1] - b[j-2]*C[1,j-2])/b[j-1]
    for i = 2:j-1
      C[i,j] = (d[i-1]*C[i-1,j-1] + (c[i]-a[j-1])*C[i,j-1] + d[i]*C[i+1,j-1] - b[j-2]*C[i,j-2])/b[j-1]
    end
    C[j,j] = d[j-1]*C[j-1,j-1]/b[j-1]
  end
  C
end

#Makes the matrix C which transforms the coefficients of an expansion in "right" matrix orthonormal polynomial with recurrence coeffs a,b (J) to those of an expansion in the "right" MOPs with recurrence coeffs c,d (D). So C is block-upper-triangular and CJ = DC (modulo the final column). The types of a,b,c,d must be BlockArrays which are like column vectors of length N whose entries are kxk blocks. The blocks of a and c should be symmetric.
function connectioncoeffsmatrix(a::AbstractBlockArray, b::AbstractBlockArray, c::AbstractBlockArray, d::AbstractBlockArray)
  k = blocksize(a,1,1)[1]
  N = nblocks(a,1)
  C = BlockArray(zeros(k*N,k*N),k*ones(Int64,N),k*ones(Int64,N))
  C[Block(1,1)] = eye(k,k)
  C[Block(1,2)] = (c[Block(1,1)]-a[Block(1,1)])*pinv(b[Block(1,1)]')
  C[Block(2,2)] = d[Block(1,1)]'/b[Block(1,1)]'
  for j = 3:N
    C[Block(1,j)] = (c[Block(1,1)]*C[Block(1,j-1)]-C[Block(1,j-1)]*a[Block(j-1,1)] + d[Block(1,1)]*C[Block(2,j-1)]-C[Block(1,j-2)]*b[Block(j-2,1)])/b[Block(j-1,1)]'
    for i = 2:j-1
        C[Block(i,j)] = (d[Block(i-1,1)]'*C[Block(i-1,j-1)] + c[Block(i,1)]*C[Block(i,j-1)]-C[Block(i,j-1)]*a[Block(j-1,1)] +d[Block(i,1)]*C[Block(i+1,j-1)] - C[Block(i,j-2)]*b[Block(j-2,1)])/b[Block(j-1,1)]'
    end
    C[Block(j,j)] = (d[Block(j-1,1)]'*C[Block(j-1,j-1)])/b[Block(j-1,1)]'
  end
  C
end

# This is for Chebyshev U
connectioncoeffsmatrix(a,b,N) = connectioncoeffsmatrix(a,b,[],[],N)

# Converts coefficients a^J to coefficients a^D using Clenshaw
function applyconversion(J::SymTriPertToeplitz,D::SymTriPertToeplitz,v::Vector)
  N = length(v)
  T = eltype(b)
  b = zeros(T,N); b1 = zeros(T,N); b2 = zeros(T,N)
  for k = N:-1:1
    # before: b = b_k+1, b1 = b_k+2, (and b2 = b_k+3 is to be forgotten)
    b2 = pad((D-J[k,k]*I)*b,N)/J[k,k+1]-b1*(J[k,k+1]/J[k+1,k+2])
    b2[1] += v[k]
    b2, b1, b = b1, b, b2
    # after: b = b_k, b1 = b_k+1, b2 = b_k+2
  end
  b
end


text:
include("./src/ShellModel.jl")
using .ShellModel

##########################
### solve EC for approx. eigenpairs
##########################
"""
calc_moment=true   # calc. mu&Q-moment
write_appwav0false # to write .appwav you need to make sample w.f.
"""
function solve(write_appwav=false)
    sntpath = "snts/random_input/validsnts_sigma1/sdshl/"
    Hs = [  sntpath*"tmp_$i"*".snt" for i = 0:99 ]
    #tpath = "/Users/Sota/Desktop/develop_ShellModel.jl/"
    tpath = "/Users/sotauu/Desktop/develop_ShellModel.jl/"
    elogf = "runlog_sigma1_valid.dat"

    num_ev_target = 1
    verbose = false # to show all Energy(EC)
    tJNs = [ [tJ,num_ev_target] for tJ in 0:2:8]
    for target_nuc in ["Mg24","Al26","Si28"]
        solveEC(Hs,target_nuc,tJNs;wpath=tpath,verbose=verbose,exact_logf=elogf)
    end

    tJNs = [ [tJ,num_ev_target] for tJ in 1:2:9]
    for target_nuc in ["Mg25"]
        solveEC(Hs,target_nuc,tJNs;wpath=tpath,verbose=verbose,exact_logf=elogf)
    end
end

solve()


text:
"""
    apply_legend_spec!(...)

Internal helper functions to apply the legend element corresponding to a
specific drawing handle.
"""
function apply_legend_spec!(
            g::GS,
            h::DrawingHandle{<:Scatter2D},
            labels::Union{String, Vector{String}},
            figid::String
        )::Nothing

    scatter = h.drawing
    labels isa Vector || (labels = fill(labels, scatter.nobj))
    for k ∈ 1:scatter.nobj
        "\n\ttext \"$(labels[k])\"" |> g
        if scatter.linestyles[k].lstyle != -1
            # line plot
            "line" |> g; apply_linestyle!(g, scatter.linestyles[k]; nosmooth=true)
            mcol = false
            if isdef(scatter.markerstyles[k].color) &&
                    (scatter.markerstyles[k].color != scatter.linestyles[k].color)
                mcol = true
                add_sub_marker!(Figure(figid; _noreset=true), scatter.markerstyles[k])
            end
            # apply markerstyle
            apply_markerstyle!(g, scatter.markerstyles[k], mcol=mcol)
        else
            # scatter plot
            if !isdef(scatter.markerstyles[k].color)
                scatter.markerstyles[k].color = scatter.linestyles[k].color
            end
            apply_markerstyle!(g, scatter.markerstyles[k])
        end
    end
    return
end
function apply_legend_spec!(
            g::GS,
            h::DrawingHandle{<:Fill2D},
            label::String,
            ::String
        )::Nothing

    fill = h.drawing
    "\n\ttext \"$label\" fill $(col2str(fill.fillstyle.fill))" |> g
    return
end
function apply_legend_spec!(
            g::GS,
            h::DrawingHandle{<:Hist2D},
            label::String,
            ::String
        )::Nothing

    hist = h.drawing
    "\n\ttext \"$label\"" |> g
    # precedence of fill over color
    if hist.barstyle.fill != c"white"
        "fill $(col2str(hist.barstyle.fill))" |> g
    else
        "marker square color $(col2str(hist.barstyle.color))" |> g
    end
    return
end
function apply_legend_spec!(
            g::GS,
            h::DrawingHandle{<:Bar2D},
            labels::Union{String,Vector{String}},
            ::String
        )::Nothing

    bar = h.drawing
    labels isa Vector || (labels = fill(labels, bar.nobj))
    for (k, barstyle) ∈ enumerate(bar.barstyles)
        "\n\ttext \"$(labels[k])\"" |> g
        if barstyle.fill != c"white"
            "fill $(col2str(barstyle.fill))" |> g
        else
            "marker square color $(col2str(barstyle.color))" |> g
        end
    end
    return
end

"""
    apply_legend!(g, leg, parent_font)

Internal function to apply a `Legend` object `leg` in a GLE context with
entries `entries` (constructed through the `apply_drawings` process).
"""
function apply_legend!(
            g::GS,
            l::Legend,
            parent_font::String,
            figid::String
        )::Nothing

    l.off && return
    isempty(l.handles) && return
    # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    "\nbegin key"  |> g
    apply_textstyle!(g, l.textstyle, parent_font; addset=true)

    # global commands
    "\n\tcompact"      |> g
    l.nobox && "nobox" |> g
    isdef(l.bgcolor)   && "background $(col2str(l.bgcolor))"        |> g
    isdef(l.margins)   && "margins $(l.margins[1]) $(l.margins[2])" |> g
    sum(l.offset) > 0  && "offset $(l.offset[1]) $(l.offset[2])"    |> g
    isdef(l.position)  && "\n\tposition $(l.position)"              |> g

    for (handle, label) ∈ zip(l.handles, l.labels)
        apply_legend_spec!(g, handle, label, figid)
    end
    #
    "\nend key"    |> g
    # <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    return
end


text:
function solve_007()
    # init sieve vector
    numbers = Set(2:150000)
    primes = Set()

    # run sieve
    for candidate = 2:150000
        for num in numbers
            if num%candidate == 0 && num ≠ candidate delete!(numbers, num) end
        end
        if candidate % 100 == 0 println(candidate) end
    end

    return numbers
end


text:
""" Sometimes useful given `@confcall`'s tendency to grow methods endlessly.
Useful special case: `delete_all_methods(fn, types=Tuple{Val{whatev},Vararg{Any}})` """
delete_all_methods(fn, types=Tuple{Vararg{Any}}) =
    foreach(Base.delete_method, methods(fn, types).ms)


text:
# Autogenerated wrapper script for Dynare_preprocessor_jll for x86_64-w64-mingw32
export dynare_preprocessor

JLLWrappers.@generate_wrapper_header("Dynare_preprocessor")
JLLWrappers.@declare_executable_product(dynare_preprocessor)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        dynare_preprocessor,
        "bin\\dynare-preprocessor.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
```Write a function that takes a list of numbers, a starting base b1 and a target
base b2 and interprets the list as a number in base b1 and converts it into a number
in base b2 (in the form of a list-of-digits). So for example [2,1,0] in base 3 gets
converted to base 10 as [2,1].```

list = [2,1,0]

function base_3_10(x::Vector{Int64})
    num = zeros(Int64, size(x,1))
    for (i,j) in enumerate(x)
        num[i] += x[i] * abs(3 ^ i-1)
    end
    return sum(num)
end

println(" $(base_3_10(list)) base 10")

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

""" Coupling elements for JuliaFEM. """
module FEMCoupling

using Reexport
@reexport using FEMBase

include("get_C.jl")
include("distributed_coupling.jl")
export Coupling, add_reference_node!, add_coupling_nodes!

end


text:
using StructuralEquationModels, ProximalSEM, ProximalOperators, CSV, DataFrames, StenoGraphs

# load data
dat = DataFrame(CSV.File("data_dem.csv"))

############################################################################
### define models
############################################################################

observed_vars = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8]
latent_vars = [:ind60, :dem60, :dem65]

graph = @StenoGraph begin

    ind60 → fixed(1)*x1 + x2 + x3
    dem60 → fixed(1)*y1 + y2 + y3 + y4
    dem65 → fixed(1)*y5 + y6 + y7 + y8

    dem60 ← ind60
    dem65 ← dem60
    dem65 ← ind60

    _(observed_vars) ↔ _(observed_vars)
    _(latent_vars) ↔ _(latent_vars)

    y1 ↔ label(:cov_15)*y5
    y2 ↔ label(:cov_24)*y4 + label(:cov_26)*y6
    y3 ↔ label(:cov_37)*y7
    y4 ↔ label(:cov_48)*y8
    y6 ↔ label(:cov_68)*y8

end

partable = ParameterTable(;
    graph = graph,
    latent_vars = latent_vars,
    observed_vars = observed_vars)

ram_mat = RAMMatrices(partable)

# use lasso from ProximalSEM
ind = get_identifier_indices([:cov_15, :cov_24, :cov_26, :cov_37, :cov_48, :cov_68], ram_mat)
λ = zeros(31); λ[16:20] .= 0.002

model_prox = Sem(
    specification = partable,
    data = dat,
    loss = (SemML,),
    diff = SemDiffProximal,
    operator_g = NormL1(λ)
)

solution_prox = sem_fit(model_prox)

@testset "lasso_converged" begin
    @test solution_prox.solution.optimization_result[:iterations] < 1000
end


text:
using MultivariateStats
using Test
import Random

@testset "Factor Rotations" begin

    Random.seed!(37273)

    X = randn(5, 2)

    ## Varimax rotation
    # Comparison with R's `varimax` function called as
    # using Printf
    # Base.show(io::IO, f::Float64) = @printf(io, "%1.8f", f)
    # using RCall
    # R"vm <- varimax($X, normalize = FALSE, eps = 1e-12)"
    # R"F <- qm$loadings"
    # R"R <- qm$rotmat"
    # @rget F
    # @rget R

    F = [-2.56254778 -0.05712524;
         -0.37387139  0.53683094;
         -1.76868624  0.67513474;
          0.08057810  0.34826689;
         -0.06196338  0.20433233]
    R = [ 0.96256370  0.27105556;
         -0.27105556  0.96256370]

    FR = fit(FactorRotation, X, alg = Varimax())
    @test isapprox(FR.F, F, rtol = √eps())
    @test isapprox(FR.R, R, rtol = √eps())

    ## Quartimax rotation
    # Comparison with R's `GPArotation::quartimax` function called as
    # R"qm <- GPArotation::quartimax($X, eps = 1e-6, maxit = 10000)"
    # R"F <- qm$loadings"
    # R"R <- qm$Th"
    # @rget F
    # @rget R

    F = [-2.55665706 -0.18280898;
         -0.39976499  0.51783707;
         -1.79968636  0.58752613;
          0.06339043  0.35180152;
         -0.07191598  0.20104540]
    R = [ 0.94810241  0.31796513;
         -0.31796513  0.94810241]

    FR = fit(FactorRotation, X, alg = Quartimax())
    @test isapprox(FR.F, F, rtol = 1e-6)
    @test isapprox(FR.R, R, rtol = 1e-6)

    # ## Equamax
    # X = randn(10, 5)

    # # Comparison with Matlab's `rotatefactors` called as
    # # using MATLAB
    # # mat"rotatefactors($X, 'Method', 'equamax', 'Normalize', 'off')"

    # F = [ 0.20295895  2.07838653 -0.54122597 -0.12487534  0.98358184;
    #      -1.03815211  0.17239352 -1.27261240 -0.41485812  1.44036802;
    #      -1.32942778  0.75883602  0.53246754 -0.39153390  2.81316627;
    #       0.17943030  0.10200669 -0.02036786 -2.05727725  0.36179405;
    #      -0.57150075  0.83976222  0.13897091 -0.58537812  1.78231347;
    #      -0.24047953 -0.14609212 -0.10808767 -0.31099008  0.53277677;
    #       1.63020394 -0.61318538 -0.77879612 -0.44644772 -0.30331055;
    #       2.92004862  0.09873379 -0.22817883 -0.00487858 -1.05015088;
    #       0.36654224  0.21341754 -2.75662510 -0.02498979 -0.47774001;
    #      -0.28002665 -0.92746802 -0.32811579 -0.52183972  0.15056833]

    # # Equamax for n x p matrix means γ = p / 2
    # FR = fit(FactorRotation, X, alg = Equamax(size(X, 2)))
    # @test isapprox(FR.F, F, rtol = √eps())

    # ## Parsimax

    # # Comparison with Matlab's `rotatefactors` called as
    # # using MATLAB
    # # mat"rotatefactors($X, 'Method', 'parsimax', 'Normalize', 'off')"

    # F = [ 0.20088755  2.13833957 -0.51952937 -0.15258917  0.85486221;
    #      -1.04430793  0.26969187 -1.25045141 -0.47113722  1.42298177;
    #      -1.32771913  0.92728472  0.57264722 -0.48321345  2.74040511;
    #       0.17856831  0.12930037 -0.00621592 -2.06837501  0.28192375;
    #      -0.57130614  0.94788128  0.16710829 -0.64273775  1.70426783;
    #      -0.24126103 -0.11168423 -0.10096713 -0.33074538  0.53024057;
    #       1.62645119 -0.62482282 -0.79024454 -0.44334505 -0.27314675;
    #       2.91932555  0.03504886 -0.25094949  0.03018747 -1.05060055;
    #       0.35443644  0.20160255 -2.76135234 -0.02259016 -0.46464651;
    #      -0.28200304 -0.91296776 -0.32805047 -0.53293768  0.19126898]

    # # Equamax for n x p matrix means γ = p / 2
    # FR = fit(FactorRotation, X, alg = Parsimax(size(X)..., maxiter = 2000, ϵ = 1e-5))
    # @test isapprox(FR.F, F, rtol = √eps())

end


text:
using AbstractAlgebra

abstract type DiffRing <: AbstractAlgebra.Ring end
abstract type DiffField <: AbstractAlgebra.Field end
abstract type DiffRingElem <: AbstractAlgebra.RingElem end
abstract type DiffFieldElem <: DiffRingElem end

max_stars_repo_name: pogudingleb/DifferentialAlgebra

max_issues_repo_name: pogudingleb/DifferentialAlgebra

max_forks_repo_name: pogudingleb/DifferentialAlgebra


text:
using ParameterSchedulers
using Test

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