text:
module FEMMShellT3DSGMTModule

using LinearAlgebra: norm, Transpose, mul!, diag, eigen, I, dot
using Statistics: mean
using FinEtools
import FinEtools.FESetModule: gradN!, nodesperelem, manifdim
using FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FEMMBaseModule: associategeometry!
using FinEtoolsDeforLinear.MatDeforLinearElasticModule: tangentmoduli!, update!, thermalstrain!
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, locjac!, add_nnt_ut_only!, add_btsigma!
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.TransformerModule: QTEQTransformer

const __nn = 3 # number of nodes
const __ndof = 6 # number of degrees of freedom per node

function local_frame!(F0, J0)
    # This is the tangent to the coordinate curve 1
    a = @view J0[:, 1]
    L0 = norm(a);
    F0[:,1] = a/L0;
    # This is the tangent to the coordinate curve 2
    b = J0[:, 2]
    #     F0(:,3)=skewmat(F0(:,1))*b;
    b = @view J0[:, 2]
    # Now compute the normal
    F0[:, 3] .= (-F0[3,1]*b[2]+F0[2,1]*b[3],
                  F0[3,1]*b[1]-F0[1,1]*b[3],
                 -F0[2,1]*b[1]+F0[1,1]*b[2]);
    F0[:, 3] /= norm(@view F0[:, 3]);
    #     F0(:,2)=skewmat(F0(:,3))*F0(:,1);
    F0[:, 2] .= (-F0[3,3]*F0[2,1]+F0[2,3]*F0[3,1],
                  F0[3,3]*F0[1,1]-F0[1,3]*F0[3,1],
                 -F0[2,3]*F0[1,1]+F0[1,3]*F0[2,1]);
    return F0
end

"""
    FEMMShellT3DSGMTMT{S<:AbstractFESet, F<:Function} <: AbstractFEMM

Class for Discrete Shear Gap shell finite element modeling machine. With
averaging of the transverse strain-displacement matrix to provide isotropic
transverse shear response. Also, the formulation is developed to correctly
handle the coupling of twisting moments and transverse shear (such as in the
twisted beam or the Raasch hook problems) by incorporating "nodal" normals. The
MacNeal transformation between element and nodal ordinates is used (hence MT).

Programming developed consistently with the paper
[1] Cui et al, Analysis of plates and shells using an edge-based smoothed finite
element method, Comput Mech (2010) 45:141–156
DOI 10.1007/s00466-009-0429-9

In this reference, the sign next to Ae in equation (44) is wrong.
[2] A superconvergent alpha finite element method (S a FEM) for static and
free vibration analysis of shell structures
Chai et al. (2017).

The treatment of the transformation between the element and nodal coordinates
is carried out as described in the paper:

Finite Elements in Analysis and Design 30 (1998) 235—242
The treatment of shell normals in ﬁnite element analysis
Richard H. MacNeal, Charles T. Wilson, Robert L. Harder, Claus C. Hoﬀ
The MacNeal-Schwendler Corporation, 815 Colorado Blvd., Los Angeles, CA 90041, USA

The following features are incorporated to deal with nodal normals:
- Nodal normals are averages of the normals of elements that meet at a node.
- A crease in the surface is take into account. In that case the normal are not
  averaged across the crease. At the nodes along the crease every element uses
  the normal to its surface instead of the nodal normal.
"""
mutable struct FEMMShellT3DSGMT{S<:AbstractFESet, F<:Function, M} <: AbstractFEMM
    integdomain::IntegDomain{S, F} # integration domain data
    material::M # material object.
    drilling_stiffness_scale::FFlt
    _associatedgeometry::Bool
    _normals::FFltMat
    _normal_valid::Vector{Bool}
    # The attributes below are buffers used in various operations.
    _loc::FFltMat
    _J::FFltMat
    _J0::FFltMat
    _ecoords::FFltMat
    _edisp::FFltMat
    _ecoords_e::FFltMat
    _dofnums::FIntMat
    _e_g::FFltMat
    _n_e::Vector{FFltMat} # transformation nodal-element matrices
    _nvalid::Vector{Bool}
    _Te::FFltMat
    _elmat::FFltMat
    _loc_e::FFltMat
    _J_e::FFltMat
    _gradN_e::FFltMat
    _Bm::FFltMat
    _Bb::FFltMat
    _Bs::FFltMat
    _DpsBmb::FFltMat
    _DtBs::FFltMat
end

function FEMMShellT3DSGMT(integdomain::IntegDomain{S, F}, material::M) where {S<:AbstractFESet, F<:Function, M}
    _nnmax = 0
    for j in 1:count(integdomain.fes)
        for k in eachindex(integdomain.fes.conn[j])
            _nnmax = max(_nnmax, integdomain.fes.conn[j][k])
        end
    end
    _normals = fill(0.0, _nnmax, 3)
    _normal_valid = fill(true, _nnmax)
    # Alocate_e the buffers2
    _loc = fill(0.0, 1, 3)
    _J = fill(0.0, 3, 2)
    _J0 = fill(0.0, 3, 2)
    _ecoords = fill(0.0, __nn, 3)
    _edisp = fill(0.0, __nn, 6);
    _ecoords_e = fill(0.0, __nn, 2)
    _dofnums = zeros(FInt, 1, __nn*__ndof);
    _e_g = fill(0.0, 3, 3);
    _n_e = [fill(0.0, 3, 3), fill(0.0, 3, 3), fill(0.0, 3, 3)];
    _nvalid = fill(false, 3)
    _Te = fill(0.0, __nn*__ndof, __nn*__ndof)
    _elmat = fill(0.0, __nn*__ndof, __nn*__ndof);
    _loc_e = fill(0.0, 1, 2)
    _J_e = fill(0.0, 2, 2)
    _gradN_e = fill(0.0, __nn, 2)
    _Bm = fill(0.0, 3, __nn*__ndof)
    _Bb = fill(0.0, 3, __nn*__ndof)
    _Bs = fill(0.0, 2, __nn*__ndof)
    _DpsBmb = similar(_Bm)
    _DtBs = similar(_Bs)

    return FEMMShellT3DSGMT(integdomain, material, 1.0,
        false,
        _normals, _normal_valid,
        _loc, _J, _J0,
        _ecoords, _edisp, _ecoords_e,
        _dofnums,
        _e_g, _n_e, _nvalid, _Te,
        _elmat,
        _loc_e, _J_e, _gradN_e,
        _Bm, _Bb, _Bs, _DpsBmb, _DtBs)
end

function make(integdomain, material)
    return FEMMShellT3DSGMT(integdomain, material)
end

function _compute_J0!(J0, ecoords)
    # Compute the Jacobian: the Jacobian is constant within the triangle
    x, y, z = ecoords[2, :].-ecoords[1, :]
    J0[:, 1] .= (x, y, z)
    x, y, z = ecoords[3, :].-ecoords[1, :]
    J0[:, 2] .= (x, y, z)
end

function _shell_material_stiffness(material)
    # Compute the two material stiffness matrices: the plane stress matrix, and
    # the transverse shear matrix
    D = fill(0.0, 6, 6)
    t::FFlt, dt::FFlt, loc::FFltMat, label::FInt = 0.0, 0.0, [0.0 0.0 0.0], 0
    tangentmoduli!(material,  D,  t, dt, loc, label)
    Dps = fill(0.0, 3, 3)
    Dps[1:2, 1:2] = D[1:2, 1:2] -  (reshape(D[1:2,3], 2, 1) * reshape(D[3,1:2], 1, 2))/D[3, 3]
    ix=[1, 2, 4];
    for i = 1:3
        Dps[3,i] = Dps[i,3] = D[4, ix[i]];
    end
    Dt = fill(0.0, 2, 2)
    ix=[5, 6];
    for i = 1:2
        Dt[i,i] = D[ix[i], ix[i]];
    end
    return Dps, Dt
end

function _nodal_triads_e!(n_e, nvalid, e_g, normals, normal_valid, c)
    # Components of nodal cartesian ordinate systems such that the third
    # direction is the direction of the nodal normal, and the angle to rotate
    # the element normal into the nodal normal is as short as possible; these
    # components are given on the local element coordinate system basis
    # vectors. If the angle of rotation is excessively small, the nodal matrix
    # is the identity.
    #
    # These triads are the nodal-to-element transformation
    # matrices.

    # TO DO Get rid of the temporaries
    r = fill(0.0, 3)
    nk_e = fill(0.0, 3)
    nk = fill(0.0, 3)
    f3_e = [0.0, 0.0, 1.0]
    for k in 1:length(n_e)
        nk .= vec(view(normals, c[k], :))
        nvalid[k] = normal_valid[c[k]]
        if nvalid[k]
            # If the nodal normal is valid, pull it back into the element frame.
            nk_e .= e_g'*nk
        else
            # Otherwise the element normal replaces the nodal normal.
            nk_e .= 0.0; nk_e[3] = 1.0
        end
        cross3!(r, f3_e, nk_e)
        if norm(r) > 1.0e-6
            rotmat3!(n_e[k], r)
        else
            n_e[k] .= I(3)
        end
    end
    return n_e, nvalid
end

function _transfmat_g_to_n!(Te, n_e, e_g)
    # Global-to-nodal transformation matrix.

    # The 3x3 blocks consist of the nodal triad expressed on the global basis.
    # The nodal basis vectors in `n_e[i]` are expressed on the element basis,
    # and are then rotated with `e_g` into the global coordinate system.

    # Output
    # - `Te` = transformation matrix, input in the global basis, output in the
    #   nodal basis
    Te .= 0.0
    for i in 1:__nn
        Teblock = transpose(n_e[i]) * transpose(e_g) # TO DO remove temporary
        offset = (i-1)*__ndof
        r = offset+1:offset+3
        @. Te[r, r] = Teblock
        r = offset+4:offset+6
        @. Te[r, r] = Teblock
    end
    return Te
end

function _transfmat_n_to_e!(Te, n_e, gradN_e)
# Nodal-to-element transformation matrix.
# lTn = matrix with the nodal triad vectors in columns, components
# on the element basis
# TO DO avoid a temporary
    Te .= 0.0
# Translation degrees of freedom
    for i in 1:__nn
        roffset = (i-1)*__ndof
        r = roffset+1:roffset+3
        Te[r, r] .= n_e[i] # this needs to be inverse
    end
# Rotation degrees of freedom. The drilling rotation of the mid surface
# produced by the 1/2*(v,x - u,y) effect is linked to the out of plane
# rotations.
    for i in 1:__nn
        roffst = (i-1)*__ndof
        n_ei = n_e[i] # TO DO avoid temporary
        invn_e = inv(n_ei[1:2, 1:2])
        Te[roffst.+(4:5), roffst.+(4:5)] .= invn_e
        Te[roffst+6, roffst+6] = 1/n_ei[3, 3]
        a1, a2 = invn_e * vec(n_ei[1:2, 3]) # TO DO avoid the temporary
        for j in 1:__nn
            coffst = (j-1)*__ndof
            Te[roffst+4, coffst+1] = (-a1 * 1/2 * gradN_e[j, 2])
            Te[roffst+4, coffst+2] = (+a1 * 1/2 * gradN_e[j, 1])
            Te[roffst+5, coffst+1] = (-a2 * 1/2 * gradN_e[j, 2])
            Te[roffst+5, coffst+2] = (+a2 * 1/2 * gradN_e[j, 1])
            # Te[roffst+4, coffst+1] += (-a1 * 1/2) * (n_ei[1, 1] * gradN_e[j, 2])
            # Te[roffst+4, coffst+2] += (-a1 * 1/2) * (n_ei[1, 2] * gradN_e[j, 2])
            # Te[roffst+4, coffst+3] += (-a1 * 1/2) * (n_ei[1, 3] * gradN_e[j, 2])
            # Te[roffst+4, coffst+1] += (+a1 * 1/2) * (n_ei[2, 1] * gradN_e[j, 1])
            # Te[roffst+4, coffst+2] += (+a1 * 1/2) * (n_ei[2, 2] * gradN_e[j, 1])
            # Te[roffst+4, coffst+3] += (+a1 * 1/2) * (n_ei[2, 3] * gradN_e[j, 1])
            # Te[roffst+5, coffst+1] += (-a2 * 1/2) * (n_ei[1, 1] * gradN_e[j, 2])
            # Te[roffst+5, coffst+2] += (-a2 * 1/2) * (n_ei[1, 2] * gradN_e[j, 2])
            # Te[roffst+5, coffst+3] += (-a2 * 1/2) * (n_ei[1, 3] * gradN_e[j, 2])
            # Te[roffst+5, coffst+1] += (+a2 * 1/2) * (n_ei[2, 1] * gradN_e[j, 1])
            # Te[roffst+5, coffst+2] += (+a2 * 1/2) * (n_ei[2, 2] * gradN_e[j, 1])
            # Te[roffst+5, coffst+3] += (+a2 * 1/2) * (n_ei[2, 3] * gradN_e[j, 1])
        end
    end
    return Te
end

"""
    _Bsmat!(Bs, gradN, N)

Compute the linear transverse shear strain-displacement matrix.
"""
function _Bsmat!(Bs, ecoords_e)
    Bs .= 0.0

    # Orientation 1
    s, p, q = 3, 1, 2
    a, b = ecoords_e[p, :] .- ecoords_e[s, :]
    c, d = ecoords_e[q, :] .- ecoords_e[s, :]
    Ae = (a*d - b*c)/2

    # The first node in the triangle
    m = (1/2/Ae) * (1/3) # multiplier
    # Node s
    co = (s - 1) * 6 # column offset
    Bs[1, co+3] += m*(b-d);                             Bs[1, co+5] += m*(Ae)
    Bs[2, co+3] += m*(c-a); Bs[2, co+4] += m*(-Ae);
    # The other two nodes
    # Node p
    co = (p - 1) * 6 # column offset
    Bs[1, co+3] += m*(d);   Bs[1, co+4] += m*(-b*d/2);  Bs[1, co+5] += m*(a*d/2)
    Bs[2, co+3] += m*(-c);  Bs[2, co+4] += m*(b*c/2);   Bs[2, co+5] += m*(-a*c/2)
    # Node q
    co = (q - 1) * 6 # column offset
    Bs[1, co+3] += m*(-b);  Bs[1, co+4] += m*(b*d/2);   Bs[1, co+5] += m*(-b*c/2)
    Bs[2, co+3] += m*(a);   Bs[2, co+4] += m*(-a*d/2);  Bs[2, co+5] += m*(a*c/2)

    # Orientation 2
    s, p, q = 1, 2, 3
    a, b = ecoords_e[p, :] .- ecoords_e[s, :]
    c, d = ecoords_e[q, :] .- ecoords_e[s, :]

    # The first node in the triangle
    m = (1/2/Ae) * (1/3) # multiplier
    # Node s
    co = (s - 1) * 6 # column offset
    Bs[1, co+3] += m*(b-d);                             Bs[1, co+5] += m*(Ae)
    Bs[2, co+3] += m*(c-a); Bs[2, co+4] += m*(-Ae);
    # The other two nodes
    # Node p
    co = (p - 1) * 6 # column offset
    Bs[1, co+3] += m*(d);   Bs[1, co+4] += m*(-b*d/2);  Bs[1, co+5] += m*(a*d/2)
    Bs[2, co+3] += m*(-c);  Bs[2, co+4] += m*(b*c/2);   Bs[2, co+5] += m*(-a*c/2)
    # Node q
    co = (q - 1) * 6 # column offset
    Bs[1, co+3] += m*(-b);  Bs[1, co+4] += m*(b*d/2);   Bs[1, co+5] += m*(-b*c/2)
    Bs[2, co+3] += m*(a);   Bs[2, co+4] += m*(-a*d/2);  Bs[2, co+5] += m*(a*c/2)

    # Orientation 3
    s, p, q = 2, 3, 1
    a, b = ecoords_e[p, :] .- ecoords_e[s, :]
    c, d = ecoords_e[q, :] .- ecoords_e[s, :]

    # The first node in the triangle
    m = (1/2/Ae) * (1/3) # multiplier
    # Node s
    co = (s - 1) * 6 # column offset
    Bs[1, co+3] += m*(b-d);                             Bs[1, co+5] += m*(Ae)
    Bs[2, co+3] += m*(c-a); Bs[2, co+4] += m*(-Ae);
    # The other two nodes
    # Node p
    co = (p - 1) * 6 # column offset
    Bs[1, co+3] += m*(d);   Bs[1, co+4] += m*(-b*d/2);  Bs[1, co+5] += m*(a*d/2)
    Bs[2, co+3] += m*(-c);  Bs[2, co+4] += m*(b*c/2);   Bs[2, co+5] += m*(-a*c/2)
    # Node q
    co = (q - 1) * 6 # column offset
    Bs[1, co+3] += m*(-b);  Bs[1, co+4] += m*(b*d/2);   Bs[1, co+5] += m*(-b*c/2)
    Bs[2, co+3] += m*(a);   Bs[2, co+4] += m*(-a*d/2);  Bs[2, co+5] += m*(a*c/2)

    return Bs
end

"""
    _Bmmat!(Bm, gradN)

Compute the linear membrane strain-displacement matrix.
"""
function _Bmmat!(Bm, gradN)
    for i in 1:__nn
        Bm[1,6*(i-1)+1] = gradN[i,1];
        Bm[2,6*(i-1)+2] = gradN[i,2];
        Bm[3,6*(i-1)+1] = gradN[i,2];
        Bm[3,6*(i-1)+2] = gradN[i,1];
    end
end

"""
    _Bbmat!(Bb, gradN)

Compute the linear, displacement independent, curvature-displacement/rotation matrix for a shell quadrilateral element with nfens=3 nodes. Displacements and rotations are in a local coordinate system.
"""
function _Bbmat!(Bb, gradN)
    for i in 1:__nn
        Bb[1,6*(i-1)+5] = gradN[i,1];
        Bb[2,6*(i-1)+4] = -gradN[i,2];
        Bb[3,6*(i-1)+4] = -gradN[i,1];
        Bb[3,6*(i-1)+5] = gradN[i,2];
    end
end

function associategeometry!(self::FEMMShellT3DSGMT,  geom::NodalField{FFlt})
    J0 = self._J0
    e_g = self._e_g
    normals, normal_valid = self._normals, self._normal_valid
    # Compute the normals at the nodes
    for el in 1:count(self.integdomain.fes)
        i, j, k = self.integdomain.fes.conn[el]
        J0[:, 1] = geom.values[j, :] - geom.values[i, :]
        J0[:, 2] = geom.values[k, :] - geom.values[i, :]
        local_frame!(e_g, J0)
        normals[i, :] .+= e_g[:, 3]
        normals[j, :] .+= e_g[:, 3]
        normals[k, :] .+= e_g[:, 3]
    end
    # Normalize to unit length
    for j in 1:size(normals, 1)
        nn = norm(normals[j, :])
        if nn > 0.0
            normals[j, :] ./= nn
        end
    end
    # Now perform a second pass. If the nodal normal differs by a substantial
    # amount from the element normals, zero out the nodal normal to indicate
    # that the nodal normal should not be used at that vertex. TO DO: Should we
    # make this a configurable parameter?
    ntolerance = 1 - sqrt(1 - sin(30/180*pi)^2)
    for el in 1:count(self.integdomain.fes)
        i, j, k = self.integdomain.fes.conn[el]
        J0[:, 1] = geom.values[j, :] - geom.values[i, :]
        J0[:, 2] = geom.values[k, :] - geom.values[i, :]
        local_frame!(e_g, J0)
        en = view(e_g, :, 3)
        nd = dot(normals[i, :], en)
        if nd < 1 - ntolerance
            normal_valid[i] = false
        end
        nd = dot(normals[j, :], en)
        if nd < 1 - ntolerance
            normal_valid[j] = false
        end
        nd = dot(normals[k, :], en)
        if nd < 1 - ntolerance
            normal_valid[k] = false
        end
    end
    nz = 0
    for j in 1:size(normals, 1)
        if !normal_valid[j]
            nz += 1
        end
    end
    # @show size(normals, 1) , nz
    self._associatedgeometry = true
    return self
end

"""
    stiffness(self::FEMMShellT3DSGMT, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{T}) where {ASS<:AbstractSysmatAssembler, T<:Number}

Compute the material stiffness matrix.
"""
function stiffness(self::FEMMShellT3DSGMT, assembler::ASS, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {ASS<:AbstractSysmatAssembler, T<:Number, TI<:Number}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    normals, normal_valid = self._normals, self._normal_valid
    loc, J, J0 = self._loc, self._J, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e = self._ecoords_e
    e_g, n_e, nvalid, Te = self._e_g, self._n_e, self._nvalid, self._Te
    elmat = self._elmat
    transformwith = QTEQTransformer(elmat)
    loc_e, J_e, gradN_e = self._loc_e, self._J_e, self._gradN_e
    Bm, Bb, Bs, DpsBmb, DtBs = self._Bm, self._Bb, self._Bs, self._DpsBmb, self._DtBs
    Br = fill(0.0, 1, size(elmat, 1))
    Dps, Dt = _shell_material_stiffness(self.material)
    scf=5/6;  # shear correction factor
    Dt .*= scf
    drilling_stiffness_scale = self.drilling_stiffness_scale
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes), dchi.nfreedofs, dchi.nfreedofs);
    for i in 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        _compute_J0!(J0, ecoords)
        local_frame!(e_g, J0)
        fill!(elmat,  0.0); # Initialize element matrix
        for j in 1:npts
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            t = self.integdomain.otherdimension(loc, fes.conn[i], Ns[j])
            mul!(ecoords_e, ecoords, view(e_g, :, 1:2))
            locjac!(loc_e, J_e, ecoords_e, Ns[j], gradNparams[j])
            gradN!(fes, gradN_e, gradNparams[j], J_e);
            _Bmmat!(Bm, gradN_e)
            _Bbmat!(Bb, gradN_e)
            _Bsmat!(Bs, ecoords_e)
            add_btdb_ut_only!(elmat, Bm, t*Jac*w[j], Dps, DpsBmb)
            add_btdb_ut_only!(elmat, Bb, (t^3)/12*Jac*w[j], Dps, DpsBmb)
            # TO DO The stabilization expression has a significant effect
            # (at least for the pinched cylinder). What is the recommended
            # multiplier of he^2?
            he = sqrt(Jac)
            # add_btdb_ut_only!(elmat, Bs, (t)*Jac*w[j], Dt, DtBs)
            add_btdb_ut_only!(elmat, Bs, (t^3/(t^2+0.2*he^2))*Jac*w[j], Dt, DtBs)
        end
        # Complete the elementwise matrix by filling in the lower triangle
        complete_lt!(elmat)
        # Now treat the transformation from the element to the nodal triad
        _nodal_triads_e!(n_e, nvalid, e_g, normals, normal_valid, fes.conn[i])
        _transfmat_n_to_e!(Te, n_e, gradN_e)
         # Transform the elementwise matrix into the nodal coordinates
        transformwith(elmat, Te)
        # Bending diagonal stiffness coefficients
        kavg4 = mean((elmat[4, 4], elmat[10, 10], elmat[16, 16]))
        kavg5 = mean((elmat[5, 5], elmat[11, 11], elmat[17, 17]))
        kavg = (kavg4 + kavg5) * drilling_stiffness_scale
        # Add the artificial drilling stiffness in the nodal cordinates, but
        # only if the nodal normal is valid.
        nvalid[1] && (elmat[6,6] += kavg)
        nvalid[2] && (elmat[12,12] += kavg)
        nvalid[3] && (elmat[18,18] += kavg)
        # Transform from nodal into global coordinates
        _transfmat_g_to_n!(Te, n_e, e_g)
        transformwith(elmat, Te)
        # Assembly
        gatherdofnums!(dchi, dofnums, fes.conn[i]);
        assemble!(assembler, elmat, dofnums, dofnums);
    end # Loop over elements
    return makematrix!(assembler);
end

function stiffness(self::FEMMShellT3DSGMT, geom0::NodalField{FFlt}, u1::NodalField{T}, Rfield1::NodalField{T}, dchi::NodalField{TI}) where {T<:Number, TI<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return stiffness(self, assembler, geom0, u1, Rfield1, dchi);
end

"""
    mass(self::FEMMShellT3DSGMT,  assembler::A,  geom::NodalField{FFlt}, dchi::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}

Compute the consistent mass matrix

This is a general routine for the shell FEMM.
"""
function mass(self::FEMMShellT3DSGMT,  assembler::A,  geom0::NodalField{FFlt}, dchi::NodalField{T}) where {A<:AbstractSysmatAssembler, T<:Number}
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    normals, normal_valid = self._normals, self._normal_valid
    loc, J, J0 = self._loc, self._J, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e = self._ecoords_e
    e_g, n_e, nvalid, Te = self._e_g, self._n_e, self._nvalid, self._Te
    elmat = self._elmat
    transformwith = QTEQTransformer(elmat)
    loc_e, J_e, gradN_e = self._loc_e, self._J_e, self._gradN_e
    rho::FFlt = massdensity(self.material); # mass density
    npe = nodesperelem(fes)
    tmss = fill(0.0, npe);# basis f. matrix -- buffer
    rmss = fill(0.0, npe);# basis f. matrix -- buffer
    ndn = ndofs(dchi)
    startassembly!(assembler,  size(elmat,1),  size(elmat,2),  count(fes), dchi.nfreedofs,  dchi.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        _compute_J0!(J0, ecoords)
        local_frame!(e_g, J0)
        fill!(tmss, 0.0)
        fill!(rmss, 0.0)
        # Compute the translational and rotational masses corresponding to nodes
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            t = self.integdomain.otherdimension(loc, fes.conn[i], Ns[j])
            mul!(ecoords_e, ecoords, view(e_g, :, 1:2))
            tfactor = rho*(t*Jac*w[j]);
            rfactor = rho*(t^3/12*Jac*w[j]);
            for k in 1:npe
                tmss[k] += tfactor*(Ns[j][k])
                rmss[k] += rfactor*(Ns[j][k])
            end
        end # Loop over quadrature points
        # Now treat the transformation from the element to the nodal triad
        _nodal_triads_e!(n_e, nvalid, e_g, normals, normal_valid, fes.conn[i])
        # Fill the elementwise matrix in the nodal basis: the matrix is lumped,
        # hence diagonal, and no transformation from element to nodal is needed.
        fill!(elmat,  0.0); # Initialize element matrix
        for k in 1:npe
            # Translation degrees of freedom
            for d in 1:3
                c = (k - 1) * __ndof + d
                elmat[c, c] += tmss[k]
            end
            # Bending degrees of freedom
            for d in 4:5
                c = (k - 1) * __ndof + d
                elmat[c, c] += rmss[k]
            end
            # Drilling rotations
            if nvalid[k]
                d = 6
                c = (k - 1) * __ndof + d
                elmat[c, c] += rmss[k] / 1e6
            end
        end
        # Transform into global coordinates
        _transfmat_n_to_g!(Te, n_e, e_g)
        transformwith(elmat, Te)
        # Assemble
        gatherdofnums!(dchi,  dofnums,  fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler,  elmat,  dofnums,  dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function mass(self::FEMMShellT3DSGMT,  geom::NodalField{FFlt},  u::NodalField{T}) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return mass(self, assembler, geom, u);
end

"""
    inspectintegpoints(self::AbstractFEMMDeforLinear,
      geom::NodalField{FFlt},  u::NodalField{T},
      dT::NodalField{FFlt},
      felist::FIntVec,
      inspector::F,  idat, quantity=:Cauchy;
      context...) where {T<:Number, F<:Function}

Inspect integration point quantities.

- `geom` - reference geometry field
- `u` - displacement field
- `dT` - temperature difference field
- `felist` - indexes of the finite elements that are to be inspected:
     The fes to be included are: `fes[felist]`.
- `context`    - structure: see the update!() method of the material.
- `inspector` - function with the signature
        `idat = inspector(idat, j, conn, x, out, loc);`
   where
    `idat` - a structure or an array that the inspector may
           use to maintain some state,  for instance minimum or maximum of
           stress, `j` is the element number, `conn` is the element connectivity,
           `out` is the output of the update!() method,  `loc` is the location
           of the integration point in the *reference* configuration.
### Return
The updated inspector data is returned.
"""
function inspectintegpoints(self::FEMMShellT3DSGMT, geom::NodalField{FFlt},  u::NodalField{T}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {T<:Number, F<:Function}
    fes = self.integdomain.fes
    @assert self._associatedgeometry == true
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    normals, normal_valid = self._normals, self._normal_valid
    loc, J, J0 = self._loc, self._J, self._J0
    ecoords, edisp, dofnums = self._ecoords, self._edisp, self._dofnums
    ecoords_e = self._ecoords_e
    e_g, n_e, n, ln, nvalid, Te = self._e_g, self._n_e, self._n, self._nv_e, self._nvalid, self._Te
    elmat = self._elmat
    loc_e, J_e, gradN_e = self._loc_e, self._J_e, self._gradN_e
    Bm, Bb, Bs, DpsBmb, DtBs = self._Bm, self._Bb, self._Bs, self._DpsBmb, self._DtBs
    Br = fill(0.0, 1, size(elmat, 1))
    Dps, Dt = _shell_material_stiffness(self.material)
    scf=5/6;  # shear correction factor
    Dt .*= scf
    # Sort out  the output requirements
    outputcsys = self.mcsys; # default: report the stresses in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = val
        end
    end
    # t= 0.0
    # dt = 0.0
    # dTe = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    # ue = fill(zero(FFlt), size(elmat, 1)); # array of node displacements -- buffer
    # nne = nodesperelem(fes); # number of nodes for element
    # sdim = ndofs(geom);            # number of space dimensions
    # xe = fill(zero(FFlt), nne, sdim); # array of node coordinates -- buffer
    # qpdT = 0.0; # node temperature increment
    # qpstrain = fill(zero(FFlt), nstressstrain(self.mr), 1); # total strain -- buffer
    # qpthstrain = fill(zero(FFlt), nthermstrain(self.mr)); # thermal strain -- buffer
    # qpstress = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    # out1 = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    # out =  fill(zero(FFlt), nstressstrain(self.mr));# output -- buffer
    # Loop over  all the elements and all the quadrature points within them
    for ilist = 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom0, ecoords, fes.conn[i]);
        gathervalues_asmat!(u1, edisp, fes.conn[i]);
        ecoords1 .= ecoords .+ edisp
        _compute_J0!(J0, ecoords)
        local_frame!(e_g, J0)
        fill!(elmat,  0.0); # Initialize element matrix
        for j in 1:npts
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            Jac = Jacobiansurface(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            t = self.integdomain.otherdimension(loc, fes.conn[i], Ns[j])
            mul!(ecoords_e, ecoords, view(e_g, :, 1:2))
            locjac!(loc_e, J_e, ecoords_e, Ns[j], gradNparams[j])
            gradN!(fes, gradN_e, gradNparams[j], J_e);
            _Bmmat!(Bm, gradN_e)
            _Bbmat!(Bb, gradN_e)
            _Bsmat!(Bs, ecoords_e)
            # add_btdb_ut_only!(elmat, Bm, t*Jac*w[j], Dps, DpsBmb)
            # add_btdb_ut_only!(elmat, Bb, (t^3)/12*Jac*w[j], Dps, DpsBmb)
            # TO DO The stabilization expression has a significant effect
            # (at least for the pinched cylinder). What is the recommended
            # multiplier of he^2?
            he = sqrt(Jac)
            # add_btdb_ut_only!(elmat, Bs, (t^3/(t^2+0.2*he^2))*Jac*w[j], Dt, DtBs)
            updatecsmat!(outputcsys, loc, J, fes.label[i]);
            # Call the inspector
            #     idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
        end # Loop over quadrature points

        # gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        # gathervalues_asvec!(u, ue, fes.conn[i]);# retrieve element displacements
        # # gathervalues_asvec!(dT, dTe, fes.conn[i]);# retrieve element temp. increments
        # for j = 1:npts # Loop over quadrature points
        #     locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
        #     Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
        #     updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        #     At_mul_B!(csmatTJ,  self.mcsys.csmat,  J); # local Jacobian matrix
        #     gradN!(fes, gradN, gradNparams[j], csmatTJ);
        #     Blmat!(self.mr, B, Ns[j], gradN, loc, self.mcsys.csmat);
        #     updatecsmat!(outputcsys, loc, J, fes.label[i]);
        #     # Quadrature point quantities
        #     A_mul_B!(qpstrain, B, ue); # strain in material coordinates
        #     qpdT = dot(vec(dTe), vec(Ns[j]));# Quadrature point temperature increment
        #     thermalstrain!(self.material, qpthstrain, qpdT)
        #     # Material updates the state and returns the output
        #     out = update!(self.material, qpstress, out, vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
        #     if (quantity == :Cauchy)   # Transform stress tensor,  if that is "out"
        #         (length(out1) >= length(out)) || (out1 = zeros(length(out)))
        #         rotstressvec!(self.mr, out1, out, transpose(self.mcsys.csmat))# To global coord sys
        #         rotstressvec!(self.mr, out, out1, outputcsys.csmat)# To output coord sys
        #     end
        #     # Call the inspector
        #     idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
        # end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function inspectintegpoints(self::FEMMShellT3DSGMT, geom::NodalField{FFlt},  u::NodalField{T}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {T<:Number, F<:Function}
    dT = NodalField(fill(zero(FFlt), nnodes(geom), 1)) # zero difference in temperature
    return inspectintegpoints(self, geom, u, dT, felist, inspector, idat, quantity; context...);
end

end # module


text:
resolve_flat(tab, x::Symbol) = (x in colnames(tab)) ? x => bycol(identity) : throw(KeyError("column `$(x)` is not present in this $(typeof(tab))."))

function resolve_flat(tab, x::Pair{<:Any,<:Callable})
    column_names = _resolve(tab, selection(first(x)))
    if length(column_names) > 1
        column_names => byrow.(last(x), Tuple(column_names))
    else
        column_names => transformation(last(x))
    end
end

function resolve_flat(tab, x::Pair{<:Any,<:Union{RowwiseTrans,TabwiseTrans}})
    column_names = _resolve(tab, selection(first(x)))
    column_names => typeof(last(x))(last(x).f, Tuple(column_names))
end

function resolve_flat(tab, x::Pair{<:Any,<:ColwiseTrans})
    _resolve(tab, selection(first(x))) => last(x)
end

resolve_nested(tab, s::SelectionQuery) = SelectionResult(_resolve(tab, s.s), s.r, s.t)
resolve_nested(tab, s::AbstractMultiSelection) = _resolve(tab, s)

positive_selection(tab, symbols::Tuple{Vararg{Symbol}}, b::Bool) = b ? collect(symbols) : setdiff(colnames(tab), symbols)
positive_selection(tab, symbols::Vector{Symbol}, b::Bool) = b ? symbols : setdiff(colnames(tab), symbols)
positive_selection(tab, symbol::Symbol, b::Bool) = b ? [symbol] : setdiff(colnames(tab), [symbol])
positive_selection(tab, s::AbstractArray{Bool}, b::Bool) = b ? @inbounds(colnames(tab)[s]) : @inbounds(colnames(tab)[.-s])

function _resolve(tab, s::SymbolSelection)
    if !(s.s in colnames(tab))
        throw(KeyError(s.s))
    else
        positive_selection(tab, s.s, bool(s))
    end
end

function _resolve(tab, s::IntSelection)
    positive_selection(tab, colnames(tab)[s.s], bool(s))
end

function _resolve(tab, s::BoolSelection)
    c = colnames(tab)
    if length(c) != length(s.s)
        throw(BoundsError(c, s.s))
    else
        @inbounds(c[bool(s) ? s.s : .!s.s])
    end
end

function _resolve(tab, s::RangeSelection{<:Symbol})
    s1, s2 = s.s1, s.s2
    c = colnames(tab)

    i = findfirst(x -> x == s1, c)
    i == nothing && throw(KeyError(s1))

    j = findfirst(x -> x == s2, c)
    j == nothing && throw(KeyError(s2))

    idx = i > j ? reverse(j:s.step:i) : i:s.step:j
    positive_selection(tab, @inbounds(c[idx]), bool(s))
end

function _resolve(tab, s::RangeSelection{<:Int})
    s1, s2 = s.s1, s.s2
    c = colnames(tab)
    idx = s1 > s2 ? reverse(s2:abs(s.step):s1) : s1:s.step:s2
    positive_selection(tab, c[idx], bool(s))
end

function _resolve(tab, s::ArraySelection{<:Symbol})
    sd = setdiff(s.s, colnames(tab))
    !(isempty(sd)) && throw(ArgumentError("column$(length(sd) == 1 ? "" : "s") `$(join(sd, "`, `"))` not found."))
    positive_selection(tab, s.s, bool(s))
end

function _resolve(tab, s::ArraySelection{<:Int})
    positive_selection(tab, colnames(tab)[s.s], bool(s))
end

function _resolve(tab, s::Union{Tuple,AbstractArray{<:AbstractSelection}})
    if length(s) == 1
        _resolve(tab, selection(@inbounds(s[1])))
    else
        _resolve(tab, cols(s...))
    end
end

function _resolve(tab, s::PredicateSelection)
    [k for (k, v) in pairs(map(s.f, columntable(tab))) if bool(s) ? v : !v]
end
function _resolve(tab, s::KeyPredicateSelection)
    c = colnames(tab)
    [k for (k, v) in (c .=> map(s.f, string.(c))) if bool(s) ? v : !v]
end
function _resolve(tab, s::PairPredicateSelection)
    [k for (k, v) in pairs(columntable(tab)) if bool(s) ? s.f(string(k), v) : !s.f(string(k), v)]
end
_resolve(tab, s::AllSelection) = collect(colnames(tab))
# These needs to be resolved at later stage
_resolve(tab, s::OtherSelection) = s
_resolve(tab, s::ElseSelection) = s

ext: jl


text:
function _update_primal_constitutive!(e, d, C::Array{T,5}, nx, ny, nz) where {T}
    @turbo for k = 2:nz-1
        for j = 2:ny-1
            for i = 2:nx-1
                dy = 0.25 * (d[i, j, k, 2] + d[i+1, j, k, 2] + d[i+1, j-1, k, 2] + d[i, j-1, k, 2])
                dz = 0.25 * (d[i, j, k, 3] + d[i+1, j, k, 3] + d[i+1, j, k-1, 3] + d[i, j, k-1, 3])
                e[i, j, k, 1] = C[i, j, k, 1, 1] * d[i, j, k, 1] + C[i, j, k, 1, 2] * dy + C[i, j, k, 1, 3] * dz

                dx = 0.25 * (d[i, j, k, 1] + d[i, j+1, k, 1] + d[i-1, j, k, 1] + d[i-1, j+1, k, 1])
                dz = 0.25 * (d[i, j, k, 3] + d[i+1, j, k, 3] + d[i+1, j, k-1, 3] + d[i, j, k-1, 3])

                e[i, j, k, 2] = C[i, j, k, 2, 1] * dx + C[i, j, k, 2, 2] * d[i, j, k, 2] + C[i, j, k, 2, 3] * dz

                dx = 0.25 * (d[i, j, k, 1] + d[i, j, k+1, 1] + d[i-1, j, k+1, 1] + d[i-1, j, k, 1])
                dy = 0.25 * (d[i, j, k, 2] + d[i, j, k+1, 2] + d[i, j-1, k+1, 2] + d[i, j-1, k, 2])

                e[i, j, k, 3] = C[i, j, k, 3, 1] * dx + C[i, j, k, 3, 2] * dy + C[i, j, k, 3, 3] * d[i, j, k, 3]

            end
        end
    end

    @turbo for k = 2:nz-1
        for j = 2:ny-1
            for i = 1:1

                dy = 0.25 * (d[i, j, k, 2] + d[i+1, j, k, 2] + d[i+1, j-1, k, 2] + d[i, j-1, k, 2])
                dz = 0.25 * (d[i, j, k, 3] + d[i+1, j, k, 3] + d[i+1, j, k-1, 3] + d[i, j, k-1, 3])
                e[i, j, k, 1] = C[i, j, k, 1, 1] * d[i, j, k, 1] + C[i, j, k, 1, 2] * dy + C[i, j, k, 1, 3] * dz

            end
        end
    end
    @turbo for k = 2:nz-1
        for j = 1:1
            for i = 2:nx-1

                dx = 0.25 * (d[i, j, k, 1] + d[i, j+1, k, 1] + d[i-1, j, k, 1] + d[i-1, j+1, k, 1])
                dz = 0.25 * (d[i, j, k, 3] + d[i+1, j, k, 3] + d[i+1, j, k-1, 3] + d[i, j, k-1, 3])
                e[i, j, k, 2] = C[i, j, k, 2, 1] * dx + C[i, j, k, 2, 2] * d[i, j, k, 2] + C[i, j, k, 2, 3] * dz

            end
        end
    end
    @turbo for k = 1:1
        for j = 2:ny-1
            for i = 2:nx-1
                dx = 0.25 * (d[i, j, k, 1] + d[i, j, k+1, 1] + d[i-1, j, k+1, 1] + d[i-1, j, k, 1])
                dy = 0.25 * (d[i, j, k, 2] + d[i, j, k+1, 2] + d[i, j-1, k+1, 2] + d[i, j-1, k, 2])
                e[i, j, k, 3] = C[i, j, k, 3, 1] * dx + C[i, j, k, 3, 2] * dy + C[i, j, k, 3, 3] * d[i, j, k, 3]

            end
        end
    end
end

function _update_primal_constitutive!(e, d, C::Array{T,4}, nx, ny, nz) where {T}
    @turbo for p in 1:3
        for k = 2:nz-1
            for j = 2:ny-1
                for i = 2:nx-1
                    e[i, j, k, p] = C[i, j, k, p] * d[i, j, k, p]
                end
            end
        end
    end
    @turbo for k = 2:nz-1
        for j = 2:ny-1
            for i = 1:1
                e[i, j, k, 1] = C[i, j, k, 1] * d[i, j, k, 1]
            end
        end
    end
    @turbo for k = 2:nz-1
        for j = 1:1
            for i = 2:nx-1
                e[i, j, k, 2] = C[i, j, k, 2] * d[i, j, k, 2]
            end
        end
    end
    @turbo for k = 1:1
        for j = 2:ny-1
            for i = 2:nx-1
                e[i, j, k, 3] = C[i, j, k, 3] * d[i, j, k, 3]
            end
        end
    end
end


text:
__precompile__()
module ElectronDisplay

using Electron

struct ElectronDisplayType <: Base.Display end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("image/png")}, x)
    img = stringmime(MIME("image/png"), x)

    payload = string("")

    w = Electron.Window(payload)
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("image/png")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)

    w = Electron.Window(payload)
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("image/svg+xml")}) = true

asset(url...) = replace(normpath(joinpath(@__DIR__, "..", "assets", "vega", url...)), "\\", "/")

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v2+json"), x)

    html_page = """

    """

    w = Electron.Window(html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true

function Base.display(d::ElectronDisplayType, x)
    if mimewritable("application/vnd.vegalite.v2+json", x)
        display(d,MIME("application/vnd.vegalite.v2+json"), x)
    elseif mimewritable("image/svg+xml", x)
        display(d,"image/svg+xml", x)
    elseif mimewritable("image/png", x)
        display(d,"image/png", x)
    else
        throw(MethodError(Base.display,(d,x)))
    end
end

function __init__()
    Base.Multimedia.pushdisplay(ElectronDisplayType())
end

end # module


text:
module Agent

import ..SnakeAI

import DataStructures: CircularBuffer
import Flux
import StatsBase: sample
import Zygote: Buffer

export SnakeAgent,
    train!,
    CircularBufferMemory

include("memory.jl")
include("agents.jl")
include("functions.jl")

end # module


text:
using ControlSystemIdentification, ControlSystems, Optim, Plots, DSP, TotalLeastSquares
using Test, Random, LinearAlgebra, Statistics

using ControlSystemIdentification: time1, time2
using MonteCarloMeasurements

function ⟂(x)
    u, s, v = svd(x)
    u * v
end
function generate_system(nx, ny, nu)
    U, S = ⟂(randn(nx, nx)), diagm(0 => 0.2 .+ 0.78rand(nx))
    A    = S * U
    B    = randn(nx, nu)
    C    = randn(ny, nx)
    sys  = ss(A, B, C, 0, 1)
end

wtest = exp10.(LinRange(-3, log10(pi), 30))
freqresptest(G, model) =
    maximum(abs, log10.(abs2.(freqresp(model, wtest))) - log10.(abs2.(freqresp(G, wtest))))

freqresptest(G, model, tol) = freqresptest(G, model) < tol

@testset "ControlSystemIdentification.jl" begin

    @testset "basis functions" begin
        @info "Testing basis functions"
        include("test_basis_functions.jl")
    end

    @testset "FRD arma fit" begin
        @info "Testing FRD arma fit"
        include("test_frq_tf.jl")
    end

    @testset "Frequency weights" begin
        @info "Testing Frequency weights"
        include("test_frequency_weights.jl")
    end

    @testset "Utils" begin
        @info "Testing Utils"
        v = zeros(2)
        M = zeros(2, 4)
        vv = fill(v, 4)
        @test time1(v) == v
        @test time1(M) == M'
        @test time1(vv) == M'
        @test time2(v) == v'
        @test time2(M) == M
        @test time2(vv) == M

        v = zeros(2)
        M = zeros(1, 2)
        vv = [[0.0], [0.0]]
        for x in (v, M, vv), t in (M, vv)
            # @show typeof(t), typeof(x)
            @test oftype(t, t) == t
            @test typeof(oftype(t, x)) == typeof(t)
        end

        y = randn(1, 3)
        @test ControlSystemIdentification.modelfit(y, y) == [100]

    end

    @testset "iddata" begin
        @info "Testing iddata"
        @testset "vectors" begin
            T = 100
            y = randn(1, T)
            @show d = iddata(y)
            @test d isa ControlSystemIdentification.OutputData
            @test length(d) == T
            @test output(d) == y
            @test !hasinput(d)
            @test ControlSystemIdentification.time2(y) == y
            @test sampletime(d) == 1
            @test d[1:10] isa typeof(d)
            @test length(d[1:10]) == 10
            @test length(timevec(d)) == length(d)

            @test_nowarn plot(d)

            u = randn(1, T)
            @show d = iddata(y, u)
            @test d isa ControlSystemIdentification.InputOutputData
            @test length(d) == T
            @test output(d) == y
            @test hasinput(d)
            @test input(d) == u
            @test ControlSystemIdentification.time2(y) == y
            @test sampletime(d) == 1
            @test d[1, 1] == d
            @test d[1:10] isa typeof(d)

            @test length([d d]) == 2length(d)

            @test oftype(Matrix, output(d)) == y

            yr, A = getARXregressor(d, 2, 2)
            @test size(A, 2) == 4

            @test_nowarn plot(d)
        end

        @testset "matrices" begin
            T = 10
            ny, nu = 2, 3
            y = randn(ny, T)
            @show d = iddata(y)
            @test d isa ControlSystemIdentification.OutputData
            @test length(d) == T
            @test output(d) == y
            @test !hasinput(d)
            @test ControlSystemIdentification.time1(y) == y'
            @test sampletime(d) == 1

            @test_nowarn plot(d)

            u = randn(nu, T)
            @show d = iddata(y, u)
            @test d isa ControlSystemIdentification.InputOutputData
            @test length(d) == T
            @test output(d) == y
            @test hasinput(d)
            @test input(d) == u
            @test sampletime(d) == 1

            @test_nowarn plot(d)

            u = randn(T, nu)
            @show d = iddata(y, u, 2)
            @test d isa ControlSystemIdentification.InputOutputData
            @test length(d) == T
            @test output(d) == y
            @test hasinput(d)
            @test input(d) == u'
            @test ControlSystemIdentification.time1(input(d)) == u
            @test sampletime(d) == 2

            @test_nowarn plot(d)

        end

        @testset "vectors of vectors" begin
            T = 100
            ny, nu = 2, 3
            y = [randn(ny) for _ = 1:T]
            @show d = iddata(y)
            @test d isa ControlSystemIdentification.OutputData
            @test length(d) == T
            @test output(d) == y
            @test !hasinput(d)
            @test ControlSystemIdentification.time1(y) == reduce(hcat, y)'
            @test sampletime(d) == 1

            @test_nowarn plot(d)

            u = [randn(nu) for _ = 1:T]
            @show d = iddata(y, u)
            @test d isa ControlSystemIdentification.InputOutputData
            @test length(d) == T
            @test output(d) == y
            @test hasinput(d)
            @test input(d) == u
            @test sampletime(d) == 1

            u = randn(T, nu)
            @show d = iddata(y, u, 2)
            @test d isa ControlSystemIdentification.InputOutputData
            @test length(d) == T
            @test output(d) == y
            @test hasinput(d)
            @test input(d) == u'
            @test ControlSystemIdentification.time1(y) == reduce(hcat, y)'
            @test ControlSystemIdentification.time1(input(d)) == u
            @test sampletime(d) == 2

            @test_nowarn plot(d)

        end

    end

    @testset "n4sid" begin
        @info "Testing n4sid"

        N = 200
        r = 2
        m = 2
        l = 2
        for r = 1:5, m = 1:2, l = 1:2
            Random.seed!(0)
            @show r, m, l
            A = Matrix{Float64}(I(r))
            A[1, 1] = 1.01
            G = ss(A, randn(r, m), randn(l, r), 0 * randn(l, m), 1)
            u = randn(N, m)
            x0 = randn(r)
            y, t, x = lsim(G, u, 1:N, x0 = x0)
            @test sum(!isfinite, y) == 0
            yn = y + 0.1randn(size(y))
            d = iddata(yn, u, 1)
            res = n4sid(d, r, γ = 0.99)
            @test maximum(abs, pole(res.sys)) <= 1.00001 * 0.99

            ys = simulate(res, d, res.x[:, 1])
            # @show mean(abs2,y-ys') / mean(abs2,y)
            ys = simulate(res, d, res.x[:, 1], stochastic = true)
            # @show mean(abs2,y-ys') / mean(abs2,y)

            yp = predict(res, d, res.x[:, 1])
            # @show mean(abs2,y-yp') / mean(abs2,y)
            @test mean(abs2, y - yp') / mean(abs2, y) < 0.01

            res = n4sid(d, r, i = r + 1)
            freqresptest(G, res.sys) < 0.2 * m * l
            w = exp10.(LinRange(-5, log10(pi), 600))
            bodeplot(G, w)
            bodeplot!(res.sys, w)

            G = ss(0.2randn(r, r), randn(r, m), randn(l, r), 0 * randn(l, m), 1)
            u = randn(N, m)

            y, t, x = lsim(G, u, 1:N, x0 = randn(r))
            @assert sum(!isfinite, y) == 0
            ϵ = 0.01
            yn = y + ϵ * randn(size(y))
            d = iddata(yn, u)

            res = n4sid(d, r)
            @test res.sys.nx == r
            @test sqrt.(diag(res.R)) ≈ ϵ * ones(l) rtol = 0.5
            @test norm(res.S) < ϵ
            @test freqresptest(G, res.sys) < 0.2 * m * l

            if m == l # == 1
                iy, it, ix = impulse(G, 20)
                H = okid(d, r, 20)
                @test norm(iy - permutedims(H, (3, 1, 2))) < 0.05
                # plot([iy vec(H)])

                sys = era(d, r)
                @test sys.nx == r
                @test freqresptest(G, sys) < 0.2 * m * l
            end

            res = ControlSystemIdentification.n4sid(d)
            @test res.sys.nx <= r # test that auto rank selection don't choose too high rank when noise is low
            kf = KalmanFilter(res)
            @test kf isa KalmanFilter
            @test kf.A == res.A
            @test kf.B == res.B
            @test kf.C == res.C

        end

        m = l = r = 1
        for a in LinRange(-0.99, 0.99, 10), b in LinRange(-5, 5, 10)
            G = tf(b, [1, -a], 1)
            u = randn(N, m)
            y, t, x = lsim(G, u, 1:N, x0 = randn(r))
            yn = y + 0.01randn(size(y))
            d = iddata(yn, u)
            res = n4sid(d, r)
            @test res.sys.A[1] ≈ a atol = 0.01
            @test numvec(tf(res.sys))[1][2] ≈ b atol = 0.01
            @test abs(numvec(tf(res.sys))[1][1]) < 1e-2
            @test freqresptest(G, res.sys) < 0.01
        end

    end

    @testset "pem" begin
        Random.seed!(1)
        T = 1000
        nx = 3
        nu = 1
        ny = 1
        x0 = randn(nx)
        sim(sys, u, x0 = x0) = lsim(sys, u', 1:T, x0 = x0)[1]'
        sys = generate_system(nx, nu, ny)
        sysn = generate_system(nx, nu, ny)

        σu = 0
        σy = 0

        u  = randn(nu, T)
        un = u + sim(sysn, σu * randn(size(u)), 0 * x0)
        y  = sim(sys, un, x0)
        yn = y + sim(sysn, σy * randn(size(u)), 0 * x0)
        d  = iddata(yn, un, 1)

        # using BenchmarkTools
        # @btime begin
        # Random.seed!(0)
        sysh, x0h, opt = pem(d, nx = nx, focus = :prediction, iterations=5000)
        # bodeplot([sys,ss(sysh)], exp10.(range(-3, stop=log10(pi), length=150)), legend=false, ylims=(0.01,100))
        # end
        # 462ms 121 29
        # 296ms
        # 283ms
        # 173ms
        @test sysh.C * x0h ≈ sys.C * x0 atol = 0.1
        @test freqresptest(sys, StateSpace(sysh)) < 1e-7
        yh = sim(convert(StateSpace, sysh), u, x0h)
        @test Optim.minimum(opt) < 1 # Should reach 0

        # Test with some noise
        # Only measurement noise
        σu = 0.0
        σy = 0.1
        u = randn(nu, T)
        un = u + sim(sysn, σu * randn(size(u)), 0 * x0)
        y = sim(sys, un, x0)
        yn = y + sim(sysn, σy * randn(size(u)), 0 * x0)
        d = iddata(yn, un, 1)
        sysh, x0h, opt = pem(d, nx = nx, focus = :prediction)
        @test sysh.C * x0h ≈ sys.C * x0 atol = 0.1
        @test Optim.minimum(opt) < 2σy^2 * T # A factor of 2 margin

        # Only input noise
        σu = 0.1
        σy = 0.0
        u = randn(nu, T)
        un = u + sim(sysn, σu * randn(size(u)), 0 * x0)
        y = sim(sys, un, x0)
        yn = y + sim(sysn, σy * randn(size(u)), 0 * x0)
        d = iddata(yn, un, 1)
        @time sysh, x0h, opt = pem(d, nx = nx, focus = :prediction)
        @test sysh.C * x0h ≈ sys.C * x0 atol = 0.1
        @test Optim.minimum(opt) < 1 # Should depend on system gramian, but too lazy to figure out

        # Both noises
        σu = 0.2
        σy = 0.2

        u = randn(nu, T)
        un = u + sim(sysn, σu * randn(size(u)), 0 * x0)
        y = sim(sys, un, x0)
        yn = y + sim(sysn, σy * randn(size(u)), 0 * x0)
        d = iddata(yn, un, 1)
        sysh, x0h, opt = pem(d, nx = 3nx, focus = :prediction, iterations = 400)
        @test sysh.C * x0h ≈ sys.C * x0 atol = 0.1
        @test Optim.minimum(opt) < 2σy^2  # A factor of 2 margin

        # Simulation error minimization
        σu = 0.01
        σy = 0.01

        u = randn(nu, T)
        un = u + sim(sysn, σu * randn(size(u)), 0 * x0)
        y = sim(sys, un, x0)
        yn = y + sim(sysn, σy * randn(size(u)), 0 * x0)
        d = iddata(yn, un, 1)
        @time sysh, x0h, opt = pem(d, nx = nx, focus = :simulation)
        @test sysh.C * x0h ≈ sys.C * x0 atol = 0.3
        @test Optim.minimum(opt) < 0.01

        # L1 error minimization
        σu = 0.01
        σy = 0.01

        u = randn(nu, T)
        un = u + sim(sysn, σu * randn(size(u)), 0 * x0)
        y = sim(sys, un, x0)
        yn = y + sim(sysn, σy * randn(size(u)), 0 * x0)
        d = iddata(yn, un, 1)
        sysh, x0h, opt = pem(
            d,
            nx = nx,
            focus = :prediction,
            metric = e -> sum(abs, e),
            regularizer = p -> (0.1 / T) * norm(p),
        )
        # 409ms
        @test sysh.C * x0h ≈ sys.C * x0 atol = 0.1
        @test Optim.minimum(opt) < 0.01

        yh = ControlSystemIdentification.predict(sysh, yn, u, x0h)
        @test mean(abs2, y - yh) < 0.01

        yh = ControlSystemIdentification.simulate(sysh, u, x0h)
        @test mean(abs2, y - yh) < 0.01

        yh = ControlSystemIdentification.predict(sysh, iddata(yn, u), x0h)
        @test mean(abs2, y - yh) < 0.01
    end

    """
    Compare a tf to a tf based on particles .. maybe replaye by correct dispatched isapprox
    """
    function compareTFs(tf, tfStochastic, atol = 10e-5)
        aok = bok = false
        try
            bok = isapprox(numvec(tf), mean.(numvec(tfStochastic)), atol = atol)
            aok = isapprox(denvec(tf), mean.(denvec(tfStochastic)), atol = atol)
        catch
            return false
        end
            return aok && bok
    end

    @testset "arx" begin
        @info "Testing arx"

        # Test params2poly2
        y = randn(10)
        u = randn(10)
        for na in 1:3, nb in 1:3, inputdelay in 0:2
            Y, A = getARXregressor(y, u, na, nb, inputdelay = inputdelay)
            w = \(A, Y)
            a, b = ControlSystemIdentification.params2poly2(w, na, nb)
            @test length(a) == length(b[1])
            @test length(w) == na + nb
        end

        unsafe_comparisons(true)

        N = 20
        t = 1:N
        u = randn(N)
        G = tf(0.8, [1, -0.9], 1)
        y = lsim(G, u, t)[1][:]

        pars = ControlSystemIdentification.params(G)
        @test pars == ([0.9, 0.8], [0.9], [[0.8]], [1])
        @test ControlSystemIdentification.params2poly(pars[1], 1, 1) == ([1, -0.9], [[0.8]])

        na, nb = 1, 1
        yr, A = getARXregressor(y, u, na, nb)
        @test length(yr) == N - na
        @test size(A) == (N - na, na + nb)

        @test yr == y[na+1:end]
        @test A[:, 1] == y[1:end-na]
        @test A[:, 2] == u[1:end-1]

        na = 1
        d = iddata(y, u, 1)
        Gh = arx(d, na, nb)
        @test Gh ≈ G # Should recover the original transfer function exactly
        @test freqresptest(G, Gh, 0.0001)
        ω = exp10.(range(-2, stop = 1, length = 200))
        # bodeplot(G,ω)
        # bodeconfidence!(Gh, Σ, ω=ω, color=:blue, lab="Est")

        # Test MISO estimation
        u2 = randn(N)
        G2 = [G tf(0.5, [1, -0.9], 1)]
        y2 = lsim(G2, [u u2], t)[1][:]

        nb = [1, 1]
        d = iddata(y2, [u u2], 1)
        Gh2 = arx(d, na, nb)
        @test Gh2 ≈ G2
        Gh2s = arx(d, na, nb, stochastic = true)
        @test compareTFs(G2, Gh2s)

        # SISO
        ## Test na < nb
        na, nb = 1, 2
        G1 =  tf([1, -2], [1, -0.5, 0], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb)
        @test G1 ≈ Gest

        Gests = arx(d, na, nb, stochastic = true)
        @test compareTFs(G1, Gests)

        ## Test na > nb
        na, nb = 2, 1
        G1  = tf([0.8, 0], [1,0.5,0.25], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb)
        @test G1 ≈ Gest

        # MISO nb1 != nb2
        na, nb = 1, [2, 3]
        G1 =  tf([1, -2], [1, -0.5, 0], 1)
        G2 =  tf([1, -2, 3], [1, -0.5, 0, 0], 1)
        G = [G1 G2]
        u1 = randn(N)
        u2 = randn(N)
        u = [u1 u2]
        y = lsim(G, u, t)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb)
        @test Gest ≈ G

        Gests = arx(d, na, nb, stochastic = true)
        @test compareTFs(G, Gests)

        # Test na = 0
        na, nb = 0, 1
        G1 =  tf([1], [1, 0], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb)
        @test Gest ≈ G1

        Gests = arx(d, na, nb, stochastic = true)
        @test compareTFs(G1, Gests)

        # Test inputdelay
        ## SISO
        na, nb, inputdelay = 1, 1, 2
        G1 =  tf([0,1], [1, -0.5,0], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb, inputdelay = inputdelay)
        @test Gest ≈ G1
        Gests = arx(d, na, nb, inputdelay = inputdelay, stochastic = true)
        @test compareTFs(G1, Gests)

        ## Test na > nb
        na, nb, inputdelay = 2, 1, 2
        G1   = tf(0.8, [1,0.5,0.25], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb, inputdelay = inputdelay)
        @test G1 ≈ Gest

        ## MISO
        na, nb, inputdelay = 1, [1, 1], [2, 3]
        G1 =  tf([0,1], [1, -0.5,0], 1)
        G2 =  tf([0,0,1], [1, -0.5,0,0], 1)
        G = [G1 G2]
        u1 = randn(N)
        u2 = randn(N)
        u = [u1 u2]
        y = lsim(G, u, 1:N)[1][:]
        d = iddata(y, u, 1)
        Gest = arx(d, na, nb, inputdelay = inputdelay)
        @test Gest ≈ G

        Gests = arx(d, na, nb, inputdelay = inputdelay, stochastic =  true)
        @test compareTFs(G, Gests)

        # direct input
        G1 =  tf([0.3, 1], [1, -0.5], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        na, nb, inputdelay = 1,2,0
        Gest = arx(d, na, nb, inputdelay = inputdelay)
        @test Gest ≈ G1

        Gests = arx(d, na, nb, inputdelay = inputdelay, stochastic = true)
        @test compareTFs(G1, Gests)

        ## with inputdelay
        G1 =  tf([0.3,0, 1], [1, -0.5, 0], 1)
        u = randn(N)
        y = lsim(G1, u, t)[1][:]
        d = iddata(y, u, 1)
        na, nb, inputdelay = 1,3,0
        Gest = arx(d, na, nb, inputdelay = inputdelay)
        @test Gest ≈ G1

        Gests = arx(d, na, nb, inputdelay = inputdelay, stochastic = true)
        @test compareTFs(G1, Gests)

    end

    @testset "ar" begin
        N = 10000
        t = 1:N
        y = zeros(N)
        y[1] = -0.2
        u = copy(y)
        for i = 2:N
            y[i] = 0.9y[i-1]
        end
        # G = tf(1, [1, -0.9], 1)
        G = tf([1, 0], [1, -0.9], 1)
        y2 = lsim(G, u, t)[1][:]
        @test all(y .== y2)

        na = 1
        yr, A = getARregressor(y, na)
        @test length(yr) == N - na
        @test size(A) == (N - na, na)

        @test yr == y[na+1:end]
        @test A[:, 1] == y[1:end-na]

        d = iddata(y, 1)
        Gh = ar(d, na)
        @test Gh ≈ G # We should be able to recover this transfer function
        uh = lsim(1/Gh, y, t)[1][:]
        @test u ≈ uh atol = eps()

        N = 10000
        t = 1:N
        y = zeros(N)
        y[1] = 5randn()
        for i = 2:N
            y[i] = 0.9y[i-1] + 0.01randn()
        end
        d = iddata(y, 1)
        Gh = ar(d, na)
        @test Gh ≈ G atol = 0.02 # We should be able to recover this transfer function
        @test freqresptest(G, Gh, 0.05)
        yh = predict(Gh, y)
        @test rms(y[2:end] - yh) < 0.0102

        Gh2 = ar(d, na, stochastic = true)
        @test denvec(Gh2)[1][end] ≈ denvec(Gh)[1][end]
    end

    @testset "arxar" begin
        @info "Testing arxar"
        # Test examples are taken from Söderstöms paper and compared against it (this is also where the (high) tolerances come from)
        Random.seed!(0)
        N = 500 # Number of samples used for simulation by Söderström
        time = 1:N
        sim(G, u) = lsim(G, u, time)[1][:]

        #### S1 ####
        A = tf([1, -0.8], [1, 0], 1)
        B = tf([0, 1], [1, 0], 1)
        G = minreal(B / A)
        D = tf([1, 0.7], [1, 0], 1)
        H = minreal(1 / (D * A))

        u = rand(Normal(0, 1), N)
        e = rand(Normal(0, 1), N)
        y = sim(G, u)
        v = sim(H, e)
        yv = y.+ v
        d = iddata(yv, u, 1)
        ###########
        na, nb , nd = 1, 1, 1
        Gest, Hest, res = arxar(d, na, nb, nd, iterations = 10, verbose = true, δmin = 1e-3)
        @test isapprox(Gest, G, atol = 1e-1)
        @test isapprox(Hest, 1/D, atol = 1e-1)
        @test var(res .- e[2:end]) < 1e-1

        #### S12 ####
        A = tf([1, -0.7], [1, 0], 1)
        B = tf([0, 1], [1, 0], 1)
        G = minreal(B / A)
        D = tf([1, 0.9], [1, 0], 1)
        H = minreal(1 / (D * A))

        u = rand(Normal(0, 1), N)
        e = rand(Normal(0, sqrt(1.2)), N)
        y = sim(G, u)
        v = sim(H, e)
        yv = y.+ v
        d = iddata(yv, u, 1)
        ############
        na, nb , nd = 1, 1, 1
        Gest, Hest, res = arxar(d, na, nb, nd, iterations = 10, verbose = true, δmin = 1e-3)
        @test isapprox(Gest, G, atol = 2e-1)
        @test isapprox(Hest, 1/D, atol = 1e-1)

        #### S2 #### structure of Ay = Bu + Ce
        A = tf([1, -0.8], [1, 0], 1)
        B = tf([0, 1], [1, 0], 1)
        G = minreal(B / A)
        D = tf([1, 0.7], [1, 0], 1)
        H = minreal((D / A))

        u = rand(Normal(0, 1), N)
        e = rand(Normal(0, 0.1), N)
        y = sim(G, u)
        v = sim(H, e)
        yv = y .+ v
        d = iddata(yv, u, 1)
        ############
        na, nb , nd = 1, 1, 1
        Gest, Hest, res = arxar(d, na, nb, nd, iterations = 10, verbose = true, δmin = 1e-3)
        @test isapprox(Gest, G, atol = 1e-1)
        @test isapprox(Hest, 1/tf([1, -0.49], [1, 0], 1), atol = 1e-1)

        #### S10 #### prior knowledge neccessary for identification
        N = 2000
        A = tf([1, -0.5], [1, 0], 1)
        B = tf([0, 1], [1, 0], 1)
        G = minreal(B / A)
        D = tf([1, 0.5], [1,0], 1)
        H = minreal(1 / (D * A))

        u = rand(Normal(0, 1), N)
        e = rand(Normal(0, 5), N)
        y = lsim(G, u, 1:N)[1][:]
        v = lsim(H, e, 1:N)[1][:]
        yv = y.+ v
        d = iddata(yv, u, 1)
        ############
        na, nb , nd = 1, 1, 1
        Gest, Hest, res = arxar(d, na, nb, nd, H = 1/D, iterations = 10, verbose = true, δmin = 1e-3)
        @test isapprox(Gest, G, atol = 4e-1)
        @test isapprox(Hest, 1/D, atol = 1e-1)
        @test freqresptest(G, Gest) < 0.5

        # MISO
        N = 500
        A = tf([1, -0.8], [1, 0], 1)
        B1 = tf([0, 1], [1, 0], 1)
        B2 = tf([0, -1], [1, 0], 1)
        G1 = minreal(B1 / A)
        G2 = minreal(B2 / A)
        G = [G1 G2]
        D = tf([1, 0.7], [1, 0], 1)
        H = minreal(1 / (D * A))

        u1 = rand(Normal(0, 1), N)
        u2 = rand(Normal(0, 1), N)
        u = [u1 u2]
        e = rand(Normal(0, 1), N)
        y = sim(G, u)
        v = sim(H, e)
        yv = y.+ v
        d = iddata(yv, u', 1)
        ###########
        na, nb , nd = 1, [1, 1], 1
        Gest, Hest, res = arxar(d, na, nb, nd, iterations = 10, verbose = true, δmin = 1e-3)
        @test isapprox(Gest, G, atol = 3e-1)
        @test isapprox(Hest, 1/D, atol = 3e-1)

        # inputdelay
        A = tf([1, -0.8], [1, 0], 1)
        B = tf([0, 1], [1, 0, 0], 1)
        G = minreal(B / A)
        D = tf([1, 0.7], [1, 0], 1)
        H = minreal(1 / (D * A))

        u = rand(Normal(0, 1), N)
        e = rand(Normal(0, 1), N)
        y = sim(G, u)
        v = sim(H, e)
        yv = y.+ v
        d = iddata(yv, u, 1)
        ###########
        na, nb , nd, inputdelay = 1, 1, 1, 2
        Gest, Hest, res = arxar(d, na, nb, nd, inputdelay = inputdelay, iterations = 10, verbose = true, δmin = 1e-3)
        @test isapprox(Gest, G, atol = 1e-1)
        @test isapprox(Hest, 1/D, atol = 1e-1)
        @test var(res .- e[3:end]) < 1e-1
    end

    @testset "plr" begin

        N = 2000
        t = 1:N
        u = randn(N)
        G = tf(0.8, [1, -0.9], 1)
        y = lsim(G, u, t)[1][:]
        e = randn(N)
        yn = y + e

        na, nb, nc = 1, 1, 1
        d = iddata(yn, u, 1)
        find_na(y, 6)
        find_nanb(d, 6, 6)
        Gls = arx(d, na, nb)
        Gtls = arx(d, na, nb, estimator = tls)
        Gwtls = arx(d, na, nb, estimator = wtls_estimator(y, na, nb))
        Gplr, Gn = ControlSystemIdentification.plr(d, na, nb, nc, initial_order = 10)

        # @show Gplr, Gn

        @test freqresptest(G, Gls) < 1.5
        @test freqresptest(G, Gtls) < 1
        @test freqresptest(G, Gwtls) < 0.1
        @test freqresptest(G, Gplr) < 0.1

        d = iddata(y, u, 1)
        Gls = arx(d, na, nb)
        Gtls = arx(d, na, nb, estimator = tls)
        Gwtls = arx(d, na, nb, estimator = wtls_estimator(y, na, nb))
        Gplr, Gn = ControlSystemIdentification.plr(d, na, nb, nc, initial_order = 10)

        @test freqresptest(G, Gls) < sqrt(eps())
        @test freqresptest(G, Gtls) < sqrt(eps())
        @test freqresptest(G, Gwtls) < sqrt(eps())
        @test freqresptest(G, Gplr) < sqrt(eps())

    end

    @testset "arma" begin
        @info "Testing arma"

        N  = 2000     # Number of time steps
        t  = 1:N
        Δt = 1        # Sample time
        u  = randn(N) # A random control input
        a  = 0.9

        G = tf([1, 0.1], [1, -a, 0], 1)
        y, t, x = lsim(G, u, 1:N) .|> vec

        na, nc = 2, 2   # Number of polynomial coefficients
        e = 0.0001randn(N) #+ 20randn(N) .* (rand(N) .< 0.01)
        yn = y + e    # Measurement signal with noise

        d = iddata(yn, Δt)
        model = arma(d, na, nc, initial_order = 20)

        @test numvec(model)[1] ≈ numvec(G)[1] atol = 0.5
        @test denvec(model)[1] ≈ denvec(G)[1] atol = 0.5
        @test freqresptest(G, model) < 0.3

        uh = ControlSystemIdentification.estimate_residuals(model, yn)
        @show mean(abs2, uh - u) / mean(abs2, u)
        @test mean(abs2, uh - u) / mean(abs2, u) < 0.01

    end

    @testset "arma ssa" begin
        @info "Testing arma ssa"

        T = 1000
        G = tf(1, [1, 2 * 0.1 * 1, 1])
        G = c2d(G, 1)
        u = randn(T)
        y = lsim(G, u, 1:T)[1][:]
        d = iddata(y)
        model = ControlSystemIdentification.arma_ssa(d, 2, 2, L = 200)

        @test numvec(model)[1] ≈ numvec(G)[1] atol = 0.7
        @test denvec(model)[1] ≈ denvec(G)[1] atol = 0.2

        @test freqresptest(G, model) < 2
    end

    # end
    @testset "frd" begin
        Random.seed!(1)
        ##
        T           = 100000
        h           = 1
        t           = range(0, step = h, length = T)
        sim(sys, u) = lsim(sys, u, t)[1][:]
        σy          = 0.5
        sys         = tf(1, [1, 2 * 0.1, 0.1])
        ωn          = sqrt(0.3)
        sysn        = tf(σy * ωn, [1, 2 * 0.1 * ωn, ωn^2])

        u  = randn(T)
        y  = sim(sys, u)
        yn = y + sim(sysn, randn(size(u)))
        d  = iddata(y, u, h)
        dn = iddata(yn, u, 1)

        # using BenchmarkTools
        # @btime begin
        # Random.seed!(0)
        k = coherence(d)
        @test all(k.r .> 0.99)
        k = coherence(dn)
        @test all(k.r[1:10] .> 0.9)
        @test k.r[end] .> 0.8
        @test k.r[findfirst(k.w .> ωn)] < 0.6
        G, N = tfest(dn, 0.02)
        noisemodel = innovation_form(ss(sys), syse = ss(sysn))
        noisemodel.D .*= 0
        bodeplot(
            [sys, sysn],
            exp10.(range(-3, stop = log10(pi), length = 200)),
            layout    = (1, 4),
            plotphase = false,
            subplot   = [1, 2],
            size      = (3 * 800, 600),
            linecolor = :blue,
        )#, ylims=(0.1,300))

        coherenceplot!(dn, subplot = 3)
        crosscorplot!(dn, -10:100, subplot = 4)
        plot!(G, subplot = 1, lab = "G Est", alpha = 0.3, title = "Process model")
        plot!(√N, subplot = 2, lab = "N Est", alpha = 0.3, title = "Noise model")

        for op in (+, -, *)
            @test op(G, G) isa FRD
        end

        ω = exp10.(LinRange(-2, 2, 100))
        P, C = tf(1.0, [1, 1]), pid(kp = 1, ki = 1)
        S, D, N, T = gangoffour(P, C)
        S2, D2, N2, T2 = gangoffour(FRD(ω, P), FRD(ω, C))

        @test ninputs(S2) == noutputs(S2) == 1
        @test size(S2) == (1, 1)
        @test S2 == S2
        @test FRD(ω, S) ≈ S2
        @test FRD(ω, D) ≈ D2
        @test FRD(ω, N) ≈ N2
        @test FRD(ω, T) ≈ T2

        @test FRD(ω, P)*P ≈ FRD(ω, P)*FRD(ω, P)
        @test FRD(ω, P)-P ≈ FRD(ω, P)-FRD(ω, P)
        @test FRD(ω, P)+P ≈ FRD(ω, P)+FRD(ω, P)

    end

    @testset "plots and difficult" begin

        Random.seed!(1)
        ##
        T = 200
        nx = 2
        nu = 1
        ny = 1
        x0 = randn(nx)
        σy = 0.5
        sim(sys, u) = lsim(sys, u', 1:T)[1]'
        sys = tf(1, [1, 2 * 0.1, 0.1])
        sysn = tf(σy, [1, 2 * 0.1, 0.3])

        u  = randn(nu, T)
        un = u + 0.1randn(size(u))
        y  = sim(sys, u)
        yn = y + sim(sysn, σy * randn(size(u)))
        dd = iddata(yn, un, 1)

        uv  = randn(nu, T)
        yv  = sim(sys, uv)
        ynv = yv + sim(sysn, σy * randn(size(uv)))
        dv  = iddata(yv, uv, 1)
        dnv = iddata(ynv, uv, 1)
        ##

        res = [
            pem(dnv, nx = nx, iterations = 1000, difficult = true, focus = :prediction)
            for nx in [1, 3, 4]
        ]

        ω = exp10.(range(-2, stop = log10(pi), length = 150))
        fig = plot(layout = 4, size = (1000, 600))
        for i in eachindex(res)
            (sysh, x0h, opt) = res[i]
            ControlSystemIdentification.simplot!(
                sysh,
                dnv,
                x0h;
                subplot = 1,
                ploty = i == 1,
            )
            ControlSystemIdentification.predplot!(
                sysh,
                dnv,
                x0h;
                subplot = 2,
                ploty = i == 1,
            )
        end
        bodeplot!(
            ss.(getindex.(res, 1)),
            ω,
            plotphase = false,
            subplot = 3,
            title = "Process",
            linewidth = 2 * [4 3 2 1],
        )
        bodeplot!(
            ControlSystems.innovation_form.(getindex.(res, 1)),
            ω,
            plotphase = false,
            subplot = 4,
            linewidth = 2 * [4 3 2 1],
        )
        bodeplot!(
            sys,
            ω,
            plotphase = false,
            subplot = 3,
            lab = "True",
            linecolor = :blue,
            l = :dash,
            legend = :bottomleft,
            title = "System model",
        )
        bodeplot!(
            ControlSystems.innovation_form(ss(sys), syse = ss(sysn), R2 = σy^2 * I),
            ω,
            plotphase = false,
            subplot = 4,
            lab = "True",
            linecolor = :blue,
            l = :dash,
            ylims = (0.1, 100),
            legend = :bottomleft,
            title = "Noise model",
        )
        display(fig)

    end

    @testset "impulseest" begin
        T = 200
        h = 0.1
        t = h:h:T
        sim(sys, u) = lsim(sys, u, t)[1][:]
        sys = c2d(tf(1, [1, 2 * 0.1, 0.1]), h)

        u = randn(length(t))
        y = sim(sys, u) + 0.1randn(length(t))

        d = iddata(y, u, h)
        impulseestplot(d, Int(50 / h), λ = 0)
        impulseplot!(sys, 50, l = (:dash, :blue))
    end

    @testset "Spectrogram" begin
        @info "Testing Spectrogram"

        T = 1000
        s = sin.((1:T) .* 2pi / 10)

        S1 = spectrogram(s, window = hanning)

        estimator = model_spectrum(ar, 1, 2)
        S2 = spectrogram(s, estimator, window = hanning)
        @test maximum(
            findmax(S1.power, dims = 1)[2][:] - findmax(S2.power, dims = 1)[2][:],
        ) <= CartesianIndex(1, 0)
        @test minimum(
            findmax(S1.power, dims = 1)[2][:] - findmax(S2.power, dims = 1)[2][:],
        ) >= CartesianIndex(-1, 0)

        estimator = model_spectrum(arma, 1, 2, 1)
        S2 = spectrogram(s, estimator, window = hanning)
        @test maximum(
            findmax(S1.power, dims = 1)[2][:] - findmax(S2.power, dims = 1)[2][:],
        ) <= CartesianIndex(1, 0)
        @test minimum(
            findmax(S1.power, dims = 1)[2][:] - findmax(S2.power, dims = 1)[2][:],
        ) >= CartesianIndex(-1, 0)
    end
end

##
# using TotalLeastSquares
# using LinearAlgebra: QRIteration
# using Random
# plotly(size=(600,400))
# N = 500
# ivec = 1:5:50
# res = tmap(ivec) do i
#     Random.seed!(i)
#     r = 10
#     A = Matrix{Float64}(I(r))
#     A[1, 1] = 1.01
#     G = ss(A, randn(r, m), randn(l, r), 0 * randn(l, m), 1)
#     r = G.nx
#     res = map(1:20) do _
#         u = randn(N, G.nu)
#         x0 = randn(r)
#         y, t, x = lsim(G, u, 1:N, x0 = x0)
#         @test sum(!isfinite, y) == 0
#         yn = y + 0.1randn(size(y))
#         d = iddata(yn, u, 1)
#
#         # res = n4sid(d, r, i = r+i, svd = x->svd(x, alg = QRIteration()))
#         res = n4sid2(d, r, i = r+i)
#         # freqresptest(G, res.sys) < 0.2 * m * l
#         w = exp10.(LinRange(-2, log10(pi), 600))
#         b1,_ = bode(G, w)
#         b2,_ = bode(res.sys, w)
#         mean(abs2, log.(b1)-log.(b2))
#     end
#     mean(res)
# end
#
# plot(ivec, res)


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

function flattenParfors(function_name, ast)
    #assert(isfunctionhead(ast))
    flatten_start = time_ns()

    @dprintln(1,"Starting flattenParfors.  function = ", function_name, " ast = ", ast)

    #bt = backtrace() ;
    #s = sprint(io->Base.show_backtrace(io, bt))
    #@dprintln(3, "from_root backtrace ")
    #@dprintln(3, s)

    start_time = time_ns()

    if isa(ast, Tuple)
        (LambdaVarInfo, body) = ast
    else
        LambdaVarInfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast)
    end

    expanded_args = Any[]
    flattenParfors(expanded_args, body.args, LambdaVarInfo)
    args = expanded_args

    if print_times
    @dprintln(1,"Flattening parfor bodies time = ", ns_to_sec(time_ns() - flatten_start))
    end

    @dprintln(3, "After flattening")
    for j = 1:length(args)
        @dprintln(3, args[j])
    end

    if shortcut_array_assignment != 0
        new_lives = CompilerTools.LivenessAnalysis.from_lambda(LambdaVarInfo, args, pir_live_cb, LambdaVarInfo)

        for i = 1:length(args)
            node = args[i]
            if isAssignmentNode(node)
                lhs = node.args[1]
                rhs = node.args[2]
                @dprintln(3,"shortcut_array_assignment = ", node)
                if isa(lhs, TypedVar) && isArrayType(lhs) && isa(rhs, TypedVar)
                    @dprintln(3,"shortcut_array_assignment to array detected")
                    live_info = CompilerTools.LivenessAnalysis.find_top_number(i, new_lives)
                    if !in(toLHSVar(rhs), live_info.live_out)
                        @dprintln(3,"rhs is dead")
                        # The RHS of the assignment is not live out so we can do a special assignment where the j2c_array for the LHS takes over the RHS and the RHS is nulled.
                        push!(node.args, RhsDead())
                    end
                end
            end
        end
    end

    return LambdaVarInfo, CompilerTools.LambdaHandling.getBody(args, CompilerTools.LambdaHandling.getReturnType(LambdaVarInfo))
end

function flattenParfors(out_body :: Array{Any,1}, in_body :: Array{Any,1}, linfo :: LambdaVarInfo)
    for i = 1:length(in_body)
        @dprintln(3,"Flatten index ", i, " ", in_body[i], " type = ", typeof(in_body[i]))
        if isBareParfor(in_body[i])
            # flattern nested DelayedFunc if any
            for red in in_body[i].args[1].reductions
                if isa(red.reductionVarInit, DelayedFunc)
                    func = red.reductionVarInit
                    body = Any[]
                    flattenParfors(body, func.args[1], linfo)
                    func.args[1] = body
                end
                if isa(red.reductionFunc, DelayedFunc)
                    func = red.reductionFunc
                    body = Any[]
                    flattenParfors(body, func.args[1], linfo)
                    func.args[1] = body
                end
            end
            flattenParfor(out_body, in_body[i].args[1], linfo)
        else
            push!(out_body, in_body[i])
        end
    end
end

"""
Takes a new array of body statements in the process of construction in "new_body" and takes a parfor to add to that
body.  This parfor is in the nested (parfor code is in the parfor node itself) temporary form we use for fusion although
pre-statements and post-statements are already elevated by this point.  We replace this nested form with a non-nested
form where we have a parfor_start and parfor_end to delineate the parfor code.
"""
function flattenParfor(new_body, the_parfor :: ParallelAccelerator.ParallelIR.PIRParForAst, linfo :: LambdaVarInfo)
    @dprintln(2,"Flattening ", the_parfor)

    private_set = getPrivateSet(the_parfor.body, linfo)
    private_array = collect(private_set)

    # Output to the new body that this is the start of a parfor.
    push!(new_body, TypedExpr(Int64, :parfor_start, PIRParForStartEnd(the_parfor.loopNests, the_parfor.reductions, the_parfor.instruction_count_expr, private_array)))
    # Output the body of the parfor as top-level statements in the new function body and convert any other parfors we may find.
    flattenParfors(new_body, the_parfor.body, linfo)
    # Output to the new body that this is the end of a parfor.
    push!(new_body, TypedExpr(Int64, :parfor_end, PIRParForStartEnd(deepcopy(the_parfor.loopNests), deepcopy(the_parfor.reductions), deepcopy(the_parfor.instruction_count_expr), deepcopy(private_array))))
    nothing
end


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

"""
    LearningProblem(sourcedata, targetdata, task)

A spatial learning problem with `sourcedata`, `targetdata`,
and learning `task`.

## Examples

Create a clustering problem based on a set of soil features:

```julia
julia> LearningProblem(sourcedata, targetdata,
                       ClusteringTask((:moisture,:mineral)))
```
"""
struct LearningProblem{DΩₛ<:AbstractData,
                       DΩₜ<:AbstractData,
                       T<:AbstractLearningTask} <: AbstractProblem
  sourcedata::DΩₛ
  targetdata::DΩₜ
  task::T

  function LearningProblem{DΩₛ,DΩₜ,T}(sourcedata, targetdata,
                                      task) where {DΩₛ<:AbstractData,
                                                   DΩₜ<:AbstractData,
                                                   T<:AbstractLearningTask}
    sourcevars = keys(variables(sourcedata))
    targetvars = keys(variables(targetdata))

    # assert spatial configuration
    @assert ndims(sourcedata) == ndims(targetdata) "source and target data must have the same number of dimensions"
    @assert coordtype(sourcedata) == coordtype(targetdata) "source and target data must have the same coordinate type"

    # assert task is compatible with the data
    if iscomposite(task)
      # TODO
    else
      @assert features(task) ⊆ sourcevars "features must be present in source data"
      @assert features(task) ⊆ targetvars "features must be present in target data"
      if issupervised(task)
        @assert label(task) ∈ sourcevars "label must be present in source data"
      end
    end

    new(sourcedata, targetdata, task)
  end
end

LearningProblem(sdata::DΩₛ, tdata::DΩₜ, task::T) where {DΩₛ<:AbstractData,
                                                        DΩₜ<:AbstractData,
                                                        T<:AbstractLearningTask} =
  LearningProblem{DΩₛ,DΩₜ,T}(sdata, tdata, task)

"""
    sourcedata(problem)

Return the source data of the learning `problem`.
"""
sourcedata(problem::LearningProblem) = problem.sourcedata

"""
    targetdata(problem)

Return the target data of the learning `problem`.
"""
targetdata(problem::LearningProblem) = problem.targetdata

"""
    task(problem)

Return the learning task of the learning `problem`.
"""
task(problem::LearningProblem) = problem.task

# ------------
# IO methods
# ------------
function Base.show(io::IO, problem::LearningProblem)
  dim = ndims(problem.sourcedata)
  print(io, "$(dim)D LearningProblem")
end

function Base.show(io::IO, ::MIME"text/plain", problem::LearningProblem)
  println(io, problem)
  println(io, "  source")
  println(io, "    └─data: ", problem.sourcedata)
  println(io, "  target")
  println(io, "    └─data: ", problem.targetdata)
  println(io, "  task: ", problem.task)
end


text:
# LANGUAGE: Julia
# AUTHOR: Akram Rameez
# GITHUB: https://github.com/akram-rameez

println("Hello World!")


text:
players=477;
maxn=7085100;
turn=0;
circle=[];
i=0;
position=2;
score=zeros(Int64,players);
push!(circle,0);
i=i+1;
turn=turn+1;
push!(circle,i);
i=i+1;
turn=turn+1;
splice!(circle,2,[2,1])
i=i+1;
position=2;
while i<=maxn
    global turn= turn+1;
    if (i%23)==0
        global position = position-7;
        if position<0
            global position= length(circle)+position;
        end
        if (turn%players) ==0
            global score[players]= score[players]+i+circle[position];
            splice!(circle,position);
        else
            global score[(turn%players)]= score[(turn%players)]+i+circle[position];
            splice!(circle,position)
        end
    else
        global position=position+2;
        if length(circle)+1<position
           position=position-length(circle);
        end
        if length(circle)+1==position
            push!(circle,i);
        else
            splice!(circle,position,[i,circle[position]]);
        end
    end
    global i=i+1;
end
max = 0;
for i in score
    if i>max
        global max = i;
    end
end
println(max);


text:
module SemiclassicalOrthogonalPolynomials
using ClassicalOrthogonalPolynomials, FillArrays, LazyArrays, ArrayLayouts, QuasiArrays, InfiniteArrays, ContinuumArrays, LinearAlgebra, BandedMatrices,
        SpecialFunctions, HypergeometricFunctions

import Base: getindex, axes, size, \, /, *, +, -, summary, ==, copy, sum, unsafe_getindex, convert, OneTo

import ArrayLayouts: MemoryLayout, ldiv, diagonaldata, subdiagonaldata, supdiagonaldata
import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedLayout, _BandedMatrix
import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector, AccumulateAbstractVector, LazyVector, AbstractCachedMatrix
import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight, NormalizedBasisLayout,
                                        Bidiagonal, Tridiagonal, SymTridiagonal, symtridiagonalize, normalizationconstant, LanczosPolynomial,
                                        OrthogonalPolynomialRatio, Weighted, Expansion, UnionDomain, oneto, Hilbert, WeightedOrthogonalPolynomial, HalfWeighted,
                                        Associated
import InfiniteArrays: OneToInf, InfUnitRange
import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid
import FillArrays: SquareEye
import HypergeometricFunctions: _₂F₁general2

export LanczosPolynomial, Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, OrthogonalPolynomialRatio, TwoBandJacobi, TwoBandWeight

""""
    SemiclassicalJacobiWeight(t, a, b, c)

is a quasi-vector corresponding to the weight `x^a * (1-x)^b * (t-x)^c` on `x in 0..1`
"""
struct SemiclassicalJacobiWeight{T} <: Weight{T}
    t::T
    a::T
    b::T
    c::T
end

function SemiclassicalJacobiWeight(t, a, b, c)
    T = promote_type(eltype(t), eltype(a), eltype(b), eltype(c))
    SemiclassicalJacobiWeight{T}(t, a, b, c)
end

copy(w::SemiclassicalJacobiWeight) = w

axes(P::SemiclassicalJacobiWeight{T}) where T = (Inclusion(UnitInterval{T}()),)
function getindex(P::SemiclassicalJacobiWeight, x::Real)
    t,a,b,c = P.t,P.a,P.b,P.c
    @boundscheck checkbounds(P, x)
    x^a * (1-x)^b * (t-x)^c
end

function sum(P::SemiclassicalJacobiWeight{T}) where T
    (t,a,b,c) = BigFloat.((P.t,P.a,P.b,P.c))
    return convert(T, t^c * beta(1+a,1+b) * _₂F₁general2(1+a,-c,2+a+b,1/t))
end

function summary(io::IO, P::SemiclassicalJacobiWeight)
    t,a,b,c = P.t,P.a,P.b,P.c
    print(io, "x^$a * (1-x)^$b * ($t-x)^$c")
end

function ==(A::SemiclassicalJacobiWeight, B::SemiclassicalJacobiWeight)
    A.a == B.a && A.b == B.b &&
        ((iszero(A.c) && iszero(B.c)) ||
         (A.t == B.t && A.c == B.c))
end

function Expansion(w::SemiclassicalJacobiWeight{T}) where T
    t,a,b,c = w.t,w.a,w.b,w.c
    P = jacobi(b, a, UnitInterval{T}())
    x = axes(P,1)
    # TODO: simplify
    LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), P).w
end

==(A::SemiclassicalJacobiWeight, B::Expansion) = Expansion(A) == B
==(A::Expansion, B::SemiclassicalJacobiWeight) = A == Expansion(B)

"""
   RaisedOP(P, y)

gives the OPs w.r.t. (y - x) .* w based on lowering to Q.
"""

struct RaisedOP{T, QQ, LL<:OrthogonalPolynomialRatio} <: OrthogonalPolynomial{T}
    Q::QQ
    ℓ::LL
end

RaisedOP(Q, ℓ::OrthogonalPolynomialRatio) = RaisedOP{eltype(Q),typeof(Q),typeof(ℓ)}(Q, ℓ)
RaisedOP(Q, y::Number) = RaisedOP(Q, OrthogonalPolynomialRatio(Q,y))

function jacobimatrix(P::RaisedOP{T}) where T
    ℓ = P.ℓ
    X = jacobimatrix(P.Q)
    a,b = diagonaldata(X), supdiagonaldata(X)
    # non-normalized lower diag of Jacobi
    v = Vcat(zero(T),b .* ℓ)
    c = BroadcastVector((ℓ,a,b,sa,v) -> ℓ*a + b - b*ℓ^2 - sa*ℓ + ℓ*v, ℓ, a, b, a[2:∞], v)
    Tridiagonal(c, BroadcastVector((ℓ,a,b,v) -> a - b * ℓ + v, ℓ,a,b,v), b)
end

"""
  the bands of the Jacobi matrix
"""
mutable struct SemiclassicalJacobiBand{dv,T} <: AbstractCachedVector{T}
    data::Vector{T}
    a::AbstractVector{T}
    b::AbstractVector{T}
    ℓ::AbstractVector{T}
    datasize::Tuple{Int}
end

function LazyArrays.cache_filldata!(r::SemiclassicalJacobiBand{:dv}, inds::AbstractUnitRange)
    rℓ = r.ℓ[inds[1]-1:inds[end]]; ℓ,sℓ = rℓ[2:end], rℓ[1:end-1]
    ra = r.a[inds[1]:(inds[end]+1)]; a = ra[1:end-1]; sa = ra[2:end]
    rb = r.b[inds[1]-1:inds[end]]; b = rb[2:end]; sb = rb[1:end-1]
    r.data[inds] .= @.(a - b * ℓ + sb*sℓ)
end

function LazyArrays.cache_filldata!(r::SemiclassicalJacobiBand{:ev}, inds::AbstractUnitRange)
    rℓ = r.ℓ[inds[1]-1:inds[end]]; ℓ,sℓ = rℓ[2:end], rℓ[1:end-1]
    ra = r.a[inds[1]:(inds[end]+1)]; a = ra[1:end-1]; sa = ra[2:end]
    rb = r.b[inds[1]-1:inds[end]]; b = rb[2:end]; sb = rb[1:end-1]
    r.data[inds] .= @.(sqrt((ℓ*a + b - b*ℓ^2 - sa*ℓ + ℓ*sb*sℓ)*b))
end

size(::SemiclassicalJacobiBand) = (ℵ₀,)

SemiclassicalJacobiBand{:dv,T}(a,b,ℓ) where T = SemiclassicalJacobiBand{:dv,T}(T[a[1] - b[1]ℓ[1]], a, b, ℓ, (1,))
SemiclassicalJacobiBand{:ev,T}(a,b,ℓ) where T = SemiclassicalJacobiBand{:ev,T}(T[sqrt((ℓ[1]*a[1] + b[1] - b[1]*ℓ[1]^2 - a[2]*ℓ[1])b[1])], a, b, ℓ, (1,))
SemiclassicalJacobiBand{dv}(a,b,ℓ) where dv = SemiclassicalJacobiBand{dv,promote_type(eltype(a),eltype(b),eltype(ℓ))}(a,b,ℓ)

copy(r::SemiclassicalJacobiBand) = r # immutable

"""
   SemiclassicalJacobi(t, a, b, c)

is a quasi-matrix for the  orthogonal polynomials w.r.t. x^a * (1-x)^b * (t-x)^c on [0,1]
"""
struct SemiclassicalJacobi{T} <: OrthogonalPolynomial{T}
    t::T
    a::T
    b::T
    c::T
    X::AbstractMatrix{T}
    SemiclassicalJacobi{T}(t::T,a::T,b::T,c::T,X::AbstractMatrix{T}) where T = new{T}(t,a,b,c,X)
end

const WeightedSemiclassicalJacobi{T} = WeightedBasis{T,<:SemiclassicalJacobiWeight,<:SemiclassicalJacobi}

function SemiclassicalJacobi(t, a, b, c, X::AbstractMatrix)
    T = float(promote_type(typeof(t), typeof(a), typeof(b), typeof(c), eltype(X)))
    SemiclassicalJacobi{T}(T(t), T(a), T(b), T(c), convert(AbstractMatrix{T},X))
end

SemiclassicalJacobi(t, a, b, c, P::SemiclassicalJacobi) = SemiclassicalJacobi(t, a, b, c, semiclassical_jacobimatrix(P, a, b, c))
SemiclassicalJacobi(t, a, b, c) = SemiclassicalJacobi(t, a, b, c, semiclassical_jacobimatrix(t, a, b, c))
SemiclassicalJacobi{T}(t, a, b, c) where T = SemiclassicalJacobi(convert(T,t), convert(T,a), convert(T,b), convert(T,b))

WeightedSemiclassicalJacobi(t, a, b, c, P...) = SemiclassicalJacobiWeight(t, a, b, c) .* SemiclassicalJacobi(t, a, b, c, P...)

function semiclassical_jacobimatrix(t, a, b, c)
    T = float(promote_type(typeof(t), typeof(a), typeof(b), typeof(c)))
    if a == 0 && b == 0 && c == -1
        # for this special case we can generate the Jacobi operator explicitly
        N = (1:∞)
        α = T.(neg1c_αcfs(t)) # cached coefficients
        A = Vcat((α[1]+1)/2 , -N./(N.*4 .- 2).*α .+ (N.+1)./(N.*4 .+ 2).*α[2:end].+1/2)
        C = -(N)./(N.*4 .- 2)
        B = Vcat((α[1]^2*3-α[1]*α[2]*2-1)/6 , -(N)./(N.*4 .+ 2).*α[2:end]./α)
        # if J is Tridiagonal(c,a,b) then for norm. OPs it becomes SymTridiagonal(a, sqrt.( b.* c))
        return SymTridiagonal(A, sqrt.(B.*C))
    end
    P = jacobi(b, a, UnitInterval{T}())
    iszero(c) && return symtridiagonalize(jacobimatrix(P))
    x = axes(P,1)
    jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), P))
end

function symraised_jacobimatrix(Q, y)
    ℓ = OrthogonalPolynomialRatio(Q,y)
    X = jacobimatrix(Q)
    a,b = diagonaldata(X), supdiagonaldata(X)
    SymTridiagonal(SemiclassicalJacobiBand{:dv}(a,b,ℓ), SemiclassicalJacobiBand{:ev}(a,b,ℓ))
end

function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
    if  iszero(a) && iszero(b) && c == -1 # (a,b,c) = (0,0,-1) special case
        semiclassical_jacobimatrix(Q.t, a, b, c)
    elseif a == Q.a+1 && b == Q.b && c == Q.c  # raising by 1
        symraised_jacobimatrix(Q, 0)
    elseif a == Q.a && b == Q.b+1 && c == Q.c
        symraised_jacobimatrix(Q, 1)
    elseif a == Q.a && b == Q.b && c == Q.c+1
        symraised_jacobimatrix(Q, Q.t)
    elseif a == Q.a && b == Q.b && c == Q.c-1 # lowering by 1
        symlowered_jacobimatrix(Q, :c)
    elseif a == Q.a-1 && b == Q.b && c == Q.c
        symlowered_jacobimatrix(Q, :a)
    elseif a == Q.a && b == Q.b-1 && c == Q.c
        symlowered_jacobimatrix(Q, :b)
    elseif a > Q.a  # iterative raising
        semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a+1, Q.b, Q.c, Q), a, b, c)
    elseif b > Q.b
        semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b+1, Q.c, Q), a, b, c)
    elseif c > Q.c
        semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1, Q), a, b, c)
    elseif b < Q.b  # iterative lowering
        semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c)
    elseif c < Q.c
        semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c-1, Q), a, b, c)
    elseif a < Q.a
        semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c)
    elseif a == Q.a && b == Q.b && c == Q.c # same basis
        jacobimatrix(Q)
    else
        error("Not Implemented")
    end
end

LanczosPolynomial(P::SemiclassicalJacobi{T}) where T =
    LanczosPolynomial(orthogonalityweight(P), Normalized(jacobi(P.b, P.a, UnitInterval{T}())), P.X.dv.data)

"""
    toclassical(P::SemiclassicalJacobi)

gives either a mapped `Jacobi` or `LanczosPolynomial` version of `P`.
"""
toclassical(P::SemiclassicalJacobi{T}) where T = iszero(P.c) ? Normalized(jacobi(P.b, P.a, UnitInterval{T}())) : LanczosPolynomial(P)

copy(P::SemiclassicalJacobi) = P
axes(P::SemiclassicalJacobi{T}) where T = (Inclusion(UnitInterval{T}()),OneToInf())

==(A::SemiclassicalJacobi, B::SemiclassicalJacobi) = A.t == B.t && A.a == B.a && A.b == B.b && A.c == B.c
==(::AbstractQuasiMatrix, ::SemiclassicalJacobi) = false
==(::SemiclassicalJacobi, ::AbstractQuasiMatrix) = false

orthogonalityweight(P::SemiclassicalJacobi) = SemiclassicalJacobiWeight(P.t, P.a, P.b, P.c)

function summary(io::IO, P::SemiclassicalJacobi)
    t,a,b,c = P.t,P.a,P.b,P.c
    print(io, "SemiclassicalJacobi with weight x^$a * (1-x)^$b * ($t-x)^$c")
end

jacobimatrix(P::SemiclassicalJacobi) = P.X

"""
    op_lowering(Q, y)
Gives the Lowering operator from OPs w.r.t. (x-y)*w(x) to Q
as constructed from Chistoffel–Darboux
"""
function op_lowering(Q, y)
    # we first use Christoff-Darboux with d = 1
    # But we want the first OP to be 1 so we rescale
    P = RaisedOP(Q, y)
    A,_,_ = recurrencecoefficients(Q)
    d = -inv(A[1]*_p0(Q)*P.ℓ[1])
    κ = d * normalizationconstant(1, P)
    T = eltype(κ)
    # hide array type for compilation
    Bidiagonal(κ, -(κ .* P.ℓ), :L)
end

function semijacobi_ldiv(Q, P::SemiclassicalJacobi)
    if P.a ≤ 0 && P.b ≤ 0 && P.c ≤ 0
        P̃ = toclassical(P)
        (Q \ P̃)/_p0(P̃)
    else
        error("Implement")
    end
end

function semijacobi_ldiv(P::SemiclassicalJacobi, Q)
    R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1))
    R̃ = toclassical(R)
    (P \ R) * _p0(R̃) * (R̃ \ Q)
end

struct SemiclassicalJacobiLayout <: AbstractBasisLayout end
MemoryLayout(::Type{<:SemiclassicalJacobi}) = SemiclassicalJacobiLayout()

copy(L::Ldiv{<:NormalizedBasisLayout,SemiclassicalJacobiLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},SemiclassicalJacobiLayout}(L.A, L.B))
copy(L::Ldiv{SemiclassicalJacobiLayout,<:NormalizedBasisLayout}) = copy(Ldiv{SemiclassicalJacobiLayout,ApplyLayout{typeof(*)}}(L.A, L.B))

copy(L::Ldiv{ApplyLayout{typeof(*)},SemiclassicalJacobiLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},BasisLayout}(L.A, L.B))
copy(L::Ldiv{SemiclassicalJacobiLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{BasisLayout,ApplyLayout{typeof(*)}}(L.A, L.B))

copy(L::Ldiv{MappedBasisLayout,SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
copy(L::Ldiv{SemiclassicalJacobiLayout,MappedBasisLayout}) = semijacobi_ldiv(L.A, L.B)

copy(L::Ldiv{SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
copy(L::Ldiv{<:Any,SemiclassicalJacobiLayout}) = semijacobi_ldiv(L.A, L.B)
function copy(L::Ldiv{SemiclassicalJacobiLayout,SemiclassicalJacobiLayout})
    Q,P = L.A,L.B
    @assert Q.t == P.t
    Q == P && return SquareEye{eltype(L)}((axes(P,2),))
    M_Q = massmatrix(Q)
    M_P = massmatrix(P)
    L = P \ (SemiclassicalJacobiWeight(Q.t, Q.a-P.a, Q.b-P.b, Q.c-P.c) .* Q)
    inv(M_Q) * L' * M_P
end

\(A::LanczosPolynomial, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
\(A::SemiclassicalJacobi, B::LanczosPolynomial) = semijacobi_ldiv(A, B)
function \(w_A::WeightedSemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi)
    wA,A = w_A.args
    wB,B = w_B.args
    @assert wA.t == wB.t == A.t == B.t

    if wA.a == wB.a && wA.b == wB.b && wA.c == wB.c
        A \ B
    elseif wA.a+1 == wB.a && wA.b == wB.b && wA.c == wB.c
        @assert A.a+1 == B.a && A.b == B.b && A.c == B.c
        op_lowering(A,0)
    elseif wA.a == wB.a && wA.b+1 == wB.b && wA.c == wB.c
        @assert A.a == B.a && A.b+1 == B.b && A.c == B.c
        -op_lowering(A,1)
    elseif wA.a == wB.a && wA.b == wB.b && wA.c+1 == wB.c
        @assert A.a == B.a && A.b == B.b && A.c+1 == B.c
        -op_lowering(A,A.t)
    elseif wA.a+1 ≤ wB.a
        C = SemiclassicalJacobi(A.t, A.a+1, A.b, A.c, A)
        w_C = SemiclassicalJacobiWeight(wA.t, wA.a+1, wA.b, wA.c) .* C
        L_2 = w_C \ w_B
        L_1 = w_A \ w_C
        L_1 * L_2
    elseif wA.b+1 ≤ wB.b
        C = SemiclassicalJacobi(A.t, A.a, A.b+1, A.c, A)
        w_C = SemiclassicalJacobiWeight(wA.t, wA.a, wA.b+1, wA.c) .* C
        L_2 = w_C \ w_B
        L_1 = w_A \ w_C
        L_1 * L_2
    else
        error("Implement")
    end
end

\(A::SemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi) = (SemiclassicalJacobiWeight(A.t,0,0,0) .* A) \ w_B

massmatrix(P::SemiclassicalJacobi) = Diagonal(Fill(sum(orthogonalityweight(P)),∞))

@simplify function *(Ac::QuasiAdjoint{<:Any,<:SemiclassicalJacobi}, wB::WeightedBasis{<:Any,<:SemiclassicalJacobiWeight,<:SemiclassicalJacobi})
    A = parent(Ac)
    w,B = arguments(wB)
    P = SemiclassicalJacobi(w.t, w.a, w.b, w.c)
    (P\A)' * massmatrix(P) * (P \ B)
end

function ldiv(Q::SemiclassicalJacobi, f::AbstractQuasiVector)
    if iszero(Q.a) && iszero(Q.b) && Q.c == -1
        # todo: due to a stdlib error this won't work with bigfloat as is
        T = typeof(Q.t)
        R = Legendre{T}()[affine(Inclusion(zero(T)..one(T)), axes(Legendre{T}(),1)), :]
        B = neg1c_tolegendre(Q.t)
        return (B \ (R \ f))
    end
    R = SemiclassicalJacobi(Q.t, mod(Q.a,-1), mod(Q.b,-1), mod(Q.c,-1))
    R̃ = toclassical(R)
    return (Q \ R̃) * (R̃ \ f)
end
function ldiv(Qn::SubQuasiArray{<:Any,2,<:SemiclassicalJacobi,<:Tuple{<:Inclusion,<:Any}}, C::AbstractQuasiArray)
    _,jr = parentindices(Qn)
    Q = parent(Qn)
    if iszero(Q.a) && iszero(Q.b) && Q.c == -1
        # todo: due to a stdlib error this won't work with bigfloat as is
        T = typeof(Q.t)
        R = Legendre{T}()[affine(Inclusion(zero(T)..one(T)), axes(Legendre{T}(),1)), :]
        B = neg1c_tolegendre(Q.t)
        return (B[jr,jr] \ (R[:,jr] \ C))
    end
    R = SemiclassicalJacobi(Q.t, mod(Q.a,-1), mod(Q.b,-1), mod(Q.c,-1))
    R̃ = toclassical(R)
    (Q \ R̃)[jr,jr] * (R̃[:,jr] \ C)
end

# sqrt(1-(1-x)^2) == sqrt(2x-x^2) == sqrt(x)*sqrt(2-x)
# sqrt(1-(1-x)^2) == sqrt(2x-x^2) == sqrt(x)*sqrt(2-x)

convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:a,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, Q.P.a,zero(T),zero(T)) .* Q.P
convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:b,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, zero(T),Q.P.b,zero(T)) .* Q.P
convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:c,T,<:SemiclassicalJacobi}) where T = SemiclassicalJacobiWeight(Q.P.t, zero(T),zero(T),Q.P.c) .* Q.P

include("derivatives.jl")

###
# Hierarchy
###

function Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, ar::AbstractUnitRange, b::Number, c::Number)
    Ps = [SemiclassicalJacobi(t, first(ar), b, c)]
    for a in ar[2:end]
        push!(Ps, SemiclassicalJacobi(t, a, b, c, Ps[end]))
    end
    Ps
end

function Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, br::AbstractUnitRange, c::Number)
    Ps = [SemiclassicalJacobi(t, a, first(br), c)]
    for b in br[2:end]
        push!(Ps, SemiclassicalJacobi(t, a, b, c, Ps[end]))
    end
    Ps
end

function Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, cr::AbstractUnitRange)
    Ps = [SemiclassicalJacobi(t, a, b, first(cr))]
    for c in cr[2:end]
        push!(Ps, SemiclassicalJacobi(t, a, b, c, Ps[end]))
    end
    Ps
end

include("twoband.jl")
include("lowering.jl")

end


text:
# module Jive

struct EndError <: Exception
end

function isenderror(ex::LoadError)
    if ex.error isa EndError
        return true
    elseif ex.error isa LoadError
        return isenderror(ex.error)
    else
        return false
    end
end

"""
    @__END__

`throw(Jive.EndError())`
"""
macro __END__()
    @eval function Base.showerror(io::IO, ex::LoadError, bt; backtrace=true)
        if isenderror(ex)
            printstyled(io, "@__END__", color=:cyan)
            print(io, " at ", basename(ex.file), ":", ex.line)
        else
            print(io, "Error while loading expression starting at ", ex.file, ":", ex.line)
        end
    end
    throw(EndError())
end

# module Jive


text:
module TDD

function neighbor(n)
    a=reshape(graph[n],(1,length(graph[n])))
    b=hcat(n,a)
    println(sort(b,dims=2))
end
end