text:
@testset "Testing ResDense" begin
	@test size(FluxExtensions.ResDense(3,3,NNlib.relu)(randn(3,4))) == (3,4)
	@test size(FluxExtensions.ResDense(5,3,NNlib.relu)(randn(5,4))) == (3,4)
	@test length(Flux.params(FluxExtensions.ResDense(5,3,NNlib.relu))) == 4
	@test length(Flux.params(FluxExtensions.ResDense(3,3,NNlib.relu))) == 3
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Shapes (Modia3D/Shapes/_module.jl)
#

import JSON

# Shadow definitions
@enum ShadowMask NoShadows=0 CastsShadows=1 ReceivesShadows=2 CastsAndReceivesShadows=3

# Material for visual elements
"""
    VisualMaterial(; color="Green",
                     wireframe=false, transparency=0.5,
                     reflectslight=true, shininess=0.7,
                     shadowMask=CastsAndReceivesShadows))

Return a material object that defines attributes for the visualization of an Object3D
that has visual or solid properties.

# Arguments
- `color` defines the material color as a String or RGB values. In case of Strings the color
          is constructed using [Colors.jl](https://github.com/JuliaGraphics/Colors.jl). So
          [named colors](https://juliagraphics.github.io/Colors.jl/stable/namedcolors/) and
          many other kinds of [color specifications](https://juliagraphics.github.io/Colors.jl/stable/constructionandconversion/#Color-Parsing)
          are available. In addition, colors can be directly defined by a vector of RGB values.
          E.g. `color="MediumVioletRed"` or `color=[255, 0, 0]`.
- `wireframe`: = false, if solid, otherwise wireframe representation.
- `transparency`: = 0.0 (opaque) ... 1.0 (fully transparent).
- `reflectslight`: = true if it reflects light and false, if it does not reflect light.
- `shininess`: = 0.0 (matte surface) ... 1.0 (very shiny surface).
- `shadowMask`: defines whether or not an object casts or receives shadows. Possible values:
   NoShadows, CastsShadows, ReceivesShadows, CastsAndReceivesShadows.
"""
mutable struct VisualMaterial
  color::MVector{3,Cint}  # RGB color
  wireframe::Bool         # if true, the wireframe of object is displayed
  transparency::Float64   # [0,1]... 0.0 is opaque, 1.0 is transparent
  reflectslight::Bool     # if true, specular Highlights are active
  shininess::Float64      # [0,1]... 0.0 is matte surface, 1.0 is very shiny
  shadowMask::ShadowMask  # defines whether or not an object casts or receives shadows 0: No Shadows, 1: Casts Shadows, 2: Receives Shadows, 3: Casts and Recieves Shadows

  function VisualMaterial(;color=defaultColor(), wireframe=false, transparency=0.0, reflectslight=true, shininess=0.7, shadowMask=CastsAndReceivesShadows)
    color2 = rgb(color)
    @assert(color2[1]>=0 && color2[1]<=255)
    @assert(color2[2]>=0 && color2[2]<=255)
    @assert(color2[3]>=0 && color2[3]<=255)
    @assert(transparency >= 0.0 && transparency <= 1.0)
    @assert(shininess    >= 0.0 && shininess    <= 1.0)
  new(color2, wireframe, transparency, reflectslight, shininess, shadowMask)
  end
end

"""
    readDictOfStructsFromJSONVisualMaterial(fileName, StructType)

Read a JSON file from `fileName` and return a `Dict{String, StructType}` dictionary.
`StructType` must be a mutable struct type with a constructor `StructType()`.
"""

function readDictOfStructsFromJSONVisualMaterial(fileName, StructType)
    dict1 = JSON.parsefile(fileName)
    palette = Dict{String, StructType}()
    for (key1,value1) in dict1
        obj = StructType()
        for (key2,value2) in value1
            if Symbol(key2) == Symbol("color")
              value2 = Shapes.rgb(value2)
            end
            setfield!(obj, Symbol(key2), value2)
        end
        palette[key1] = obj
    end
    return palette
end

"""
    readVisualMaterialFromJSON(fileName)

Read a JSON file consisting of a dictionary of VisualMaterial instances from `fileName` and
return a `Dict{String, VisualMaterial}` dictionary.
"""
readDictOfStructsFromJSONVisualMaterial(fileName::AbstractString) = readDictOfStructsFromJSONVisualMaterial(fileName, VisualMaterial)

"""
    const visualMaterialPalette

Dictionary of visual material data, see [`Modia3D.VisualMaterial`](@ref)
"""
const visualMaterialPalette = readDictOfStructsFromJSONVisualMaterial( joinpath(Modia3D.path, "palettes", "visualMaterials.json") )

ext: jl


text:
include("weeds.jl")

#run the inverse analysis repeatedly for a set of parameters
sigma, lambda, mu, deltah, sqrtnumobs, alpha, beta, gamma, iterations = 1.0, 50.0, 0.0, 5.0, 16, 1.0, 1.0, 1.0, 100#here are the parameters we will run the inverse analysis with
@time results = pmap(x->inversemodel(sigma, lambda, mu, deltah, sqrtnumobs, alpha, beta, gamma, iterations), 1:100)#run the inverse analysis 100 times (in parallel if julia is run like, e.g., "julia -p 8" where 8 processors will be used)

#postprocess the results
eig_residuals = map(x->x[1], results)
eig_rmse = sum(map(x->x .^ 2, eig_residuals)) / length(eig_residuals)
logK_rmse = map(x->x[2], results)
of_vals = map(x->x[3], results)

#make a couple plots based on the results
fig, ax = PyPlot.subplots()
ax.semilogy(eig_rmse, ".")
ax.set_xlabel("Eigenvector Index")
ax.set_ylabel("Eigenvector Coefficient RMSE")
PyPlot.savefig("eigenvector_coefficient_fit.pdf")
display(fig)
println()
PyPlot.close(fig)

fig, ax = PyPlot.subplots()
for i = 1:length(logK_rmse)
	ax.loglog(of_vals[i], logK_rmse[i], ".", alpha=0.1)
end
ax.set_xlabel("Objective Function")
ax.set_ylabel("L2 Error in logK")
PyPlot.savefig("logK_vs_objective_function.pdf")
display(fig)
println()
PyPlot.close(fig)


text:
## Caches

function CachedOperator(io::InterlaceOperator{T,1};padding::Bool=false) where T
    ds=domainspace(io)
    rs=rangespace(io)

    ind=findall(op->isinf(size(op,1)), io.ops)
    if length(ind) ≠ 1  || !isbanded(io.ops[ind[1]])  # is almost banded
        return default_CachedOperator(io;padding=padding)
    end
    i=ind[1]
    bo=io.ops[i]
    lin,uin=bandwidths(bo)

    # calculate number of rows interlaced
    # each each row in the rangespace increases the lower bandwidth by 1
    nds=0
    md=0
    for k=1:length(io.ops)
        if k ≠ i
            d=dimension(rs[k])
            nds+=d
            md=max(md,d)
        end
    end

    isend=true
    for k=i+1:length(io.ops)
        if dimension(rs[k]) == md
            isend=false
        end
    end

    numoprows=isend ? md-1 : md
    n=nds+numoprows

    (l,u) = (max(lin+nds,n-1),max(0,uin+1-ind[1]))

    # add extra rows for QR
    if padding
        u+=l
    end

    ret=AlmostBandedMatrix(Zeros{T}(n,n+u),(l,u),nds)

    # populate the finite rows
    jr=1:n+u
    ioM=io[1:n,jr]

    bcrow=1
    oprow=0
    for k=1:n
        K,J=io.rangeinterlacer[k]

        if K ≠ i
            # fill the fill matrix
            ret.fill.V[:,bcrow] = Matrix(view(io.ops[K],J:J,jr))
            ret.fill.U[k,bcrow] = 1
            bcrow += 1
        else
            oprow+=1
        end

        for j=rowrange(ret.bands,k)
            ret[k,j] = ioM[k,j]
        end
    end

    CachedOperator(io,ret,(n,n+u),ds,rs,(l,∞))
end

function CachedOperator(io::InterlaceOperator{T,2};padding::Bool=false) where T
    ds=domainspace(io)
    rs=rangespace(io)
    di=io.domaininterlacer
    ri=io.rangeinterlacer
    ddims=dimensions(di.iterator)
    rdims=dimensions(ri.iterator)

    # we are only almost banded if every operator is either finite
    # range or banded, and if the # of  ∞ spaces is the same
    # for the domain and range
    isab=all(op->isfinite(size(op,1)) || isbanded(op),io.ops) &&
        count(isinf,ddims) == count(isinf,rdims)

    if !isab
        return default_CachedOperator(io;padding=padding)
    end

    # these are the bandwidths if we only had ∞-dimensional operators
    d∞=findall(isinf,collect(ddims))
    r∞=findall(isinf,collect(rdims))
    p=length(d∞)

    # we only support block size 1 for now
    for k in d∞
        bl = blocklengths(ds[k])
        if !(bl isa AbstractFill) || getindex_value(bl) ≠ 1
            return default_CachedOperator(io;padding=padding)
        end
    end
    for k in r∞
        bl = blocklengths(rs[k])
        if !(bl isa AbstractFill) || getindex_value(bl) ≠ 1
            return default_CachedOperator(io;padding=padding)
        end
    end

    l∞,u∞ = 0,0
    for k=1:p,j=1:p
        l∞=max(l∞,p*bandwidth(io.ops[r∞[k],d∞[j]],1)+k-j)
    end
    for k=1:p,j=1:p
        u∞=max(u∞,p*bandwidth(io.ops[r∞[k],d∞[j]],2)+j-k)
    end

    # now we move everything by the finite rank
    ncols=mapreduce(d->isfinite(d) ? d : 0,+,ddims)
    nbcs=mapreduce(d->isfinite(d) ? d : 0,+,rdims)
    shft=ncols-nbcs
    l∞,u∞=l∞-shft,u∞+shft

    # iterate through finite rows to find worst case bandwidth
    l,u=l∞,u∞
        for k=1+nbcs+p,j=1:ncols+p
            N,n=ri[k]
            M,m=di[j]
            l=max(l,bandwidth(io.ops[N,M],1)+m-n+k-j)
            u=max(u,bandwidth(io.ops[N,M],2)+n-m+j-k)
        end

    # add extra rows for QR
    if padding
        u += l
    end

    n=1+nbcs+p
    ret=AlmostBandedMatrix(Zeros{T}(n,n+u),(l,u),nbcs)

    # populate entries and fill functionals
    bcrow=1
    oprow=0
    jr=1:n+u
    for k=1:n
        K,J=io.rangeinterlacer[k]

        if isfinite(rdims[K] )
            # fill the fill matrix
            ret.fill.V[:,bcrow] = Matrix(view(io,k:k,jr))
            ret.fill.U[k,bcrow] = 1
            bcrow += 1
        else
            oprow+=1
        end

        for j=rowrange(ret.bands,k)
            ret[k,j] = io[k,j]
        end
    end

    CachedOperator(io,ret,(n,n+u),ds,rs,(l,∞))
end

# Grow cached interlace operator

function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T},
                              InterlaceOperator{T,1,DS,RS,DI,RI,BI}},
           n::Integer,::Colon) where {T<:Number,DS,RS,DI,RI,BI}
    if n ≤ co.datasize[1]
        return co
    end

    (l,u)=bandwidths(co.data.bands)
    pad!(co.data,n,n+u)

    r=rank(co.data.fill)
    ind=findfirst(op->isinf(size(op,1)),co.op.ops)

    k=1
    for (K,J) in co.op.rangeinterlacer
        if K ≠ ind
            co.data.fill.V[co.datasize[2]:end,k] = co.op.ops[K][J,co.datasize[2]:n+u]
            k += 1
            if k > r
                break
            end
        end
    end

    kr=co.datasize[1]+1:n
    jr=max(1,kr[1]-l):n+u
    BLAS.axpy!(1.0,view(co.op.ops[ind],kr .- r,jr),
                    view(co.data.bands,kr,jr))

    co.datasize=(n,n+u)
    co
end

function resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T},
                              InterlaceOperator{T,2,DS,RS,DI,RI,BI}},
           n::Integer,::Colon) where {T<:Number,DS,RS,DI,RI,BI}
    if n ≤ co.datasize[1]
        return co
    end

    io=co.op
    ds=domainspace(io)
    rs=rangespace(io)
    di=io.domaininterlacer
    ri=io.rangeinterlacer
    ddims=dimensions(di.iterator)
    rdims=dimensions(ri.iterator)

    d∞=findall(isinf,collect(ddims))
    r∞=findall(isinf,collect(rdims))
    p=length(d∞)

    (l,u)=bandwidths(co.data.bands)
    pad!(co.data,n,n+u)
    co.data
    # r is number of extra rows, ncols is number of extra columns
    r=rank(co.data.fill)
    ncols=mapreduce(d->isfinite(d) ? d : 0,+,ddims)

    # fill rows
    K=k=1
    while k ≤ r
        if isfinite(dimension(rs[ri[K][1]]))
            co.data.fill.V[co.datasize[2]:end,k] = co.op[K,co.datasize[2]:n+u]
            k += 1
        end
        K += 1
    end

    kr=co.datasize[1]+1:n
    jr=max(ncols+1,kr[1]-l):n+u
    io∞=InterlaceOperator(io.ops[r∞,d∞])

    BLAS.axpy!(1.0,view(io∞,kr.-r,jr.-ncols),view(co.data.bands,kr,jr))

    co.datasize=(n,n+u)
    co
end

resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T},
                              InterlaceOperator{T,1,DS,RS,DI,RI,BI}},
n::Integer,m::Integer) where {T<:Number,DS,RS,DI,RI,BI} = resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:)

resizedata!(co::CachedOperator{T,AlmostBandedMatrix{T},
                               InterlaceOperator{T,2,DS,RS,DI,RI,BI}},
n::Integer,m::Integer) where {T<:Number,DS,RS,DI,RI,BI} = resizedata!(co,max(n,m+bandwidth(co.data.bands,1)),:)

##
# These give highly optimized routines for delaying with Cached
#

## QR

function QROperator(R::CachedOperator{T,AlmostBandedMatrix{T}}) where T
    M = R.data.bands.l+1   # number of diag+subdiagonal bands
    H = Matrix{T}(undef,M,100)
    QROperator(R,H,0)
end

function resizedata!(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T},
                                                  MM,DS,RS,BI}},
         ::Colon,col) where {T,MM,DS,RS,BI}
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R_cache
    W=QR.H

    R=MO.data.bands
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    F=MO.data.fill.U

    for k=QR.ncols+1:col
        W[:,k] = view(R.data,R.u+1:R.u+R.l+1,k) # diagonal and below
        wp=view(W,:,k)
        W[1,k]+= flipsign(norm(wp),W[1,k])
        normalize!(wp)

        # scale banded entries
        for j=k:k+R.u
            dind=R.u+1+k-j
            v=view(R.data,dind:dind+M-1,j)
            dt=dot(wp,v)
            LinearAlgebra.axpy!(-2*dt,wp,v)
        end

        # scale banded/filled entries
        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=view(R.data,1:M-p,j)  # shift down each time
            wp2=view(wp,p+1:M)
            dt=dot(wp2,v)
            for ℓ=k:k+p-1
                @inbounds dt=muladd(conj(W[ℓ-k+1,k]),
                                    unsafe_getindex(MO.data.fill,ℓ,j),dt)
            end
            LinearAlgebra.axpy!(-2*dt,wp2,v)
        end

        # scale filled entries

        for j=1:size(F,2)
            v=view(F,k:k+M-1,j) # the k,jth entry of F
            dt=dot(wp,v)
            LinearAlgebra.axpy!(-2*dt,wp,v)
        end
    end
    QR.ncols=col
    QR
end

# BLAS versions, requires BlasFloat

function resizedata!(QR::QROperator{CachedOperator{T,AlmostBandedMatrix{T},
                                       MM,DS,RS,BI}},
::Colon,col) where {T<:BlasFloat,MM,DS,RS,BI}
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R_cache
    W=QR.H

    R=MO.data.bands
    M=R.l+1   # number of diag+subdiagonal bands

    if col+M-1 ≥ MO.datasize[1]
        resizedata!(MO,(col+M-1)+100,:)  # double the last rows
    end

    if col > size(W,2)
        W=QR.H=unsafe_resize!(W,:,2col)
    end

    F=MO.data.fill.U

    f=pointer(F)
    m,n=size(R)
    w=pointer(W)
    r=pointer(R.data)
    sz=sizeof(T)
    st=stride(R.data,2)
    stw=stride(W,2)

    for k=QR.ncols+1:col
        v=r+sz*(R.u + (k-1)*st)    # diagonal entry
        wp=w+stw*sz*(k-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W[1,k]+= flipsign(BLAS.nrm2(M,wp,1),W[1,k])
        normalize!(M,wp)

        for j=k:k+R.u
            v=r+sz*(R.u + (k-1)*st + (j-k)*(st-1))
            dt = BandedMatrices.dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end

        for j=k+R.u+1:k+R.u+M-1
            p=j-k-R.u
            v=r+sz*((j-1)*st)  # shift down each time
            dt = BandedMatrices.dot(M-p,wp+p*sz,1,v,1)
            for ℓ=k:k+p-1
                @inbounds dt=muladd(conj(W[ℓ-k+1,k]),
                                    unsafe_getindex(MO.data.fill,ℓ,j),dt)
            end
            BLAS.axpy!(M-p,-2*dt,wp+p*sz,1,v,1)
        end

        fp=f+(k-1)*sz
        fst=stride(F,2)
        for j=1:size(F,2)
            v=fp+fst*(j-1)*sz   # the k,jth entry of F
            dt = BandedMatrices.dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end
    end
    QR.ncols=col
    QR
end

## back substitution
# loop to avoid ambiguity with AbstractTRiangular
for ArrTyp in (:AbstractVector, :AbstractMatrix)
    @eval function ldiv!(U::UpperTriangular{T, SubArray{T, 2, AlmostBandedMatrix{T}, Tuple{UnitRange{Int}, UnitRange{Int}}, false}},
                       u::$ArrTyp{T}) where T
        n = size(u,1)
        n == size(U,1) || throw(DimensionMismatch())

        V = parent(U)
        @assert parentindices(V)[1][1] == 1
        @assert parentindices(V)[2][1] == 1

        B = parent(V)

        A = B.bands
        F = B.fill
        b=bandwidth(A,2)
        nbc = rank(B.fill)

        pk = zeros(T,nbc)

        for c=1:size(u,2)
            fill!(pk,zero(T))

            # before we get to filled rows
            for k=n:-1:max(1,n-b)
                @simd for j=k+1:n
                    @inbounds u[k,c] = muladd(-A.data[k-j+A.u+1,j],u[j,c],u[k,c])
                end

                @inbounds u[k,c] /= A.data[A.u+1,k]
            end

           #filled rows
            for k=n-b-1:-1:1
                @simd for j=1:nbc
                    @inbounds pk[j] = muladd(u[k+b+1,c],F.V[k+b+1,j],pk[j])
                end

                @simd for j=k+1:k+b
                    @inbounds u[k,c]=muladd(-A.data[k-j+A.u+1,j],u[j,c],u[k,c])
                end

                @simd for j=1:nbc
                    @inbounds u[k,c] = muladd(-F.U[k,j],pk[j],u[k,c])
                end

                @inbounds u[k,c] /= A.data[A.u+1,k]
            end
        end
        u
    end
end


text:
# reference from https://www.bioinfor.com/amino-acid/
const RESNAME2MASS = Dict{String, Float32}(
"ALA" =>  71.03711,
"ARG" => 156.10111,
"ASN" => 114.04293,
"ASP" => 115.02694,
"CYS" => 103.00919,
"GLU" => 129.04259,
"GLN" => 128.05858,
"GLY" =>  57.02146,
"HIS" => 137.05895,
"ILE" => 113.08406,
"LEU" => 113.08406,
"LYS" => 128.09496,
"MET" => 131.04049,
"PHE" => 147.06841,
"PRO" =>  97.05276,
"SER" =>  87.03203,
"THR" => 101.04768,
"TRP" => 186.07931,
"TYR" => 163.06333,
"VAL" =>  99.06841,
"HOH" =>  18.01528
                                    )

const ATOMNAME2MASS = Dict{String, Float32}(
"H"  => 1.00798,
"C"  => 12.0106,
"N"  => 14.0069,
"O"  => 15.9994,
"P"  => 30.9738,
"S"  => 32.0680,
"CL" => 35.4520,
"NA" => 22.9898
                                     )


text:
using ApproxFun, MultivariateOrthogonalPolynomials, Test, StaticArrays
import MultivariateOrthogonalPolynomials: DirichletTriangle
import ApproxFunBase: testbandedblockbandedoperator, testblockbandedoperator, Vec
x,y=0.1,0.2

@testset "Dirchlet conversions" begin
    C=Conversion(DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        u=C*Fun(domainspace(C),[1.0])
        @test u(0.1,0.2) ≈ 1
        u=C*Fun(domainspace(C),[0,1.0])
        @test u(0.1,0.2) ≈ 0.1*Fun(JacobiTriangle(1,0,0),[1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[0,0,1.0])
        @test u(0.1,0.2) ≈ Fun(JacobiTriangle(0,0,0),[0,0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[0,0,0,1.0])
        @test u(0.1,0.2) ≈ 0.1*Fun(JacobiTriangle(1,0,0),[0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[0,0,0,0,1.0])
        @test u(0.1,0.2) ≈ 0.1*Fun(JacobiTriangle(1,0,0),[0,0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[0,0,0,0,0,1.0])
        @test u(0.1,0.2) ≈ Fun(JacobiTriangle(0,0,0),[0,0,0,0,0,1.])(0.1,0.2)

    C=Conversion(DirichletTriangle{0,1,0}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        u=C*Fun(domainspace(C),[1.0])
        @test u(0.1,0.2) ≈ 1
        u=C*Fun(domainspace(C),[0,1.0])
        @test u(0.1,0.2) ≈ Fun(Legendre(0..1),[0,1.])(0.1)
        u=C*Fun(domainspace(C),[zeros(2);1.0])
        @test u(0.1,0.2) ≈ 0.2Fun(JacobiTriangle(0,1,0),[1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[zeros(3);1.0])
        @test u(0.1,0.2) ≈ Fun(Legendre(0..1),[0,0,1.])(0.1)
        u=C*Fun(domainspace(C),[zeros(4);1.0])
        @test u(0.1,0.2) ≈ 0.2*Fun(JacobiTriangle(0,1,0),[0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[zeros(5);1.0])
        @test u(0.1,0.2) ≈ 0.2*Fun(JacobiTriangle(0,1,0),[0,0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[zeros(6);1.0])
        @test u(0.1,0.2) ≈ Fun(Legendre(0..1),[0,0,0,1.])(0.1)

    C=Conversion(DirichletTriangle{0,0,1}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        u=C*Fun(domainspace(C),[1.0])
        @test u(0.1,0.2) ≈ 1
        u=C*Fun(domainspace(C),[0,1.0])
        @test u(0.1,0.2) ≈ Fun(Legendre(0..1),[0,1.])(0.1)
        u=C*Fun(domainspace(C),[zeros(2);1.0])
        @test u(0.1,0.2) ≈ (1-x-y)*Fun(JacobiTriangle(0,0,1),[1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[zeros(3);1.0])
        @test u(0.1,0.2) ≈ Fun(Legendre(0..1),[0,0,1.])(0.1)
        u=C*Fun(domainspace(C),[zeros(4);1.0])
        @test u(0.1,0.2) ≈ (1-x-y)*Fun(JacobiTriangle(0,0,1),[0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[zeros(5);1.0])
        @test u(0.1,0.2) ≈ (1-x-y)*Fun(JacobiTriangle(0,0,1),[0,0,1.])(0.1,0.2)
        u=C*Fun(domainspace(C),[zeros(6);1.0])
        @test u(0.1,0.2) ≈ Fun(Legendre(0..1),[0,0,0,1.])(0.1)

    C=Conversion(DirichletTriangle{1,1,0}(),DirichletTriangle{0,1,0}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        n,k=0,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ 1
        n,k=1,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*Fun(Jacobi(1,0,0..1),[zeros(n-1);1.0])(x)
        n,k=1,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*Fun(JacobiTriangle(0,1,0),[zeros(sum(0:(n-1))+n-1);1.0])(x,y)
        n,k=2,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*Fun(Jacobi(1,0,0..1),[zeros(n-1);1.0])(x)
        n,k=2,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*y*Fun(JacobiTriangle(1,1,0),[zeros(sum(0:(n-2))+k-1);1.0])(x,y)
        n,k=2,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*Fun(JacobiTriangle(0,1,0),[zeros(sum(0:(n-1))+n-1);1.0])(x,y)
        n,k=3,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*Fun(Jacobi(1,0,0..1),[zeros(n-1);1.0])(x)
        n,k=3,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*y*Fun(JacobiTriangle(1,1,0),[zeros(sum(0:(n-2))+k-1);1.0])(x,y)
        n,k=3,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*y*Fun(JacobiTriangle(1,1,0),[zeros(sum(0:(n-2))+k-1);1.0])(x,y)
        n,k=3,3;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*Fun(JacobiTriangle(0,1,0),[zeros(sum(0:(n-1))+n-1);1.0])(x,y)

    C̃=Conversion(DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C̃)
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]

    C=Conversion(DirichletTriangle{1,0,1}(),DirichletTriangle{0,0,1}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        n,k=0,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ 1
        n,k=1,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*Fun(Jacobi(1,0,0..1),[zeros(n-1);1.0])(x)
        n,k=1,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x-y)*Fun(JacobiTriangle(0,0,1),[zeros(sum(0:(n-1))+n-1);1.0])(x,y)
        n,k=2,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*Fun(Jacobi(1,0,0..1),[zeros(n-1);1.0])(x)
        n,k=2,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x-y)*Fun(JacobiTriangle(1,0,1),[zeros(sum(0:(n-2))+k-1);1.0])(x,y)
        n,k=2,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x-y)*Fun(JacobiTriangle(0,0,1),[zeros(sum(0:(n-1))+n-1);1.0])(x,y)
        n,k=3,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*Fun(Jacobi(1,0,0..1),[zeros(n-1);1.0])(x)
        n,k=3,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x-y)*Fun(JacobiTriangle(1,0,1),[zeros(sum(0:(n-2))+k-1);1.0])(x,y)
        n,k=3,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x-y)*Fun(JacobiTriangle(1,0,1),[zeros(sum(0:(n-2))+k-1);1.0])(x,y)
        n,k=3,3;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x-y)*Fun(JacobiTriangle(0,0,1),[zeros(sum(0:(n-1))+n-1);1.0])(x,y)

    C̃=Conversion(DirichletTriangle{1,0,1}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C̃)
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]

    C=Conversion(DirichletTriangle{0,1,1}(),DirichletTriangle{0,1,0}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        n,k=0,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ 1
        n,k=1,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x)*Fun(Jacobi(0,1,0..1),[zeros(n-1);1.0])(x)
        n,k=1,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x-2y)*Fun(Jacobi(0,1,0..1),[zeros(n-1);1.0])(x)
        n,k=2,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x)*Fun(Jacobi(0,1,0..1),[zeros(n-1);1.0])(x)
        n,k=2,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x-2y)*Fun(Jacobi(0,1,0..1),[zeros(n-1);1.0])(x)
        n,k=2,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*(1-x-y)*Fun(JacobiTriangle(0,1,1),[zeros(sum(0:(n-2))+k-2);1.0])(x,y)
        n,k=3,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x)*Fun(Jacobi(0,1,0..1),[zeros(n-1);1.0])(x)
        n,k=3,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ (1-x-2y)*Fun(Jacobi(0,1,0..1),[zeros(n-1);1.0])(x)
        n,k=3,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*(1-x-y)*Fun(JacobiTriangle(0,1,1),[zeros(sum(0:(n-2))+k-2);1.0])(x,y)
        n,k=3,3;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*(1-x-y)*Fun(JacobiTriangle(0,1,1),[zeros(sum(0:(n-2))+k-2);1.0])(x,y)

    C̃=Conversion(DirichletTriangle{0,1,1}(),DirichletTriangle{0,0,1}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C̃)
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]

    C=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{0,1,1}(),DirichletTriangle{0,0,1}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C)
        n,k=0,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ 1
        n,k=1,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ 1-2x
        n,k=1,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ 1-x-2y
        n,k=2,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x)*Fun(Jacobi(1,1,0..1),[zeros(n-2);1.0])(x)
        n,k=2,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x-2y)*Fun(Jacobi(1,1,0..1),[zeros(n-2);1.0])(x)
        n,k=2,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*(1-x-y)*Fun(JacobiTriangle(0,1,1),[zeros(sum(0:(n-2))+n-2);1.0])(x,y)
        n,k=3,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x)*Fun(Jacobi(1,1,0..1),[zeros(n-2);1.0])(x)
        n,k=3,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x-2y)*Fun(Jacobi(1,1,0..1),[zeros(n-2);1.0])(x)
        n,k=3,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*y*(1-x-y)*Fun(JacobiTriangle(1,1,1),[zeros(sum(0:(n-3))+k-2);1.0])(x,y)
        n,k=3,3;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*(1-x-y)*Fun(JacobiTriangle(0,1,1),[zeros(sum(0:(n-2))+n-2);1.0])(x,y)
        n,k=4,0;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x)*Fun(Jacobi(1,1,0..1),[zeros(n-2);1.0])(x)
        n,k=4,1;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*(1-x-2y)*Fun(Jacobi(1,1,0..1),[zeros(n-2);1.0])(x)
        n,k=4,2;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*y*(1-x-y)*Fun(JacobiTriangle(1,1,1),[zeros(sum(0:(n-3))+k-2);1.0])(x,y)
        n,k=4,3;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ x*y*(1-x-y)*Fun(JacobiTriangle(1,1,1),[zeros(sum(0:(n-3))+k-2);1.0])(x,y)
        n,k=4,4;u=C*Fun(domainspace(C),[zeros(sum(0:n)+k);1.0])
        @test u(x,y) ≈ y*(1-x-y)*Fun(JacobiTriangle(0,1,1),[zeros(sum(0:(n-2))+n-2);1.0])(x,y)

    C̃=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,0,1}(),DirichletTriangle{0,0,1}(),JacobiTriangle(0,0,0))
        testbandedblockbandedoperator(C̃)
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]

    C̃=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,0,1}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]

    C̃=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{0,1,0}(),JacobiTriangle(0,0,0))
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]

    C̃=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
        @test C[1:20,1:20] ≈ C̃[1:20,1:20]
end

@testset "Dirichlet evaluation" begin
    f=Fun((x,y)->exp(x-0.12*cos(y)),JacobiTriangle(0,0,0))

    C=Conversion(DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
    R=Conversion(domainspace(C),Legendre(Segment(Vec(0.,0.),Vec(0.,1.))))
    u=C\f
    @test (R*u)(0.,0.3) ≈ f(0.,0.3)
    @test u(0.1,0.2) ≈ f(0.1,0.2)

    C=Conversion(DirichletTriangle{0,1,0}(),JacobiTriangle(0,0,0))
    R=Conversion(domainspace(C),Legendre(Segment(Vec(0.,0.),Vec(1.,0.))))
    u=C\f
    @test (R*u)(0.3,0.) ≈ f(0.3,0.)
    @test u(0.1,0.2) ≈ f(0.1,0.2)

    C=Conversion(DirichletTriangle{0,0,1}(),JacobiTriangle(0,0,0))
    R=Conversion(domainspace(C),Legendre(Segment(Vec(0.,1.),Vec(1.,0.))))
    u=C\f
    @test (R*u)(0.3,1-0.3) ≈ f(0.3,1-0.3)
    @test u(0.1,0.2) ≈ f(0.1,0.2)

    C=Conversion(DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
    Rx=Conversion(DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),Legendre(Segment(Vec(0.,0.),Vec(0.,1.))))
    Ry=Conversion(DirichletTriangle{1,1,0}(),DirichletTriangle{0,1,0}(),Legendre(Segment(Vec(0.,0.),Vec(1.,0.))))
    u=C\f
    @test (Rx*u)(0.,0.3) ≈ f(0.,0.3)
    @test (Ry*u)(0.3,0.) ≈ f(0.3,0.)
    @test u(0.1,0.2) ≈ f(0.1,0.2)

    C=Conversion(DirichletTriangle{1,0,1}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
    Rx=Conversion(DirichletTriangle{1,0,1}(),DirichletTriangle{1,0,0}(),Legendre(Segment(Vec(0.,0.),Vec(0.,1.))))
    Rz=Conversion(DirichletTriangle{1,0,1}(),DirichletTriangle{0,0,1}(),Legendre(Segment(Vec(0.,1.),Vec(1.,0.))))
    u=C\f
    @test (Rx*u)(0.,0.3) ≈ f(0.,0.3)
    @test (Rz*u)(0.3,1-0.3) ≈ f(0.3,1-0.3)
    @test u(0.1,0.2) ≈ f(0.1,0.2)

    C=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
    Rx=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),Legendre(Segment(Vec(0.,0.),Vec(0.,1.))))
    Ry=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{0,1,0}(),Legendre(Segment(Vec(0.,0.),Vec(1.,0.))))
    Rz=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{0,1,1}(),DirichletTriangle{0,0,1}(),Legendre(Segment(Vec(0.,1.),Vec(1.,0.))))
    u=C\f
    @test (Rx*u)(0.,0.3) ≈ f(0.,0.3)
    @test (Ry*u)(0.3,0.) ≈ f(0.3,0.)
    @test (Rz*u)(0.3,1-0.3) ≈ f(0.3,1-0.3)
    @test u(0.1,0.2) ≈ f(0.1,0.2)
end

@testset "Dirichlet restriction" begin
    f=Fun((x,y)->exp(x-0.12*cos(y)),JacobiTriangle(0,0,0))
    C=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
    u=C\f

    B=Dirichlet(DirichletTriangle{1,1,1}())
    testblockbandedoperator(B)
    ∂u=B*u
    @test ∂u(0.3,0.)≈ f(0.3,0.)
    @test ∂u(0.,0.3)≈ f(0.,0.3)
    @test ∂u(0.3,1-0.3)≈ f(0.3,1-0.3)

    d = Triangle(Vec(2,3),Vec(3,4),Vec(1,6))
    f=Fun((x,y)->exp(x-0.12*cos(y)),JacobiTriangle(0,0,0,d))
    C=Conversion(DirichletTriangle{1,1,1}(d),DirichletTriangle{1,1,0}(d),DirichletTriangle{1,0,0}(d),JacobiTriangle(0,0,0,d))
    u=C\f

    B=Dirichlet(DirichletTriangle{1,1,1}(d))
    testblockbandedoperator(B)
    ∂u=B*u
    @test ∂u(2.5,3.5)≈ f(2.5,3.5)
    @test ∂u(1.5,4.5)≈ f(1.5,4.5)
    @test ∂u(2,5)≈ f(2,5)

end

@testset "Triangle Dirichlet Derivatives" begin
    @testset "Triangle() Derivative" begin
        S = DirichletTriangle{1,0,1}()
        Dx = Derivative(S, [1,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle()
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,0])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        Dx = Derivative(S, [2,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(1,0,1)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([2,0])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        Dx = Derivative(S, [2,1])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(1,1,2)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([2,1])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        S = DirichletTriangle{0,1,1}()
        Dy = Derivative(S, [0,1])
        testbandedblockbandedoperator(Dy)
        @test rangespace(Dy) == JacobiTriangle()
        f = Fun(S, randn(20))
        @test (Dy*f)(0.1,0.2) ≈ (Derivative([0,1])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        Dx = Derivative(S, [1,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(1,0,1)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,0])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        Dx = Derivative(S, [2,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(2,0,2)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([2,0])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        Dx = Derivative(S, [1,1])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(1,0,1)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,1])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)

        S = DirichletTriangle{1,1,1}()
        Dx = Derivative(S,[1,0])
        @test rangespace(Dx) == JacobiTriangle()
        testbandedblockbandedoperator(Dx)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,0])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)
        Dy = Derivative(S, [0,1])
        testbandedblockbandedoperator(Dy)
        @test rangespace(Dy) == JacobiTriangle()
        f = Fun(S, randn(20))
        @test (Dy*f)(0.1,0.2) ≈ (Derivative([0,1])*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)
    end

    @testset "Other triangle Derivative" begin
        d = Triangle(Vec(2,3),Vec(3,4),Vec(1,6))
        S = DirichletTriangle{1,0,1}(d)
        Dx = Derivative(S, [1,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(1,1,1,d)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,0])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)

        Dx = Derivative(S, [2,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(2,2,2,d)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([2,0])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)

        Dx = Derivative(S, [2,1])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(3,3,3,d)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([2,1])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)

        S = DirichletTriangle{0,1,1}(d)
        Dy = Derivative(S, [0,1])
        testbandedblockbandedoperator(Dy)
        @test rangespace(Dy) == JacobiTriangle(1,1,1,d)
        f = Fun(S, randn(20))
        @test (Dy*f)(0.1,0.2) ≈ (Derivative([0,1])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)

        Dx = Derivative(S, [1,0])
        testbandedblockbandedoperator(Dx)
        @test rangespace(Dx) == JacobiTriangle(1,1,1,d)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,0])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)

        S = DirichletTriangle{1,1,1}(d)
        Dx = Derivative(S,[1,0])
        @test rangespace(Dx) == JacobiTriangle(0,0,0,d)
        testbandedblockbandedoperator(Dx)
        f = Fun(S, randn(20))
        @test (Dx*f)(0.1,0.2) ≈ (Derivative([1,0])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)
        Dy = Derivative(S, [0,1])
        testbandedblockbandedoperator(Dy)
        @test rangespace(Dy) == JacobiTriangle(0,0,0,d)
        f = Fun(S, randn(20))
        @test (Dy*f)(0.1,0.2) ≈ (Derivative([0,1])*Fun(f,JacobiTriangle(0,0,0,d)))(0.1,0.2)
    end
end

@testset "Triangle Dirichlet Laplacian" begin
    @testset "Triangle Laplace" begin
        S = DirichletTriangle{1,1,1}()
        Δ=Laplacian(S)
        testbandedblockbandedoperator(Δ)
        @test rangespace(Δ) == JacobiTriangle(1,1,1)
        f = Fun(S, randn(20))
        @test (Δ*f)(0.1,0.2) ≈ (Laplacian()*Fun(f,JacobiTriangle(0,0,0)))(0.1,0.2)
        B=Dirichlet(S)
        f=Fun((x,y)->real(exp(x+im*y)),rangespace(B))
        u=\([B;Δ],[f;0.];tolerance=1E-13)
        @test u(0.1,0.2) ≈ real(exp(0.1+0.2im))
    end

    @testset "Triangle Laplace via JacobiTriangle(0,0,0)" begin
        C=Conversion(DirichletTriangle{1,1,1}(),DirichletTriangle{1,1,0}(),DirichletTriangle{1,0,0}(),JacobiTriangle(0,0,0))
        B=Dirichlet(DirichletTriangle{1,1,1}())
        Δ=Laplacian(rangespace(C))
        f=Fun((x,y)->real(exp(x+im*y)),rangespace(B))
        u=\([B;Δ*C],[f;0.];tolerance=1E-13)

        @test u(0.1,0.2) ≈ real(exp(0.1+0.2im))
    end

    @testset "Other triangle Laplace" begin
        d = Triangle(Vec(2,3),Vec(3,4),Vec(1,6))
        S = DirichletTriangle{1,1,1}(d)
        Δ=Laplacian(S)
        @test rangespace(Δ) == JacobiTriangle(1,1,1,d)
        testbandedblockbandedoperator(Δ)
        f = Fun(S, randn(5))
        x,y = fromcanonical(S, 0.1,0.2)
        @test (Δ*f)(x,y) ≈ (Laplacian()*Fun(f, JacobiTriangle(0,0,0,d)))(0.1,0.2)

        B = Dirichlet(d)
        f=Fun((x,y)->real(exp(x+im*y)),rangespace(B))
        u=\([B;Δ],[f;0.];tolerance=1E-13)

        @test u(x,y) ≈ real(exp(x+y*im))
    end
end


text:
export CalorEspecifico
typealias DanaCalorEspecifico DanaRealParametric
type _CalorEspecifico
	function _CalorEspecifico(_::Dict{Symbol,Any})
		fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
		fields[:Default]=1e-3
		fields[:Lower]=0
		fields[:Upper]=1
		fields[:Unit]="MJ/kg/K"
		drive!(fields,_)
		new(fields)
	end
	value::Dict{Symbol,Any}
end
typealias CalorEspecifico DanaRealParametric{_CalorEspecifico}


text:
function HDF5._getindex(dset::HDF5Dataset, T::Type, indices::Union{AbstractRange{Int},Int}...)
    dsel_id = hyperslab(dset, indices...)
    ret = Array{T}(undef,map(length, indices))
    memtype = datatype(ret)
    memspace = dataspace(ret)
    try
        h5d_read(dset.id, memtype.id, memspace.id, dsel_id, dset.xfer, ret)
    finally
        close(memtype)
        close(memspace)
        h5s_close(dsel_id)
    end
    ret
end

function HDF5._setindex!(dset::HDF5Dataset,T::Type, X::Array, indices::Union{AbstractRange{Int},Int}...)
    if !(T <: Array)
        error("Dataset indexing (hyperslab) is available only for arrays")
    end
    ET = eltype(T)
    if length(X) != prod(map(length, indices))
        error("number of elements in range and length of array must be equal")
    end
    if eltype(X) != ET
        X = convert(Array{ET}, X)
    end
    dsel_id = hyperslab(dset, indices...)
    memtype = datatype(X)
    memspace = dataspace(X)
    try
        h5d_write(dset.id, memtype.id, memspace.id, dsel_id, dset.xfer, X)
    finally
        close(memtype)
        close(memspace)
        h5s_close(dsel_id)
    end
    X
end


text:
# The structure follows https://developer.spotify.com/documentation/web-api/reference/

export Albums, Artists, Browse, Episodes, Follow, Library
export Markets, Personalization, Player, Playlists, Search
export Shows, Tracks, UsersProfile, Objects
@info "Exported Library & etc."
"import export shorthand"
macro _ie(foo...)
    expr = :()
    for f in foo
        push!(expr.args, esc(
            quote
               import ..($f)
               export $f
            end))
    end
    expr
end
macro _ie()
    nothing
end

baremodule Albums
    import ..@_ie
    @_ie album_get album_get_tracks
end
baremodule Artists
    import ..@_ie
    @_ie artist_get artist_get_albums artist_top_tracks artist_get_related_artists
end
baremodule Browse
    import ..@_ie
    @_ie category_get    category_get_playlist    category_get_several
    @_ie category_get_new_releases    category_get_featured_playlist
    @_ie recommendations_get  recommendations_dict_parser
    # single parameter validation - integer or float
    @_ie min_acousticness  min_danceability  min_duration_ms  min_energy  min_instrumentalness
    @_ie min_key  min_liveness  min_loudness  min_mode  min_popularity  min_speechiness
    @_ie min_tempo  min_time_signature  min_valence  max_acousticness  max_danceability
    @_ie max_duration_ms  max_energy  max_instrumentalness  max_key  max_liveness
    @_ie max_loudness  max_mode  max_popularity  max_speechiness  max_tempo
    @_ie max_time_signature  max_valence  target_acousticness  target_danceability
    @_ie target_duration_ms  target_energy  target_instrumentalness  target_key
    @_ie target_liveness  target_loudness  target_mode  target_popularity
    @_ie target_speechiness  target_tempo  target_time_signature  target_valence
end
baremodule Episodes
    import ..@_ie
    @_ie episodes_get
end
baremodule Follow
    import ..@_ie
    @_ie follow_check follow_check_playlist follow_artists follow_artists_users follow_playlist unfollow_artists_users unfollow_playlist
end
baremodule Library
    import ..@_ie
    @_ie library_get_saved_tracks library_get_saved_shows library_get_saved_albums library_check_saved_tracks
    @_ie library_check_saved_shows library_check_saved_albums library_remove_albums library_remove_shows
    @_ie library_remove_tracks library_save_album library_save_show library_save_track
end
baremodule Markets
    import ..@_ie
    @_ie
end

baremodule Personalization
    import ..@_ie
    @_ie top_tracks top_artists
end
baremodule Player
    import ..@_ie
    @_ie show_get show_get_episodes
end
baremodule Playlists
    import ..@_ie
    @_ie tracks_get_audio_analysis tracks_get_audio_features tracks_get
end
baremodule Search
    import ..@_ie
    @_ie
end
baremodule Shows
    import ..@_ie
    @_ie
end
baremodule Tracks
    import ..@_ie
    @_ie tracks_get
end
baremodule UsersProfile
    import ..@_ie
    @_ie
end
baremodule Objects
    import ..@_ie
    @_ie
end


text:
module MuJoCo

using MacroTools: combinestructdef
using Base: RefValue
using AxisArrays: Axis, AxisArray
using
    UnsafeArrays,
    StaticArrays,
    Reexport,
    BangBang

const TESTMODELXML = joinpath(@__DIR__, "../test/testmodel.xml")

export
    # Sugar
    jlData,
    jlModel,
    namify,
    jl_enabled,
    jl_disabled,
    jl_name2id,
    jl_id2name,

    # Core
    MJCore,
    mj_activate,
    mj_resetData,
    mj_resetCallbacks,
    mj_step,
    mj_step1,
    mj_step2,
    mj_forward,
    mj_inverse,
    mj_forwardSkip,
    mj_inverseSkip

include("MJCore/MJCore.jl")
using .MJCore
using .MJCore: AbstractMJModel,
                AbstractMJData,
                UnsafeVector,
                UnsafeMatrix,
                MuJoCoException,
                @mjerror,
                MuJoCoException,
                union_types

####
#### Sugar
####

function build_structdef(
    name::Symbol,
    fields::Vector;
    supertype = Any,
    constructors = Any[],
    mutable = false,
    params = Any[],
)
    Dict{Symbol,Any}(
        :name => name,
        :fields => fields,
        :supertype => supertype,
        :constructors => constructors,
        :mutable => mutable,
        :params => params,
    )
end

function build_funcdef(
    name::Symbol,
    body::Expr;
    args = Any[],
    kwargs = Any[],
    whereparams = (),
)
    Dict{Symbol,Any}(
        :name => name,
        :args => args,
        :kwargs => kwargs,
        :body => body,
        :whereparams => whereparams,
    )
end

include("metadata.jl")
include("jlmodel.jl")
include("jldata.jl")
include("jlmujoco.jl")
include("namify.jl")

const JLTypes = Union{jlModel,jlData}

const STATE_FIELDS = (
    :time,
    :qpos,
    :qvel,
    :act,
    :qacc_warmstart,
    :mocap_pos,
    :mocap_quat,
    :userdata,
)

const CONTROL_FIELDS = (
    :ctrl,
    :xfrc_applied,
    :qfrc_applied
)

# For @set!! macro to work
@inline function BangBang.setproperties!(value::JLTypes, patch::NamedTuple)
    ntuple(length(patch)) do i
        setproperty!(value, keys(patch)[i], patch[i])
    end
    return value
end
BangBang.possible(::typeof(BangBang.setproperties!), obj::JLTypes, patch) = true

function printstate(io, jd::jlData)
    header = "jlData"
    n = 10
    println(io, repeat('-', n), header, repeat('-', n))
    println(io, "MuJoCo State:")
    for name in STATE_FIELDS
        x = getproperty(jd, name)
        println(io, "$(name):")
        show(io, x)
        println(io)
    end
    println(io, "MuJoCo Controls:")
    for name in CONTROL_FIELDS
        x = getproperty(jd, name)
        println(io, "$(name):")
        show(io, x)
        println(io)
    end
    println(io, repeat('-', n * 2 + length(header)))
    io
end

function Base.show(io::IO, ::MIME"text/plain", d::jlData)
    printstate(io, d)
    mjshow(io, d, 0, filter= (n, T) -> T <: Number)
    io
end

function Base.show(io::IO, ::MIME"text/plain", m::jlModel)
    mjshow(io, m, 0, filter=(n, T) -> T <: Number)
    io
end

Base.show(io::IO, ::MIME"text/plain", x::Union{MJCore.MJTypes,JLTypes}) = mjshow(io, x, 0)

function mjshow(io::IO, x::Union{MJCore.MJTypes,JLTypes}, level::Int; filter=(x, T)->true)
    alltypes = Union{MJCore.MJTypes,JLTypes}
    for name in propertynames(x)
        field = getproperty(x, name)
        T = typeof(field)
        !filter(name, T) && continue

        if (T <: AbstractArray{<:alltypes} ||
            T <: Ptr || name === :ref || (x isa jlData && name in STATE_FIELDS))
            continue
        elseif T <: alltypes
            mjshow(io, field, level + 1)
        elseif !(T <: Ptr)
            for _ = 1:level
                print(io, " ")
            end
            print(io, "$name: ")
            show(io, field)
            println(io)
        end
    end
end

end # module


text:
# Geometry-related routines, including distances and transformations

module Geometry

using LinearAlgebra
using Statistics
using ..Dorothy.Utils
using StaticArrays

export
		Vector3D, Basis3D,

		sqnorm, dist, sqdist, dihedral, com, cog,

		Transformation, transform, Translation, Scaling, LinearTransformation,
		AffineTransformation, translation, scaling, rotation, superposition,
		SuperpositionBuffer, FitBuffer, fitline, fitplane, rmsd, msd,

		BoundingBox, extent, dims, center, volume, isinside

struct Vector3D <: FieldVector{3,Float64}
	x::Float64
	y::Float64
	z::Float64
end

StaticArrays.similar_type(::Type{Vector3D}, ::Type{Float64}, ::Size{(3,)}) =
		Vector3D

const Basis3D = SMatrix{3,3,Float64,9}

sqnorm(V::AbstractVector{<:Real}) = sqnorm(Vector3D(V))

function sqnorm(V::Vector3D)
	d2 = V.x^2
	d2 += V.y^2
	d2 += V.z^2
end

LinearAlgebra.norm(V::Vector3D) = sqrt(sqnorm(V))

dist(R1::AbstractVector{<:Real}, R2::AbstractVector{<:Real}) =
		dist(Vector3D(R1), Vector3D(R2))

dist(R1::Vector3D, R2::Vector3D) = sqrt(sqdist(R1, R2))

sqdist(R1::AbstractVector{<:Real}, R2::AbstractVector{<:Real}) =
		sqdist(Vector3D(R1), Vector3D(R2))

function sqdist(R1::Vector3D, R2::Vector3D)
	d2 = (R2.x - R1.x)^2
	d2 += (R2.y - R1.y)^2
	d2 += (R2.z - R1.z)^2
end

Base.angle(V1::AbstractVector{<:Real}, V2::AbstractVector{<:Real}) =
		angle(Vector3D(V1), Vector3D(V2))

Base.angle(V1::Vector3D, V2::Vector3D) = acos(V1⋅V2 / (norm(V1)*norm(V2)))

Base.angle(R1::AbstractVector{<:Real}, R2::AbstractVector{<:Real},
		R3::AbstractVector{<:Real}) =
		angle(Vector3D(R1), Vector3D(R2), Vector3D(R3))

Base.angle(R1::Vector3D, R2::Vector3D, R3::Vector3D) = angle(R1-R2, R3-R2)

dihedral(V1::AbstractVector{<:Real}, V2::AbstractVector{<:Real},
		V3::AbstractVector{<:Real}) =
		dihedral(Vector3D(V1), Vector3D(V2), Vector3D(V3))

function dihedral(V1::Vector3D, V2::Vector3D, V3::Vector3D)
	V1xV2 = V1 × V2
	V2xV3 = V2 × V3
	y = (V1xV2 × V2xV3)⋅(V2/norm(V2))
	x = V1xV2 ⋅ V2xV3
	atan(y,x)
end

dihedral(R1::AbstractVector{<:Real}, R2::AbstractVector{<:Real},
		R3::AbstractVector{<:Real}, R4::AbstractVector{<:Real}) =
		dihedral(Vector3D(R1), Vector3D(R2), Vector3D(R3), Vector3D(R4))

dihedral(R1::Vector3D, R2::Vector3D, R3::Vector3D, R4::Vector3D) =
		dihedral(R2-R1, R3-R2, R4-R3)

function com(R::AbstractVector{Vector3D}, W::AbstractVector{<:Real})
	@boundscheck begin
		length(R) > 0 || error("com of zero positions is undefined")
		length(R) == length(W) ||
				error("size mismatch between position/weight arrays")
	end
	sumx, sumy, sumz, sumW = 0.0, 0.0, 0.0, 0.0
	for i in eachindex(W)
		Wi = W[i]
		Ri = R[i]
		sumRx += Ri.x * Wi
		sumRy += Ri.y * Wi
		sumRz += Ri.z * Wi
		sumW += Wi
	end
	Vector3D(sumRx / sumW, sumRy / sumW, sumRz / sumW)
end

function cog(R::AbstractVector{Vector3D})
	@boundscheck length(R) > 0 || error("cog of zero positions is undefined")
	sumRx, sumRy, sumRz = 0.0, 0.0, 0.0
	for Ri in R
		sumRx += Ri.x
		sumRy += Ri.y
		sumRz += Ri.z
	end
	n = length(R)
	Vector3D(sumRx / n, sumRy / n, sumRz / n)
end

abstract type Transformation end

transform(T::Transformation, R::AbstractVector{<:Real}) =
		transform(T, Vector3D(R))

transform(T::Transformation, R::Vector3D) = T * R

transform(T::Transformation, R::AbstractVector{Vector3D}) = (R .= T .* R)

Base.broadcastable(T::Transformation) = Ref(T)

struct Translation <: Transformation
	V::Vector3D
end

Base.inv(T::Translation) = Translation(-T.V)

Base.:(*)(T::Translation, R::AbstractVector{<:Real}) = T * Vector3D(R)

Base.:(*)(T::Translation, R::Vector3D) = T.V + R

Base.:(*)(T2::Translation, T1::Translation) = Translation(T2.V + T1.V)

(T::Translation)(x) = transform(T, x)

struct Scaling <: Transformation
	V::Vector3D
end

Base.inv(T::Scaling) = Scaling(1.0 ./ T.V)

Base.:(*)(T::Scaling, R::AbstractVector{<:Real}) = T * Vector3D(R)

Base.:(*)(T::Scaling, R::Vector3D) = T.V .* R

Base.:(*)(T2::Scaling, T1::Scaling) = Scaling(T2.V .* T1.V)

(T::Scaling)(x) = transform(T, x)

struct LinearTransformation <:Transformation
	M::Basis3D
end

Base.inv(T::LinearTransformation) = LinearTransformation(inv(T.M))

Base.:(*)(T::LinearTransformation, R::AbstractVector{<:Real}) = T * Vector3D(R)

Base.:(*)(T::LinearTransformation, R::Vector3D) = T.M * R

Base.:(*)(T2::LinearTransformation, T1::LinearTransformation) =
		LinearTransformation(T2.M * T1.M)

Base.:(*)(T2::LinearTransformation, T1::Scaling) =
		LinearTransformation(T2.M * Diagonal(T1.V))

Base.:(*)(T2::Scaling, T1::LinearTransformation) =
		LinearTransformation(Diagonal(T2) * T1.M)

(T::LinearTransformation)(x) = transform(T, x)

struct AffineTransformation <: Transformation
	T1::LinearTransformation
	T2::Translation
end

Base.inv(T::AffineTransformation) =
		convert(AffineTransformation, inv(convert(Matrix{Float64}, T)))

Base.:(*)(T::AffineTransformation, R::AbstractVector{<:Real}) = T * Vector3D(R)

Base.:(*)(T::AffineTransformation, R::Vector3D) = T.T2 * (T.T1 * R)

Base.:(*)(T2::Translation, T1::LinearTransformation) =
		AffineTransformation(T1, T2)

Base.:(*)(T2::Translation, T1::AffineTransformation) =
		AffineTransformation(T1.T1, T2 * T1.T2)

function Base.:(*)(T2::AffineTransformation, T1::AffineTransformation)
	M2 = convert(Matrix{Float64}, T2)
	M1 = convert(Matrix{Float64}, T1)
	convert(AffineTransformation, M2 * M1)
end

(T::AffineTransformation)(x) = transform(T, x)

function Base.convert(::Type{Matrix{Float64}}, T::AffineTransformation)
	M = zeros(Float64, 4,4)
	for i = 1:3
		for j = 1:3
			M[j,i] = T.T1.M[j,i]
		end
	end
	for j = 1:3
		M[j,4] = T.T2.V[j]
	end
	M[4,4] = 1.0
	M
end

function Base.convert(::Type{AffineTransformation}, M::AbstractMatrix{<:Real})
	@boundscheck size(M) == (4,4) && M[4,:] == [0.0, 0.0, 0.0, 1.0] ||
			error("expected 4×4 affine transformation matrix")
	AffineTransformation(LinearTransformation(M[1:3,1:3]),
	Translation(M[1:3,4]))
end

Base.convert(::Type{AffineTransformation}, T::Translation) =
		AffineTransformation(LinearTransformation(Matrix{Float64}(I, 3,3)),
		Translation(copy(T.V)))

Base.convert(::Type{AffineTransformation}, T::LinearTransformation) =
		AffineTransformation(LinearTransformation(copy(T.M)),
		Translation(zeros(3)))

function Base.:(*)(T2::Transformation, T1::Transformation)
	T2p, T1p = promote(T2, T1)
	T2p * T1p
end

Base.promote_rule(::Type{<:Transformation}, ::Type{<:Transformation}) =
		AffineTransformation

translation(V::AbstractVector{<:Real}) = Translation(V)

translation(; x::Real = 0.0, y::Real = 0.0, z::Real = 0.0) =
		translation(Vector3D(x, y, z))

function scaling(V::AbstractVector{<:Real};
		O::Union{AbstractVector{<:Real},Nothing} = nothing)
	T = LinearTransformation(Diagonal(V))
	if O == nothing
		T
	else
		translation(O) * T * translation(-O)
	end
end

function scaling(c::Real; O::Union{AbstractVector{<:Real},Nothing} = nothing)
	T = scaling(Vector3D(c,c,c))
	if O == nothing
		T
	else
		translation(O) * T * translation(-O)
	end
end

function scaling(; x::Real = 1.0, y::Real = 1.0, z::Real = 1.0,
		O::Union{AbstractVector{<:Real},Nothing} = nothing)
	T = scaling(Vector3D(x,y,z))
	if O == nothing
		T
	else
		translation(O) * T * translation(-O)
	end
end

function rotation(θ::Real, V::AbstractVector{<:Real};
		O::Union{AbstractVector{<:Real},Nothing} = nothing)
	l, m, n = normalize(V)
	cosθ = cos(θ)
	sinθ = sin(θ)
	M = @SMatrix [l*l*(1-cosθ)+cosθ   m*l*(1-cosθ)-n*sinθ n*l*(1-cosθ)+m*sinθ;
		          l*m*(1-cosθ)+n*sinθ m*m*(1-cosθ)+cosθ   n*m*(1-cosθ)-l*sinθ;
		          l*n*(1-cosθ)-m*sinθ m*n*(1-cosθ)+l*sinθ n*n*(1-cosθ)+cosθ  ]
	T = LinearTransformation(M)
	if O == nothing
		T
	else
		translation(O) * T * translation(-O)
	end
end

function rotation(θ::Real; R1::AbstractVector{<:Real},
		R2::AbstractVector{<:Real},
		R3::Union{AbstractVector{<:Real},Nothing} = nothing)
	if R3 == nothing
		rotation(θ, R2-R1, O = R1)
	else
		V1 = R1 - R2
		V2 = R3 - R2
		rotation(θ, V1×V2, O = R2)
	end
end

function rotation(; x::Union{Real,Nothing} = nothing,
		y::Union{Real,Nothing} = nothing, z::Union{Real,Nothing} = nothing,
		O::Union{Real,Nothing} = nothing)
	if x != nothing
		if y != nothing || z != nothing
			error("expected only one of x, y or z as keyword argument")
		end
		if O == nothing
			rotation(x, Vector3D(1.0, 0.0, 0.0))
		else
			rotation(x, Vector3D(1.0, 0.0, 0.0); O = O)
		end
	elseif y != nothing
		if x != nothing || z != nothing
			error("expected only one of x, y or z as keyword argument")
		end
		if O == nothing
			rotation(y, Vector3D(0.0, 1.0, 0.0))
		else
			rotation(y, Vector3D(0.0, 1.0, 0.0); O = O)
		end
	elseif z != nothing
		if x != nothing || y != nothing
			error("expected only one of x, y or z as keyword argument")
		end
		if O == nothing
			rotation(z, Vector3D(0.0, 0.0, 1.0))
		else
			rotation(z, Vector3D(0.0, 0.0, 1.0); O = O)
		end
	else
		error("expected one of x, y or z as keyword argument")
	end
end

mutable struct SuperpositionBuffer
	R::Matrix{Float64}
	Rref::Matrix{Float64}
	H::Matrix{Float64}
	T::Matrix{Float64}
	center::Bool

	SuperpositionBuffer(; center::Bool = true) =
			new(Matrix{Float64}(undef, (0,0)), Matrix{Float64}(undef, (0,0)),
			Matrix{Float64}(undef, (3,3)), Matrix{Float64}(undef, (3,3)),
			center)
end

superposition(R::AbstractVector{Vector3D}, Rref::AbstractVector{Vector3D}) =
		superposition(R, Rref, SuperpositionBuffer())

superposition(R::AbstractVector{Vector3D}, Rref::AbstractVector{Vector3D},
		W::AbstractVector{<:Real}) =
		superposition(R, Rref, W, SuperpositionBuffer())

superposition(R::AbstractVector{Vector3D}, Rref::AbstractVector{Vector3D},
		buffer::SuperpositionBuffer) =
		superposition(R, Rref, Repeated(1.0, length(R)), buffer)

function superposition(R::AbstractVector{Vector3D},
		Rref::AbstractVector{Vector3D}, W::AbstractVector{<:Real},
		buffer::SuperpositionBuffer)
	n = length(R)
	@boundscheck begin
		n > 0 || error("fit of zero positions is undefined")
		length(Rref) == length(W) == n ||
				error("size mismatch between position/weight arrays")
	end
	if ncols(buffer.R) != n
		buffer.R = Matrix{Float64}(undef, (3,n))
		buffer.Rref = Matrix{Float64}(undef, (3,n))
	end
	for i = 1:n
		buffer.R[:,i] = R[i]
		buffer.Rref[:,i] = Rref[i]
	end
	if buffer.center
		cogR = cog(R)
		cogRref = cog(Rref)
		for i = 1:n
			buffer.R[:,i] -= cogR
			buffer.Rref[:,i] -= cogRref
		end
	end
	kabsch!(buffer.T, buffer.H, buffer.R, buffer.Rref, W)
	rotation = LinearTransformation(buffer.T)
	if buffer.center
		translation(cogRref) * rotation * translation(-cogR)
	else
		rotation
	end
end

function kabsch!(T::AbstractMatrix{<:Real}, H::AbstractMatrix{<:Real},
		R::AbstractMatrix{<:Real}, Rref::AbstractMatrix{<:Real},
		W::AbstractVector{<:Real})
	mul!(H, Rref, lmul!(Diagonal(W), R'))
	U, _, V = svd!(H)
	if sign(det(U) * det(V')) < 0
		V[:,3] .= .- V[:,3]
	end
	mul!(T, U, V')
end

mutable struct FitBuffer
	R::Matrix{Float64}
	Radj::Matrix{Float64}
	center::Bool

	FitBuffer(; center::Bool = true) = new(Matrix{Float64}(undef, (0,0)),
			Matrix{Float64}(undef, (0,0)),center)
end

fitline(R::AbstractVector{Vector3D}) = fitline(R, FitBuffer())

fitline(R::AbstractVector{Vector3D}, W::AbstractVector{<:Real}) =
		fitline(R, W, FitBuffer())

fitline(R::AbstractVector{Vector3D}, buffer::FitBuffer) =
		fitline(R, Repeated(1.0, length(R)), buffer)

function fitline(R::AbstractVector{Vector3D}, W::AbstractVector{<:Real},
		buffer::FitBuffer)
	n = length(R)
	@boundscheck begin
		n >= 2 || error("line fit requires two or more positions")
		length(W) == n || error("size mismatch between position/weight arrays")
	end
	if ncols(buffer.R) != n
		buffer.R = Matrix{Float64}(undef, (3,n))
		buffer.Radj = Matrix{Float64}(undef, (n,3))
	end
	for i = 1:n
		buffer.R[:,i] = R[i]
	end
	if buffer.center
		cogR = cog(R)
		for i = 1:n
			buffer.R[:,i] -= cogR
		end
	end
	T = svdlinefit!(buffer.Radj, buffer.R, W)
	Vector3D(T)
end

function svdlinefit!(Radj::AbstractMatrix{<:Real}, R::AbstractMatrix{<:Real},
		W::AbstractVector{<:Real})
	mul!(Radj, Diagonal(W), R')
	T = svd!(Radj).V[:,1]
	Rn = Radj'[:,end]
	if sign(Rn ⋅ T) < 0
		T .= .- T
	end
	T
end

fitplane(R::AbstractVector{Vector3D}) = fitplane(R, FitBuffer())

fitplane(R::AbstractVector{Vector3D}, W::AbstractVector{<:Real}) =
		fitplane(R, W, FitBuffer())

fitplane(R::AbstractVector{Vector3D}, buffer::FitBuffer) =
		fitplane(R, Repeated(1.0, length(R)), buffer)

function fitplane(R::AbstractVector{Vector3D}, W::AbstractVector{<:Real},
		buffer::FitBuffer)
	n = length(R)
	@boundscheck begin
		n >= 3 || error("plane fit requires three or more positions")
		length(W) == n || error("size mismatch between position/weight arrays")
	end
	if ncols(buffer.R) != n
		buffer.R = Matrix{Float64}(undef, (3,n))
		buffer.Radj = Matrix{Float64}(undef, (n,3))
	end
	for i = 1:n
		buffer.R[:,i] = R[i]
	end
	if buffer.center
		cogR = cog(R)
		for i = 1:n
			buffer.R[:,i] -= cogR
		end
	end
	T = svdplanefit!(buffer.Radj, buffer.R, W)
	Vector3D(T)
end

function svdplanefit!(Radj::AbstractMatrix{<:Real}, R::AbstractMatrix{<:Real},
		W::AbstractVector{<:Real})
	mul!(Radj, Diagonal(W), R')
	T = svd!(Radj).V[:,3]
	Rn = Radj'[:,end]
	if sign(Rn ⋅ T) < 0
		T .= .- T
	end
	T
end

rmsd(R1::AbstractVector{Vector3D}, R2::AbstractVector{Vector3D}) =
		sqrt(msd(R1, R2))

rmsd(R1::AbstractVector{Vector3D}, R2::AbstractVector{Vector3D},
		W::AbstractVector{<:Real}) = sqrt(msd(R1, R2, W))

function msd(R1::AbstractVector{Vector3D}, R2::AbstractVector{Vector3D})
	@boundscheck begin
		length(R1) > 0 || error("deviation of zero positions is undefined")
		length(R1) == length(R2) ||
				error("size mismatch between position arrays")
	end
	sd = 0.0
	for i in eachindex(R1)
		R1i = R1[i]
		R2i = R2[i]
		sd += (R2i.x - R1i.x)^2
		sd += (R2i.y - R1i.y)^2
		sd += (R2i.z - R1i.z)^2
	end
	sd / length(R1)
end

function msd(R1::AbstractVector{Vector3D}, R2::AbstractVector{Vector3D},
		W::AbstractVector{<:Real})
	@boundscheck begin
		length(R1) > 0 || error("deviation of zero positions is undefined")
		length(R1) == length(R2) == length(W) ||
				error("size mismatch between position/weight arrays")
	end
	sd = 0.0
	meanW = mean(W)
	for i in eachindex(W)
		R1i = R1[i]
		R2i = R2[i]
		Wi = W[i]
		sd += (R2i.x - R1i.x)^2 * Wi / meanW
		sd += (R2i.y - R1i.y)^2 * Wi / meanW
		sd += (R2i.z - R1i.z)^2 * Wi / meanW
	end
	sd / length(W)
end

struct BoundingBox
	min::Vector3D
	max::Vector3D

	BoundingBox(min::AbstractVector{<:Real}, max::AbstractVector{<:Real}) =
			new(min, max)
end

function extent(R::AbstractVector{Vector3D})
	@boundscheck length(R) > 0 || error("extent of zero positions is undefined")
	xmin, ymin, zmin = Inf, Inf, Inf
	xmax, ymax, zmax = -Inf, -Inf, -Inf
	for Ri in R
		xmin = min(xmin, Ri.x)
		ymin = min(ymin, Ri.y)
		zmin = min(zmin, Ri.z)
		xmax = max(xmax, Ri.x)
		ymax = max(ymax, Ri.y)
		zmax = max(zmax, Ri.z)
	end
	BoundingBox(Vector3D(xmin, ymin, zmin), Vector3D(xmax, ymax, zmax))
end

extent(R::AbstractVector{<:Real}...) = extent([Vector3D(Ri) for Ri in R])

extent(bb::BoundingBox) = bb

Base.minimum(bb::BoundingBox) = bb.min

Base.maximum(bb::BoundingBox) = bb.max

dims(bb::BoundingBox) = maximum(bb) - minimum(bb)

Base.extrema(bb::BoundingBox) = minimum(bb), maximum(bb)

center(R1::Vector3D, R2::Vector3D) = R1 + (R2 - R1) / 2

center(bb::BoundingBox) = center(minimum(bb), maximum(bb))

volume(bb::BoundingBox) = prod(dims(bb))

isinside(R::AbstractVector{<:Real}, bb::BoundingBox) = isinside(Vector3D(R), bb)

isinside(R::Vector3D, bb::BoundingBox) =
		! (R.x < minimum(bb).x || R.y < minimum(bb).y || R.z < minimum(bb).z ||
		R.x > maximum(bb).x || R.y > maximum(bb).y || R.z > maximum(bb).z)

isinside(bb1::BoundingBox, bb2::BoundingBox) =
		isinside(minimum(bb1), bb2) && isinside(maximum(bb1), bb2)

end # module


text:
mutable struct SchellingAgent <: AbstractAgent
  id::Int # The identifier number of the agent
  pos::Tuple{Int,Int} # The x, y location of the agent
  mood::Bool # whether the agent is happy in its position. (true = happy)
  group::Int # The group of the agent,
             # determines mood as it interacts with neighbors
end

function instantiate_modelS(;numagents=320, griddims=(20, 20), min_to_be_happy=3)

  space = GridSpace(griddims, periodic = false)

  properties = Dict(:min_to_be_happy => min_to_be_happy)
  schelling = ABM(SchellingAgent, space; properties = properties, scheduler=random_activation)

  for n in 1:numagents
    agent = SchellingAgent(n, (1,1), false, n < numagents/2 ? 1 : 2)
    add_agent_single!(agent, schelling)
  end
  return schelling
end

function agent_step!(agent, model)
  agent.mood == true && return # do nothing if already happy
  minhappy = model.properties[:min_to_be_happy]
  # while agent.mood == false
      neighbor_positions = nearby_positions(agent, model)
      count_neighbors_same_group = 0
      # For each neighbor, get group and compare to current agent's group
      # and increment count_neighbors_same_group as appropriately.
      for neighbor_pos in neighbor_positions
          pos_contents = ids_in_position(neighbor_pos, model)
          # Skip iteration if the position is empty.
          length(pos_contents) == 0 && continue
          # Otherwise, get the first agent in the position...
          agent_id = pos_contents[1]
          # ...and increment count_neighbors_same_group if the neighbor's group is
          # the same.
          neighbor_agent_group = model.agents[agent_id].group
          if neighbor_agent_group == agent.group
              count_neighbors_same_group += 1
          end
      end

      # After counting the neighbors, decide whether or not to move the agent.
      # If count_neighbors_same_group is at least the min_to_be_happy, set the
      # mood to true. Otherwise, move the agent to a random position.
      if count_neighbors_same_group ≥ minhappy
          agent.mood = true
      else
          move_agent_single!(agent, model)
      end
  # end
  # return
end


text:
const AbstractNumOrArray = Union{AbstractVecOrMat,Number}
"""
    DescriptorStateSpace{T}(A::Matrix{T}, E::Union{Matrix{T},UniformScaling},
                            B::Matrix{T}, C::Matrix{T}, D::Matrix{T},
                            Ts::Real) where T <: Number

Construct a descriptor state-space model from a quintuple of matrices `(A,E,B,C,D)` and a sampling time `Ts`.

If `SYS::DescriptorStateSpace{T}` is a descriptor system model object
defined by the 4-tuple `SYS = (A-λE,B,C,D)`, then:

`SYS.A` is the `nx × nx` state matrix `A` with elements of type `T`.

`SYS.E` is the `nx × nx` descriptor matrix `E` with elements of type `T`.
 For a standard state-space system `SYS.E = I`, the `UniformScaling` of type `Bool`.

`SYS.B` is the `nx × nu` system input matrix `B` with elements of type `T`.

`SYS.C` is the `ny × nx` system output matrix `C` with elements of type `T`.

`SYS.D` is the `ny × nu` system feedthrough matrix `D` with elements of type `T`.

`SYS.Ts` is the real sampling time `Ts`, where `Ts = 0` for a continuous-time system,
    and `Ts > 0` or `Ts = -1` for a discrete-time system.
    `Ts = -1` indicates a discrete-time system with an unspecified sampling time.

The dimensions `nx`, `ny` and `nu` can be obtained as `SYS.nx`, `SYS.ny` and `SYS.nu`, respectively.
"""
struct DescriptorStateSpace{T} <: AbstractDescriptorStateSpace
    A::Matrix{T}
    E::Union{Matrix{T},UniformScaling}
    B::Matrix{T}
    C::Matrix{T}
    D::Matrix{T}
    Ts::Float64
    function DescriptorStateSpace{T}(A::Matrix{T}, E::Union{Matrix{T},UniformScaling},
                                     B::Matrix{T}, C::Matrix{T}, D::Matrix{T},  Ts::Real) where T
        dss_validation(A, E, B, C, D, Ts)
        new{T}(A, E, B, C, D, Float64(Ts))
    end
end
function dss_validation(A::Matrix{T}, E::Union{Matrix{T},UniformScaling},
                        B::Matrix{T}, C::Matrix{T}, D::Matrix{T},  Ts::Real) where T
    nx = LinearAlgebra.checksquare(A)
    (ny, nu) = size(D)

    # Validate dimensions
    if typeof(E) <: Matrix
       nx == LinearAlgebra.checksquare(E) || error("A and E must have the same size")
    end
    size(B, 1) == nx ||  error("B must have the same row size as A")
    size(C, 2) == nx ||  error("C must have the same column size as A")
    nu == size(B, 2) ||  error("D must have the same column size as B")
    ny == size(C, 1) ||  error("D must have the same row size as C")

    # Validate sampling time
    Ts >= 0 || Ts == -1 || error("Ts must be either a positive number, 0
                                (continuous system), or -1 (unspecified)")
    return nx, nu, ny
end
function Base.getproperty(sys::DescriptorStateSpace, d::Symbol)
    if d === :nx
        return size(sys.A,1)
    elseif d === :ny
        return size(sys.C,1)
    elseif d === :nu
        return size(sys.B,2)
    else
        getfield(sys, d)
    end
end
Base.propertynames(sys::DescriptorStateSpace) =
    (:nx, :ny, :nu, fieldnames(typeof(sys))...)

# Base.ndims(::AbstractDescriptorStateSpace) = 2
# Base.size(sys::AbstractDescriptorStateSpace) = size(sys.D)
# Base.length(sys::AbstractDescriptorStateSpace) = length(sys.D)
# Base.size(sys::AbstractDescriptorStateSpace, d::Integer) = d <= 2 ? size(sys)[d] : 1
Base.eltype(sys::DescriptorStateSpace) = eltype(sys.A)

function Base.getindex(sys::DST, inds...) where DST <: DescriptorStateSpace
    size(inds, 1) != 2 &&
        error("Must specify 2 indices to index descriptor state-space model")
    rows, cols = index2range(inds...)
    return DescriptorStateSpace{eltype(sys)}(copy(sys.A), copy(sys.E), sys.B[:, cols], sys.C[rows, :], sys.D[rows, cols], sys.Ts)
end
index2range(ind1, ind2) = (index2range(ind1), index2range(ind2))
index2range(ind::T) where {T<:Number} = ind:ind
index2range(ind::T) where {T<:AbstractArray} = ind
index2range(ind::Colon) = ind
function Base.lastindex(sys::DST, dim::Int) where DST <: DescriptorStateSpace
    lastindex(sys.D,dim)
end

# Basic Operations
function ==(sys1::DST1, sys2::DST2) where {DST1<:DescriptorStateSpace, DST2<:DescriptorStateSpace}
    # fieldnames(DST1) == fieldnames(DST2) || (return false)
    return all(getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(DST1))
end

function isapprox(sys1::DST1, sys2::DST2; atol = zero(real(eltype(sys1))),
                  rtol = rtol::Real =  ((max(size(sys1.A)...))+1)*eps(real(float(one(real(eltype(sys1))))))*iszero(atol)) where
                  {DST1<:DescriptorStateSpace,DST2<:DescriptorStateSpace}
    #fieldnames(DST1) == fieldnames(DST2) || (return false)
    return all(isapprox(getfield(sys1, f), getfield(sys2, f); atol = atol, rtol = rtol) for f in fieldnames(DST1))
end

# sum sys1+sys2
function +(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2}) where {T1,T2}
    #Ensure systems have same dimensions and sampling times
    size(sys1) == size(sys2) || error("The systems have different shapes.")
    Ts = promote_Ts(sys1.Ts,sys2.Ts)

    T = promote_type(T1, T2)
    n1 = sys1.nx
    n2 = sys2.nx
    A = [sys1.A  zeros(T,n1,n2);
         zeros(T,n2,n1) sys2.A]
    B = [sys1.B ; sys2.B]
    C = [sys1.C sys2.C;]
    D = [sys1.D + sys2.D;]
    if sys1.E == I && sys2.E == I
        E = I
    else
        E = [sys1.E  zeros(T,n1,n2);
             zeros(T,n2,n1) sys2.E]
    end

    return DescriptorStateSpace{T}(A, E, B, C, D, Ts)
end
# difference sys1-sys2
function -(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2}) where {T1,T2}
    #Ensure systems have same dimensions and sampling times
    size(sys1) == size(sys2) || error("The systems have different shapes.")
    Ts = promote_Ts(sys1.Ts,sys2.Ts)

    T = promote_type(T1, T2)
    n1 = sys1.nx
    n2 = sys2.nx

    A = [sys1.A  zeros(T,n1,n2);
         zeros(T,n2,n1) sys2.A]
    B = [sys1.B ; sys2.B]
    C = [sys1.C -sys2.C;]
    D = [sys1.D - sys2.D;]
    if sys1.E == I && sys2.E == I
        E = I
    else
        E = [sys1.E  zeros(T,n1,n2);
             zeros(T,n2,n1) sys2.E]
    end

    return DescriptorStateSpace{T}(A, E, B, C, D, Ts)
end
# negation -sys
function -(sys::DescriptorStateSpace{T}) where T
    return DescriptorStateSpace{T}(sys.A, sys.E, sys.B, -sys.C, -sys.D, sys.Ts)
end
# sys+mat and mat+sys
function +(sys::DescriptorStateSpace{T1}, mat::VecOrMat{T2}) where {T1,T2}
    p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
    size(sys) == (p, m) || error("The input-output dimensions of system does not match the shape of matrix.")
    T = promote_type(T1, T2)
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T), copy_oftype(sys.C,T),
                                   copy_oftype(sys.D,T) + mat, sys.Ts)
end
+(mat::VecOrMat{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(sys,mat)
# sys+I and I+sys
function +(sys::DescriptorStateSpace{T1}, mat::UniformScaling{T2}) where {T1,T2}
    size(sys,1) == size(sys,2) || error("The system must have the same number of inputs and outputs")
    T = promote_type(T1, T2)
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T), copy_oftype(sys.C,T),
                                   copy_oftype(sys.D,T) + mat, sys.Ts)
end

# sys-mat and mat-sys
function -(sys::DescriptorStateSpace{T1}, mat::VecOrMat{T2}) where {T1,T2}
    p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
    size(sys) == (p, m) || error("The input-output dimensions of system does not match the shape of matrix.")
    T = promote_type(T1, T2)
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T), copy_oftype(sys.C,T),
                                   copy_oftype(sys.D,T) - mat, sys.Ts)
end
-(mat::VecOrMat{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(-sys,mat)

# sys-I and I-sys
function -(sys::DescriptorStateSpace{T1}, mat::UniformScaling{T2}) where {T1,T2}
    size(sys,1) == size(sys,2) || error("The system must have the same number of inputs and outputs")
    T = promote_type(T1, T2)
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T), copy_oftype(sys.C,T),
                                   copy_oftype(sys.D,T) - mat, sys.Ts)
end
+(mat::UniformScaling{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(sys,mat)
-(mat::UniformScaling{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2} = +(-sys,mat)
# sys ± n and n ± sys
function +(sys::DescriptorStateSpace{T1}, n::Number) where T1
    T = promote_type(T1, eltype(n))
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T), sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T), copy_oftype(sys.C,T),
                                   copy_oftype(sys.D,T) .+ n, sys.Ts)
end
+(n::Number, sys::DescriptorStateSpace{T1}) where T1 = +(sys, n)
-(n::Number, sys::DescriptorStateSpace{T1}) where T1 = +(-sys, n)
-(sys::DescriptorStateSpace{T1},n::Number) where T1 = +(sys, -n)

# multiplication sys1*sys2
function *(sys1::DescriptorStateSpace{T1}, sys2::DescriptorStateSpace{T2}) where {T1,T2}
    sys1.nu == sys2.ny || error("sys1 must have same number of inputs as sys2 has outputs")
    Ts = promote_Ts(sys1.Ts, sys2.Ts)
    T = promote_type(T1, T2)
    n1 = sys1.nx
    n2 = sys2.nx

    A = [sys1.A    sys1.B*sys2.C;
         zeros(T,n2,n1)   sys2.A]
    B = [sys1.B*sys2.D ; sys2.B]
    C = [sys1.C   sys1.D*sys2.C;]
    D = [sys1.D*sys2.D;]
    if sys1.E == I && sys2.E == I
        E = I
    else
        E = [sys1.E  zeros(T,n1,n2);
             zeros(T,n2,n1) sys2.E]
    end

    return DescriptorStateSpace{T}(A, E, B, C, D, Ts)
end
# sys*mat
function *(sys::DescriptorStateSpace{T1}, mat::VecOrMat{T2}) where {T1,T2}
    p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
    sys.nu == p || error("The input dimension of system does not match the number of rows of the matrix.")
    T = promote_type(T1, T2)
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T),
                                   sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T)*to_matrix(T,mat), copy_oftype(sys.C,T),
                                   copy_oftype(sys.D,T)*to_matrix(T,mat), sys.Ts)
end
# mat*sys
function *(mat::VecOrMat{T1}, sys::DescriptorStateSpace{T2}) where {T1,T2}
    p, m = typeof(mat) <: Vector ? (length(mat),1) : size(mat)
    sys.ny == m || error("The output dimension of system does not match the number of columns of the matrix.")
    T = promote_type(T1, T2)
    return DescriptorStateSpace{T}(copy_oftype(sys.A,T),
                                   sys.E == I ? I : copy_oftype(sys.E,T),
                                   copy_oftype(sys.B,T), to_matrix(T,mat)*copy_oftype(sys.C,T),
                                   to_matrix(T,mat)*copy_oftype(sys.D,T), sys.Ts)
end
# sI*sys
function *(s::Union{UniformScaling,Number}, sys::DescriptorStateSpace{T}) where T
    T1 = promote_type(eltype(s),T)
    return DescriptorStateSpace{T1}(copy_oftype(sys.A,T1), sys.E == I ? I : copy_oftype(sys.E,T1), copy_oftype(sys.B,T1),
                                   lmul!(s,copy_oftype(sys.C,T1)), lmul!(s,copy_oftype(sys.D,T1)), sys.Ts)

end
# sys*sI
function *(sys::DescriptorStateSpace{T},s::Union{UniformScaling,Number}) where T
    T1 = promote_type(eltype(s),T)
    return DescriptorStateSpace{T1}(copy_oftype(sys.A,T1), sys.E == I ? I : copy_oftype(sys.E,T1), rmul!(copy_oftype(sys.B,T1),s),
                                   copy_oftype(sys.C,T1), rmul!(copy_oftype(sys.D,T1),s), sys.Ts)

end

# right division sys1/sys2 = sys1*inv(sys2)
/(n::Union{UniformScaling,Number}, sys::DescriptorStateSpace) = n*inv(sys)
/(sys::DescriptorStateSpace, n::Number) = sys*(1/n)
/(sys::DescriptorStateSpace, n::UniformScaling) = sys*(1/n.λ)

# left division sys1\sys2 = inv(sys1)*sys2
\(n::Number, sys::DescriptorStateSpace) = (1/n)*sys
\(n::UniformScaling, sys::DescriptorStateSpace) = (1/n.λ)*sys
\(sys::DescriptorStateSpace, n::Union{UniformScaling,Number}) = inv(sys)*n

# display sys
Base.print(io::IO, sys::DescriptorStateSpace) = show(io, sys)
Base.show(io::IO, sys::DescriptorStateSpace) = show(io, MIME("text/plain"), sys)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::DescriptorStateSpace)
    summary(io, sys); println(io)
    n = size(sys.A,1)
    typeof(sys.D) <: Vector ? ((p,m) = (length(sys.D),1)) : ((p, m) = size(sys.D))
    if n > 0
       println(io, "\nState matrix A:")
       show(io, mime, sys.A)
       if sys.E != I
          println(io, "\n\nDescriptor matrix E:")
          show(io, mime, sys.E)
       end
       m > 0 ? println(io, "\n\nInput matrix B:") : println(io, "\n\nEmpty input matrix B.")
       show(io, mime, sys.B)
       p > 0 ? println(io, "\n\nOutput matrix C:") : println(io, "\n\nEmpty output matrix C.")
       show(io, mime, sys.C)
       (m > 0 && p > 0) ? println(io, "\n\nFeedthrough matrix D:") : println(io, "\n\nEmpty feedthrough matrix D.")
       show(io, mime, sys.D)
       sys.Ts < 0 && println(io, "\n\nSample time: unspecified.")
       sys.Ts > 0 && println(io, "\n\nSample time: $(sys.Ts) second(s).")
       iszero(sys.Ts) ? println(io, "\n\nContinuous-time state-space model.") :
                        println(io, "Discrete-time state-space model.")
    elseif m > 0 && p > 0
       println(io, "\nFeedthrough matrix D:")
       show(io, mime, sys.D)
       println(io, "\n\nStatic gain.")
    else
       println(io, "\nEmpty state-space model.")
    end
end


text:
module AtomicOpt

using LinearAlgebra
using LinearMaps
using Printf
using Random
using SparseArrays
using Arpack
using IterativeSolvers
using Distributed
using FFTW
using LBFGSB
using BlockArrays

import Base: show, vec, size
import LinearAlgebra: dot, rank

export AbstractAtomicSet, AbstractAtom, AbstractFace
export MappedAtomicSet, MappedAtom, MappedFace
export SumAtomicSet, SumAtom, SumFace
export OneBall, OneBallAtom, OneBallFace
export NucBall, NucBallAtom, NucBallFace
export BlkNucBall, BlkNucBallAtom, BlkNucBallFace
# export PosSimplex, PosSimplexAtom, PosSimplexFace
# export TwoBall, TwoBallAtom, TwoBallFace
export TraceBall, TraceBallAtom, TraceBallFace
export face, face!, gauge, rank, support, expose, expose!
export face_project!
export level_set, conditional_graident, coordinate_descent
export Solution, constructPrimal
export MaskOP, TMaskOP
export dual_obj_gap

"Abstract atomic set."
abstract type AbstractAtomicSet end

"Abstract atoms."
abstract type AbstractAtom end

"Abstract face."
abstract type AbstractFace end

"Abstract operators."
abstract type AbstractOperator end

"Abstract Linear Operator"
LinearOp = Union{LinearMap, AbstractMatrix}
AbstractLinearOp = Union{LinearOp, AbstractOperator}

include("./BasicSets/OneBall.jl")
include("./BasicSets/NucBall.jl")
include("./BasicSets/BlkNucBall.jl")
include("./BasicSets/TraceBall.jl")
# include("src/BasicSets/PosSimplex.jl")
# include("src/BasicSets/TwoBall.jl")
include("./SetOperations/mapped.jl")
include("./SetOperations/sum.jl")
include("./SetOperations/operators.jl")
include("./SetOperations/facialprojection.jl")
include("./SetOperations/utils.jl")
include("./Optimization/solution.jl")
include("./Optimization/align.jl")
include("./Optimization/levelset.jl")
include("./Optimization/boxls.jl")
include("./Optimization/conditionalGradient.jl")
include("./Optimization/liftedCoordinateDescent.jl")

dot(a1::AbstractAtom, a2::Vector) = dot(vec(a1), a2)
dot(a1::Vector, a2::AbstractAtom) = dot(vec(a2), a1)
dot(a1::AbstractAtom, a2::AbstractAtom) = dot(vec(a1), vec(a2))

function Base.show(io::IO, A::AbstractAtomicSet)
    println(io, "atomic set  : ", atom_name(A))
    println(io, "description : ", atom_description(A))
    print(  io, "parameters  : ", atom_parameters(A))
end

function Base.show(io::IO, F::AbstractFace)
    println(io, "face of atomic set  : ", face_name(F))
    print(  io, "parameters          : ", face_parameters(F))
end

end # module


text:
using DataFrames
using Test

m = page_model()
include("../src/components/SeaLevelRise.jl")

SLR = add_comp!(m, SeaLevelRise)

set_param!(m, :SeaLevelRise, :rt_g_globaltemperature, readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv"))
set_param!(m, :SeaLevelRise, :y_year_0, 2008.)
set_param!(m, :SeaLevelRise, :y_year, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])

run(m)

es_equilibriumSL = m[:SeaLevelRise, :es_equilibriumSL]
es_equilibriumSL_compare = readpagedata(m, "test/validationdata/es_equilibriumSL.csv")
@test ones(10) ≈ es_equilibriumSL ./ es_equilibriumSL rtol=.01

s_sealevel = m[:SeaLevelRise, :s_sealevel]
s_sealevel_compare = readpagedata(m, "test/validationdata/s_sealevel.csv")
@test ones(10) ≈ s_sealevel ./ s_sealevel_compare rtol=.01

expfs_exponential = m[:SeaLevelRise, :expfs_exponential]
expfs_exponential_compare = readpagedata(m, "test/validationdata/expfs_exponential.csv")
@test ones(10) ≈ expfs_exponential ./ expfs_exponential_compare rtol=.01


text:
# This file defines the Attrition type. This type represents an attrition scheme
#   to make it possible for personnel members to leave service before their
#   mandatory retirement.

# The attrition type requires no additional types. [ XXX true? ]
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export Attrition
type Attrition
    # The name of the attrition scheme.
    name::String

    # The attrition period.
    attrPeriod::Float64

    # The time points of the attrition curve. This curve is always sorted, and
    #   the first element is 0.0.
    attrCurvePoints::Vector{Float64}

    # The attrition rates at the points of the attrition curve.
    attrRates::Vector{Float64}

    # The lambdas of each part of the attrition time distribution.
    lambdas::Vector{Float64}

    # An auxiliary vector of parameters.
    betas::Vector{Float64}

    # The exceeding probabilities for each time node of the attrition time
    #   distribution.
    gammas::Vector{Float64}

    function Attrition( name::String, rate::T1 = 0.0, period::T2 = 1.0 ) where T1 <: Real where T2 <: Real

        if ( rate < 0.0 ) || ( rate >= 1.0 )
            error( "Attrition rate must be between 0.0 and 1.0." )
        end  # if ( rate < 0.0 ) || ( rate > 100.0 )

        if period <= 0.0
            error( "Attrition period must be > 0.0." )
        end  # if period <= 0.0

        newAttr = new()
        newAttr.name = name
        # If there's no attrition, the period doesn't matter.
        newAttr.attrPeriod = rate == 0.0 ? 1.0 : period
        newAttr.attrCurvePoints = [ 0.0 ]
        newAttr.attrRates = [ rate ]
        computeDistPars( newAttr )
        return newAttr

    end  # Attrition( name, rate, period )

    function Attrition( rate::T1 = 0.0, period::T2 = 1.0 ) where T1 <: Real where T2 <: Real

        newAttr = Attrition( "default", rate, period )
        return newAttr

    end  # Attrition( rate, period )

    function Attrition( name::String, curve::T1, period::T2 = 1.0 ) where T1 <: Union{Dict{Float64, Float64}, Array{Float64, 2}} where T2 <: Real

        if period <= 0.0
            error( "Attrition period must be > 0.0." )
        end  # if period <= 0.0

        newAttr = Attrition( name )
        setAttritionCurve( newAttr, curve )
        setAttritionPeriod( newAttr, period )
        return newAttr

    end  # Attrition( name, curve, period )

    function Attrition( curve::T1, period::T2 = 1.0 ) where T1 <: Union{Dict{Float64, Float64}, Array{Float64, 2}} where T2 <: Real

        if period <= 0.0
            error( "Attrition period must be > 0.0." )
        end  # if period <= 0.0

        newAttr = Attrition( "default" )
        setAttritionCurve( newAttr, curve )
        setAttritionPeriod( newAttr, period )
        return newAttr

    end  # Attrition( curve, period )

end  # type Attrition


text:
struct ModelPosterior{M,A,O} <: AbstractConditionalModel{M,A,O}
    closure::ModelClosure{M,A}
    obs::O
end

model(post::ModelPosterior) = model(post.closure)

function Base.show(io::IO, cm::ModelPosterior)
    println(io, "ModelPosterior given")
    println(io, "    arguments    ", keys(argvals(cm)))
    println(io, "    observations ", keys(observations(cm)))
    println(io, model(cm))
end

type2model(::Type{MP}) where {M, MP<:ModelPosterior{M}} = type2model(M)
type2model(::ModelPosterior{M}) where {M} = type2model(M)

export argvals
argvals(c::ModelPosterior) = argvals(c.closure)

argvalstype(mp::ModelPosterior{M,A}) where {M,A} = A
argvalstype(::Type{MP}) where {M,A,MP<:ModelPosterior{M,A}} = A

obstype(mp::ModelPosterior{M,A,O}) where {M,A,O} = O
obstype(::Type{MP}) where {M,A,O,MP<:ModelPosterior{M,A,O}} = O

export observations
observations(c::ModelPosterior) = c.obs

export observed
function observed(cm::ModelPosterior{M,A,O}) where {M,A,O}
    keys(schema(O))
end

ModelPosterior(m::AbstractModel) = ModelPosterior(m,NamedTuple(), NamedTuple())

(post::ModelPosterior)(nt::NamedTuple) = ModelPosterior(model(post)(merge(argvals(post), nt)), post.obs)

Base.:|(post::ModelPosterior, nt::NamedTuple) = ModelPosterior(post.closure, merge(post.obs, nt))


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

# ## Setup
using Dojo

# ## Visualizer
vis= Visualizer()
open(vis)
set_camera!(vis,
    cam_pos=[-0.01, 0.0, 90.0],
    zoom=30)
set_floor!(vis,
    x=20.0,
    y=20.0,
    color=RGBA(1.0, 1.0, 1.0, 1.0))
set_light!(vis,
    ambient=0.65,
    direction="Positive")

################################################################################
# Nonlinear Friction Cone vs. Linearized Friction Cone
################################################################################
timestep = 0.01
gravity = -9.81
friction_coefficient = 0.25
x0 = [-1.5, -0.50, 0.25]
v0 = [4, 0.80, 0.0]
ω0 = [0.0, 0.0, 0.0]
opts = SolverOptions(rtol=1.0e-6, btol=1.0e-6)

# ## Linear cone
color_lc = orange;
mech_lc = get_mechanism(:block,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    contact_type=:linear,
    mode=:box,
    color=color_lc);

# ## Simulate
initialize!(mech_lc, :block,
    x=x0,
    q=one(Quaternion),
    v=v0,
    ω=ω0)

storage_lc = simulate!(mech_lc, 4.0,
    record=true,
    opts=opts)

# ## Visualize
via, anim = visualize(mech_lc, storage_lc,
    vis=vis,
    name=:lc)

# ## Nonlinear cone
color_nc = cyan;
mech_nc = get_mechanism(:block,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    contact_type=:nonlinear,
    mode=:box,
    color=color_nc);

# ## Simulate
initialize!(mech_nc, :block,
    x=x0,
    q=one(Quaternion),
    v=v0,
    ω=ω0)

storage_nc = simulate!(mech_nc, 4.0,
    record=true,
    opts=opts)

# ## Visualize
for (i, x) in enumerate(storage_nc.x[1])
    storage_nc.x[1][i] += [0.0; 0.0; 0.1]
end

vis, anim = visualize(mech_nc, storage_nc,
    vis=vis,
    name=:nc,
    animation=anim)

# ## MuJoCo pyramidal cone
color_mjlc = magenta;

mech_mjlc = get_mechanism(:block,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    contact_type=:linear,
    mode=:box, color=color_mjlc);

# ## Load
initialize!(mech_mjlc, :block,
    x=x0,
    q=one(Quaternion),
    v=v0,
    ω=ω0)
file = jldopen(joinpath(@__DIR__, "../mujoco_benchmark/results/cone_compare_pyramidal.jld2"))
storage_mjlc = generate_storage(mech_mjlc, [get_maximal_state(mech_mjlc), file["ztraj"]...])

# ## Visualize
vis, anim = visualize(mech_mjlc, storage_mjlc,
    vis=vis,
    name=:mjlc,
    animation=anim)

line_mat_mjlc = LineBasicMaterial(color=color_mjlc, linewidth=25.0)
points_mjlc = Vector{Point{3,Float64}}()
for (i, xt) in enumerate(storage_mjlc.x[1])
    k = xt
    push!(points_mjlc, Point(k[1], k[2], k[3]))
end
setobject!(vis[:path_mjlc], MeshCat.Line(points_mjlc, line_mat_mjlc))

# ## MuJoCo elliptic cone
color_mjnc = RGBA(0,0,0);

mech_mjnc = get_mechanism(:block,
    timestep=timestep,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    contact_type=:linear,
    mode=:box, color=color_mjnc);

# ## Load
initialize!(mech_mjnc, :block,
    x=x0,
    q=one(Quaternion),
    v=v0,
    ω=ω0)
file = jldopen(joinpath(@__DIR__, "../mujoco_benchmark/results/cone_compare_elliptic.jld2"))
storage_mjnc = generate_storage(mech_mjnc, [get_maximal_state(mech_mjnc), file["ztraj"]...])

# ## Visualize
vis, anim = visualize(mech_mjnc, storage_mjnc,
    vis=vis,
    name=:mjnc,
    animation=anim)


text:
import ACE.SphericalHarmonics: SHBasis, index_y

@doc raw"""
`struct Ylm1pBasis <: OneParticleBasis`

One-particle basis of the form
```math
\phi_{lm}({\bm r}) = Y_l^m(\hat{\br r})
```
Fundamental building block of ACE basis sets of the form
```math
   R_{nl}^{\mu_i, \mu_j}(r_{ij}) Y_l^m(\hat{\bm r})
```
This type basically just translates the `SHBasis` into a valid one-particle
basis.
"""
mutable struct Ylm1pBasis{T} <: OneParticleBasis{Complex{T}}
   SH::SHBasis{T}  # SH = Ylm
end

# TODO: allow for the possibility that the symbol where the
#       position is stored is not `rr` but something else!
#
# TODO: Should we drop this type altogether and just replace it with SHBasis?

# ---------------------- Implementation of Ylm1pBasis

Ylm1pBasis(maxL::Integer, T = Float64) = Ylm1pBasis(SHBasis(maxL, T))

Base.length(basis::Ylm1pBasis) = length(basis.SH)

# -> TODO : figure out how to do this well!!!
# Base.rand(basis::Ylm1pBasis) =
#       AtomState(rand(basis.zlist.list), ACE.Random.rand_vec(basis.J))

function get_spec(basis::Ylm1pBasis)
   spec = NamedTuple{(:n, :l), (Int, Int)}[]
   for l = 0:basis.SH.maxL, m = -l:l
      spec[index_y(l, m)] = Ylm1pBasisFcn(l, m)
   end
   return spec
end

==(P1::Ylm1pBasis, P2::Ylm1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::Ylm1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_Ylm1pBasis",
          "SH" => write_dict(basis.SH) )

read_dict(::Val{:ACE_Ylm1pBasis}, D::Dict) = Ylm1pBasis(read_dict(D["SH"]))

fltype(basis::Ylm1pBasis{T}) where T = Complex{T}
rfltype(basis::Ylm1pBasis{T}) where T = T

gradtype(::Ylm1pBasis{T}) where {T} = SVector{3, Complex{T}}

symbols(::Ylm1pBasis) = [:l, :m]

function indexrange(basis::Ylm1pBasis)
   maxl = basis.SH.maxL
   # note we create a stupid tensor product domain and then make up for it
   # by using an index filter during the basis generation process
   return Dict( :l => 0:maxl, :m => -maxl:maxl )
end

isadmissible(b, basis::Ylm1pBasis) =
      (  (0 <= b.l <= basis.SH.maxL) &&
         (-b.l <= b.m <= b.l) )

# ---------------------------  Evaluation code
#

alloc_B(basis::Ylm1pBasis, args...) = alloc_B(basis.SH)

alloc_dB(basis::Ylm1pBasis, args...) = alloc_dB(basis.SH)

alloc_temp(basis::Ylm1pBasis, args...) = alloc_temp(basis.SH)

alloc_temp_d(basis::Ylm1pBasis, args...) = alloc_temp_d(basis.SH, args...)

evaluate!(B, tmp, basis::Ylm1pBasis, X::AbstractState) =
      evaluate!(B, tmp, basis.SH, X.rr)

evaluate_d!(dB, tmpd, basis::Ylm1pBasis, X::AbstractState) =
      evaluate_d!(dB, tmpd, basis.SH, X.rr)

evaluate_ed!(B, dB, tmpd, basis::Ylm1pBasis, X::AbstractState) =
      evaluate_ed!(B, dB, tmpd, basis.SH, X.rr)

degree(b, basis::Ylm1pBasis) = b.l

get_index(basis::Ylm1pBasis, b) = index_y(b.l, b.m)

#
# function add_into_A_dA!(A, dA, tmpd, basis::RnYlm1pBasis, R, iz::Integer, iz0::Integer)
#    r = norm(R)
#    R̂ = R / r
#    # evaluate the r-basis and the R̂-basis for the current neighbour at R
#    evaluate_d!(tmpd.BJ, tmpd.dBJ, tmpd.tmpdJ, basis.J, r)
#    evaluate_d!(tmpd.BY, tmpd.dBY, tmpd.tmpdY, basis.SH, R)
#    # add the contributions to the A_zklm, ∇A
#    @inbounds for (i, nlm) in enumerate(basis.spec)
#       iY = index_y(nlm.l, nlm.m)
#       A[i] += tmpd.BJ[nlm.n] * tmpd.BY[iY]
#       dA[i] = (tmpd.dBJ[nlm.n] * tmpd.BY[iY]) * R̂ + tmpd.BJ[nlm.n] * tmpd.dBY[iY]
#    end
#    return nothing
# end


text:
using Colors

export noborder,
       dotted,
       dashed,
       solid,
       border,
       bordercolor,
       borderwidth,
       borderstyle,
       hline, vline,
       roundcorner,
       shadow,
       fillcolor

## Borders

const allsides = Side[]

@api bordercolor => (BorderColor <: Tile) begin
    doc(md"""Set the border color. Note that you need to set the `borderwidth` as
    well for the border to actually appear.""")
    typedarg(
        sides::AbstractArray=allsides,
        doc=md"""An array of sides to set the border color for. Valid values are
               `left`, `right`, `top` and `bottom`. By default the border color
               is set for all sides."""
    )
    arg(color::Color, doc="The color.")
    curry(tile::Tile, doc="The tile to border.")
end
render(t::BorderColor, state) =
    render(t.tile, state) & (
        mapparts(
            allsides, t.sides, "border", "Color", render_color(t.color)
        ) |> style
    )

@api borderwidth => (BorderWidth <: Tile) begin
    doc("Set the border width.")
    typedarg(
        sides::AbstractArray=allsides,
        doc=md"""An array of sides to set the border width for. Valid values are
               `left`, `right`, `top` and `bottom`. By default the border width
               is set for all sides."""
    )
    arg(width::Length, doc="The width.")
    curry(tile::Tile, doc="The tile to border.")
end
render(t::BorderWidth, state) =
    render(t.tile, state) &
        style(mapparts(allsides, t.sides, "border", "Width", t.width))

abstract StrokeStyle

@terms StrokeStyle begin
    noborder => NoStroke
    dotted => Dotted
    dashed => Dashed
    solid => Solid
end

name(::NoStroke) = "none"
name(::Solid) = "solid"
name(::Dotted) = "dotted"
name(::Dashed) = "dashed"

@api borderstyle => (BorderStyle <: Tile) begin
    doc(md"""Set the border style. Note that you need to set the `borderwidth`
             as well for the border to actually appear.""")
    typedarg(
        sides::AbstractArray=allsides,
        doc=md"""An array of sides to set the border style for. Valid values are
               `left`, `right`, `top` and `bottom`. By default the border style
               is set for all sides."""
        )
    arg(
        style::StrokeStyle,
        doc=md"""The border style. Valid values are `noborder`, `dotted`, `dashed`
               and `solid`."""
    )
    curry(tile::Tile, doc="The tile to border.")
end
render(t::BorderStyle, state) =
    render(t.tile, state) &
        style(mapparts(allsides, t.sides, "border", "Style", name(t.style)))

border(
    sides::AbstractArray,
    style::StrokeStyle,
    width::Length,
    color::Color,
    tile
) =
    borderstyle(sides, style,
        borderwidth(sides, width,
            bordercolor(sides, color, tile)))

border(sides::AbstractArray, style, width, color) =
    tile -> border(sides, style, width, color, tile)

border(style::StrokeStyle, width::Length, color::Color, tile) =
    borderstyle(allsides, style,
        borderwidth(allsides, width,
            bordercolor(allsides, color, tile)))

border(style::StrokeStyle, width, color) =
    tile -> border(style, width, color, tile)

@apidoc border => (BorderStyle <: Tile) begin
    doc("A helper function for setting border properties.")
    typedarg(sides::AbstractArray=allsides, doc="An array of sides to set the border for. Valid values are `left`, `right`, `top` and `bottom`. By default the border is set for all sides.")
    arg(style::StrokeStyle, doc="The border style. Valid values are `noborder`, `dotted`, `dashed` and `solid`.")
    arg(width::Length, doc="The width.")
    arg(color::Color, doc="The color.")
    curry(tile::Tile, doc="The tile to border.")
end
const default_border_color = RGB(0.6, 0.6, 0.6)

hline(;style=solid, width=1px, color=default_border_color) =
    border([bottom], style, width, color, empty)

@apidoc hline => (BorderStyle <: Tile) begin
    doc("Create a horizontal line. Returns a bordered tile of height 0.")
    kwarg(style::StrokeStyle, doc="The line style. Valid values are `noborder`, `dotted`, `dashed` and `solid`.")
    kwarg(width::Length, doc="The width.")
    kwarg(color::Color, doc="The color.")
end

vline(;style=solid, width=1px, color=default_border_color) =
    border([left], style, width, color, empty)

@apidoc vline => (BorderStyle <: Tile) begin
    doc("Create a vertical line. Returns a bordered tile of width 0.")
    kwarg(style::StrokeStyle, doc="The line style. Valid values are `noborder`, `dotted`, `dashed` and `solid`.")
    kwarg(width::Length, doc="The width.")
    kwarg(color::Color, doc="The color.")
end

## RoundRects

const allcorners = Corner[]
@api roundcorner => (RoundedRect <: Tile) begin
    doc("Round the corner of a tile.")
    typedarg(corners::AbstractArray=allcorners, doc=md"An array of corners to set the rounding for. By default, all corners are rounded. Valid corners are `topleft`,`midtop`,`topright`,`midleft`,`middle`,`midright`,`bottomleft`,`midbottom` and `bottomright`.")
    arg(radius::Length, doc="The radius.")
    curry(tile::Tile, doc="A tile.")
end

render(t::RoundedRect, state) =
    render(t.tile, state) &
        style(mapparts(allcorners, t.corners, "border", "Radius", t.radius))

## Box shadow

@api shadow => (Shadow <: Tile) begin
    doc("Show a shadow around a tile.")
    curry(tile::Tile, doc="A tile.")
    kwarg(inset::Bool=false, doc="Set whether the shadow falls inward or outward")
    kwarg(offset::(@compat Tuple{Length, Length})=(0px, 0px), doc="The offset. Displaces the shadow relative to the tile.")
    kwarg(blur_radius::Length=5px, doc="The radius for blurring.")
    kwarg(spread_radius::Length=5px, doc="The radius for spreading.")
    kwarg(color::Color=parse(Colorant,"black"), doc="Base color of the shadow.")
end

## Fill color

@api fillcolor => (FillColor <: Tile) begin
    doc("Fill a tile with a color.")
    arg(color::Color, doc="The color.")
    curry(tile::Tile, doc="The tile.")
end

render(t::FillColor, state) =
    render(t.tile, state) & style(@d(:backgroundColor => render_color(t.color)))