text:
# PART 5 – ASSESSMENT OF LOCAL METAL LOSS

# Determine Asessment Applicability
#Determine the assessment applicability
# @doc DesignCodeCriteria
# @doc MaterialToughness
# @doc CyclicService
# @doc Part5ComponentType
print("Begin -- Assessment Applicability and Component Type Checks\n")
creep_range = CreepRangeTemperature("Carbon Steel (UTS ≤ 414MPa (60 ksi))"; design_temperature=100.0, units="nmm-mm-mpa")
design = DesignCodeCriteria("ASME B&PV Code, Section VIII, Division 1")
toughness = MaterialToughness("Certain")
cyclic = CyclicService(100, "Meets Part 14")
x = Part5ComponentType("Cylindrical Vessel", vessel_orientation="horizontal", material="Carbon and Low Alloy Steels", D=24.75,Lss=120.0,H=0.0, NPS=3.0, design_temperature=500.0, units="lbs-in-psi")
part5_applicability = Part5AsessmentApplicability(x,design,toughness,cyclic,creep_range)

# For all assessments - determine the inspection data grid
M1 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M2 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.100, 0.220, 0.280, 0.250, 0.240, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M3 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.215, 0.255, 0.215, 0.145, 0.275, 0.170, 0.240, 0.250, 0.250, 0.280, 0.290, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M4 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.170, 0.270, 0.190, 0.190, 0.285, 0.250, 0.225, 0.275, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M5 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M6 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
CTPGrid = hcat(M6,M5,M4,M3,M2,M1) # build in descending order
CTPGrid = rotl90(CTPGrid) # rotate to correct orientation

# Level 1 fit for service
    annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
    # "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
    # "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]
    equipment_group = "piping" # "vessel", "tank"
    flaw_location = "external" # "External","Internal"
    metal_loss_categorization = "LTA" # "LTA" or "Groove-Like Flaw"
    units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"
    tnom = .3 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.
    trd = .3 # uniform thickness away from the local metal loss location established by thickness measurements at the time of the assessment.
    FCAml = 0.0 # Future Corrosion Allowance applied to the region of metal loss.
    FCA = 0.0 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).
    LOSS = 0.0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.
    Do = 3.5 # Outside Diameter
    D = Do - 2*(tnom)
    P = 1480.0 # internal design pressure.
    S = 20000.0 # allowable stress.
    E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.
    MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.
    Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .
    t = trd # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.
    tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).
    spacings = 0.5 # spacings determine by visual inspection to adequately ccategorizse the corrosion -----------+ may add to CTP_Grid function for plotting purposes
    # Flaw dimensions
    s = 6.0 # longitudinal extent or length of the region of local metal loss based on future corroded thickness,
    c = 2.0 # circumferential extent or length of the region of local metal loss (see Figure 5.2 and Figure 5.10), based on future corroded thickness, tc .
    Ec = 1.0 # circumferential weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
    El = 1.0 # longitudinal weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
    RSFa = 0.9 # remaining strength factor - consult API 579 is go lower than 0.9

    # For all assessments determine far enough from structural discontinuity
    # Flaw-To-Major Structural Discontinuity Spacing
    L1msd = [12.0] # distance to the nearest major structural discontinuity.
    L2msd = [12.0] # distance to the nearest major structural discontinuity.
    L3msd = [12.0] # distance to the nearest major structural discontinuity.
    L4msd = [12.0] # distance to the nearest major structural discontinuity.
    L5msd = [12.0] # distance to the nearest major structural discontinuity.
    Lmsd = minimum([L1msd,L2msd,L3msd,L4msd,L5msd])
    if (Lmsd[1] >= (1.8*(sqrt(D*(trd - LOSS - FCA)))))
        print("Satisfied - Flaw is located far enough from structural discontinuity\n")
        lmsd_satisfied = 1
    else
        print("Not satisfied - Flaw is too close to a structural discontinuity - Conduct a level 3 assessment\n")
        lmsd_satisfied = 0
    end

    # Groove Like Flaw dimensions
    gl = .05 # length of the Groove-Like Flaw based on future corroded condition.
    gw = .4 # width of the Groove-Like Flaw based on future corroded condition.
    gr = 0.1 # radius at the base of a Groove-Like Flaw based on future corroded condition.
    β = 40.0 # see (Figure 5.4) :: orientation of the groove-like flaw with respect to the longitudinal axis or a parameter to compute an effective fracture toughness for a groove being evaluated as a crack-like flaw, as applicable.

# Perform level 1 assessment
if (part5_applicability[1] == 1 && lmsd_satisfied == 1) # begin level 1 assessment
    #let part_5_lta_output = Array{Any,2},
    part_5_lta_output = Part5LTALevel1(annex2c_tmin_category; equipment_group=equipment_group, flaw_location=flaw_location, metal_loss_categorization=metal_loss_categorization, units=units, tnom=tnom,
        trd=trd, FCA=FCA, FCAml=FCAml, LOSS=LOSS, Do=Do, D=D, P=P, S=S, E=E, MA=MA, Yb31=Yb31, tsl=tsl, spacings=spacings, s=s, c=c, El=El, Ec=Ec, RSFa=RSFa, gl=gl, gw=gw, gr=gr,β=β)
    #end # let end
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 1 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 1)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
end

#=
# Level 2 Assessment
# STEP 1 – Determine the CTP (see paragraph 5.3.3.2).
# Conducted at the top of the script
# STEP 2 – Determine the wall thickness to be used in the assessment using Equation (5.3) or Equation (5.4), as applicable.
tc = trd - LOSS - FCA # wall thickness away from the damaged area adjusted for LOSS and FCA , as applicable. # eq (5.3)
tc = trd - FCA # wall thickness away from the damaged area adjusted for LOSS and FCA , as applicable. # eq (5.4)
# STEP 3 – Determine the minimum measured thickness, tmm , and the flaw dimensions s and c (see paragraph 5.3.3.2.b).
tmm = CTP_Grid(CTPGrid) # minimum measured thickness determined at the time of the inspection.
# STEP 4 – Determine the remaining thickness ratio, Rt , using Equation (5.5) and the longitudinal flaw length parameter, λ , using Equation (5.6).
Rt = (tmm-FCA) / tc # remaining thickness ratio. # (5.5)
lambda = (1.285*s)/(sqrt(D*tc)) # longitudinal flaw length parameter eq (5.6)
# STEP 5 – Check the limiting flaw size criteria in paragraph 5.4.2.2.e. If all of these requirements are satisfied, then proceed to STEP 6; otherwise, the flaw is not acceptable per the Level 2 Assessment procedure.
x = FlawSizeLimitCriteria("piping","lbs-in-psi")
flaw_size_accept = FlawSizeLevel1Acceptance(x,"piping")

RSFa = .9 # allowable remaining strength factor (see Part 2).
if (Rt < RSFa)
Q = round(1.123*((((1-Rt)/(1-Rt/RSFa))^2-1)^.5),digits=2) # factor used to determine the length for thickness averaging based on an allowable Remaining Strength Factor (see Part 2) and the remaining thickness ratio, Rt (see Table 4.8).
elseif (Rt >= RSFa)
    Q = 50.0
end
L = Q*(sqrt(Dml*tc)) # length for thickness averaging along the shell.

# If visual inspection or NDE methods are utilized to quantify the metal loss, an alternative spacing can be used as long as the metal loss on the component can be adequately characterized.
if (flaw_location == "Internal")
    Ls = minimum([L,(2*trd)]) # recommended spacing of thickness readings
    print("Recommended spacing of thickness readings = ",Ls)
elseif (flaw_location == "External")
    print("Can Determine alterante spacing providing component can be adequately characterized - Using Spacing = ", spacings)
end

# Piping
if (part5_applicability[1] == 1) # begin level 1 assessment
    print("Begin part 5 - Level 1 assessment - applicability criteria has been met")

end # end piping level 1
=#


text:
module Lib # module H3

include("../deps/deps.jl")
include("../gen/libh3_types.jl")
include("../gen/libh3_functions.jl")

using CEnum # @cenum

###
#
# * descriptions are taken from
#   - https://github.com/uber/h3/tree/master/docs/api
#   - https://github.com/uber/h3/blob/master/src/h3lib/include/h3api.h.in
#
###

### Constants

const H3_INIT               = UInt64(35184372088831)  # H3 index with mode 0, res 0, base cell 0, and 7 for all index digits.  0x00001fffffffffff
const H3_MODE_OFFSET        = 59  # The bit offset of the mode in an H3 index.
const H3_MODE_MASK          = UInt64(15) << H3_MODE_OFFSET  # 1's in the 4 mode bits, 0's everywhere else.  0x7800000000000000
const H3_MODE_MASK_NEGATIVE = ~H3_MODE_MASK
const MAX_H3_RES            = 15  # max H3 resolution; H3 version 1 has 16 resolutions, numbered 0 through 15
const H3_PER_DIGIT_OFFSET   = 3   # The number of bits in a single H3 resolution digit.
const H3_DIGIT_MASK         = UInt64(7)  # 1's in the 3 bits of res 15 digit bits, 0's everywhere else.
const H3_HEXAGON_MODE       = 1  # H3 index modes
const H3_UNIEDGE_MODE       = 2  # H3 index modes
const H3_INVALID_INDEX      = UInt64(0)  # Invalid index used to indicate an error from geoToH3 and related functions.
const NUM_BASE_CELLS        = 122  # The number of H3 base cells

### Types

"""
    const H3Index = UInt64

the H3Index fits within a 64-bit unsigned integer
"""
H3Index

"""
    struct GeoCoord
        lat::Cdouble
        lon::Cdouble
    end

latitude/longitude in radians
"""
GeoCoord

"""
    struct CoordIJ
        i::Cint
        j::Cint
    end

IJ hexagon coordinates
"""
CoordIJ

"""
    struct Vec2d
        x::Cdouble
        y::Cdouble
    end

2D floating-point vector
"""
struct Vec2d
    x::Cdouble
    y::Cdouble
end

"""
    struct Vec3d
        x::Cdouble
        y::Cdouble
        z::Cdouble
    end

3D floating point structure
"""
struct Vec3d
    x::Cdouble
    y::Cdouble
    z::Cdouble
end

"""
    struct CoordIJK
        i::Cint
        j::Cint
        k::Cint
    end

IJK hexagon coordinates
Each axis is spaced 120 degrees apart.
"""
struct CoordIJK
    i::Cint
    j::Cint
    k::Cint
end

"""
    struct FaceIJK
        face::Cint
        coord::CoordIJK
    end

Face number and ijk coordinates on that face-centered coordinate system
"""
struct FaceIJK
    face::Cint
    coord::CoordIJK
end

"""
    struct VertexNode
        from::GeoCoord
        to::GeoCoord
        next::Ptr{VertexNode}
    end

A single node in a vertex graph, part of a linked list
"""
struct VertexNode
    from::GeoCoord
    to::GeoCoord
    next::Ptr{VertexNode}
end

"""
    struct VertexGraph
        buckets::Ptr{Ptr{VertexNode}}
        numBuckets::Cint
        size::Cint
        res::Cint
    end

A data structure to store a graph of vertices
"""
struct VertexGraph
    buckets::Ptr{Ptr{VertexNode}}
    numBuckets::Cint
    size::Cint
    res::Cint
end

"""
    Direction

H3 digit representing ijk+ axes direction. Values will be within the lowest 3 bits of an integer.
"""
Direction

@cenum(Direction,
    CENTER_DIGIT  = 0,  # H3 digit in center
    K_AXES_DIGIT  = 1,  # H3 digit in k-axes direction
    J_AXES_DIGIT  = 2,  # H3 digit in j-axes direction
    JK_AXES_DIGIT = 3,  # H3 digit in j == k direction
    I_AXES_DIGIT  = 4,  # H3 digit in i-axes direction
    IK_AXES_DIGIT = 5,  # H3 digit in i == k direction
    IJ_AXES_DIGIT = 6,  # H3 digit in i == j direction
    INVALID_DIGIT = 7,  # H3 digit in the invalid direction
)

### H3Index functions

"""
    setH3Index(refh::Ref{H3Index}, res::Int, baseCell::Int, initDigit::Union{Int, Direction})

Initializes an H3 index.
"""
function setH3Index(refh::Ref{H3Index}, res::Int, baseCell::Int, initDigit::Direction)
    ccall((:setH3Index, libh3), Cvoid, (Ptr{H3Index}, Cint, Cint, Direction), refh, res, baseCell, initDigit)
end

function setH3Index(refh::Ref{H3Index}, res::Int, baseCell::Int, initDigit::Int)
    setH3Index(refh, res, baseCell, Direction(initDigit))
end

"""
    h3GetIndexDigit(h::H3Index, res::Int)::Direction

Gets the resolution res integer digit (0-7) of h3.
`#define H3_GET_INDEX_DIGIT(h3, res)`
"""
function h3GetIndexDigit(h::H3Index, res::Int)::Direction
    Direction((h >> ((MAX_H3_RES - res) * H3_PER_DIGIT_OFFSET)) & H3_DIGIT_MASK)
end

"""
    h3GetMode(h::H3Index)::Int

Gets the integer mode of h3.
`#define H3_GET_MODE(h3)`
"""
function h3GetMode(h::H3Index)::Int
    Int((h & H3_MODE_MASK) >> H3_MODE_OFFSET)
end

"""
    h3SetMode(refh::Ref{H3Index}, v::Int)

Sets the integer mode of h3 to v.
`#define H3_SET_MODE(h3, v)`
"""
function h3SetMode(refh::Ref{H3Index}, v::Int)
    refh[] = (refh[] & H3_MODE_MASK_NEGATIVE) | (UInt64(v) << H3_MODE_OFFSET)
end

### Region functions

"""
    polyfill(refpolygon::Ref{GeoPolygon}, res::Int)::Vector{H3Index}

polyfill takes a given GeoJSON-like data structure and preallocated, zeroed memory, and fills it with the hexagons that are contained by the GeoJSON-like data structure.
"""
function polyfill(refpolygon::Ref{GeoPolygon}, res::Int)::Vector{H3Index}
    numHexagons = Lib.maxPolyfillSize(refpolygon, res)
    out = Vector{H3Index}(undef, numHexagons)
    Lib.polyfill(refpolygon, res, out)
    out
end

"""
    maxPolyfillSize(refpolygon::Ref{GeoPolygon}, res::Int)::Int

maxPolyfillSize returns the number of hexagons to allocate space for when performing a polyfill on the given GeoJSON-like data structure.
"""
maxPolyfillSize

"""
    h3SetToLinkedGeo(h3Set::Vector{H3Index})::Ref{LinkedGeoPolygon}

Create a LinkedGeoPolygon describing the outline(s) of a set of hexagons. Polygon outlines will follow GeoJSON MultiPolygon order: Each polygon will have one outer loop, which is first in the list, followed by any holes.
"""
function h3SetToLinkedGeo(h3Set::Vector{H3Index})::Ref{LinkedGeoPolygon}
    refpolygon = Ref{LinkedGeoPolygon}(LinkedGeoPolygon(C_NULL,C_NULL,C_NULL))
    Lib.h3SetToLinkedGeo(C_NULL, 0, refpolygon)
    refpolygon
end

"""
    destroyLinkedPolygon(refpolygon::Ref{LinkedGeoPolygon})

Free all allocated memory for a linked geo structure. The caller is responsible for freeing memory allocated to the input polygon struct.
"""
destroyLinkedPolygon

### Vertex Graph

"""
    initVertexGraph(refgraph::Ref{VertexGraph}, numBuckets::Int, res::Int)

Initialize a new VertexGraph
"""
function initVertexGraph(refgraph::Ref{VertexGraph}, numBuckets::Int, res::Int)
    ccall((:initVertexGraph, libh3), Cvoid, (Ptr{VertexGraph}, Cint, Cint), refgraph, numBuckets, res)
end

"""
    destroyVertexGraph(refgraph::Ref{VertexGraph})

Destroy a VertexGraph's sub-objects, freeing their memory. The caller is responsible for freeing memory allocated to the VertexGraph struct itself.
"""
function destroyVertexGraph(refgraph::Ref{VertexGraph})
    ccall((:destroyVertexGraph, libh3), Cvoid, (Ptr{VertexGraph},), refgraph)
end

end # module H3.Lib

ext: jl


text:
include(joinpath(pwd(), "examples/implicit_dynamics/utils.jl"))

"""
    Hopper2D
    	model inspired by "Dynamically Stable Legged Locomotion"
		s = (x, z, t, r)
			x - lateral position
			z - vertical position
			t - body orientation
			r - leg length
"""
mutable struct Hopper2D{T}
    dim::Dimensions

    mb::T # mass of body
    ml::T # mass of leg
    Jb::T # inertia of body
    Jl::T # inertia of leg

    μ_world::T  # coefficient of friction
    μ_joint::T  # gravity
	g::T

	joint_friction::SVector
end

lagrangian(model::Hopper2D, q, q̇) = 0.0

# Kinematics
function kinematics(::Hopper2D, q)
	[q[1] + q[4] * sin(q[3]),
	 q[2] - q[4] * cos(q[3])]
end

# Methods
function M_func(model::Hopper2D, q)
	Diagonal([model.mb + model.ml,
    					   model.mb + model.ml,
    					   model.Jb + model.Jl,
    					   model.ml])
 end

function C_func(model::Hopper2D, q, q̇)
	[0.0,
	  (model.mb + model.ml) * model.g,
	  0.0,
	  0.0]
end

function ϕ_func(model::Hopper2D, q)
    [q[2] - q[4] * cos(q[3])]
end

function J_func(::Hopper2D, q)
    [1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
		          0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]))]
end

function B_func(::Hopper2D, q)
	[0.0 0.0 1.0 0.0;
                -sin(q[3]) cos(q[3]) 0.0 1.0]
end

# Working Parameters
gravity = 9.81 # gravity
μ_world = 0.8 # coefficient of friction
μ_joint = 0.0

# TODO: change to Raibert parameters
mb = 3.0 # body mass
ml = 0.3  # leg mass
Jb = 0.75 # body inertia
Jl = 0.075 # leg inertia

# Dimensions
nq = 4
nu = 2
nw = 2
nc = 1
nb = 2

model = Hopper2D(Dimensions(nq, nu, nw, nc),
			   mb, ml, Jb, Jl,
			   μ_world, μ_joint, gravity,
			   SVector{4}(zeros(4)))

function lagrangian_derivatives(model::Hopper2D, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

lagrangian_derivatives(model, rand(nq), rand(nq))

function dynamics(model::Hopper2D, h, q0, q1, u1, λ1, q2)
    # evalutate at midpoint
    qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

    D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
    D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

    return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2
    	+ transpose(B_func(model, qm2)) * u1
        + transpose(J_func(model, q2)) * λ1)
        # -h[1] * model.joint_friction .* vm2)
end

function residual(model, z, θ, κ)
    nq = model.dim.q
    nu = model.dim.u
    nc = model.dim.c
    nb = nc * 2

    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    h = θ[2nq + nu .+ (1:1)]

    μ = model.μ_world

    q2 = z[1:nq]
    γ1 = z[nq .+ (1:nc)]
    b1 = z[nq + nc .+ (1:nb)]
    ψ1 = z[nq + nc + nb .+ (1:nc)]
    η1 = z[nq + nc + nb + nc .+ (1:nb)]
    s1 = z[nq + nc + nb + nc + nb .+ (1:nc)]
    s2 = z[nq + nc + nb + nc + nb + nc .+ (1:nc)]

	ϕ = ϕ_func(model, q2)

	λ1 = [b1[1] - b1[2]; γ1]
    vT = (J_func(model, q2) * (q2 - q1) / h[1])[1]
	vT_stack = [vT; -vT]
	ψ_stack = ψ1 .* ones(nb)

	[
     dynamics(model, h, q0, q1, u1, λ1, q2);
	 s1 - ϕ;
	 vT_stack + ψ_stack - η1;
	 s2 .- (μ[1] * γ1 .- sum(b1));
	 γ1 .* s1 .- κ;
	 b1 .* η1 .- κ;
	 ψ1 .* s2 .- κ
    ]
end

nz = nq + nc + nb + nc + nb + nc + nc
nθ = nq + nq + nu + 1

idx_ineq = collect(nq .+ (1:(nc + nb + nc + nb + nc + nc)))
z_subset_init = 0.1 * ones(nc + nb + nc + nb + nc + nc)

# Declare variables
@variables z[1:nz]
@variables θ[1:nθ]
@variables κ[1:1]

# Residual
r = residual(model, z, θ, κ)
r = Symbolics.simplify.(r)
rz = Symbolics.jacobian(r, z, simplify = true)
rθ = Symbolics.jacobian(r, θ, simplify = true) # TODO: sparse version

# Build function
r_func = eval(build_function(r, z, θ, κ)[2])
rz_func = eval(build_function(rz, z, θ)[2])
rθ_func = eval(build_function(rθ, z, θ)[2])

rz_array = similar(rz, Float64)
rθ_array = similar(rθ, Float64)

@save joinpath(@__DIR__, "dynamics/residual.jl") r_func rz_func rθ_func rz_array rθ_array
@load joinpath(@__DIR__, "dynamics/residual.jl") r_func rz_func rθ_func rz_array rθ_array

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays, Base64
import LazyArrays: MulAdd, MemoryLayout, DenseColumnMajor, DiagonalLayout, SymTridiagonalLayout, Add, AddArray,
                    MulStyle, MulAddStyle, Applied, ApplyStyle, Lmul, ApplyArrayBroadcastStyle, DefaultArrayApplyStyle,
                    Rmul, ApplyLayout, arguments, colsupport, rowsupport, lazymaterialize
import Base.Broadcast: materialize, materialize!, broadcasted
import MatrixFactorizations: QRCompactWYQLayout, AdjQRCompactWYQLayout

@testset "Matrix * Vector" begin
    @testset "eltype" begin
        @test @inferred(eltype(Applied(*, zeros(Int,2,2), zeros(Float64,2)))) == Float64
        @test @inferred(eltype(Applied(*, zeros(ComplexF16,2,2),zeros(Int,2,2),zeros(Float64,2)))) == ComplexF64

        v = Applied(*, zeros(Int,2,2), zeros(Float64,2))
        A = Applied(*, zeros(Int,2,2), zeros(Float64,2,2))
        @test @inferred(axes(v)) == (@inferred(axes(v,1)),) == (Base.OneTo(2),)
        @test @inferred(size(v)) == (@inferred(size(v,1)),) == (2,)
        @test @inferred(axes(A)) == (@inferred(axes(A,1)),@inferred(axes(A,2))) == (Base.OneTo(2),Base.OneTo(2))
        @test @inferred(size(A)) == (@inferred(size(A,1)),size(A,2)) == (2,2)

        Ã = materialize(A)
        @test Ã isa Matrix{Float64}
        @test Ã == zeros(2,2)

        c = similar(v)
        @test c isa Vector
        fill!(c,NaN)
        @test copyto!(c,v) == zeros(2)
        fill!(c,NaN)
        c .= v
        @test c == zeros(2)

        A = randn(6,5); B = randn(5,5); C = randn(5,6)
        M = Applied(*, A,B,C)
        @test @inferred(eltype(M)) == Float64
        @test materialize(M) == A*B*C
    end

    @testset "gemv Float64" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
            b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))
            c = similar(b);

            c .= applied(*,A,b)
            @test applied(*,A,b) isa Applied{MulStyle}
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))
            copyto!(c, applied(*,A,b))
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))
            c .= @~ A*b
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))

            b̃ = copy(b)
            copyto!(b̃, applied(*,A,b̃))
            @test all(c .=== b̃)

            M = applied(*, 2.0, A,b)
            @test M isa Applied{MulAddStyle}
            c .= M
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))
            copyto!(c, M)
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))
            c .= @~ 2.0*A*b
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))

            c = copy(b)
            M = applied(+, applied(*,A,b), c)
            @test M isa Applied{MulAddStyle}
            copyto!(c, M)
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))
            c = copy(b)
            c .= applied(+, applied(*,A,b), c)
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))
            c = copy(b)
            c .= @~ A*b + c
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))

            c = copy(b)
            @test (@~ A*b + 2.0*c) isa Applied{MulAddStyle}
            c .= @~ A*b + 2.0*c
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 2.0, copy(b)))

            c = copy(b)
            @test (@~ 2.0*A*b + c) isa Applied{MulAddStyle}
            c .= @~ 2.0*A*b + c
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 1.0, copy(b)))

            c = copy(b)
            @test ( @~ 3.0*A*b + 2.0*c) isa Applied{MulAddStyle}
            c .= @~ 3.0*A*b + 2.0*c
            @test all(c .=== BLAS.gemv!('N', 3.0, A, b, 2.0, copy(b)))

            d = similar(c)
            c = copy(b)
            d .= @~ 3.0*A*b + 2.0*c
            @test all(d .=== BLAS.gemv!('N', 3.0, A, b, 2.0, copy(b)))
        end
    end

    @testset "gemv mixed array types" begin
        (A, b, c) = (randn(5,5), randn(5), 1.0:5.0)
        d = similar(b)
        d .= @~ A*b + c
        @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 1.0, Vector{Float64}(c)))

        d .= @~ A*b + 2.0*c
        @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 2.0, Vector{Float64}(c)))

        d .= @~ 2.0*A*b + c
        @test all(d .=== BLAS.gemv!('N', 2.0, A, b, 1.0, Vector{Float64}(c)))

        d .= @~ 3.0*A*b + 2.0 * c
        @test all(d .=== BLAS.gemv!('N', 3.0, A, b, 2.0, Vector{Float64}(c)))

        @test (similar(b) .= @~ 1 * A * b) ≈
            (similar(b) .= @~ 1 * A *b + 0 * b) ≈
            A*b
    end

    @testset "gemv Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                b in (randn(T,5), view(randn(T,5),:))
            c = similar(b);

            c .= @~ A*b
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, zero(T), similar(c)))

            b .= @~ A*b
            @test all(c .=== b)

            c .= applied(*, 2one(T), A, b)
            @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, zero(T), similar(c)))

            c = copy(b)
            c .= @~ A*b + c
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+, @~(A*b), applied(*, 2one(T), c))
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, 2one(T), copy(b)))

            c = copy(b)
            c .= applied(+, applied(*, 2one(T), A, b), c)
            @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+, applied(*, 3one(T), A, b), applied(*, 2one(T), c))
            @test all(c .=== BLAS.gemv!('N', 3one(T), A, b, 2one(T), copy(b)))

            d = similar(c)
            c = copy(b)
            d .= applied(+, applied(*, 3one(T), A, b), applied(*, 2one(T), c))
            @test all(d .=== BLAS.gemv!('N', 3one(T), A, b, 2one(T), copy(b)))
        end
    end
end

@testset "Matrix * Matrix" begin
    @testset "gemm" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
            B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                      view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))
            C = similar(B);

            C .= @~ A*B
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, similar(C)))

            B .= @~ A*B
            @test all(C .=== B)

            C .= @~ 2.0 * A*B
            @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 0.0, similar(C)))

            C = copy(B)
            C .= @~ A*B + C
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, copy(B)))

            C = copy(B)
            C .= @~ A*B + 2.0 * C
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 2.0, copy(B)))

            C = copy(B)
            C .= @~ 2.0 * A*B + C
            @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, copy(B)))

            C = copy(B)
            C .= @~ 3.0 * A*B + 2.0 * C
            @test all(C .=== BLAS.gemm!('N', 'N', 3.0, A, B, 2.0, copy(B)))

            d = similar(C)
            C = copy(B)
            d .= @~ 3.0 * A*B + 2.0 * C
            @test all(d .=== BLAS.gemm!('N', 'N', 3.0, A, B, 2.0, copy(B)))
        end
    end

    @testset "gemm Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                B in (randn(T,5,5), view(randn(T,5,5),:,:))
            C = similar(B);

            C .= @~ A*B
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, zero(T), similar(C)))

            B .= @~ A*B
            @test all(C .=== B)

            C .= applied(*, 2one(T), A, B)
            @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, zero(T), similar(C)))

            C = copy(B)
            C .= @~ A*B + C
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, one(T), copy(B)))

            C = copy(B)
            C .= applied(+, @~(A*B), applied(*, 2one(T), C))
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, 2one(T), copy(B)))

            C = copy(B)
            C .= applied(+, applied(*, 2one(T), A,B), C)
            @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, one(T), copy(B)))

            C = copy(B)
            C .= applied(+, applied(*, 3one(T),A,B), applied(*, 2one(T), C))
            @test all(C .=== BLAS.gemm!('N', 'N', 3one(T), A, B, 2one(T), copy(B)))

            d = similar(C)
            C = copy(B)
            d .= applied(+, applied(*,3one(T),A,B), applied(*,2one(T), C))
            @test all(d .=== BLAS.gemm!('N', 'N', 3one(T), A, B, 2one(T), copy(B)))
        end
    end

    @testset "gemm mixed array types" begin
        (A, B, C) = (randn(5,5), randn(5,5), reshape(1.0:25.0,5,5))
        D = similar(B)
        D .= @~ A*B + C
        @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, Matrix{Float64}(C)))

        D .= @~ A*B + 2.0 * C
        @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 2.0, Matrix{Float64}(C)))

        D .= @~ 2.0 * A*B + C
        @test all(D .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, Matrix{Float64}(C)))

        D .= @~ 3.0 * A*B + 2.0 * C
        @test all(D .=== BLAS.gemm!('N', 'N', 3.0, A, B, 2.0, Matrix{Float64}(C)))

        @test (similar(B) .= 1 * A*B) ≈
            (similar(B) .= 1 * A*B + 0 * B) ≈
            A*B
    end
end

@testset "adjtrans" begin
    @testset "gemv adjtrans" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
                    b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)),
                    Ac in (transpose(A), A')
            c = similar(b);

            c .= @~ Ac*b
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 0.0, similar(c)))

            b .= @~ Ac*b
            @test all(c .=== b)

            c .= @~ 2.0 * Ac*b
            @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 0.0, similar(c)))

            c = copy(b)
            c .= @~ Ac*b + c
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 1.0, copy(b)))

            c = copy(b)
            c .= @~ Ac*b + 2.0 * c
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 2.0, copy(b)))

            c = copy(b)
            c .= @~ 2.0 * Ac*b + c
            @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 1.0, copy(b)))

            c = copy(b)
            c .= @~ 3.0 * Ac*b + 2.0 * c
            @test all(c .=== BLAS.gemv!('T', 3.0, A, b, 2.0, copy(b)))

            d = similar(c)
            c = copy(b)
            d .= @~ 3.0 * Ac*b + 2.0 * c
            @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, copy(b)))
        end
    end

    @testset "gemv adjtrans mixed types" begin
        (A, b, c) = (randn(5,5), randn(5), 1.0:5.0)
        for Ac in (transpose(A), A')
            d = similar(b)
            d .= @~ Ac*b + c
            @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 1.0, Vector{Float64}(c)))

            d .= @~ Ac*b + 2.0 * c
            @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 2.0, Vector{Float64}(c)))

            d .= @~ 2.0 * Ac*b + c
            @test all(d .=== BLAS.gemv!('T', 2.0, A, b, 1.0, Vector{Float64}(c)))

            d .= @~ 3.0 * Ac*b + 2.0 * c
            @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, Vector{Float64}(c)))
        end
    end

    @testset "gemv adjtrans Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                b in (randn(T,5), view(randn(T,5),:)),
                (Ac,trans) in ((transpose(A),'T'), (A','C'))
            c = similar(b);

            c .= @~ Ac*b
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, zero(T), similar(c)))

            b .= @~ Ac*b
            @test all(c .=== b)

            c .= applied(*,2one(T),Ac,b)
            @test all(c .=== BLAS.gemv!(trans, 2one(T), A, b, zero(T), similar(c)))

            c = copy(b)
            c .= @~ Ac*b + c
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+, applied(*,Ac,b), applied(*,2one(T), c))
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, 2one(T), copy(b)))

            c = copy(b)
            c .= applied(+,applied(*,2one(T),Ac,b), c)
            @test all(c .=== BLAS.gemv!(trans, 2one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+,applied(*,3one(T),Ac,b), applied(*,2one(T), c))
            @test all(c .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b)))

            d = similar(c)
            c = copy(b)
            d .= applied(+,applied(*,3one(T),Ac,b), applied(*,2one(T), c))
            @test all(d .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b)))
        end
    end

    @testset "gemm adjtrans" begin
        for A in (randn(5,5), view(randn(5,5),:,:)),
            B in (randn(5,5), view(randn(5,5),1:5,:))
            for Ac in (transpose(A), A')
                C = similar(B)
                C .= @~ Ac*B
                @test all(C .=== BLAS.gemm!('T', 'N', 1.0, A, B, 0.0, similar(C)))

                B .= @~ Ac*B
                @test all(C .=== B)

                C .= @~ 2.0 * Ac*B
                @test all(C .=== BLAS.gemm!('T', 'N', 2.0, A, B, 0.0, similar(C)))

                C = copy(B)
                C .= @~ Ac*B + C
                @test all(C .=== BLAS.gemm!('T', 'N', 1.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ Ac*B + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'N', 1.0, A, B, 2.0, copy(B)))

                C = copy(B)
                C .= @~ 2.0 * Ac*B + C
                @test all(C .=== BLAS.gemm!('T', 'N', 2.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ 3.0 * Ac*B + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B)))

                d = similar(C)
                C = copy(B)
                d .= @~ 3.0 * Ac*B + 2.0 * C
                @test all(d .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B)))
            end
            for Bc in (transpose(B), B')
                C = similar(B)
                C .= @~ A*Bc
                @test all(C .=== BLAS.gemm!('N', 'T', 1.0, A, B, 0.0, similar(C)))

                B .= @~ A*Bc
                @test all(C .=== B)

                C .= @~ 2.0 * A*Bc
                @test all(C .=== BLAS.gemm!('N', 'T', 2.0, A, B, 0.0, similar(C)))

                C = copy(B)
                C .= @~ A*Bc + C
                @test all(C .=== BLAS.gemm!('N', 'T', 1.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ A*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('N', 'T', 1.0, A, B, 2.0, copy(B)))

                C = copy(B)
                C .= @~ 2.0 * A*Bc + C
                @test all(C .=== BLAS.gemm!('N', 'T', 2.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ 3.0 * A*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B)))

                d = similar(C)
                C = copy(B)
                d .= @~ 3.0 * A*Bc + 2.0 * C
                @test all(d .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B)))
            end
            for Ac in (transpose(A), A'), Bc in (transpose(B), B')
                C = similar(B)
                C .= @~ Ac*Bc
                @test all(C .=== BLAS.gemm!('T', 'T', 1.0, A, B, 0.0, similar(C)))

                B .= @~ Ac*Bc
                @test all(C .=== B)

                C .= @~ 2.0 * Ac*Bc
                @test all(C .=== BLAS.gemm!('T', 'T', 2.0, A, B, 0.0, similar(C)))

                C = copy(B)
                C .= @~ Ac*Bc + C
                @test all(C .=== BLAS.gemm!('T', 'T', 1.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ Ac*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'T', 1.0, A, B, 2.0, copy(B)))

                C = copy(B)
                C .= @~ 2.0 * Ac*Bc + C
                @test all(C .=== BLAS.gemm!('T', 'T', 2.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ 3.0 * Ac*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B)))

                d = similar(C)
                C = copy(B)
                d .= @~ 3.0 * Ac*Bc + 2.0 * C
                @test all(d .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B)))
            end
        end
    end

    @testset "symv adjtrans" begin
        A = randn(100,100)
        x = randn(100)

        @test all( (similar(x) .= applied(*,Symmetric(A),x)) .===
                    (similar(x) .= applied(*,Hermitian(A),x)) .===
                    (similar(x) .= applied(*,Symmetric(A)',x)) .===
                    (similar(x) .= applied(*,Symmetric(view(A,:,:)',:L),x)) .===
                    (similar(x) .= applied(+,applied(*,1.0,Symmetric(A),x), applied(*,0.0,similar(x)))) .===
                    (similar(x) .= applied(*,view(Symmetric(A),:,:),x)) .===
                    BLAS.symv!('U', 1.0, A, x, 0.0, similar(x)) )

        y = copy(x)
        y .= applied(*,Symmetric(A), y)
        @test all( (similar(x) .= applied(*, Symmetric(A),x)) .=== y)

        @test all( (similar(x) .= applied(*, Symmetric(A,:L),x)) .===
                    (similar(x) .= applied(*, Hermitian(A,:L),x)) .===
                    (similar(x) .= applied(*, Symmetric(A,:L)',x)) .===
                    (similar(x) .= applied(*, Symmetric(view(A,:,:)',:U),x)) .===
                    (fill!(similar(x),NaN) .= applied(+,applied(*,1.0,Symmetric(A,:L),x), applied(*,0.0,similar(x)))) .===
                    (fill!(similar(x),NaN) .= applied(*, view(Hermitian(A,:L),:,:),x)) .===
                    BLAS.symv!('L', 1.0, A, x, 0.0, similar(x)) )
        T = ComplexF64
        A = randn(T,100,100)
        x = randn(T,100)
        @test all( (similar(x) .= applied(*, Symmetric(A),x)) .===
                    (similar(x) .= applied(*, transpose(Symmetric(A)),x)) .===
                    (similar(x) .= applied(*, Symmetric(transpose(view(A,:,:)),:L),x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Symmetric(A),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= applied(*, view(Symmetric(A),:,:),x)) .===
                    BLAS.symv!('U', one(T), A, x, zero(T), similar(x)) )

        @test all( (similar(x) .= applied(*, Symmetric(A,:L),x)) .===
                    (similar(x) .= applied(*, transpose(Symmetric(A,:L)),x)) .===
                    (similar(x) .= applied(*, Symmetric(transpose(view(A,:,:)),:U),x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Symmetric(A,:L),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= applied(*, view(Symmetric(A,:L),:,:),x)) .===
                    BLAS.symv!('L', one(T), A, x, zero(T), similar(x)) )

        @test all( (similar(x) .= applied(*, Hermitian(A),x)) .===
                    (similar(x) .= applied(*, Hermitian(A)',x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Hermitian(A),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= applied(*, view(Hermitian(A),:,:),x)) .===
                    BLAS.hemv!('U', one(T), A, x, zero(T), similar(x)) )

        @test all( (similar(x) .= applied(*, Hermitian(A,:L),x)) .===
                    (similar(x) .= applied(*, Hermitian(A,:L)',x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Hermitian(A,:L),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= applied(*, view(Hermitian(A,:L),:,:),x)) .===
                    BLAS.hemv!('L', one(T), A, x, zero(T), similar(x)) )

        y = copy(x)
        y .= applied(*, Hermitian(A), y)
        @test all( (similar(x) .= applied(*, Hermitian(A),x)) .=== y)
    end
end

@testset "Applied(*, ...)" begin
    @testset "Mixed types" begin
        A = randn(5,6)
        b = rand(Int,6)
        c = Array{Float64}(undef, 5)
        c .= applied(*,A,b)
        d = similar(c)
        mul!(d, A, b)
        @test c ≈ d

        @test_throws DimensionMismatch (similar(c,3) .= applied(*,A,b))
        @test_throws DimensionMismatch (c .= applied(*,A,similar(b,2)))

        copyto!(d, MulAdd(1, A, b, 0.0, d))
        @test d == copyto!(similar(d), MulAdd(1, A, b, 0.0, d)) ≈ A*b
        @test copyto!(similar(d), MulAdd(1, A, b, 1.0, d)) ≈ A*b + d

        @test all((similar(d) .= MulAdd(1, A, b, 1.0, d)) .=== copyto!(similar(d), MulAdd(1, A, b, 1.0, d)))

        B = rand(Int,6,4)
        C = Array{Float64}(undef, 5, 4)
        C .= applied(*,A,B)

        @test_throws DimensionMismatch materialize!(MulAdd(1,A,B,0,similar(C,4,4)))
        @test_throws DimensionMismatch (similar(C,4,4) .= applied(*,A,B))
        @test_throws DimensionMismatch (C .= applied(*,A,similar(B,2,2)))

        D = similar(C)
        mul!(D, A, B)
        @test all(C .=== D)

        A = randn(Float64,20,22)
        B = randn(ComplexF64,22,24)
        C = similar(B,20,24)
        @test all((C .= applied(*,A,B)  ) .=== copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)))
        @test copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)) ≈ A*B
    end

    @testset "no allocation" begin
        function blasnoalloc(c, α, A, x, β, y)
            c .= @~ A*x
            c .= @~ α*A*x
            c .= @~ A*x + y
            c .= @~ α * A*x + y
            c .= @~ A*x + β * y
            c .= @~ α * A*x + β * y
        end

        A = randn(5,5); x = randn(5); y = randn(5); c = similar(y);

        @inferred(MulAdd(@~ A*x + y))
        @test blasnoalloc(c, 2.0, A, x, 3.0, y) === c
        @test @allocated(blasnoalloc(c, 2.0, A, x, 3.0, y)) == 0
        Ac = A'
        blasnoalloc(c, 2.0, Ac, x, 3.0, y)
        @test @allocated(blasnoalloc(c, 2.0, Ac, x, 3.0, y)) == 0
        Aa = ApplyArray(+, A, Ac)
        blasnoalloc(c, 2.0, Aa, x, 3.0, y)
        @test_broken @allocated(blasnoalloc(c, 2.0, Aa, x, 3.0, y)) == 0
    end

    @testset "multi-argument mul" begin
        A = randn(5,5)
        B = apply(*,A,A,A)
        @test B isa Matrix{Float64}
        @test all(B .=== (A*A)*A)
    end

    @testset "#14" begin
        A = ones(1,1) * 1e200
        B = ones(1,1) * 1e150
        C = ones(1,1) * 1e-300

        @test apply(*, A, applied(*,B,C)) == A*(B*C)
        @test apply(*, A , applied(*,B,C), C) == A * (B*C) * C
    end

    @testset "#15" begin
        N = 10
        A = randn(N,N); B = randn(N,N); C = randn(N,N); R1 = similar(A); R2 = similar(A)
        M = applied(*, A, applied(*, B, C))
        @test axes(M) == (Base.OneTo(N),Base.OneTo(N))
        @test ndims(M) == ndims(typeof(M)) == 2
        @test eltype(M) == Float64
        @test all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M))
    end

    @testset "broadcasting" begin
        A = randn(5,5); B = randn(5,5); C = randn(5,5)
        C .= NaN
        C .= @~ 1.0 * A*B + 0.0 * C
        @test C == A*B
    end

    @testset "BigFloat" begin
        A = BigFloat.(randn(5,5))
        x = BigFloat.(randn(5))
        @test A*x == apply(*,A,x) == copyto!(similar(x), applied(*,A,x))
        @test_throws UndefRefError materialize!(MulAdd(1.0,A,x,0.0,similar(x)))
    end

    @testset "Scalar * Vector" begin
        A, x =  [1 2; 3 4] , [[1,2],[3,4]]
        @test apply(*,A,x) == A*x
    end

    @testset "Complex broadcast" begin
        A = randn(5,5) .+ im*randn(5,5)
        x = randn(5) .+ im*randn(5)
        y = randn(5) .+ im*randn(5)
        z = similar(x)
        @test all((z .= @~ (2.0+0.0im)*A*x + (3.0+0.0im)*y) .=== BLAS.gemv!('N',2.0+0.0im,A,x,3.0+0.0im,y))
    end
end

@testset "Lmul/Rmul" begin
    @testset "tri Lmul" begin
        @testset "Float * Float vector" begin
            A = randn(Float64, 100, 100)
            x = randn(Float64, 100)

            L = Lmul(UpperTriangular(A),x)
            @test size(L) == (size(L,1),) == (100,)
            @test axes(L) == (axes(L,1),) == (Base.OneTo(100),)
            @test eltype(L) == Float64
            @test length(L) == 100

            @test similar(L) isa Vector{Float64}
            @test similar(L,Int) isa Vector{Int}

            @test applied(*, UpperTriangular(A), x) isa Applied{MulStyle}
            @test similar(applied(*, UpperTriangular(A), x), Float64) isa Vector{Float64}

            @test ApplyStyle(*, typeof(UpperTriangular(A)), typeof(x)) isa MulStyle

            @test all((y = copy(x); y .= applied(*, UpperTriangular(A),y) ) .===
                        (similar(x) .= applied(*, UpperTriangular(A),x)) .===
                        copy(Lmul(UpperTriangular(A),x)) .===
                        materialize!(Lmul(UpperTriangular(A),copy(x))) .===
                        copyto!(similar(x),Lmul(UpperTriangular(A),x)) .===
                        UpperTriangular(A)*x .===
                        BLAS.trmv!('U', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitUpperTriangular(A),y) ) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(A),x)) .===
                        UnitUpperTriangular(A)*x .===
                        BLAS.trmv!('U', 'N', 'U', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, LowerTriangular(A),y) ) .===
                        (similar(x) .= applied(*, LowerTriangular(A),x)) .===
                        LowerTriangular(A)*x .===
                        BLAS.trmv!('L', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitLowerTriangular(A),y) ) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(A),x)) .===
                        UnitLowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'U', A, copy(x)))

            @test all((y = copy(x); y .= applied(*, UpperTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, UpperTriangular(A)',x)) .===
                        (similar(x) .= applied(*, LowerTriangular(A'),x)) .===
                        UpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitUpperTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(A)',x)) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(A'),x)) .===
                        UnitUpperTriangular(A)'*x .===
                        BLAS.trmv!('U', 'T', 'U', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, LowerTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, LowerTriangular(A)',x)) .===
                        (similar(x) .= applied(*, UpperTriangular(A'),x)) .===
                        LowerTriangular(A)'*x .===
                        BLAS.trmv!('L', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitLowerTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(A)',x)) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(A'),x)) .===
                        UnitLowerTriangular(A)'*x .===
                        BLAS.trmv!('L', 'T', 'U', A, copy(x)))
        end

        @testset "Float * Complex vector"  begin
            T = ComplexF64
            A = randn(T, 100, 100)
            x = randn(T, 100)

            @test all((y = copy(x); y .= applied(*, UpperTriangular(A),y) ) .===
                        (similar(x) .= applied(*, UpperTriangular(A),x)) .===
                        UpperTriangular(A)x .===
                        BLAS.trmv!('U', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitUpperTriangular(A),y) ) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(A),x)) .===
                        UnitUpperTriangular(A)x .===
                        BLAS.trmv!('U', 'N', 'U', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, LowerTriangular(A),y) ) .===
                        (similar(x) .= applied(*, LowerTriangular(A),x)) .===
                        LowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitLowerTriangular(A),y) ) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(A),x)) .===
                        UnitLowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'U', A, copy(x)))
            LowerTriangular(A')  == UpperTriangular(A)'

            @test all((y = copy(x); y .= applied(*, transpose(UpperTriangular(A)),y) ) .===
                        (similar(x) .= applied(*, transpose(UpperTriangular(A)),x)) .===
                        (similar(x) .= applied(*, LowerTriangular(transpose(A)),x)) .===
                        transpose(UpperTriangular(A))x .===
                        BLAS.trmv!('U', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, transpose(UnitUpperTriangular(A)),y) ) .===
                        (similar(x) .= applied(*, transpose(UnitUpperTriangular(A)),x)) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(transpose(A)),x)) .===
                        transpose(UnitUpperTriangular(A))x .===
                        BLAS.trmv!('U', 'T', 'U', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, transpose(LowerTriangular(A)),y) ) .===
                        (similar(x) .= applied(*, transpose(LowerTriangular(A)),x)) .===
                        (similar(x) .= applied(*, UpperTriangular(transpose(A)),x)) .===
                        transpose(LowerTriangular(A))x .===
                        BLAS.trmv!('L', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, transpose(UnitLowerTriangular(A)),y) ) .===
                        (similar(x) .= applied(*, transpose(UnitLowerTriangular(A)),x)) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(transpose(A)),x)) .===
                        transpose(UnitLowerTriangular(A))x .===
                        BLAS.trmv!('L', 'T', 'U', A, copy(x)))

            @test all((y = copy(x); y .= applied(*, UpperTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, UpperTriangular(A)',x)) .===
                        (similar(x) .= applied(*, LowerTriangular(A'),x)) .===
                        UpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'C', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitUpperTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(A)',x)) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(A'),x)) .===
                        UnitUpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'C', 'U', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, LowerTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, LowerTriangular(A)',x)) .===
                        (similar(x) .= applied(*, UpperTriangular(A'),x)) .===
                        LowerTriangular(A)'x .===
                        BLAS.trmv!('L', 'C', 'N', A, copy(x)))
            @test all((y = copy(x); y .= applied(*, UnitLowerTriangular(A)',y) ) .===
                        (similar(x) .= applied(*, UnitLowerTriangular(A)',x)) .===
                        (similar(x) .= applied(*, UnitUpperTriangular(A'),x)) .===
                        UnitLowerTriangular(A)'x .===
                        BLAS.trmv!('L', 'C', 'U', A, copy(x)))
        end

        @testset "Float * Float Matrix" begin
            A = randn(Float64, 100, 100)
            x = randn(Float64, 100, 100)

            @test UpperTriangular(A)*x ≈ (similar(x) .= applied(*, UpperTriangular(A), x))
        end

        @testset "adjtrans" begin
            for T in (Float64, ComplexF64)
                A = randn(T,100,100)
                b = randn(T,100)

                @test all(apply(*, UpperTriangular(A)', b) .=== UpperTriangular(A)'b)
                @test all(apply(*, UnitUpperTriangular(A)', b) .=== UnitUpperTriangular(A)'b)
                @test all(apply(*, LowerTriangular(A)', b) .=== LowerTriangular(A)'b)
                @test all(apply(*, UnitLowerTriangular(A)', b) .=== UnitLowerTriangular(A)'b)

                @test all(apply(*, transpose(UpperTriangular(A)), b) .=== transpose(UpperTriangular(A))b)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), b) .=== transpose(UnitUpperTriangular(A))b)
                @test all(apply(*, transpose(LowerTriangular(A)), b) .=== transpose(LowerTriangular(A))b)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), b) .=== transpose(UnitLowerTriangular(A))b)

                B = randn(T,100,100)

                @test all(apply(*, UpperTriangular(A)', B) .=== UpperTriangular(A)'B)
                @test all(apply(*, UnitUpperTriangular(A)', B) .=== UnitUpperTriangular(A)'B)
                @test all(apply(*, LowerTriangular(A)', B) .=== LowerTriangular(A)'B)
                @test all(apply(*, UnitLowerTriangular(A)', B) .=== UnitLowerTriangular(A)'B)

                @test all(apply(*, transpose(UpperTriangular(A)), B) .=== transpose(UpperTriangular(A))B)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), B) .=== transpose(UnitUpperTriangular(A))B)
                @test all(apply(*, transpose(LowerTriangular(A)), B) .=== transpose(LowerTriangular(A))B)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), B) .=== transpose(UnitLowerTriangular(A))B)
            end

            for T in (Float64, ComplexF64)
                A = big.(randn(T,100,100))
                b = big.(randn(T,100))

                @test all(apply(*, UpperTriangular(A)', b) ≈ UpperTriangular(A)'b)
                @test all(apply(*, UnitUpperTriangular(A)', b) ≈ UnitUpperTriangular(A)'b)
                @test all(apply(*, LowerTriangular(A)', b) ≈ LowerTriangular(A)'b)
                @test all(apply(*, UnitLowerTriangular(A)', b) ≈ UnitLowerTriangular(A)'b)

                @test all(apply(*, transpose(UpperTriangular(A)), b) ≈ transpose(UpperTriangular(A))b)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), b) ≈ transpose(UnitUpperTriangular(A))b)
                @test all(apply(*, transpose(LowerTriangular(A)), b) ≈ transpose(LowerTriangular(A))b)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), b) ≈ transpose(UnitLowerTriangular(A))b)

                B = big.(randn(T,100,100))

                @test all(apply(*, UpperTriangular(A)', B) ≈ UpperTriangular(A)'B)
                @test all(apply(*, UnitUpperTriangular(A)', B) ≈ UnitUpperTriangular(A)'B)
                @test all(apply(*, LowerTriangular(A)', B) ≈ LowerTriangular(A)'B)
                @test all(apply(*, UnitLowerTriangular(A)', B) ≈ UnitLowerTriangular(A)'B)

                @test all(apply(*, transpose(UpperTriangular(A)), B) ≈ transpose(UpperTriangular(A))B)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), B) ≈ transpose(UnitUpperTriangular(A))B)
                @test all(apply(*, transpose(LowerTriangular(A)), B) ≈ transpose(LowerTriangular(A))B)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), B) ≈ transpose(UnitLowerTriangular(A))B)
            end
        end
    end

    @testset "tri Rmul" begin
        for T in (Float64, ComplexF64)
            A = randn(T, 100,100)
            B = randn(T, 100,100)
            R = Rmul(copy(A), UpperTriangular(B))
            @test size(R) == (size(R,1),size(R,2)) == (100,100)
            @test axes(R) == (axes(R,1),axes(R,2)) == (Base.OneTo(100),Base.OneTo(100))
            @test eltype(R) == T
            @test length(R) == 100^2

            @test similar(R) isa Matrix{T}
            @test similar(R,Int) isa Matrix{Int}

            @test applied(*, A, UpperTriangular(B)) isa Applied{MulStyle}
            @test similar(applied(*, A, UpperTriangular(B)), Float64) isa Matrix{Float64}

            R2 = deepcopy(R)
            Ap = applied(*, copy(A), UpperTriangular(B))
            @test all(BLAS.trmm('R', 'U', 'N', 'N', one(T), B, A) .=== apply(*, A, UpperTriangular(B)) .===
                    copyto!(similar(Ap),Ap) .=== copyto!(similar(R2), R2) .=== materialize!(R))
            @test R.A ≠ A
            @test all(BLAS.trmm('R', 'U', 'T', 'N', one(T), B, A) .=== apply(*, A, transpose(UpperTriangular(B))) .=== A*transpose(UpperTriangular(B)))
            @test all(BLAS.trmm('R', 'U', 'N', 'U', one(T), B, A) .=== apply(*, A, UnitUpperTriangular(B)) .=== A*UnitUpperTriangular(B))
            @test all(BLAS.trmm('R', 'U', 'T', 'U', one(T), B, A) .=== apply(*, A, transpose(UnitUpperTriangular(B))) .=== A*transpose(UnitUpperTriangular(B)))
            @test all(BLAS.trmm('R', 'L', 'N', 'N', one(T), B, A) .=== apply(*, A, LowerTriangular(B)) .=== A*LowerTriangular(B))
            @test all(BLAS.trmm('R', 'L', 'T', 'N', one(T), B, A) .=== apply(*, A, transpose(LowerTriangular(B))) .=== A*transpose(LowerTriangular(B)))
            @test all(BLAS.trmm('R', 'L', 'N', 'U', one(T), B, A) .=== apply(*, A, UnitLowerTriangular(B)) .=== A*UnitLowerTriangular(B))
            @test all(BLAS.trmm('R', 'L', 'T', 'U', one(T), B, A) .=== apply(*, A, transpose(UnitLowerTriangular(B))) .=== A*transpose(UnitLowerTriangular(B)))

            if T <: Complex
                @test all(BLAS.trmm('R', 'U', 'C', 'N', one(T), B, A) .=== apply(*, A, UpperTriangular(B)') .=== A*UpperTriangular(B)')
                @test all(BLAS.trmm('R', 'U', 'C', 'U', one(T), B, A) .=== apply(*, A, UnitUpperTriangular(B)') .=== A*UnitUpperTriangular(B)')
                @test all(BLAS.trmm('R', 'L', 'C', 'N', one(T), B, A) .=== apply(*, A, LowerTriangular(B)') .=== A*LowerTriangular(B)')
                @test all(BLAS.trmm('R', 'L', 'C', 'U', one(T), B, A) .=== apply(*, A, UnitLowerTriangular(B)') .=== A*UnitLowerTriangular(B)')
            end
        end

        T = Float64
        A = big.(randn(100,100))
        B = big.(randn(100,100))
        materialize!(Rmul(A,UpperTriangular(B)))
    end

    @testset "Diagonal and SymTridiagonal" begin
        A = randn(5,5)
        B = Diagonal(randn(5))
        @test MemoryLayout(typeof(B)) == DiagonalLayout{DenseColumnMajor}()
        @test ApplyStyle(*, typeof(A), typeof(B)) == MulStyle()
        @test apply(*,A,B) == A*B == materialize!(Rmul(copy(A),B))

        @test ApplyStyle(*, typeof(B), typeof(A)) == MulStyle()
        @test apply(*,B,A) == B*A

        @test ApplyStyle(*, typeof(B), typeof(B)) == MulStyle()
        @test apply(*,B,B) == B*B
        @test apply(*,B,B) isa Diagonal

        A = randn(5,5)
        B = SymTridiagonal(randn(5),randn(4))
        @test MemoryLayout(typeof(B)) == SymTridiagonalLayout{DenseColumnMajor,DenseColumnMajor}()
        @test apply(*,A,B) ≈ A*B
    end
end

@testset "Factorizations" begin
    @testset "QR" begin
        A = randn(5,3)
        b = randn(3)
        B = randn(3,3)
        Q,R = qr(A)
        @test MemoryLayout(typeof(Q)) isa QRCompactWYQLayout
        @test all(Q*b .=== apply(*,Q,b))
        @test all(Q*B .=== copyto!(similar(B,5,3),Lmul(Q,B)) .=== copyto!(similar(B,5,3),applied(*,Q,B)) .=== apply(*,Q,B))
        @test all(Q*B*b .=== apply(*,Q,B,b))
        @test_throws DimensionMismatch apply(*, Q, randn(4))
        @test_throws DimensionMismatch apply(*, Q, randn(4,3))
        dest = fill(NaN,5)
        @test copyto!(dest, applied(*,Q,b)) == Q*b

        b = randn(5)
        B = randn(5,5)
        @test all(Q*b .=== apply(*,Q,b))
        @test all(Q*b .=== apply(*,Q,b))
        @test all(Q*B .=== apply(*,Q,B))
        @test all(Q*B*b .=== apply(*,Q,B,b))

        @test MemoryLayout(typeof(Q')) isa AdjQRCompactWYQLayout
        @test all(Q'b .=== apply(*,Q',b))
        @test all(Q'B .=== apply(*,Q',B))
        @test all(*(Q',B,b) .=== apply(*,Q',B,b))
        @test_throws DimensionMismatch apply(*,Q',randn(3))
        @test_throws DimensionMismatch apply(*,Q',randn(3,3))
        dest = fill(NaN,5)
        @test all(copyto!(dest, applied(*,Q',b)) .=== Q'*b)
    end

    @testset "applied axes" begin
        T = Float64
        A = randn(5,5)
        x = randn(5)
        M =  applied(+, applied(*,T(1.0),A,x), applied(*,T(0.0),x))
        @test axes(M) == axes(M.args[1]) == (Base.OneTo(5),)
    end

    @testset "Broadcast" begin
        A = randn(5,5)
        b = randn(5)
        M = applied(*,A,b)
        @test Base.BroadcastStyle(typeof(M)) isa ApplyArrayBroadcastStyle{1}
        @test M .+ 1 ≈ A*b .+ 1

        @test Base.BroadcastStyle(ApplyArrayBroadcastStyle{1}(), Broadcast.DefaultArrayStyle{1}()) == Broadcast.DefaultArrayStyle{1}()
        @test Base.BroadcastStyle(ApplyArrayBroadcastStyle{1}(), Broadcast.DefaultArrayStyle{2}()) == Broadcast.DefaultArrayStyle{2}()
    end

    @testset "Diagonal" begin
       @test apply(*, Diagonal(Fill(2,10)), Fill(3,10)) ≡ Diagonal(Fill(2,10))  * Fill(3,10) ≡ Fill(6,10)
       @test apply(*, Diagonal(Fill(2,10)), Fill(3,10,3)) ≡ Diagonal(Fill(2,10))  * Fill(3,10,3) ≡ Fill(6,10,3)
       @test apply(*,Fill(3,10,10),Fill(3,10)) ≡ Fill(3,10,10) * Fill(3,10) ≡ Fill(90,10)
       @test apply(*, Eye(10), Ones(10)) ≡ Ones(10)
       @test apply(*, Eye(10), Eye(10)) ≡ Eye(10)
    end

    @testset "ApplyArray MulTest" begin
        A = ApplyArray(*,randn(2,2), randn(2,2))
        @test ApplyStyle(*,typeof(A),typeof(randn(2,2))) isa MulStyle
        @test ApplyArray(*,Diagonal(Fill(2,10)), Fill(3,10,10))*Fill(3,10) ≡ Fill(180,10)
        @test ApplyArray(*,Diagonal(Fill(2,10)), Fill(3,10,10))*ApplyArray(*,Diagonal(Fill(2,10)), Fill(3,10,10)) == Fill(360,10,10)
        @test A' isa ApplyArray
        @test A' ≈ Matrix(A)'
        @test transpose(A) isa ApplyArray
        @test transpose(A) ≈ transpose(Matrix(A))

        @testset "vector adjoint acts like dual vector" begin
            a = ApplyArray(*, randn(2,2), randn(2))
            @test a ≈ *(a.args...)
            @test a' isa Adjoint
            @test a'a ≈ Vector(a)'a ≈ a'Vector(a) ≈ Vector(a)'Vector(a)
            @test transpose(a) isa Transpose
            @test transpose(a)a ≈ transpose(Vector(a))a ≈ transpose(a)Vector(a) ≈ transpose(Vector(a))Vector(a)
        end
    end

    @testset "copyto!" begin
        A = ApplyArray(*,randn(2,2), randn(2,2))
        copyto!(similar(A), A)
    end
end

@testset "MulAdd" begin
    A = randn(5,5)
    B = randn(5,4)
    C = randn(5,4)
    b = randn(5)
    c = randn(5)

    M = MulAdd(2.0,A,B,3.0,C)
    @test size(M) == size(C)
    @test size(M,1) == size(C,1)
    @test size(M,2) == size(C,2)
    @test_broken size(M,3) == size(C,3)
    @test length(M) == length(C)
    @test axes(M) == axes(C)
    @test eltype(M) == Float64
    @test materialize(M) ≈ 2.0A*B + 3.0C

    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,randn(3),1.0,B))
    @test_throws DimensionMismatch materialize(MulAdd([1,2],A,B,[1,2],C))
    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(3,4)))
    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(5,5)))

    B = randn(5,5)
    C = randn(5,5)
    @test materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) ==
          materialize!(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(copy(C)))) == 2.0Diagonal(A)*Diagonal(B) + 3.0*Diagonal(C)
    @test_broken materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) isa Diagonal

    @test materialize(MulAdd(1.0, Eye(5), A, 3.0, C)) == materialize!(MulAdd(1.0, Eye(5), A, 3.0, copy(C))) == A + 3.0C
    @test materialize(MulAdd(1.0, A, Eye(5), 3.0, C)) == materialize!(MulAdd(1.0, A, Eye(5), 3.0, copy(C))) == A + 3.0C
end

@testset "MulArray" begin
  @testset "Basic" begin
        A = randn(5,5)
        M = ApplyArray(*,A,A)
        @test M[5] == M[5,1]
        @test M[6] == M[1,2]
        @test Matrix(M) ≈ A^2
        @test M^2 ≈ A^4
        @test M^2 isa Matrix{Float64}
        x = randn(5)
        @test x'M ≈ transpose(x)*M ≈ x'Matrix(M)
        @test_throws DimensionMismatch materialize(applied(*, randn(5,5), randn(4)))
        @test_throws DimensionMismatch ApplyArray(*, randn(5,5), randn(4))
    end

    @testset "Bug in getindex" begin
        M = ApplyArray(*,[1,2,3],Ones(1,20))
        @test M[1,1] == 1
        @test M[2,1] == 2
        M = Applied(*,[1 2; 3 4], [1 2; 3 4])
        @test M[1] == 7
    end

    @testset "Views" begin
        A = randn(500,500)
        b = randn(500)
        M = ApplyArray(*,A,b)

        V = view(M,2:300)
        @test MemoryLayout(typeof(V)) isa ApplyLayout{typeof(*)}
        @test arguments(V) == (view(A,2:300,Base.OneTo(500)),view(b, Base.OneTo(500)))
        @test Applied(V) isa Applied{MulStyle}
        @test ApplyArray(V) ≈ (A*b)[2:300]
        c = similar(V)
        copyto!(c,Applied(V))
        @test @allocated(copyto!(c,Applied(V))) ≤ 200
        copyto!(c, V)
        @test @allocated(copyto!(c, V)) ≤ 500
        @test all(c .=== apply(*, arguments(V)...))

        B = randn(500,500)
        M = ApplyArray(*,A,B)
        V = view(M,2:300,3:400)
        @test MemoryLayout(typeof(V)) isa ApplyLayout{typeof(*)}
        @test arguments(V) == (view(A,2:300,Base.OneTo(500)),view(B, Base.OneTo(500),3:400))
        @test Applied(V) isa Applied{MulStyle}
        c = similar(V)
        copyto!(c,Applied(V))
        @test @allocated(copyto!(c,Applied(V))) ≤ 1000
        copyto!(c, V)
        @test @allocated(copyto!(c, V)) ≤ 1000
        @test all(c .=== apply(*, arguments(V)...))
    end

    @testset "* algebra" begin
        A = ApplyArray(*,[1 2; 3 4], Vcat(Fill(1,1,3),Fill(2,1,3)))
        @test 2.0A isa ApplyArray
        @test 2.0\A isa ApplyArray
        @test A/2 isa ApplyArray
        @test (2.0A) == 2.0Array(A)
        @test (2.0\A) == 2.0\Array(A)
        @test A/2.0 == Array(A)/2.0
    end

    @testset "* trans" begin
        A = ApplyArray(*,[1 2; 3 4], Vcat(Fill(1,1,3),Fill(2,1,3)))
        V = view(A,1:2,1:2)
        @test MemoryLayout(typeof(V)) isa ApplyLayout{typeof(*)}
        @test V == ApplyArray(V)
        @test MemoryLayout(typeof(V')) isa ApplyLayout{typeof(*)}
        @test V' == ApplyArray(V')
        @test MemoryLayout(typeof(transpose(V))) isa ApplyLayout{typeof(*)}
        @test transpose(V) == ApplyArray(transpose(V))
    end

    @testset "row-vec * fix" begin
        A = ApplyArray(*,[1 2; 3 4], Vcat(Fill(1,1,3),Fill(2,1,3)))
        V = view(A, 2, 1:2)
        @test arguments(V) == ([1 2; 1 2], [3,4])
        @test ApplyArray(V) == Array(V) == A[2,1:2] == Array(A)[2,1:2]

        V = view(A, 1:2, 2)
        @test arguments(V) == ([1 2; 3 4], [1,2])
        @test ApplyArray(V) == Array(V) == A[1:2,2] == Array(A)[1:2,2]
    end

    @testset "argument type inferrence" begin
        n = 10
        L = ApplyArray(*,fill(3.0,n), ones(1,n))
        A,B = @inferred(arguments(L))
        V = view(L,1:3,1:3)
        a,b = @inferred(arguments(V))
        @test a == view(A,1:3,Base.OneTo(1))
        @test b == view(B,Base.OneTo(1),1:3)
    end

    @testset "Mul colsupport" begin
        D = ApplyArray(*,Diagonal(randn(5)),Diagonal(randn(5)),Diagonal(randn(5)))
        D̃ = Applied(*,Diagonal(randn(5)),Diagonal(randn(5)),Diagonal(randn(5)))
        @test colsupport(D,3) == colsupport(D̃,3) == 3
        @test rowsupport(D,3) == rowsupport(D̃,3) == 3
    end

    @testset "Constructors" begin
        A = randn(5,5)
        B = randn(5,5)
        M = ApplyArray(*,A, B)
        b = randn(5)
        @test ApplyArray(Mul(A,b)) ≈ lazymaterialize(Mul(A,b)) ≈ A*b
        @test ApplyArray(Lmul(UpperTriangular(A),b)) ≈ UpperTriangular(A)b
        @test ApplyArray(Rmul(A,UpperTriangular(B))) ≈ A*UpperTriangular(B)
        @test ApplyArray(Mul(M,M)) ≈ A*B*A*B
        @test ApplyArray(Mul(M,B)) ≈ A*B*B
        @test ApplyArray(Mul(A,M)) ≈ A*A*B

        @test materialize!(Lmul(applied(*,UpperTriangular(A),copy(b)))) ≈ UpperTriangular(A)b
        @test materialize!(Rmul(applied(*,copy(A),UpperTriangular(B)))) ≈ A*UpperTriangular(B)
    end

    @testset "show" begin
        A = ApplyArray(*, [1 2; 3 4], 1:2)
        @test stringmime("text/plain", A) == "(2×2 Matrix{Int64}) * (2-element UnitRange{Int64}):\n  5\n 11"
    end

    @testset "mat-mul-vec matrix indexing" begin
        A = ApplyArray(*, [1 2; 3 4], 1:2)
        # following needed for ContinuumArrays.jl
        g = [1 1; 2 2]
        V = view(A,g)
        @test MemoryLayout(V) isa UnknownLayout
        @test A[g] == layout_getindex(A, g) == [A[1] A[1]; A[2] A[2]]
    end

    @testset "permutedims" begin
        A = ApplyArray(*, [1 2; 3 4], 1:2)
        B = ApplyArray(*, [1 2; 3 4], [1 2; 3 4])
        @test permutedims(A) == A'
        @test permutedims(B) == B'
    end

    @testset "Colon indexing" begin
        A = ApplyArray(*, [1 2; 3 4], [1 2; 3 4], [1 2; 3 4])
        @test A[:,1] == Matrix(A)[:,1]
    end
end


text:
module day24

dirs = Dict("e" => :e, "se" => :se, "sw" => :sw, "w" => :w, "ne" => :ne, "nw" => :nw)
offsets = Dict(:e => (2, 0), :se => (1, 1), :sw => (-1, 1), :w => (-2, 0), :ne => (1, -1), :nw => (-1, -1))

function parseline(line)
    result = []
    while length(line) > 0
        for dir in keys(dirs)
            if startswith(line, dir)
                push!(result, dirs[dir])
                line = line[length(dir)+1:end]
                break
            end
        end
    end
    result
end

function indextocoord(index)
    coord = CartesianIndex(0, 0)
    for i in index
        coord += CartesianIndex(offsets[i]...)
    end
    coord
end

function fliptile!(tiles, index)
    coord = indextocoord(index)
    if coord ∈ keys(tiles)
        tiles[coord] = !tiles[coord]
    else
        tiles[coord] = true
    end
end

input = parseline.(readlines("day24_input.txt"))
tiles = Dict()
map(i -> fliptile!(tiles, i), input)
println("Part 1: black = $(count(values(tiles)))")

# Part 2
adjacent(coord) = [coord + CartesianIndex(o) for o in values(offsets)]

function runstep!(tiles)
    candidates = Set(keys(tiles))
    for (k, v) in tiles
        if v
            union!(candidates, Set(adjacent(k)))
        end
    end
    oldtiles = copy(tiles)
    for c in candidates
        nblack = count([oldtiles[n] for n in adjacent(c) if n ∈ keys(oldtiles)])
        if c ∈ keys(oldtiles) && oldtiles[c]
            if nblack == 0 || nblack > 2 tiles[c] = false end
        else
            if nblack == 2 tiles[c] = true end
        end
    end
end

[runstep!(tiles) for i in 1:100]
println("Part 2: black = $(count(values(tiles)))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
using Plots
using LinearAlgebra
using ProgressBars
using StaticArrays
using ArraysOfArrays

function split_vec_of_arrays(u)
    vec.(u) |>
    x ->
    VectorOfSimilarVectors(x).data |>
    transpose |>
    VectorOfSimilarVectors
end

"""ルンゲクッタ法（4次）"""
function solve_RungeKutta(dx, x₀::Vector{T}, t_span, Δt::T) where T
    t = range(t_span..., step = Δt)  # 時間軸
    x = Vector{typeof(x₀)}(undef, length(t))  # 解を格納する1次元配列

    x[1] = x₀  # 初期値
    for i in tqdm(1:length(x)-1)
        k₁ = dx(x[i])
        k₂ = dx(x[i]+k₁*Δt/2)
        k₃ = dx(x[i]+k₂*Δt/2)
        k₄ = dx(x[i]+k₃*Δt)
        x[i+1] = x[i] .+ (k₁ .+ 2k₂ .+ 2k₃ .+k₄) .* Δt/6
    end

    t, x
end

const N = 1  # セクションの数

"""Phi0のラッパー"""
function Phi0!(q::Vector{Float64}, ξ::Float64, out::Vector{Float64})
    ccall(
        (:Phi0, "o/soft_robot/derived/ikko_dake/eqs/c_so/Phi0.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], ξ, out
    )
end

"""グローバル位置"""
function Phi0(q::Vector{Float64}, ξ)
    Z = Vector{Float64}(undef, 3*N)
    Phi0!(q, ξ, Z)
    return Z
end

"""アームのディスク位置を計算"""
function Arm(q::Vector{T}, Ξ::Vector{T}) where T
    Phis = Vector{typeof(q)}(undef, length(Ξ))
    for (i, ξ) in enumerate(Ξ)
        Phis[i] = Phi0(q, ξ)
    end
    return Phis
end

"""Mのラッパー"""
function M!(q::Vector{Float64}, out::Matrix{Float64})
    ccall(
        (:M, "o/soft_robot/derived/ikko_dake/eqs/c_so/M.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], out
    )
end

"""慣性行列"""
function M(q::Vector{Float64})
    Z = Matrix{Float64}(undef, 3*N, 3*N)
    M!(q, Z)
    return Z
end

function C!(q::Vector{Float64}, q_dot::Vector{Float64}, out::Matrix{Float64})
    ccall(
        (:C, "o/soft_robot/derived/ikko_dake/eqs/c_so/C.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], q_dot[1], q_dot[2], q_dot[3], out
    )
end

"""コリオリ＋遠心力行列"""
function C(q::Vector{Float64}, q_dot::Vector{Float64})
    Z = Matrix{Float64}(undef, 3*N, 3*N)

    C!(q, q_dot, Z)

    return Z
end

"""Gのラッパー"""
function G!(q::Vector{Float64}, out::Vector{Float64})
    ccall(
        (:G, "o/soft_robot/derived/ikko_dake/eqs/c_so/G.so"),
        Cvoid,
        (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
        q[1], q[2], q[3], out
    )
end

"""重力行列"""
function G(q::Vector{Float64})
    Z = Vector{Float64}(undef, 3*N)

    G!(q, Z)

    return Z
end

const K = diagm([1700.0, 1700.0, 1700.0])
const D = diagm([110.0, 110.0, 110.0])

"""状態方程式"""
function X_dot(X::Vector{T}) where T
    q = X[1:3]
    q_dot = X[4:6]

    τ = [10.0, 0.0, 0.0]

    x1_dot = q_dot
    x2_dot = inv(M(q)) * (τ .- G(q) .- (C(q, q_dot) .+ D)*q_dot .- K*q)
    #x2_dot = M(q) \ (τ .- G(q) .- (C(q, q_dot) .+ D)*q_dot .- K*q)

    #println(norm([x1_dot; x2_dot]))
    return [x1_dot; x2_dot]

end

const Ξ = Vector(0:0.01:1)

"""1フレームを描写"""
function draw_frame(t::T, q::Vector{T}, q_dot::Vector{T}, fig_shape) where T

    arm = Arm(q, Ξ)

    x, y, z = split_vec_of_arrays(arm)
    fig = plot(
        x, y, z,
        #marker=:circle,
        aspect_ratio = 1,
        #markersize=2,
        label="arm",
        xlabel = "X[m]", ylabel = "Y[m]", zlabel = "Z[m]",
    )

    # scatter!(
    #     fig,
    #     [xd_true[1]], [xd_true[2]],
    #     label="xd_true",
    #     markershape=:star6,
    # )

    # scatter!(
    #     fig,
    #     [xd[1]], [xd[2]],
    #     label="xd_true",
    #     markershape=:star6,
    # )

    plot!(
        fig,
        xlims=(fig_shape.xl, fig_shape.xu),
        ylims=(fig_shape.yl, fig_shape.yu),
        zlims=(fig_shape.zl, fig_shape.zu),
        legend = true,
        size=(600, 600),
        title = string(round(t, digits=2)) * "[s]"
    )

    return fig
end

"""アニメ作成"""
function make_animation(t, x)
    println("アニメ作成中...")
    # 枚数決める
    #println(data.t)
    epoch_max = 100
    epoch = length(t)
    if epoch < epoch_max
        step = 1
    else
        step = div(epoch, epoch_max)
    end

    #println(step)

    x_max = 0.18
    x_min = -0.18
    y_max = 0.18
    y_min = -0.18
    z_max = 0.15
    z_min = 0.0
    max_range = max(x_max-x_min, y_max-y_min, z_max-z_min)*0.5
    x_mid = (x_max + x_min) / 2
    y_mid = (y_max + y_min) / 2
    z_mid = (z_max + z_min) / 2
    fig_shape = (
        xl = x_mid-max_range, xu = x_mid+max_range,
        yl = y_mid-max_range, yu = y_mid+max_range,
        zl = z_mid-max_range, zu = z_mid+max_range,
    )

    anim = Animation()
    @gif for i in tqdm(1:step:length(t))
        _fig = draw_frame(t[i], x[i][1:3], x[i][4:6], fig_shape)
        frame(anim, _fig)
    end

    gif(anim, "julia.gif", fps = 60)

    println("アニメ作成完了")
end

"""合ってるかテスト"""
function test()
    println("計算中...")
    TIME_SPAN = 10.0
    TIME_INTERVAL = 0.0001  # これより大きいと発散．方程式が硬すぎる?

    q = [0.0, 0.0, 0.0]
    xi = 1.0
    q_dot = [0.0, 0.0, 0.0]

    t, X = solve_RungeKutta(
        X_dot,
        [q; q_dot],
        (0.0, TIME_SPAN),
        TIME_INTERVAL
    )

    l1, l2, l3, l1_dot, l2_dot, l3_dot = split_vec_of_arrays(X)
    fig = plot(xlim=(0, TIME_SPAN))
    plot!(fig, t, l1, label="l1")
    plot!(fig, t, l2, label="l2")
    plot!(fig, t, l3, label="l3")

    fig2 = plot(xlim=(0, TIME_SPAN))
    plot!(fig2, t, l1_dot, label="l1_dot")
    plot!(fig2, t, l2_dot, label="l2_dot")
    plot!(fig2, t, l3_dot, label="l3_dot")

    fig_I = plot(fig, fig2, layout=(1, 2), size=(800, 400))
    savefig(fig_I, "julia.png")

    # アニメ制作
    make_animation(t, X)

end

#@time X = Arm([0.0, 0.0, 0.0], Ξ)
@time test()


text:
type TreeSwaptionEngine{S <: ShortRateModel} <: LatticeShortRateModelEngine{S}
  model::S
  timeSteps::Int
  common::LatticeShortRateModelEngineCommon
  # ts::Y

  # function call{S, I}(::Type{TreeSwaptionEngine}, m::S, tsteps::I)
  #   t = new{S, I, YieldTermStructure}(m, tsteps)
  #   add_observer!(m, t)
  #
  #   return t
  # end
  #
  # call{S, I}(::Type{TreeSwaptionEngine}, m::S, tsteps::I, l::LatticeShortRateModelEngineCommon) = new{S, I, YieldTermStructure}(m, tsteps, l)
  #
  # call{S, I, Y}(::Type{TreeSwaptionEngine}, m::S, tsteps::I, l::LatticeShortRateModelEngineCommon, ts::Y) = new{S, I, T, Y}(m, tsteps, l, ts)
  function TreeSwaptionEngine{S}(model::S, timeSteps::Int, common::LatticeShortRateModelEngineCommon)
    ts = new{S}(model, timeSteps, common)
    add_observer!(model, ts)

    return ts
  end

  function TreeSwaptionEngine{S}(model::S, timeSteps::Int)
    ts = new{S}(model, timeSteps)
    add_observer!(model, ts)

    return ts
  end
end

function TreeSwaptionEngine{S <: ShortRateModel}(model::S, tg::TimeGrid)
  lattice = tree(model, tg)
  return TreeSwaptionEngine{S}(model, 0, LatticeShortRateModelEngineCommon(tg, lattice))
end

TreeSwaptionEngine{S <: ShortRateModel}(model::S, timeSteps::Int) = TreeSwaptionEngine{S}(model, timeSteps)

# methods
function _calculate!(pe::TreeSwaptionEngine, swaption::Swaption)
  tsmodel = pe.model

  refDate = reference_date(tsmodel.ts)
  dc = tsmodel.ts.dc

  dSwaption = DiscretizedSwaption(swaption, refDate, dc)

  if isdefined(pe, :common)
    lattice = pe.common.lattice
  else
    times = mandatory_times(dSwaption)
    tg = TimeGrid(times, pe.timeSteps)
    lattice = tree(pe.model, tg)
    pe.common = LatticeShortRateModelEngineCommon(tg, lattice)
  end

  stoppingTimes = zeros(length(swaption.exercise.dates))
  @simd for i = 1:length(stoppingTimes)
    @inbounds stoppingTimes[i] = year_fraction(dc, refDate, swaption.exercise.dates[i])
  end

  initialize!(dSwaption, lattice.treeLattice, stoppingTimes[end])

  nextExerciseIdx = findnext(greater_than_or_equal_to, stoppingTimes, 1, 0.0)

  nextExercise = nextExerciseIdx != 0 ? stoppingTimes[nextExerciseIdx] : stoppingTimes[end]

  rollback!(dSwaption, nextExercise)
  swaption.results.value = present_value(dSwaption)
end


text:
using MultivariateStats
using LinearAlgebra
using Test
using StableRNGs
import Statistics: mean, cov
using StatsBase

"""
    UnitCovarianceEstimator

Covariance estimator that always returns a unit matrix of appropriate size.
For testing purposes only.
"""
struct UnitCovarianceEstimator <: StatsBase.CovarianceEstimator
end

function cov(::UnitCovarianceEstimator, X::AbstractMatrix; dims=1, mean=nothing)
    n = size(X, 3-dims)
    return Matrix{eltype(X)}(I, n, n)
end

@testset "LDA" begin

    rng = StableRNG(34568)

    ## LinearDiscriminant

    @assert LinearDiscriminant <: RegressionModel

    w = [1., 2., 3., 4., 5.]
    b = 2.5
    x = [4., 5., 2., 3., 1.]

    f = LinearDiscriminant(w, b)
    @test length(f) == 5
    @test dof(f) == 6
    @test coef(f) == (b, w)
    @test weights(f) == w
    @test evaluate(f, x) == 39.5
    @test evaluate(f, -x) == -34.5

    @test predict(f, x) == true
    @test predict(f, -x) == false

    X = rand(rng, 5, 8)
    Y = evaluate(f, X)
    @test size(Y) == (8,)
    for i = 1:8
        @test Y[i] ≈ evaluate(f, X[:,i])
    end
    @test predict(f, X) == (Y .> 0)

    ## prepare data

    t1 = deg2rad(30)
    t2 = deg2rad(75)

    R1 = [cos(t1) -sin(t1); sin(t1) cos(t1)]
    R2 = [cos(t2) -sin(t2); sin(t2) cos(t2)]

    n = 20
    Xp = Diagonal([1.2, 3.6]) * randn(rng, 2, n) .+ [1.0, -3.0]
    Xn = Diagonal([2.8, 1.8]) * randn(rng, 2, n) .+ [-5.0, 2.0]

    up = vec(mean(Xp, dims=2))
    un = vec(mean(Xn, dims=2))
    Cp = cov(Xp, dims=2)
    Cn = cov(Xn, dims=2)
    C = 0.5 * (Cp + Cn)

    w_gt = C \ (up - un)
    w_gt .*= (2 / (dot(w_gt, up) - dot(w_gt, un)))
    b_gt = 1.0 - dot(w_gt, up)

    @test dot(w_gt, up) + b_gt ≈ 1.0
    @test dot(w_gt, un) + b_gt ≈ -1.0

    ## LDA

    f = ldacov(C, up, un)
    @test f.w ≈ w_gt
    @test f.b ≈ b_gt

    f = ldacov(Cp, Cn, up, un)
    @test f.w ≈ w_gt
    @test f.b ≈ b_gt

    f = fit(LinearDiscriminant, Xp, Xn)
    @test f.w ≈ w_gt
    @test f.b ≈ b_gt

    covestimator = SimpleCovariance(corrected=true)
    fs = fit(LinearDiscriminant, Xp, Xn; covestimator=covestimator)
    @test fs.w ≈ w_gt
    @test fs.b ≈ b_gt

    let unitcovestimator = UnitCovarianceEstimator(),
        up = vec(mean(Xp, dims=2)),
        un = vec(mean(Xn, dims=2)),
        C = cov(unitcovestimator, Xn, dims=2)
        w_gt = C \ (up - un)
        w_gt .*= (2 / (dot(w_gt, up) - dot(w_gt, un)))
        b_gt = 1.0 - dot(w_gt, up)

        fs = fit(LinearDiscriminant, Xp, Xn; covestimator=unitcovestimator)
        @test fs.w ≈ w_gt
        @test fs.b ≈ b_gt
    end

    # input of Float32 type
    L = fit(LinearDiscriminant, [1 2 3 4f0], [5 6 7 8f0])
    @test isa(L.b, Float32)

    Lce = fit(LinearDiscriminant, [1 2 3 4f0], [5 6 7 8f0], covestimator=covestimator)
    @test isa(Lce.b, Float32)
end