text:
"""
    μ ± σ

Creates $DEFAULT_NUM_PARTICLES `Particles` with mean `μ` and std `σ`.
If `μ` is a vector, the constructor `MvNormal` is used, and `σ` is thus treated as std if it's a scalar, and variances if it's a matrix or vector.
See also [`∓`](@ref), [`..`](@ref)
"""
±

"""
    μ ∓ σ

Creates $DEFAULT_STATIC_NUM_PARTICLES `StaticParticles` with mean `μ` and std `σ`.
If `μ` is a vector, the constructor `MvNormal` is used, and `σ` is thus treated as std if it's a scalar, and variances if it's a matrix or vector.
See also [`±`](@ref), [`⊗`](@ref)
"""
∓

±(μ::Real,σ) = Particles{promote_type(float(typeof(μ)),float(typeof(σ))),DEFAULT_NUM_PARTICLES}(systematic_sample(DEFAULT_NUM_PARTICLES,Normal(μ,σ); permute=true))
±(μ::AbstractVector,σ) = Particles(DEFAULT_NUM_PARTICLES, MvNormal(μ, σ))
∓(μ::Real,σ) = StaticParticles{promote_type(float(typeof(μ)),float(typeof(σ))),DEFAULT_STATIC_NUM_PARTICLES}(systematic_sample(DEFAULT_STATIC_NUM_PARTICLES,Normal(μ,σ); permute=true))
∓(μ::AbstractVector,σ) = StaticParticles(DEFAULT_STATIC_NUM_PARTICLES, MvNormal(μ, σ))

"""
    a .. b

Creates $DEFAULT_NUM_PARTICLES `Particles` with a `Uniform` distribution between `a` and `b`.
See also [`±`](@ref), [`⊗`](@ref)
"""
(..)(a,b) = Particles{float(promote_type(eltype(a), eltype(b))), DEFAULT_NUM_PARTICLES}(Random.GLOBAL_RNG, Uniform(a,b))

"""
    a ⊠ Distribution()

Multiplies `a` by $DEFAULT_NUM_PARTICLES `Particles` sampled from a specified `::Distribution`.
Shorthand for `a * Particles(Distribution())`, e.g., `a ⊠ Gamma(1)`.
"""
⊠(a,d::Distribution) = a * Particles{eltype(d), DEFAULT_NUM_PARTICLES}(Random.GLOBAL_RNG, d)

"""
    a ⊞ Distribution()

Adds $DEFAULT_NUM_PARTICLES `Particles` sampled from a specified `::Distribution` to `a`.
Shorthand for `a + Particles(Distribution())`, e.g., `1 ⊞ Binomial(3)`.
"""
⊞(a,d::Distribution) = a + Particles{eltype(d), DEFAULT_NUM_PARTICLES}(Random.GLOBAL_RNG, d)

"""
    ⊗(μ,σ) = outer_product(Normal.(μ,σ))

See also [`outer_product`](@ref), [`±`](@ref)
"""
⊗(μ,σ) = outer_product(Normal.(μ,σ))

"""
    p = outer_product([rng::AbstractRNG,] dists::Vector{<:Distribution}, N=100_000)

Creates a multivariate systematic sample where each dimension is sampled according to the corresponding univariate distribution in `dists`. Returns `p::Vector{Particles}` where each Particles has a length approximately equal to `N`.
The particles form the outer product between `d` systematically sampled vectors with length given by the d:th root of N, where `d` is the length of `dists`, All particles will be independent and have marginal distributions given by `dists`.

See also `MonteCarloMeasurements.⊗`
"""
function outer_product(rng::AbstractRNG, dists::AbstractVector{<:Distribution}, N=100_000)
    d = length(dists)
    N = floor(Int,N^(1/d))
    dims = map(dists) do dist
        v = systematic_sample(rng,N,dist; permute=true)
    end
    cart_prod = vec(collect(Iterators.product(dims...)))
    p = map(1:d) do i
        Particles(getindex.(cart_prod,i))
    end
end
function outer_product(dists::AbstractVector{<:Distribution}, N=100_000)
    return outer_product(Random.GLOBAL_RNG, dists, N)
end

# StaticParticles(N::Integer = DEFAULT_NUM_PARTICLES; permute=true) = StaticParticles{Float64,N}(SVector{N,Float64}(systematic_sample(N, permute=permute)))

function print_functions_to_extend()
    excluded_functions = [fill, |>, <, display, show, promote, promote_rule, promote_type, size, length, ndims, convert, isapprox, ≈, <, (<=), (==), zeros, zero, eltype, getproperty, fieldtype, rand, randn]
    functions_to_extend = setdiff(names(Base), Symbol.(excluded_functions))
    for fs in functions_to_extend
        ff = @eval $fs
        ff isa Function || continue
        isempty(methods(ff)) && continue # Sort out intrinsics and builtins
        f = nameof(ff)
        if !isempty(methods(ff, (Real,Real)))
            println(f, ",")
        end
    end
end
"""
    shortform(p::AbstractParticles)
Return a short string describing the type
"""
shortform(p::Particles) = "Part"
shortform(p::StaticParticles) = "SPart"
function to_num_str(p::AbstractParticles{T}, d=3, ds=d-1) where T
    s = pstd(p)
    # TODO: be smart and select sig digits based on s
    if T <: AbstractFloat && s < eps(p)
        string(round(pmean(p), sigdigits=d))
    else
        string(round(pmean(p), sigdigits=d), " ± ", round(s, sigdigits=ds))
    end
end

function Base.show(io::IO, p::AbstractParticles{T,N}) where {T,N}
    print(io, to_num_str(p, 3))
end

function Base.show(io::IO, ::MIME"text/plain", p::AbstractParticles{T,N}) where {T,N}
    sPT = MonteCarloMeasurements.shortform(p)
    compact = get(io, :compact, false)
    if compact
        print(io, MonteCarloMeasurements.to_num_str(p, 6, 3))
    else
        print(io, "$(typeof(p))\n ", MonteCarloMeasurements.to_num_str(p, 6, 3))
    end
end

function Base.show(io::IO, z::Complex{<:AbstractParticles})
    r, i = reim(z)
    compact = get(io, :compact, false)
    print(io, "(")
    show(io, r)
    print(io, ")")
    if pmaximum(i) < 0
        i = -i
        print(io, compact ? "-" : " - ")
    else
        print(io, compact ? "+" : " + ")
    end
    print(io, "(")
    show(io, i)
    print(io, ")")
    print(io, "im")
end

# function Base.show(io::IO, p::MvParticles)
#     sPT = shortform(p)
#     print(io, "(", N, " $sPT with mean ", round.(mean(p), sigdigits=3), " and std ", round.(sqrt.(diag(cov(p))), sigdigits=3),")")
# end
for mime in (MIME"text/x-tex", MIME"text/x-latex")
    @eval function Base.show(io::IO, ::$mime, p::AbstractParticles)
        print(io, "\$"); show(io, p); print("\$")
    end

    @eval function Base.show(io::IO, ::$mime, z::Complex{<:AbstractParticles})
        print(io, "\$")
        r, i = reim(z)
        compact = get(io, :compact, false)
        print(io, "(")
        show(io, r)
        print(io, ")")
        if pmaximum(i) < 0
            i = -i
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        print(io, "(")
        show(io, i)
        print(io, ")")
        print(io, "i")
        print("\$")
    end
end

# Two-argument functions
# foreach(register_primitive_binop, [+,-,*,/,//,^])
foreach(register_primitive_multi, [+,-,*,/,//,^,max,min,mod,mod1,atan,atand,add_sum,hypot])
# One-argument functions
foreach(register_primitive_single, [+,-,
exp,exp2,exp10,expm1,
log,log10,log2,log1p,
sin,cos,tan,sind,cosd,tand,sinh,cosh,tanh,
asin,acos,atan,asind,acosd,atand,asinh,acosh,atanh,
zero,sign,abs,sqrt,rad2deg,deg2rad])

MvParticles(x::AbstractVector{<:AbstractArray{<:Number}}) = Particles(copy(reduce(hcat, x)'))
MvParticles(v::AbstractVector{<:Number}) = Particles(v)

function MvParticles(v::AbstractVector{<:Tuple})
    Particles.([getindex.(v,i) for i in 1:length(v[1])])
end

function MvParticles(s::Vector{NamedTuple{vs, T}}) where {vs, T}
    nt = NamedTuple()
    for k in keys(s[1])
        nt = merge(nt, [k => MvParticles(getproperty.(s,k))])
    end
    nt
end

function _finish_individuals(::Type{PT}, N, individuals::AbstractArray{<:Tuple}, p) where PT
    ntuple(length(individuals[1])) do ti
        RT = eltype(first(individuals)[ti])
        PRT = PT{RT,N}
        no = length(individuals[1][ti])
        out = Vector{PRT}(undef, no)
        for i = 1:no
            out[i] = PRT(getindex.(getindex.(individuals,ti), i))
        end
        reshape(out, size(individuals[1][ti]))
    end
end

function _finish_individuals(::Type{PT}, ::Val{N}, individuals, p) where {PT, N}
    RT = eltype(eltype(individuals))
    PRT = PT{RT,N}
    out = similar(p, PRT)
    for i = 1:length(p)
        out[i] = PRT(getindex.(individuals,i))
    end
    reshape(out, size(p))
end

for PT in ParticleSymbols
    # Constructors
    @eval begin

        """
            ℝⁿ2ℝⁿ_function(f::Function, p::AbstractArray{T})

        Helper function for performing uncertainty propagation through vector-valued functions with vector inputs.
        Applies  `f : ℝⁿ → ℝⁿ` to an array of particles. E.g., `Base.log(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(log,p)`
        """
        function ℝⁿ2ℝⁿ_function(f::F, p::AbstractArray{$PT{T,N}}) where {F,T,N}
            individuals = map(1:nparticles(p[1])) do i
                f(vecindex.(p,i))
            end
            _finish_individuals($PT, Val{N}(), individuals, p)
        end

        function ℝⁿ2ℝⁿ_function(f::F, p::AbstractArray{$PT{T,N}}, p2::AbstractArray{$PT{T,N}}) where {F,T,N}
            individuals = map(1:nparticles(p[1])) do i
                f(vecindex.(p,i), vecindex.(p2,i))
            end
            _finish_individuals($PT, Val{N}(), individuals, p)
        end

        """
            ℝⁿ2ℂⁿ_function(f::Function, p::AbstractArray{T})

        Helper function for performing uncertainty propagation through complex-valued functions with vector inputs.
        Applies  `f : ℝⁿ → Cⁿ` to an array of particles. E.g., `LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℂⁿ_function(eigvals,p)`
        """
        function ℝⁿ2ℂⁿ_function(f::F, p::AbstractArray{$PT{T,N}}; kwargs...) where {F,T,N}
            individuals = map(1:nparticles(p[1])) do i
                f(vecindex.(p,i); kwargs...)
            end
            PRT = $PT{T,N}
            RT = eltype(eltype(individuals))
            if RT <: Complex
                CRT = Complex{PRT}
            else
                CRT = PRT
            end
            out = Array{CRT}(undef, size(individuals[1]))
            for i = eachindex(out)
                ind = getindex.(individuals,i)
                if RT <: Complex
                    out[i] = complex(PRT(real.(ind)), PRT(imag.(ind)))
                else
                    out[i] = PRT(ind)
                end
            end
            out
        end
        #
        # function ℝⁿ2ℂⁿ_function(f::F, p::AbstractArray{$PT{T,N}}, p2::AbstractArray{$PT{T,N}}) where {F,T,N}
        #     individuals = map(1:nparticles(p[1])) do i
        #         f(getindex.(p,i), getindex.(p2,i))
        #     end
        #     RT = eltype(eltype(individuals))
        #     @assert RT <: Complex
        #     PRT = $PT{T,N}
        #     CRT = Complex{PRT}
        #     out = similar(p, CRT)
        #     for i = 1:length(p)
        #         ind = getindex.(individuals,i)
        #         out[i] = complex(PRT(real.(ind)), PRT(imag.(ind)))
        #     end
        #     reshape(out, size(p))
        # end
    end
end

for ff in (var, std)
    f = nameof(ff)
    @eval function (Statistics.$f)(p::ParticleDistribution{<:AbstractParticles{T,N}}, args...; kwargs...) where {N,T}
        N == 1 && (return zero(T))
        $f(p.p.particles, args...;kwargs...)
    end
    pname = Symbol("p"*string(f))
    @eval function ($pname)(p::AbstractParticles{T,N}, args...; kwargs...) where {N,T}
        N == 1 && (return zero(T))
        $f(p.particles, args...;kwargs...)
    end
end
# Instead of @forward
# TODO: convert all these to operate on ParticleDistribution
for ff in [Statistics.mean, Statistics.cov, Statistics.median, Statistics.quantile, Statistics.middle, Base.iterate, Base.extrema, Base.minimum, Base.maximum]
    f = nameof(ff)
    m = Base.parentmodule(ff)
    @eval ($m.$f)(p::ParticleDistribution, args...; kwargs...) = ($m.$f)(p.p.particles, args...; kwargs...)
    pname = Symbol("p"*string(f))
    @eval ($pname)(p::AbstractParticles, args...; kwargs...) = ($m.$f)(p.particles, args...; kwargs...)
    @eval ($pname)(p::Number, args...; kwargs...) = ($m.$f)(p, args...; kwargs...)
end

for PT in ParticleSymbols

    @eval begin
        # Base.length(::Type{$PT{T,N}}) where {T,N} = N
        # Base.eltype(::Type{$PT{T,N}}) where {T,N} = $PT{T,N} # TODO: remove

        Base.convert(::Type{StaticParticles{T,N}}, p::$PT{T,N}) where {T,N} = StaticParticles(p.particles)
        Base.convert(::Type{$PT{T,N}}, f::Real) where {T,N} = $PT{T,N}(fill(T(f),N))
        Base.convert(::Type{$PT{T,N}}, f::$PT{S,N}) where {T,N,S} = $PT{promote_type(T,S),N}(convert.(promote_type(T,S),f.particles))
        function Base.convert(::Type{S}, p::$PT{T,N}) where {S<:ConcreteFloat,T,N}
            N == 1 && (return S(p.particles[1]))
            pstd(p) < eps(S) || throw(ArgumentError("Cannot convert a particle distribution to a float if not all particles are the same."))
            return S(p.particles[1])
        end
        function Base.convert(::Type{S}, p::$PT{T,N}) where {S<:ConcreteInt,T,N}
            isinteger(p) || throw(ArgumentError("Cannot convert a particle distribution to an int if not all particles are the same."))
            return S(p.particles[1])
        end
        Base.zeros(::Type{$PT{T,N}}, dim::Integer) where {T,N} = [$PT{T,N}(zeros(eltype(T),N)) for d = 1:dim]
        Base.zero(::Type{$PT{T,N}}) where {T,N} = $PT{T,N}(zeros(eltype(T),N))
        Base.isfinite(p::$PT{T,N}) where {T,N} = isfinite(pmean(p))
        Base.round(p::$PT{T,N}, r::RoundingMode, args...; kwargs...) where {T,N} = round(pmean(p), r, args...; kwargs...)
        Base.round(::Type{S}, p::$PT{T,N}, args...; kwargs...) where {S,T,N} = round(S, pmean(p), args...; kwargs...)
        function Base.AbstractFloat(p::$PT{T,N}) where {T,N}
            N == 1 && (return p.particles[1])
            pstd(p) < eps(T) || throw(ArgumentError("Cannot convert a particle distribution to a number if not all particles are the same."))
            return p.particles[1]
        end

        """
            union(p1::AbstractParticles, p2::AbstractParticles)

        A `Particles` containing all particles from both `p1` and `p2`. Note, this will be twice as long as `p1` or `p2` and thus of a different type.
        `pu = Particles([p1.particles; p2.particles])`
        """
        function Base.union(p1::$PT{T,NT},p2::$PT{T,NS}) where {T,NT,NS}
            $PT{T,NT+NS}([p1.particles; p2.particles])
        end

        """
            intersect(p1::AbstractParticles, p2::AbstractParticles)

        A `Particles` containing all particles from the common support of `p1` and `p2`. Note, this will be of undetermined length and thus undetermined type.
        """
        function Base.intersect(p1::$PT,p2::$PT)
            mi = max(pminimum(p1),pminimum(p2))
            ma = min(pmaximum(p1),pmaximum(p2))
            f = x-> mi <= x <= ma
            $PT([filter(f, p1.particles); filter(f, p2.particles)])
        end

        function Base.:^(p::$PT{T,N}, i::Integer) where {T,N} # Resolves ambiguity
            res = p.particles.^i
            $PT{eltype(res),N}(res)
        end
        Base.:\(p::Vector{<:$PT}, p2::Vector{<:$PT}) = Matrix(p)\Matrix(p2) # Must be here to be most specific

        function LinearAlgebra.eigvals(p::Matrix{$PT{T,N}}; kwargs...) where {T,N} # Special case to propte types differently
            individuals = map(1:N) do i
                eigvals(vecindex.(p,i); kwargs...)
            end
            PRT = Complex{$PT{T,N}}
            out = Vector{PRT}(undef, length(individuals[1]))
            for i = eachindex(out)
                c = getindex.(individuals,i)
                out[i] = complex($PT{T,N}(real(c)),$PT{T,N}(imag(c)))
            end
            out
        end

    end

    # for XT in (:T, :($PT{T,N})), YT in (:T, :($PT{T,N})), ZT in (:T, :($PT{T,N}))
    #     XT == YT == ZT == :T && continue
    #     @eval function Base.muladd(x::$XT,y::$YT,z::$ZT) where {T<:Number,N}
    #          res = muladd.(maybe_particles(x),maybe_particles(y),maybe_particles(z))
    #          $PT{eltype(res),N}(res)
    #     end
    # end

    @eval function Base.muladd(x::$PT{T,N},y::$PT{T,N},z::$PT{T,N}) where {T<:Number,N}
        res = muladd.(x.particles,y.particles,z.particles)
        $PT{T,N}(res)
    end
    @eval function Base.muladd(x::T,y::$PT{T,N},z::$PT{T,N}) where {T<:Number,N}
        res = muladd.(x,y.particles,z.particles)
        $PT{T,N}(res)
    end
    @eval function Base.muladd(x::T,y::T,z::$PT{T,N}) where {T<:Number,N}
        res = muladd.(x,y,z.particles)
        $PT{T,N}(res)
    end

    @eval Base.promote_rule(::Type{S}, ::Type{$PT{T,N}}) where {S<:Number,T,N} = $PT{promote_type(S,T),N} # This is hard to hit due to method for real 3 lines down
    @eval Base.promote_rule(::Type{Bool}, ::Type{$PT{T,N}}) where {T,N} = $PT{promote_type(Bool,T),N}

    for PT2 in ParticleSymbols
        if PT == PT2
            @eval Base.promote_rule(::Type{$PT{S,N}}, ::Type{$PT{T,N}}) where {S,T,N} = $PT{promote_type(S,T),N}
        elseif any(==(:StaticParticles), (PT, PT2))
            @eval Base.promote_rule(::Type{$PT{S,N}}, ::Type{$PT2{T,N}}) where {S,T,N} = StaticParticles{promote_type(S,T),N}
        else
            @eval Base.promote_rule(::Type{$PT{S,N}}, ::Type{$PT2{T,N}}) where {S,T,N} = Particles{promote_type(S,T),N}
        end
    end

    @eval Base.promote_rule(::Type{<:AbstractParticles}, ::Type{$PT{T,N}}) where {T,N} = Union{}
end

# Base.length(p::AbstractParticles{T,N}) where {T,N} = N
Base.ndims(p::AbstractParticles{T,N}) where {T,N} = ndims(T)
Base.:\(H::MvParticles,p::AbstractParticles) = Matrix(H)\p.particles
# Base.:\(p::AbstractParticles, H) = p.particles\H
# Base.:\(p::MvParticles, H) = Matrix(p)\H
# Base.:\(H,p::MvParticles) = H\Matrix(p)

Base.Broadcast.broadcastable(p::AbstractParticles) = Ref(p)
# Base.setindex!(p::AbstractParticles, val, i::Integer) = setindex!(p.particles, val, i)
# Base.getindex(p::AbstractParticles, i::Integer) = getindex(p.particles, i)
# Base.getindex(v::MvParticles, i::Int, j::Int) = v[j][i] # Defining this methods screws with show(::MvParticles)

Base.Array(p::AbstractParticles) = p.particles
Base.Vector(p::AbstractParticles) = Array(p)

function Base.Array(v::Array{<:AbstractParticles})
    m = reduce(hcat, Array.(v))
    return reshape(m, size(m, 1), size(v)...)
end
Base.Matrix(v::MvParticles) = Array(v)

# function Statistics.var(v::MvParticles,args...;kwargs...) # Not sure if it's a good idea to define this. Is needed for when var(v::AbstractArray) is used
#     s2 = map(1:length(v[1])) do i
#         var(getindex.(v,i))
#     end
#     eltype(v)(s2)
# end

pmean(v::MvParticles) = pmean.(v)
pcov(v::MvParticles,args...;kwargs...) = cov(Matrix(v), args...; kwargs...)
pcor(v::MvParticles,args...;kwargs...) = cor(Matrix(v), args...; kwargs...)
pvar(v::MvParticles,args...; corrected = true, kwargs...) = sum(abs2, v)/(nparticles(v) - corrected)
Distributions.fit(d::Type{<:MultivariateDistribution}, p::MvParticles) = fit(d,Matrix(p)')
Distributions.fit(d::Type{<:Distribution}, p::AbstractParticles) = fit(d,p.particles)

Distributions.Normal(p::AbstractParticles) = Normal(pmean(p), pstd(p))
Distributions.MvNormal(p::MvParticles) = MvNormal(pmean(p), pcov(p))

meanstd(p::AbstractParticles) = pstd(p)/sqrt(nparticles(p))
meanvar(p::AbstractParticles) = pvar(p)/nparticles(p)

Base.:(==)(p1::AbstractParticles{T,N},p2::AbstractParticles{T,N}) where {T,N} = p1.particles == p2.particles
Base.:(!=)(p1::AbstractParticles{T,N},p2::AbstractParticles{T,N}) where {T,N} = p1.particles != p2.particles
Base.hash(p::AbstractParticles) = hash(p.particles) + hash(typeof(p))

function zip_longest(a_,b_)
    a,b = maybe_particles(a_), maybe_particles(b_)
    l = max(length(a), length(b))
    Iterators.take(zip(Iterators.cycle(a), Iterators.cycle(b)), l)
end

function safe_comparison(a_, b_, op::F) where F
    a,b = maybe_particles(a_), maybe_particles(b_)
    all(((a,b),)->op(a,b), Iterators.product(extrema(a),extrema(b))) && (return true)
    !any(((a,b),)->op(a,b), Iterators.product(extrema(a),extrema(b))) && (return false)
    _comparison_error()
end

function do_comparison(a,b,op::F) where F
    mode = COMPARISON_MODE[]
    if mode === :reduction
        op(COMPARISON_FUNCTION[](a), COMPARISON_FUNCTION[](b))
    elseif mode === :montecarlo
        all(((a,b),)->op(a,b), zip_longest(a,b)) && return true
        !any(((a,b),)->op(a,b), zip_longest(a,b)) && return false
        _comparison_error()
    elseif mode === :safe
        safe_comparison(a,b,op)
    else
        error("Got unsupported comparison mode.")
    end
end

function _comparison_error()
    msg = "Comparison of uncertain values using comparison mode $(COMPARISON_MODE[]) failed. Comparison operators are not well defined for uncertain values. Call `unsafe_comparisons(true)` to enable comparison operators for particles using the current reduction function $(COMPARISON_FUNCTION[]). Change this function using `set_comparison_function(f)`. "
    if COMPARISON_MODE[] === :safe
        msg *= "For safety reasons, the default safe comparison function is maximally conservative and tests if the extreme values of the distributions fulfil the comparison operator."
    elseif COMPARISON_MODE[] === :montecarlo
        msg *= "For safety reasons, montecarlo comparison is conservative and tests if pairwise particles fulfil the comparison operator. If some do *and* some do not, this error is thrown. Consider if you can define a primitive function ([docs](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/overloading/#Overloading-a-new-function-1)) or switch to `unsafe_comparisons(:reduction)`"
    end

    error(msg)
end

function Base.:<(a::Real,p::AbstractParticles)
    do_comparison(a,p,<)
end
function Base.:<(p::AbstractParticles,a::Real)
    do_comparison(p,a,<)
end
function Base.:<(p::AbstractParticles, a::AbstractParticles)
    do_comparison(p,a,<)
end
function Base.:(<=)(p::AbstractParticles{T,N}, a::AbstractParticles{T,N}) where {T,N}
    do_comparison(p,a,<=)
end

"""
    p1 ≈ p2

Determine if two particles are not significantly different
"""
Base.:≈(p::AbstractParticles, a::AbstractParticles, lim=2) = abs(pmean(p)-pmean(a))/(2sqrt(pstd(p)^2 + pstd(a)^2)) < lim
function Base.:≈(a::Real,p::AbstractParticles, lim=2)
    m = pmean(p)
    s = pstd(p, mean=m)
    s == 0 && (return m == a)
    abs(pmean(p)-a)/pstd(p) < lim
end
function Base.:≈(p::AbstractParticles, a::Real, lim=2)
    m = pmean(p)
    s = pstd(p, mean=m)
    s == 0 && (return m == a)
    abs(pmean(p)-a)/pstd(p) < lim
end
Base.:≈(p::MvParticles, a::AbstractVector) = all(a ≈ b for (a,b) in zip(a,p))
Base.:≈(a::AbstractVector, p::MvParticles) = all(a ≈ b for (a,b) in zip(a,p))
Base.:≈(a::MvParticles, p::MvParticles) = all(a ≈ b for (a,b) in zip(a,p))
Base.:≉(a,b::AbstractParticles,lim=2) = !(≈(a,b,lim))
Base.:≉(a::AbstractParticles,b,lim=2) = !(≈(a,b,lim))
"""
    p1 ≉ p2

Determine if two particles are significantly different
"""
Base.:≉(a::AbstractParticles,b::AbstractParticles,lim=2) = !(≈(a,b,lim))

Base.sincos(x::AbstractParticles) = sin(x),cos(x)
Base.minmax(x::AbstractParticles,y::AbstractParticles) = (min(x,y), max(x,y))

Base.:!(p::AbstractParticles) = all(p.particles .== 0)

Base.isinteger(p::AbstractParticles) = all(isinteger, p.particles)
Base.iszero(p::AbstractParticles) = all(iszero, p.particles)
Base.iszero(p::AbstractParticles, tol) = abs(mean(p.particles)) < tol

≲(a,b,args...) = a < b
≲(a::Real,p::AbstractParticles,lim=2) = (pmean(p)-a)/pstd(p) > lim
≲(p::AbstractParticles,a::Real,lim=2) = (a-pmean(p))/pstd(p) > lim
≲(p::AbstractParticles,a::AbstractParticles,lim=2) = (pmean(a)-pmean(p))/(2sqrt(pstd(p)^2 + pstd(a)^2)) > lim
≳(a::Real,p::AbstractParticles,lim=2) = ≲(p,a,lim)
≳(p::AbstractParticles,a::Real,lim=2) = ≲(a,p,lim)
≳(p::AbstractParticles,a::AbstractParticles,lim=2) = ≲(a,p,lim)
Base.eps(p::Type{<:AbstractParticles{T,N}}) where {T,N} = eps(T)
Base.eps(p::AbstractParticles{T,N}) where {T,N} = eps(T)
Base.eps(p::AbstractParticles{<:Complex{T},N}) where {T,N} = eps(T)

Base.rtoldefault(::Type{<:AbstractParticles{T,N}}) where {T,N} = sqrt(eps(T))

LinearAlgebra.norm(x::AbstractParticles, args...) = abs(x)

Base.log(p::Matrix{<:AbstractParticles}) = ℝⁿ2ℂⁿ_function(log,p) # Matrix more specific than StridedMatrix used in Base.log
# LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles}; kwargs...) = ℝⁿ2ℂⁿ_function(eigvals,p; kwargs...) # Replaced with implementation below
Base.exp(p::AbstractMatrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(exp, p)
LinearAlgebra.lyap(p1::Matrix{<:AbstractParticles}, p2::Matrix{<:AbstractParticles}) = ℝⁿ2ℝⁿ_function(lyap, p1, p2)

## Particle BLAS

# pgemv is up to twice as fast as the naive way already for A(2,2)-A(20,20)
"""
    _pgemv(A, p::Vector{StaticParticles{T, N}}) where {T, N}

Perform `A*p::Vector{StaticParticles{T,N}` using BLAS matrix-matrix multiply. This function is automatically used when applicable and there is no need to call it manually.
"""
function _pgemv(
    A,
    p::Vector{StaticParticles{T,N}},
) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N}
    pm = reinterpret(T, p)
    M = reshape(pm, N, :)'
    AM = A * M
    reinterpret(StaticParticles{T,N}, vec(AM'))
end

Base.:*(A::Matrix{T}, p::Vector{StaticParticles{T,N}}) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N} = _pgemv(A,p)

"""
    _pdot(v::Vector{T}, p::Vector{StaticParticles{T, N}}) where {T, N}

Perform `v'p::Vector{StaticParticles{T,N}` using BLAS matrix-vector multiply. This function is automatically used when applicable and there is no need to call it manually.
"""
function _pdot(
    v::AbstractVector{T},
    p::Vector{StaticParticles{T,N}},
) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N}
    pm = reinterpret(T, p)
    M = reshape(pm, N, :)
    Mv = M*v
    StaticParticles{T,N}(Mv)
end

LinearAlgebra.dot(v::AbstractVector{T}, p::Vector{StaticParticles{T,N}}) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N} = _pdot(v,p)
LinearAlgebra.dot(p::Vector{StaticParticles{T,N}}, v::AbstractVector{T}) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N} = _pdot(v,p)

function _paxpy!(
    a::T,
    x::Vector{StaticParticles{T,N}},
    y::Vector{StaticParticles{T,N}},
) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N}
    X = reinterpret(T, x)
    Y = reinterpret(T, y)
    LinearAlgebra.axpy!(a,X,Y)
    reinterpret(StaticParticles{T,N}, Y)
end

LinearAlgebra.axpy!(
    a::T,
    x::Vector{StaticParticles{T,N}},
    y::Vector{StaticParticles{T,N}},
) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N} = _paxpy!(a,x,y)

function LinearAlgebra.mul!(
    y::Vector{StaticParticles{T,N}},
    A::AbstractMatrix{T},
    b::Vector{StaticParticles{T,N}},
) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N}
    Bv = reinterpret(T, b)
    B = reshape(Bv, N, :)'
    # Y0 = A*B
    # reinterpret(StaticParticles{T,N}, vec(Y0'))
    Yv = reinterpret(T, y)
    Y = reshape(Yv, :, N)
    mul!(Y,A,B)
    reinterpret(StaticParticles{T,N}, vec(Y'))
end

function LinearAlgebra.mul!(
    y::Vector{Particles{T,N}},
    A::AbstractMatrix{T},
    b::Vector{Particles{T,N}},
) where {T<:Union{Float32,Float64,ComplexF32,ComplexF64},N}

    B = Matrix(b)
    # Y = A*B
    Y = B*A' # This order makes slicing below more efficient
    @inbounds if isdefined(y, 1)
        for i in eachindex(y)
            @views y[i].particles .= Y[:,i]
        end
    else
        for i in eachindex(y)
            y[i] = Particles(Y[:,i])
        end
    end
    y
end


text:
using Distributed
using DelimitedFiles
using Printf

@everywhere using RAFF

# This script runs the parallel version of RAFF

"""

    run_praff()

Load and run the parallel/distributed version of RAFF. It assumes that
there is a problem file `/tmp/output.txt`.

"""
function run_praff(maxms=1, initguess=nothing; model_str="logistic", foutliers=0.5)

    n, model, modelstr = RAFF.model_list[model_str]

    open("/tmp/output.txt") do fp

        global N = parse(Int, readline(fp))

        global data = readdlm(fp)

    end

    if initguess == nothing

        initguess = zeros(Float64, n)

    end

    rsol = praff(model, data[:, 1:end - 1], n; MAXMS=maxms, initguess=initguess,
                 noutliers=Int(round(foutliers * size(data)[1])))

    @printf("Solution found:
            fbest = %f
            p     = %d\n", rsol.f, rsol.p)
    println(rsol.solution)

    return rsol

end

ext: jl


text:
# # Huber regression

# This example can be found here: <https://web.stanford.edu/~boyd/papers/pdf/cvx_applications.pdf>.
# Here we set `big_example = false` to only generate a small example which takes less time to run.
big_example = false
if big_example
    n = 300
    number_tests = 50
else
    n = 100
    number_tests = 20
end

# Generate data for Huber regression.
using Random
Random.seed!(1);
number_samples = round(Int,1.5*n);
beta_true = 5*randn(n);
X = randn(n, number_samples);
Y = zeros(number_samples);
v = randn(number_samples);

#-

## Generate data for different values of p.
## Solve the resulting problems.
using Convex, SCS, Distributions
lsq_data = zeros(number_tests);
huber_data = zeros(number_tests);
prescient_data = zeros(number_tests);
p_vals = range(0, stop=0.15, length=number_tests);
for i=1:length(p_vals)
    p = p_vals[i];
    ## Generate the sign changes.
    factor = 2 * rand(Binomial(1, 1-p), number_samples) .- 1;
    Y = factor .* X' * beta_true + v;

    ## Form and solve a standard regression problem.
    beta = Variable(n);
    fit = norm(beta - beta_true) / norm(beta_true);
    cost = norm(X' * beta - Y);
    prob = minimize(cost);
    solve!(prob, SCS.Optimizer(verbose=0));
    lsq_data[i] = evaluate(fit);

    ## Form and solve a prescient regression problem,
    ## i.e., where the sign changes are known.
    cost = norm(factor .* (X'*beta) - Y);
    solve!(minimize(cost), SCS.Optimizer(verbose=0))
    prescient_data[i] = evaluate(fit);

    ## Form and solve the Huber regression problem.
    cost = sum(huber(X' * beta - Y, 1));
    solve!(minimize(cost), SCS.Optimizer(verbose=0))
    huber_data[i] = evaluate(fit);
end

#-
using Plots

plot(p_vals, huber_data, label="Huber", xlabel="p", ylabel="Fit" )
plot!(p_vals, lsq_data, label="Least squares")
plot!(p_vals, prescient_data, label="Prescient")
#-

## Plot the relative reconstruction error for Huber and prescient regression,
## zooming in on smaller values of p.
indices = findall(p_vals .<= 0.08);
plot(p_vals[indices], huber_data[indices], label="Huber")
plot!(p_vals[indices], prescient_data[indices], label="Prescient")


text:
module V5RPC
export Strategy,EmptyStrategy,V5Server

module PB
(filename->include(joinpath(@__DIR__, "..", "deps", "julia_pb", filename))).(["DataStructures_pb.jl","Events_pb.jl","API_pb.jl"])
end
using .PB

abstract type Strategy end

function on_event(::Strategy, ::Int32, ::EventArguments)
    @warn "on_event not implemented"
    nothing
end

function get_team_info(::Strategy)
    @warn "get_team_info not implemented"
    ""
end

function get_team_info(s::Strategy, ::typeof(Version.V1_1))
    @warn "Fallback to v1.0 API"
    get_team_info(s)
end

function get_instruction(::Strategy, ::Field)
    @warn "get_instruction not implemented"
    (fill((0, 1), 5), ControlType.Continue)
end

function get_placement(::Strategy, ::Field)
    @warn "get_placement not implemented"
    fill((0, 0, 0), 5)
end

struct EmptyStrategy <: Strategy end
on_event(::EmptyStrategy, ::Int32, ::EventArguments) = nothing
get_team_info(::EmptyStrategy, ::typeof(Version.V1_1)) = "Empty Strategy"
get_instruction(::EmptyStrategy, ::Field) = fill((0, 0), 5)
get_placement(::EmptyStrategy, ::Field) = fill((0, 0, 0), 5)

destructure(v::Vector2) = (v.x, v.y)
destructure(b::Ball) = destructure(b.position)
destructure(w::Wheel) = (w.left_speed, w.right_speed)
destructure(r::Robot) = (
    position = destructure(r.position),
    rotation = r.rotation,
    wheel = destructure(r.wheel)
)
destructure(f::Field) = (
    self_robots = tuple(destructure.(f.self_robots)...),
    opponent_robots = tuple(destructure.(f.opponent_robots)...),
    ball = destructure(f.ball),
    tick = f.tick
)

include("V5Packet.jl")
using .ModV5Packet
include("V5Server.jl")
using .ModV5Server

end # module


text:
using MichiganAutonomousVehicles
using NLOptControl
using Plots

case_name = "s12"
c = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/planner/","RTPP",".yaml")))
c["vehicle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/vehicle/","hmmwv",".yaml")))
c["goal"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["goal"]
c["X0"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["X0"]
c["obstacle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["obstacle"]
setConfig(c,"misc";(:N=>40),(:model=>:ThreeDOFv2),(:solver=>:Ipopt),(:integrationScheme=>:trapezoidal))
setConfig(c,"X0";(:ux=>0.1))# the type needs to be a float
# otherwise, this will pop up: WARNING: Ipopt finished with status Invalid_Number_Detected
fixYAML(c)   # fix messed up data types
n = initializeAutonomousControl(c);
simMPC!(n;updateFunction=updateAutoParams!,checkFunction=checkCrash)

include("postProcess.jl")


text:
module Pastebin

import HTTP

const API = "https://pastebin.com/api/api_post.php"

struct Client
    devKey::String
    userKey::String

    """
        Client(devKey, userKey = "")

    Creates a `Pastebin.Client`. If `userKey` is empty, accesses are as guest.
    """
    Client(devKey::String, userKey::String = "") = new(devKey, userKey)

    """
        Client(devKey, username, password)

    Login, generating a `userKey` for the session.
    """
    function Client(devKey::String, username::String, password::String)
        data = Dict(
            "api_dev_key" => devKey,
            "api_user_name" => username,
            "api_user_password" => password,
        )
        response = String(HTTP.post("https://pastebin.com/api/api_login.php", [], HTTP.Form(data)).body)
        return new(devKey, response)
    end
end

function Base.show(io::IO, client::Client)
    if isempty(client.userKey)
        print(io, "Pastebin.Client(...) -- no userKey")
    else
        print(io, "Pastebin.Client(...) -- userKey set")
    end
end

struct Response
    successful::Bool
    response::String
end

function Response(response::AbstractString)
    successful = !contains("Bad API Request", response)
    return Response(successful, response)
end

"""
    paste(client, title, text; expire = "N", access = 0)

Submits `text` to pastebin. If `client` has a `userKey`, the paste will be created by
that user.

`expire` date is one of: `"N"` (never, default), `"10M"` (10 minutes), `"1H"` (1 hour),
`"1D"` (1 day), `"1W"` (1 week), `"2W"` (2 weeks), `"1M"` (1 month), `"6M"` (6 months),
`"1Y"` (1 year).

`access` is an integer: `0` (public), `1` (unlisted, default) or `2` (private, requires a `userKey`).
`format` determines syntax highlight (plain text by default).

See <https://pastebin.com/doc_api>.
"""
function paste(
    client::Client, title::String, text::String;
    expire::String="N", access::Int=1, format::String="text"
)
    if access == 2 && isempty(client.userKey)
        throw(ArgumentError("for a private paste, providing a user key is required"))
    end
    data = Dict(
        "api_option"            => "paste",
        "api_paste_private"     => string(access),
        "api_paste_name"        => title,
        "api_paste_expire_date" => expire,
        "api_paste_format"      => format,
        "api_dev_key"           => client.devKey,
        "api_user_key"          => client.userKey,
        "api_paste_code"        => text
    )
    response = String(HTTP.post(API, [], HTTP.Form(data)).body)
    return Response(response)
end

"""
    delete(client, key)

Deletes paste `key` from the server. Requires a `userKey`.
"""
function delete(client::Client, key::String)
    if isempty(client.userKey)
        throw(ArgumentError("providing a userKey is required to delete a paste"))
    end
    data = Dict(
        "api_option"    => "delete",
        "api_dev_key"   => client.devKey,
        "api_user_key"  => client.userKey,
        "api_paste_key" => key,
    )
    response = String(HTTP.post(API, [], HTTP.Form(data)).body)
    return Response(response)
end

"""
    getraw(key)

Fetch a raw paste by `key`. Returns a `response = Pastebin.Response(...)` object.
The raw paste string can be accessed as `Pastebin.content(response)`.
"""
function getraw(key::AbstractString)
    response = String(HTTP.get("https://pastebin.com/raw/$key").body)
    return Response(response)
end

content(resp::Response) = resp.response
success(resp::Response) = resp.successful
parsekey(url::String) = split(url, "/")[end]
parsekey(resp::Response) = parsekey(content(resp))

end # end of module


text:
import Term.Style: apply_style
import Term: tprint, tprintln

@testset "\e[34mStyle\e[0m" begin
    @test apply_style("test") == "test"

    @test apply_style("{red}test{/red}") == "\e[31mtest\e[39m"

    @test apply_style("{red}te{blue on_green}s{/blue on_green}t{/red}") ==
          "\e[31mte\e[34m\e[42ms\e[39m\e[49m\e[31mt\e[39m"

    @test apply_style("{red}\ntest{bold} sdfsfsd{/bold}sdfsdf{/red}") ==
          "\e[31m\ntest\e[1m sdfsfsd\e[22m\e[31msdfsdf\e[39m"

    @test apply_style("""
    test {red} sdfsdf
    fdsf{/red} {bold} sfsdfp{green} sdfsdp{/green}sdsdfs
    pdfsdp{/bold}""") ==
          "test \e[31m sdfsdf\nfdsf\e[39m \e[1m sfsdfp\e[32m sdfsdp\e[39m\e[1msdsdfs\npdfsdp\e[22m\e[39m"

    # check that parentheses are escaped correctly
    @test apply_style("This and that {{something}} for") ==
          "This and that {{something}} for"
end

@testset "\e[34mTprint\e[0m" begin
    stprint(x) = chomp(sprint(tprint, x; context = stdout))
    stprintln(x) = chomp(sprint(tprintln, x; context = stdout))

    # ------------------------------- with strings ------------------------------- #
    @test stprint("test") == "test"
    # @test stprintln("test") == "test "

    @test stprint("{red}test{/red}") == "\e[31mtest\e[39m"
    # @test stprintln("{red}test{/red}") == "\e[31mtest\e[39m "

    @test stprint("{red}te{blue on_green}s{/blue on_green}t{/red}") ==
          "\e[31mte\e[34m\e[42ms\e[39m\e[49m\e[31mt\e[39m"
    # @test stprintln("{red}te{blue on_green}s{/blue on_green}t{/red}") == "\e[31mte\e[34m\e[42ms\e[39m\e[49m\e[31mt\e[39m "

    @test stprint("{red}\ntest{bold} sdfsfsd{/bold}sdfsdf{/red}") ==
          "\e[31m\ntest\e[1m sdfsfsd\e[22m\e[31msdfsdf\e[39m"
    # @test stprintln("{red}\ntest{bold} sdfsfsd{/bold}sdfsdf{/red}") == "\e[31m\ntest\e[1m sdfsfsd\e[22m\e[31msdfsdf\e[39m "

    @test stprint("""
    test {red} sdfsdf
    fdsf{/red} {bold} sfsdfp{green} sdfsdp{/green}sdsdfs
    pdfsdp{/bold}""") ==
          "test \e[31m sdfsdf\nfdsf\e[39m \e[1m sfsdfp\e[32m sdfsdp\e[39m\e[1msdsdfs\npdfsdp\e[22m\e[39m"

    @test stprint("This and that {{something}} for") == "This and that {{something}} for"
end


text:
function mdl_pss_ode!(
    device_states,
    output_ode,
    ω_sys,
    device::PSY.DynamicGenerator{M, S, A, TG, PSY.PSSFixed},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov}

    #Update V_pss on inner vars
    get_inner_vars(device)[V_pss_var] = PSY.get_V_pss(PSY.get_pss(device))

    return
end

#Currently not working properly.
#=
function mdl_pss_ode!(
    device_states,
    output_ode,
    ω_sys,
    device::PSY.DynamicGenerator{M, S, A, TG, PSY.PSSSimple},
) where {M <: PSY.Machine, S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov}

    #Get references
    static_gen = PSY.get_static_injector(device)
    P0 = PSY.get_active_power(static_gen)

    #Obtain external states for device
    external_ix = get_input_port_ix(device, PSY.PSSSimple)
    ω = @view device_states[external_ix]

    #Define external inner vars for component
    V_th = sqrt(get_inner_vars(device)[VR_gen_var]^2 + get_inner_vars(device)[VI_gen_var]^2)
    τe = get_inner_vars(device)[τe_var]

    #Get parameters
    pss = PSY.get_pss(device)
    K_ω = PSY.get_K_ω(pss)
    K_p = PSY.get_K_p(pss)

    #Update V_pss on inner vars
    get_inner_vars(device)[V_pss_var] = K_ω * (ω[1] - ω_sys) + K_p * (ω[1] * τe - P0)

    return
end
=#


text:
mutable struct Statevector
    xhat::Vector
    P::Matrix

    function Statevector(xhat::Vector, P::Matrix)
        m = size(xhat)[1]
        n,o = size(P)

        m==n==o ? new(xhat, P) : error("P must be NxN and xhat Nx1")

    end
end

Statevector(xhat::Matrix, P::Matrix) = Statevector(xhat[:], P)


text:
using Revise
using FunctionWrappers: FunctionWrapper
using BenchmarkTools

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Expressions, with function wrappers
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#
struct RheoModelClass

    name::String
    params::Vector{Symbol}

    G::FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}}
    Ga::FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}}
    J::FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}}
    Ja::FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}}
    Gp::FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}}
    Gpa::FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}}
    Gpp::FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}}
    Gppa::FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}}

    constraint::FunctionWrapper{Bool,Tuple{Vector{Float64}}}

    info::String
    expressions::NamedTuple

end

function RheoModelClass(;name::String="Custom model",
                         p::Array{Symbol}=[],
                         G::Expr = quote nothing end,
                         J::Expr = quote nothing end,
                         Gp::Expr = quote nothing end,
                         Gpp::Expr = quote nothing end,
                         constraint::Expr = quote true end,
                         info=name)

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))
    expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint)

    @eval return(RheoModelClass($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end)                 |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end)  |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((t,params) -> begin $unpack_expr; $J; end)                 |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end)  |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((ω,params) -> begin $unpack_expr; $Gp; end)                |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((ω,params) -> begin $unpack_expr; $Gpp; end)               |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        (params -> begin $unpack_expr; $constraint; end)            |> FunctionWrapper{Bool,Tuple{Vector{Float64}}},
        $info, $expressions) )
end

expr_replace!(ex, s, v) = ex
expr_replace!(ex::Symbol, s, v) = s == ex ? v : ex

function expr_replace!(ex::Expr, s, v)
    for i=1:length(ex.args)
        ex.args[i] = expr_replace!(ex.args[i], s, v)
    end
    return ex
end

expr_replace(ex, s, v) = expr_replace!(copy(ex), s, v)

function expr_replace(ex, nt)
    e = copy(ex)
    k=keys(nt)
    v=values(nt)
    for i in 1:length(nt)
        expr_replace!(e, k[i], v[i])
    end
    return e
end

function model_parameters(nt::NamedTuple, params::Vector{Symbol}, err_string::String)
    # check that every parameter in m exists in the named tuple nt
    @assert all(i->i in keys(nt),params) "Missing parameter(s) in " * err_string
    # check that no extra parameters have been provided
    @assert length(params) == length(nt) "Mismatch number of model parameters and parameters provided in " * err_string

    p = map(i->Float64(nt[i]), params)
end

function freeze_params(m::RheoModelClass, nt0::NamedTuple)

    # check that every parameter in m exists in the named tuple nt
    @assert all( i-> i in m.params,keys(nt0)) "A parameter to freeze is not involved in the model"
    # convert values format for consistency
    nt=NamedTuple{keys(nt0)}([Float64(i) for i in nt0])

    # create array of remaining variables
    p = filter(s -> !(s in keys(nt)),m.params)

    name="$(m.name) with set parameters: $nt"
    info = m.info

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))

    # This section creates moduli expressions with material parameters
    # replaced by specific values.

    G = expr_replace(m.expressions.G, nt)
    J = expr_replace(m.expressions.J, nt)
    Gp = expr_replace(m.expressions.Gp, nt)
    Gpp = expr_replace(m.expressions.Gpp, nt)
    constraint = expr_replace(m.expressions.constraint, nt)

    expressions=NamedTuple{(:G,:J,:Gp,:Gpp)}( ( G, J, Gp, Gpp ) )

    @eval return( RheoModelClass($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end) |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((t,params) -> begin $unpack_expr; $J; end) |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((ω,params) -> begin $unpack_expr; $Gp; end) |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((ω,params) -> begin $unpack_expr; $Gpp; end) |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        (params -> begin $unpack_expr; $constraint; end) |> FunctionWrapper{Bool,Tuple{Vector{Float64}}},
        $info, $expressions)   )
end

SLS_exp_fw = RheoModelClass(
          # Model name
          name="SLS",
          # Model parameters,
          p = [ :η, :kᵦ, :kᵧ],
          # Relaxation modulus
          G = quote
                 kᵧ + kᵦ*exp(-t*kᵦ/η)
              end,
          # Creep modulus
          J = quote
                  c₀ = 1/kᵧ
                  c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
                  τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
                  c₀ - c₁*exp(-t/τᵣ)
              end,
          # Storage modulus
          Gp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω^2*τ^2*kᵦ
                  numerator/denominator + kᵧ
               end,
         # Loss modulus
          Gpp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω*τ*kᵦ
                  numerator/denominator
                end,
          # Network
          info= " "
          )

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Expressions, no function wrappers
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#

struct RheoModelClassNOFW

    name::String
    params::Vector{Symbol}

    G::Function
    Ga::Function
    J::Function
    Ja::Function
    Gp::Function
    Gpa::Function
    Gpp::Function
    Gppa::Function

    constraint::Function

    info::String
    expressions::NamedTuple

end

function RheoModelClassNOFW(;name::String="Custom model",
                         p::Array{Symbol}=[],
                         G::Expr = quote nothing end,
                         J::Expr = quote nothing end,
                         Gp::Expr = quote nothing end,
                         Gpp::Expr = quote nothing end,
                         constraint::Expr = quote true end,
                         info=name)

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))
    expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint)

    @eval return(RheoModelClassNOFW($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end),
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end),
        ((t,params) -> begin $unpack_expr; $J; end),
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end),
        ((ω,params) -> begin $unpack_expr; $Gp; end),
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end),
        ((ω,params) -> begin $unpack_expr; $Gpp; end),
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end),
        (params -> begin $unpack_expr; $constraint; end),
        $info, $expressions) )
end

function freeze_paramsNOFW(m::RheoModelClassNOFW, nt0::NamedTuple)
    # check that every parameter in m exists in the named tuple nt
    @assert all( i-> i in m.params,keys(nt0)) "A parameter to freeze is not involved in the model"
    # convert values format for consistency
    nt=NamedTuple{keys(nt0)}([Float64(i) for i in nt0])

    # create array of remaining variables
    p = filter(s -> !(s in keys(nt)),m.params)

    name="$(m.name) with set parameters: $nt"
    info = m.info

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))

    # This section creates moduli expressions with material parameters
    # replaced by specific values.
    G = expr_replace(m.expressions.G, nt)
    J = expr_replace(m.expressions.J, nt)
    Gp = expr_replace(m.expressions.Gp, nt)
    Gpp = expr_replace(m.expressions.Gpp, nt)
    constraint = expr_replace(m.expressions.constraint, nt)

    expressions=NamedTuple{(:G,:J,:Gp,:Gpp)}( ( G, J, Gp, Gpp ) )

    @eval return( RheoModelClass($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end),
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end),
        ((t,params) -> begin $unpack_expr; $J; end),
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end),
        ((ω,params) -> begin $unpack_expr; $Gp; end),
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end),
        ((ω,params) -> begin $unpack_expr; $Gpp; end),
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end),
        (params -> begin $unpack_expr; $constraint; end),
        $info,
        $expressions))
end

SLS_exp_NOFW = RheoModelClassNOFW(
          # Model name
          name="SLS",
          # Model parameters,
          p = [ :η, :kᵦ, :kᵧ],
          # Relaxation modulus
          G = quote
                 kᵧ + kᵦ*exp(-t*kᵦ/η)
              end,
          # Creep modulus
          J = quote
                  c₀ = 1/kᵧ
                  c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
                  τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
                  c₀ - c₁*exp(-t/τᵣ)
              end,
          # Storage modulus
          Gp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω^2*τ^2*kᵦ
                  numerator/denominator + kᵧ
               end,
         # Loss modulus
          Gpp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω*τ*kᵦ
                  numerator/denominator
                end,
          # Network
          info= " "
          )

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Functions, no freezing capability
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#
struct RheoModelClassNOFWNOEXP

    name::String
    params::Vector{Symbol}

    G::Function
    Ga::Function
    J::Function
    Ja::Function
    Gp::Function
    Gpa::Function
    Gpp::Function
    Gppa::Function

    constraint::Function

    info::String
    expressions::NamedTuple

end

function RheoModelClassNOFWNOEXP(;name::String="Custom model",
                         p::Array{Symbol}=[],
                         G::Function = (t, params) -> nothing,
                         J::Function = (t, params) -> nothing,
                         Gp::Function = (t, params) -> nothing,
                         Gpp::Function = (t, params) -> nothing,
                         constraint::Function = (params) -> true,
                         info=name)

    return(RheoModelClassNOFWNOEXP(name, p,
        G,
        (ta, params) -> G.(ta, (params,)),
        J,
        (ta, params) -> J.(ta, (params,)),
        Gp,
        (ωa, params) -> Gp.(ωa, (params,)),
        Gpp,
        (ωa, params) -> Gpp.(ωa, (params,)),
        constraint,
        info,
        (DefinedAsFunctions = 1,)))
end

function SLSFunkNOFWNOEXPG(t, params)
    η, kᵦ, kᵧ = params
    kᵧ + kᵦ*exp(-t*kᵦ/η)
end

function SLSFunkNOFWNOEXPJ(t, params)
    η, kᵦ, kᵧ = params
    c₀ = 1/kᵧ
    c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
    τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
    c₀ - c₁*exp(-t/τᵣ)
end

function SLSFunkNOFWNOEXPGp(ω, params)
    η, kᵦ, kᵧ = params
    τ = η/kᵦ
    denominator = 1 + τ^2*ω^2
    numerator = ω^2*τ^2*kᵦ
    numerator/denominator + kᵧ
end

function SLSFunkNOFWNOEXPGpp(ω, params)
    η, kᵦ, kᵧ = params
    τ = η/kᵦ
    denominator = 1 + τ^2*ω^2
    numerator = ω*τ*kᵦ
    numerator/denominator
end

SLS_exp_NOFWNOEXP = RheoModelClassNOFWNOEXP(
          # Model name
          name="SLS",
          # Model parameters,
          p = [ :η, :kᵦ, :kᵧ],
          # Relaxation modulus
          G = SLSFunkNOFWNOEXPG,
          # Creep modulus
          J = SLSFunkNOFWNOEXPJ,
          # Storage modulus
          Gp = SLSFunkNOFWNOEXPGp,
         # Loss modulus
          Gpp = SLSFunkNOFWNOEXPGpp,
          # Network
          info= " "
          )

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Parametrised seems would be messy due to all the differently parameterised FunctionWrappers, trying generic
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#
struct RheoModelClassGENERIC

    name::String
    params::Vector{Symbol}

    G
    Ga
    J
    Ja
    Gp
    Gpa
    Gpp
    Gppa

    constraint

    info::String
    expressions::NamedTuple

end

function RheoModelClassGENERICfromfunk(;name::String="Custom model",
                         p::Array{Symbol}=[],
                         G::Function = (t, params) -> nothing,
                         J::Function = (t, params) -> nothing,
                         Gp::Function = (t, params) -> nothing,
                         Gpp::Function = (t, params) -> nothing,
                         constraint::Function = (params) -> true,
                         info=name)

    return(RheoModelClassGENERIC(name, p,
        G,
        (ta, params) -> G.(ta, (params,)),
        J,
        (ta, params) -> J.(ta, (params,)),
        Gp,
        (ωa, params) -> Gp.(ωa, (params,)),
        Gpp,
        (ωa, params) -> Gpp.(ωa, (params,)),
        constraint,
        info,
        (DefinedAsFunctions = 1,)))
end

function RheoModelClassGENERICfromexp(;name::String="Custom model",
                         p::Array{Symbol}=[],
                         G::Expr = quote nothing end,
                         J::Expr = quote nothing end,
                         Gp::Expr = quote nothing end,
                         Gpp::Expr = quote nothing end,
                         constraint::Expr = quote true end,
                         info=name)

    # Expression to unpack parameter array into suitably names variables in the moduli expressions
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))
    expressions = (G=G,J=J,Gp=Gp,Gpp=Gpp,constraint=constraint)

    @eval return(RheoModelClassGENERIC($name, $p,
        ((t,params) -> begin $unpack_expr; $G; end)                 |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ta,params) -> begin $unpack_expr; [$G for t in ta]; end)  |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((t,params) -> begin $unpack_expr; $J; end)                 |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ta,params) -> begin $unpack_expr; [$J for t in ta]; end)  |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((ω,params) -> begin $unpack_expr; $Gp; end)                |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ωa,params) -> begin $unpack_expr; [$Gp for ω in ωa]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        ((ω,params) -> begin $unpack_expr; $Gpp; end)               |> FunctionWrapper{Float64,Tuple{Float64,Vector{Float64}}},
        ((ωa,params) -> begin $unpack_expr; [$Gpp for ω in ωa]; end) |> FunctionWrapper{Vector{Float64},Tuple{Vector{Float64},Vector{Float64}}},
        (params -> begin $unpack_expr; $constraint; end)            |> FunctionWrapper{Bool,Tuple{Vector{Float64}}},
        $info, $expressions) )
end

function SLSFunkGENERICG(t, params)
    η, kᵦ, kᵧ = params
    kᵧ + kᵦ*exp(-t*kᵦ/η)
end

function SLSFunkGENERICJ(t, params)
    η, kᵦ, kᵧ = params
    c₀ = 1/kᵧ
    c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
    τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
    c₀ - c₁*exp(-t/τᵣ)
end

function SLSFunkGENERICGp(ω, params)
    η, kᵦ, kᵧ = params
    τ = η/kᵦ
    denominator = 1 + τ^2*ω^2
    numerator = ω^2*τ^2*kᵦ
    numerator/denominator + kᵧ
end

function SLSFunkGENERICGpp(ω, params)
    η, kᵦ, kᵧ = params
    τ = η/kᵦ
    denominator = 1 + τ^2*ω^2
    numerator = ω*τ*kᵦ
    numerator/denominator
end

SLS_exp_GENERIC = RheoModelClassGENERICfromfunk(
          # Model name
          name="SLS",
          # Model parameters,
          p = [ :η, :kᵦ, :kᵧ],
          # Relaxation modulus
          G = SLSFunkNOFWNOEXPG,
          # Creep modulus
          J = SLSFunkNOFWNOEXPJ,
          # Storage modulus
          Gp = SLSFunkNOFWNOEXPGp,
         # Loss modulus
          Gpp = SLSFunkNOFWNOEXPGpp,
          # Network
          info= " "
          )

SLS_exp_GENERIC_fw = RheoModelClassGENERICfromexp(
          # Model name
          name="SLS",
          # Model parameters,
          p = [ :η, :kᵦ, :kᵧ],
          # Relaxation modulus
          G = quote
                 kᵧ + kᵦ*exp(-t*kᵦ/η)
              end,
          # Creep modulus
          J = quote
                  c₀ = 1/kᵧ
                  c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
                  τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
                  c₀ - c₁*exp(-t/τᵣ)
              end,
          # Storage modulus
          Gp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω^2*τ^2*kᵦ
                  numerator/denominator + kᵧ
               end,
         # Loss modulus
          Gpp = quote
                  τ = η/kᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω*τ*kᵦ
                  numerator/denominator
                end,
          # Network
          info= " "
          )

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Defined as functions, named tuple arguements with wrapper function for optimisation
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#
struct RheoModelFunkNT
    name::String
    free_params::Tuple
    fixed_params::NamedTuple
    G::Function
    J::Function
    Gp::Function
    Gpp::Function
end

function RheoModelFunkNT(;name::String="Custom model", free::Tuple=(), fixed::NamedTuple=NamedTuple{}(), G, J, Gp, Gpp)
   return(RheoModelFunkNT(name, free, fixed, G, J, Gp, Gpp))
end

function SLSFunk1G(t; η, kᵦ, kᵧ)
    kᵧ + kᵦ*exp(-t*kᵦ/η)
end

function SLSFunk1J(t; η, kᵦ, kᵧ)
    c₀ = 1/kᵧ
    c₁ = kᵦ/(kᵧ*(kᵧ + kᵦ))
    τᵣ = η*(kᵧ + kᵦ)/(kᵧ*kᵦ)
    c₀ - c₁*exp(-t/τᵣ)
end

function SLSFunk1Gp(ω; η, kᵦ, kᵧ)
    τ = η/kᵦ
    denominator = 1 + τ^2*ω^2
    numerator = ω^2*τ^2*kᵦ
    numerator/denominator + kᵧ
end

function SLSFunk1Gpp(ω; η, kᵦ, kᵧ)
    τ = η/kᵦ
    denominator = 1 + τ^2*ω^2
    numerator = ω*τ*kᵦ
    numerator/denominator
end

function freeze_paramsNT(orig::RheoModelFunkNT, tofix::NamedTuple)
    # some error checking, could be neatened up to show all incorrect parameters,
    # not just first incorrect one it hits
    for p in keys(tofix)
        @assert !(p in keys(orig.fixed_params)) "Parameter $p is already fixed."
        @assert (p in orig.free_params) "Parameter $p does not belong to the provided model."
    end

    # create new free symbols list without symbols about to be fixed
    newfree =  Tuple([i for i in orig.free_params if !(i in keys(tofix))])
    # combine previously fixed and newly fixed
    newfixed = merge(orig.fixed_params, tofix)

    return RheoModelFunkNT(orig.name, newfree, newfixed, orig.G, orig.J, orig.Gp, orig.Gpp)
end

function moduluswitharrayinput(t_or_ω::Vector{Float64}, numparams::Vector{Float64}, modulusfunc::Function, freeparams::Tuple, fixedargs::NamedTuple)
    freeargs = NamedTuple{freeparams}(numparams)
    args = merge(freeargs, fixedargs)

    modulusfunc.(t_or_ω; args...)
end

SLSFunk = RheoModelFunkNT(
          name="SLSFunk",
          free = (:η, :kᵦ, :kᵧ),
          G = SLSFunk1G,
          J = SLSFunk1J,
          Gp = SLSFunk1Gp,
          Gpp = SLSFunk1Gpp)

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Defined as functions, named tuple arguements with wrapper function for optimisation, WITH CLOSURES ON FUNCTIONS
SO THAT THEY BECOME SPECIALISED WITH EACH FIXING!
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#
function freeze_paramsNTclosure(orig::RheoModelFunkNT, tofix::NamedTuple)
    # some error checking, could be neatened up to show all incorrect parameters,
    # not just first incorrect one it hits
    for p in keys(tofix)
        @assert !(p in keys(orig.fixed_params)) "Parameter $p is already fixed."
        @assert (p in orig.free_params) "Parameter $p does not belong to the provided model."
    end

    # create new free symbols list without symbols about to be fixed
    newfree =  Tuple([i for i in orig.free_params if !(i in keys(tofix))])
    # combine previously fixed and newly fixed
    newfixed = merge(orig.fixed_params, tofix)

    G = let tofix=tofix; orig=orig;
            (t0; kwargs0...) -> orig.G(t0; kwargs0..., tofix...)
    end
    J = let tofix=tofix; orig=orig;
            (t0; kwargs0...) -> orig.J(t0; kwargs0..., tofix...)
    end
    Gp = let tofix=tofix; orig=orig;
            (t0; kwargs0...) -> orig.Gp(t0; kwargs0..., tofix...)
    end
    Gpp = let tofix=tofix; orig=orig;
            (t0; kwargs0...) -> orig.Gpp(t0; kwargs0..., tofix...)
    end

    return RheoModelFunkNT(orig.name, newfree, newfixed, G, J, Gp, Gpp)
end

function moduluswitharrayinputNTclosure(t_or_ω::Vector{Float64}, numparams::Vector{Float64}, modulusfunc::Function, freeparams::Tuple, fixedargs::NamedTuple)
    freeargs = NamedTuple{freeparams}(numparams)
    modulusfunc.(t_or_ω; freeargs...)
end

SLSFunkNTclosure = RheoModelFunkNT(
          name="SLSFunk",
          free = (:η, :kᵦ, :kᵧ),
          G = SLSFunk1G,
          J = SLSFunk1J,
          Gp = SLSFunk1Gp,
          Gpp = SLSFunk1Gpp)

#=
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
Benchmarking
=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~=~-~#
=#
test_time = collect(0.0:0.1:1000.0);
test_params = [1.0, 1.0, 1.0];
test_params_tuple = (η=1.0, kᵦ=1.0, kᵧ=1.0)
test_params_symbols = (:η, :kᵦ, :kᵧ)
no_fixed = NamedTuple{}()

less_named_params=(kᵦ=1.0, kᵧ=1.0)
less_params=[1.0,1.0]

# println("Expressions and Function wrappers")
# @btime SLS_exp_fw.Ga($test_time, $test_params);

# println("\nExpressions and Functions")
# @btime SLS_exp_NOFW.Ga($test_time, $test_params);

# println("\nJust functions predefined (no freezing capability)")
# @btime SLS_exp_NOFWNOEXP.Ga($test_time, $test_params);

# println("\nGeneric struct that can accept functions or functionwrappers (functions tested first then FW)")
# @btime SLS_exp_GENERIC.Ga($test_time, $test_params);
# @btime SLS_exp_GENERIC_fw.Ga($test_time, $test_params);
# println("note another option is to have two different structs and all functions accept both... They interact with RHEOS in the same way so should be fine...")

# println("\nNamed tuple approach, testing function directly...")
# @btime SLSFunk.G.($test_time; $test_params_tuple...)

# println("\nName tuple approach with calling wrapper function that would be used from optimiser")
# @btime moduluswitharrayinput($test_time, $test_params, $SLSFunk.G, $test_params_symbols, $no_fixed)

# println("\nSame as above but with one parameter 'fixed'")
# SLSFunk_1fixed = freeze_paramsNT(SLSFunk, (η=1.0,))
# @btime moduluswitharrayinput($test_time, $less_params, $SLSFunk_1fixed.G, $SLSFunk_1fixed.free_params, $SLSFunk_1fixed.fixed_params)

####################################################################################################
println("\nNAMED TUPLE APPROACH WITH CLOSURES TO SPECIALISE MODULI - BEST SOLUTION - JUST MODULI")
@btime SLSFunkNTclosure.G.($test_time; (η=1.0, kᵦ=1.0, kᵧ=1.0)...)

println("\nNAMED TUPLE APPROACH WITH CLOSURES TO SPECIALISE MODULI - BEST SOLUTION - JUST MODULI FIXED")
SLSFunkNTclosure_1fixed = freeze_paramsNTclosure(SLSFunkNTclosure, (η=1.0,))
@btime SLSFunkNTclosure_1fixed.G.($test_time; (kᵦ=1.0, kᵧ=1.0)...)

# println("\nNAMED TUPLE APPROACH WITH CLOSURES TO SPECIALISE MODULI - BEST SOLUTION - FROM ARRAY")
# @btime moduluswitharrayinput($test_time, $test_params, $SLSFunkNTclosure.G, $test_params_symbols, $no_fixed)

# println("\nNAMED TUPLE APPROACH WITH CLOSURES TO SPECIALISE MODULI - BEST SOLUTION - FROM ARRAY FIXED")
# @btime moduluswitharrayinput($test_time, $less_params, $SLSFunkNTclosure_1fixed.G, $SLSFunkNTclosure_1fixed.free_params, $SLSFunkNTclosure_1fixed.fixed_params)
######################################################################################################

println("\n===============\n[Test Complete]\n===============")

# note there is also one more option similar to named tuples but with an array wrapper function like seen on Julia discourse

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
@testset "Sets" begin
    @testset "Copy" begin
        @testset "for $S" for S in [MOI.SOS1, MOI.SOS2]
            s = S([1.0])
            s_copy = copy(s)
            s_copy.weights[1] = 2.0
            @test s.weights[1] == 1.0
        end
    end
end


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

@doc raw"""IDRange provides a min/max of an allowed range of IDs. Deprecated: use IDRange from policy API Group instead.

    IoK8sApiExtensionsV1beta1IDRange(;
        max=nothing,
        min=nothing,
    )

    - max::Int64 : max is the end of the range, inclusive.
    - min::Int64 : min is the start of the range, inclusive.
"""
mutable struct IoK8sApiExtensionsV1beta1IDRange <: SwaggerModel
    max::Any # spec type: Union{ Nothing, Int64 } # spec name: max
    min::Any # spec type: Union{ Nothing, Int64 } # spec name: min

    function IoK8sApiExtensionsV1beta1IDRange(;max=nothing, min=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1IDRange, Symbol("max"), max)
        setfield!(o, Symbol("max"), max)
        validate_property(IoK8sApiExtensionsV1beta1IDRange, Symbol("min"), min)
        setfield!(o, Symbol("min"), min)
        o
    end
end # type IoK8sApiExtensionsV1beta1IDRange

const _property_map_IoK8sApiExtensionsV1beta1IDRange = Dict{Symbol,Symbol}(Symbol("max")=>Symbol("max"), Symbol("min")=>Symbol("min"))
const _property_types_IoK8sApiExtensionsV1beta1IDRange = Dict{Symbol,String}(Symbol("max")=>"Int64", Symbol("min")=>"Int64")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1IDRange }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1IDRange))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1IDRange }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1IDRange[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1IDRange }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1IDRange[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1IDRange)
    (getproperty(o, Symbol("max")) === nothing) && (return false)
    (getproperty(o, Symbol("min")) === nothing) && (return false)
    true
end

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


text:
using DataDeps
using DataDepsGenerators
using Test
using MD5
using ReferenceTests

@testset "Figshare test" begin

    @testset "Chile Activism" begin
        registration_code = generate(Figshare(), "https://figshare.com/articles/Youth_Activism_in_Chile_from_urban_educational_inequalities_to_experiences_of_living_together_and_solidarity/6504206")

        @testset "Integration Test" begin
            eval(Meta.parse(registration_code))
            @test length(collect(readdir(datadep"Youth Activism in Chile from urban educational inequalities to experiences of living together and solidarity"))) > 0
        end

        @test_reference "references/Figshare Chile Activism.txt" registration_code
    end

    @test_reference "references/Figshare from doi.txt" generate(Figshare(), "10.1371/journal.pbio.2001414")
end


text:
# This file was generated, do not modify it. # hide
rf_mdl = RFR()
rf = machine(rf_mdl, X, y)
fit!(rf, rows=train)

rms(y[test], MLJ.predict(rf, rows=test))

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module AbstractMergers

import Flux: children, mapchildren
import Base: show

export AbstractMerger, inputname

"""
    AbstractMerger

Abstract base type for mergers.

# Interface

Any subtype should support to combine a forward stream with other streams that
can be accessed through a state dictionary via their `Splitter` name.
"""
abstract type AbstractMerger end

end # module AbstractMergers


text:
using Test
using MechGluecode
using MechanicalUnits
using MechanicalUnits: Unitfu.promote_to_derived
@import_expand(~m, s, N)
using DifferentialEquations
import MechGluecode: norm # should work without

@testset "Initial checks ArrayPartition" begin
    r0 = [1131.340, -2282.343, 6672.423]∙km
    v0 = [-5.64305, 4.30333, 2.42879]∙km/s
    Δt = 86400.0*365∙s
    μ = 398600.4418∙km³/s²
    rv0 = ArrayPartition(r0,v0)

    function fo(dy, y, μ, t)
        r = norm(y.x[1])
        dy.x[1] .= y.x[2]
        dy.x[2] .= -μ .* y.x[1] / r^3
    end
    prob = ODEProblem(fo, rv0, (0.0s, 0.035Δt), μ)
    @time sol = solve(prob, Tsit5(), dt = 0.1s);
    # Plot all variables vs time
    plot(sol)
    plot(sol, vars = (1))
    plot(sol, vars = (2))
    plot(sol, vars = (3))
    plot(sol, vars = (4))
    plot(sol, vars = (5))
    plot(sol, vars = (6))

    # Plot in phase space
    plot(sol, vars=(1,4))
    plot(sol, vars=(2,5))
    plot(sol, vars=(3,6))
    plot(sol, vars=(2,5))
    # Three distances vs time (0)
    plot(sol, vars=[(0,1), (0,2),(0,3)])
    # Three velocities vs time (0)
    plot(sol, vars=[(0,4), (0,5),(0,6)])
end
#=
@testset "Plot ODE system single pendulum" begin
    # Second order non-linear, undamped pendulum:
    # mLθ'' + mg∙sin(θ)
    # where m is mass, unit kg
    #       L is length, unit m
    #       θ is angle to vertical, strictly no unit
    #       g is gravity's acceleration
    # Acceleration
    # θ'' + k * sin(θ) = 0
    #    where k = g / L , unit(k) = s⁻²
    # Reduced to first order non-linear system by defining:
    # θ = y[1]
    # θ' = y[2]
    # It follows that
    # θ'' = y'[2]
    # and
    # y'[1] = y[2]
    # y'[2] = -k * sin(y[1])
    # To solve numerically, we define the derivative.
    # For efficiency, it's evaluated in place, i.e. dy changes at every call.
    function f!(dy, y, k, t)
        dy[1] = y[2]
        dy[2] = -k ∙ sin(y[1])
    end
    # Fixed parameter
    k = upreferred(g / 1.0m)
    # Initial conditions
    θ0 = 0.1
    θ0_vel = 0.0s⁻¹
    θ0_acc = -k *sin(θ0)

    y0 = [θ0, θ0_vel]
    y0 = ArrayPartition([θ0_vel, θ0_acc])
    tspan = (0.0, 100.0)s
    prob = ODEProblem(f!, y0, tspan, k)
    sol = solve(prob)
    # Plot all variables vs time
    plot(sol)
    # Plot in phase space
    plot(solve(prob), vars=(1,2))
end
=#


text:
# makes calling C functions a bit easier
macro cpx_ccall(func, args...)
    f = "CPX$(func)"
    if is_unix()
        return quote
            ccall(($f,libcplex), $(args...))
        end
    end
    if is_windows()
        return quote
            ccall(($f,libcplex), stdcall, $(args...))
        end
    end
end

macro cpx_ccall_intercept(model, func, args...)
    f = "CPX$(func)"
    quote
        ccall(:jl_exit_on_sigint, Void, (Cint,), convert(Cint,0))
        ret = try
            $(@static is_windows() ? :(ccall(($f,libcplex), stdcall, $(args...))) : :(ccall(($f,libcplex), $(args...))) )
        catch ex
            println("Caught exception")
            if !isinteractive()
                ccall(:jl_exit_on_sigint, Void, (Cint,), convert(Cint,1))
            end
            if isa(ex, InterruptException)
                model.terminator[1] = 1
            end
            rethrow(ex)
        end
        if !isinteractive()
            ccall(:jl_exit_on_sigint, Void, (Cint,), convert(Cint,1))
        end
        ret
    end
end

typealias GChars Union{Cchar, Char}
typealias IVec Vector{Cint}
typealias FVec Vector{Cdouble}
typealias CVec Vector{Cchar}
typealias CoeffMat Union{Matrix{Cdouble}, SparseMatrixCSC{Cdouble}}
typealias Bounds{T<:Real} Union{T, Vector{T}}

typealias GCharOrVec Union{Cchar, Char, Vector{Cchar}, Vector{Char}}

# empty vector & matrix (for the purpose of supplying default arguments)
const emptyfvec = Array(Float64, 0)
const emptyfmat = Array(Float64, 0, 0)

cchar(c::Cchar) = c
cchar(c::Char) = convert(Cchar, c)

ivec(v::IVec) = v
fvec(v::FVec) = v
cvec(v::CVec) = v

ivec(v::Vector) = convert(IVec, v)
fvec(v::Vector) = convert(FVec, v)
cvec(v::Vector) = convert(CVec, v)

# cvecx(v, n) and fvecx(v, n)
# converts v into a vector of Cchar or Float64 of length n,
# where v can be either a scalar or a vector of length n.

_chklen(v, n::Integer) = (length(v) == n || error("Inconsistent argument dimensions."))

cvecx(c::GChars, n::Integer) = fill(cchar(c), n)
cvecx(c::Vector{Cchar}, n::Integer) = (_chklen(c, n); c)
cvecx(c::Vector{Char}, n::Integer) = (_chklen(c, n); convert(Vector{Cchar}, c))

fvecx(v::Real, n::Integer) = fill(Float64(v), n)
fvecx(v::Vector{Float64}, n::Integer) = (_chklen(v, n); v)
fvecx{T<:Real}(v::Vector{T}, n::Integer) = (_chklen(v, n); convert(Vector{Float64}, v))


text:
"""
A package for dealing with the estimation of (co)variance matrices.

    1. Define CovarianceEstimator(s) which are robust to heteroskedasticity, time-dependence, amd spatial dependence (cluster)
    2. Optionally, interface with GLM.jl to provide estimates of coefficients' standard errors which are robust

These use cases utilize different parts of this package, make sure you read the
documentation.
"""

module CovarianceMatrices

using CategoricalArrays
using LinearAlgebra
using Requires: @require
using StatsBase: CovarianceEstimator
using StatsModels
using Statistics
include("types.jl")
include("HAC.jl")
include("HC.jl")
include("CRHC.jl")
include("VARHAC.jl")
include("lrvar.jl")
include("CovarianceMatrix.jl")
# using GLM
# include("glm.jl")

function __init__()
    @require GLM="38e38edf-8417-5370-95a0-9cbb8c7f171a" include("glm.jl")
end

export Andrews,
       BartlettKernel,
       CRHC0,
       CRHC1,
       CRHC2,
       CRHC3,
       CovarianceMatrix,
       HC0,
       HC1,
       HC2,
       HC3,
       HC4,
       HC4m,
       HC5,
       NeweyWest,
       ParzenKernel,
       QuadraticSpectralKernel,
       TruncatedKernel,
       TukeyHanningKernel,
       lrvar,
       lrvarmatrix,
       optimalbandwidth,
       #stderror,
       #vcov,
       vcovmatrix

end # module


text:
function L1TransformProvider(arg0::CelestialBody, arg1::CelestialBody)
    return L1TransformProvider((CelestialBody, CelestialBody), arg0, arg1)
end

function get_transform(obj::L1TransformProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::L1TransformProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
uni_rand = (r) -> rand(Uniform(r...));

function rand_hermitian(x_range, y_range, scale)
    a = uni_rand(x_range)
    b = uni_rand(y_range)
    m = min(a, b)
    diag_range = -m:scale:m
    c = rand(diag_range)
    [a c; c b]
end

function collapse(p, scale)
    p = round(p)
    i = convert(Int, p)
    i - i % scale
end

function generate_clients(
    counts, map_size, min_var=1, max_var=2, scale=10)
    x = 1:map_size;
    y = 1:map_size;

    var_inter = (x) -> [min_var  max_var] .* x;

    clients = zeros(map_size, map_size)
    centers = rand(1:map_size, (length(counts), 2))

    for (c, μ) in zip(counts, eachrow(centers))
        Σ = rand_hermitian(var_inter(c), var_inter(c), 0.1)
        nd = MvNormal(μ, Σ)
        points = collapse.(rand(nd, c), scale)
        for (p1, p2) in eachcol(points)
            if p1 ∈ x && p2 ∈ y
                clients[p1, p2] += 1
            end
        end

    end

    clients
end


text:
#Reference
#岑松, 龙志飞, 龙驭球. 对转角场和剪应变场进行合理插值的厚薄板通用四边形单元[J]. 工程力学, 1999,16(4): 1―15.
function K_TMQ(elm::Quad)::Matrix{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    G=elm.material.G
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    x₄,y₄,z₄=T*(elm.node4.loc-center)
    b₁=y₂-y₃
    b₂=y₃-y₄
    b₃=y₄-y₁
    b₄=y₁-y₂

    c₁=x₃-x₂
    c₂=x₄-x₃
    c₃=x₁-x₄
    c₄=x₂-x₁

    d₁=norm(elm.node3.loc-elm.node2.loc)
    d₂=norm(elm.node4.loc-elm.node3.loc)
    d₃=norm(elm.node1.loc-elm.node4.loc)
    d₄=norm(elm.node2.loc-elm.node1.loc)

    δ₁=(t/d₁)^2/(5/6*(1-ν₀)+2*(t/d₁)^2)
    δ₂=(t/d₂)^2/(5/6*(1-ν₀)+2*(t/d₂)^2)
    δ₃=(t/d₃)^2/(5/6*(1-ν₀)+2*(t/d₃)^2)
    δ₄=(t/d₄)^2/(5/6*(1-ν₀)+2*(t/d₄)^2)

    Γ=[0 0 0 -2δ₁ -c₁*δ₁ b₁*δ₁ 2δ₁ -c₁*δ₁ b₁*δ₁ 0 0 0;
       0 0 0 0 0 0 -2δ₂ -c₂*δ₂ b₂*δ₂ 2δ₂ -c₂*δ₂ b₂*δ₂;
       2δ₃ -c₃*δ₃ b₃*δ₃ 0 0 0 0 0 0 -2δ₃ -c₃*δ₃ b₃*δ₃;
      -2δ₄ -c₄*δ₄ b₄*δ₄ 2δ₄ -c₄*δ₄ b₄*δ₄ 0 0 0 0 0 0]

    Xₛ=[0 0 -b₄/(b₃*c₄-b₄*c₃) b₃/(b₃*c₄-b₄*c₃);
       b₄/(b₄*c₁-b₁*c₄) 0 0 -b₁/(b₄*c₁-b₁*c₄);
       -b₂/(b₁*c₂-b₂*c₁) b₁/(b₁*c₂-b₂*c₁) 0 0;
       0 -b₃/(b₂*c₃-b₃*c₂) b₂/(b₂*c₃-b₃*c₂) 0]

    Yₛ=[0 0 -c₄/(b₃*c₄-b₄*c₃) c₃/(b₃*c₄-b₄*c₃);
       c₄/(b₄*c₁-b₁*c₄) 0 0 -c₁/(b₄*c₁-b₁*c₄);
       -c₂/(b₁*c₂-b₂*c₁) c₁/(b₁*c₂-b₂*c₁) 0 0;
       0 -c₃/(b₂*c₃-b₃*c₂) c₂/(b₂*c₃-b₃*c₂) 0]

    α=[0 0 0 -3c₁/2/d₁^2*(1-2δ₁) 1/2/d₁^2*(b₁^2-c₁^2/2*(1-6δ₁)) 3b₁*c₁/4/d₁^2*(1-2δ₁) 3c₁/2/d₁^2*(1-2δ₁) 1/2/d₁^2*(b₁^2-c₁^2/2*(1-6δ₁)) 3b₁*c₁/4/d₁^2*(1-2δ₁) 0 0 0;
       0 0 0 0 0 0 -3c₂/2/d₂^2*(1-2δ₂) 1/2/d₂^2*(b₂^2-c₂^2/2*(1-6δ₂)) 3b₂*c₂/4/d₂^2*(1-2δ₂) 3c₂/2/d₂^2*(1-2δ₂) 1/2/d₂^2*(b₂^2-c₂^2/2*(1-6δ₂)) 3b₂*c₂/4/d₂^2*(1-2δ₂);
       3c₃/2/d₃^2*(1-2δ₃) 1/2/d₃^2*(b₃^2-c₃^2/2*(1-6δ₃)) 3b₃*c₃/4/d₃^2*(1-2δ₃) 0 0 0 0 0 0 -3c₃/2/d₃^2*(1-2δ₃) 1/2/d₃^2*(b₃^2-c₃^2/2*(1-6δ₃)) 3b₃*c₃/4/d₃^2*(1-2δ₃);
       -3c₄/2/d₄^2*(1-2δ₄) 1/2/d₄^2*(b₄^2-c₄^2/2*(1-6δ₄)) 3b₄*c₄/4/d₄^2*(1-2δ₄) 3c₄/2/d₄^2*(1-2δ₄) 1/2/d₄^2*(b₄^2-c₄^2/2*(1-6δ₄)) 3b₄*c₄/4/d₄^2*(1-2δ₄) 0 0 0 0 0 0]

    β=[0 0 0 3b₁/2/d₁^2*(1-2δ₁) 3b₁*c₁/4/d₁^2*(1-2δ₁) 1/2/d₁^2*(c₁^2-b₁^2/2*(1-6δ₁)) -3b₁/2/d₁^2*(1-2δ₁) 3b₁*c₁/4/d₁^2*(1-2δ₁) 1/2/d₁^2*(c₁^2-b₁^2/2*(1-6δ₁)) 0 0 0;
       0 0 0 0 0 0 3b₂/2/d₂^2*(1-2δ₂) 3b₂*c₂/4/d₂^2*(1-2δ₂) 1/2/d₂^2*(c₂^2-b₂^2/2*(1-6δ₂)) -3b₂/2/d₂^2*(1-2δ₂) 3b₂*c₂/4/d₂^2*(1-2δ₂) 1/2/d₂^2*(c₂^2-b₂^2/2*(1-6δ₂));
       -3b₃/2/d₃^2*(1-2δ₃) 3b₃*c₃/4/d₃^2*(1-2δ₃) 1/2/d₃^2*(c₃^2-b₃^2/2*(1-6δ₃)) 0 0 0 0 0 0 3b₃/2/d₃^2*(1-2δ₃) 3b₃*c₃/4/d₃^2*(1-2δ₃) 1/2/d₃^2*(c₃^2-b₃^2/2*(1-6δ₃));
       3b₄/2/d₄^2*(1-2δ₄) 3b₄*c₄/4/d₄^2*(1-2δ₄) 1/2/d₄^2*(c₄^2-b₄^2/2*(1-6δ₄)) -3b₄/2/d₄^2*(1-2δ₄) 3b₄*c₄/4/d₄^2*(1-2δ₄) 1/2/d₄^2*(c₄^2-b₄^2/2*(1-6δ₄)) 0 0 0 0 0 0]

    D₀=E₀*t^3/12/(1-ν₀^2)
    D=D₀*[1 ν₀ 0;
          ν₀ 1 0;
          0  0 (1-ν₀)/2]
    k=5/6
    C=k*G*t*[1 0;0 1]

    K=Matrix{Float64}(undef,12,12)
    J=Matrix{Float64}(undef,2,2)
    function BtDB(x)
        ξ,η=x[1],x[2]
        J[1,1]=-0.25*(1 - η)*x₁ + 0.25*(1 - η)*x₂ + 0.25*(1 + η)*x₃ - 0.25*(1 + η)*x₄
        J[1,2]=-0.25*(1 - η)*y₁ + 0.25*(1 - η)*y₂ + 0.25*(1 + η)*y₃ - 0.25*(1 + η)*y₄
        J[2,1]=-0.25*(1 - ξ)*x₁ + 0.25*(1 - ξ)*x₄ - 0.25*(1 + ξ)*x₂ + 0.25*(1 + ξ)*x₃
        J[2,2]=-0.25*(1 - ξ)*y₁ + 0.25*(1 - ξ)*y₄ - 0.25*(1 + ξ)*y₂ + 0.25*(1 + ξ)*y₃

        N₁⁰=1/4*(1-ξ)*(1-η)
        N₂⁰=1/4*(1+ξ)*(1-η)
        N₃⁰=1/4*(1+ξ)*(1+η)
        N₄⁰=1/4*(1-ξ)*(1+η)

        Nₛ⁰=[N₁⁰ N₂⁰ N₃⁰ N₄⁰]
        Bs=[Nₛ⁰*Xₛ*Γ;
            Nₛ⁰*Yₛ*Γ]

        dN₁dξ = -0.25*(1 - ξ)*(1 - η) + 0.25*(1 + η + ξ)*(1 - η)
        dN₁dη = -0.25*(1 - ξ)*(1 - η) + 0.25*(1 + η + ξ)*(1 - ξ)
        dN₂dξ = 0.25*(1 + ξ)*(1 - η) - 0.25*(1 + η - ξ)*(1 - η)
        dN₂dη = -0.25*(1 + ξ)*(1 - η) + 0.25*(1 + η - ξ)*(1 + ξ)
        dN₃dξ = 0.25*(1 + ξ)*(1 + η) - 0.25*(1 - η - ξ)*(1 + η)
        dN₃dη = 0.25*(1 + ξ)*(1 + η) - 0.25*(1 - η - ξ)*(1 + ξ)
        dN₄dξ = -0.25*(1 - ξ)*(1 + η) + 0.25*(1 - η + ξ)*(1 + η)
        dN₄dη = 0.25*(1 - ξ)*(1 + η) - 0.25*(1 - η + ξ)*(1 - ξ)
        dN₅dξ = 0.5*(1 - η^2)
        dN₅dη = -1.0*η*(1 + ξ)
        dN₆dξ = -1.0*ξ*(1 + η)
        dN₆dη = 0.5*(1 - ξ^2)
        dN₇dξ = -0.5*(1 - η^2)
        dN₇dη = -1.0*η*(1 - ξ)
        dN₈dξ = -1.0*ξ*(1 - η)
        dN₈dη = -0.5*(1 - ξ^2)
        # J=Matrix(1.0I,2,2)
        dN₁dx,dN₁dy=inv(J)*[dN₁dξ,dN₁dη]
        dN₂dx,dN₂dy=inv(J)*[dN₂dξ,dN₂dη]
        dN₃dx,dN₃dy=inv(J)*[dN₃dξ,dN₃dη]
        dN₄dx,dN₄dy=inv(J)*[dN₄dξ,dN₄dη]
        dN₅dx,dN₅dy=inv(J)*[dN₅dξ,dN₅dη]
        dN₆dx,dN₆dy=inv(J)*[dN₆dξ,dN₆dη]
        dN₇dx,dN₇dy=inv(J)*[dN₇dξ,dN₇dη]
        dN₈dx,dN₈dy=inv(J)*[dN₈dξ,dN₈dη]

        H₀=[0 dN₁dx     0 0 dN₂dx     0 0     dN₃dx 0 0     dN₄dx 0;
            0     0 dN₁dy 0     0 dN₂dy 0     0 dN₃dy 0     0 dN₄dy;
            0 dN₁dy dN₁dx 0 dN₂dy dN₂dx 0 dN₃dy dN₃dx 0 dN₄dy dN₄dx]

        H₁=[dN₅dx dN₆dx dN₇dx dN₈dx;
            0 0 0 0;
            dN₅dy dN₆dy dN₇dy dN₈dy]

        H₂=[0 0 0 0;
            dN₅dy dN₆dy dN₇dy dN₈dy;
            dN₅dx dN₆dx dN₇dx dN₈dx]

        Bb=-(H₀+H₁*α+H₂*β)

        Kb=transpose(Bb)*D*Bb
        Ks=transpose(Bs)*C*Bs
        K=Kb+Ks

        return K*det(J)
    end
    Kb=hcubature(BtDB,[-1,-1],[1,1])[1]
    #left-hand system to right-hand system
    I=1:12
    J=[1,3,2,4,6,5,7,9,8,10,12,11]
    L=sparse(I,J,1.,12,12)
    for i in [1,3,4,6,7,9,10,12]
        for j in [2,5,8,11]
            Kb[i,j]=Kb[j,i]=-Kb[i,j]
        end
    end
    K=L'*Kb*L
    #12x12 to 24x24
    I=1:12
    J=[3,4,5,9,10,11,15,16,17,21,22,23]
    L=sparse(I,J,1.,12,24)
    return L'*K*L
end

function K2_TMQ(elm::Quad)::Matrix{Float64}
    spzeros(24,24)
end