text:
"""
    FEMMDeforLinearMSModule

Module for operations on interiors of domains to construct system matrices and
system vectors for linear deformation models:  mean-strain  formulation.
"""
module FEMMDeforLinearMSModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, FESetH8, FESetT10, manifdim, nodesperelem, gradN!
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FEMMDeforLinearBaseModule: FEMMDeforLinearAbstract
import FinEtools.DeforModelRedModule: DeforModelRed, DeforModelRed3D
import FinEtools.MatDeforModule: MatDefor
import FinEtools.MatDeforElastIsoModule: MatDeforElastIso
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtools.DeforModelRedModule: nstressstrain, nthermstrain, Blmat!
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
using FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, add_btv!, loc!, jac!, locjac!
import FinEtools.FEMMDeforLinearBaseModule: stiffness, nzebcloadsstiffness, mass, thermalstrainloads, inspectintegpoints
import FinEtools.FEMMBaseModule: associategeometry!
import FinEtools.MatDeforModule: rotstressvec
import LinearAlgebra: mul!, Transpose, UpperTriangular
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: norm, qr, diag, dot, cond
import Statistics: mean

abstract type FEMMDeforLinearAbstractMS <: FEMMDeforLinearAbstract end

mutable struct FEMMDeforLinearMSH8{MR<:DeforModelRed, S<:FESetH8, F<:Function, M<:MatDefor} <: FEMMDeforLinearAbstractMS
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabilization_material::MatDeforElastIso
    phis::Vector{FFlt}
end

function FEMMDeforLinearMSH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:DeforModelRed,  S<:FESetH8, F<:Function, M<:MatDefor}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = make_stabilization_material(material)
    return FEMMDeforLinearMSH8(mr, integdomain, mcsys, material, stabilization_material, fill(zero(FFlt), 1))
end

function FEMMDeforLinearMSH8(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:DeforModelRed,  S<:FESetH8, F<:Function, M<:MatDefor}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = make_stabilization_material(material)
    return FEMMDeforLinearMSH8(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabilization_material, fill(zero(FFlt), 1))
end

mutable struct FEMMDeforLinearMST10{MR<:DeforModelRed, S<:FESetT10, F<:Function, M<:MatDefor} <: FEMMDeforLinearAbstractMS
    mr::Type{MR}
    integdomain::IntegDomain{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
    stabilization_material::MatDeforElastIso
    phis::Vector{FFlt}
end

function FEMMDeforLinearMST10(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:DeforModelRed,  S<:FESetT10, F<:Function, M<:MatDefor}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = make_stabilization_material(material)
    return FEMMDeforLinearMST10(mr, integdomain, mcsys, material, stabilization_material, fill(zero(FFlt), 1))
end

function FEMMDeforLinearMST10(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:DeforModelRed,  S<:FESetT10, F<:Function, M<:MatDefor}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (mr == DeforModelRed3D) "3D model required"
    stabilization_material = make_stabilization_material(material)
    return FEMMDeforLinearMST10(mr, integdomain, CSys(manifdim(integdomain.fes)), material, stabilization_material, fill(zero(FFlt), 1))
end

function make_stabilization_material(material::M) where {M}
    ns = fieldnames(typeof(material))
    E = 0.0; nu = 0.0
    if :E in ns
        E = material.E
        if material.nu < 0.3
            nu = material.nu
        else
            nu = 0.3 + (material.nu - 0.3) / 2.0
        end
    else
        if :E1 in ns
            E = mean([material.E1, material.E2, material.E3])
            nu = min(material.nu12, material.nu13, material.nu23)
        else
            error("No clues on how to construct the stabilization material")
        end
    end
    return  MatDeforElastIso(material.mr, 0.0, E, nu, 0.0)
end

function buffers1(self::FEMMDeforLinearAbstractMS, geom::NodalField, npts::FInt)
    fes = self.integdomain.fes
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    # Prepare buffers
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    gradN = fill(zero(FFlt), nne, mdim);
    return loc, J, csmatTJ, gradN
end

function buffers2(self::FEMMDeforLinearAbstractMS, geom::NodalField, u::NodalField, npts::FInt)
    fes = self.integdomain.fes
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    nstrs = nstressstrain(self.mr);  # number of stresses
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare buffers
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    dofnums = zeros(FInt, 1, elmatdim); # degree of freedom array -- buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    AllgradN = FFltMat[];
    for ixxxx = 1:npts
        push!(AllgradN, fill(zero(FFlt), nne, mdim));
    end
    Jac = fill(zero(FFlt), npts);
    MeangradN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    D = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    Dstab = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    DB = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    Bbar = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    elvecfix = fill(zero(FFlt), elmatdim); # vector of prescribed displ. -- buffer
    elvec = fill(zero(FFlt), elmatdim); # element vector -- buffer
    return dofnums, loc, J, csmatTJ, AllgradN, MeangradN, Jac, D, Dstab, B, DB, Bbar, elmat, elvec, elvecfix
end

function centroid!(self::F,  loc, X::FFltMat, conn::C) where {F<:FEMMDeforLinearMSH8, C}
    weights = [0.125
                0.125
                0.125
                0.125
                0.125
                0.125
                0.125
                0.125]
    return loc!(loc, X, conn, reshape(weights, 8, 1))
end

function centroid!(self::F, loc, X::FFltMat, conn::C) where {F<:FEMMDeforLinearMST10, C}
    weights = [ -0.125
                -0.125
                -0.125
                -0.125
                0.250
                0.250
                0.250
                0.250
                0.250
                0.250]
    return loc!(loc, X, conn, reshape(weights, 10, 1))
end

"""
    associategeometry!(self::FEMMAbstractBase,  geom::NodalField{FFlt})

Associate geometry field with the FEMM.

Compute the  correction factors to account for  the shape of the  elements.
"""
function associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:FEMMDeforLinearMSH8}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    loc, J, csmatTJ, gradN = buffers1(self, geom, npts)
    self.phis = fill(zero(FFlt), count(fes))
    for i = 1:count(fes) # Loop over elements
        # NOTE: the coordinate system should be evaluated at a single point within the
        # element in order for the derivatives to be consistent at all quadrature points
        loc = centroid!(self,  loc, geom.values, fes.conn[i])
        updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        for j = 1:npts # Loop over quadrature points
            jac!(J, geom.values, fes.conn[i], gradNparams[j])
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, gradN, gradNparams[j], csmatTJ);
            h2 = diag(transpose(csmatTJ)*csmatTJ)
            cap_phi = (2 * (1 + self.stabilization_material.nu) * (minimum(h2) / maximum(h2)))  # Plane stress
            phi = cap_phi / (1 + cap_phi)
            self.phis[i] = max(self.phis[i], phi)
        end # Loop over quadrature points
    end # Loop over elements
    return self
end

"""
    associategeometry!(self::FEMMAbstractBase,  geom::NodalField{FFlt})

Associate geometry field with the FEMM.

Compute the  correction factors to account for  the shape of the  elements.
"""
function associategeometry!(self::F,  geom::NodalField{FFlt}) where {F<:FEMMDeforLinearMST10}
    gamma = 2.6; C = 1e4;
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    loc, J, csmatTJ, gradN = buffers1(self, geom, npts)
    self.phis = fill(zero(FFlt), count(fes))
    for i = 1:count(fes) # Loop over elements
        # NOTE: the coordinate system should be evaluated at a single point within the
        # element in order for the derivatives to be consistent at all quadrature points
        loc = centroid!(self,  loc, geom.values, fes.conn[i])
        updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        for j = 1:npts # Loop over quadrature points
            jac!(J, geom.values, fes.conn[i], gradNparams[j])
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            condJ = cond(csmatTJ);
            cap_phi = C*(1.0/condJ)^(gamma);
            phi = cap_phi / (1 + cap_phi)
            self.phis[i] = max(self.phis[i], phi)
        end # Loop over quadrature points
    end # Loop over elements
    return self
end

"""
stiffness(self::FEMMDeforLinearAbstractMS, assembler::A,
      geom::NodalField{FFlt},
      u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}

Compute and assemble  stiffness matrix.
"""
function stiffness(self::FEMMDeforLinearAbstractMS, assembler::A,
    geom::NodalField{FFlt},
    u::NodalField{T}) where {A<:SysmatAssemblerBase, T<:Number}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, csmatTJ, AllgradN, MeangradN, Jac, D, Dstab, B, DB, Bbar, elmat, elvec, elvecfix = buffers2(self, geom, u, npts)
    realmat = self.material
    stabmat = self.stabilization_material
    realmat.tangentmoduli!(realmat, D, 0.0, 0.0, loc, 0)
    stabmat.tangentmoduli!(stabmat, Dstab, 0.0, 0.0, loc, 0)
    startassembly!(assembler, size(elmat, 1), size(elmat, 2), count(fes),
    u.nfreedofs, u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        # NOTE: the coordinate system should be evaluated at a single point within the
        # element in order for the derivatives to be consistent at all quadrature points
        loc = centroid!(self,  loc, geom.values, fes.conn[i])
        updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        vol = 0.0; # volume of the element
        fill!(MeangradN, 0.0) # mean basis function gradients
        for j = 1:npts # Loop over quadrature points
            jac!(J, geom.values, fes.conn[i], gradNparams[j])
            Jac[j] = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, AllgradN[j], gradNparams[j], csmatTJ);
            dvol = Jac[j]*w[j]
            MeangradN .= MeangradN .+ AllgradN[j]*dvol
            vol = vol + dvol
        end # Loop over quadrature points
        MeangradN .= MeangradN/vol
        Blmat!(self.mr, Bbar, Ns[1], MeangradN, loc, self.mcsys.csmat);
        fill!(elmat,  0.0); # Initialize element matrix
        add_btdb_ut_only!(elmat, Bbar, vol, D, DB)
        add_btdb_ut_only!(elmat, Bbar, -self.phis[i]*vol, Dstab, DB)
        for j = 1:npts # Loop over quadrature points
            Blmat!(self.mr, B, Ns[j], AllgradN[j], loc, self.mcsys.csmat);
            add_btdb_ut_only!(elmat, B, self.phis[i]*Jac[j]*w[j], Dstab, DB)
        end # Loop over quadrature points
        complete_lt!(elmat)
        gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
        assemble!(assembler, elmat, dofnums, dofnums); # assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

"""
nzebcloadsstiffness(self::FEMMDeforLinearAbstract,  assembler::A,
  geom::NodalField{FFlt},
  u::NodalField{T}) where {A<:SysvecAssemblerBase, T<:Number}

Compute load vector for nonzero EBC for fixed displacement.
"""
function nzebcloadsstiffness(self::FEMMDeforLinearAbstractMS,  assembler::A, geom::NodalField{FFlt}, u::NodalField{T}) where {A<:SysvecAssemblerBase, T<:Number}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, csmatTJ, AllgradN, MeangradN, Jac, D, Dstab, B, DB, Bbar, elmat, elvec, elvecfix = buffers2(self, geom, u, npts)
    realmat = self.material
    stabmat = self.stabilization_material
    realmat.tangentmoduli!(realmat, D, 0.0, 0.0, loc, 0)
    stabmat.tangentmoduli!(stabmat, Dstab, 0.0, 0.0, loc, 0)
    startassembly!(assembler,  u.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gatherfixedvalues_asvec!(u, elvecfix, fes.conn[i]);# retrieve element displacement vector
        if norm(elvecfix, Inf) != 0.0   # Is the load nonzero?
            # NOTE: the coordinate system should be evaluated at a single point within the
            # element in order for the derivatives to be consistent at all quadrature points
            loc = centroid!(self,  loc, geom.values, fes.conn[i])
            updatecsmat!(self.mcsys, loc, J, fes.label[i]);
            vol = 0.0; # volume of the element
            fill!(MeangradN, 0.0) # mean basis function gradients
            for j = 1:npts # Loop over quadrature points
                jac!(J, geom.values, fes.conn[i], gradNparams[j])
                Jac[j] = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
                At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
                gradN!(fes, AllgradN[j], gradNparams[j], csmatTJ);
                dvol = Jac[j]*w[j]
                MeangradN .= MeangradN .+ AllgradN[j]*dvol
                vol = vol + dvol
            end # Loop over quadrature points
            MeangradN .= MeangradN/vol
            Blmat!(self.mr, Bbar, Ns[1], MeangradN, loc, self.mcsys.csmat);
            fill!(elmat,  0.0); # Initialize element matrix
            add_btdb_ut_only!(elmat, Bbar, vol, D, DB)
            add_btdb_ut_only!(elmat, Bbar, -self.phis[i]*vol, Dstab, DB)
            for j = 1:npts # Loop over quadrature points
                Blmat!(self.mr, B, Ns[j], AllgradN[j], loc, self.mcsys.csmat);
                add_btdb_ut_only!(elmat, B, self.phis[i]*Jac[j]*w[j], Dstab, DB)
            end # Loop over quadrature points
            complete_lt!(elmat)
            gatherdofnums!(u, dofnums, fes.conn[i]); # retrieve degrees of freedom
            A_mul_B!(elvec, elmat, elvecfix)
            assemble!(assembler,  -elvec,  dofnums); # assemble element load vector
        end
    end # Loop over elements
    return makevector!(assembler);
end

function nzebcloadsstiffness(self::FEMMDeforLinearAbstractMS, geom::NodalField{FFlt}, u::NodalField{T}) where {T<:Number}
    assembler = SysvecAssembler()
    return  nzebcloadsstiffness(self, assembler, geom, u);
end

function _iip_meanonly(self::FEMMDeforLinearAbstractMS, 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
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, csmatTJ, AllgradN, MeangradN, Jac, D, Dstab, B, DB, Bbar, elmat, elvec, elvecfix = buffers2(self, geom, u, npts)
    MeanN = deepcopy(Ns[1])
    realmat = self.material
    stabmat = self.stabilization_material
    # Sort out  the output requirements
    outputcsys = deepcopy(self.mcsys); # default: report the stresses in the material coord system
    for apair in pairs(context)
        sy, val = apair
        if sy == :outputcsys
            outputcsys = deepcopy(val)
        end
    end
    t= 0.0
    dt = 0.0
    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
    dTe = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    ue = fill(zero(FFlt), size(elmat, 1)); # array of node displacements -- 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!(geom, xe, fes.conn[i]);# retrieve element coords
        gathervalues_asvec!(u, ue, fes.conn[i]);# retrieve element displacements
        gathervalues_asvec!(dT, dTe, fes.conn[i]);# retrieve element temperature increments
        # NOTE: the coordinate system should be evaluated at a single point within the
        # element in order for the derivatives to be consistent at all quadrature points
        loc = centroid!(self,  loc, geom.values, fes.conn[i])
        updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        updatecsmat!(outputcsys, loc, J, fes.label[i]);
        vol = 0.0; # volume of the element
        fill!(MeangradN, 0.0) # mean basis function gradients
        fill!(MeanN, 0.0) # mean basis function gradients
        for j = 1:npts # Loop over quadrature points
            jac!(J, geom.values, fes.conn[i], gradNparams[j])
            Jac[j] = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, AllgradN[j], gradNparams[j], csmatTJ);
            dvol = Jac[j]*w[j]
            MeangradN .= MeangradN .+ AllgradN[j]*dvol
            MeanN .= MeanN .+ Ns[j]*dvol
            vol = vol + dvol
        end # Loop over quadrature points
        MeangradN .= MeangradN/vol
        Blmat!(self.mr, Bbar, MeanN, MeangradN, loc, self.mcsys.csmat);
        MeanN .= MeanN/vol
        qpdT = dot(vec(dTe), vec(MeanN));# Quadrature point temperature increment
        # Quadrature point quantities
        A_mul_B!(qpstrain, Bbar, ue); # strain in material coordinates
        realmat.thermalstrain!(realmat, qpthstrain, qpdT)
        # Material updates the state and returns the output
        out = realmat.update!(realmat, qpstress, out, vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
        if (quantity == :Cauchy)   # Transform stress tensor,  if that is "quantity"
            (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], xe, out, loc);
    end # Loop over elements
    return idat; # return the updated inspector data
end

function _iip_extrapmean(self::FEMMDeforLinearAbstractMS, 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
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, csmatTJ, AllgradN, MeangradN, Jac, D, Dstab, B, DB, Bbar, elmat, elvec, elvecfix = buffers2(self, geom, u, npts)
    MeanN = deepcopy(Ns[1])
    realmat = self.material
    stabmat = self.stabilization_material
    # Sort out  the output requirements
    outputcsys = deepcopy(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
    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
    dTe = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    ue = fill(zero(FFlt), size(elmat, 1)); # array of node displacements -- 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!(geom, xe, fes.conn[i]);# retrieve element coords
        gathervalues_asvec!(u, ue, fes.conn[i]);# retrieve element displacements
        gathervalues_asvec!(dT, dTe, fes.conn[i]);# retrieve element temperature increments
        # NOTE: the coordinate system should be evaluated at a single point within the
        # element in order for the derivatives to be consistent at all quadrature points
        loc = centroid!(self,  loc, geom.values, fes.conn[i])
        updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        updatecsmat!(outputcsys, loc, J, fes.label[i]);
        vol = 0.0; # volume of the element
        fill!(MeangradN, 0.0) # mean basis function gradients
        fill!(MeanN, 0.0) # mean basis function gradients
        for j = 1:npts # Loop over quadrature points
            jac!(J, geom.values, fes.conn[i], gradNparams[j])
            Jac[j] = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, AllgradN[j], gradNparams[j], csmatTJ);
            dvol = Jac[j]*w[j]
            MeangradN .= MeangradN .+ AllgradN[j]*dvol
            MeanN .= MeanN .+ Ns[j]*dvol
            vol = vol + dvol
        end # Loop over quadrature points
        MeangradN .= MeangradN/vol
        Blmat!(self.mr, Bbar, MeanN, MeangradN, loc, self.mcsys.csmat);
        MeanN .= MeanN/vol
        qpdT = dot(vec(dTe), vec(MeanN));# Quadrature point temperature increment
        # Quadrature point quantities
        A_mul_B!(qpstrain, Bbar, ue); # strain in material coordinates
        realmat.thermalstrain!(realmat, qpthstrain, qpdT)
        # Material updates the state and returns the output
        out = realmat.update!(realmat, qpstress, out,
            vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
        if (quantity == :Cauchy)   # Transform stress tensor,  if that is "quantity"
            (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 for each node location
        for nod = 1:size(xe, 1)
            idat = inspector(idat, i, fes.conn[i], xe, out, xe[nod, :]);
        end
    end # Loop over elements
    return idat; # return the updated inspector data
end

function _iip_extraptrend(self::FEMMDeforLinearAbstractMS, 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
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    dofnums, loc, J, csmatTJ, AllgradN, MeangradN, Jac, D, Dstab, B, DB, Bbar, elmat, elvec, elvecfix = buffers2(self, geom, u, npts)
    MeanN = deepcopy(Ns[1])
    realmat = self.material
    stabmat = self.stabilization_material
    # Sort out  the output requirements
    outputcsys = deepcopy(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
    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
    dTe = fill(zero(FFlt), nodesperelem(fes)) # nodal temperatures -- buffer
    ue = fill(zero(FFlt), size(elmat, 1)); # array of node displacements -- 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
    rout1 = fill(zero(FFlt), nstressstrain(self.mr)); # stress -- buffer
    rout =  fill(zero(FFlt), nstressstrain(self.mr));# output -- buffer
    sbout = deepcopy(rout)
    sbout1 = deepcopy(sbout)
    sout = deepcopy(rout)
    sout1 = deepcopy(sout)
    sqploc = deepcopy(loc)
    sstoredout = fill(zero(FFlt), npts, length(sout))
    A = ones(FFlt, npts, 4)
    # 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!(geom, xe, fes.conn[i]);# retrieve element coords
        gathervalues_asvec!(u, ue, fes.conn[i]);# retrieve element displacements
        gathervalues_asvec!(dT, dTe, fes.conn[i]);# retrieve element temperature increments
        # NOTE: the coordinate system should be evaluated at a single point within the
        # element in order for the derivatives to be consistent at all quadrature points
        loc = centroid!(self,  loc, geom.values, fes.conn[i])
        updatecsmat!(self.mcsys, loc, J, fes.label[i]);
        updatecsmat!(outputcsys, loc, J, fes.label[i]);
        vol = 0.0; # volume of the element
        fill!(MeangradN, 0.0) # mean basis function gradients
        fill!(MeanN, 0.0) # mean basis function gradients
        for j = 1:npts # Loop over quadrature points
            jac!(J, geom.values, fes.conn[i], gradNparams[j])
            Jac[j] = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            At_mul_B!(csmatTJ, self.mcsys.csmat, J); # local Jacobian matrix
            gradN!(fes, AllgradN[j], gradNparams[j], csmatTJ);
            dvol = Jac[j]*w[j]
            MeangradN .= MeangradN .+ AllgradN[j]*dvol
            MeanN .= MeanN .+ Ns[j]*dvol
            vol = vol + dvol
        end # Loop over quadrature points
        MeangradN .= MeangradN/vol
        Blmat!(self.mr, Bbar, MeanN, MeangradN, loc, self.mcsys.csmat);
        MeanN .= MeanN/vol
        qpdT = dot(vec(dTe), vec(MeanN));# Quadrature point temperature increment
        # Quadrature point quantities
        A_mul_B!(qpstrain, Bbar, ue); # strain in material coordinates
        realmat.thermalstrain!(realmat, qpthstrain, qpdT)
        # REAL Material updates the state and returns the output
        rout = realmat.update!(realmat, qpstress, rout,
            vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
        if (quantity == :Cauchy)   # Transform stress tensor,  if that is "quantity"
            (length(rout1) >= length(rout)) || (rout1 = zeros(length(rout)))
            rotstressvec(self.mr, rout1, rout, transpose(self.mcsys.csmat))# To global coord sys
            rotstressvec(self.mr, rout, rout1, outputcsys.csmat)# To output coord sys
        end
        # STABILIZATION Material updates the state and returns the output
        sbout = stabmat.update!(stabmat, qpstress, sbout,
            vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
        if (quantity == :Cauchy)   # Transform stress tensor,  if that is "quantity"
            (length(sbout1) >= length(sbout)) || (sbout1 = zeros(length(sbout)))
            rotstressvec(self.mr, sbout1, sbout, transpose(self.mcsys.csmat))# To global coord sys
            rotstressvec(self.mr, sbout, sbout1, outputcsys.csmat)# To output coord sys
        end
        for j = 1:npts # Loop over quadrature points (STABILIZATION material)
            At_mul_B!(sqploc, Ns[j], xe);# Quadrature point location
            A[j, 1:3] .= vec(sqploc - loc);
            Blmat!(self.mr, B, Ns[j], AllgradN[j], sqploc, self.mcsys.csmat);
            qpdT = dot(vec(dTe), vec(Ns[j]));# Quadrature point temperature increment
            #  Quadrature point quantities
            A_mul_B!(qpstrain, B, ue); # strain in material coordinates
            stabmat.thermalstrain!(stabmat, qpthstrain, qpdT)
            # Material updates the state and returns the output
            sout = stabmat.update!(stabmat, qpstress, sout,
                vec(qpstrain), qpthstrain, t, dt, loc, fes.label[i], quantity)
            if (quantity == :Cauchy)   # Transform stress tensor,  if that is "quantity"
                (length(sout1) >= length(sout)) || (sout1 = zeros(length(sout)))
                rotstressvec(self.mr, sout1, sout, transpose(self.mcsys.csmat))# To global coord sys
                rotstressvec(self.mr, sout, sout1, outputcsys.csmat)# To output coord sys
            end
            sstoredout[j, :] .= sout # store  the output for this quadrature point
        end # Loop over quadrature points
        #  Solve for the least-square fit parameters
        fact = qr(A);
        p = UpperTriangular(fact.R)\(transpose(Array(fact.Q))*sstoredout) # R \ (transpose(Q) * sstoredout)
        for nod = 1:size(xe, 1)
            #  Predict the value  of the output quantity at the node
            xdel = vec(@view xe[nod, :]) - vec(loc)
            nout = rout + self.phis[i] * (- sbout + vec(reshape(xdel, 1, 3) * p[1:3, :]) + p[4, :])
            # Call the inspector for the node location
            idat = inspector(idat, i, fes.conn[i], xe, nout, xe[nod, :]);
        end
    end # Loop over elements
    return idat; # return the updated inspector data
end

"""
   inspectintegpoints(self::FEMMDeforLinear,
     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` - functionwith 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::FEMMDeforLinearAbstractMS, geom::NodalField{FFlt},  u::NodalField{T}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F,  idat, quantity=:Cauchy; context...) where {T<:Number, F<:Function}
    reportat = :meanonly
    for apair in pairs(context)
        sy, val = apair
        if sy == :reportat
            reportat = val
        end
    end
    if reportat == :extraptrend
        return _iip_extraptrend(self, geom, u, dT, felist,
            inspector, idat, quantity; context...);
    elseif reportat == :extrapmean
        return _iip_extrapmean(self, geom, u, dT, felist,
            inspector, idat, quantity; context...);
    elseif reportat == :meanonly || true # this is the default
        return _iip_meanonly(self, geom, u, dT, felist,
            inspector, idat, quantity; context...);
    end
    return idat
end

function inspectintegpoints(self::FEMMDeforLinearAbstractMS, 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


text:
using BenchmarkTools
using Folds
using FoldsCatwalk
using Transducers

type_instability(x) = Val(round(Int, 2x))
asint(::Val{x}) where {x} = x::Int

function sum_baseline(xs)
    itr = xs |> Map(type_instability) |> Map(asint)
    Folds.sum(itr, SequentialEx())
end

function sum_catwalk(xs)
    itr = xs |> Map(type_instability) |> OptimizeInner() |> Map(asint)
    Folds.sum(itr, CatwalkEx())
end

function demo_sum()
    xs = randn(1000_000)
    @assert sum_baseline(xs) == sum_catwalk(xs)
    @btime sum_baseline($xs)
    @btime sum_catwalk($xs)
end

ext: jl


text:
export soundspeed, absorption, density, reflectioncoef, surfaceloss, doppler, bubbleresonance

function soundspeed(temperature=27.0, salinity=35.0, depth=10.0; voidfrac=0.0, cgas=340.0, reldensity=1000.0)
  # based on Mackenzie (1981)
  c = 1448.96 + 4.591*temperature - 5.304e-2*temperature^2 + 2.374e-4*temperature^3
  c += 1.340*(salinity-35) + 1.630e-2*depth + 1.675e-7*depth^2
  c += -1.025e-2*temperature*(salinity-35) - 7.139e-13*temperature*depth^3
  if voidfrac > 0.0
    # based on Wood (1964) or Buckingham (1997)
    m = √reldensity
    c = 1.0/(1.0/c*√((voidfrac*(c/cgas)^2*m+(1.0-voidfrac)/m)*(voidfrac/m+(1-voidfrac)*m)))
  end
  return c
end

function absorption(frequency, distance=1000.0, temperature=27.0, salinity=35.0, depth=10.0, pH=8.1)
  # based on Francois-Garrison model
  f = frequency/1000.0
  d = distance/1000.0
  c = 1412.0 + 3.21*temperature + 1.19*salinity + 0.0167*depth
  A1 = 8.86/c * 10.0^(0.78*pH-5.0)
  P1 = 1.0
  f1 = 2.8*_np.sqrt(salinity/35.0) * 10.0^(4.0-1245.0/(temperature+273.0))
  A2 = 21.44*salinity/c*(1.0+0.025*temperature)
  P2 = 1.0 - 1.37e-4*depth + 6.2e-9*depth*depth
  f2 = 8.17 * 10^(8.0-1990.0/(temperature+273.0)) / (1.0+0.0018*(salinity-35.0))
  P3 = 1.0 - 3.83e-5*depth + 4.9e-10*depth*depth
  if temperature < 20.0
    A3 = 4.937e-4 - 2.59e-5*temperature + 9.11e-7*temperature^2 - 1.5e-8*temperature^3
  else
    A3 = 3.964e-4 - 1.146e-5*temperature + 1.45e-7*temperature^2 - 6.5e-10*temperature^3
  end
  a = A1*P1*f1*f*f/(f1*f1+f*f) + A2*P2*f2*f*f/(f2*f2+f*f) + A3*P3*f*f
  10^(-a*d/20.0)
end

function density(temperature=27, salinity=35)
  # based on Fofonoff (1985 - IES 80)
  t = temperature
  A = 1.001685e-04 + t * (-1.120083e-06 + t * 6.536332e-09)
  A = 999.842594 + t * (6.793952e-02 + t * (-9.095290e-03 + t * A))
  B = 7.6438e-05 + t * (-8.2467e-07 + t * 5.3875e-09)
  B = 0.824493 + t * (-4.0899e-03 + t * B)
  C = -5.72466e-03 + t * (1.0227e-04 - t * 1.6546e-06)
  D = 4.8314e-04
  A + salinity * (B + C*_np.sqrt(salinity) + D*salinity)
end

function reflectioncoef(θ, ρ1, c1, α, ρ=density(), c=soundspeed())
  # based on Brekhovskikh & Lysanov
  n = c/c1*(1.0+1.0im*α)
  m = ρ1/ρ
  t1 = m*cos(θ)
  t2 = √(n^2 - sin(θ)^2)
  V = (t1-t2)/(t1+t2)
  imag(V) ≈ 0.0 ? real(V) : V
end

function surfaceloss(windspeed, frequency, θ)
  # based on APL model (1994)
  β = π/2 - θ
  if windspeed >= 6.0
    a = 1.26e-3/sin(β) * windspeed^1.57 * frequency^0.85
  else
    a = 1.26e-3/sin(β) * 6^1.57 * frequency^0.85 * exp(1.2*(windspeed-6.0))
  end
  db2amp(-a)
end

doppler(speed, frequency, c=soundspeed()) = (1.0+speed/c)*frequency
bubbleresonance(radius, depth=0) = 3.25/radius * √(1+0.1*depth)     # Medwin & Clay (1998)


text:
using EquilibriumMeasures, StaticArrays, ClassicalOrthogonalPolynomials, FillArrays, BlockArrays, LazyBandedMatrices, LinearAlgebra, IntervalSets, Test
import EquilibriumMeasures: EquilibriumMeasureMoment
using ForwardDiff: jacobian
using DomainSets: components

@testset "one-interval" begin
    μ = equilibriummeasure(x -> x^2)
    @test first(axes(μ,1)) ≈ -sqrt(2) atol=1E-13
    @test last(axes(μ,1)) ≈ sqrt(2) atol=1E-13
    @test sum(μ) ≈ 1 atol=1E-13
    @test μ[0.1] ≈ sqrt(2-0.1^2)/π atol=1E-13

    μ = equilibriummeasure(x -> x^4)
    @test first(axes(μ,1)) ≈ -1.0745699318235422 atol=1E-13
    @test last(axes(μ,1)) ≈ 1.0745699318235422 atol=1E-13
    @test sum(μ) ≈ 1 atol=1E-13
    @test μ[0.1] ≈ 0.40005825716887744

    μ = equilibriummeasure(x -> x^4 + x^3 - x^2)
    @test sum(μ) ≈ 1 atol=1E-13
    @test μ[0.1] ≈ 0.19674192408704289

    μ = equilibriummeasure(x -> (x-3)*(x-2)*(1+x)*(2+x)*(3+x)*(2x-1)/250; a=SVector(-3,3))
    @test sum(μ) ≈ 1 atol=1E-13

    # non-unique minimiser
    μ, b = equilibriummeasure(x -> (x-3)*(x-2)*(1+x)*(2+x)*(3+x)*(2x-1)/20; a=SVector(-3,-2), returnendpoint=true)
    @test sum(μ) ≈ 1 atol=1E-13

    # Deflate found solution and find another one whilst damping the Newton step
    μ = equilibriummeasure(x -> (x-3)*(x-2)*(1+x)*(2+x)*(3+x)*(2x-1)/20; a=SVector(-3,-2), knownsolutions=[b], dampening=0.3)
    @test sum(μ) ≈ 1 atol=1E-6
end

@testset "two-interval" begin
    # see EquilibriumMeasureExamples.nb
    V = x -> 0.7*(x^4 - 2x^3 - x^2 + 2x)
    a_ex = SVector(-1.0637226766068189, 0.2659671162729329, 0.7340328837270674, 2.063722676606819)
    @test norm(EquilibriumMeasureMoment(V)(a_ex)) ≤ 1E-13

    @time μ = equilibriummeasure(V; a=SVector(-1,0.25,0.7,2))
    @test all(components(axes(μ,1).domain) .≈ (a_ex[1]..a_ex[2], a_ex[3]..a_ex[4]))

    V = x -> (x-3)*(x-2)*(1+x)*(2+x)*(3+x)*(2x-1)/150
    a_ex = SVector(-3.0782937003244046, 0.9082915776561011, 1.8384672185745663, 3.0223539638196737)
    @test norm(EquilibriumMeasureMoment(V)(a_ex)) ≤ 1E-13

    @time μ = equilibriummeasure(V; a=SVector(-3,1,2,3))
    @test vcat(leftendpoint.(components(axes(μ,1).domain))...) ≈ a_ex[[1,3]]
    @test vcat(rightendpoint.(components(axes(μ,1).domain))...) ≈ a_ex[[2,4]]
end

max_stars_repo_name: ApproxFun/EquilbriumMeasures

max_issues_repo_name: ApproxFun/EquilbriumMeasures

max_forks_repo_name: ApproxFun/EquilbriumMeasures


text:
@testset "PTBLM" begin
    x, y = PTBLM.traindata()
    x, y = PTBLM.testdata()

end

@testset "UD_English" begin
    x = UD_English.traindata()
    x = UD_English.devdata()
    x = UD_English.testdata()
end

@testset "SMS Spam Collection" begin
    X = SMSSpamCollection.features()
    y = SMSSpamCollection.targets()
    @test X isa Vector
    @test y isa Vector
    @test size(X) == (5574,)
    @test size(y) == (5574,)
end


text:
struct PosteriorMean <: Estimator
    prior::Prior
end

string(estimator::PosteriorMean) = @sprintf "PosteriorMean<%s>" string(estimator.prior)

function (estimator::PosteriorMean)(x1::TI, x2::TI, design::TD) where {TI<:Integer,TD<:AbstractDesign}
    mean(update(estimator.prior, x1 + x2, n(design, x1)))
end

# TODO: implement this as cached() / CachedEstimator for all!
struct PosteriorMeanPrecalculated{TI<:Integer,TR<:Real,TD<:AbstractDesign,TP<:Prior} <: Estimator
    design::TD
    prior::TP
    XX::Array{TI,2} # sample space
    estimates::Array{TR,1}
end

string(estimator::PosteriorMeanPrecalculated) = @sprintf "PosteriorMean(pre-calculated)<%s;%s>" string(estimator.prior) string(estimator.design)

function (estimator::PosteriorMeanPrecalculated)(x1::TI, x2::TI, design::TD) where {TI<:Integer,TD<:AbstractDesign}

    for i in 1:size(estimator.XX, 1)
        if [x1, x2] == estimator.XX[i,:]
            return estimator.estimates[i]
        end
    end
    error("(x1,x2) not found in sample space, valid observation?")
end

function PosteriorMeanPrecalculated(prior::Prior, design::AbstractDesign)

    XX        = sample_space(design)
    estimates = zeros(size(XX, 1))
    for i in 1:size(XX, 1)
        x1, x2 = XX[i,:]
        estimates[i] = mean(update(prior, x1 + x2, n(design, x1)))
    end
    return PosteriorMeanPrecalculated{eltype(XX),eltype(estimates),typeof(design),typeof(prior)}(
        design, prior, XX, estimates
    )
end


text:
module test_julia_ccall

using Test

# https://github.com/JuliaLang/julia/blob/master/base/essentials.jl#L339

end # module test_julia_ccall


text:
using MPI, Test

include("../testhelpers.jl")

@testset "Driver" begin
    runmpi(joinpath(@__DIR__, "cr_unit_tests.jl"), ntasks = 1)
    runmpi(joinpath(@__DIR__, "les_driver_test.jl"), ntasks = 1)
    runmpi(joinpath(@__DIR__, "mms3.jl"), ntasks = 2)
end


text:
# Perfect power function

@doc raw"""
```julia
isperfectpower(n::Integer) -> Bool
```

Given an integer `n`, returns `true` or `false` depending on whether or not it is a perfect power.

Here, a perfect power is some number of the form ``a^b``, where ``a, b \in \mathbb{N}``, and ``b > 1``.

This function is a wrapper around [Hecke.jl's really excellent and efficient `ispower` function](https://github.com/thofma/Hecke.jl/blob/master/src/Misc/Integer.jl#L413-L443).

!!! note

	It should be noted that there is another defintion for "perfect powers":
		> some number of the form ``p^b``, where ``p, b \in \mathbb{N}``, *and ``p`` is a prime number``*.
	Here, we call numbers of this form *prime powers*, or *proper perfect powers*.  By this definition, 36 is *not* a perfect power, because 6 is not prime.

See also: `CodingTheory.isperfectpower`.
"""
function isperfectpower(n::Integer)
    (n == 0) && return false
    e = first(ispower(n))
    return (e != 1 && e != 0)
end

@doc raw"""
```julia
isprimepower(n::Integer) -> Bool
```

Given an integer `n`, returns `true` or `false` depending on whether or not it is a prime power.

A prime power is some number of the form ``p^b``, where ``p, b \in \mathbb{N}``, and ``p`` is a prime number``.

This function is a wrapper around [Hecke.jl's really excellent and effcient `ispower` function](https://github.com/thofma/Hecke.jl/blob/master/src/Misc/Integer.jl#L756-L769).

See also: `CodingTheory.isperfectpower`.
"""
function isprimepower(n::Integer)
    return isprime_power(n)
end

## Prime power function

#=
function __isprime(n)
	return isprime(fmpz(n))
end

function __factors(n)
	return n == 0 ? [] : factor(fmpz(n))
end
# function __factors(n)
# 	first(i) for i in Primes.factor(n)
# end

function primedivisors(n)
    n == 0 && return []
    __isprime(n) && return [fmpz(n)]
    f = __factors(n)
    return sort!([p for (p, e) ∈ f])
end

function ω(n)
	nprimedivisors = 0
	if n == 0
		nprimedivisors = 0
	elseif __isprime(n)
		nprimedivisors = 1
    else
        nprimedivisors = length(__factors(n))
	end

	return fmpz(nprimedivisors)

    # # original definition:
	# return fmpz(length(primedivisors(n)))
end

function Ω(n)
    n == fmpz(0) && return 0
    __isprime(n) && return fmpz(1)

    return sum(e for (__, e) ∈ __factors(n))
end

# this is a slow version compared to the Hecke one
function isperfectpower(n)
	return ω(n) == 1 && Ω(n) != 1
end

# this is slower than the Hecke verson
function isprimepower(n)
    return ω(n) == 1
end

=#


text:
include("algorithm.jl")

graphNodes = 3
#              A  B   C   D    E
adjMatrix = [ -1  500 400 500 800;
              500 -1  200 200 300;
              400 200 -1  300 500;
              500 200 300 -1  200;
              800 300 500 200 -1  ]
foodNode = 3
maxIteration = 200
antsCounts = 20
pheromoneExponent = 1
desirabilityExponent = 0.1
pheromoneDepositFactor = 1
evaporationRate = 0.1

graph = generateGraph(graphNodes,adjMatrix,foodNode)

bestPath, bestFitness, graph = algorithmAntColonyOptimization(graph, maxIteration, antsCounts, pheromoneExponent,
                                            desirabilityExponent, pheromoneDepositFactor, evaporationRate)

print("Best solution $(bestPath) with length $(bestFitness)")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
*.jl.cov
*.jl.*.cov
*.jl.mem

Manifest.toml
docs/build/
docs/site/


text:
# beehive larger test

using Test
using RoME
using Distributed

while nprocs() < 3
  addprocs(1)
end

using RoME
@everywhere using RoME

##

@testset "sequential beehive, repeated growing and solving" begin

##

fg = RoME.generateCanonicalFG_Honeycomb!(7, graphinit=true, useMsgLikelihoods = true)
tree = solveTree!(fg);

fg = RoME.generateCanonicalFG_Honeycomb!(14, graphinit=true, dfg=fg)
tree = solveTree!(fg, tree);

tree_ = deepcopy(tree); fg_ = deepcopy(fg);
fg = RoME.generateCanonicalFG_Honeycomb!(21, graphinit=true, dfg=fg)
tree = solveTree!(fg  , tree);

# fg = RoME.generateCanonicalFG_Honeycomb!(28, graphinit=true, dfg=fg)
# tree = solveTree!(fg, tree);

# fg = RoME.generateCanonicalFG_Honeycomb!(35, graphinit=true, dfg=fg)
# tree = solveTree!(fg, tree);

# fg = RoME.generateCanonicalFG_Honeycomb!(42, graphinit=true, dfg=fg)
# tree = solveTree!(fg, tree);

## test numerical results

# working before IIF 1010

@warn("Test for beehive graph is using loose bounds until IIF #1010 is resolved.")
@test isapprox( getPPE(fg, :l11).suggested , [5;10*sin(pi/3)], atol=6)
@test isapprox( getPPE(fg, :l0).suggested , [20;0], atol=4); @error "degraded numerical performance"
@test isapprox( getPPE(fg, :l7).suggested , [20;-20*sin(pi/3)], atol=6)

# likely to fail until IIF 1010 is completed
@test_skip  isapprox( getPPE(fg, :x21).suggested[1:2] , [10;-20*sin(pi/3)], atol=4)

##

end

@testset "Test stochastic Beehive trajectory generator" begin

##

fg = generateCanonicalFG_Beehive!(8, graphinit=false, postpose_cb=(g,l)->@show l)

@test isapprox( getPPE(fg, :x0, :simulated).suggested[1:2], [0.0;0.0], atol = 1e-8)

##

end

##

# using RoMEPlotting
# Gadfly.set_default_plot_size(40cm,25cm)

# ##

# plotSLAM2D(fg, drawContour=false, drawPoints=false, drawEllipse=false)

#


text:
module EmbeddingsAnalysis

    using LinearAlgebra
    using Statistics
    using Languages
    using Word2Vec
    using ConceptnetNumberbatch
    using StatsBase
    using MultivariateStats
    using Distances
    using QuantizedArrays

    import Base: size
    import Word2Vec: analogy_words,
                     WordVectors,
                     get_vector,
                     vocabulary

    export WordVectors,
           get_vector,
           cosine_vec,
           conceptnet2wv,
           CompressedWordVectors,
           compressedwordvectors,
           compress,
           analogy_words,
           vocabulary,
           write2disk,
           similarity_order,
           pca_reduction,
           vocab_reduction,
           translate

    include("defaults.jl")          # defaults
    include("utils.jl")             # useful utilities
    include("conceptnet2wv.jl")     # convert ConceptNet to WordVectors
    include("cwv.jl")               # CompressedWordVectors
    include("write2disk.jl")        # save WordVectors to disk
    include("similarity_order.jl")  # preprocess WordVectors
    include("pca_reduction.jl")     # preprocess/reduce dimensionality of WordVectors
    include("vocab_reduction.jl")   # reduce vocabulary size of WordVectors using a seed vocabulary
    include("translator.jl")        # word-vector based translation

end # module

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
import Base.iterate

function iterate(a::RingList)
    if length(a)==0
        return nothing
    end
    f = first(a)
    return f,f
end

function iterate(a::RingList, st)
    y = a[st]
    if y==first(a)
        return nothing
    end
    return y,y
end


text:
function __msc0(D::GrB_Matrix{T}, cols) where T
    r = GrB_Matrix_nrows(D)
    Ires, Xres = ZeroBasedIndex[], T[]
    vec = SparseMM.gbv_new(T, r)

    for j in cols
        GrB_Col_extract(vec, GrB_NULL, GrB_NULL, D, GrB_ALL, 0, ZeroBasedIndex(j.x-1), GrB_NULL)
        I, X = GrB_Vector_extractTuples(vec)
        append!(Ires, I)
        append!(Xres, X)
        GrB_Vector_clear(vec)
    end

    GrB_free(vec)
    return Ires, Xres

end

"""
    msc(D::GrB_Matrix{T}, v)

Sums D's columns which indexes are in V

# Arguments
 - `D`: a GraphBLAS matrix
 - `v`: list of indices
"""
function msc(D::GrB_Matrix{T}, v) where T
    r = GrB_Matrix_nrows(D)
    res = SparseMM.gbm_new(T, r, length(v))

    X = T[]
    J = ZeroBasedIndex[]
    I = ZeroBasedIndex[]

    for (i, c) in enumerate(v)
        _I, _X = __msc0(D, c)

        append!(I, _I)
        append!(X, _X)
        append!(J, fill(i-1, length(_X)))
    end

    GrB_Matrix_build(res, I, J, X, length(X), SparseMM.GrB_op("PLUS", T))

    return res

end

"""
    cellCongruence(Delta::GrB_Matrix{T}, vclasses)

Performs the Cell Congruence Enabling algorithm.

# Arguments
 - `Delta`: the Cochain operator
 - `classes`: list that indentifies old cells
"""
function cellCongruence(Delta::GrB_Matrix{T}, classes) where T
    copEV = msc(Delta, classes)

    function __congruence(copEV::GrB_Matrix{T}) where T
        ec = Dict{ZeroBasedIndex, Array{ZeroBasedIndex, 1}}()
        foreach(zip(GrB_Matrix_extractTuples(copEV)...)) do (i, j, _)
            push!(get!(ec, i, []), j)
        end
        c = Dict{Array{ZeroBasedIndex, 1}, Array{ZeroBasedIndex, 1}}()
        for (k, v) in ec
            push!(get!(c, v, []), k)
        end
        eclasses = [sort(v) for (k, v) in c]
        sort!(eclasses)
        return eclasses
    end
    eclasses = __congruence(copEV)
    newedges = first.(eclasses)

    res = SparseMM.gbm_new(T, GrB_Matrix_ncols(copEV), length(newedges))

    I, J, X = GrB_Matrix_extractTuples(copEV)
    _I, _J, _X = ZeroBasedIndex[], ZeroBasedIndex[], T[]

    for (i, j, x) in zip(I, J, X)
        if i in newedges
            push!(_J, findfirst(x -> x==i ,newedges)-1)
            push!(_I, j)
            push!(_X, x)
        end
    end

    GrB_Matrix_build(res, _I, _J, _X, length(_X), SparseMM.GrB_op("FIRST", T))
    return res, eclasses
end

"""
    chainCongruenceGB(G, Top)

Perform the Geometry ``G`` congruence and coherently reshape the Topology ``T`` (GraphBLAS implementation)
"""
function chainCongruenceGB(W, Top)
	V, cls = vertCongruence(W)
	cls = [[ZeroBasedIndex(c) for c in cl] for cl in cls]

	Topn = Array{GrB_Matrix{Int8}}(undef, length(Top))
	for i in eachindex(Topn)
		Topn[i], cls = cellCongruence(Top[i], cls)
	end

	return V, Topn
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ContainerServiceSshPublicKey <: SwaggerModel
    keyData::Any # spec type: Union{ Nothing, String } # spec name: keyData

    function ContainerServiceSshPublicKey(;keyData=nothing)
        o = new()
        validate_property(ContainerServiceSshPublicKey, Symbol("keyData"), keyData)
        setfield!(o, Symbol("keyData"), keyData)
        o
    end
end # type ContainerServiceSshPublicKey

const _property_map_ContainerServiceSshPublicKey = Dict{Symbol,Symbol}(Symbol("keyData")=>Symbol("keyData"))
const _property_types_ContainerServiceSshPublicKey = Dict{Symbol,String}(Symbol("keyData")=>"String")
Base.propertynames(::Type{ ContainerServiceSshPublicKey }) = collect(keys(_property_map_ContainerServiceSshPublicKey))
Swagger.property_type(::Type{ ContainerServiceSshPublicKey }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ContainerServiceSshPublicKey[name]))}
Swagger.field_name(::Type{ ContainerServiceSshPublicKey }, property_name::Symbol) =  _property_map_ContainerServiceSshPublicKey[property_name]

function check_required(o::ContainerServiceSshPublicKey)
    (getproperty(o, Symbol("keyData")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ContainerServiceSshPublicKey }, name::Symbol, val)
end


text:
struct Stream{T_<:Real}
    subst_::Substance
    ṁ::T_
    V̇::T_
    p::T_
    T::T_
    ρ::T_
    Ĥ::T_
    H::T_
    Cp::T_
    μ::T_
    k::T_
    ν::T_
    α::T_
    Pr::T_
    β::T_
end

function Stream(subst_::Substance, ṁ, p, T)
    ρ = subst_.density(p, T)
    V̇ = ṁ / ρ
    Ĥ = subst_.enthalpy(p, T)
    H = Ĥ * ṁ
    Cp = subst_.heat_capacity(p, T)
    μ = subst_.viscosity(p, T)
    k = subst_.thermal_conductivity(p, T)
    ν = subst_.visc_kin(p, T)
    α = subst_.therm_diff(p, T)
    Pr = subst_.Pr(p, T)
    β = subst_.exp_coeff(p, T)
    return Stream(subst_, ṁ, V̇, p, T, ρ, Ĥ, H, Cp, μ, k, ν, α, Pr, β)
end

# function Stream(subst_::String, ṁ, p, T)
#     return Stream(Stream_substances[subst_], ṁ, p, T)
# end

struct FlowStream{T<:Real}
    tpipe::Pipe
    tstream::Stream
    Re::T
    Pr::T
    f::T
    ΔP::T
    Δh::T
end

"""
FlowStream: Basically takes in a Stream and calculates relevant parameters for
its flow through a pipe, including the Reynolds and Prandtl numbers, the
friction factor, and pressure/head losses. Assumes negligible change in fluid
properties (notably density and viscosity) down the length of the pipe.

Uses well-known equations that can be found in (or derived from) any good book
on fluid mechanics.
"""
function FlowStream(tpipe::Pipe, tStream::Stream)
    d = tpipe.diameter
    Re = tStream.ṁ * d / (tStream.μ * tpipe.flowarea)
    f = getf2(Re, tpipe)
    vel = tStream.V̇ / tpipe.flowarea
    Δh = (f*tpipe.length / d) * vel^2 / (2*9.81) + tpipe.elevation_change
    ΔP = Δh * (9.81 * tStream.ρ)
    return FlowStream(tpipe, tStream, Re, tStream.Pr, f, ΔP, Δh)
end


text:
# fitting

"""
    struct FittedPlane{A<:AbstractVector,B<:AbstractVector} <: FittedShape

Plane primitive, defined by one of its point, and its normalvector.
"""
struct FittedPlane{A<:AbstractVector,B<:AbstractVector} <: FittedShape
    point::A
    normal::B
end

Base.show(io::IO, x::FittedPlane) =
    print(io, """plane""")

Base.show(io::IO, ::MIME"text/plain", x::FittedPlane{A,B}) where {A,B} =
    print(io, """FittedPlane{$A,$B}\nnormal: $(x.normal), point: $(x.point) """)

strt(x::FittedPlane) = "plane"

function defaultshapeparameters(::Type{FittedPlane})
    return (plane=(ϵ=0.3, α=deg2rad(5)),)
end

"""
    fit(::Type{FittedPlane}, p, n, pc, params)

Fit a plane to 3 points. Their and additional point's normals are used to validate the fit.
Return `nothing` if points do not fit to a plane.

A collinearity check is used to not filter out points on one line.
"""
function fit(::Type{FittedPlane}, p, n, pc, params)
    #@unpack α_plane, collin_threshold = params
    @extract params : params_plane=plane
    @extract params_plane : α_plane=α
    @extract params.common : collin_threshold
    lp = length(p)
    @assert length(p) > 2 "At least 3 point is needed."
    @assert lp == length(n) "Size must be the same."
    crossv = normalize(cross(p[2]-p[1], p[3]-p[1]))
    # how many point's normal must be checked
    norm(crossv) < collin_threshold && return nothing
    # here we have the normal of the theoretical plane
    norm_ok = falses(lp)
    invnorm_ok = falses(lp)

    thr = cos(α_plane)
    for i in 1:lp
        dotp = dot(crossv, normalize(n[i]))
        norm_ok[i] = dotp > thr
        invnorm_ok[i] = dotp < -thr
    end
    norm_ok == trues(lp) && return FittedPlane(p[1], crossv)
    invnorm_ok == trues(lp) && return FittedPlane(p[1], -1*crossv)
    return nothing
end

# bitmapping

function scorecandidate(pc, candidate::FittedPlane, subsetID, params)
    ps = @view pc.vertices[pc.subsets[subsetID]]
    ns = @view pc.normals[pc.subsets[subsetID]]
    ens = @view pc.isenabled[pc.subsets[subsetID]]

    cp = compatiblesPlane(candidate, ps, ns, params)
    inder = cp.&ens
    inpoints = (pc.subsets[subsetID])[inder]
    score = estimatescore(length(pc.subsets[subsetID]), pc.size, length(inpoints))
    return (score, inpoints)
end

"""
    project2plane(plane, points)

Project `points` on to the `plane`.

# Arguments:
- 'plane::FittedPlane': a plane.
- 'points::AbstractArray': an array of points (array-of-arrays).
"""
function project2plane(plane, points)
    # create a coordinate frame
    # z is the plane's normal
    o_z = normalize(plane.normal)
    # x is a random orthogonal vector (in the plane)
    o_x = normalize(arbitrary_orthogonal(o_z))
    # y is created so it's a right hand coord. frame
    o_y = normalize(cross(o_z, o_x))
    #answer = similar(points)
    # get the coordinates of the points in the prev. created coord. frame
    function proj_plane(p)
        v = p-plane.point
        return eltype(points)(dot(o_x,v), dot(o_y,v), dot(o_z,v))
    end

    #for i in eachindex(points)
    #    v = points[i]-plane.point
    #    answer[i] = eltype(answer)(dot(o_x,v), dot(o_y,v), dot(o_z,v))
    #end
    #answer
    return (proj_plane(ps) for ps in points)
end

"""
    compatiblesPlane(plane, points, normals, eps, alpharad)

Create a bool-indexer array for those points that are compatible to the plane.
Give back the projected points too for parameter space magic.

Compatibility is measured with an `eps` distance to the plane
and an `alpharad` angle to it's normal.
"""
function compatiblesPlane(plane, points, normals, params)
    #@unpack ϵ_plane, α_plane = params
    @extract params : params_plane=plane
    @extract params_plane : α_plane=α ϵ_plane=ϵ
    @assert length(points) == length(normals) "Size must be the same."
    projecteds = project2plane(plane, points)
    # eps check
    #c1 = [abs(a[3]) < ϵ_plane for a in projecteds]
    # alpha check
    #c2 = [isparallel(plane.normal, normals[i], α_plane) && c1[i] for i in eachindex(normals)]
    # projecteds[comp] are the compatible points

    # zip projected points and normals
    zpn = zip(projecteds, normals)
    comp = [isparallel(plane.normal, n, α_plane) && (abs(p[3]) < ϵ_plane) for (p,n) in zpn]
    return comp
end

"""
    refit(s::T, pc, params) where {T<:FittedPlane}

Refit plane.
"""
function refit(s::T, pc, params) where {T<:FittedPlane}
    # TODO: use octree for that
    pcv = @view pc.vertices[pc.isenabled]
    pcn = @view pc.normals[pc.isenabled]
    cp = compatiblesPlane(s, pcv, pcn, params)
    return ExtractedShape(s, ((1:pc.size)[pc.isenabled])[cp])
end


text:
module MultiWayAggregation

@assert VERSION>=v"1.4"

using DataFrames
##using CSV
using IterTools

export multiwayaggregation
function multiwayaggregation(df::DataFrame,v::Symbol,cs::Union{Pair, typeof(nrow), DataFrames.ColumnIndex, DataFrames.MultiColumnIndex}...)
    res = multiwayaggregation(df,vcat(v),cs...)
    return res
end

function multiwayaggregation(df::DataFrame,v::Vector{Symbol},cs::Union{Pair, typeof(nrow), DataFrames.ColumnIndex, DataFrames.MultiColumnIndex}...)
    res=DataFrame()

    for c in v
        @assert !(any(ismissing,df[!,c])) "Column $(c) has missing values. This is currently not supported by multiwayaggregation" #otherwise the appending will not be meaningful, as we set the values to missing for columns which are not considered in the multi way summary
    end

    for subsetlength=length(v):-1:0
        for subs in IterTools.subsets(v,subsetlength)
            #@show subs
            if subsetlength == 0
                agg = DataFrames.combine(df,cs...)
            else
                agg = DataFrames.combine(DataFrames.groupby(df,subs),cs...)
            end
            nonAggregatedVars=setdiff(v,subs)

            k = 1
            DataFrames.insertcols!(agg,k,:_TYPE_ => repeat(vcat(subsetlength),size(agg,1)))
            k += 1
            for addcol in nonAggregatedVars
                DataFrames.insertcols!(agg,k,addcol => repeat(vcat(missing),size(agg,1)))
                k += 1
            end

            DataFrames.allowmissing!(agg)
            DataFrames.append!(res,agg)
        end
    end

    sort!(res,vcat(:_TYPE_,v))
    return res
end

#multiwayaggregationkw
#same as multiwayaggregation but instead of 'missing' a keyword (e.g. "ALL") is used to denote subtotal summaries
export multiwayaggregationkw
function multiwayaggregationkw(df::DataFrame,v::Symbol,subtotalkw,cs::Union{Pair, typeof(nrow), DataFrames.ColumnIndex, DataFrames.MultiColumnIndex}...)
    res = multiwayaggregationkw(df,vcat(v),subtotalkw,cs...)
    return res
end

function multiwayaggregationkw(df::DataFrame,v::Vector{Symbol},subtotalkw,cs::Union{Pair, typeof(nrow), DataFrames.ColumnIndex, DataFrames.MultiColumnIndex}...)
    res = DataFrame()

    for c in v
        if eltype(df[!,c]) != typeof(subtotalkw)
            @show c
            @show eltype(df[!,c]),typeof(subtotalkw)
            @assert eltype(df[!,c]) == typeof(subtotalkw)
        end
        @assert !(any(isequal(subtotalkw),df[!,c])) #otherwise the appending will not be meaningful, as we set the values to missing for columns which are not considered in the multi way summary
    end

    for subsetlength = length(v):-1:0
        for subs in IterTools.subsets(v,subsetlength)
            #@show subs
            if subsetlength == 0
                agg = DataFrames.combine(df,cs...)
            else
                agg = DataFrames.combine(DataFrames.groupby(df,subs),cs...)
            end
            nonAggregatedVars = setdiff(v,subs)

            k = 1
            DataFrames.insertcols!(agg,k,:_TYPE_ => repeat(vcat(subsetlength),size(agg,1)))
            k += 1
            for addcol in nonAggregatedVars
                DataFrames.insertcols!(agg,k,addcol => repeat(vcat(subtotalkw),size(agg,1)))
                k += 1
            end

            DataFrames.allowmissing!(agg)
            DataFrames.append!(res,agg)
        end
    end

    sort!(res,vcat(:_TYPE_,v))
    return res
end

export addkey!

#adds comma separated row identifier to each row
function addkey!(df,v::Vector{Symbol};keyname=:key)
    @assert issubset(v,propertynames(df))
    @assert !in(keyname,propertynames(df))
    @assert !in(keyname,v)

    df[!,keyname]=repeat(vcat(""),size(df,1))
    for i=1:size(df,1)
        df[i,keyname]  = join(df[i,v],",")
    end

    have=vcat(keyname)
    donthave=setdiff(propertynames(df),have)
    select!(df,vcat(have,donthave))

    return nothing
end

end # module


text:
using BinaryBuilder

name = "Eigen"
version = v"3.3.9"

sources = [
    ArchiveSource("https://gitlab.com/libeigen/eigen/-/archive/$(version)/eigen-$(version).tar.bz2",
                  "0fa5cafe78f66d2b501b43016858070d52ba47bd9b1016b0165a7b8e04675677")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/eigen-*

mkdir build
cd build/

cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN} \
    -DCMAKE_Fortran_COMPILER=/opt/${MACHTYPE}/bin/${MACHTYPE}-gfortran \
    ..
make -j${nproc}
make install
cd ..
install_license COPYING.*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
# No products: Eigen is a pure header library
products = Product[
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using OASA
using Base.Test

# Synthesize random data
X = hcat(ones(100, 1), rand(100, 3))
truth_theta = [.5 1.2 3.3 5.0]
Y = truth_theta * X'

# Run some basic tests
function TestA()
    yh = infer(myModel, X)
    @test size(yh) == size(Y)
    @test assess(myModel, Y, yh) > 0.0
    @test size(assessgrad(myModel, Y, yh, X)) == (1, size(X,2))
end

# Ensure that all of the objective functions are operational
@testset "Objective Functions" begin
    @testset "L1 Regression" begin
        myModel = Model(L1_Obj(),
                    Linear(),
                    GradientDescent(),
                    SingleLayer(zeros(1, size(X,2))),
                    Dict())
        TestA()
    end
    @testset "L2 Regression" begin
        myModel.obj = L2_Obj()
        TestA()
    end
    @testset "Logistic Regression" begin
        myModel.obj = Logit_Obj()
        TestA()
    end
end


text:
# This file was generated, do not modify it. # hide
df, df_test, df_dumped = partition(df0, 0.07, 0.03, # in ratios 7:3:90
                                   stratify=df0.Churn,
                                   rng=123);

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module DCEMRI

using ArgParse
using Calculus    #.jacobian
using MAT, PyPlot
using Pkg, LinearAlgebra, Random, Statistics, Distributed, Printf, LsqFit

export ser, r1eff, tissueconc, fitr1, fitdce, fitdata,
  defaults, ccc, nlsfit, makeplots, demo, validate

const verbose = true
const version = v"0.2.2"

include("util.jl")
include("fitting.jl")
include("models.jl")
include("science.jl")
include("plotting.jl")
include("demo.jl")
include("validate.jl")

end


text:
# system/z-list,jl

include("Afft.jl")


text:
include("src/Fabrics.jl")
using Makie

function main()
    # Get Eulerian mesh
    M = EulerianMesh()
    # Get Lagrangian particles
    P = LagrangianParticles(M.Δx,M.Δy)
    # Save initial particles
    P0 = deepcopy(P)
    #Eulerian velocity field
    U = EulerianVelocity(M)
    ∇U = VelocityGradient(M)
    # Allocate gradiend of deformation
    F = [Matrix(1.0I,3,3) for i in 1:P.n]

    # Prepare Loop
    t = 0.0
    sec2year = 365.25*86400
    timemax = 1e6*sec2year
    Δt = min(M.Δx/2/maximum(abs.(U.x)), M.Δy/2/maximum(abs.(U.y)))
    cyclemax = div(timemax, Δt);

    # Start loop
    l = fill(0.0,3,3)
    RK4 = RK(fill(0.0,3,3), fill(0.0,3,3), fill(0.0,3,3), fill(0.0,3,3))
    for cycle in 1:cyclemax
        for m in 1:P.n
            interp∇U!(l,m,∇U,P,M)
            F[m] = updateFSE(F[m],l,Δt,RK4)
            advectparticle!(P,M,U,Δt,m)
        end
        t += Δt
    end

    return U,∇U,F,P,P0
end

U,∇U,F3D,P,P0 = main()
scatter(P.x/1e3, P.y/1e3, color=P0.y, markersize=2)

max_stars_repo_name: albert-de-montserrat/ForwardFSE

max_issues_repo_name: albert-de-montserrat/ForwardFSE

max_forks_repo_name: albert-de-montserrat/ForwardFSE


text:
# Defines CellScalarValues and CellVectorValues and common methods
"""
    CellScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
    CellVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])

A `CellValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions,
values of nodal functions, gradients and divergences of nodal functions etc. in the finite element cell. There are
two different types of `CellValues`: `CellScalarValues` and `CellVectorValues`. As the names suggest, `CellScalarValues`
utilizes scalar shape functions and `CellVectorValues` utilizes vectorial shape functions. For a scalar field, the
`CellScalarValues` type should be used. For vector field, both subtypes can be used.

**Arguments:**
* `T`: an optional argument (default to `Float64`) to determine the type the internal data is stored as.
* `quad_rule`: an instance of a [`QuadratureRule`](@ref)
* `func_interpol`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
* `geom_interpol`: an optional instance of a [`Interpolation`](@ref) which is used to interpolate the geometry

**Common methods:**
* [`reinit!`](@ref)
* [`getnquadpoints`](@ref)
* [`getdetJdV`](@ref)

* [`shape_value`](@ref)
* [`shape_gradient`](@ref)
* [`shape_symmetric_gradient`](@ref)
* [`shape_divergence`](@ref)

* [`function_value`](@ref)
* [`function_gradient`](@ref)
* [`function_symmetric_gradient`](@ref)
* [`function_divergence`](@ref)
* [`spatial_coordinate`](@ref)
"""
CellValues

# CellScalarValues
struct CellScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: CellValues{dim,T,refshape}
    N::Matrix{T}
    dNdx::Matrix{Vec{dim,T}}
    dNdξ::Matrix{Vec{dim,T}}
    detJdV::Vector{T}
    M::Matrix{T}
    dMdξ::Matrix{Vec{dim,T}}
    qr_weights::Vector{T}
end

function CellScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol)
    CellScalarValues(Float64, quad_rule, func_interpol, geom_interpol)
end

function CellScalarValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol) where {dim,T,shape<:AbstractRefShape}

    @assert getdim(func_interpol) == getdim(geom_interpol)
    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
    n_qpoints = length(getweights(quad_rule))

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol)
    N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
    dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)

    # Geometry interpolation
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)

    for (qp, ξ) in enumerate(quad_rule.points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
        end
        for i in 1:n_geom_basefuncs
            dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
        end
    end

    detJdV = fill(T(NaN), n_qpoints)

    CellScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
end

# CellVectorValues
struct CellVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: CellValues{dim,T,refshape}
    N::Matrix{Vec{dim,T}}
    dNdx::Matrix{Tensor{2,dim,T,M}}
    dNdξ::Matrix{Tensor{2,dim,T,M}}
    detJdV::Vector{T}
    M::Matrix{T}
    dMdξ::Matrix{Vec{dim,T}}
    qr_weights::Vector{T}
end

function CellVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
    CellVectorValues(Float64, quad_rule, func_interpol, geom_interpol)
end

function CellVectorValues(::Type{T}, quad_rule::QuadratureRule{dim,shape}, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol) where {dim,T,shape<:AbstractRefShape}

    @assert getdim(func_interpol) == getdim(geom_interpol)
    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
    n_qpoints = length(getweights(quad_rule))

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol) * dim
    N    = fill(zero(Vec{dim,T})      * T(NaN), n_func_basefuncs, n_qpoints)
    dNdx = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)

    # Geometry interpolation
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)

    for (qp, ξ) in enumerate(quad_rule.points)
        basefunc_count = 1
        for basefunc in 1:getnbasefunctions(func_interpol)
            dNdξ_temp, N_temp = gradient(ξ -> value(func_interpol, basefunc, ξ), ξ, :all)
            for comp in 1:dim
                N_comp = zeros(T, dim)
                N_comp[comp] = N_temp
                N[basefunc_count, qp] = Vec{dim,T}((N_comp...,))

                dN_comp = zeros(T, dim, dim)
                dN_comp[comp, :] = dNdξ_temp
                dNdξ[basefunc_count, qp] = Tensor{2,dim,T}((dN_comp...,))
                basefunc_count += 1
            end
        end
        for basefunc in 1:n_geom_basefuncs
            dMdξ[basefunc, qp], M[basefunc, qp] = gradient(ξ -> value(geom_interpol, basefunc, ξ), ξ, :all)
        end
    end

    detJdV = fill(T(NaN), n_qpoints)
    MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))

    CellVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, M, dMdξ, quad_rule.weights)
end

function reinit!(cv::CellValues{dim}, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_geom_basefuncs = getngeobasefunctions(cv)
    n_func_basefuncs = getnbasefunctions(cv)
    @assert length(x) == n_geom_basefuncs

    @inbounds for i in 1:length(cv.qr_weights)
        w = cv.qr_weights[i]
        fecv_J = zero(Tensor{2,dim})
        for j in 1:n_geom_basefuncs
            fecv_J += x[j] ⊗ cv.dMdξ[j, i]
        end
        detJ = det(fecv_J)
        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
        cv.detJdV[i] = detJ * w
        Jinv = inv(fecv_J)
        for j in 1:n_func_basefuncs
            cv.dNdx[j, i] = cv.dNdξ[j, i] ⋅ Jinv
        end
    end
end


text:
# this small script takes as an input an example
# file and outputs a formatted file with the outputs

include(ARGS[1])


text:
# neural network factor implementation

cd(@__DIR__)
using Pkg
pkg"activate ."

# using Revise

using Flux
using DifferentialEquations, OrdinaryDiffEq
using DiffEqFlux
using Optim
using Zygote, DiffEqSensitivity

import Base: convert
import IncrementalInference: getSample

# export NeuralFactor

struct NeuralSetBinaryFactor{T<:Vector{Union{Chain,FastChain}},D} <: FunctorPairwiseMinimize
  Zij::T
  data::D
end

function getSample(nfb::NeuralSetBinaryFactor, N::Int=1)::Tuple
  # evaluate
  len = length(nfg.Zij)
  spmls = []
  lbls = Vector{Int}(undef, len)
  sel = Categorical(ones(len)*(1/len))
  for i in 1:N
    lb = rand(sel)
    lbls[i] = lb
    val = 0 # evaluate nfb.Zij[lb]
    push!(smpls, val)
  end
  return (smpls, lbls)
end

function (nfb::NeuralSetBinaryFactor)(
            res::Array{Float64},
            userdata::FactorMetadata,
            idx::Int,
            meas::Tuple,
            xi::Array{Float64,2},
            xj::Array{Float64,2}  )
  #
  len = size(xi, 1)
  fill!(res, 0.0)
  for i in 1:len
    res[1] += meas[1][i,idx] - (xj[i,idx] - xi[i,idx])
  end
  res[1]
end

#


text:
# Solve generalized assignment problem with column generation
using StudentHousing, StatsBase, Distributions, JuMP, Gurobi

srand(32)
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.5, 0.5])
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.5, 0.5])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0]
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)
budget = 1e6
problem_data = StudentHousingData(market_data, nhouses = 2, budget = budget, demand_distribution = Uniform())

d = problem_data

nhouses = length(d.houses)
npatterns = length(d.patterns_allow)
println(npatterns) # 13

h_init = collect(1:nhouses)
V = zeros(npatterns, nhouses)
for i = 1:nhouses
    fsb_pattern = house_allowedby(i,d)[1]
    V[fsb_pattern, i] = 1.0
end

# find one feasible assignment to house 1
# julia> house_allowedby(1,d)
# 5-element Array{Int64,1}:
#   6
#   7
#  11
#  12
#  13
# V1 = zeros(npatterns)
# V1[6] = 1.0
# we could have started with 1 column for each house consisting of zeroes
# everywhere and 1 in the first pattern allowed

c = ones(npatterns)

# =============================================================================
# Master problem
# =============================================================================
m = Model(solver = GurobiSolver(OutputFlag=0))
@variable(m, 0 <= λ[i in h_init] <= 1)
@objective(m, Max, sum(λ[i] * dot(c, V[:, i]) for i in h_init))
@constraints(m, begin
    capsupply[p = 1:npatterns], sum(λ[i] * V[p, i] for i in h_init) <= 1
    convexity[i in h_init], sum(λ[i]) <= 1
end)

solve(m)
print(m)
println(getobjectivevalue(m))
π = getdual(m[:capsupply])
w = getdual(m[:convexity])
println("Dual vector = ", π)
println("Convexity dual = ", w)

# =============================================================================
# Sub-problems
# =============================================================================
best_rc = -Inf
best_sp  = Model()
best_i = 0

# The knapsack problem for each house
function get_kp(i::Int)
    sp = Model(solver = GurobiSolver(OutputFlag=0))
    @variable(sp, v[p = 1:npatterns], Bin)
    @constraints(sp, begin
        sum(v[p] * beds_needed(p) for p = 1:npatterns) <= beds_avail(i, d.houses)
        [p = 1:npatterns], v[p] <= house_fits_pattern(i, p, d)
    end)
    @objective(sp, Max, dot((c - π), v))
    sp
end

# Get the best reduced cost over all houses
for i = 1:nhouses
    sp = get_kp(i)
    @assert solve(sp) == :Optimal
    col_rc = getobjectivevalue(sp) - w[i]
    if col_rc > best_rc
        best_rc = col_rc
        best_sp = sp
        best_i = i
    end
end

@show best_rc
@assert best_rc > 0
@show best_i
@show find(getvalue(best_sp[:v]) .> 0)
# We chose to generate a column for house 1, which accommodates students 7
# and 11

# sanity check
StudentHousing.house_allowedby(1, d)
# 5-element Array{Int64,1}:
#   6
#   7
#  11
#  12
#  13
beds_avail(1, d.houses)
# 2
beds_needed(7)
# 1.0
beds_needed(11)
# 1.0


text:
"""
    walk_tree!(g, labels, ex, show_call=true)

Walk the abstract syntax tree (AST) of the given expression `ex`.
Builds up the graph `g` and the set of `labels`
for each node, both modified in place

`show_call` specifies whether to include `call` nodes in the graph.
Including them represents the Julia AST more precisely, but adds visual noise.

Returns the number of the top vertex.
"""
function walk_tree!(g, labels, ex, show_call=true)

    top_vertex = add_numbered_vertex!(g)

    where_start = 1  # which argument to start with

    if !(show_call) && ex.head == :call
        f = ex.args[1]   # the function name
        push!(labels, f)

        where_start = 2   # drop "call" from tree

    else
        push!(labels, ex.head)
    end

    for i in where_start:length(ex.args)

        if isa(ex.args[i], Expr)

            child = walk_tree!(g, labels, ex.args[i], show_call)
            add_edge!(g, top_vertex, child)

        else
            n = add_numbered_vertex!(g)
            add_edge!(g, top_vertex, n)

            push!(labels, ex.args[i])

        end
    end

    return top_vertex

end

function walk_tree(ex::Expr, show_call=false)
    g = DiGraph()
    labels = Any[]

    walk_tree!(g, labels, ex, show_call)

    return LabelledTree(g, labels)

end

"""
Make a tree from a Julia `Expr`ession.
Omits `call`.
"""
macro tree(ex::Expr)
    walk_tree(ex)
end

"""
Make a tree from a Julia `Expr`ession.
Includes `call`.
"""
macro tree_with_call(ex::Expr)
    walk_tree(ex, true)
end


text:
mutable struct SimpleCell <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    length::Float64
    orientation::Float64
    growthprog::Float64
    growthrate::Float64

    ## node positions/forces
    p1::NTuple{2,Float64}
    p2::NTuple{2,Float64}
    f1::NTuple{2,Float64}
    f2::NTuple{2,Float64}
end

function SimpleCell(id, pos, l, φ, g, γ)
    a = SimpleCell(id, pos, l, φ, g, γ, (0.0, 0.0), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
    update_nodes!(a)
    return a
end

"""
``` julia
growing_bacteria()
```
Same as in [Bacterial Growth](@ref).
"""
function growing_bacteria()
    space = ContinuousSpace((14, 9), 1.0; periodic = false)
    model = ABM(
        SimpleCell,
        space,
        properties = Dict(:dt => 0.005, :hardness => 1e2, :mobility => 1.0),
    )

    add_agent!((6.5, 4.0), model, 0.0, 0.3, 0.0, 0.1)
    add_agent!((7.5, 4.0), model, 0.0, 0.0, 0.0, 0.1)

    return model, growing_bacteria_agent_step!, growing_bacteria_model_step!
end

function update_nodes!(a::SimpleCell)
    offset = 0.5 * a.length .* unitvector(a.orientation)
    a.p1 = a.pos .+ offset
    a.p2 = a.pos .- offset
end

# Some geometry convenience functions
unitvector(φ) = reverse(sincos(φ))
cross2D(a, b) = a[1] * b[2] - a[2] * b[1]

# ## Stepping functions

function growing_bacteria_model_step!(model)
    for a in allagents(model)
        if a.growthprog ≥ 1
            ## When a cell has matured, it divides into two daughter cells on the
            ## positions of its nodes.
            add_agent!(a.p1, model, 0.0, a.orientation, 0.0, 0.1 * rand(model.rng) + 0.05)
            add_agent!(a.p2, model, 0.0, a.orientation, 0.0, 0.1 * rand(model.rng) + 0.05)
            kill_agent!(a, model)
        else
            ## The rest lengh of the internal spring grows with time. This causes
            ## the nodes to physically separate.
            uv = unitvector(a.orientation)
            internalforce = model.hardness * (a.length - a.growthprog) .* uv
            a.f1 = -1 .* internalforce
            a.f2 = internalforce
        end
    end
    ## Bacteria can interact with more than on other cell at the same time, therefore,
    ## we need to specify the option `:all` in `interacting_pairs`
    for (a1, a2) in interacting_pairs(model, 2.0, :all)
        interact!(a1, a2, model)
    end
end

function growing_bacteria_agent_step!(agent::SimpleCell, model::ABM)
    fsym, compression, torque = transform_forces(agent)
    new_pos = agent.pos .+ model.dt * model.mobility .* fsym
    move_agent!(agent, new_pos, model)
    agent.length += model.dt * model.mobility .* compression
    agent.orientation += model.dt * model.mobility .* torque
    agent.growthprog += model.dt * agent.growthrate
    update_nodes!(agent)
    return agent.pos
end

# ### Helper functions
function interact!(a1::SimpleCell, a2::SimpleCell, model)
    n11 = noderepulsion(a1.p1, a2.p1, model)
    n12 = noderepulsion(a1.p1, a2.p2, model)
    n21 = noderepulsion(a1.p2, a2.p1, model)
    n22 = noderepulsion(a1.p2, a2.p2, model)
    a1.f1 = @. a1.f1 + (n11 + n12)
    a1.f2 = @. a1.f2 + (n21 + n22)
    a2.f1 = @. a2.f1 - (n11 + n21)
    a2.f2 = @. a2.f2 - (n12 + n22)
end

function noderepulsion(p1::NTuple{2,Float64}, p2::NTuple{2,Float64}, model::ABM)
    delta = p1 .- p2
    distance = norm(delta)
    if distance ≤ 1
        uv = delta ./ distance
        return (model.hardness * (1 - distance)) .* uv
    end
    return (0, 0)
end

function transform_forces(agent::SimpleCell)
    ## symmetric forces (CM movement)
    fsym = agent.f1 .+ agent.f2
    ## antisymmetric forces (compression, torque)
    fasym = agent.f1 .- agent.f2
    uv = unitvector(agent.orientation)
    compression = dot(uv, fasym)
    torque = 0.5 * cross2D(uv, fasym)
    return fsym, compression, torque
end


text:
"""
    ModelNet40(;kwargs...)

Returns ModelNet40 dataset.

### Optional Key Arguments:
    * `root::String=default_root`   - Root directory of dataset
    * `train::Bool=true`            - Specifies the trainset
    * `transform=nothing`           - Transform to be applied to data point.
    * `categories::Vector{String}`  - Specifies the categories to be used in dataset.

### Examples:

```jldoctest
julia> dset = ModelNet40(train=false)
julia> typeof(dset[1].data) == TriMesh
```
"""
ModelNet40(;
    root::String = default_root,
    train::Bool = true,
    download::Bool = false,
    transform = nothing,
    categories::Vector{String} = MN10_classes,
) = _mn_dataset(root, train, download, transform, categories, 40)

const MN40_classes = [
    "airplane",
    "bathtub",
    "bed",
    "bench",
    "bookshelf",
    "bottle",
    "bowl",
    "car",
    "chair",
    "cone",
    "cup",
    "curtain",
    "desk",
    "door",
    "dresser",
    "flower_pot",
    "glass_box",
    "guitar",
    "keyboard",
    "lamp",
    "laptop",
    "mantel",
    "monitor",
    "night_stand",
    "person",
    "piano",
    "plant",
    "radio",
    "range_hood",
    "sink",
    "sofa",
    "stairs",
    "stool",
    "table",
    "tent",
    "toilet",
    "tv_stand",
    "vase",
    "wardrobe",
    "xbox",
]


text:
"""
    module Simulate_FallingBall4

Models and simulates a sphere that is falling downwards as 6 dof body
(described by quaternions), **with** collision handling.

See also models

- [`Simulate_FallingBall1`](@ref) (falling ball constrained with prismatic joint),
- [`Simulate_FallingBall2`](@ref) (falling ball constrained with prismatic joint and collision handling).
"""
module Simulate_FallingBall4

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial  = Modia3D.Material(color="DarkGreen", transparency=0.5)
ballMaterial    = Modia3D.Material(color="LightBlue", transparency=0.5)
contactMaterial = Modia3D.ContactMaterialElastic(mu1=0.5)

@assembly FallingBall4(;h_start = 1.0, D=0.1) begin
   world  = Object3D(Modia3D.CoordinateSystem(0.5))
   ground = Object3D(world, Modia3D.Solid(Modia3D.SolidBox(2.0,0.1,2.0), nothing, groundMaterial,
                                          contactMaterial=contactMaterial), r = [1.0, -0.05, 0.0])
   ball   = Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(D), "Aluminium", ballMaterial,
                                          contactMaterial=contactMaterial),
                     fixed=false, r=[0.5,h_start,0.0], v_start=[0.2,0,0], visualizeFrame=true)
end

fallingBall = FallingBall4()

# Modia3D.visualizeAssembly!( fallingBall )

simulationModel = Modia3D.SimulationModel( fallingBall, stopTime=2.5, scaleConstraintsAtEvents = true)
#Modia3D.print_ModelVariables(simulationModel)
result = ModiaMath.simulate!(simulationModel, log=true)

ModiaMath.plot(result, ["ball.r[2]"; "ball.v[1]"; "ball.v[2]"; "ball.v[3]"])

println("... success of Simulate_FallingBall4.jl!")

end


text:
export OLSModel, OLSData, OLSParams, add_intercept

"""
Model data with a scalar- or vector-valued ordinary least squares regression.
See [`OLSData`](@ref) for wrapping data.
"""
struct OLSModel end

"""
    OLSData(Y, X)

Ordinary least squares with dependent variable `Y` and design matrix `X`.

Either

1. ``Y`` is an ``n×m`` matrix, then ``Y = X B + E`` where ``X`` is a ``n×k``
    matrix, ``B`` is a ``k×m`` parameter matrix, and ``Eᵢ ∼ N(0, Σ)`` is IID
    with ``m×m`` variance matrix ``Σ`` (multivariate linear regression), or

2. ``Y`` is a length ``n`` vector, then ``Y = X β + ϵ``, where ``X`` is a
    ``n×k`` matrix, ``β`` is a parameter vector of ``k`` elements, and
    ``ϵᵢ∼N(0,σ)`` where ``σ`` is the variance of the normal error.

See also [`add_intercept`](@ref).
"""
struct OLSData{TY <: Union{AbstractMatrix, AbstractVector},
               TX <: AbstractMatrix}
    Y::TY
    X::TX
    function OLSData(Y::TY, X::TX) where {TY, TX}
        @argcheck size(Y, 1) == size(X, 1)
        new{TY, TX}(Y, X)
    end
end

function show(io::IO, data::OLSData{<:AbstractVector})
    n, k = size(data.X)
    print(io, "OLS regression data, $(n) scalar observations, $(k) covariates")
end

function show(io::IO, data::OLSData{<:AbstractMatrix})
    n, k = size(data.X)
    m = size(data.Y, 2)
    print(io, "OLS regression data, $(n) vector observations of length $(m), $(k) covariates")
end

"""
    OLSParams(B, Σ)

Maximum likelihood estimated parameters for an OLS regression. See
[`OLSData`](@ref).
"""
struct OLSParams{TB <: Union{AbstractMatrix, AbstractVector},
                  TΣ <: Union{<:Real, AbstractMatrix}}
    B::TB
    Σ::TΣ
    function OLSParams(B::TB, Σ::TΣ) where {TB <: AbstractVector, TΣ <: Real}
        new{TB, TΣ}(B, Σ)
    end
    function OLSParams(B::TB, Σ::TΣ) where {TB <: AbstractMatrix,
                                             TΣ <: AbstractMatrix}
        m = size(B, 2)
        @argcheck size(Σ) == (m, m)
        new{TB, TΣ}(B, Σ)
    end
end

function show(io::IO, params::OLSParams{<:AbstractVector})
    k = length(params.B)
    print(io, "OLS regression parameters, $(k) covariates")
end

function show(io::IO, params::OLSParams{<:AbstractMatrix})
    k, m = size(params.B)
    print(io, "OLS regression parameters, vector observations of length $(m), $(k) covariates")
end

function MLE(::OLSModel, data::OLSData)
    @unpack Y, X = data
    B = qrfact(X, Val{true}) \ Y
    E = Y - X*B
    Σ = E'*E / size(X, 1)
    OLSParams(B, Σ)
end

function loglikelihood(::OLSModel, data::OLSData{<: AbstractVector},
                       params::OLSParams{<: AbstractVector})
    @unpack Y, X = data
    @unpack B, Σ = params
    n = length(Y)
    @argcheck size(X, 1) == n
    @argcheck size(X, 2) == length(B)
    E = Y - X * B
    -0.5 * (n*(log(2*π) + log(Σ)) + sum(abs2, E)/Σ)
end

"""
    $SIGNATURES

Sum of the log pdf for observations ``xᵢ ∼ MultivariateNormal(0, Σ)``, where
``xᵢ`` are stored in the rows of `X`.

Handles `ForwardDiff.Dual` types and singular covariance matrices. Helper
function, not part of the API.
"""
function _logpdf_normal(X::AbstractMatrix{TX},
                        Σ::AbstractMatrix{TΣ}) where {TX, TΣ}
    n, m = size(X)
    try
        U = chol(Σ)
        A = X / U
        -0.5 * (n*(m*log(2*π) + 2*logdet(U)) + sum(abs2, A))
    catch
        # NOTE: technically, should check for X=0 in this case, but that almost
        # never happens and would need to throw an error, -Inf rejects without a
        # hassle too
        convert(promote_type(TX, TΣ), -Inf)
    end
end

function loglikelihood(::OLSModel, data::OLSData{<: AbstractMatrix},
                       params::OLSParams{<: AbstractMatrix})
    @unpack Y, X = data
    @unpack B, Σ = params
    n, m = size(Y)
    @argcheck size(B, 1) == size(X, 2) # k
    @argcheck size(B, 2) == m
    _logpdf_normal(Y - X * B, Σ)
end

"""
    $SIGNATURES

Add an intercept to a matrix or vector of covariates.
"""
add_intercept(X::AbstractMatrix{T}) where T = hcat(ones(T, size(X, 1)), X)

add_intercept(x::AbstractVector{T}) where T = hcat(ones(T, length(x)), x)

add_intercept(xs...) = add_intercept(hcat(xs...))

vec_parameters(ϕ::OLSParams) = vec_parameters((ϕ.B, ϕ.Σ))


text:
module RobotPattern

include("./../../constraints/geometric/constants.jl")
include("./1_equalities.jl")
include("./2_bounds.jl")
include("./3_constraints.jl")
using JuMP
using .Constants
using .RPEqualities
using .RPBounds
using .RPConstraints

export cons_rp_robot_pattern

function cons_rp_robot_pattern(m)

    m = rp_equalities(m)
    m = rp_bounds(m)
    m = rp_constraints(m)

    return m
end

end


text:
# Unit tests for implementations/nbodysimulator/nbodysimulator_base.jl

@testset "implementations/nbodysimulator/nbodysimulator_base.jl" begin
    body1 = ElementMassBody((@SVector [1, 2, 3])u"nm", (@SVector [4, 5, 6])u"nm/s", :Ar; meta="data")

    @test body1 isa ElementMassBody{Float64,Float64}
    @test body1.r ≈ (@SVector [1, 2, 3])austrip(1u"nm")
    @test body1.v ≈ (@SVector [4, 5, 6])austrip(1u"nm/s")
    @test body1.m ≈ 39.9481austrip(1u"u")
    @test body1.symbol == :Ar
    @test body1.data == Dict(:meta => "data")

    atom1 = AtomsBase.Atom(:Ar, [-2, -1, 3]u"bohr"; hello=:world)
    body2 = ElementMassBody(atom1)

    @test body2 isa ElementMassBody{Float64,Float64}
    @test body2.r ≈ @SVector [-2, -1, 3]
    @test body2.v ≈ @SVector [0, 0, 0]
    @test body2.m ≈ 39.9481austrip(1u"u")
    @test body2.symbol == :Ar
    @test body2.data == Dict(:hello => :world)

    nbs_boundary_conditions = CubicPeriodicBoundaryConditions(1.5)
    atom2 = AtomsBase.Atom(body2, nbs_boundary_conditions)

    @test position(atom2) ≈ (@SVector [1.0, 0.5, 0.0])u"bohr"
    @test velocity(atom2) ≈ velocity(atom1)
    @test atomic_mass(atom2) ≈ atomic_mass(atom1)
    @test atomic_symbol(atom2) == atomic_symbol(atom1)
    @test atomic_number(atom2) == atomic_number(atom1)
    @test atom2.data == atom1.data

    atom3 = AtomsBase.Atom(body2, (@SVector [2, 3, 4]), (@SVector [5, 6, 7]), nbs_boundary_conditions)

    @test position(atom3) ≈ (@SVector [0.5, 0.0, 1.0])u"bohr"
    @test velocity(atom3) ≈ (@SVector [5, 6, 7])u"bohr * hartree / ħ_au"
    @test atomic_mass(atom3) ≈ atomic_mass(atom1)
    @test atomic_symbol(atom3) == atomic_symbol(atom1)
    @test atomic_number(atom3) == atomic_number(atom1)
    @test atom3.data == atom1.data

    atom4 = AtomsBase.Atom(body2, (@SVector [1.0, 0.5, 0.0])u"bohr", (@SVector [1.0, 0.5, 0.0])u"bohr/s")

    @test position(atom4) ≈ (@SVector [1.0, 0.5, 0.0])u"bohr"
    @test velocity(atom4) ≈ (@SVector [1.0, 0.5, 0.0])u"bohr/s"
    @test atomic_mass(atom4) ≈ atomic_mass(atom1)
    @test atomic_symbol(atom4) == atomic_symbol(atom1)
    @test atomic_number(atom4) == atomic_number(atom1)
    @test atom4.data == atom1.data

    particles = [
        AtomsBase.Atom(:Ar, [-2, -1, 3]u"bohr", [3, 5, 21]u"bohr * hartree / ħ_au"; a=:b),
        AtomsBase.Atom(:Ar, [-7, 2, 13]u"bohr", [-3, 7, 1]u"bohr * hartree / ħ_au"; c=:d)
    ]
    box = [[1.5, 0.0, 0.0], [0.0, 1.5, 0.0], [0.0, 0.0, 1.5]]u"bohr"
    system = periodic_system(particles, box)

    @test Atomistic.nbs_boundary_conditions(system) == nbs_boundary_conditions
    @test get_boundary_conditions(nbs_boundary_conditions) == [Periodic(), Periodic(), Periodic()]
    @test get_bounding_box(nbs_boundary_conditions) == box

    @test Atomistic.bound_position((@SVector [1, -1, 10]), nbs_boundary_conditions) == (@SVector [1, 0.5, 1])
end


text:
const TESTDATA = Dict("number_of_females_1"=> 2,
             "number_of_females_2"=> 3,
             "number_of_females_3"=> 3,
             "number_of_males_1"=> 11,
             "number_of_males_2"=> 13,
             "number_of_males_3"=> 42,
             "number_of_initial_vacancies_1"=> 0,
             "number_of_initial_vacancies_2"=> 0,
             "number_of_initial_vacancies_3"=> 0,
             "hiring_rate_women_1"=> 0.2059,
             "hiring_rate_women_2"=> 0.0588,
             "hiring_rate_women_3"=> 0.00,
             "hiring_rate_men_1"=> 0.5294,
             "hiring_rate_men_2"=> 0.1176,
             "hiring_rate_men_3"=> 0.0882,
             "attrition_rate_women_1"=> 0.0615,
             "attrition_rate_women_2"=> 0.00,
             "attrition_rate_women_3"=> 0.06,
             "attrition_rate_men_1"=> 0.0859,
             "attrition_rate_men_2"=> 0.0473,
             "attrition_rate_men_3"=> 0.0414,
             "female_promotion_probability_1"=> 0.0769,
             "female_promotion_probability_2"=> 0.1111,
             "male_promotion_probability_1"=> 0.0707,
             "male_promotion_probability_2"=> 0.0946,
             "t_fpct"=> 0.15,
             "upperbound"=> 94,
             "lowerbound"=> 74,
             "duration"=> 20,
             "annotation_text"=> "this is an annotation",
             "growth_rate"=> [0.01, 0.01, -0.01, -0.01, 0.01, -0.01, 0.01],
             "model_name" => "model 3 with drift")


text:
# Lucas Ondel, 2021

using SemifieldAlgebra
using Test

const TropicalSemifield{T} = Semifield{T, max, +, -, -Inf, 0} where T
const _x = 1.6
const _y = -11.3

for T in [Float32, Float64]
    x = TropicalSemifield{T}(_x)
    y = TropicalSemifield{T}(_y)

    @test (x + y).val ≈ _x
    @test (x * y).val ≈ _x + _y
    @test (x / y).val ≈ _x - _y
    @test zero(typeof(x)).val == -Inf
    @test zero(TropicalSemifield{T}).val == -Inf
    @test one(typeof(x)).val == 0
    @test one(TropicalSemifield{T}).val == 0
end


text:
pos1 = Real2D(1.0,1.0)
pos2 = Real2D(1.0,1.0)
pos3 = Real2D(10.0,10.0)

@test (pos1 - pos2) == Real2D(0.0,0.0)
@test (pos1 + pos2) == Real2D(2.0,2.0)
@test pos1 == pos2

@test toroidal(Real2D(11.2,11.2),10.0,10.0) == Real2D(1.2,1.2)

@test toroidal(Real2D(11.2,11.2),10.0,10.0) != Real2D(1.3,1.3)

@test toroidal(Real2D(11.2,11.2),3.0,3.0) == Real2D(2.2,2.2)

@test distance(Real2D(1.0,1.0), Real2D(1.0,5.0), 10.0, 10.0, false) == 4.0

@test distance(Real2D(1.0,1.0), Real2D(1.0,5.0), 10.0, 10.0, true) == 4.0

@test distance(Real2D(1.0,1.0), Real2D(1.0,5.0), 4.0, 4.0, true) == 0.0