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

# Matrix-matrix multiplication

AdjOrTransStridedMat{T} = Union{Adjoint{T, <:StridedMatrix}, Transpose{T, <:StridedMatrix}}
StridedMaybeAdjOrTransMat{T} = Union{StridedMatrix{T}, Adjoint{T, <:StridedMatrix}, Transpose{T, <:StridedMatrix}}

# matmul.jl: Everything to do with dense matrix multiplication

matprod(x, y) = x*y + x*y

# dot products

dot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasReal} = BLAS.dot(x, y)
dot(x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where {T<:BlasComplex} = BLAS.dotc(x, y)

function dot(x::Vector{T}, rx::AbstractRange{TI}, y::Vector{T}, ry::AbstractRange{TI}) where {T<:BlasReal,TI<:Integer}
    if length(rx) != length(ry)
        throw(DimensionMismatch(lazy"length of rx, $(length(rx)), does not equal length of ry, $(length(ry))"))
    end
    if minimum(rx) < 1 || maximum(rx) > length(x)
        throw(BoundsError(x, rx))
    end
    if minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError(y, ry))
    end
    GC.@preserve x y BLAS.dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end

function dot(x::Vector{T}, rx::AbstractRange{TI}, y::Vector{T}, ry::AbstractRange{TI}) where {T<:BlasComplex,TI<:Integer}
    if length(rx) != length(ry)
        throw(DimensionMismatch(lazy"length of rx, $(length(rx)), does not equal length of ry, $(length(ry))"))
    end
    if minimum(rx) < 1 || maximum(rx) > length(x)
        throw(BoundsError(x, rx))
    end
    if minimum(ry) < 1 || maximum(ry) > length(y)
        throw(BoundsError(y, ry))
    end
    GC.@preserve x y BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end

function *(transx::Transpose{<:Any,<:StridedVector{T}}, y::StridedVector{T}) where {T<:BlasComplex}
    x = transx.parent
    return BLAS.dotu(x, y)
end

# Matrix-vector multiplication
function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S<:Real}
    TS = promote_op(matprod, T, S)
    y = isconcretetype(TS) ? convert(AbstractVector{TS}, x) : x
    mul!(similar(x, TS, size(A,1)), A, y)
end
function (*)(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(similar(x,TS,axes(A,1)),A,x)
end

# these will throw a DimensionMismatch unless B has 1 row (or 1 col for transposed case):
(*)(a::AbstractVector, tB::Transpose{<:Any,<:AbstractMatrix}) = reshape(a, length(a), 1) * tB
(*)(a::AbstractVector, adjB::Adjoint{<:Any,<:AbstractMatrix}) = reshape(a, length(a), 1) * adjB
(*)(a::AbstractVector, B::AbstractMatrix) = reshape(a, length(a), 1) * B

@inline mul!(y::StridedVector{T}, A::StridedVecOrMat{T}, x::StridedVector{T},
             alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemv!(y, 'N', A, x, alpha, beta)

# Complex matrix times real vector. Reinterpret the matrix as a real matrix and do real matvec compuation.
for elty in (Float32, Float64)
    @eval begin
        @inline function mul!(y::StridedVector{Complex{$elty}}, A::StridedVecOrMat{Complex{$elty}}, x::StridedVector{$elty},
                alpha::Real, beta::Real)
            Afl = reinterpret($elty, A)
            yfl = reinterpret($elty, y)
            mul!(yfl, Afl, x, alpha, beta)
            return y
        end
    end
end

# Real matrix times complex vector.
# Multiply the matrix with the real and imaginary parts separately
@inline mul!(y::StridedVector{Complex{T}}, A::StridedMaybeAdjOrTransMat{T}, x::StridedVector{Complex{T}},
        alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemv!(y, A isa StridedArray ? 'N' : 'T', A isa StridedArray ? A : parent(A), x, alpha, beta)

@inline mul!(y::AbstractVector, A::AbstractVecOrMat, x::AbstractVector,
             alpha::Number, beta::Number) =
    generic_matvecmul!(y, 'N', A, x, MulAddMul(alpha, beta))

function *(tA::Transpose{<:Any,<:StridedMatrix{T}}, x::StridedVector{S}) where {T<:BlasFloat,S}
    TS = promote_op(matprod, T, S)
    mul!(similar(x, TS, size(tA, 1)), tA, convert(AbstractVector{TS}, x))
end
function *(tA::Transpose{<:Any,<:AbstractMatrix{T}}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(similar(x, TS, size(tA, 1)), tA, x)
end
@inline mul!(y::StridedVector{T}, tA::Transpose{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T},
                      alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemv!(y, 'T', tA.parent, x, alpha, beta)
@inline mul!(y::AbstractVector, tA::Transpose{<:Any,<:AbstractVecOrMat}, x::AbstractVector,
                      alpha::Number, beta::Number) =
    generic_matvecmul!(y, 'T', tA.parent, x, MulAddMul(alpha, beta))

function *(adjA::Adjoint{<:Any,<:StridedMatrix{T}}, x::StridedVector{S}) where {T<:BlasFloat,S}
    TS = promote_op(matprod, T, S)
    mul!(similar(x, TS, size(adjA, 1)), adjA, convert(AbstractVector{TS}, x))
end
function *(adjA::Adjoint{<:Any,<:AbstractMatrix{T}}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(similar(x, TS, size(adjA, 1)), adjA, x)
end

@inline mul!(y::StridedVector{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T},
                      alpha::Number, beta::Number) where {T<:BlasReal} =
    mul!(y, transpose(adjA.parent), x, alpha, beta)
@inline mul!(y::StridedVector{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, x::StridedVector{T},
                      alpha::Number, beta::Number) where {T<:BlasComplex} =
    gemv!(y, 'C', adjA.parent, x, alpha, beta)
@inline mul!(y::AbstractVector, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, x::AbstractVector,
                      alpha::Number, beta::Number) =
    generic_matvecmul!(y, 'C', adjA.parent, x, MulAddMul(alpha, beta))

# Vector-Matrix multiplication
(*)(x::AdjointAbsVec,   A::AbstractMatrix) = (A'*x')'
(*)(x::TransposeAbsVec, A::AbstractMatrix) = transpose(transpose(A)*transpose(x))

_parent(A) = A
_parent(A::Adjoint) = parent(A)
_parent(A::Transpose) = parent(A)

"""
    *(A::AbstractMatrix, B::AbstractMatrix)

Matrix multiplication.

# Examples
```jldoctest
julia> [1 1; 0 1] * [1 0; 1 1]
2×2 Matrix{Int64}:
 2  1
 1  1
```
"""
function (*)(A::AbstractMatrix, B::AbstractMatrix)
    TS = promote_op(matprod, eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A, 1), size(B, 2))), A, B)
end
# optimization for dispatching to BLAS, e.g. *(::Matrix{Float32}, ::Matrix{Float64})
# but avoiding the case *(::Matrix{<:BlasComplex}, ::Matrix{<:BlasReal})
# which is better handled by reinterpreting rather than promotion
function (*)(A::StridedMaybeAdjOrTransMat{<:BlasReal}, B::StridedMaybeAdjOrTransMat{<:BlasReal})
    TS = promote_type(eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A, 1), size(B, 2))),
         wrapperop(A)(convert(AbstractArray{TS}, _parent(A))),
         wrapperop(B)(convert(AbstractArray{TS}, _parent(B))))
end
function (*)(A::StridedMaybeAdjOrTransMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMat{<:BlasComplex})
    TS = promote_type(eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A, 1), size(B, 2))),
         wrapperop(A)(convert(AbstractArray{TS}, _parent(A))),
         wrapperop(B)(convert(AbstractArray{TS}, _parent(B))))
end

@inline function mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, B::StridedVecOrMat{T},
                      alpha::Number, beta::Number) where {T<:BlasFloat}
    return gemm_wrapper!(C, 'N', 'N', A, B, MulAddMul(alpha, beta))
end
# Complex Matrix times real matrix: We use that it is generally faster to reinterpret the
# first matrix as a real matrix and carry out real matrix matrix multiply
function (*)(A::StridedMatrix{<:BlasComplex}, B::StridedMaybeAdjOrTransMat{<:BlasReal})
    TS = promote_type(eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A, 1), size(B, 2))),
         convert(AbstractArray{TS}, A),
         wrapperop(B)(convert(AbstractArray{real(TS)}, _parent(B))))
end
function (*)(A::AdjOrTransStridedMat{<:BlasComplex}, B::StridedMaybeAdjOrTransMat{<:BlasReal})
    TS = promote_type(eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A, 1), size(B, 2))),
         copy_oftype(A, TS), # remove AdjOrTrans to use reinterpret trick below
         wrapperop(B)(convert(AbstractArray{real(TS)}, _parent(B))))
end
# the following case doesn't seem to benefit from the translation A*B = (B' * A')'
function (*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasComplex})
    temp = real(B)
    R = A * temp
    temp .= imag.(B)
    I = A * temp
    Complex.(R, I)
end
(*)(A::AdjOrTransStridedMat{<:BlasReal}, B::StridedMatrix{<:BlasComplex}) = copy(transpose(transpose(B) * parent(A)))
(*)(A::StridedMaybeAdjOrTransMat{<:BlasReal}, B::AdjOrTransStridedMat{<:BlasComplex}) = copy(wrapperop(B)(parent(B) * transpose(A)))

for elty in (Float32,Float64)
    @eval begin
        @inline function mul!(C::StridedMatrix{Complex{$elty}}, A::StridedVecOrMat{Complex{$elty}}, B::StridedVecOrMat{$elty},
                              alpha::Real, beta::Real)
            Afl = reinterpret($elty, A)
            Cfl = reinterpret($elty, C)
            mul!(Cfl, Afl, B, alpha, beta)
            return C
        end
    end
end

"""
    muladd(A, y, z)

Combined multiply-add, `A*y .+ z`, for matrix-matrix or matrix-vector multiplication.
The result is always the same size as `A*y`, but `z` may be smaller, or a scalar.

!!! compat "Julia 1.6"
     These methods require Julia 1.6 or later.

# Examples
```jldoctest
julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; z=[0, 100];

julia> muladd(A, B, z)
2×2 Matrix{Float64}:
   3.0    3.0
 107.0  107.0
```
"""
function Base.muladd(A::AbstractMatrix, y::AbstractVecOrMat, z::Union{Number, AbstractArray})
    Ay = A * y
    for d in 1:ndims(Ay)
        # Same error as Ay .+= z would give, to match StridedMatrix method:
        size(z,d) > size(Ay,d) && throw(DimensionMismatch("array could not be broadcast to match destination"))
    end
    for d in ndims(Ay)+1:ndims(z)
        # Similar error to what Ay + z would give, to match (Any,Any,Any) method:
        size(z,d) > 1 && throw(DimensionMismatch(string("dimensions must match: z has dims ",
            axes(z), ", must have singleton at dim ", d)))
    end
    Ay .+ z
end

function Base.muladd(u::AbstractVector, v::AdjOrTransAbsVec, z::Union{Number, AbstractArray})
    if size(z,1) > length(u) || size(z,2) > length(v)
        # Same error as (u*v) .+= z:
        throw(DimensionMismatch("array could not be broadcast to match destination"))
    end
    for d in 3:ndims(z)
        # Similar error to (u*v) + z:
        size(z,d) > 1 && throw(DimensionMismatch(string("dimensions must match: z has dims ",
            axes(z), ", must have singleton at dim ", d)))
    end
    (u .* v) .+ z
end

Base.muladd(x::AdjointAbsVec, A::AbstractMatrix, z::Union{Number, AbstractVecOrMat}) =
    muladd(A', x', z')'
Base.muladd(x::TransposeAbsVec, A::AbstractMatrix, z::Union{Number, AbstractVecOrMat}) =
    transpose(muladd(transpose(A), transpose(x), transpose(z)))

function Base.muladd(A::StridedMaybeAdjOrTransMat{<:Number}, y::AbstractVector{<:Number}, z::Union{Number, AbstractVector})
    T = promote_type(eltype(A), eltype(y), eltype(z))
    C = similar(A, T, axes(A,1))
    C .= z
    mul!(C, A, y, true, true)
end

function Base.muladd(A::StridedMaybeAdjOrTransMat{<:Number}, B::StridedMaybeAdjOrTransMat{<:Number}, z::Union{Number, AbstractVecOrMat})
    T = promote_type(eltype(A), eltype(B), eltype(z))
    C = similar(A, T, axes(A,1), axes(B,2))
    C .= z
    mul!(C, A, B, true, true)
end

"""
    mul!(Y, A, B) -> Y

Calculates the matrix-matrix or matrix-vector product ``AB`` and stores the result in `Y`,
overwriting the existing value of `Y`. Note that `Y` must not be aliased with either `A` or
`B`.

# Examples
```jldoctest
julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); mul!(Y, A, B);

julia> Y
2×2 Matrix{Float64}:
 3.0  3.0
 7.0  7.0
```

# Implementation
For custom matrix and vector types, it is recommended to implement
5-argument `mul!` rather than implementing 3-argument `mul!` directly
if possible.
"""
@inline function mul!(C, A, B)
    return mul!(C, A, B, true, false)
end

"""
    mul!(C, A, B, α, β) -> C

Combined inplace matrix-matrix or matrix-vector multiply-add ``A B α + C β``.
The result is stored in `C` by overwriting it.  Note that `C` must not be
aliased with either `A` or `B`.

!!! compat "Julia 1.3"
    Five-argument `mul!` requires at least Julia 1.3.

# Examples
```jldoctest
julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; C=[1.0 2.0; 3.0 4.0];

julia> mul!(C, A, B, 100.0, 10.0) === C
true

julia> C
2×2 Matrix{Float64}:
 310.0  320.0
 730.0  740.0
```
"""
@inline mul!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat,
             alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'N', 'N', A, B, MulAddMul(alpha, beta))

"""
    rmul!(A, B)

Calculate the matrix-matrix product ``AB``, overwriting `A`, and return the result.
Here, `B` must be of special matrix type, like, e.g., [`Diagonal`](@ref),
[`UpperTriangular`](@ref) or [`LowerTriangular`](@ref), or of some orthogonal type,
see [`QR`](@ref).

# Examples
```jldoctest
julia> A = [0 1; 1 0];

julia> B = LinearAlgebra.UpperTriangular([1 2; 0 3]);

julia> LinearAlgebra.rmul!(A, B);

julia> A
2×2 Matrix{Int64}:
 0  3
 1  2

julia> A = [1.0 2.0; 3.0 4.0];

julia> F = qr([0 1; -1 0]);

julia> rmul!(A, F.Q)
2×2 Matrix{Float64}:
 2.0  1.0
 4.0  3.0
```
"""
rmul!(A, B)

"""
    lmul!(A, B)

Calculate the matrix-matrix product ``AB``, overwriting `B`, and return the result.
Here, `A` must be of special matrix type, like, e.g., [`Diagonal`](@ref),
[`UpperTriangular`](@ref) or [`LowerTriangular`](@ref), or of some orthogonal type,
see [`QR`](@ref).

# Examples
```jldoctest
julia> B = [0 1; 1 0];

julia> A = LinearAlgebra.UpperTriangular([1 2; 0 3]);

julia> LinearAlgebra.lmul!(A, B);

julia> B
2×2 Matrix{Int64}:
 2  1
 3  0

julia> B = [1.0 2.0; 3.0 4.0];

julia> F = qr([0 1; -1 0]);

julia> lmul!(F.Q, B)
2×2 Matrix{Float64}:
 3.0  4.0
 1.0  2.0
```
"""
lmul!(A, B)

@inline function mul!(C::StridedMatrix{T}, tA::Transpose{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T},
                 alpha::Number, beta::Number) where {T<:BlasFloat}
    A = tA.parent
    if A === B
        return syrk_wrapper!(C, 'T', A, MulAddMul(alpha, beta))
    else
        return gemm_wrapper!(C, 'T', 'N', A, B, MulAddMul(alpha, beta))
    end
end
@inline mul!(C::AbstractMatrix, tA::Transpose{<:Any,<:AbstractVecOrMat}, B::AbstractVecOrMat,
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'T', 'N', tA.parent, B, MulAddMul(alpha, beta))

@inline function mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, tB::Transpose{<:Any,<:StridedVecOrMat{T}},
                 alpha::Number, beta::Number) where {T<:BlasFloat}
    B = tB.parent
    if A === B
        return syrk_wrapper!(C, 'N', A, MulAddMul(alpha, beta))
    else
        return gemm_wrapper!(C, 'N', 'T', A, B, MulAddMul(alpha, beta))
    end
end
# Complex matrix times transposed real matrix. Reinterpret the first matrix to real for efficiency.
for elty in (Float32,Float64)
    @eval begin
        @inline function mul!(C::StridedMatrix{Complex{$elty}}, A::StridedVecOrMat{Complex{$elty}}, tB::Transpose{<:Any,<:StridedVecOrMat{$elty}},
                         alpha::Real, beta::Real)
            Afl = reinterpret($elty, A)
            Cfl = reinterpret($elty, C)
            mul!(Cfl, Afl, tB, alpha, beta)
            return C
        end
    end
end
# collapsing the following two defs with C::AbstractVecOrMat yields ambiguities
@inline mul!(C::AbstractVector, A::AbstractVecOrMat, tB::Transpose{<:Any,<:AbstractVecOrMat},
             alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'N', 'T', A, tB.parent, MulAddMul(alpha, beta))
@inline mul!(C::AbstractMatrix, A::AbstractVecOrMat, tB::Transpose{<:Any,<:AbstractVecOrMat},
             alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'N', 'T', A, tB.parent, MulAddMul(alpha, beta))

@inline mul!(C::StridedMatrix{T}, tA::Transpose{<:Any,<:StridedVecOrMat{T}}, tB::Transpose{<:Any,<:StridedVecOrMat{T}},
                 alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemm_wrapper!(C, 'T', 'T', tA.parent, tB.parent, MulAddMul(alpha, beta))
@inline mul!(C::AbstractMatrix, tA::Transpose{<:Any,<:AbstractVecOrMat}, tB::Transpose{<:Any,<:AbstractVecOrMat},
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'T', 'T', tA.parent, tB.parent, MulAddMul(alpha, beta))

@inline mul!(C::StridedMatrix{T}, tA::Transpose{<:Any,<:StridedVecOrMat{T}}, adjB::Adjoint{<:Any,<:StridedVecOrMat{T}},
                 alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemm_wrapper!(C, 'T', 'C', tA.parent, adjB.parent, MulAddMul(alpha, beta))
@inline mul!(C::AbstractMatrix, tA::Transpose{<:Any,<:AbstractVecOrMat}, tB::Adjoint{<:Any,<:AbstractVecOrMat},
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'T', 'C', tA.parent, tB.parent, MulAddMul(alpha, beta))

@inline mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T},
                 alpha::Real, beta::Real) where {T<:BlasReal} =
    mul!(C, transpose(adjA.parent), B, alpha, beta)
@inline function mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, B::StridedVecOrMat{T},
                 alpha::Number, beta::Number) where {T<:BlasComplex}
    A = adjA.parent
    if A === B
        return herk_wrapper!(C, 'C', A, MulAddMul(alpha, beta))
    else
        return gemm_wrapper!(C, 'C', 'N', A, B, MulAddMul(alpha, beta))
    end
end
@inline mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, B::AbstractVecOrMat,
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'C', 'N', adjA.parent, B, MulAddMul(alpha, beta))

@inline mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:StridedVecOrMat{<:BlasReal}},
                 alpha::Number, beta::Number) where {T<:BlasFloat} =
    mul!(C, A, transpose(adjB.parent), alpha, beta)
@inline function mul!(C::StridedMatrix{T}, A::StridedVecOrMat{T}, adjB::Adjoint{<:Any,<:StridedVecOrMat{T}},
                 alpha::Number, beta::Number) where {T<:BlasComplex}
    B = adjB.parent
    if A === B
        return herk_wrapper!(C, 'N', A, MulAddMul(alpha, beta))
    else
        return gemm_wrapper!(C, 'N', 'C', A, B, MulAddMul(alpha, beta))
    end
end
@inline mul!(C::AbstractMatrix, A::AbstractVecOrMat, adjB::Adjoint{<:Any,<:AbstractVecOrMat},
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'N', 'C', A, adjB.parent, MulAddMul(alpha, beta))

@inline mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, adjB::Adjoint{<:Any,<:StridedVecOrMat{T}},
                 alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemm_wrapper!(C, 'C', 'C', adjA.parent, adjB.parent, MulAddMul(alpha, beta))
@inline mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, adjB::Adjoint{<:Any,<:AbstractVecOrMat},
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'C', 'C', adjA.parent, adjB.parent, MulAddMul(alpha, beta))

@inline mul!(C::StridedMatrix{T}, adjA::Adjoint{<:Any,<:StridedVecOrMat{T}}, tB::Transpose{<:Any,<:StridedVecOrMat{T}},
                 alpha::Number, beta::Number) where {T<:BlasFloat} =
    gemm_wrapper!(C, 'C', 'T', adjA.parent, tB.parent, MulAddMul(alpha, beta))
@inline mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, tB::Transpose{<:Any,<:AbstractVecOrMat},
                 alpha::Number, beta::Number) =
    generic_matmatmul!(C, 'C', 'T', adjA.parent, tB.parent, MulAddMul(alpha, beta))

# Supporting functions for matrix multiplication

# copy transposed(adjoint) of upper(lower) side-digonals. Optionally include diagonal.
@inline function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false, diag::Bool=false)
    n = checksquare(A)
    off = diag ? 0 : 1
    if uplo == 'U'
        for i = 1:n, j = (i+off):n
            A[j,i] = conjugate ? adjoint(A[i,j]) : transpose(A[i,j])
        end
    elseif uplo == 'L'
        for i = 1:n, j = (i+off):n
            A[i,j] = conjugate ? adjoint(A[j,i]) : transpose(A[j,i])
        end
    else
        throw(ArgumentError(lazy"uplo argument must be 'U' (upper) or 'L' (lower), got $uplo"))
    end
    A
end

function gemv!(y::StridedVector{T}, tA::AbstractChar, A::StridedVecOrMat{T}, x::StridedVector{T},
               α::Number=true, β::Number=false) where {T<:BlasFloat}
    mA, nA = lapack_size(tA, A)
    if nA != length(x)
        throw(DimensionMismatch(lazy"second dimension of A, $nA, does not match length of x, $(length(x))"))
    end
    if mA != length(y)
        throw(DimensionMismatch(lazy"first dimension of A, $mA, does not match length of y, $(length(y))"))
    end
    if mA == 0
        return y
    end
    if nA == 0
        return _rmul_or_fill!(y, β)
    end

    alpha, beta = promote(α, β, zero(T))
    if alpha isa Union{Bool,T} && beta isa Union{Bool,T} && stride(A, 1) == 1 && stride(A, 2) >= size(A, 1)
        return BLAS.gemv!(tA, alpha, A, x, beta, y)
    else
        return generic_matvecmul!(y, tA, A, x, MulAddMul(α, β))
    end
end

function gemv!(y::StridedVector{Complex{T}}, tA::AbstractChar, A::StridedVecOrMat{T}, x::StridedVector{Complex{T}},
    α::Number = true, β::Number = false) where {T<:BlasFloat}
    mA, nA = lapack_size(tA, A)
    nA != length(x) &&
        throw(DimensionMismatch(lazy"second dimension of A, $nA, does not match length of x, $(length(x))"))
    mA != length(y) &&
        throw(DimensionMismatch(lazy"first dimension of A, $mA, does not match length of y, $(length(y))"))
    mA == 0 && return y
    nA == 0 && return _rmul_or_fill!(y, β)
    alpha, beta = promote(α, β, zero(T))
    @views if alpha isa Union{Bool,T} && beta isa Union{Bool,T} && stride(A, 1) == 1 && stride(A, 2) >= size(A, 1)
        xfl = reinterpret(reshape, T, x) # Use reshape here.
        yfl = reinterpret(reshape, T, y)
        BLAS.gemv!(tA, alpha, A, xfl[1, :], beta, yfl[1, :])
        BLAS.gemv!(tA, alpha, A, xfl[2, :], beta, yfl[2, :])
        return y
    else
        return generic_matvecmul!(y, tA, A, x, MulAddMul(α, β))
    end
end

function syrk_wrapper!(C::StridedMatrix{T}, tA::AbstractChar, A::StridedVecOrMat{T},
        _add = MulAddMul()) where {T<:BlasFloat}
    nC = checksquare(C)
    if tA == 'T'
        (nA, mA) = size(A,1), size(A,2)
        tAt = 'N'
    else
        (mA, nA) = size(A,1), size(A,2)
        tAt = 'T'
    end
    if nC != mA
        throw(DimensionMismatch(lazy"output matrix has size: $(nC), but should have size $(mA)"))
    end
    if mA == 0 || nA == 0 || iszero(_add.alpha)
        return _rmul_or_fill!(C, _add.beta)
    end
    if mA == 2 && nA == 2
        return matmul2x2!(C, tA, tAt, A, A, _add)
    end
    if mA == 3 && nA == 3
        return matmul3x3!(C, tA, tAt, A, A, _add)
    end

    # BLAS.syrk! only updates symmetric C
    # alternatively, make non-zero β a show-stopper for BLAS.syrk!
    if iszero(_add.beta) || issymmetric(C)
        alpha, beta = promote(_add.alpha, _add.beta, zero(T))
        if (alpha isa Union{Bool,T} &&
            beta isa Union{Bool,T} &&
            stride(A, 1) == stride(C, 1) == 1 &&
            stride(A, 2) >= size(A, 1) &&
            stride(C, 2) >= size(C, 1))
            return copytri!(BLAS.syrk!('U', tA, alpha, A, beta, C), 'U')
        end
    end
    return gemm_wrapper!(C, tA, tAt, A, A, _add)
end

function herk_wrapper!(C::Union{StridedMatrix{T}, StridedMatrix{Complex{T}}}, tA::AbstractChar, A::Union{StridedVecOrMat{T}, StridedVecOrMat{Complex{T}}},
        _add = MulAddMul()) where {T<:BlasReal}
    nC = checksquare(C)
    if tA == 'C'
        (nA, mA) = size(A,1), size(A,2)
        tAt = 'N'
    else
        (mA, nA) = size(A,1), size(A,2)
        tAt = 'C'
    end
    if nC != mA
        throw(DimensionMismatch(lazy"output matrix has size: $(nC), but should have size $(mA)"))
    end
    if mA == 0 || nA == 0 || iszero(_add.alpha)
        return _rmul_or_fill!(C, _add.beta)
    end
    if mA == 2 && nA == 2
        return matmul2x2!(C, tA, tAt, A, A, _add)
    end
    if mA == 3 && nA == 3
        return matmul3x3!(C, tA, tAt, A, A, _add)
    end

    # Result array does not need to be initialized as long as beta==0
    #    C = Matrix{T}(undef, mA, mA)

    if iszero(_add.beta) || issymmetric(C)
        alpha, beta = promote(_add.alpha, _add.beta, zero(T))
        if (alpha isa Union{Bool,T} &&
            beta isa Union{Bool,T} &&
            stride(A, 1) == stride(C, 1) == 1 &&
            stride(A, 2) >= size(A, 1) &&
            stride(C, 2) >= size(C, 1))
            return copytri!(BLAS.herk!('U', tA, alpha, A, beta, C), 'U', true)
        end
    end
    return gemm_wrapper!(C, tA, tAt, A, A, _add)
end

function gemm_wrapper(tA::AbstractChar, tB::AbstractChar,
                      A::StridedVecOrMat{T},
                      B::StridedVecOrMat{T}) where {T<:BlasFloat}
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    C = similar(B, T, mA, nB)
    gemm_wrapper!(C, tA, tB, A, B)
end

function gemm_wrapper!(C::StridedVecOrMat{T}, tA::AbstractChar, tB::AbstractChar,
                       A::StridedVecOrMat{T}, B::StridedVecOrMat{T},
                       _add = MulAddMul()) where {T<:BlasFloat}
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)

    if nA != mB
        throw(DimensionMismatch(lazy"A has dimensions ($mA,$nA) but B has dimensions ($mB,$nB)"))
    end

    if C === A || B === C
        throw(ArgumentError("output matrix must not be aliased with input matrix"))
    end

    if mA == 0 || nA == 0 || nB == 0 || iszero(_add.alpha)
        if size(C) != (mA, nB)
            throw(DimensionMismatch(lazy"C has dimensions $(size(C)), should have ($mA,$nB)"))
        end
        return _rmul_or_fill!(C, _add.beta)
    end

    if mA == 2 && nA == 2 && nB == 2
        return matmul2x2!(C, tA, tB, A, B, _add)
    end
    if mA == 3 && nA == 3 && nB == 3
        return matmul3x3!(C, tA, tB, A, B, _add)
    end

    alpha, beta = promote(_add.alpha, _add.beta, zero(T))
    if (alpha isa Union{Bool,T} &&
        beta isa Union{Bool,T} &&
        stride(A, 1) == stride(B, 1) == stride(C, 1) == 1 &&
        stride(A, 2) >= size(A, 1) &&
        stride(B, 2) >= size(B, 1) &&
        stride(C, 2) >= size(C, 1))
        return BLAS.gemm!(tA, tB, alpha, A, B, beta, C)
    end
    generic_matmatmul!(C, tA, tB, A, B, _add)
end

# blas.jl defines matmul for floats; other integer and mixed precision
# cases are handled here

lapack_size(t::AbstractChar, M::AbstractVecOrMat) = (size(M, t=='N' ? 1 : 2), size(M, t=='N' ? 2 : 1))

function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int})
    if tM == 'N'
        copyto!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    else
        LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src)
        tM == 'C' && conj!(@view B[ir_dest, jr_dest])
    end
    B
end

function copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int})
    if tM == 'N'
        LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src)
    else
        copyto!(B, ir_dest, jr_dest, M, jr_src, ir_src)
        tM == 'C' && conj!(@view B[ir_dest, jr_dest])
    end
    B
end

# TODO: It will be faster for large matrices to convert to float,
# call BLAS, and convert back to required type.

# NOTE: the generic version is also called as fallback for
#       strides != 1 cases

function generic_matvecmul!(C::AbstractVector{R}, tA, A::AbstractVecOrMat, B::AbstractVector,
                            _add::MulAddMul = MulAddMul()) where R
    require_one_based_indexing(C, A, B)
    mB = length(B)
    mA, nA = lapack_size(tA, A)
    if mB != nA
        throw(DimensionMismatch(lazy"matrix A has dimensions ($mA,$nA), vector B has length $mB"))
    end
    if mA != length(C)
        throw(DimensionMismatch(lazy"result C has length $(length(C)), needs length $mA"))
    end

    Astride = size(A, 1)

    @inbounds begin
    if tA == 'T'  # fastest case
        if nA == 0
            for k = 1:mA
                _modify!(_add, false, C, k)
            end
        else
            for k = 1:mA
                aoffs = (k-1)*Astride
                s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
                for i = 1:nA
                    s += transpose(A[aoffs+i]) * B[i]
                end
                _modify!(_add, s, C, k)
            end
        end
    elseif tA == 'C'
        if nA == 0
            for k = 1:mA
                _modify!(_add, false, C, k)
            end
        else
            for k = 1:mA
                aoffs = (k-1)*Astride
                s = zero(A[aoffs + 1]*B[1] + A[aoffs + 1]*B[1])
                for i = 1:nA
                    s += A[aoffs + i]'B[i]
                end
                _modify!(_add, s, C, k)
            end
        end
    else # tA == 'N'
        for i = 1:mA
            if !iszero(_add.beta)
                C[i] *= _add.beta
            elseif mB == 0
                C[i] = false
            else
                C[i] = zero(A[i]*B[1] + A[i]*B[1])
            end
        end
        for k = 1:mB
            aoffs = (k-1)*Astride
            b = _add(B[k], false)
            for i = 1:mA
                C[i] += A[aoffs + i] * b
            end
        end
    end
    end # @inbounds
    C
end

function generic_matmatmul(tA, tB, A::AbstractVecOrMat{T}, B::AbstractMatrix{S}) where {T,S}
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    C = similar(B, promote_op(matprod, T, S), mA, nB)
    generic_matmatmul!(C, tA, tB, A, B)
end

const tilebufsize = 10800  # Approximately 32k/3

function generic_matmatmul!(C::AbstractMatrix, tA, tB, A::AbstractMatrix, B::AbstractMatrix,
                            _add::MulAddMul=MulAddMul())
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    mC, nC = size(C)

    if iszero(_add.alpha)
        return _rmul_or_fill!(C, _add.beta)
    end
    if mA == nA == mB == nB == mC == nC == 2
        return matmul2x2!(C, tA, tB, A, B, _add)
    end
    if mA == nA == mB == nB == mC == nC == 3
        return matmul3x3!(C, tA, tB, A, B, _add)
    end
    _generic_matmatmul!(C, tA, tB, A, B, _add)
end

generic_matmatmul!(C::AbstractVecOrMat, tA, tB, A::AbstractVecOrMat, B::AbstractVecOrMat, _add::MulAddMul) =
    _generic_matmatmul!(C, tA, tB, A, B, _add)

function _generic_matmatmul!(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVecOrMat{T}, B::AbstractVecOrMat{S},
                             _add::MulAddMul) where {T,S,R}
    require_one_based_indexing(C, A, B)
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    if mB != nA
        throw(DimensionMismatch(lazy"matrix A has dimensions ($mA,$nA), matrix B has dimensions ($mB,$nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch(lazy"result C has dimensions $(size(C)), needs ($mA,$nB)"))
    end

    if iszero(_add.alpha) || isempty(A) || isempty(B)
        return _rmul_or_fill!(C, _add.beta)
    end

    tile_size = 0
    if isbitstype(R) && isbitstype(T) && isbitstype(S) && (tA == 'N' || tB != 'N')
        tile_size = floor(Int, sqrt(tilebufsize / max(sizeof(R), sizeof(S), sizeof(T), 1)))
    end
    @inbounds begin
    if tile_size > 0
        sz = (tile_size, tile_size)
        Atile = Array{T}(undef, sz)
        Btile = Array{S}(undef, sz)

        z1 = zero(A[1, 1]*B[1, 1] + A[1, 1]*B[1, 1])
        z = convert(promote_type(typeof(z1), R), z1)

        if mA < tile_size && nA < tile_size && nB < tile_size
            copy_transpose!(Atile, 1:nA, 1:mA, tA, A, 1:mA, 1:nA)
            copyto!(Btile, 1:mB, 1:nB, tB, B, 1:mB, 1:nB)
            for j = 1:nB
                boff = (j-1)*tile_size
                for i = 1:mA
                    aoff = (i-1)*tile_size
                    s = z
                    for k = 1:nA
                        s += Atile[aoff+k] * Btile[boff+k]
                    end
                    _modify!(_add, s, C, (i,j))
                end
            end
        else
            Ctile = Array{R}(undef, sz)
            for jb = 1:tile_size:nB
                jlim = min(jb+tile_size-1,nB)
                jlen = jlim-jb+1
                for ib = 1:tile_size:mA
                    ilim = min(ib+tile_size-1,mA)
                    ilen = ilim-ib+1
                    fill!(Ctile, z)
                    for kb = 1:tile_size:nA
                        klim = min(kb+tile_size-1,mB)
                        klen = klim-kb+1
                        copy_transpose!(Atile, 1:klen, 1:ilen, tA, A, ib:ilim, kb:klim)
                        copyto!(Btile, 1:klen, 1:jlen, tB, B, kb:klim, jb:jlim)
                        for j=1:jlen
                            bcoff = (j-1)*tile_size
                            for i = 1:ilen
                                aoff = (i-1)*tile_size
                                s = z
                                for k = 1:klen
                                    s += Atile[aoff+k] * Btile[bcoff+k]
                                end
                                Ctile[bcoff+i] += s
                            end
                        end
                    end
                    if isone(_add.alpha) && iszero(_add.beta)
                        copyto!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                    else
                        C[ib:ilim, jb:jlim] .= @views _add.(Ctile[1:ilen, 1:jlen], C[ib:ilim, jb:jlim])
                    end
                end
            end
        end
    else
        # Multiplication for non-plain-data uses the naive algorithm
        if tA == 'N'
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    z2 = zero(A[i, 1]*B[1, j] + A[i, 1]*B[1, j])
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += A[i, k]*B[k, j]
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    z2 = zero(A[i, 1]*transpose(B[j, 1]) + A[i, 1]*transpose(B[j, 1]))
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += A[i, k] * transpose(B[j, k])
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            else
                for i = 1:mA, j = 1:nB
                    z2 = zero(A[i, 1]*B[j, 1]' + A[i, 1]*B[j, 1]')
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += A[i, k]*B[j, k]'
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            end
        elseif tA == 'T'
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    z2 = zero(transpose(A[1, i])*B[1, j] + transpose(A[1, i])*B[1, j])
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += transpose(A[k, i]) * B[k, j]
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    z2 = zero(transpose(A[1, i])*transpose(B[j, 1]) + transpose(A[1, i])*transpose(B[j, 1]))
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += transpose(A[k, i]) * transpose(B[j, k])
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            else
                for i = 1:mA, j = 1:nB
                    z2 = zero(transpose(A[1, i])*B[j, 1]' + transpose(A[1, i])*B[j, 1]')
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += transpose(A[k, i]) * adjoint(B[j, k])
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            end
        else
            if tB == 'N'
                for i = 1:mA, j = 1:nB
                    z2 = zero(A[1, i]'*B[1, j] + A[1, i]'*B[1, j])
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += A[k, i]'B[k, j]
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            elseif tB == 'T'
                for i = 1:mA, j = 1:nB
                    z2 = zero(A[1, i]'*transpose(B[j, 1]) + A[1, i]'*transpose(B[j, 1]))
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += adjoint(A[k, i]) * transpose(B[j, k])
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            else
                for i = 1:mA, j = 1:nB
                    z2 = zero(A[1, i]'*B[j, 1]' + A[1, i]'*B[j, 1]')
                    Ctmp = convert(promote_type(R, typeof(z2)), z2)
                    for k = 1:nA
                        Ctmp += A[k, i]'B[j, k]'
                    end
                    _modify!(_add, Ctmp, C, (i,j))
                end
            end
        end
    end
    end # @inbounds
    C
end

# multiply 2x2 matrices
function matmul2x2(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S}) where {T,S}
    matmul2x2!(similar(B, promote_op(matprod, T, S), 2, 2), tA, tB, A, B)
end

function matmul2x2!(C::AbstractMatrix, tA, tB, A::AbstractMatrix, B::AbstractMatrix,
                    _add::MulAddMul = MulAddMul())
    require_one_based_indexing(C, A, B)
    if !(size(A) == size(B) == size(C) == (2,2))
        throw(DimensionMismatch(lazy"A has size $(size(A)), B has size $(size(B)), C has size $(size(C))"))
    end
    @inbounds begin
    if tA == 'T'
        # TODO making these lazy could improve perf
        A11 = copy(transpose(A[1,1])); A12 = copy(transpose(A[2,1]))
        A21 = copy(transpose(A[1,2])); A22 = copy(transpose(A[2,2]))
    elseif tA == 'C'
        # TODO making these lazy could improve perf
        A11 = copy(A[1,1]'); A12 = copy(A[2,1]')
        A21 = copy(A[1,2]'); A22 = copy(A[2,2]')
    else
        A11 = A[1,1]; A12 = A[1,2]; A21 = A[2,1]; A22 = A[2,2]
    end
    if tB == 'T'
        # TODO making these lazy could improve perf
        B11 = copy(transpose(B[1,1])); B12 = copy(transpose(B[2,1]))
        B21 = copy(transpose(B[1,2])); B22 = copy(transpose(B[2,2]))
    elseif tB == 'C'
        # TODO making these lazy could improve perf
        B11 = copy(B[1,1]'); B12 = copy(B[2,1]')
        B21 = copy(B[1,2]'); B22 = copy(B[2,2]')
    else
        B11 = B[1,1]; B12 = B[1,2];
        B21 = B[2,1]; B22 = B[2,2]
    end
    _modify!(_add, A11*B11 + A12*B21, C, (1,1))
    _modify!(_add, A11*B12 + A12*B22, C, (1,2))
    _modify!(_add, A21*B11 + A22*B21, C, (2,1))
    _modify!(_add, A21*B12 + A22*B22, C, (2,2))
    end # inbounds
    C
end

# Multiply 3x3 matrices
function matmul3x3(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S}) where {T,S}
    matmul3x3!(similar(B, promote_op(matprod, T, S), 3, 3), tA, tB, A, B)
end

function matmul3x3!(C::AbstractMatrix, tA, tB, A::AbstractMatrix, B::AbstractMatrix,
                    _add::MulAddMul = MulAddMul())
    require_one_based_indexing(C, A, B)
    if !(size(A) == size(B) == size(C) == (3,3))
        throw(DimensionMismatch(lazy"A has size $(size(A)), B has size $(size(B)), C has size $(size(C))"))
    end
    @inbounds begin
    if tA == 'T'
        # TODO making these lazy could improve perf
        A11 = copy(transpose(A[1,1])); A12 = copy(transpose(A[2,1])); A13 = copy(transpose(A[3,1]))
        A21 = copy(transpose(A[1,2])); A22 = copy(transpose(A[2,2])); A23 = copy(transpose(A[3,2]))
        A31 = copy(transpose(A[1,3])); A32 = copy(transpose(A[2,3])); A33 = copy(transpose(A[3,3]))
    elseif tA == 'C'
        # TODO making these lazy could improve perf
        A11 = copy(A[1,1]'); A12 = copy(A[2,1]'); A13 = copy(A[3,1]')
        A21 = copy(A[1,2]'); A22 = copy(A[2,2]'); A23 = copy(A[3,2]')
        A31 = copy(A[1,3]'); A32 = copy(A[2,3]'); A33 = copy(A[3,3]')
    else
        A11 = A[1,1]; A12 = A[1,2]; A13 = A[1,3]
        A21 = A[2,1]; A22 = A[2,2]; A23 = A[2,3]
        A31 = A[3,1]; A32 = A[3,2]; A33 = A[3,3]
    end

    if tB == 'T'
        # TODO making these lazy could improve perf
        B11 = copy(transpose(B[1,1])); B12 = copy(transpose(B[2,1])); B13 = copy(transpose(B[3,1]))
        B21 = copy(transpose(B[1,2])); B22 = copy(transpose(B[2,2])); B23 = copy(transpose(B[3,2]))
        B31 = copy(transpose(B[1,3])); B32 = copy(transpose(B[2,3])); B33 = copy(transpose(B[3,3]))
    elseif tB == 'C'
        # TODO making these lazy could improve perf
        B11 = copy(B[1,1]'); B12 = copy(B[2,1]'); B13 = copy(B[3,1]')
        B21 = copy(B[1,2]'); B22 = copy(B[2,2]'); B23 = copy(B[3,2]')
        B31 = copy(B[1,3]'); B32 = copy(B[2,3]'); B33 = copy(B[3,3]')
    else
        B11 = B[1,1]; B12 = B[1,2]; B13 = B[1,3]
        B21 = B[2,1]; B22 = B[2,2]; B23 = B[2,3]
        B31 = B[3,1]; B32 = B[3,2]; B33 = B[3,3]
    end

    _modify!(_add, A11*B11 + A12*B21 + A13*B31, C, (1,1))
    _modify!(_add, A11*B12 + A12*B22 + A13*B32, C, (1,2))
    _modify!(_add, A11*B13 + A12*B23 + A13*B33, C, (1,3))

    _modify!(_add, A21*B11 + A22*B21 + A23*B31, C, (2,1))
    _modify!(_add, A21*B12 + A22*B22 + A23*B32, C, (2,2))
    _modify!(_add, A21*B13 + A22*B23 + A23*B33, C, (2,3))

    _modify!(_add, A31*B11 + A32*B21 + A33*B31, C, (3,1))
    _modify!(_add, A31*B12 + A32*B22 + A33*B32, C, (3,2))
    _modify!(_add, A31*B13 + A32*B23 + A33*B33, C, (3,3))
    end # inbounds
    C
end

const RealOrComplex = Union{Real,Complex}

# Three-argument *
"""
    *(A, B::AbstractMatrix, C)
    A * B * C * D

Chained multiplication of 3 or 4 matrices is done in the most efficient sequence,
based on the sizes of the arrays. That is, the number of scalar multiplications needed
for `(A * B) * C` (with 3 dense matrices) is compared to that for `A * (B * C)`
to choose which of these to execute.

If the last factor is a vector, or the first a transposed vector, then it is efficient
to deal with these first. In particular `x' * B * y` means `(x' * B) * y`
for an ordinary column-major `B::Matrix`. Unlike `dot(x, B, y)`, this
allocates an intermediate array.

If the first or last factor is a number, this will be fused with the matrix
multiplication, using 5-arg [`mul!`](@ref).

See also [`muladd`](@ref), [`dot`](@ref).

!!! compat "Julia 1.7"
    These optimisations require at least Julia 1.7.
"""
*(A::AbstractMatrix, B::AbstractMatrix, x::AbstractVector) = A * (B*x)

*(tu::AdjOrTransAbsVec, B::AbstractMatrix, v::AbstractVector) = (tu*B) * v
*(tu::AdjOrTransAbsVec, B::AdjOrTransAbsMat, v::AbstractVector) = tu * (B*v)

*(A::AbstractMatrix, x::AbstractVector, γ::Number) = mat_vec_scalar(A,x,γ)
*(A::AbstractMatrix, B::AbstractMatrix, γ::Number) = mat_mat_scalar(A,B,γ)
*(α::RealOrComplex, B::AbstractMatrix{<:RealOrComplex}, C::AbstractVector{<:RealOrComplex}) =
    mat_vec_scalar(B,C,α)
*(α::RealOrComplex, B::AbstractMatrix{<:RealOrComplex}, C::AbstractMatrix{<:RealOrComplex}) =
    mat_mat_scalar(B,C,α)

*(α::Number, u::AbstractVector, tv::AdjOrTransAbsVec) = broadcast(*, α, u, tv)
*(u::AbstractVector, tv::AdjOrTransAbsVec, γ::Number) = broadcast(*, u, tv, γ)
*(u::AbstractVector, tv::AdjOrTransAbsVec, C::AbstractMatrix) = u * (tv*C)

*(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix) = _tri_matmul(A,B,C)
*(tv::AdjOrTransAbsVec, B::AbstractMatrix, C::AbstractMatrix) = (tv*B) * C

function _tri_matmul(A,B,C,δ=nothing)
    n,m = size(A)
    # m,k == size(B)
    k,l = size(C)
    costAB_C = n*m*k + n*k*l  # multiplications, allocations n*k + n*l
    costA_BC = m*k*l + n*m*l  #                              m*l + n*l
    if costA_BC < costAB_C
        isnothing(δ) ? A * (B*C) : A * mat_mat_scalar(B,C,δ)
    else
        isnothing(δ) ? (A*B) * C : mat_mat_scalar(A*B, C, δ)
    end
end

# Fast path for two arrays * one scalar is opt-in, via mat_vec_scalar and mat_mat_scalar.

mat_vec_scalar(A, x, γ) = A * (x * γ)  # fallback
mat_vec_scalar(A::StridedMaybeAdjOrTransMat, x::StridedVector, γ) = _mat_vec_scalar(A, x, γ)
mat_vec_scalar(A::AdjOrTransAbsVec, x::StridedVector, γ) = (A * x) * γ

function _mat_vec_scalar(A, x, γ)
    T = promote_type(eltype(A), eltype(x), typeof(γ))
    C = similar(A, T, axes(A,1))
    mul!(C, A, x, γ, false)
end

mat_mat_scalar(A, B, γ) = (A*B) * γ # fallback
mat_mat_scalar(A::StridedMaybeAdjOrTransMat, B::StridedMaybeAdjOrTransMat, γ) =
    _mat_mat_scalar(A, B, γ)

function _mat_mat_scalar(A, B, γ)
    T = promote_type(eltype(A), eltype(B), typeof(γ))
    C = similar(A, T, axes(A,1), axes(B,2))
    mul!(C, A, B, γ, false)
end

mat_mat_scalar(A::AdjointAbsVec, B, γ) = (γ' * (A * B)')' # preserving order, adjoint reverses
mat_mat_scalar(A::AdjointAbsVec{<:RealOrComplex}, B::StridedMaybeAdjOrTransMat{<:RealOrComplex}, γ::RealOrComplex) =
    mat_vec_scalar(B', A', γ')'

mat_mat_scalar(A::TransposeAbsVec, B, γ) = transpose(γ * transpose(A * B))
mat_mat_scalar(A::TransposeAbsVec{<:RealOrComplex}, B::StridedMaybeAdjOrTransMat{<:RealOrComplex}, γ::RealOrComplex) =
    transpose(mat_vec_scalar(transpose(B), transpose(A), γ))

# Four-argument *, by type
*(α::Number, β::Number, C::AbstractMatrix, x::AbstractVector) = (α*β) * C * x
*(α::Number, β::Number, C::AbstractMatrix, D::AbstractMatrix) = (α*β) * C * D
*(α::Number, B::AbstractMatrix, C::AbstractMatrix, x::AbstractVector) = α * B * (C*x)
*(α::Number, vt::AdjOrTransAbsVec, C::AbstractMatrix, x::AbstractVector) = α * (vt*C*x)
*(α::RealOrComplex, vt::AdjOrTransAbsVec{<:RealOrComplex}, C::AbstractMatrix{<:RealOrComplex}, D::AbstractMatrix{<:RealOrComplex}) =
    (α*vt*C) * D # solves an ambiguity

*(A::AbstractMatrix, x::AbstractVector, γ::Number, δ::Number) = A * x * (γ*δ)
*(A::AbstractMatrix, B::AbstractMatrix, γ::Number, δ::Number) = A * B * (γ*δ)
*(A::AbstractMatrix, B::AbstractMatrix, x::AbstractVector, δ::Number, ) = A * (B*x*δ)
*(vt::AdjOrTransAbsVec, B::AbstractMatrix, x::AbstractVector, δ::Number) = (vt*B*x) * δ
*(vt::AdjOrTransAbsVec, B::AbstractMatrix, C::AbstractMatrix, δ::Number) = (vt*B) * C * δ

*(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, x::AbstractVector) = A * B * (C*x)
*(vt::AdjOrTransAbsVec, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix) = (vt*B) * C * D
*(vt::AdjOrTransAbsVec, B::AbstractMatrix, C::AbstractMatrix, x::AbstractVector) = vt * B * (C*x)

# Four-argument *, by size
*(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, δ::Number) = _tri_matmul(A,B,C,δ)
*(α::RealOrComplex, B::AbstractMatrix{<:RealOrComplex}, C::AbstractMatrix{<:RealOrComplex}, D::AbstractMatrix{<:RealOrComplex}) =
    _tri_matmul(B,C,D,α)
*(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix) =
    _quad_matmul(A,B,C,D)

function _quad_matmul(A,B,C,D)
    c1 = _mul_cost((A,B),(C,D))
    c2 = _mul_cost(((A,B),C),D)
    c3 = _mul_cost(A,(B,(C,D)))
    c4 = _mul_cost((A,(B,C)),D)
    c5 = _mul_cost(A,((B,C),D))
    cmin = min(c1,c2,c3,c4,c5)
    if c1 == cmin
        (A*B) * (C*D)
    elseif c2 == cmin
        ((A*B) * C) * D
    elseif c3 == cmin
        A * (B * (C*D))
    elseif c4 == cmin
        (A * (B*C)) * D
    else
        A * ((B*C) * D)
    end
end
@inline _mul_cost(A::AbstractMatrix) = 0
@inline _mul_cost((A,B)::Tuple) = _mul_cost(A,B)
@inline _mul_cost(A,B) = _mul_cost(A) + _mul_cost(B) + *(_mul_sizes(A)..., last(_mul_sizes(B)))
@inline _mul_sizes(A::AbstractMatrix) = size(A)
@inline _mul_sizes((A,B)::Tuple) = first(_mul_sizes(A)), last(_mul_sizes(B))

max_stars_repo_name: adharshkamath/julia

max_issues_repo_name: adharshkamath/julia

max_forks_repo_name: adharshkamath/julia


text:
mutable struct Variable <: Expression
    name::Symbol
    subtype::Symbol
    diff::Union{Function,Nothing}  # FIXME
    dependents::Vector{Variable}
end

Variable(name; subtype::Symbol, dependents::Vector{Variable} = Variable[]) =
    Variable(name, subtype, nothing, dependents)

Parameter(name; kwargs...) = Variable(name; subtype=:Parameter, kwargs...)
Unknown(name, ;kwargs...) = Variable(name; subtype=:Unknown, kwargs...)

export Variable, Unknown, Parameter, Constant, @Unknown, @Param, @Const

Base.copy(x::Variable) = Variable(x.name, x.subtype, x.diff, x.dependents)

struct Constant <: Expression
    value::Number
end
Base.get(c::Constant) = c.value

Base.iszero(ex::Expression) = isa(ex, Constant) && iszero(ex.value)
Base.isone(ex::Expression)  = isa(ex, Constant) && isone(ex.value)

# Variables use isequal for equality since == is an Operation
Base.:(==)(x::Variable, y::Variable) = (x.name, x.subtype, x.diff) == (y.name, y.subtype, y.diff)
Base.:(==)(::Variable, ::Number) = false
Base.:(==)(::Number, ::Variable) = false
Base.:(==)(::Variable, ::Constant) = false
Base.:(==)(::Constant, ::Variable) = false
Base.:(==)(c::Constant, n::Number) = c.value == n
Base.:(==)(n::Number, c::Constant) = c.value == n
Base.:(==)(a::Constant, b::Constant) = a.value == b.value

function Base.convert(::Type{Expr}, x::Variable)
    x.diff === nothing && return x.name
    return Symbol("$(x.name)_$(x.diff.x.name)")
end
Base.convert(::Type{Expr}, c::Constant) = c.value

function Base.show(io::IO, x::Variable)
    print(io, x.subtype, '(', x.name, ')')
    x.diff === nothing || print(io, ", diff = ", x.diff)
end

# Build variables more easily
function _parse_vars(macroname, fun, x)
    ex = Expr(:block)
    var_names = Symbol[]
    # if parsing things in the form of
    # begin
    #     x
    #     y
    #     z
    # end
    x = flatten_expr!(x)
    for _var in x
        iscall = isa(_var, Expr) && _var.head == :call
        issym    = _var isa Symbol
        @assert iscall || issym "@$macroname expects a tuple of expressions (`@$macroname x y z(t)`)"

        if iscall
            dependents = :(Variable[$(_var.args[2:end]...)])
            var_name = _var.args[1]
        else
            dependents = Variable[]
            var_name = _var
        end

        push!(var_names, var_name)
        expr = :($var_name = $fun($(Meta.quot(var_name)), dependents = $dependents))
        push!(ex.args, expr)
    end
    push!(ex.args, build_expr(:tuple, var_names))
    return ex
end
macro Unknown(xs...)
    esc(_parse_vars(:Unknown, Unknown, xs))
end
macro Param(xs...)
    esc(_parse_vars(:Param, Parameter, xs))
end

function _const_assign(x)
    ex = Expr(:block)
    lhss = Symbol[]
    x = flatten_expr!(x)
    for eq in x
        @assert eq isa Expr && eq.head == :(=) "@Const expects a tuple of assignments!\nE.g. `@Const D=t W=g`"
        lhs = eq.args[1]
        push!(lhss, lhs)
        rhs = eq.args[2]
        expr = :($lhs = Constant($rhs))
        push!(ex.args,  expr)
    end
    push!(ex.args, Expr(:tuple, lhss...))
    ex
end

macro Const(x...)
    esc(_const_assign(x))
end

ext: jl


text:
# Autogenerated wrapper script for GAP_pkg_io_jll for x86_64-unknown-freebsd
export io

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_io")
JLLWrappers.@declare_file_product(io)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        io,
        "lib/gap/io.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, JLD, LinearAlgebra
mutable struct problem
    data::Array{Float64}
    Sigma::Array{Float64}
end
###############################################################################
### This test script aims to verify correctness of our code,
### by verifying that the sparse principal components found here match those
### found by (Berk+Bertsimas, 2017) on the UCI dataset
###
### We have hard-coded the test data here, in order to avoid Julia dependency issues
###############################################################################
include(normpath(joinpath(@__FILE__,"..",".."))*"Algorithm/branchAndBound.jl")
include(normpath(joinpath(@__FILE__,"..",".."))*"Algorithm/utilities.jl")
include(normpath(joinpath(@__FILE__,"..",".."))*"Algorithm/multiComponents.jl")
################################################################################
######## Test pitprops data set for each cardinality level
################################################################################
@testset "Berk and Bertsimas SPCA Table 8" begin

pitprops=[[1,0.954,0.364,0.342,-0.129,0.313,0.496,0.424,0.592,0.545,0.084,-0.019,0.134];
       [0.954,1,0.297,0.284,-0.118,0.291,0.503,0.419,0.648,0.569,0.076,-0.036,0.144];
       [0.364,0.297,1,0.882,-0.148,0.153,-0.029,-0.054,0.125,-0.081,0.162,0.22,0.126];
       [0.342,0.284,0.882,1,0.22,0.381,0.174,-0.059,0.137,-0.014,0.097,0.169,0.015];
       [-0.129,-0.118,-0.148,0.22,1,0.364,0.296,0.004,-0.039,0.037,-0.091,-0.145,-0.208];
       [0.313,0.291,0.153,0.381,0.364,1,0.813,0.09,0.211,0.274,-0.036,0.024,-0.329];
       [0.496,0.503,-0.029,0.174,0.296,0.813,1,0.372,0.465,0.679,-0.113,-0.232,-0.424];
       [0.424,0.419,-0.054,-0.059,0.004,0.09,0.372,1,0.482,0.557,0.061,-0.357,-0.202];
       [0.592,0.648,0.125,0.137,-0.039,0.211,0.465,0.482,1,0.526,0.085,-0.127,-0.076];
       [0.545,0.569,-0.081,-0.014,0.037,0.274,0.679,0.557,0.526,1,-0.319,-0.368,-0.291];
       [0.084,0.076,0.162,0.097,-0.091,-0.036,-0.113,0.061,0.085,-0.319,1,0.029,0.007];
       [-0.019,-0.036,0.22,0.169,-0.145,0.024,-0.232,-0.357,-0.127,-0.368,0.029,1,0.184];
       [0.134,0.144,0.126,0.015,-0.208,-0.329,-0.424,-0.202,-0.076,-0.291,0.007,0.184,1];]
spca_pitrops_value=[1.000, 1.954, 2.475, 2.937, 3.406, 3.771, 3.996, 4.069, 4.139, 4.173, 4.208, 4.218, 4.219]#source: Berk+B paper (unpublished manuscript)
pitprops=reshape(pitprops, (13,13));
B=sqrt(pitprops);
theProb=problem(B, pitprops)
    for k=1:13
           println("Running pitprops data set k= ", k)
           obj, xVal, timetoBound, timetoConverge, timeOut, explored, toPrint=branchAndBound(theProb, k)
           println("Testing pitprops data set k= ", k)
           @show @test(abs(obj-spca_pitrops_value[k])<1e-3)
    end

end;

max_stars_repo_name: ryancorywright/Optimal-SPCA

max_issues_repo_name: ryancorywright/Optimal-SPCA

max_forks_repo_name: ryancorywright/Optimal-SPCA


text:
using Surrogates
using LinearAlgebra
using Flux
using Flux: @epochs
using PolyChaos

#######SRBF############
##### 1D #####

lb = 0.0
ub = 15.0
objective_function = x -> 2*x+1
x = [2.5,4.0,6.0]
y = [6.0,9.0,13.0]

# In 1D values of p closer to 2 make the det(R) closer and closer to 0,
#this does not happen in higher dimensions because p would be a vector and not
#all components are generally C^inf
p = 1.99
a = 2
b = 6

#Using Kriging
my_k_SRBF1 = Kriging(x,y,lb,ub)
surrogate_optimize(objective_function,SRBF(),a,b,my_k_SRBF1,UniformSample())

#Using RadialBasis
my_rad_SRBF1 = RadialBasis(x,y,a,b,rad = linearRadial())
surrogate_optimize(objective_function,SRBF(),a,b,my_rad_SRBF1,UniformSample())

my_wend_1d = Wendland(x,y,lb,ub)
surrogate_optimize(objective_function,SRBF(),a,b,my_wend_1d,UniformSample())

my_earth1d = EarthSurrogate(x,y,lb,ub)
surrogate_optimize(objective_function,SRBF(),a,b,my_earth1d,LowDiscrepancySample(2))

##### ND #####
objective_function_ND = z -> 3*norm(z)+1
lb = [1.0,1.0]
ub = [6.0,6.0]
x = sample(5,lb,ub,SobolSample())
y = objective_function_ND.(x)
p = [1.5,1.5]
theta = [1.0,1.0]

#Kriging

my_k_SRBFN = Kriging(x,y,lb,ub)
#Every optimization method now returns the y_min and its position
x_min, y_min = surrogate_optimize(objective_function_ND,SRBF(),lb,ub,my_k_SRBFN,UniformSample())

#Radials
lb = [1.0,1.0]
ub = [6.0,6.0]
x = sample(5,lb,ub,SobolSample())
objective_function_ND = z -> 3*norm(z)+1
y = objective_function_ND.(x)
my_rad_SRBFN = RadialBasis(x,y,lb,ub,rad = linearRadial())
surrogate_optimize(objective_function_ND,SRBF(),lb,ub,my_rad_SRBFN,UniformSample())

# Lobachevsky
x = sample(5,lb,ub,UniformSample())
y = objective_function_ND.(x)
alpha = [2.0,2.0]
n = 4
my_loba_ND = LobachevskySurrogate(x,y,lb,ub)
surrogate_optimize(objective_function_ND,SRBF(),lb,ub,my_loba_ND,UniformSample())
#Linear
lb = [1.0,1.0]
ub = [6.0,6.0]
x = sample(500,lb,ub,SobolSample())
objective_function_ND = z -> 3*norm(z)+1
y = objective_function_ND.(x)
my_linear_ND = LinearSurrogate(x,y,lb,ub)
surrogate_optimize(objective_function_ND,SRBF(),lb,ub,my_linear_ND,SobolSample(),maxiters=15)

#=
#SVM
lb = [1.0,1.0]
ub = [6.0,6.0]
x = sample(5,lb,ub,SobolSample())
objective_function_ND = z -> 3*norm(z)+1
y = objective_function_ND.(x)
my_SVM_ND = SVMSurrogate(x,y,lb,ub)
surrogate_optimize(objective_function_ND,SRBF(),lb,ub,my_SVM_ND,SobolSample(),maxiters=15)
=#

#Inverse distance surrogate
lb = [1.0,1.0]
ub = [6.0,6.0]
x = sample(5,lb,ub,SobolSample())
objective_function_ND = z -> 3*norm(z)+1
my_p = 2.5
y = objective_function_ND.(x)
my_inverse_ND = InverseDistanceSurrogate(x,y,lb,ub,p=my_p)
surrogate_optimize(objective_function_ND,SRBF(),lb,ub,my_inverse_ND,SobolSample(),maxiters=15)

#SecondOrderPolynomialSurrogate
lb = [0.0,0.0]
ub = [10.0,10.0]
obj_ND = x -> log(x[1])*exp(x[2])
x = sample(15,lb,ub,UniformSample())
y = obj_ND.(x)
my_second_order_poly_ND = SecondOrderPolynomialSurrogate(x,y,lb,ub)
surrogate_optimize(obj_ND,SRBF(),lb,ub,my_second_order_poly_ND,SobolSample(),maxiters=15)

####### LCBS #########
######1D######
objective_function = x -> 2*x+1
lb = 0.0
ub = 15.0
x = [2.0,4.0,6.0]
y = [5.0,9.0,13.0]
p = 1.8
a = 2.0
b = 6.0
my_k_LCBS1 = Kriging(x,y,lb,ub)
surrogate_optimize(objective_function,LCBS(),a,b,my_k_LCBS1,UniformSample())

##### ND #####
objective_function_ND = z -> 3*norm(z)+1
x = [(1.2,3.0),(3.0,3.5),(5.2,5.7)]
y = objective_function_ND.(x)
p = [1.2,1.2]
theta = [2.0,2.0]
lb = [1.0,1.0]
ub = [6.0,6.0]

#Kriging
my_k_LCBSN = Kriging(x,y,lb,ub)
surrogate_optimize(objective_function_ND,LCBS(),lb,ub,my_k_LCBSN,UniformSample())

##### EI ######

###1D###
objective_function = x -> (x+1)^2 - x + 2# Minimum of this function is at x = -0.5, y = -2.75
true_min_x = -0.5
true_min_y = objective_function(true_min_x)
lb = -5
ub = 5
x = sample(5, lb, ub, SobolSample())
y = objective_function.(x)
my_k_EI1 = Kriging(x,y,lb,ub; p = 2)
surrogate_optimize(objective_function,EI(),lb,ub,my_k_EI1, SobolSample(),maxiters=200,num_new_samples=155)

# Check that EI is correctly minimizing the objective
y_min, index_min = findmin(my_k_EI1.y)
x_min = my_k_EI1.x[index_min]
@test norm(x_min - true_min_x) < 0.05 * norm(ub .- lb)
@test abs(y_min - true_min_y) < 0.05 * (objective_function(ub) - objective_function(lb))

###ND###
objective_function_ND = z -> 3*norm(z)+1 # this is minimized at x = (0, 0), y = 1
true_min_x = (0.0, 0.0)
true_min_y = objective_function_ND(true_min_x)
x = [(1.2,3.0),(3.0,3.5),(5.2,5.7)]
y = objective_function_ND.(x)
min_y = minimum(y)
p = [1.2,1.2]
theta = [2.0,2.0]
lb = [-1.0,-1.0]
ub = [6.0,6.0]

#Kriging
my_k_EIN = Kriging(x,y,lb,ub)
surrogate_optimize(objective_function_ND,EI(),lb,ub,my_k_EIN,SobolSample())

# Check that EI is correctly minimizing instead of maximizing
y_min, index_min = findmin(my_k_EIN.y)
x_min = my_k_EIN.x[index_min]
@test norm(x_min .- true_min_x) < 0.05 * norm(ub .- lb)
@test abs(y_min - true_min_y) < 0.05 * (objective_function_ND(ub) - objective_function_ND(lb))

###ND with SectionSampler###
# We will make sure the EI function finds the minimum when constrained to a specific slice of 3D space
objective_function_section = x -> x[1]^2 + x[2]^2 + x[3]^2 # this is minimized at x = (0, 0, 0), y = 0

# We will constrain x[2] to some value
x2_constraint = 2.0
true_min_x = (0.0, x2_constraint, 0.0)
true_min_y = objective_function_section(true_min_x)

sampler = SectionSample([NaN64, x2_constraint, NaN64], SobolSample())
lb = [-1.0,x2_constraint,-1.0]
ub = [6.0,x2_constraint,6.0]
x = sample(5, lb, ub, sampler)
y = objective_function_section.(x)
p = [1.2,1.2,1.2]
theta = [2.0,2.0,2.0]

#Kriging
my_k_EIN_section = Kriging(x,y,lb,ub; p = p)
# Constrain our sampling to the plane where x[2] = 1
surrogate_optimize(objective_function_section,EI(),lb,ub,my_k_EIN_section, sampler)

# Check that EI is correctly minimizing instead of maximizing
y_min, index_min = findmin(my_k_EIN_section.y)
x_min = my_k_EIN_section.x[index_min]
@test norm(x_min .- true_min_x) < 0.05 * norm(ub .- lb)
@test abs(y_min - true_min_y) < 0.05 * (objective_function_section(ub) - objective_function_section(lb))

## DYCORS ##

#1D#
objective_function = x -> 3*x+1
x = [2.1,2.5,4.0,6.0]
y = objective_function.(x)
p = 1.9
lb = 2.0
ub = 6.0

my_k_DYCORS1 = Kriging(x,y,lb,ub,p=1.9)
surrogate_optimize(objective_function,DYCORS(),lb,ub,my_k_DYCORS1,UniformSample())

my_rad_DYCORS1 = RadialBasis(x,y,lb,ub,rad = linearRadial())
surrogate_optimize(objective_function,DYCORS(),lb,ub,my_rad_DYCORS1,UniformSample())

#ND#
objective_function_ND = z -> 2*norm(z)+1
x = [(2.3,2.2),(1.4,1.5)]
y = objective_function_ND.(x)
p = [1.5,1.5]
theta = [2.0,2.0]
lb = [1.0,1.0]
ub = [6.0,6.0]

my_k_DYCORSN = Kriging(x,y,lb,ub)
surrogate_optimize(objective_function_ND,DYCORS(),lb,ub,my_k_DYCORSN,UniformSample(),maxiters=30)

my_rad_DYCORSN = RadialBasis(x,y,lb,ub,rad = linearRadial())
surrogate_optimize(objective_function_ND,DYCORS(),lb,ub,my_rad_DYCORSN,UniformSample(),maxiters=30)

my_wend_ND = Wendland(x,y,lb,ub)
surrogate_optimize(objective_function_ND,DYCORS(),lb,ub,my_wend_ND,UniformSample(),maxiters=30)

### SOP ###
# 1D
objective_function = x -> 3*x+1
x = sample(20,1.0,6.0,SobolSample())
y = objective_function.(x)
p = 1.9
lb = 1.0
ub = 6.0
num_centers = 2
my_k_SOP1 = Kriging(x,y,lb,ub,p=1.9)
surrogate_optimize(objective_function,SOP(num_centers),lb,ub,my_k_SOP1,SobolSample(),maxiters=60)
#ND
objective_function_ND = z -> 2*norm(z)+1
x = [(2.3,2.2),(1.4,1.5)]
y = objective_function_ND.(x)
p = [1.5,1.5]
theta = [2.0,2.0]
lb = [1.0,1.0]
ub = [6.0,6.0]
my_k_SOPND = Kriging(x,y,lb,ub)
num_centers = 2
surrogate_optimize(objective_function_ND,SOP(num_centers),lb,ub,my_k_SOPND,SobolSample(),maxiters=20)

#multi optimization
#=
f  = x -> [x^2, x]
lb = 1.0
ub = 10.0
x  = sample(100, lb, ub, SobolSample())
y  = f.(x)
my_radial_basis_smb = RadialBasis(x, y, lb, ub, rad = linearRadial())
surrogate_optimize(f,SMB(),lb,ub,my_radial_basis_ego,SobolSample())

f  = x -> [x^2, x]
lb = 1.0
ub = 10.0
x  = sample(100, lb, ub, SobolSample())
y  = f.(x)
my_radial_basis_rtea = RadialBasis(x, y, lb, ub, rad = linearRadial())
Z = 0.8 #percentage
K = 2 #number of revaluations
p_cross = 0.5 #crossing vs copy
n_c = 1.0 # hyperparameter for children creation
sigma = 1.5 # mutation
surrogate_optimize(f,RTEA(Z,K,p_cross,n_c,sigma),lb,ub,my_radial_basis_rtea,SobolSample())
=#

max_stars_repo_name: ludoro/Surrogates

max_issues_repo_name: ludoro/Surrogates

max_forks_repo_name: ludoro/Surrogates


text:
load("utils/req.jl")
req("prettyshow/prettyshow.jl")
req("dag/transforms.jl")

# -- node type specific stuff -------------------------------------------------

pshow_nodeval(io::PrettyIO, node::NoNode) = pprint(io, "NoNode()")
pshow_nodeval(io::PrettyIO, node::LiteralNode) = pprint(io,
                                             "LiteralNode($(node.val.value))")
pshow_nodeval(io::PrettyIO, node::SymNode) = pprint(io,
                              "SymNode(:$(node.val.name), :$(node.val.kind))")

repargname(argname, n) = (n == 1 ? [argname] : [argname*string(k) for k=1:n])

get_signature(node::CallNode    ) = ("CallNode"    , [".op", ".arg"])
get_signature(node::RefNode     ) = ("RefNode"     , [".ref",  ".ind"])
get_signature(node::TupleNode   ) = ("TupleNode"   , [".arg"])
get_signature(node::KnotNode    ) = ("KnotNode"    , [
                           repargname(".pre",length(node.args)-1), ".value"])
get_signature(node::EllipsisNode) = ("EllipsisNode", [".arg"])
get_signature(node::AssignNode  ) = ("AssignNode"  , [".lhs", ".rhs", ".dep"])

pshow_nodeval(io::PrettyIO, node::Node) = pprint(io, "Node(", node.val, ")")

# -- pshow -------------------------------------------------------------------

function pshow(io::PrettyIO, sink::Node)
    numassigns::Int = 0
    rewrite(node::Node, args::Vector) = Node(node, args)
    function rewrite(node::ActionNode, args::Vector)
        node = Node(node, args)
        node.name = symbol("assign"*string(numassigns+=1))
        node
    end
    sink = rewrite_dag(sink, rewrite)
    sink = name_fanout_nodes(sink)
    sink = name_nodes_uniquely(sink)

    firstnode = true
    for node in forward(sink)
        if has_name(node)
            pprintln(io)
            pshow_tree(io, node)
            firstnode = false
        end
    end
    if !has_name(sink)
        pprintln(io)
        pshow_tree(io, sink)
    end
end

function pshow_tree(io::PrettyIO, node::Node)
    if has_name(node)
        pprint(io, get_name(node), "=")
    end
    if isa(node, TerminalNode)
        pshow_nodeval(io, node)
    else
        name, argnames = get_signature(node)
        numargs = length(node.args)

        dlength = numargs - length(argnames)
        argnames = [argnames[1:end-1], repargname(argnames[end], 1+dlength)]

        new_args_on_next_line::Int = 0
        on_last_arg::Bool=true
        on_last_node_line::Bool=true
        function newline_hook()
            nal, new_args_on_next_line = new_args_on_next_line, 0
#             split  = (nal >= 2 ? "#" : "+")
#             split  = on_last_arg ? " " : ((nal >= 2 ? "#" : "+"))
            split  = on_last_node_line ? "\\" : (nal >= 2 ? "#" : "+")
#             extend = (nal >= 2 ? "=" : (on_last_arg ? "\\" : "-"))
            extend = (nal >= 2 ? "=" : "-")
            if nal > 0
                return " "*split*extend*" "
            end
            return on_last_arg ? "    " : " |  "
        end

        verbs = [!has_name(arg) for arg in node.args]

        pprint(io, name, "(")
        let io=PrettyChild(io, newline_hook)
            lastverbose = true
            for ((arg, argname), k) in enumerate(zip(node.args, argnames))
                on_last_arg = k==numargs
                # argindex = k
                verbose = verbs[k]
                on_last_node_line = on_last_arg || !any(verbs[k:end])
                if lastverbose || verbose
                    if on_last_arg || verbose || verbs[k+1]
                        new_args_on_next_line = 1
                    else
                        new_args_on_next_line = 2
                    end
                    pprint(io, '\n')
                end
                pprint(io, argname, "=")
                if verbose;  pshow_tree(io, arg);
                else         pshow(io, get_name(arg));  end
                if !on_last_arg
                    pprint(io, ", ")
                end
                lastverbose = verbose
            end
        end
        pprint(io, ")")
    end
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
struct Trade
    price::Int64
    quantity::Int64
    timestamp::Int64

    from::AbstractAsset
    to::AbstractAsset

    seller::Union{String, AbstractInvestor}
    buyer::Union{String, AbstractInvestor}
end

"Trade limit orders"
function trade!(buy::BidLimitOrder, sell::AskLimitOrder; price::Int64, from::AbstractAsset, to::AbstractAsset)

    quantity = get_clearable_quantity(buy, sell)
    # Both order's quantity is reduced by trade quantity
    buy.quantity -= quantity
    sell.quantity -= quantity
    #println("Trading $quantity with $(buy.dealer.name) ($(from.name)) <--> $(sell.dealer.name) ($(to.name)) with price $price ")
    # Must be released before transfer
    # as the reserve cannot be higher
    # than the amount at any given time
    release!(buy.dealer, quantity * buy.price, from) # We must reduce buy price for reserve accounting
    release!(sell.dealer, quantity, to)

    #println("Trading $price * $quantity ($from --> $to)")
    # TODO: How to specify which currency? (asset to seller)
    transfer!(buy.dealer, sell.dealer, asset=from, quantity=quantity * price)
    transfer!(sell.dealer, buy.dealer, asset=to, quantity=quantity)

    return Trade(price, quantity, get_time(), from, to, sell.dealer, buy.dealer)
end

# Allow argument order insensitive
trade!(sell::AskLimitOrder, buy::BidLimitOrder; side=:equal, asset::Union{Nothing, AbstractAsset}=nothing) = trade!(buy, sell, asset, side=side, asset=asset)


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 54e6bd88-a6a9-11eb-3380-49b79430f6cf
html"""

Section 3.2

 ODEs and parameterized types

Lecture Video

"""

# ╔═╡ 5997bb40-b508-4647-8b12-a7cc5152e550
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="DifferentialEquations", version="6"),
        Pkg.PackageSpec(name="Plots", version="1"),
		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="MappedArrays", version="0.4"),

    ])

    using DifferentialEquations, Plots, PlutoUI, MappedArrays
end

# ╔═╡ 38c74e63-13e9-49a9-8eae-4a898197647a
TableOfContents(aside=true)

# ╔═╡ 7a1d476b-1302-4b98-b9fa-01b990979985
md"""
# Julia concepts
- Using libraries
- Libraries are built around types and objects
- Plot recipes: Defining how to plot our objects
- Making objects into callable functions
- Parameterized types
"""

# ╔═╡ 3004336a-0250-483f-be00-60db592c9fc8
md"""
# ODEs and parameterized types
"""

# ╔═╡ a1cfd00b-825f-494e-91a9-192815f06ad6
md"""
In the previous lecture we started to look at Ordinary Differential Equations (ODEs). These are equations relating a function to its derivatives, telling us how to move at each moment in time.

For example, exponential decay is modelled by

$$\dot{u} = -p \, u,$$

in other words

$$\dot{u}(t) = - p \, u(t)$$

for all real $t \ge 0$,

with initial condition

$$u(0) = u_0$$

where $u(t)$ is the number of light bulbs still functioning at time $t$, the number of radioactive nuclei that have not yet decayed, or the number of people who are still infectious and have not yet recovered. The parameter $p$ is the decay rate.

We saw that **time stepping** methods provide a class of numerical methods for solving such initial-value problems, of which the Euler method is the simplest example.

Unfortunately the Euler method is *not* a good choice, since it is not accurate enough. In numerical analysis classes (like [18.330](https://github.com/mitmath/18330)) various better methods are derived.

Fortunately those methods are available in a Julia **library** or package, `DifferentialEquations.jl`.

"""

# ╔═╡ b2d2b8f0-c81d-437a-83cf-726b5caa77ab
md"""
# The `DifferentialEquations.jl` library
"""

# ╔═╡ d27ec502-ea0c-4d15-a223-7cbf6b784d1c
md"""
[DifferentialEquations.jl](https://diffeq.sciml.ai/v2.0/) (often abbreviated as DiffEq) is a Julia library for solving many types of differential equations. It has become the best (most complete and fastest) general library for **solving** ordinary differential equations (and many other types) in any programming environment. It is also a very good example of the way that Julia libraries are often structured.

We will see some basic usage examples and how Julia **types** are used both to control the library and to store the input and output to the **solvers**.
"""

# ╔═╡ 6f87059a-dbdc-46dd-933b-7bb36980868e
md"""
## Basic usage of DifferentialEquations.jl
"""

# ╔═╡ 557cedbc-6e44-40a8-b0ab-37ae11b63f55
md"""
Let's try to solve the above equation that we derived in the previous lecture, $\dot{u} = -p \, u$. (We are using $u$ as the name of the dependent variable for consistency with the `DifferentialEquations.jl` documentation and internal structure.)

We write a general ODE in 1 variable as

$\dot{u} = f_p(t, u).$

The right-hand side of the ODE is given by a function $f$, which may depend on a **parameter** $p$, i.e. a constant which does not change for any given run of the ODE.
$f$ may also depend explicitly on time.

In our example we do indeed have a parameter $p$ that we want to change, but $f$ does not have an explicit time dependence, so we write
"""

# ╔═╡ a6ae905b-db28-4e75-b567-54840dcda9a5
f(u, p, t) = -p * u

# ╔═╡ f16644d0-e2d2-44ac-988c-aa5b178844f5
md"""
Note that you always need the three input arguments `u`, `p` and `t`, in that order, even if you don't use them.
"""

# ╔═╡ 1ea2acfe-6ba1-44d3-9a24-a904a5fac0c4
md"""
Now we define our initial condition $u_0$ and the time span to integrate over, given as a tuple (initial\_time, final\_time):
"""

# ╔═╡ 264a27cd-9ace-48cc-9556-7eb798298f17
begin
	u0 = 100.0

	time_span = (0.0, 10.0)
end

# ╔═╡ 70ebe76d-6107-492d-a3d2-84c1dd7f59fb
md"""
To set up the **problem** instance we use a type `ODEProblem` that is defined in the `DifferentialEquations.jl` package, into which we pass all the information necessary to define the problem. The parameters must go in the following order:
"""

# ╔═╡ 6e7c8e9e-bef0-4bf8-b1b1-d50c82aa203e
problem = ODEProblem(f, u0, time_span, p)

# ╔═╡ b5379dc2-d97f-47ed-8737-35e3fe59285c
md"""
(For more advanced use there are also some additional, optional, keyword arguments.)
Note that the displayed output does not currently include information about the function, nor the parameters.
"""

# ╔═╡ 0ed8e007-2137-43e0-93b8-6f5dfb9496f5
md"""
To solve the ODE we call the `solve` function:
"""

# ╔═╡ ff91515a-89a4-4423-a186-5572c712493d
solution = solve(problem)

# ╔═╡ 64b26404-1ca0-4aa4-bea0-3c9075b08298
md"""
What happened here? A suitable solver (i.e. an algorithm to calculate the solution) was chosen *automatically*, and it chose certain moments in time at which to output information about the (approximate, but very accurate) solution.

In this particular case it chose to output data at only eight points in time between $t=0$ and $t=10$.

Let's try to plot the `solution` object:
"""

# ╔═╡ b07de2ce-640b-42d6-8b60-39fcbcd116e7
plot(solution, size=(500, 300), label="solution")

# ╔═╡ 0b1855de-a8c9-492e-a249-3238e41fe84c
md"""
### Plot recipes
"""

# ╔═╡ 0190bdb0-97bf-4c35-b3c8-e97c1940b702
md"""
Two surprising things happen. Firstly, there is no reason to expect this to have worked: `solution` is some kind of Julia object, but somehow `Plots.jl` knows how to plot it. This is because `DifferentialEquations.jl` defines a **plot recipe**. This specifies a way to turn solution objects into plots! Any package can do this relatively easily.
"""

# ╔═╡ 555a6328-8aff-4f9f-87f5-51a67450002d
md"""
The second surprise is that the output looks like a smooth curve, rather than just 8 points. Let's see those points on top of the curve. We can extract the relevant data from the `solution` object:
"""

# ╔═╡ d65372df-e6ca-4165-8d9b-b21cc5c9f796
scatter!(solution.t, solution.u, label="discrete output")

# ╔═╡ b274c61b-4fab-4f4a-99df-5d39e0f56aa1
md"""
We see that the package in fact gives not only the value at those points, but it is in fact also capable of calculating an (approximate) solution at *any* intermediate point, using **interpolation**. In fact, we can access this by treating `solution` as if it were a function:
"""

# ╔═╡ e3d9e0be-6cf4-4ae1-8c7f-b444a153a9f5
begin
	tt = 3.5
	solution(tt)
end

# ╔═╡ 180bcdc1-5f51-4b32-8b9c-5000605cdf32
scatter!([tt], [solution(tt)], label="t = $(tt)", ms=5, m=:square)

# ╔═╡ 13b38f88-2ead-46b3-bc96-eae2ea10204d
md"""
For this particular ODE we know the analytical solution. Let's compare them as we vary the parameter $p$:
"""

# ╔═╡ 475ef0e5-7f3a-4bbb-b513-8be29abfd3f9
md"""
p = $(@bind p Slider(0.0:0.1:2.0, show_value=true))
"""

# ╔═╡ 0ed6d65c-707d-4666-b203-2ad0ea822687
let

	problem = ODEProblem(f, u0, time_span, p)
	solution = solve(problem)

	plot(solution,
		linewidth=3, xlabel="t", yaxis="x(t)", label="numerical", size=(500, 300))

	plot!(t -> u0 * exp(-p*t), lw=3, ls=:dash, label="exact")

	ylims!(0, 100)
	title!("p = $p")
end

# ╔═╡ cf9a293c-aa0f-40ce-960b-d4c3b82b8346
md"""
We see that the numerical and exact solutions are (to the eye) indistinguishable, and that the package is fast enough to calculate the solution basically in real time.
"""

# ╔═╡ efcfdc1a-5a80-4780-ba56-43c5e2d6ea36
md"""
## Systems of ODEs
"""

# ╔═╡ 5d7a81d7-8080-4390-924a-72c34c0a5e23
md"""
Now let's try to solve the SIR equations:

$$\begin{aligned}
\dot{s} &= -\beta s i \\
\dot{i} &= +\beta s i - \gamma i\\
\dot{r} &= + \gamma i
\end{aligned}$$
"""

# ╔═╡ 6c607126-dd0b-438d-8724-c1ef5e79c71d
md"""
We need to think of the system in a vector form,

$$\dot{\mathbf{x}} = \mathbf{f}_\mathbf{p}(t, \mathbf{x}),$$

where:
- the vector $\mathbf{x}$ is a vector of all the variables, $\mathbf{x} := (x_1, x_2, x_3) := (s, i, r)$
- the function $\mathbf{f}$ is a vector-valued function $\mathbb{R}^3 \to \mathbb{R}^3$, where $f_k$ gives the right-hand side of the equation for $\dot{x_k}$; and
- the vector $\mathbf{p} = (\beta, \gamma)$ is a vector of the parameters.

"""

# ╔═╡ d16e06f4-55c4-4c0e-bde5-c0ce4b36b043
x0 = [0.99, 0.01, 0.0]

# ╔═╡ aa0e4260-6268-4eff-a662-f3fd55e7d450
function SIR(x, p, t)

	s, i, r = x    # unpack the vectors into scalar values
	β, γ = p

	# build a new vector to return:

	return [-β * s * i,
			+β * s * i - γ * i,
					   + γ * i]

end

# ╔═╡ d0f40681-73df-4cd3-bbd5-3edb8193153e
params = [β, γ]

# ╔═╡ fd122a25-5655-410b-aa7f-34936fc97b53
SIR_problem = ODEProblem(SIR, x0, (0.0, 50.0), params)

# ╔═╡ 8ac7da6b-46e0-470b-91b3-1a61c226fa4a
sol = solve(SIR_problem)

# ╔═╡ 7aa45efe-865d-4e0e-9c71-0c032c72c40d
md"""
Now we see that the solverr has recognised that everything is a vector, and it returns a vector at each time stamp.

Again we can plot:
"""

# ╔═╡ bca8112f-cb61-43dc-ae87-383915c8a89b
gr()

# ╔═╡ a766a141-5d7b-499f-9e18-48bf926ee7ea
plot(sol)

# ╔═╡ 76cbc37d-54c8-4626-8bfe-58b63a602c38
md"""
β = $(@bind β Slider(-0.5:0.01:2.0, default=1.0, show_value=true))

γ = $(@bind γ Slider(-0.5:0.01:2.0, default=0.1, show_value=true))
"""

# ╔═╡ c2765282-bcc7-4110-9822-10557326461e
md"""
It knows to plot each variable separately.

We can instead plot combinations of variables in *phase space* or *state space*:
"""

# ╔═╡ 7671506d-d4b7-4792-9571-e003097235e1
gr()

# ╔═╡ 203e90b0-4d0c-4999-8513-d4eb43f53aac
plot(sol, vars=(1, 2), xlabel="s", ylabel="i", arrow=true, xlims=(-0.1, 1.0), size=(500, 300))

# ╔═╡ 326890ae-0675-4779-b5e8-9b3e0412a52b
md"""
And even in 3D:
"""

# ╔═╡ d7025d1f-c3b8-461a-94dc-2414fbdfd373
plotly()

# ╔═╡ 7f5c41f9-96b9-40d6-87f8-3e81c372b48e
plot(sol, vars=(1, 2, 3), xlabel="s", ylabel="i", zlabel="r")

# ╔═╡ d8f15e9b-aef1-4795-8522-85ea3564e351
md"""
The [ModelingToolkit.jl](https://mtk.sciml.ai/stable/tutorials/ode_modeling/) library  provides ways to make creating ODE models more intuitive, using symbolic equation objects.
"""

# ╔═╡ e723400e-3f8f-47d6-8944-28355a54c698
md"""
## Making objects into callable functions
"""

# ╔═╡ 15fdb430-4157-4401-85b9-665dde795a22
md"""
One of the interesting things we have seen is the way that the `solution` object can be **called** using the normal syntax for calling functions:
"""

# ╔═╡ a410e98d-dd7f-4f55-8812-8b0339a56c86
solution(3.5)

# ╔═╡ 239f39b2-53c6-4da9-ac25-b5e7d19deb5f
md"""
Not all types allow this. For example, if we define a matrix
"""

# ╔═╡ 6630cbc8-1b2f-4ea9-bcd6-e4f4fcb2e844
A = rand(2, 2)

# ╔═╡ 5c1c7016-69aa-4009-8dc8-0c68c925ee66
md"""
we cannot usually call it with a vector:
"""

# ╔═╡ 83719ef1-e8d6-44d0-98ab-a851b1082fa5
A(rand(2))

# ╔═╡ e4d36db5-352b-4fb2-9169-da5c6a3458b7
md"""
If we wanted to do so, we could, however, using the following syntax:
"""

# ╔═╡ 167015ff-1caf-42a8-85f1-070da3e5d5c5
# (A::Matrix)(x) = A * x

# ╔═╡ f8bcd82d-089c-49a2-99d0-b17a5ac5c509
md"""
Note that this makes objects (instances) of the type callable, rather than the type itself, which we can already call as a **constructor** to create objects of that type:
"""

# ╔═╡ cd17e738-089b-4934-af63-d27706c5d1c7
Matrix([1 2; 3 4])

# ╔═╡ fc525f77-ad62-43cb-92e4-0fc1354be99c
md"""
For example, if we wanted to write our own type to store the output of the Euler method, we could do something like
"""

# ╔═╡ ec90842a-a20c-4847-a803-e070d343b98d
begin
	struct SimpleEulerOutput
		times::Vector{Float64}
		values::Vector{Float64}
	end

	function (soln::SimpleEulerOutput)(t)
		# use interpolation to calculate the value at time t
		# using  `soln.times` and `soln.values`
	end
end

# ╔═╡ 9fad6f01-a3b7-4dbb-a8f5-29f2ea1f545a
md"""
The definition of the function makes it callable.

This says that "the way to call an object called `soln`, of type `EulerOutput`, is the following:". Note that inside the function you will need too access the `times` and `values` from *inside* the `soln` object, using `soln.times` and `soln.values`.

Note that at the time of writing, Pluto requires that you put this call definition in the same cell as the type definition.
"""

# ╔═╡ 658d8d15-5b3e-4d35-88b1-4e2fa0f3b17a
md"""
## Digging into objects
"""

# ╔═╡ aafeeb55-943f-4375-9cc0-e0dce659d161
md"""
A key feature of the `DifferentialEquations.jl` library, and of most other libraries in Julia, is the way they define and use new **types**.

As we have seen previously, we can define new types using `struct` which store data, and then define functions which act on those types.
"""

# ╔═╡ 3c45ad21-d18b-41fa-a0fd-64d8988621f3
md"""
Let's start by looking at the `ODEProblem` type:
"""

# ╔═╡ 96e871a8-48b9-459b-a1a0-94fc1deb941e
problem

# ╔═╡ 2622801c-4e37-42f1-ab77-939d461b8d20
md"""
To see what fields, or attributes, the object contains we can use
"""

# ╔═╡ 43e7b270-6452-4a4c-ba6f-c286a337a423
fieldnames(typeof(problem))

# ╔═╡ dbdb9998-ca3b-46fb-87cf-5b7c6e387ab9
md"""
In Pluto and other interactive environments we can get the same information using `problem.<TAB>`, i.e. by typing the TAB key after `problem.`
"""

# ╔═╡ b8fbff51-72fe-4432-9f6f-06f4628320d4
md"""
As usual we can extract the data contained in the object as usual using `.`:
"""

# ╔═╡ ef677abe-c96d-41cc-9b4a-192ff29cdb86
problem.u0

# ╔═╡ 78933b33-8a3c-492e-8a58-19ad7187b45e
problem.tspan

# ╔═╡ 9ef4aa78-2ecd-4d53-87ae-f0909bb46915
md"""
To see everything contained in the object, we can use `Dump` in Pluto, or `dump` if we are not using Pluto:
"""

# ╔═╡ 798bc482-7556-4904-8946-0d9898ce2d33
Dump(problem)

# ╔═╡ 32e9f6e5-8df6-4d40-8092-efaf2bce28d3
md"""
This seems (much) more complicated than you might expect. The library is clearly doing more than just storing the data that we provide in a type! The main thing that has happened is that the function we passed in has been processed by wrapping it into yet another type, `ODEFunction`. This contains information such as how to calculate the derivative of the function that can be provided for more advanced usage of the package.
"""

# ╔═╡ 6dfc2699-fe5e-4b55-8437-0866f715170a
md"""
Similarly we can look inside the solution object:
"""

# ╔═╡ 17d6021b-0d21-4d04-a8ae-64743cc01404
fieldnames(typeof(solution))

# ╔═╡ f4c3b174-57ac-461b-a816-869460d8896d
Dump(solution)

# ╔═╡ 500973f5-faa7-4a31-b812-d0b20dbb9d82
md"""
The solution is even more complicated, containing not only the data that was calculated, but also all information about which algorithms were used to solve the problem and tables of coefficients for interpolation.
"""

# ╔═╡ 4a575627-32c2-4842-a540-52253f79ff89
md"""
## Parameterized types
"""

# ╔═╡ bea223ef-f961-457c-8962-d0ed396f4717
md"""
Let's look now at the *type* of the DiffEq objects:
"""

# ╔═╡ 92bd82dc-eaab-4869-9d95-cb4647092352
typeof(problem)

# ╔═╡ d518828a-1972-4bfa-b96f-0ebbf86b6193
md"""
We see that indeed `problem` is an object of type `ODEProblem`, but that it also has several **type parameters**, which are listed inside the curly braces (`{` and `}`):
"""

# ╔═╡ c3570557-55bd-4189-a332-aa2e54693665
[typeof(problem).parameters...]

# ╔═╡ 1f684b1b-8e43-4e0e-8a47-1f2474cc5e4f
md"""
For example, we see that the second type parameter is of type `Tuple{Float64, Float64}`, which is the type of the variable `time_span`:
"""

# ╔═╡ 1bf84cfe-9ddc-4644-b7ea-8ac43b01d477
typeof(time_span)

# ╔═╡ d2a1f97e-2e3e-4ae9-8287-2de62065d708
md"""
One of the reasons for using parameterized types like this is for efficiency: Julia is fastest when it knows the exact types of every variable, since then it can generate efficient machine code.
"""

# ╔═╡ a71aaf0b-8ce9-4242-b31f-8f05730ce2e7
md"""
# Making parameterized types
"""

# ╔═╡ 09ddc00b-2361-44c5-b194-7b4b6d24bef4
md"""
Let's see how and why we would make a parameterized type.

Let's think again about the `EulerOutput` type that we started to write above.
What happens if, as in the `DifferentialEquations.jl` library, we would like to handle both *scalar* ODEs (one single real dependent variable), and *vector* ODEs?

We can no longer specify that the output type is a `Vector{Float64}`, so we would need to leave it unspecified:
"""

# ╔═╡ c7c16e89-3ccf-4d70-9cdb-d5ae71a5ea49
struct SimpleEulerOutput2
	times
	values
end

# ╔═╡ 1872c30d-aece-44a5-b3b6-36e82c072609
md"""
However, it turns out that we then lose efficiency (speed).

What we would need is two different types: one in which we specify the type of `times` as `Vector{Float64}`, and the other as `Vector{Vector{Float64}}`.

Julia provides a mechanism to define *both at once*: this is what parameterized types do!
"""

# ╔═╡ 5d06734e-ffbd-421e-820e-609e373ccc7b
md"""
Let's define our own `MyODEProblem` type as a simple example. We specify type parameters for each variable:
"""

# ╔═╡ 62d55814-c863-4d4d-8582-908721d37326
struct MyODEProblem{T, U}
	t0::T
	tfinal::T

	u0::U
end

# ╔═╡ 470b7808-5ff6-4cda-82cc-2ac250d51909
md"""
Note that the starting and final times have been specified to require the *same* type `T`.
"""

# ╔═╡ d392de63-137e-43fb-86ca-4f344769bc71
md"""
Now if we create objects of this type, we see that Julia matches the types we pass in to the correct type parameters:
"""

# ╔═╡ e07830be-3aab-4410-a830-786a30eb162a
newprob = MyODEProblem(3, 4, 100.0)

# ╔═╡ 6d809ff0-20c5-4910-aca0-612496955812
newprob.t0

# ╔═╡ 816a64ed-1f88-43a8-9df7-8070372e84ea
newprob.tfinal

# ╔═╡ cdd6a0dd-e188-4d37-ab6d-12e7fecd3e90
newprob.u0

# ╔═╡ c29d4796-9e52-4234-8325-5bf31099c79d
md"""
If we change the types of the inputs, the type parameters change:
"""

# ╔═╡ 7e5b0073-1bc7-41c7-ae74-c9c1bdf6c6ea
MyODEProblem(3.0, 4.0, 100.0)

# ╔═╡ 67e36f7c-0cb9-4ead-82ad-4c47639802a8
md"""
If we mix types for the first two arguments then we get an error, since Julia can't match to a single type `T`:
"""

# ╔═╡ cdc25238-ed48-42e5-9db4-2367577b215d
MyODEProblem(3.0, 5, 100.0)

# ╔═╡ f66b6078-bab5-48ef-9241-8d5fb0dc78ea
md"""
We can use vectors:
"""

# ╔═╡ 73016af9-5efe-4af2-9e99-ce1a7de6fa96
MyODEProblem(3, 4, [17.0, 18.0])

# ╔═╡ Cell order:
# ╟─54e6bd88-a6a9-11eb-3380-49b79430f6cf
# ╠═5997bb40-b508-4647-8b12-a7cc5152e550
# ╠═38c74e63-13e9-49a9-8eae-4a898197647a
# ╟─7a1d476b-1302-4b98-b9fa-01b990979985
# ╟─3004336a-0250-483f-be00-60db592c9fc8
# ╟─a1cfd00b-825f-494e-91a9-192815f06ad6
# ╟─b2d2b8f0-c81d-437a-83cf-726b5caa77ab
# ╟─d27ec502-ea0c-4d15-a223-7cbf6b784d1c
# ╟─6f87059a-dbdc-46dd-933b-7bb36980868e
# ╟─557cedbc-6e44-40a8-b0ab-37ae11b63f55
# ╠═a6ae905b-db28-4e75-b567-54840dcda9a5
# ╟─f16644d0-e2d2-44ac-988c-aa5b178844f5
# ╟─1ea2acfe-6ba1-44d3-9a24-a904a5fac0c4
# ╠═264a27cd-9ace-48cc-9556-7eb798298f17
# ╟─70ebe76d-6107-492d-a3d2-84c1dd7f59fb
# ╠═6e7c8e9e-bef0-4bf8-b1b1-d50c82aa203e
# ╟─b5379dc2-d97f-47ed-8737-35e3fe59285c
# ╟─0ed8e007-2137-43e0-93b8-6f5dfb9496f5
# ╠═ff91515a-89a4-4423-a186-5572c712493d
# ╟─64b26404-1ca0-4aa4-bea0-3c9075b08298
# ╠═b07de2ce-640b-42d6-8b60-39fcbcd116e7
# ╟─0b1855de-a8c9-492e-a249-3238e41fe84c
# ╟─0190bdb0-97bf-4c35-b3c8-e97c1940b702
# ╟─555a6328-8aff-4f9f-87f5-51a67450002d
# ╠═d65372df-e6ca-4165-8d9b-b21cc5c9f796
# ╟─b274c61b-4fab-4f4a-99df-5d39e0f56aa1
# ╠═e3d9e0be-6cf4-4ae1-8c7f-b444a153a9f5
# ╠═180bcdc1-5f51-4b32-8b9c-5000605cdf32
# ╟─13b38f88-2ead-46b3-bc96-eae2ea10204d
# ╟─475ef0e5-7f3a-4bbb-b513-8be29abfd3f9
# ╠═0ed6d65c-707d-4666-b203-2ad0ea822687
# ╟─cf9a293c-aa0f-40ce-960b-d4c3b82b8346
# ╟─efcfdc1a-5a80-4780-ba56-43c5e2d6ea36
# ╟─5d7a81d7-8080-4390-924a-72c34c0a5e23
# ╟─6c607126-dd0b-438d-8724-c1ef5e79c71d
# ╠═d16e06f4-55c4-4c0e-bde5-c0ce4b36b043
# ╠═aa0e4260-6268-4eff-a662-f3fd55e7d450
# ╠═d0f40681-73df-4cd3-bbd5-3edb8193153e
# ╠═fd122a25-5655-410b-aa7f-34936fc97b53
# ╠═8ac7da6b-46e0-470b-91b3-1a61c226fa4a
# ╟─7aa45efe-865d-4e0e-9c71-0c032c72c40d
# ╠═bca8112f-cb61-43dc-ae87-383915c8a89b
# ╠═a766a141-5d7b-499f-9e18-48bf926ee7ea
# ╟─76cbc37d-54c8-4626-8bfe-58b63a602c38
# ╟─c2765282-bcc7-4110-9822-10557326461e
# ╠═7671506d-d4b7-4792-9571-e003097235e1
# ╠═203e90b0-4d0c-4999-8513-d4eb43f53aac
# ╟─326890ae-0675-4779-b5e8-9b3e0412a52b
# ╠═d7025d1f-c3b8-461a-94dc-2414fbdfd373
# ╠═7f5c41f9-96b9-40d6-87f8-3e81c372b48e
# ╟─d8f15e9b-aef1-4795-8522-85ea3564e351
# ╟─e723400e-3f8f-47d6-8944-28355a54c698
# ╟─15fdb430-4157-4401-85b9-665dde795a22
# ╠═a410e98d-dd7f-4f55-8812-8b0339a56c86
# ╟─239f39b2-53c6-4da9-ac25-b5e7d19deb5f
# ╠═6630cbc8-1b2f-4ea9-bcd6-e4f4fcb2e844
# ╟─5c1c7016-69aa-4009-8dc8-0c68c925ee66
# ╠═83719ef1-e8d6-44d0-98ab-a851b1082fa5
# ╟─e4d36db5-352b-4fb2-9169-da5c6a3458b7
# ╠═167015ff-1caf-42a8-85f1-070da3e5d5c5
# ╟─f8bcd82d-089c-49a2-99d0-b17a5ac5c509
# ╠═cd17e738-089b-4934-af63-d27706c5d1c7
# ╟─fc525f77-ad62-43cb-92e4-0fc1354be99c
# ╠═ec90842a-a20c-4847-a803-e070d343b98d
# ╟─9fad6f01-a3b7-4dbb-a8f5-29f2ea1f545a
# ╟─658d8d15-5b3e-4d35-88b1-4e2fa0f3b17a
# ╟─aafeeb55-943f-4375-9cc0-e0dce659d161
# ╟─3c45ad21-d18b-41fa-a0fd-64d8988621f3
# ╠═96e871a8-48b9-459b-a1a0-94fc1deb941e
# ╟─2622801c-4e37-42f1-ab77-939d461b8d20
# ╠═43e7b270-6452-4a4c-ba6f-c286a337a423
# ╟─dbdb9998-ca3b-46fb-87cf-5b7c6e387ab9
# ╟─b8fbff51-72fe-4432-9f6f-06f4628320d4
# ╠═ef677abe-c96d-41cc-9b4a-192ff29cdb86
# ╠═78933b33-8a3c-492e-8a58-19ad7187b45e
# ╟─9ef4aa78-2ecd-4d53-87ae-f0909bb46915
# ╠═798bc482-7556-4904-8946-0d9898ce2d33
# ╟─32e9f6e5-8df6-4d40-8092-efaf2bce28d3
# ╟─6dfc2699-fe5e-4b55-8437-0866f715170a
# ╠═17d6021b-0d21-4d04-a8ae-64743cc01404
# ╠═f4c3b174-57ac-461b-a816-869460d8896d
# ╟─500973f5-faa7-4a31-b812-d0b20dbb9d82
# ╟─4a575627-32c2-4842-a540-52253f79ff89
# ╟─bea223ef-f961-457c-8962-d0ed396f4717
# ╠═92bd82dc-eaab-4869-9d95-cb4647092352
# ╟─d518828a-1972-4bfa-b96f-0ebbf86b6193
# ╠═c3570557-55bd-4189-a332-aa2e54693665
# ╟─1f684b1b-8e43-4e0e-8a47-1f2474cc5e4f
# ╠═1bf84cfe-9ddc-4644-b7ea-8ac43b01d477
# ╟─d2a1f97e-2e3e-4ae9-8287-2de62065d708
# ╟─a71aaf0b-8ce9-4242-b31f-8f05730ce2e7
# ╟─09ddc00b-2361-44c5-b194-7b4b6d24bef4
# ╠═c7c16e89-3ccf-4d70-9cdb-d5ae71a5ea49
# ╟─1872c30d-aece-44a5-b3b6-36e82c072609
# ╟─5d06734e-ffbd-421e-820e-609e373ccc7b
# ╠═62d55814-c863-4d4d-8582-908721d37326
# ╟─470b7808-5ff6-4cda-82cc-2ac250d51909
# ╟─d392de63-137e-43fb-86ca-4f344769bc71
# ╠═e07830be-3aab-4410-a830-786a30eb162a
# ╠═6d809ff0-20c5-4910-aca0-612496955812
# ╠═816a64ed-1f88-43a8-9df7-8070372e84ea
# ╠═cdd6a0dd-e188-4d37-ab6d-12e7fecd3e90
# ╟─c29d4796-9e52-4234-8325-5bf31099c79d
# ╠═7e5b0073-1bc7-41c7-ae74-c9c1bdf6c6ea
# ╟─67e36f7c-0cb9-4ead-82ad-4c47639802a8
# ╠═cdc25238-ed48-42e5-9db4-2367577b215d
# ╟─f66b6078-bab5-48ef-9241-8d5fb0dc78ea
# ╠═73016af9-5efe-4af2-9e99-ce1a7de6fa96