text:
module LE1NAFEMS_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
using FinEtools.MeshImportModule: import_ABAQUS
using Statistics: mean
using LinearAlgebra: Symmetric, cholesky, norm

# function LE1NAFEMS_compare_meshes()
#     Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

#     ref  = 0
#     Thickness = Thick0
#     tolerance = Thickness/2^ref/300.; # Geometrical tolerance

#     fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*5, 2^ref*6, 1; orientation = :b)
#     for i=1:count(fens)
#         t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
#         fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
#     end
#     println("$((count(fens), count(fes)))")

#     output = import_ABAQUS("LE1AbaqusExport-C3D10HS-5-6-1.inp")
#     fens1, fes1 = output["fens"], output["fesets"][1]
#     println("$((count(fens1), count(fes1[1])))")

#     fens, newfes1, fes2 = mergemeshes(fens,fes, fens1,fes1[1], tolerance)
#     # fes = cat(fes2, newfes1)
#     # println("$((count(fens), count(fes)))")

#     File =  "a1.vtk"
#     vtkexportmesh(File, newfes1.conn, fens.xyz,
#     FinEtools.MeshExportModule.T10)
#     @async run(`"paraview.exe" $File`)
#     File =  "a2.vtk"
#     vtkexportmesh(File, fes2.conn, fens.xyz,
#     FinEtools.MeshExportModule.T10)
#     @async run(`"paraview.exe" $File`)

#     #
#     # fens,fes = H8block(1.0, pi/2, Thickness, 2^ref*5, 2^ref*6, 1)
#     # for i=1:count(fens)
#     #     t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
#     #     fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
#     # end
#     # println("$((count(fens), count(fes)))")
#     #
#     # fens1,fes1 = import_ABAQUS("LE1AbaqusExport-C3D8S-80-96-1.inp")
#     # println("$((count(fens1), count(fes1[1])))")
#     #
#     #  fens, newfes1, fes2 = mergemeshes(fens,fes, fens1,fes1[1], tolerance)
#     #  fes = cat(fes2, newfes1)
#     #  println("$((count(fens), count(fes)))")
#     #
#     # File =  "a.vtk"
#     # vtkexportmesh(File, fes.conn, fens.xyz,
#     #                FinEtools.MeshExportModule.H8)
#     # @async run(`"paraview.exe" $File`)

#     true

# end # LE1NAFEMS_compare_meshes

function LE1NAFEMS_MSH8()
    println("LE1NAFEMS, 3D version."        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius = 1.0*phun("m")
    Thickness = 0.1*phun("m")
    n = 2; # number of elements per side
    tolerance = 1.0/n/1000.; # Geometrical tolerance

    fens,fes = Q4block(1.0, pi/2, n, n*2)
    fens,fes  = H8extrudeQ4(fens, fes,
    1, (xyz, layer)->[xyz[1], xyz[2], (layer)*Thickness]);

    bdryfes = meshboundary(fes);
    icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
    for i=1:count(fens)
        t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
        fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
    end

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 1, 0.0)
    l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
    setebc!(u,l1,true, 3, 0.0)

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
    function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
        pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
        forceout .=    vec(p*pt/norm(pt));
        return forceout
    end
    fi = ForceIntensity(FFlt, 3, pfun);
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
    MR = DeforModelRed3D

    material = MatDeforElastIso(MR, E, nu)

    femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    # The geometry field now needs to be associated with the FEMM
    femm = associategeometry!(femm, geom)

    K = stiffness(femm, geom, u)
    K = cholesky(K)
    U = K\(F2)
    scattersysvec!(u, U[:])

    nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, 0.0],inflate=tolerance);
    thecorneru = zeros(FFlt,1,3)
    gathervalues_asmat!(u, thecorneru, nl);
    thecorneru = thecorneru/phun("mm")
    println("$(time()-t0) [s];  displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
    nodevalmethod = :averaging, reportat = :meanonly)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
    nodevalmethod = :averaging, reportat = :extrapmean)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
    nodevalmethod = :averaging, reportat = :extraptrend)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :invdistance, reportat = :meanonly)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = :extrapmean)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = :extraptrend)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")

    println("$(n), $(fld.values[nl,1][1]/phun("MPa"))")

    File =  "a.vtk"
    vtkexportmesh(File, fes.conn, geom.values,
    FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
    scalars=[("sigmay", fld.values)])
    @async run(`"paraview.exe" $File`)
    true

end # LE1NAFEMS_MSH8

function LE1NAFEMS_MSH8_convergence()
    # Example from the "Improved Stress Recovery for Mean-strain Finite Elements" paper by Sivapuram and Krysl, 2017
    println("LE1NAFEMS, 3D version. MSH8"        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

    sigyderrs = Dict{Symbol, FFltVec}()

    nnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nnodes = []
        for ref in 0:1:4
            # Thickness = Thick0
            Thickness = Thick0/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = H8block(1.0, pi/2, Thickness, 2^ref*5, 2^ref*6, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
            end

            geom = NodalField(fens.xyz)
            u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
                forceout .=    vec(p*pt/norm(pt));
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            # Note that the material object needs to be created with the proper
            # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            thecorneru = zeros(FFlt,length(nl),3)
            gathervalues_asmat!(u, thecorneru, nl);
            thecorneru = mean(thecorneru, dims=1)[1]/phun("mm")
            println("displacement =$(thecorneru) vs -0.10215 [MM]")

            fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = extrapolation)
            sigyd = mean(fld.values[nl,1], dims=1)[1]
            println("Sigma_y =$(sigyd/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

            println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            push!(nnodes, count(fes))
            push!(sigyderrs[extrapolation], (sigyd/sigma_yD - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmay", fld.values)])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "LE1NAFEMS_MSH8_convergence.CSV"
    savecsv(File, nnodes=vec(nnodes), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrmean=vec(sigyderrs[:extrapmean]))

end # LE1NAFEMS_MSH8_convergence

function LE1NAFEMS_MSH8_export()
    println("LE1NAFEMS, 3D version."        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius = 1.0*phun("m")
    Thickness = 0.1*phun("m")
    n = 64; # number of elements per side
    tolerance = 1.0/n/1000.; # Geometrical tolerance

    fens,fes = Q4block(1.0, pi/2, n, n*2)
    fens,fes  = H8extrudeQ4(fens, fes,
    1, (xyz, layer)->[xyz[1], xyz[2], (layer)*Thickness]);

    bdryfes = meshboundary(fes);
    icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
    for i=1:count(fens)
        t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
        fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
    end

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 1, 0.0)
    l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
    setebc!(u,l1,true, 3, 0.0)

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
    function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
        pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
        forceout .=    vec(p*pt/norm(pt));
        return forceout
    end
    fi = ForceIntensity(FFlt, 3, pfun);
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
    MR = DeforModelRed3D

    material = MatDeforElastIso(MR, E, nu)

    femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    # The geometry field now needs to be associated with the FEMM
    femm = associategeometry!(femm, geom)

    K = stiffness(femm, geom, u)
    K = cholesky(K)
    U = K\(F2)
    scattersysvec!(u, U[:])

    nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, 0.0],inflate=tolerance);
    thecorneru = zeros(FFlt,1,3)
    gathervalues_asmat!(u, thecorneru, nl);
    thecorneru = thecorneru/phun("mm")
    println("$(time()-t0) [s];  displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = :extraptrend)
    println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")

    println("$(n), $(fld.values[nl,1][1]/phun("MPa"))")

    File =  "a.vtk"
    vtkexportmesh(File, fes.conn, geom.values,
    FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
    scalars=[("sigmay", fld.values)])
    @async run(`"paraview.exe" $File`)
    true

end # LE1NAFEMS_MSH8_export

function LE1NAFEMS_MST10_convergence()
    # Example from the "Improved Stress Recovery for Mean-strain Finite Elements" paper by Sivapuram and Krysl, 2017
    println("LE1NAFEMS, 3D version. MST10"        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

    sigyderrs = Dict{Symbol, FFltVec}()

    nnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nnodes = []
        for ref in 0:1:4
            # Thickness = Thick0
            Thickness = Thick0/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*5, 2^ref*6, 1; orientation = :b)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
            end

            geom = NodalField(fens.xyz)
            u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
                forceout .=    vec(p*pt/norm(pt));
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            # Note that the material object needs to be created with the proper
            # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            thecorneru = zeros(FFlt,length(nl),3)
            gathervalues_asmat!(u, thecorneru, nl);
            thecorneru = mean(thecorneru, dims=1)[1]/phun("mm")
            println("displacement =$(thecorneru) vs -0.10215 [MM]")

            fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = extrapolation)
            sigyd = mean(fld.values[nl,1], dims=1)[1]
            println("Sigma_y =$(sigyd/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

            println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            push!(nnodes, count(fens))
            push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmay", fld.values)])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "LE1NAFEMS_MST10_convergence.CSV"
    savecsv(File, nnodes=vec(nnodes), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # LE1NAFEMS_MST10_convergence

function LE1NAFEMS_MST10_one()
    println("LE1NAFEMS, 3D version."        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

    sigyderrs = Dict{Symbol, FFltVec}()

    for extrapolation in [:extraptrend]
        sigyderrs[extrapolation] = FFltVec[]
        nnodes = []
        for ref in 1:1
            Thickness = Thick0
            # Thickness = Thick0/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*5, 2^ref*6, 1; orientation = :b)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
            end

            geom = NodalField(fens.xyz)
            u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
                forceout .=    vec(p*pt/norm(pt));
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            # Note that the material object needs to be created with the proper
            # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            thecorneru = zeros(FFlt,length(nl),3)
            gathervalues_asmat!(u, thecorneru, nl);
            thecorneru = mean(thecorneru, dims=1)[1]/phun("mm")
            println("displacement =$(thecorneru) vs -0.10215 [MM]")

            fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = extrapolation)
            sigyd = mean(fld.values[nl,1], dims=1)[1]
            println("Sigma_y =$(sigyd/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

            println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            push!(nnodes, count(fens))
            push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            File =  "LE1NAFEMS_MST10_a.vtk"
            vtkexportmesh(File, fes.conn, geom.values,
            FinEtools.MeshExportModule.T10; vectors=[("u", u.values)],
            scalars=[("sigmay", fld.values)])
            @async run(`"paraview.exe" $File`)
        end
    end

    #
    #
    # using DataFrames
    # using CSV
    #
    # df = DataFrame(nnodes=vec(nnodes),
    #     sigyderrtrend=vec(sigyderrs[:extraptrend]),
    #     sigyderrdefault=vec(sigyderrs[:extrapmean]))
    # File = "LE1NAFEMS_MST10_convergence.CSV"
    # CSV.write(File, df)
    # @async run(`"paraview.exe" $File`)

end # LE1NAFEMS_MST10_one

# using SymPy
# function LE1NAFEMS_MST10_single()

#     @vars ua ub uc ud
#     @vars va vb vc vd
#     @vars wa wb wc wd
#     @vars x y z
#     u(x, y, z) = ua*x +  ub*y +  uc*z +  ud
#     v(x, y, z) = va*x +  vb*y +  vc*z +  vd
#     w(x, y, z) = wa*x +  wb*y +  wc*z +  wd

#     ex = diff(u(x, y, z), x)
#     ey = diff(v(x, y, z), y)
#     ez = diff(w(x, y, z), z)
#     gxy = diff(v(x, y, z), x) + diff(u(x, y, z), y)
#     gxz = diff(w(x, y, z), x) + diff(u(x, y, z), z)
#     gyz = diff(w(x, y, z), y) + diff(v(x, y, z), z)

#     display(ex)
#     display(ey)
#     display(ez)
#     display(gxy)
#     display(gxz)
#     display(gyz)

# end # LE1NAFEMS_MST10_single

function LE1NAFEMS_MST10_stresses_nodal()
    elementtag = "MST10"
    println("LE1NAFEMS, 3D version. Element: $(elementtag)")

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

    for extrapolation in [:extraptrend :extrapmean]
        convergencestudy = FDataDict[]
        for ref in  [0, 1, 2, 3, 4]
            Thickness = Thick0
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*2, 2^ref*3, 1; orientation = :b)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
            end

            geom = NodalField(fens.xyz)
            u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
                forceout .=    vec(p*pt/norm(pt));
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            # Note that the material object needs to be created with the proper
            # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            thecorneru = zeros(FFlt,length(nl),3)
            gathervalues_asmat!(u, thecorneru, nl);
            thecorneru = mean(thecorneru, dims=1)[1]/phun("mm")
            println("displacement =$(thecorneru) vs -0.10215 [MM]")

            fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
            nodevalmethod = :averaging, reportat = extrapolation)
            sigyd = mean(fld.values[nl,1], dims=1)[1]
            println("Sigma_y =$(sigyd/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

            stressfield = fieldfromintegpoints(femm, geom, u, :Cauchy, collect(1:6);
            nodevalmethod = :averaging, reportat = extrapolation)

            # File =  "LE1NAFEMS_MST10_sigma.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            #     FinEtools.MeshExportModule.T10; vectors=[("u", u.values)],
            #     scalars=[("sig", stressfield.values)])
            # @async run(`"paraview.exe" $File`)

            push!(convergencestudy, FDataDict(
            "elementsize"=> 1.0 / (2^ref),
            "fens"=>fens,
            "fes"=>fes,
            "geom"=>geom,
            "u"=>u,
            "femm"=>femm,
            "integrationrule"=>femm.integdomain.integration_rule,
            "stressfields"=>[stressfield],
            "tolerance"=>tolerance)
            )
        end # for ref

        # File = "LE1NAFEMS_MST10_stresses_nodal_convergence_$(elementtag)_$(extrapolation)"
        # open(File * ".jls", "w") do file
        #     serialize(file, convergencestudy)
        # end

    end # for extrapolation

end # LE1NAFEMS_MST10_stresses_nodal

function LE1NAFEMS_MST10_S_convergence()
    println("LE1NAFEMS, 3D version."        )
    t0 = time()

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

    sigyderrs = Dict{Symbol, FFltVec}()

    nnodes = []
    for extrapolation in [:extraptrend :extrapmean]
        sigyderrs[extrapolation] = FFltVec[]
        nnodes = []
        for ref in 0:1:4
            Thickness = Thick0
            # Thickness = Thick0/2^ref
            tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

            fens,fes = T4block(1.0, pi/2, Thickness, 2^ref*5, 2^ref*6, 1)

            bdryfes = meshboundary(fes);
            icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
            for i=1:count(fens)
                t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
                fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
            end
            fens,fes = T4toT10(fens,fes)

            geom = NodalField(fens.xyz)
            u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

            l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 2, 0.0)
            l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
            setebc!(u,l1,true, 1, 0.0)
            l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
            setebc!(u,l1,true, 3, 0.0)

            applyebc!(u)
            numberdofs!(u)

            el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
            function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
                pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
                forceout .=    vec(p*pt/norm(pt));
                return forceout
            end
            fi = ForceIntensity(FFlt, 3, pfun);
            F2 = distribloads(el1femm, geom, u, fi, 2);

            # Note that the material object needs to be created with the proper
            # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
            MR = DeforModelRed3D

            material = MatDeforElastIso(MR, E, nu)

            femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

            # The geometry field now needs to be associated with the FEMM
            femm = associategeometry!(femm, geom)

            K = stiffness(femm, geom, u)
            K = cholesky(K)
            U = K\(F2)
            scattersysvec!(u, U[:])

            nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
            thecorneru = zeros(FFlt,length(nl),3)
            gathervalues_asmat!(u, thecorneru, nl);
            thecorneru = mean(thecorneru, dims=1)[1]/phun("mm")
            println("displacement =$(thecorneru) vs -0.10215 [MM]")

            fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2; nodevalmethod = :averaging, reportat = extrapolation)
            sigyd = mean(fld.values[nl,1], dims=1)[1]
            println("Sigma_y =$(sigyd/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

            println("$extrapolation, $(count(fes)), $(sigyd/phun("MPa"))")
            push!(nnodes, count(fens))
            push!(sigyderrs[extrapolation], abs(sigyd/sigma_yD - 1.0))
            # File =  "a.vtk"
            # vtkexportmesh(File, fes.conn, geom.values,
            # FinEtools.MeshExportModule.H8; vectors=[("u", u.values)],
            # scalars=[("sigmay", fld.values)])
            # @async run(`"paraview.exe" $File`)
        end
    end

    File = "LE1NAFEMS_MST10_S_convergence.CSV"
    savecsv(File, nnodes=vec(nnodes), sigyderrtrend=vec(sigyderrs[:extraptrend]), sigyderrdefault=vec(sigyderrs[:extrapmean]))

end # LE1NAFEMS_MST10_S_convergence

function LE1NAFEMS_T10_stresses_nodal()
    elementtag = "T10"
    println("LE1NAFEMS, 3D version. Element: $(elementtag)")

    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Thick0 = 0.1*phun("m")/2.0 # to account for the symmetry reduction

    convergencestudy = FDataDict[]
    for ref in  [0, 1, 2, 3, 4]
        Thickness = Thick0
        tolerance = Thickness/2^ref/1000.; # Geometrical tolerance

        fens,fes = T10block(1.0, pi/2, Thickness, 2^ref*2, 2^ref*3, 1; orientation = :b)

        bdryfes = meshboundary(fes);
        icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
        for i=1:count(fens)
            t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
            fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
        end

        geom = NodalField(fens.xyz)
        u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

        l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 2, 0.0)
        l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
        setebc!(u,l1,true, 1, 0.0)
        l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
        setebc!(u,l1,true, 3, 0.0)

        applyebc!(u)
        numberdofs!(u)

        el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), TriRule(3)))
        function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
            pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
            forceout .=    vec(p*pt/norm(pt));
            return forceout
        end
        fi = ForceIntensity(FFlt, 3, pfun);
        F2 = distribloads(el1femm, geom, u, fi, 2);

        # Note that the material object needs to be created with the proper
        # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
        MR = DeforModelRed3D

        material = MatDeforElastIso(MR, E, nu)

        femm = FEMMDeforLinear(MR, IntegDomain(fes, TetRule(4)), material)

        # The geometry field now needs to be associated with the FEMM
        femm = associategeometry!(femm, geom)

        K = stiffness(femm, geom, u)
        K = cholesky(K)
        U = K\(F2)
        scattersysvec!(u, U[:])

        nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, Thickness], inflate=tolerance);
        thecorneru = zeros(FFlt,length(nl),3)
        gathervalues_asmat!(u, thecorneru, nl);
        thecorneru = mean(thecorneru, dims=1)[1]/phun("mm")
        println("displacement =$(thecorneru) vs -0.10215 [MM]")

        fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
        sigyd = mean(fld.values[nl,1], dims=1)[1]
        println("Sigma_y =$(sigyd/phun("MPa")) vs $(sigma_yD/phun("MPa")) [MPa]")

        stressfield = fieldfromintegpoints(femm, geom, u, :Cauchy, collect(1:6))

        # File =  "LE1NAFEMS_T10_sigma.vtk"
        # vtkexportmesh(File, fes.conn, geom.values,
        #     FinEtools.MeshExportModule.T10; vectors=[("u", u.values)],
        #     scalars=[("sig", stressfield.values)])
        # @async run(`"paraview.exe" $File`)

        push!(convergencestudy, FDataDict(
        "elementsize"=> 1.0 / (2^ref),
        "fens"=>fens,
        "fes"=>fes,
        "geom"=>geom,
        "u"=>u,
        "femm"=>femm,
        "integrationrule"=>femm.integdomain.integration_rule,
        "stressfields"=>[stressfield],
        "tolerance"=>tolerance)
        )
    end # for ref

    # File = "LE1NAFEMS_T10_stresses_nodal_convergence_$(elementtag)"
    # open(File * ".jls", "w") do file
    #     serialize(file, convergencestudy)
    # end

end # LE1NAFEMS_T10_stresses_nodal

function allrun()
    println("#####################################################")
    println("# LE1NAFEMS_MSH8 ")
    LE1NAFEMS_MSH8()
    println("#####################################################")
    println("# LE1NAFEMS_MSH8_convergence ")
    LE1NAFEMS_MSH8_convergence()
    println("#####################################################")
    println("# LE1NAFEMS_MSH8_export ")
    LE1NAFEMS_MSH8_export()
    println("#####################################################")
    println("# LE1NAFEMS_MST10_convergence ")
    LE1NAFEMS_MST10_convergence()
    println("#####################################################")
    println("# LE1NAFEMS_MST10_one ")
    LE1NAFEMS_MST10_one()
    # println("#####################################################")
    # println("# LE1NAFEMS_MST10_single ")
    # LE1NAFEMS_MST10_single()
    println("#####################################################")
    println("# LE1NAFEMS_MST10_stresses_nodal ")
    LE1NAFEMS_MST10_stresses_nodal()
    println("#####################################################")
    println("# LE1NAFEMS_MST10_S_convergence ")
    LE1NAFEMS_MST10_S_convergence()
    println("#####################################################")
    println("# LE1NAFEMS_T10_stresses_nodal ")
    LE1NAFEMS_T10_stresses_nodal()
    return true
end # function allrun

end # module LE1NAFEMS_examples


text:
module JamesPlots

using CairoMakie
using ColorSchemes
using Colors

# https://github.com/OrdnanceSurvey/GeoDataViz-Toolkit/tree/master/Colours
const COLORS = ColorScheme(parse.(Colorant, ["#ff1f5b", "#00cd6c", "#009ade", "#af58ba", "#ffc61e", "#f28522"]))
const COLUMN_WIDTH = 240
const TWO_COLUMN_WIDTH = 481
const RESOLUTION = (COLUMN_WIDTH, COLUMN_WIDTH/MathConstants.golden)

const NICECOLORS = ColorScheme(parse.(Colorant, ["#FCE1A4","#FABF7B","#F08F6E","#E05C5C","#D12959","#AB1866","#6E005F"]))
const DIVERGINGCOLORS = ColorScheme(parse.(Colorant,["#045275","#089099","#7CCBA2","#FCDE9C","#F0746E","#DC3977","#7C1D6F"]))

export MyAxis
export save_figure
export COLORS

function MyAxis(f; kwargs...)
    my_ax = Axis(f; kwargs...)

    extra_ax = Axis(f; xaxisposition=:top, yaxisposition=:right, kwargs...)
    linkaxes!(my_ax, extra_ax)
    hidespines!(extra_ax)
    hidedecorations!(extra_ax, ticks=false, minorticks=false)

    return my_ax
end

function choose_color(ncolors)

    selectors = [
        [1],
        [1, 3],
        [1, 3, 5],
        [1, 3, 4, 5],
        [1, 3, 4, 5, 6],
        [1, 2, 3, 4, 5, 6],
    ]
    select = selectors[ncolors]
    return [COLORS[i] for i in select]
end

function get_theme(;scale, ncolors)

    color = choose_color(ncolors)

    Theme(
        resolution=RESOLUTION .* scale,
        figure_padding=1,
        rowgap=0,
        colgap=0,
        font="Times New Roman", fontsize=8 * scale,
        Axis=(
            xminorticksvisible=true,
            yminorticksvisible=true,
            xgridvisible=false,
            ygridvisible=false,
            xtickalign=1,
            ytickalign=1,
            xminortickalign=1,
            yminortickalign=1,
            spinewidth=scale,
            xtickwidth=scale,
            xminortickwidth=scale,
            ytickwidth=scale,
            yminortickwidth=scale,
            xminorticksize=2scale,
            xticksize=3scale,
            yminorticksize=2scale,
            yticksize=3scale,
            xlabelpadding=scale,
            ylabelpadding=scale,
        ),
        palette=(patchcolor=collect(color), color=color),
        Lines=(
            linewidth=1.5scale,
        ),
        Scatter=(
            strokewidth=scale,
            markersize=4scale
        ),
        Legend=(
            framevisible=true,
            colgap=5scale,
            rowgap=0,
            patchsize=(10scale, 10scale),
            patchlabelgap=3scale,
            padding=(3scale, 3scale, 2scale, 2scale),
            merge=true,
            labelsize=8scale,
        ),
        Colorbar=(
            tickalign=scale,
            ticksize=3scale,
            lip_vertical_label=true
        ),
    )
end

function save_figure(file, plot_function; scale=1, ncolors=3)
    theme = get_theme(;scale, ncolors)
    with_theme(theme) do
        save(file, plot_function(); px_per_unit=4, pt_per_unit=1)
    end
end

end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    StdNames(spec)

Standardizes column names according to given `spec`.
Default to `:upper` case specification.

# Specs

* `:upper` - Uppercase, e.g. COLUMNNAME
* `:camel` - Camelcase, e.g. ColumnName
* `:snake` - Snakecase, e.g. column_name
"""
struct StdNames <: Stateless
  spec::Symbol
end

StdNames() = StdNames(:upper)

isrevertible(::Type{StdNames}) = true

function apply(transform::StdNames, table)
  # retrieve spec
  spec = transform.spec

  # retrieve column names
  cols = Tables.columns(table)
  oldnames = string.(Tables.columnnames(cols))

  # clean column names
  cleaned = _clean.(oldnames)

  # apply spec
  spec == :camel && (names = _camel.(cleaned))
  spec == :snake && (names = _snake.(cleaned))
  spec == :upper && (names = _upper.(cleaned))

  # make names unique
  newnames = _unique(names)

  # rename transform
  rtrans = Rename(Dict(oldnames .=> newnames))
  newtable, rcache = apply(rtrans, table)

  newtable, (rtrans, rcache)
end

function revert(::StdNames, newtable, cache)
  rtrans, rcache = cache
  revert(rtrans, newtable, rcache)
end

const delim = [' ', '\t', '-', '_']

_clean(name) = filter(c -> isdigit(c) || isletter(c) || c ∈ delim, name)

function _unique(names)
  newnames = String[]
  for name in names
    n = name
    while n ∈ newnames
      n = string(n, "_")
    end
    push!(newnames, n)
  end

  newnames
end

_camel(name) = join(uppercasefirst.(split(name, delim)))

_snake(name) = join(lowercase.(split(strip(name, delim), delim)), '_')

_upper(name) = replace(uppercase(name), delim => "")


text:
# ---------------------------------------------------------
# A backup, if no PNG generation is defined, is to try to make a PDF and use FileIO to convert

_fileio_load(@nospecialize(filename::AbstractString)) =
    FileIO.load(filename::AbstractString)
_fileio_save(@nospecialize(filename::AbstractString), @nospecialize(x)) =
    FileIO.save(filename::AbstractString, x)

function _show_pdfbackends(io::IO, ::MIME"image/png", plt::Plot)
    fn = tempname()

    # first save a pdf file
    pdf(plt, fn)

    # load that pdf into a FileIO Stream
    s = _fileio_load(fn * ".pdf")

    # save a png
    pngfn = fn * ".png"
    _fileio_save(pngfn, s)

    # now write from the file
    write(io, read(open(pngfn), String))
end

const PDFBackends =
    Union{PGFPlotsBackend,PlotlyJSBackend,PyPlotBackend,InspectDRBackend,GRBackend}


text:
using RecordedArrays: getdim

@testset "MCIndices" begin
    im = MCIndices(([1, 2, 3], [1, 3, 4], [2, 4, 5]))
    @test size(im) == (3, 3, 3)
    for i in ndims(im)
        @test getdim(im, i) == im.indices[i]
    end
    for (i, ind) in enumerate(im)
        Is = Base._ind2sub(im, i)
        @test map(getindex, im.indices, Is) == ind
    end
    for i in ndims(im)
        push!(getdim(im, i), 6)
        @test im.indices[i][end] == 6
        @test length(im.indices[i]) == 4
    end
end

@testset "DOKSparseArray" begin
    s = DOKSparseArray(Dict{Tuple{},Int}(() => 1), Size())
    @test ArrayInterface.parent_type(s) == Dict{Tuple{},Int}
    @test s[] == 1
    s[] = 2
    @test get(s, (), 1) == 2
end


text:
using Quaternions
using Test
using Random
using LinearAlgebra

@testset "Quaternions Tests" begin
  @testset "basics" begin
    for T in (Int64, Float64)
      @test real(Quaternion{T}) == T
      @test quaternion(T) == Quaternion{T}
      @test quaternion(Quaternion{T}) == Quaternion{T}
    end

    @test sprint(show, quaternion(1, 0, 0, 0)) == "1 + 0im + 0jm + 0km"
    @test sprint(show, Quaternion{Int8}(0, 0, 0, typemin(Int8))) == "0 + 0im + 0jm - 128km"

    @test real(Quaternion(1, 2, 3, 4)) == 1
    @test imagi(Quaternion(1, 2, 3, 4)) == 2
    @test imagj(Quaternion(1, 2, 3, 4)) == 3
    @test imagk(Quaternion(1, 2, 3, 4)) == 4
    @test imag(Quaternion(1, 2, 3, 4)) == Quaternion(0, 2, 3, 4)
    @test conj(Quaternion(1, 2, 3, 4)) == Quaternion(1, -2, -3, -4)
  end

  @testset "unary operator on quaternion boolean" begin
    p = [true, false]
    for components in Iterators.product(p, p, p, p)
      @test +Quaternion(components...) === Quaternion(map(+,components)...)
      @test -Quaternion(components...) === Quaternion(map(-,components)...)
    end
  end

  @testset "random values" begin
    for T in (Float16, Float32, Float64)
      @test typeof(rand(Quaternion{T})) == Quaternion{T}
      @test typeof(randn(Quaternion{T})) == Quaternion{T}
    end
    for T in (Int8, Int16, Int32, Int64, Int128)
      @test typeof(rand(Quaternion{T})) == Quaternion{T}
    end
  end

  @testset "equivalent complex matrix" begin
    x = rand(QuaternionF64, 100,100)
    @test isequal(Quaternions.equiv(Quaternions.equiv(x)), x)
  end

  @testset "linear algebra" begin
    N = 100
    M = 200
    x = rand(QuaternionF64, N, N)
    y = rand(QuaternionF64, M, N)
    atol = 1e-12

    x_inv = inv(x)
    @test x_inv isa typeof(x)
    x_mul = x_inv * x
    @test isapprox(real(x_mul), I(N); atol = atol)
    @test isapprox(imagi(x_mul), zeros(N, N); atol = atol)
    @test isapprox(imagj(x_mul), zeros(N, N); atol = atol)
    @test isapprox(imagk(x_mul), zeros(N, N); atol = atol)

    y_pinv = pinv(y)
    @test y_pinv isa typeof(y)
    y_mul = y_pinv * y
    @test isapprox(real(y_mul), I(N); atol = atol)
    @test isapprox(imagi(y_mul), zeros(N, N); atol = atol)
    @test isapprox(imagj(y_mul), zeros(N, N); atol = atol)
    @test isapprox(imagk(y_mul), zeros(N, N); atol = atol)
  end

  @testset "math functions" begin
    for T in (Float16, Float32, Float64)
      x = rand(Quaternion{T})
      @test exp(x) isa typeof(x)
      @test log(x) isa typeof(x)
      @test cos(x) isa typeof(x)
      @test sin(x) isa typeof(x)
      @test tan(x) isa typeof(x)
      @test cosh(x) isa typeof(x)
      @test sinh(x) isa typeof(x)
      @test tanh(x) isa typeof(x)
      @test sqrt(x) isa typeof(x)
      @test x^x isa typeof(x)
      @test x^2 isa typeof(x)
      @test acos(x) isa typeof(x)
      @test asin(x) isa typeof(x)
      @test atan(x) isa typeof(x)
      @test norm(x) isa T
    end
    @test round(Quaternion(0.1, 0.1, 0.1, 0.1)) == Quaternion(0.0)
    @test round(Quaternion(0.9, 0.9, 0.9, 0.9)) == Quaternion(1.0, 1.0, 1.0, 1.0)
  end

  @testset "type handling" begin
    for T in (Int8, Int16, Int32, Int64, Int128)
      @test float(Quaternion{T}) == Quaternion{float(T)}
    end
    for T in (QuaternionF16, QuaternionF32, QuaternionF64)
      @test float(T) == T
    end
  end
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "../.."))
Pkg.instantiate()

# ## visualizer
vis = Visualizer()
open(vis)

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra
using FiniteDiff
using DojoEnvironments

# ## scripts
include(joinpath(module_dir(), "examples/policy/methods/continuation.jl"))
include(joinpath(module_dir(), "examples/policy/methods/tvlqr.jl"))
include(joinpath(module_dir(), "DojoEnvironments/src",
    "quadruped/methods/template.jl"))

################################################################################
# ## system
################################################################################
gravity = -9.81
timestep = 0.02
friction_coefficient = 0.8
damper = 0.5
spring = 1.0
env = get_environment(:quadruped,
    representation=:minimal,
    timestep=timestep,
    contact_body=false,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring,
    infeasible_control=true,
    vis=vis)

# ## dimensions
n = env.num_states
m = env.num_inputs
nu_infeasible = 6

################################################################################
# ## simulation test
################################################################################
mech = get_mechanism(:quadruped,
    contact_body=false,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring)

initialize!(mech, :quadruped, body_position=[0,0,0.0])
u_hover = [0.02;0;1.9; 0;0;0; zeros(12)]
function ctrl!(m, k; u=u_hover)
    nu = input_dimension(m)
    set_input!(m, SVector{nu}(u))
end

Main.@elapsed storage = simulate!(mech, 0.6, ctrl!,
    opts=SolverOptions(rtol=1e-5, btol=1e-4, undercut=5.0),
    )
Dojo.visualize(mech, storage, vis=env.vis)

################################################################################
# ## reference trajectory
################################################################################
N = 1
initialize!(env.mechanism, :quadruped)
x1 = quadruped_trajectory(env.mechanism,
    r=0.08,
    z=0.29;
    Δx=-0.02,
    Δfront=0.10,
    width_scale=0.0,
    height_scale=1.0,
    N=1,
    Ncycles=N)[1]

thrust = -gravity * [-zeros(4); 2.69ones(8); zeros(17); 1.9ones(8)]
positions = [0.29]
velocities = [0.00]
xref = [deepcopy(x1)]
for i = 1:length(thrust)
    a = gravity + thrust[i]
    vel = velocities[end] + timestep * a
    pos = positions[end] + timestep * vel
    push!(positions, pos)
    push!(velocities, vel)
    # reference trajectory
    x = deepcopy(x1)
    x[3] = positions[end]
    x[9] = velocities[end]
    push!(xref, x)
end
plot(positions)

DojoEnvironments.visualize(env, xref)
# ## horizon
T = length(xref)

################################################################################
# ## ILQR problem
################################################################################
# ## model
dyn = IterativeLQR.Dynamics(
    (y, x, u, w) -> dynamics(y, env, x, u, w),
    (dx, x, u, w) -> dynamics_jacobian_state(dx, env, x, u, w),
    (du, x, u, w) -> dynamics_jacobian_input(du, env, x, u, w),
    n, n, m)

model = [dyn for t = 1:T-1]

# ## rollout
x1 = deepcopy(xref[1])
ū = [u_hover for t = 1:T-1]

x̄ = IterativeLQR.rollout(model, x1, ū)
DojoEnvironments.visualize(env, x̄)

# ## objective
############################################################################
qt = [0.3; 0.05; 5;
    5e-2 * ones(3);
    1e-3 * ones(3);
    1e-3 * ones(3);
    fill([2, 1e-3], 12)...]
ots = [(x, u, w) -> transpose(x - xref[t]) * Diagonal(timestep * qt) * (x - xref[t]) +
    transpose(u) * Diagonal(timestep * 0.3 * ones(m)) * u for t = 1:T-1]
oT = (x, u, w) -> transpose(x - xref[end]) * Diagonal(timestep * qt) * (x - xref[end])

cts = [IterativeLQR.Cost(ot, n, m) for ot in ots]
cT = IterativeLQR.Cost(oT, n, 0)
obj = [cts..., cT]

# ## constraints
############################################################################
ul = -1.0 * 1e-3*ones(nu_infeasible)
uu = +1.0 * 1e-3*ones(nu_infeasible)

function contt(x, u, w)
    [
        1e-1 * (ul - u[1:nu_infeasible]);
        1e-1 * (u[1:nu_infeasible] - uu);
    ]
end

function cont_mid(x, u, w)
    [
        1e-1 * (ul - u[1:nu_infeasible]);
        1e-1 * (u[1:nu_infeasible] - uu);
        1e-1 * (x[3:3] - xref[20][3:3])
    ]
end

function goal(x, u, w)
    Δ = 1e-1 * (x - xref[end])[[1:6;13:2:36]]
    return Δ
end

con_policyt = IterativeLQR.Constraint(contt, n, m, indices_inequality=collect(1:2nu_infeasible))
con_policy_mid = IterativeLQR.Constraint(cont_mid, n, m, indices_inequality=collect(1:2nu_infeasible))
con_policyT = IterativeLQR.Constraint(goal, n, 0)

cons = [[con_policyt for t = 1:19]..., con_policy_mid, [con_policyt for t = 1:17]..., con_policyT]

# ## solver
options = Options(line_search=:armijo,
        max_iterations=50,
        max_dual_updates=12,
        min_step_size=1e-2,
        objective_tolerance=1e-3,
        lagrangian_gradient_tolerance=1e-3,
        constraint_tolerance=1e-3,
        initial_constraint_penalty=1e-1,
        scaling_penalty=10.0,
        max_penalty=1e4,
        verbose=true)

s = IterativeLQR.Solver(model, obj, cons, options=options)

IterativeLQR.initialize_controls!(s, ū)
IterativeLQR.initialize_states!(s, x̄)

# ## solve
local_callback!(solver::IterativeLQR.Solver) = continuation_callback!(solver, env)
reset!(env)
@time IterativeLQR.constrained_ilqr_solve!(s, augmented_lagrangian_callback! = local_callback!)

# ## solution
x_sol, u_sol = IterativeLQR.get_trajectory(s)
K_sol = deepcopy(s.policy.K)

# ## visualize
x_view = [[x_sol[1] for t = 1:15]..., x_sol..., [x_sol[end] for t = 1:15]...]
DojoEnvironments.visualize(env, x_view)


text:
import Base: fetch, map, put!, push!, take!, isopen, close
using Sockets

include("constants.jl")
include("crc32.jl")
include("protocol.jl")
include("io.jl")
include("channels.jl")
include("client.jl")
include("requests.jl")


text:
function gen_default_preamble()
    preamble = []
    push!(preamble, "% Default preamble")
    push!(preamble, "\\usepackage{pgfplots}")
    push!(preamble, "\\pgfplotsset{compat=newest}")
    push!(preamble, "\\usepackage{lmodern}")
    push!(preamble, "\\usepackage[T1]{fontenc}")
    push!(preamble,"\\usepackage{sansmath}")
end

CUSTOM_PREAMBLE = String[]

push!(CUSTOM_PREAMBLE, gen_default_preamble()...)


text:
using BinDeps

@BinDeps.setup

libpq = library_dependency("libpq", aliases=["libpq5"])

major, minor = (9,4)
branch = "REL$(major)_$(minor)_STABLE"

postgresql_srcurl = "http://git.postgresql.org/gitweb/?p=postgresql.git;a=snapshot;h=refs/heads/$branch;sf=tgz"
postgresql_giturl = "https://github.com/postgres/postgres"  # GitHub mirror

postgresql_srcdir = joinpath(BinDeps.srcdir(libpq), "postgresql")
postgresql_usrdir = BinDeps.usrdir(libpq)

provides(AptGet, "libpq5", libpq)
provides(Yum, "libpq5", libpq)
provides(Yum, "postgresql-libs", libpq)
provides(Pacman, "postgresql-libs", libpq)

@osx_only begin
    using Homebrew
    provides(Homebrew.HB, "postgresql", libpq, os=:Darwin)
end

provides(SimpleBuild, (@build_steps begin
    BinDeps.DirectoryRule(postgresql_srcdir,  # this isn't exported for some reason
        `git clone $postgresql_giturl $postgresql_srcdir --depth 1 --branch $branch`
    )
    FileRule(joinpath(postgresql_srcdir, "GNUMakefile"), @build_steps begin
        ChangeDirectory(postgresql_srcdir)
        `./configure --prefix=$postgresql_usrdir`
    end)
    @build_steps begin
        ChangeDirectory(joinpath(postgresql_srcdir, "src", "interfaces", "libpq"))
        MakeTargets()
        MakeTargets("install")
    end
end), libpq)

@BinDeps.install Dict(:libpq => :libpq)

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
function test_algorithm(algo::BanditAlgorithm, arms::Vector, num_sims::Int64, horizon::Int64)
  chosen_arms = zeros(Int64, num_sims * horizon)
  rewards = zeros(num_sims * horizon)
  cumulative_rewards = zeros(num_sims * horizon)
  sim_nums = zeros(num_sims * horizon)
  times = zeros(num_sims * horizon)

  for sim in 1:num_sims
    initialize(algo, length(arms))

    for t = 1:horizon
      index = (sim - 1) * horizon + t
      sim_nums[index] = sim
      times[index] = t

      chosen_arm = select_arm(algo)
      chosen_arms[index] = chosen_arm

      reward = draw(arms[chosen_arms[index]])
      rewards[index] = reward

      if t == 1
        cumulative_rewards[index] = reward
      else
        cumulative_rewards[index] = cumulative_rewards[index - 1] + reward
      end

      update(algo, chosen_arm, reward)
    end
  end

  hcat(sim_nums, times, chosen_arms, rewards, cumulative_rewards)
end

function test_algorithm(algo::BanditAlgorithm, arms::Vector, num_sims::Int64, horizon::Int64)
  chosen_arms = zeros(Int64, num_sims * horizon)
  rewards = zeros(num_sims * horizon)
  cumulative_rewards = zeros(num_sims * horizon)
  sim_nums = zeros(num_sims * horizon)
  times = zeros(num_sims * horizon)

  for sim in 1:num_sims
    initialize(algo, length(arms))

    for t = 1:horizon
      index = (sim - 1) * horizon + t
      sim_nums[index] = sim
      times[index] = t

      chosen_arm = select_arm(algo)
      chosen_arms[index] = chosen_arm

      reward = draw(arms[chosen_arms[index]])
      rewards[index] = reward

      if t == 1
        cumulative_rewards[index] = reward
      else
        cumulative_rewards[index] = cumulative_rewards[index - 1] + reward
      end

      update(algo, chosen_arm, reward)
    end
  end

  hcat(sim_nums, times, chosen_arms, rewards, cumulative_rewards)
end

function test_algorithm(algo::BanditAlgorithm, env::BanditEnvironment, num_sims::Int64, horizon::Int64)
  chosen_arms = zeros(Int64, num_sims * horizon)
  rewards = zeros(num_sims * horizon)
  cumulative_rewards = zeros(num_sims * horizon)
  sim_nums = zeros(num_sims * horizon)
  times = zeros(num_sims * horizon)

  for sim in 1:num_sims
    initialize(algo, n_arms(env))

    for t = 1:horizon
      index = (sim - 1) * horizon + t
      sim_nums[index] = sim
      times[index] = t

      chosen_arm = select_arm(algo)
      chosen_arms[index] = chosen_arm

      current_arms = arms(env, t)
      reward = draw(current_arms[chosen_arms[index]])
      rewards[index] = reward

      if t == 1
        cumulative_rewards[index] = reward
      else
        cumulative_rewards[index] = cumulative_rewards[index - 1] + reward
      end

      update(algo, chosen_arm, reward)
    end
  end

  hcat(sim_nums, times, chosen_arms, rewards, cumulative_rewards)
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


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

using Markdown
using InteractiveUtils

# ╔═╡ 4ea41550-9313-11eb-1cdf-df6c5527f80f
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 4c360009-f01e-45bb-b458-29b15f8cac8f
begin
	using LinearAlgebra, SpecialMatrices, Statistics
	import Random
	Random.seed!(125)
	x=rand(-9:9,5)
end

# ╔═╡ 7d3ed16e-50ea-4ede-b968-e32baa56f76f
md"""
# Eigenvalue Decomposition - Perturbation Theory

__Prerequisites__

The reader should be familiar with basic linear algebra concepts and facts about eigenvalue decomposition.

__Competences__

The reader should be able to understand and check the facts about perturbations of eigenvalues and eigenvectors.

"""

# ╔═╡ da0d7094-fa3e-4f1e-a6f1-9d1b9316f65a
md"""
# Norms

In order to measure changes, we need to define norms. For more details and the proofs of the Facts below, see
[R. Byers and B. N. Datta, Vector and Matrix Norms, Error Analysis, Efficiency, and Stability, pp. 50.1-50.24](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

## Definitions

__Norm__ on a vector space $X$ is a real-valued function $\| \phantom{x} \| : X\to \mathbb{R}$ with the following properties:

1. __Positive definiteness.__ $\| x\|\geq 0$ and $\|x\|=0$ if and only if $x$ is the zero vector.
2. __Homogeneity.__ $\| \lambda x\|=|\lambda| \|x\|$
3. __Triangle inequality.__ $\| x+y\| \leq \|x\|+\|y\|$

Commonly encountered vector norms for $x\in\mathbb{C}^n$ are:

* __Hölder norm__ or $p$-__norm__: for $p\geq 1$, $\|x\|_p=\big(|x_1|^p+|x_2|^p+\cdots |x_n|^p)^{1/p}$,
* __Sum norm__ or $1$-__norm__: $\|x\|_1=|x_1|+|x_2|+\cdots |x_n|$,
* __Euclidean norm__ or $2$-__norm__: $\|x\|_2=\sqrt{|x_1|^2+|x_2|^2+\cdots |x_n|^2}$,
* __Sup-norm__ or $\infty$-__norm__: $\|x\|_\infty = \max\limits_{i=1,\ldots,n} |x_i|$.

Vector norm is __absolute__ if $\||x|\|=\|x\|$.

Vector norm is __monotone__ if $|x|\leq |y|$ implies $\|x\|\leq \|y\|$.

From every vector norm we can derive a corresponding __induced__ matrix norm (also, __operator norm__ or __natural norm__):

$$\|A\| = \max\limits_{x\neq 0} \frac{\|Ax\|}{\|x\|}=\max\limits_{\|x\|=1} \|Ax\|.$$

For matrix $A\in\mathbb{C}^{m\times n}$ we define:

* __Maximum absolute column sum norm__: $\|A\|_1=\max\limits_{1\leq j \leq n} \sum_{i=1}^m |a_{ij}|$,
* __Spectral norm__: $\|A\|_2=\sqrt{\rho(A^*A)}=\sigma_{\max}(A)$  (largest singular value of $A$),
* __Maximum absolute row sum norm__: $\|A\|_{\infty}=\max\limits_{1\leq i \leq m} \sum_{j=1}^n |a_{ij}|$,
* __Euclidean norm__ or __Frobenius norm__:
$\|A\|_F =\sqrt{\sum_{i,j} |a_{ij}|^2}=\sqrt{\mathop{\mathrm{tr}}(A^*A)}$.

Matrix norm is __consistent__ if $\|A\cdot B\|\leq \|A\| \cdot \| B\|$, where $A$ and $B$ are compatible for matrix multiplication.

Matrix norm is __absolute__ if $\||A|\|=\|A\|$.
"""

# ╔═╡ f8ba4230-9313-11eb-170c-ab58597826f8
md"
### Examples
"

# ╔═╡ 7bd6f898-b86c-4670-96e8-9016666ee14f
norm(x,1), norm(x), norm(x,Inf)

# ╔═╡ 8e1ced3a-d23d-4f5f-8c65-dc13ab392504
A=rand(-4:4,5,4)

# ╔═╡ d9ef7b4b-4954-464a-96ac-e3f49800e0b7
norm(A,1), norm(A), norm(A,2), norm(A,Inf),
opnorm(A), opnorm(A,1), opnorm(A,Inf)

# ╔═╡ 90fcd9f5-d98c-4db6-bf30-21e377209681
md"""
## Facts

1.  $\|x\|_1$, $\|x\|_2$, $\|x\|_\infty$ and $\|x\|_p$ are absolute and monotone vector norms.
2. A vector norm is absolute iff it is monotone.
3. __Convergence.__ $x_k\to x_*$ iff for any vector norm $\|x_k-x_*\|\to 0$.
4. Any two vector norms are equivalent in the sense that, for all $x$ and some $\alpha,\beta>0$

$$\alpha \|x\|_\mu \leq \|x\|_\nu \leq \beta \|x\|_\mu.$$

In particular:

   *  $\|x\|_2 \leq \|x\|_1\leq \sqrt{n}\|x\|_2$,
   *  $\|x\|_\infty \leq \|x\|_2\leq \sqrt{n}\|x\|_\infty$,
   *  $\|x\|_\infty \leq \|x\|_1\leq n\|x\|_\infty$.

5. __Cauchy-Schwartz inequality.__ $|x^*y|\leq \|x\|_2\|y\|_2$.
6. __Hölder inequality.__ If $p,q\geq 1$ and $\displaystyle\frac{1}{p}+\frac{1}{q}=1$, then

$$|x^*y|\leq \|x\|_p\|y\|_q.$$

7.  $\|A\|_1$, $\|A\|_2$ and $\|A\|_\infty$ are induced by the corresponding vector norms.
8.  $\|A\|_F$ is not an induced norm.
9.  $\|A\|_1$, $\|A\|_2$, $\|A\|_\infty$ and $\|A\|_F$ are consistent.
10.  $\|A\|_1$, $\|A\|_\infty$ and $\|A\|_F$ are absolute. However, $\||A|\|_2\neq \|A\|_2$.
11. Any two matrix norms are equivalent in the sense that, for all $A$ and some $\alpha,\beta>0$

$$\alpha \|A\|_\mu \leq \|A\|_\nu \leq \beta \|A\|_\mu.$$

In particular:
   *  $\frac{1}{\sqrt{n}}\|A\|_\infty \leq \|A\|_2\leq \sqrt{m}\|A\|_\infty$,
   *  $\|A\|_2 \leq \|A\|_F\leq \sqrt{n}\|A\|_2$,
   *  $\frac{1}{\sqrt{m}}\|A\|_1 \leq \|A\|_2\leq \sqrt{n}\|A\|_1$.

12.  $\|A\|_2\leq \sqrt{\|A\|_1 \|A\|_\infty}$.
13.  $\|AB\|_F\leq \|A\|_F\|B\|_2$ and $\|AB\|_F\leq \|A\|_2\|B\|_F$.
14. If $A=xy^*$, then $\|A\|_2=\|A\|_F=\|x\|_2\|y\|_2$.
15.  $\|A^*\|_2=\|A\|_2$ and $\|A^*\|_F=\|A\|_F$.
16. For a unitary matrix $U$ of compatible dimension,

$$\|AU\|_2=\|A\|_2,\quad \|AU\|_F=\|A\|_F,\quad
\|UA\|_2=\|A\|_2,\quad  \|UA\|_F=\|A\|_F.$$

17. For $A$ square, $\rho(A)\leq\|A\|$.
18. For $A$ square, $A^k\to 0$ iff $>\rho(A)<1$.
"""

# ╔═╡ 8261b5ce-9315-11eb-2378-af63c3eb4c02
md"
### Random matrix
"

# ╔═╡ bb6eee61-c3ee-4b6a-ab44-a638fba776b8
A

# ╔═╡ 060b8d2f-e6f3-4625-9cf6-d4cc504145fa
# Absolute norms. Spectral norm is not absolute.
opnorm(A,1), opnorm(abs.(A),1), norm(A,1), opnorm(A,Inf), opnorm(abs.(A),Inf), norm(A), norm(abs.(A)),  opnorm(A),opnorm(abs.(A))

# ╔═╡ f55a3c60-9315-11eb-027a-dd008031f6f9
m,n=size(A)

# ╔═╡ c4a1b888-9497-4eb4-aa1c-23671123f426
# Equivalence of norms
opnorm(A,Inf)\ √n,opnorm(A), √m*opnorm(A,Inf)

# ╔═╡ b46cbb6f-36f4-4a66-b970-8fe9b0506b6a
opnorm(A), norm(A), √n*opnorm(A)

# ╔═╡ 254b2df4-f57e-4284-a735-6fd5d4546c52
opnorm(A,1)\ √m,opnorm(A), √n*opnorm(A,1)

# ╔═╡ 09fd0923-1987-4a2e-be28-78451aa2c11a
# Fact 12
opnorm(A), √(opnorm(A,1)*opnorm(A,Inf))

# ╔═╡ dcba6892-b90e-41ad-8d58-53e72d40e882
begin
	# Fact 13
	B=rand(n,rand(1:9))
	norm(A*B), norm(A)*opnorm(B), opnorm(A)*norm(B)
end

# ╔═╡ c765c02e-9316-11eb-146c-ff5d3ba81323
x

# ╔═╡ d948741f-6748-44d7-ae8c-cf074f7549bd
begin
	# Fact 14
	y=rand(ComplexF64,10)
	opnorm(x*y'), norm(x*y'), norm(x)*norm(y)
end

# ╔═╡ ec0746a1-41bd-4d61-9deb-8e227106e592
begin
	# Fact 15
	A₁=rand(-4:4,7,5)+im*rand(-4:4,7,5)
	opnorm(A₁), opnorm(A₁'), norm(A₁), norm(A₁')
end

# ╔═╡ d92799f7-1813-450b-8353-9d03078f69f0
# Unitary invariance - generate random unitary matrix U
U=qr(rand(ComplexF64,size(A₁)));

# ╔═╡ 0d58d4c4-16c6-4a5b-bdbd-f3f599d0d648
opnorm(A₁), opnorm(U.Q*A₁), norm(A₁), norm(U.Q*A₁)

# ╔═╡ b1a90c59-8c3e-4530-83ef-6b3e6c6fa687
begin
	# Spectral radius
	A₂=rand(ComplexF64,7,7)
	maximum(abs,eigvals(A₂)), opnorm(A₂,Inf), opnorm(A₂,1), opnorm(A₂), norm(A₂)
end

# ╔═╡ 31e73f82-93e7-448c-940b-5c9f5e573e5c
begin
	# Fact 18
	B₂=A₂/(maximum(abs,eigvals(A₂))+2)
	maximum(abs,eigvals(B₂)), norm(B₂^100)
end

# ╔═╡ 909033d9-bb27-4f0e-b07d-574f7fcee700
md"""
# Errors and condition numbers

We want to answer the question:

__How much the value of a function changes with respect to the change of its argument?__

## Definitions

For function $f(x)$ and argument $x$, the __absolute error__ with respect to the __perturbation__ of the argument
$\delta x$ is

$$
\| f(x+\delta x)-f(x)\| = \frac{\| f(x+\delta x)-f(x)\|}{\| \delta x \|} \|\delta x\| \equiv \kappa \|\delta x\|.$$

The  __condition__ or  __condition number__ $\kappa$ tells how much does the perturbation of the argument increase. (Its form resembles derivative.)

Similarly, the __relative error__ with respect to the relative perturbation of the argument is

$$
\frac{\| f(x+\delta x)-f(x)\|}{\| f(x) \|}= \frac{\| f(x+\delta x)-f(x)\|\cdot  \|x\| }{\|\delta x\| \cdot\| f(x)\|}
\cdot \frac{\|\delta x\|}{\|x\|} \equiv \kappa_{rel} \frac{\|\delta x\|}{\|x\|}.$$
"""

# ╔═╡ 7ef162c0-5c80-4e96-bcc3-125b302b1d5e
md"""
# Peturbation bounds

## Definitions

Let $A\in\mathbb{C}^{n\times n}$.

Pair $(\lambda,x)\in\mathbb{C}\times\mathbb{C}^{n\times n}$ is an __eigenpair__ of $A$ if $x\neq 0$ and $Ax=\lambda x$.

Triplet $(y,\lambda,x)\in\times\mathbb{C}^{n}\times\mathbb{C}\times\mathbb{C}^{n}$ is an __eigentriplet__ of $A$ if $x,y\neq 0$ and $Ax=\lambda x$ and $y^*A=\lambda y^*$.

__Eigenvalue matrix__ is a diagonal matrix $\Lambda=\mathop{\mathrm{diag}}(\lambda_1,\lambda_2,\ldots,\lambda_n)$.

If all eigenvalues are real, they can be increasingly ordered. $\Lambda^\uparrow$ is the eigenvalue matrix of increasingly ordered eigenvalues.

 $\tau$ is a __permutation__ of $\{1,2,\ldots,n\}$.

 $\tilde A=A+\Delta A$ is a __perturbed matrix__, where $\Delta A$ is __perturbation__. $(\tilde \lambda,\tilde x)$ are the eigenpairs of $\tilde A$.

__Condition number__ of a nonsingular matrix $X$ is $\kappa(X)=\|X\| \|X^{-1}\|$.

Let $X,Y\in\mathbb{C}^{n\times k}$ with $\mathop{\mathrm{rank}}(X)=\mathop{\mathrm{rank}}(Y)=k$. The __canonical angles__ between their column spaces, $\theta_i$, are defined by $\cos \theta_i=\sigma_i$, where $\sigma_i$ are the singular values of the matrix

$$(Y^*Y)^{-1/2}Y^*X(X^*X)^{-1/2}.$$

The __canonical angle matrix__ between $X$ and $Y$ is

$$\Theta(X,Y)=\mathop{\mathrm{diag}}(\theta_1,\theta_2,\ldots,\theta_k).$$

"""

# ╔═╡ 3c8f71a5-0c2f-43a0-bf71-db9df707f84e
md"""
## Facts

Bounds become more strict as matrices have more structure.
Many bounds have versions in spectral norm and Frobenius norm.
For more details and the proofs of the Facts below, see
[R.-C. Li, Matrix Perturbation Theory, pp 21.1-21.20](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897), and the references therein.

1. There exists $\tau$ such that

$$\|\Lambda- \tilde\Lambda_\tau\|_2\leq 4(\|A\|_2+\|\tilde A\|_2)^{1-1/n}\|\Delta A\|_2^{1/n}.$$

2. __First-order perturbation bounds.__ Let $(y,\lambda,x)$ be an eigentriplet of a simple $\lambda$. $\Delta A$ changes $\lambda$ to $\tilde\lambda=\lambda+ \delta\lambda$, where _(for proof see GVL p. 359)_

$$\delta\lambda=\frac{y^*(\Delta A)x}{y^*x}+O(\|\Delta A\|_2^2).$$

3. Let $\mu$ be a semisimple eigenvalue of $A$ with multiplicitiy $k$, and let $X,Y\in \mathbb{C}^{n\times k}$ be the matrices of the corresponding right and left eigenvectors, that is, $AX=\lambda X$ and $Y^*A=\lambda Y^*$, such that $Y^*X=I_k$. $\Delta A$ changes the $k$ copies of $\mu$ to $\tilde \mu=\mu+\delta\mu_i$, where $\delta\mu_i$ are the eigenvalues of $Y^*(\Delta A) X$ up to $O(\|\Delta A\|_2^2)$.

4. Perturbations and the inverse: if $\|A\|_p<1$, then $I-A$ is nonsingular and _(for proof see GVL p. 74)_

$$(I-A)^{-1}=\sum\limits_{k=0}^\infty A^k,$$

$$\|(I-A)^{-1}\|_p \leq \frac{1}{1-\|A\|_p},\qquad
\|(I-A)^{-1}-I\|_p \leq \frac{\|A\|_p}{1-\|A\|_p}.$$

5. __Geršgorin Circle Theorem.__ If $X^{-1} A X=D+F$, where $D=\mathop{\mathrm{diag}}(d_1,\ldots,d_n)$ and $F$ has zero diagonal entries, then _(for proof see GVL p. 357)_

$$\sigma(A) \subseteq \bigcup\limits_{i=1}^n D_i,$$

where

$$D_i=\big\{z\in\mathbb{C} : |z-d_i| \leq \sum\limits_{j=1}^n |f_{ij}| \big\}.$$

Moreover, by continuity, if a connected component of $D$ consists of $k$ circles, it contains $k$ eigenvalues.

6. __Bauer-Fike Theorem.__ If $A$ is diagonalizable and $A=X\Lambda X^{-1}$ is its eigenvalue decomposition, then _(for proof see GVL p. 357)_

$$\max_i\min_j |\tilde \lambda_i -
\lambda_j|\leq \|X^{-1}(\Delta A)X\|_p\leq \kappa_p(X)\|\Delta A\|_p.$$

7. If $A$ and $\tilde A$ are diagonalizable, there exists $\tau$ such that

$$\|\Lambda-\tilde\Lambda_\tau\|_F\leq \sqrt{\kappa_2(X)\kappa_2(\tilde X)}\|\Delta A\|_F.$$

If $\Lambda$ and  $\tilde\Lambda$ are real, then

$$\|\Lambda^\uparrow-\tilde\Lambda^\uparrow\|_{2,F} \leq \sqrt{\kappa_2(X)\kappa_2(\tilde X)}\|\Delta A\|_{2,F}.$$

8. If $A$ is normal, there exists $\tau$ such that $\|\Lambda-\tilde\Lambda_\tau\|_F\leq\sqrt{n}\|\Delta A\|_F$.

9. __Hoffman-Wielandt Theorem.__ If $A$ and $\tilde A$ are normal, there exists $\tau$ such that $\|\Lambda-\tilde\Lambda_\tau\|_F\leq\|\Delta A\|_F$.

10. If $A$ and $\tilde A$ are Hermitian, for any unitarily invariant norm $\|\Lambda^\uparrow-\tilde\Lambda^\uparrow\| \leq \|\Delta A\|$. In particular,

$$\begin{aligned}
\max_i|\lambda^\uparrow_i-\tilde\lambda^\uparrow_i|&\leq \|\Delta A\|_2,\\
\sqrt{\sum_i(\lambda^\uparrow_i-\tilde\lambda^\uparrow_i)^2}&\leq \|\Delta A\|_F.
\end{aligned}$$

11. __Residual bounds.__ Let $A$ be Hermitian. For some $\tilde\lambda\in\mathbb{R}$ and $\tilde x\in\mathbb{C}^n$ with $\|\tilde x\|_2=1$, define __residual__ $r=A\tilde x-\tilde\lambda\tilde x$. Then $|\tilde\lambda-\lambda|\leq \|r\|_2$ for some $\lambda\in\sigma(A)$.

12. Let, in addition,  $\tilde\lambda=\tilde x^* A\tilde x$, let $\lambda$ be closest to $\tilde\lambda$ and $x$ be its unit eigenvector, and let

$$\eta=\mathop{\mathrm{gap}}(\tilde\lambda)= \min_{\lambda\neq\mu\in\sigma(A)}|\tilde\lambda-\mu|.$$

If $\eta>0$, then

$$|\tilde\lambda-\lambda|\leq \frac{\|r\|_2^2}{\eta},\quad \sin\theta(x,\tilde x)\leq \frac{\|r\|_2}{\eta}.$$

13. Let $A$ be Hermitian, $X\in\mathbb{C}^{n\times k}$ have full column rank, and $M\in\mathcal{H}_k$ having eigenvalues $\mu_1\leq\mu_2\leq\cdots\leq\mu_k$. Set $R=AX-XM$. Then there exist $\lambda_{i_1}\leq\lambda_{i_2}\leq\cdots\leq\lambda_{i_k}\in\sigma(A)$ such that

$$
\begin{aligned}
\max_{1\leq j\leq k} |\mu_j-\lambda_{i_j}|& \leq \frac{\|R\|_2}{\sigma_{\min}(X)},\\
\sqrt{\sum_{j=1}^k (\mu_j-\lambda_{i_j})^2}&\leq \frac{\|R\|_F}{\sigma_{\min}(X)}.
\end{aligned}$$

(The indices $i_j$ need not be the same in the above formulae.)

14. If, additionally, $X^*X=I$ and $M=X^*AX$, and if all but $k$ of $A$'s eigenvalues differ from every one of $M$'s eigenvalues by at least $\eta>0$, then

$$\sqrt{\sum_{j=1}^k (\mu_j-\lambda_{i_j})^2}\leq \frac{\|R\|_F^2}{\eta\sqrt{1-\|R\|_F^2/\eta^2}}.$$

15. Let $A=\begin{bmatrix} M & E^* \\ E & H \end{bmatrix}$ and $\tilde A=\begin{bmatrix} M & 0 \\ 0 & H \end{bmatrix}$ be Hermitian, and set $\eta=\min |\mu-\nu|$ over all $\mu\in\sigma(M)$ and $\nu\in\sigma(H)$. Then

$$\max |\lambda_j^\uparrow -\tilde\lambda_j^\uparrow| \leq \frac{2\|E\|_2^2}{\eta+\sqrt{\eta^2+4\|E\|_2^2}}.$$

16. Let

$$\begin{bmatrix} X_1^*\\ X_2^* \end{bmatrix} A \begin{bmatrix} X_1 & X_2 \end{bmatrix}=
\begin{bmatrix} A_1 &  \\ & A_2 \end{bmatrix}, \quad \begin{bmatrix} X_1 & X_2 \end{bmatrix} \quad \textrm{unitary},
\quad X_1\in\mathbb{C}^{n\times k}.$$

Let $Q\in\mathbb{C}^{n\times k}$ have orthonormal columns and for a Hermitian $k\times k$ matrix $M$ set
$R=AQ-QM$. Let $\eta=\min|\mu-\nu|$ over all $\mu\in\sigma(M)$ and $\nu\in\sigma(A_2)$. If $\eta > 0$, then

$$\|\sin\Theta(X_1,Q)\|_F\leq \frac{\|R\|_F}{\eta}.$$
"""

# ╔═╡ 98fade4c-cdcb-495e-9328-96e86ca8e90e
md"""
## Examples

### Nondiagonalizable matrix
"""

# ╔═╡ 8f064645-cf1f-4226-add4-53bc72665bf6
A₃=[-3 7 -1; 6 8 -2; 72 -28 19]

# ╔═╡ 568366c7-3eee-4948-a096-7bf87537412d
# (Right) eigenvectors
X₃=eigen(A₃)

# ╔═╡ fe46c02b-74bd-4303-8314-3b6e96b9e51a
cond(X₃.vectors)

# ╔═╡ c48aa813-8012-409f-9b08-ce9805da9cb3
# Left eigenvectors
Y₃=eigen(Matrix(A₃'))

# ╔═╡ f0bcadb5-882c-4e46-9079-b36097650435
cond(Y₃.vectors)

# ╔═╡ 21b412f3-f460-4f60-8b31-85c839566d62
begin
	# Try k=2,3
	k=1
	Y₃.vectors[:,k]'*A₃-Y₃.values[k]*Y₃.vectors[:,k]'
end

# ╔═╡ 6654b8d5-0c31-4a44-a16a-3f71e714ce31
begin
	ΔA₃=rand(3,3)/20
	B₃=A₃+ΔA₃
end

# ╔═╡ 3fbcbecd-6e69-49d0-94bf-06a4f6704858
norm(ΔA₃)

# ╔═╡ 210137a6-1b7b-4741-84d2-9d72e18cd999
Z₃=eigen(B₃)

# ╔═╡ e946019f-5e75-4195-8991-da371300e9a2
begin
	# Fact 2
	l=1
	Z₃.values[l]-X₃.values[l], Y₃.vectors[:,l]'*ΔA₃*X₃.vectors[:,l] /(Y₃.vectors[:,l]⋅X₃.vectors[:,l])
end

# ╔═╡ 4723a4f9-af98-4f2d-8628-0ffc58488e62
md"""
### Jordan form
"""

# ╔═╡ f3d63eb6-0174-458c-aeef-46ed6d62db7f
begin
	n₄=6
	c=0.5
	J=Bidiagonal(c*ones(n₄),ones(n₄-1),'U')
end

# ╔═╡ b9770791-4fa6-4d3f-9688-ac79b7ae7f22
# Accurately defined eigenvalues
λ₄=eigvals(J)

# ╔═╡ 8112bf44-3146-4b9f-953a-54de6f4c8c5b
# Only one eigenvector
eigvecs(J)

# ╔═╡ 6b26afc1-e427-496e-b9a5-4618551aeede
begin
	x₄=eigvecs(J)[:,1]
	y₄=eigvecs(J')[:,1]
end

# ╔═╡ b2a07659-cb08-4c53-99ac-32f817485062
y₄'*J-0.5*y₄'

# ╔═╡ 989b7c60-c673-4b77-ada8-909db2904fa2
# Just one perturbed element in the lower left corner
ΔJ=√eps()*[zeros(n₄-1);1]*Matrix(I,1,n₄)

# ╔═╡ 1da2bfee-6792-4aed-8798-91cc058cd540
μ₄=eigvals(J+ΔJ)

# ╔═╡ a6769e47-59f7-430c-9c20-a121a98f19bf
# Fact 2
maximum(abs,λ₄-μ₄)

# ╔═╡ fab42c40-7d01-4787-a491-3d37cae1768d
y₄'*ΔJ*x₄/(y₄⋅x₄)

# ╔═╡ 2f783e6e-5b29-477f-8b45-00461013c66f
md"""
However, since $J+\Delta J$ is diagonalizable, we can apply Bauer-Fike theorem to it:
"""

# ╔═╡ 9385e7e0-17d6-4840-98a2-fa69d53ebea3
Y₄=eigvecs(J+ΔJ)

# ╔═╡ 18e6478b-bb8e-482c-bd02-0c85bca6fc7f
cond(Y₄)

# ╔═╡ ae524f6a-3e64-49ea-92e3-671559f85a5f
opnorm(inv(Y₄)*ΔJ*Y₄), cond(Y₄)*opnorm(ΔJ)

# ╔═╡ 9c0c9fcf-9581-479f-8ec5-d885dd3e6665
md"""
### Normal matrix
"""

# ╔═╡ 300e0ddd-3ab0-4803-9528-665e754e54ca
begin
	n₅=5
	C=Circulant(rand(-5:5,n₅))
end

# ╔═╡ fdecade7-4d5c-460f-a566-d52868c75ec0
eigvals(Matrix(C))

# ╔═╡ 01250371-f8ab-4a0a-b283-ffbcffdc2929
ΔC=randn(n₅,n₅)*0.0001

# ╔═╡ 0787f637-44bf-4c5c-b8d2-92df842a91b5
opnorm(ΔC), eigvals(C+ΔC)

# ╔═╡ 5645f7ea-50e5-43f7-afc3-1c46fbdd7c5b
md"""
### Hermitian matrix
"""

# ╔═╡ 51cd2fc5-e295-4b76-9fd1-fdd2b92ce05a
begin
	# Random matrix with strong column scaling
	m₆=10
	n₆=6
	D₆=Diagonal(exp.(10*randn(n₆)))
	A₆=cor(randn(m₆,n₆)*D₆)
end

# ╔═╡ bc5a0ee0-931f-11eb-204c-b9e03584c787
D₆

# ╔═╡ 1eb12479-0662-410b-9569-c4e8939a291f
ΔA₆=cor(rand(m₆,n₆)*D₆)*1e-5

# ╔═╡ 821f6b8b-0e12-404b-b1dc-57d809ebf79b
begin
	λ₆,U₆=eigen(A₆)
	μ₆=eigvals(A₆+ΔA₆)
	[λ₆ μ₆]
end

# ╔═╡ ac3d3677-1160-4a79-a25f-8faa7fd89940
norm(ΔA₆)

# ╔═╡ f4051f63-a0d3-4b43-8c61-2bbd08338a72
# ?round

# ╔═╡ 612870ba-75df-4729-a16e-432f70707fb4
begin
	# Residual bounds - how close is μ, y to λ[2],X[:,2]
	k₆=3
	ζ₆=round(λ₆[k₆],sigdigits=2)
	y₆=round.(U₆[:,k₆],sigdigits=2)
	normalize!(y₆)
end

# ╔═╡ 3ba55b4b-5445-48a3-b327-b7d3af28397a
ζ₆

# ╔═╡ 6f3da676-c953-4d23-8b9a-bb62b8b8ac94
# Fact 9
r₆=A₆*y₆-ζ₆*y₆

# ╔═╡ 886ce4f8-73a0-4b57-8d5e-082bcfc1f9aa
minimum(abs,ζ₆.-λ₆), norm(r₆)

# ╔═╡ 82a4adaa-3663-4200-9f15-d09fd4cd0259
begin
	# Fact 10 - μ is Rayleigh quotient
	ξ₆=y₆⋅(A₆*y₆)
	ρ₆=A₆*y₆-ξ₆*y₆
end

# ╔═╡ b1025690-4244-4933-9693-9bbd3ad8aedb
η₆=min(abs(ξ₆-λ₆[k₆-1]),abs(ξ₆-λ₆[k₆+1]))

# ╔═╡ 5e69cab3-712c-4754-b6b1-579bc3c8e114
ξ₆-λ₆[k₆], norm(ρ₆)^2/η₆

# ╔═╡ 3792c478-fbaa-4683-849a-ca0d0cb751f3
begin
	# Eigenvector bound
	# cos(θ)
	cosθ=dot(y₆,U₆[:,k₆])
	# sin(θ)
	sinθ=sqrt(1-cosθ^2)
	sinθ,norm(ρ₆)/η₆
end

# ╔═╡ 0dd080f1-575d-43e0-9285-20279f390d2f
begin
	# Residual bounds - Fact 13
	Q₆=round.(U₆[:,1:3],sigdigits=2)
	# Orthogonalize
	F=qr(Q₆)
	X=Matrix(F.Q)
	# Make sure M₆ is hermitian
	M₆=Hermitian(X'*A₆*X)
	μₘ=eigvals(M₆)
	R₆=A₆*X-X*M₆
end

# ╔═╡ ea07eea3-bab7-4add-90cc-43acd7c2ec89
λ₆

# ╔═╡ 53f2860a-db33-49cd-b95e-3086f76c2592
μₘ

# ╔═╡ 8cdffcc2-c8c1-410e-95db-6327931782a1
M₆

# ╔═╡ 0e871e72-47ee-4492-9f73-fcfd2ef1aada
# The entries of μ are not ordered - which algorithm was called?
issymmetric(M₆)

# ╔═╡ 108f73bc-b1cf-47ed-80f0-281cb5ebd048
ηₘ=λ₆[4]-λ₆[3]

# ╔═╡ dbc219f9-9fcf-4aac-b7d5-feb8fcafe3e9
norm(λ₆[1:3]-μₘ), norm(R₆)^2/ηₘ

# ╔═╡ 7cae5501-8a51-4874-8b23-5825f79624e3
A₆

# ╔═╡ 7cb49041-f326-4b2c-90f8-99c4d58178eb
begin
	# Fact 15
	M₇=A₆[1:3,1:3]
	H=A₆[4:6,4:6]
	E=A₆[4:6,1:3]
	# Block-diagonal matrix
	B₇=cat(M₇,H,dims=(1,2))
end

# ╔═╡ 75c00141-e691-4f9c-8b7a-ed262b88a62d
begin
	η₇=minimum(abs,eigvals(M₇)-eigvals(H))
	μ₇=eigvals(B₇)
	[λ₆ μ₇]
end

# ╔═╡ 433a0652-912f-4f59-8494-4d0381e22115
2*norm(E)^2/(η₇+√(η₇^2+4*norm(E)^2))

# ╔═╡ ad6ca597-65b2-4403-a6aa-1394399c3ed8
begin
	# Eigenspace bounds - Fact 16
	B₆=A₆+ΔA₆
	μ,V=eigen(B₆)
end

# ╔═╡ 4727c0ba-124e-41a4-861a-57eb5a6d849d
begin
	# sin(Θ(U[:,1:3],V[:,1:3]))
	X₆=U₆[:,1:3]
	Q=V[:,1:3]
	cosθ₆=svdvals(√(Q'*Q)*Q'*X₆*√(X₆'*X₆))
	sinθ₆=sqrt.(1 .-cosθ₆.^2)
end

# ╔═╡ dd00cbef-e1f6-43bc-95b0-1b6694c23583
# Bound
M₈=Q'*A₆*Q

# ╔═╡ 99f44a10-b3e1-442c-a042-b2f0afad1797
R₈=A₆*Q-Q*M₈

# ╔═╡ b6f8f48e-2e80-4999-95a9-701d491411bd
eigvals(M₈), λ₆

# ╔═╡ 7a6554bb-0518-4068-9784-678f846d3b3f
begin
	η₈=abs(eigvals(M₈)[3]-λ₆[4])
	norm(sinθ₆), norm(R₈)/η₈
end

# ╔═╡ Cell order:
# ╟─4ea41550-9313-11eb-1cdf-df6c5527f80f
# ╟─7d3ed16e-50ea-4ede-b968-e32baa56f76f
# ╟─da0d7094-fa3e-4f1e-a6f1-9d1b9316f65a
# ╟─f8ba4230-9313-11eb-170c-ab58597826f8
# ╠═4c360009-f01e-45bb-b458-29b15f8cac8f
# ╠═7bd6f898-b86c-4670-96e8-9016666ee14f
# ╠═8e1ced3a-d23d-4f5f-8c65-dc13ab392504
# ╠═d9ef7b4b-4954-464a-96ac-e3f49800e0b7
# ╟─90fcd9f5-d98c-4db6-bf30-21e377209681
# ╟─8261b5ce-9315-11eb-2378-af63c3eb4c02
# ╠═bb6eee61-c3ee-4b6a-ab44-a638fba776b8
# ╠═060b8d2f-e6f3-4625-9cf6-d4cc504145fa
# ╠═f55a3c60-9315-11eb-027a-dd008031f6f9
# ╠═c4a1b888-9497-4eb4-aa1c-23671123f426
# ╠═b46cbb6f-36f4-4a66-b970-8fe9b0506b6a
# ╠═254b2df4-f57e-4284-a735-6fd5d4546c52
# ╠═09fd0923-1987-4a2e-be28-78451aa2c11a
# ╠═dcba6892-b90e-41ad-8d58-53e72d40e882
# ╠═c765c02e-9316-11eb-146c-ff5d3ba81323
# ╠═d948741f-6748-44d7-ae8c-cf074f7549bd
# ╠═ec0746a1-41bd-4d61-9deb-8e227106e592
# ╠═d92799f7-1813-450b-8353-9d03078f69f0
# ╠═0d58d4c4-16c6-4a5b-bdbd-f3f599d0d648
# ╠═b1a90c59-8c3e-4530-83ef-6b3e6c6fa687
# ╠═31e73f82-93e7-448c-940b-5c9f5e573e5c
# ╟─909033d9-bb27-4f0e-b07d-574f7fcee700
# ╟─7ef162c0-5c80-4e96-bcc3-125b302b1d5e
# ╟─3c8f71a5-0c2f-43a0-bf71-db9df707f84e
# ╟─98fade4c-cdcb-495e-9328-96e86ca8e90e
# ╠═8f064645-cf1f-4226-add4-53bc72665bf6
# ╠═568366c7-3eee-4948-a096-7bf87537412d
# ╠═fe46c02b-74bd-4303-8314-3b6e96b9e51a
# ╠═c48aa813-8012-409f-9b08-ce9805da9cb3
# ╠═f0bcadb5-882c-4e46-9079-b36097650435
# ╠═21b412f3-f460-4f60-8b31-85c839566d62
# ╠═6654b8d5-0c31-4a44-a16a-3f71e714ce31
# ╠═3fbcbecd-6e69-49d0-94bf-06a4f6704858
# ╠═210137a6-1b7b-4741-84d2-9d72e18cd999
# ╟─e946019f-5e75-4195-8991-da371300e9a2
# ╟─4723a4f9-af98-4f2d-8628-0ffc58488e62
# ╠═f3d63eb6-0174-458c-aeef-46ed6d62db7f
# ╠═b9770791-4fa6-4d3f-9688-ac79b7ae7f22
# ╠═8112bf44-3146-4b9f-953a-54de6f4c8c5b
# ╠═6b26afc1-e427-496e-b9a5-4618551aeede
# ╠═b2a07659-cb08-4c53-99ac-32f817485062
# ╠═989b7c60-c673-4b77-ada8-909db2904fa2
# ╠═1da2bfee-6792-4aed-8798-91cc058cd540
# ╠═a6769e47-59f7-430c-9c20-a121a98f19bf
# ╠═fab42c40-7d01-4787-a491-3d37cae1768d
# ╟─2f783e6e-5b29-477f-8b45-00461013c66f
# ╠═9385e7e0-17d6-4840-98a2-fa69d53ebea3
# ╠═18e6478b-bb8e-482c-bd02-0c85bca6fc7f
# ╠═ae524f6a-3e64-49ea-92e3-671559f85a5f
# ╟─9c0c9fcf-9581-479f-8ec5-d885dd3e6665
# ╠═300e0ddd-3ab0-4803-9528-665e754e54ca
# ╠═fdecade7-4d5c-460f-a566-d52868c75ec0
# ╠═01250371-f8ab-4a0a-b283-ffbcffdc2929
# ╠═0787f637-44bf-4c5c-b8d2-92df842a91b5
# ╟─5645f7ea-50e5-43f7-afc3-1c46fbdd7c5b
# ╠═51cd2fc5-e295-4b76-9fd1-fdd2b92ce05a
# ╠═bc5a0ee0-931f-11eb-204c-b9e03584c787
# ╠═1eb12479-0662-410b-9569-c4e8939a291f
# ╠═821f6b8b-0e12-404b-b1dc-57d809ebf79b
# ╠═ac3d3677-1160-4a79-a25f-8faa7fd89940
# ╠═f4051f63-a0d3-4b43-8c61-2bbd08338a72
# ╠═612870ba-75df-4729-a16e-432f70707fb4
# ╠═3ba55b4b-5445-48a3-b327-b7d3af28397a
# ╠═6f3da676-c953-4d23-8b9a-bb62b8b8ac94
# ╠═886ce4f8-73a0-4b57-8d5e-082bcfc1f9aa
# ╠═82a4adaa-3663-4200-9f15-d09fd4cd0259
# ╠═b1025690-4244-4933-9693-9bbd3ad8aedb
# ╠═5e69cab3-712c-4754-b6b1-579bc3c8e114
# ╠═3792c478-fbaa-4683-849a-ca0d0cb751f3
# ╠═0dd080f1-575d-43e0-9285-20279f390d2f
# ╠═ea07eea3-bab7-4add-90cc-43acd7c2ec89
# ╠═53f2860a-db33-49cd-b95e-3086f76c2592
# ╠═8cdffcc2-c8c1-410e-95db-6327931782a1
# ╠═0e871e72-47ee-4492-9f73-fcfd2ef1aada
# ╠═108f73bc-b1cf-47ed-80f0-281cb5ebd048
# ╠═dbc219f9-9fcf-4aac-b7d5-feb8fcafe3e9
# ╠═7cae5501-8a51-4874-8b23-5825f79624e3
# ╠═7cb49041-f326-4b2c-90f8-99c4d58178eb
# ╠═75c00141-e691-4f9c-8b7a-ed262b88a62d
# ╠═433a0652-912f-4f59-8494-4d0381e22115
# ╠═ad6ca597-65b2-4403-a6aa-1394399c3ed8
# ╠═4727c0ba-124e-41a4-861a-57eb5a6d849d
# ╠═dd00cbef-e1f6-43bc-95b0-1b6694c23583
# ╠═99f44a10-b3e1-442c-a042-b2f0afad1797
# ╠═b6f8f48e-2e80-4999-95a9-701d491411bd
# ╠═7a6554bb-0518-4068-9784-678f846d3b3f

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L3b Eigenvalue Decomposition - Perturbation Theory.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
using FastSplat
using Compat
using Compat.Test

@testset "vararg" begin
    splat_single(args...) = string(args...)
    fastsplat_single(args...) = @fastsplat string(args...)
    @test splat_single(1,2,3) == fastsplat_single(1,2,3)

    splat_trailing(args...) = string(0,args...)
    fastsplat_trailing(args...) = @fastsplat string(0,args...)
    @test splat_trailing(1,2,3) == fastsplat_trailing(1,2,3)

    splat_leading(args...) = string(args...,0)
    fastsplat_leading(args...) = @fastsplat string(args...,0)
    @test splat_leading(1,2,3) == fastsplat_leading(1,2,3)

    splat_middle(args...) = string(0,args...,0)
    fastsplat_middle(args...) = @fastsplat string(0,args...,0)
    @test splat_middle(1,2,3) == fastsplat_middle(1,2,3)
end

@testset "array" begin
    splat_single(args) = string(args...)
    fastsplat_single(args) = @fastsplat string(args...)
    @test splat_single([1,2,3]) == fastsplat_single([1,2,3])

    splat_trailing(args) = string(0,args...)
    fastsplat_trailing(args) = @fastsplat string(0,args...)
    @test splat_trailing([1,2,3]) == fastsplat_trailing([1,2,3])

    splat_leading(args) = string(args...,0)
    fastsplat_leading(args) = @fastsplat string(args...,0)
    @test splat_leading([1,2,3]) == fastsplat_leading([1,2,3])

    splat_middle(args) = string(0,args...,0)
    fastsplat_middle(args) = @fastsplat string(0,args...,0)
    @test splat_middle([1,2,3]) == fastsplat_middle([1,2,3])
end


text:
"""
Hybrid learner (core)
Convex combination of MB and MF Q-values
Time decaying (or increasing) MB weights
"""
mutable struct HybridConvex
    nr_states::Int
    nr_actions::Int
    w_MB::Float64
    w_MF::Float64
    t_hybrid::Float64
    decay_offset::Float64
    decay_slope::Float64
    QSA::Array{Float64, 2}
    do_debug_print::Bool
end
function HybridConvex(; nr_states = 10, nr_actions = 4,
    decay_offset = 0.1, decay_slope = 0.1,
    QSA = zeros(Float64, nr_states, nr_actions),
    do_debug_print = false)
    t_hybrid = 0.
    w_MB = get_weight(decay_offset, decay_slope, t_hybrid)
    w_MF = 1. - w_MB

    HybridConvex(nr_states, nr_actions, w_MB, w_MF, t_hybrid, decay_offset, decay_slope,
                QSA, do_debug_print)
end
export HybridConvex
function gethybridweights!(learner::HybridConvex)
    learner.w_MB = get_weight(learner.decay_offset, learner.decay_slope, learner.t_hybrid)
    learner.w_MF = 1. - learner.w_MB
end
function updatetimecounter!(learner::HybridConvex, reward)
    if reward>0 #goal state, reward, NO NEXT ACTION (in our paradigm)
        learner.t_hybrid += 1
    end
end
function get_weight(w0, T, t)
    # compute a weight that can decrease to 0 or grow to 1 (unlike original Glascher paper)
    epsilon = 0.01
    if T>epsilon # have a small margin, avoiding division by
        tau = 1.0 / T
        return w0 * exp(-t/tau)
    elseif T >= -epsilon
        return w0
    else
        tau = -1.0 / T
        return 1.0 - (1.0 -w0) * exp(-t/tau)
    end
end


text:
#==============================================================================#
# AWSCredentials.jl
#
# Load AWS Credentials from:
# - EC2 Instance Profile,
# - Environment variables, or
# - ~/.aws/credentials file.
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

export AWSCredentials,
       localhost_is_lambda,
       localhost_is_ec2,
       localhost_maybe_ec2,
       aws_user_arn,
       aws_account_number

"""
When you interact with AWS, you specify your [AWS Security Credentials](http://docs.aws.amazon.com/general/latest/gr/aws-security-credentials.html) to verify who you are and whether you have permission to access the resources that you are requesting. AWS uses the security credentials to authenticate and authorize your requests.

The fields `access_key_id` and `secret_key` hold the access keys used to authenticate API requests (see [Creating, Modifying, and Viewing Access Keys](http://docs.aws.amazon.com/IAM/latest/UserGuide/id_credentials_access-keys.html#Using_CreateAccessKey)).

[Temporary Security Credentials](http://docs.aws.amazon.com/IAM/latest/UserGuide/id_credentials_temp.html) require the extra session `token` field.

The `user_arn` and `account_number` fields are used to cache the result of the [`aws_user_arn`](@ref) and [`aws_account_number`](@ref) functions.

The `AWSCredentials()` constructor tries to load local Credentials from
environment variables, `~/.aws/credentials`, `~/.aws/config` or EC2 instance credentials.
"""
mutable struct AWSCredentials
    access_key_id::String
    secret_key::String
    token::String
    user_arn::String
    account_number::String

    function AWSCredentials(access_key_id, secret_key,
                            token="", user_arn="", account_number="")
        new(access_key_id, secret_key, token, user_arn, account_number)
    end
end

function Base.show(io::IO,c::AWSCredentials)
    println(io, string(c.user_arn,
                       c.user_arn == "" ? "" : " ",
                       "(",
                       c.account_number,
                       c.account_number == "" ? "" : ", ",
                       c.access_key_id,
                       c.secret_key == "" ? "" : ", $(c.secret_key[1:3])...",
                       c.token == "" ? "" : ", $(c.token[1:3])..."),
                       ")")
end

function Base.copyto!(dest::AWSCredentials, src::AWSCredentials)
    for f in fieldnames(typeof(dest))
        setfield!(dest, f, getfield(src, f))
    end
end
import Base: copy!
Base.@deprecate copy!(dest::AWSCredentials, src::AWSCredentials) copyto!(dest, src)

function AWSCredentials()

    if haskey(ENV, "AWS_ACCESS_KEY_ID")

        creds = env_instance_credentials()

    elseif isfile(dot_aws_credentials_file()) || isfile(dot_aws_config_file())

        creds = dot_aws_credentials()

    elseif haskey(ENV, "AWS_CONTAINER_CREDENTIALS_RELATIVE_URI")

        creds = ecs_instance_credentials()

    elseif localhost_maybe_ec2()

        creds = ec2_instance_credentials()

    else
        error("Can't find AWS credentials!")
    end

    if debug_level > 0
        display(creds)
        println()
    end

    return creds
end

"""
Is Julia running in an AWS Lambda sandbox?
"""
localhost_is_lambda() = haskey(ENV, "LAMBDA_TASK_ROOT")

"""
Is Julia running on an EC2 virtual machine?
"""
function localhost_is_ec2()

    if localhost_is_lambda()
        return false
    end

    if isfile("/sys/hypervisor/uuid") &&
       String(read("/sys/hypervisor/uuid",3)) == "ec2"
        return true
    end

    if isfile("/sys/devices/virtual/dmi/id/product_uuid")
        try
            # product_uuid is not world readable!
            # https://patchwork.kernel.org/patch/6461521/
            # https://github.com/JuliaCloud/AWSCore.jl/issues/24
            if String(read("/sys/devices/virtual/dmi/id/product_uuid")) == "EC2"
                return true
            end
        catch
        end
    end

    return false
end

localhost_maybe_ec2() = localhost_is_ec2() ||
                        isfile("/sys/devices/virtual/dmi/id/product_uuid")

"""
    aws_user_arn(::AWSConfig)

Unique
[Amazon Resource Name]
(http://docs.aws.amazon.com/IAM/latest/UserGuide/id_users.html)
for configrued user.

e.g. `"arn:aws:iam::account-ID-without-hyphens:user/Bob"`
"""
function aws_user_arn(aws::AWSConfig)

    creds = aws[:creds]

    if creds.user_arn == ""

        r = Services.sts(aws, "GetCallerIdentity", [])
        creds.user_arn = r["Arn"]
        creds.account_number = r["Account"]
    end

    return creds.user_arn
end

"""
    aws_account_number(::AWSConfig)

12-digit [AWS Account Number](http://docs.aws.amazon.com/general/latest/gr/acct-identifiers.html).
"""
function aws_account_number(aws::AWSConfig)
    creds = aws[:creds]
    if creds.account_number == ""
        aws_user_arn(aws)
    end
    return creds.account_number
end

"""
    ec2_metadata(key)

Fetch [EC2 meta-data]
(http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ec2-instance-metadata.html)
for `key`.
"""
function ec2_metadata(key)

    @assert localhost_maybe_ec2()

    String(http_get("http://169.254.169.254/latest/meta-data/$key").body)
end

"""
Load [Instance Profile Credentials]
(http://docs.aws.amazon.com/AWSEC2/latest/UserGuide/iam-roles-for-amazon-ec2.html#instance-metadata-security-credentials)
for EC2 virtual machine.
"""
function ec2_instance_credentials()

    @assert localhost_maybe_ec2()

    info  = ec2_metadata("iam/info")
    info  = LazyJSON.value(info)

    name  = ec2_metadata("iam/security-credentials/")
    creds = ec2_metadata("iam/security-credentials/$name")
    new_creds = LazyJSON.value(creds)

    if debug_level > 0
        print("Loading AWSCredentials from EC2 metadata... ")
    end

    AWSCredentials(new_creds["AccessKeyId"],
                   new_creds["SecretAccessKey"],
                   new_creds["Token"],
                   info["InstanceProfileArn"])
end

"""
Load [ECS Task Credentials]
(http://docs.aws.amazon.com/AmazonECS/latest/developerguide/task-iam-roles.html)
"""
function ecs_instance_credentials()

    @assert haskey(ENV, "AWS_CONTAINER_CREDENTIALS_RELATIVE_URI")

    uri = ENV["AWS_CONTAINER_CREDENTIALS_RELATIVE_URI"]

    new_creds = String(http_get("http://169.254.170.2$uri").body)
    new_creds = LazyJSON.value(new_creds)

    if debug_level > 0
        print("Loading AWSCredentials from ECS metadata... ")
    end

    AWSCredentials(new_creds["AccessKeyId"],
                   new_creds["SecretAccessKey"],
                   new_creds["Token"],
                   new_creds["RoleArn"])
end

"""
Load Credentials from [environment variables]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-environment.html)
`AWS_ACCESS_KEY_ID`, `AWS_SECRET_ACCESS_KEY` etc.
(e.g. in Lambda sandbox).
"""
function env_instance_credentials()

    if debug_level > 0
        print("Loading AWSCredentials from ENV[\"AWS_ACCESS_KEY_ID\"]... ")
    end

    AWSCredentials(ENV["AWS_ACCESS_KEY_ID"],
                   ENV["AWS_SECRET_ACCESS_KEY"],
                   get(ENV, "AWS_SESSION_TOKEN", ""),
                   get(ENV, "AWS_USER_ARN", ""))
end

using IniFile

dot_aws_credentials_file() = get(ENV, "AWS_SHARED_CREDENTIALS_FILE",
                                 joinpath(homedir(), ".aws", "credentials"))

dot_aws_config_file() = get(ENV, "AWS_CONFIG_FILE",
                                 joinpath(homedir(), ".aws", "config"))

function aws_get_role_details(profile::AbstractString, ini::Inifile)
    if debug_level > 0
        println("Loading \"$profile\" Profile from " *
                dot_aws_config_file() * "... ")
    end

    role_arn = get(ini, profile, "role_arn")
    source_profile = get(ini, profile, "source_profile")

    profile = "profile $profile"
    role_arn = get(ini, profile, "role_arn", role_arn)
    source_profile = get(ini, profile, "source_profile", source_profile)

    (source_profile, role_arn)
end

function aws_get_credential_details(profile::AbstractString, ini::Inifile, config::Bool)
    if debug_level > 0
        filename = config ? dot_aws_config_file() : dot_aws_credentials_file()
        println("Loading \"$profile\" AWSCredentials from " * filename
                * "... ")
    end

    key_id = get(ini, profile, "aws_access_key_id")
    key = get(ini, profile, "aws_secret_access_key")
    token = get(ini, profile, "aws_session_token", "")

    if config
        profile = "profile $profile"
        key_id = get(ini, profile, "aws_access_key_id", key_id)
        key = get(ini, profile, "aws_secret_access_key", key)
        token = get(ini, profile, "aws_session_token", token)
    end

    (key, key_id, token)
end

function aws_get_region(profile::AbstractString, ini::Inifile)
    region = get(ENV, "AWS_DEFAULT_REGION", "us-east-1")

    region = get(ini, profile, "region", region)
    region = get(ini, "profile $profile", "region", region)
end

function aws_get_role(role::AbstractString, ini::Inifile)
    source_profile, role_arn = aws_get_role_details(role, ini)

    if source_profile == :notfound
        error("Can't find AWS credentials!")
    end

    if debug_level > 0
        println("Assuming \"$source_profile\"... ")
    end
    credentials = dot_aws_credentials(source_profile)

    config = AWSConfig(:creds=>credentials, :region=>aws_get_region(source_profile, ini))

    role = Services.sts(
        config,
        "AssumeRole",
        RoleArn=role_arn,
        RoleSessionName=replace(role, r"[^\w+=,.@-]" => s"-"),
    )
    role_creds = role["Credentials"]

    credentials = AWSCredentials(role_creds["AccessKeyId"],
        role_creds["SecretAccessKey"],
        role_creds["SessionToken"]
    )
end

"""
Load Credentials from [AWS CLI ~/.aws/credentials file] or [AWS CLI ~/.aws/config file]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html).
"""
function dot_aws_credentials(profile = nothing)
    @assert isfile(dot_aws_credentials_file()) || isfile(dot_aws_config_file())

    if profile == nothing
        profile = get(ENV, "AWS_DEFAULT_PROFILE",
                get(ENV, "AWS_PROFILE", "default"))
    end

    # According to the docs the order of precedence is:
    # 1. credentials in the credential file
    # 2. credentials in the config file
    # 3. roles in the config file
    credential_file = dot_aws_credentials_file()
    ini = nothing
    if isfile(credential_file)
        ini = read(Inifile(), credential_file)
        key, key_id, token = aws_get_credential_details(profile, ini, false)
        if key != :notfound
            return AWSCredentials(key_id, key, token)
        end
    end

    config_file = dot_aws_config_file()
    if isfile(config_file)
        ini = read(Inifile(), config_file)
        key, key_id, token = aws_get_credential_details(profile, ini, true)
        if key != :notfound
            AWSCredentials(key_id, key, token)
        else
            aws_get_role(profile, ini)
        end
    end
end

#==============================================================================#
# End of file.
#==============================================================================#