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

mutable struct ConnectionMonitorEndpoint <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    address::Any # spec type: Union{ Nothing, String } # spec name: address
    filter::Any # spec type: Union{ Nothing, ConnectionMonitorEndpointFilter } # spec name: filter

    function ConnectionMonitorEndpoint(;name=nothing, resourceId=nothing, address=nothing, filter=nothing)
        o = new()
        validate_property(ConnectionMonitorEndpoint, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ConnectionMonitorEndpoint, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(ConnectionMonitorEndpoint, Symbol("address"), address)
        setfield!(o, Symbol("address"), address)
        validate_property(ConnectionMonitorEndpoint, Symbol("filter"), filter)
        setfield!(o, Symbol("filter"), filter)
        o
    end
end # type ConnectionMonitorEndpoint

const _property_map_ConnectionMonitorEndpoint = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("resourceId")=>Symbol("resourceId"), Symbol("address")=>Symbol("address"), Symbol("filter")=>Symbol("filter"))
const _property_types_ConnectionMonitorEndpoint = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("resourceId")=>"String", Symbol("address")=>"String", Symbol("filter")=>"ConnectionMonitorEndpointFilter")
Base.propertynames(::Type{ ConnectionMonitorEndpoint }) = collect(keys(_property_map_ConnectionMonitorEndpoint))
Swagger.property_type(::Type{ ConnectionMonitorEndpoint }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorEndpoint[name]))}
Swagger.field_name(::Type{ ConnectionMonitorEndpoint }, property_name::Symbol) =  _property_map_ConnectionMonitorEndpoint[property_name]

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

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


text:
using Revise
using NNFEM
using PoreFlow
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function small_continuum_stiffness(k)
    small_continuum_stiffness_ = load_op_and_grad("./build/libSmallContinuumStiffness","small_continuum_stiffness", multiple=true)
    k = convert_to_tensor([k], [Float64]); k = k[1]
    ii, jj, vv = small_continuum_stiffness_(k)
    SparseTensor(ii+1, jj+1, vv, 2*domain.nnodes, 2*domain.nnodes)
end

m = 10
n = 10
h = 0.1
domain = example_domain(m, n, h)
# TODO: specify your input parameters
k = rand(domain.neles * length(domain.elements[1].weights),3,3)
q = copy(k)
for i = 1:4
    q[i,:,:] = k[i,:,:]'
end
init_nnfem(domain)
u = small_continuum_stiffness(k)
sess = Session(); init(sess)
@show run(sess, u)-run(sess, compute_fem_stiffness_matrix(constant(k), m, n, h))
@show

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    return sum(values(small_continuum_stiffness(m))^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(domain.neles * length(domain.elements[1].weights),3,3))
v_ = rand(domain.neles * length(domain.elements[1].weights),3,3)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

ext: jl


text:
pts = [rand(3) for i = 1:1000]

@test invariantmeasure(pts, RectangularBinning(3)) isa RectangularInvariantMeasure

#@test invariantmeasure(pts, InducedRectangularBinning(3, 4)) isa InducedRectangularInvariantMeasure
#@test invariantmeasure(pts, RectangularBinning(3), RectangularBinning(4)) isa InducedRectangularInvariantMeasure

#@test invariantmeasure(pts, AverageRectangularBinning(3, [4, 5])) isa AverageRectangularInvariantMeasure
#@test invariantmeasure(pts, RectangularBinning(3), [RectangularBinning(4), RectangularBinning(5)]) isa AverageRectangularInvariantMeasure


text:
using Test, bad

prior1 = Beta(mean = .4, sd = .1)
@test expectation(p -> 1, prior1) ≈ 1
@test (prior1.a ≈ 9.2) & (prior1.b ≈ 13.8)
@test mean(prior1) ≈ .4
prior_sd = expectation(p -> (p - mean(prior1))^2, prior1) |> sqrt
@test prior_sd ≈ .1

prior2 = condition(prior1; low = .3)
@test minimum(prior2.pivots) > .3
@test expectation(p -> 1, prior2) ≈ 1
@test (prior2.a ≈ 9.2) & (prior2.b ≈ 13.8)
@test (prior2.low ≈ .3) & (prior2.high ≈ 1.0)
@test mean(prior2) > mean(prior1)

prior3 = condition(prior1; high = .5)
@test maximum(prior3.pivots) < .5
@test expectation(p -> 1, prior3) ≈ 1
@test (prior3.a ≈ 9.2) & (prior3.b ≈ 13.8)
@test (prior3.low ≈ .0) & (prior3.high ≈ .5)
@test mean(prior3) < mean(prior1)

prior4 = Beta(2, 3)
@test expectation(p -> 1, prior4) ≈ 1
@test (prior4.a ≈ 2) & (prior4.b ≈ 3)
@test (prior4.low ≈ .0) & (prior4.high ≈ 1.)
@test mean(prior4) ≈ 2/5

post1 = update(prior4, 1, 1)
@test expectation(p -> 1, post1) ≈ 1
@test (post1.a ≈ 3) & (post1.b ≈ 3)
@test (post1.low ≈ .0) & (post1.high ≈ 1.)

post2 = update(prior3, 0, 1)
@test expectation(p -> 1, post2) ≈ 1
@test (post2.a ≈ 9.2) & (post2.b ≈ 14.8)
@test (post2.low ≈ .0) & (post2.high ≈ .5)

# mixtures

prior1  = condition(Beta(mean = .4, sd = .1), low = .3)
prior2  = Beta(1, 1)

mprior1 = .2prior1 + .3prior2
@test !is_proper(mprior1)
p = .01:.01:.99
@test all( .2*pdf.(p, prior1) + .3*pdf.(p, prior2) .== pdf.(p, mprior1) )

mprior2    = .9prior1 + .1prior2
@test mean(prior1) < mean(mprior2) < mean(prior2)
prior1_sd  = expectation(p -> (p - mean(prior1))^2, prior1) |> sqrt
mprior2_sd = expectation(p -> (p - mean(mprior2))^2, mprior2) |> sqrt
@test prior1_sd < mprior2_sd

mprior3 = .5prior1 + .5prior2
@test mean(mprior2) < mean(mprior3) < mean(prior2)
prior2_sd  = expectation(p -> (p - mean(prior2))^2, prior2) |> sqrt
mprior3_sd = expectation(p -> (p - mean(mprior3))^2, mprior3) |> sqrt
@test prior2_sd > mprior3_sd

mpost = update(mprior2, 1, 1)
@test (mpost.priors[1].a ≈ 10.2) & (mpost.priors[1].b ≈ 13.8)
@test (mpost.priors[2].a ≈ 2) & (mpost.priors[2].b ≈ 1)
@test mean(mpost) > mean(mprior2)

mpost2 = update(mprior2, 2, 8)
@test mean(mpost2) < mean(mprior2)
mpost2_sd = expectation(p -> (p - mean(mpost2))^2, mpost2) |> sqrt
@test mpost2_sd < mprior2_sd


text:
# KNITRO environment & context

#--------------------------------------------------
# KN_context
#--------------------------------------------------
"Wrapper for KNITRO KN_context."
mutable struct Env
    ptr_env::Ptr{Cvoid}

    function Env()
        ptrptr_env = Ref{Ptr{Cvoid}}()
        res = @kn_ccall(new, Cint, (Ptr{Cvoid},), ptrptr_env)
        if res != 0
            error("Fail to retrieve a valid KNITRO KN_context. Error $res")
        end
        new(ptrptr_env[])
    end

    Env(ptr::Ptr{Cvoid}) = new(ptr)
end

Base.unsafe_convert(ptr::Type{Ptr{Cvoid}}, env::Env) = env.ptr_env::Ptr{Cvoid}

is_valid(env::Env)= env.ptr_env != C_NULL

"""
Free all memory and release any Knitro license acquired by calling KN_new.
"""
function free_env(env::Env)
    if env.ptr_env != C_NULL
        ptrptr_env = Ref{Ptr{Cvoid}}(env.ptr_env)
        @kn_ccall(free, Cint, (Ptr{Cvoid},), ptrptr_env)
        env.ptr_env = C_NULL
    end
    return
end


text:
module PropCheck

using Base: Signed
using Test
using InteractiveUtils: subtypes

export generate, shrink
export check
export Integrated, Generator

# Trees
export root, subtrees, unfold, interleave
export igen

"""
    shrink(x)

Function to be used during shrinking of `x`. Returns an iterable of shrunk values.

Must _never_ return a previously input value, i.e. no value `x` used as input should ever be produced by `shrink(x)` or subsequent applications of `shrink` on the produced elements. This _will_ lead to infinite looping during shrinking.
"""
function shrink end

include("util.jl")
include("iteratorextras.jl")
include("config.jl")
include("generators.jl")
include("shrinkers.jl")
include("tree.jl")
include("integrated.jl")
include("checking.jl")

const igen = (Integrated ∘ Generator)

end # module


text:
using Plots, Distributions, LaTeXStrings, Statistics, StatsPlots, JLD, ProgressBars, StatsBase

function Metropolis(p :: Function, steps :: Integer, dx :: Real, X₀ :: Real)
    success=0
    x= X₀
    xList=Float64[i for i in 1:steps]
    n=0
    for n in 1:steps
        xList[n]=x
        y= x + rand(Uniform(-dx, dx))
        if rand()< p(y)/p(x)
            x=y
            success+=1
        end
    end
    return xList, success/steps
end

#histogram
xList, successful_steps=Metropolis(x -> ℯ^(-x^2/2), 10^5, 3 , 0.0)
save("../../chapter8/8.1/HistoData.jld","xList", xList, "SuccessfulSteps", successful_steps)
histogram(xList, normalize=true, alpha=0.7, label="Normal line")
plot!(Normal(0,1), ylabel="Probability Density", xlabel="x",
    title=L"Guassian\ Distribution\ for\ 100000\ steps\ and\ a_r=%$successful_steps", label="Metropolis Series", dpi=600)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter8\\HistoPlot.png")
#acceptence rate per unit length
N=100
stepsList= range(0.1, 50, length=N)
aᵣData=[]
for n in ProgressBar(1:N)
    xList, successful_steps=Metropolis(x -> ℯ^(-x^2/2), 10^5, stepsList[n] , 0.0)
    push!(aᵣData, successful_steps)
end
plot(stepsList,aᵣData , xlabel=L"dx", ylabel=L"a_r(dx)", title="Acceptance per Length of Steps", label=nothing, dpi=600)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter8\\acceptance_ratio.png")
#C(j)_j plot
stepsList=[15.9, 7.95, 5.27, 3.88, 2.94, 2.2, 1.57, 1.03, 0.5]
aᵣList=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
CⱼList=[]
for n in ProgressBar(1:9)
    xList, successful_steps=Metropolis(x -> ℯ^(-x^2/2), 10^5, stepsList[n] , 0.0)
    push!(CⱼList, autocor(xList))
end
save("../../chapter8/8.1/CjList.jld","CjList",CⱼList)
selectedNums=[1, 3, 6, 9 ]

plot()
for i in selectedNums
    plot!(1:length(CⱼList[i]), CⱼList[i], linewidth=3, label=L"a_r\ =\ %$(aᵣList[i])")
end
plot!(ylabe="Cⱼ", xlabel="j", title=L"C_j \_\ j", dpi=600)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter8\\Cj_j.png")
#length correlation
logCl=[]
ξList=Float64[]
x=[i for i in 1:10]
for i in CⱼList
    push!(logCl,log.(i[1:10]))
end
hcat(x, ones(10))
for y in logCl
    ξ⁽⁻¹⁾, Const= hcat(x, ones(10)) \ y
    ξ= -1/ξ⁽⁻¹⁾
    push!(ξList, ξ)
end
plot(aᵣList, ξList, title="Correlation Length", xlabel=L"a_r", ylabel=L"\xi", label=nothing)
scatter!(aᵣList, ξList, dpi=600, label=nothing)
save("../../chapter8/8.1/xiList.jld","ξList",ξList)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter8\\correlationLength.png")


text:
function _rec_get_fr(oct, lfs, areas, ndim, n_aux)

    if isa(oct, OctLeaf)
        push!(lfs, oct)
        push!(areas, leaf_area(oct, ndim))
    else
        for cinds in CartesianIndices(oct.subnodes)
            if cinds[ndim]==n_aux
                _rec_get_fr(oct.subnodes[cinds], lfs, areas,  ndim, n_aux)
            end
        end
    end

end

"""
Function to get vectors for all leaves in the frontier of an `AbstractOctNode`

* `oct`: `AbstractOctNode` at hand
* `ndim`: number of the searched dimension
* `n_aux`: 1 if looking for rear interfaces, 2 if otherwise

* return: a vector to all interfaced leaves and a vector with their respective interface area
"""
function get_frontier_interfaces(oct::AbstractOctNode, ndim::Int64, n_aux::Int64)

    lfs=Vector{AbstractOctNode}()
    areas=Vector{Float64}()

    _rec_get_fr(oct, lfs, areas, ndim, n_aux)

    return lfs, areas

end

"""
Function to get interfaces from neighboring leaves

* `oct`: `AbstractOctNode` at hand
* `ndim`: number of the searched dimension
* `n_aux`: 1 if looking for rear interfaces, 2 if otherwise

* return: a vector to all interfaced leaves and a vector with their respective interface area
"""
function get_neighbor_interfaces(oct::AbstractOctNode, ndim::Int64, n_aux::Int64)

    if !isnothing(oct.adjacent[ndim, n_aux])
        lfs, areas=get_frontier_interfaces(
            oct.adjacent[ndim, n_aux],
            ndim,
            (n_aux==1 ? 2 : 1)
        )

        if length(areas)==1
            areas[1]=leaf_area(oct, ndim)
        end

        return lfs, areas
    end

    lfs=Vector{AbstractOctNode}()
    areas=Vector{Float64}()

    return lfs, areas

end

function _rec_maxind(oct, prev)

    if isa(oct, OctLeaf)
        return (prev > oct.ind ? prev : oct.ind)
    end

    mx=prev

    for sd in oct.subnodes
        mx=_rec_maxind(sd, mx)
    end

    mx

end

function _rec_setlvint(oct, ndim, n_aux, m)

    if isa(oct, OctLeaf)
        lfs, areas=get_neighbor_interfaces(oct, ndim, n_aux)

        for (lf, a) in zip(lfs, areas)
            m[oct.ind, lf.ind]+=a
        end
    else
        for sd in oct.subnodes
            _rec_setlvint(sd, ndim, n_aux, m)
        end
    end

end

"""
Function to get interfaces to get rear or aft interfaces to leaves in the form of
a linear transformation.

**Pressuposes that `set_indexing!` has already been ran**

* `oct`: `AbstractOctNode` at hand
* `ndim`: number for the dimension for which the linear transformation is desired
* `n_aux`: 1 if expecting rear interfaces, 2 if otherwise
* `nlvs`: number of leaves, for matrix order. Surmised from leaves in `oct` if not provided

* return: a sparse matrix mapping values to `∫u_adj dS` surface integrals
"""
function neighbor_lintransf(oct::AbstractOctNode, ndim::Int64, n_aux::Int64;
    nlvs::Union{Int64, Nothing}=nothing)

    _nlvs=(isnothing(nlvs) ? _rec_maxind(oct, 0) : nlvs)

    m=spzeros(_nlvs, _nlvs)

    _rec_setlvint(oct, ndim, n_aux, m)

    m

end

function _rec_volset(oct, v)

    if isa(oct, OctLeaf)
        v[oct.ind]=leaf_volume(oct)
    else
        for sd in oct.subnodes
            _rec_volset(sd, v)
        end
    end

end

"""
Function to generate a vector of cell volumes through recursion through the tree

* `oct`: `AbstractOctNode` at hand
* `nlvs`: number of leaves, for vector order. Surmised from leaves in `oct` if not provided
"""
function volume_vector(oct::AbstractOctNode; nlvs::Union{Int64, Nothing}=nothing)

    _nlvs=(isnothing(nlvs) ? _rec_maxind(oct, 0) : nlvs)

    v=zeros(Float64, _nlvs)

    _rec_volset(oct, v)

    v

end

"""
Get a vector that implies on the additional fluxes to all cells in a given direction
when a unit Neumann condition is applied at a frontier

* `oct`: `AbstractOctNode` at hand
* `ndim`: number for the dimension for which the flux is desired
* `n_aux`: 1 if expecting rear interfaces, 2 if otherwise
* `nlvs`: number of leaves, for vector order. Surmised from leaves in `oct` if not provided
"""
function get_BC_flux(oct::AbstractOctNode, ndim::Int64, n_aux::Int64;
    nlvs::Union{Int64, Nothing}=nothing)

    _nlvs=(isnothing(nlvs) ? _rec_maxind(oct, 0) : nlvs)

    v=spzeros(_nlvs)

    lfs, areas=get_frontier_interfaces(oct, ndim, n_aux)

    for (lf, a) in zip(lfs, areas)
        v[lf.ind]=a
    end

    v

end

export get_frontier_interfaces, get_neighbor_interfaces, neighbor_lintransf, volume_vector, get_BC_flux


text:
# AbstractDict interface for JSON Object

Base.IteratorSize(::Type{JSON.Object{W, T}}) where W where T = Base.SizeUnknown()
Base.IteratorEltype(::Type{JSON.Object{W, T}}) where W where T = Base.EltypeUnknown()

# Access

function Base.getindex(o::JSON.Object, key)
    v = get(o, key, :not_found)
    if v === :not_found
        throw(KeyError(key))
    end
    return v
end

Base.get(o::JSON.Object, key, default) = _get(o, key, default)
function _get(o::JSON.Object{W}, key, default) where W
    i, c = get_ic(o, key, (0, 0x00))
    if i == 0
        return default
    end
    return getvalue(W, o.s, i, c)
end

# Dimensions

Base.length(a::JSON.Object) = div(collection_length(a), 2)

# Iterate

Base.iterate(j::JSON.Object, i = (j.i, 0x00)) = _iterate(j, i)
function _iterate(j::JSON.Object{W}, (i, c)) where W

    i, c = nextindex(j, i, c)

    if c == '}'
        return nothing
    end

    k = getvalue(W, j.s, i, c)
    i, c = nextindex(j, i, c)
    v = getvalue(W, j.s, i, c)

    return k => v, (i, c)
end

# Conversion to Dict Types

Base.convert(::Type{Any}, o::JSON.Object) = o

Base.convert(::Type{JSON.Object{W, T}},
                 o::JSON.Object{W, T}) where W where T <: AbstractString = o

Base.convert(::Type{T}, o::JSON.Object) where T <: AbstractDict = T(o)

# Conversion to Struct Types

#https://github.com/JuliaLang/julia/issues/26090
isreserved(s::Symbol) = s in (
    :while, :if, :for, :try, :return, :break, :continue,
    :function, :macro, :quote, :let, :local, :global, :const, :do,
    :struct,
    :type, :immutable, :importall,  # to be deprecated
    :module, :baremodule, :using, :import, :export,
    :end, :else, :catch, :finally, :true, :false)

unmangled_fieldnames(T) = [n[1] == '_' &&
                           isreserved(Symbol(n[2:end])) ?  n[2:end] : n
                           for n in map(Base.String, fieldnames(T))]

"""
Convert a `JSON.Object` to a `struct T`.
Optimised for the case where the order of the JSON fields matches the struct.
"""
@generated function Base.convert(::Type{T}, o::JSON.Object) where T

    fn = fieldnames(T)
    fk = unmangled_fieldnames(T)

    Expr(:block,
         :(    i = o.i                          ),
        (:(    (i, $n) = get_field(o, $k, i)    ) for (n,k) in zip(fn, fk))...,
         :(    $T($(fn...))                     ))
end

"""
Get a `field` from a `JSON.Object` starting at `start_i`.
Returns `start_i` for next field and field value.
"""
function get_field(o::JSON.Object{W}, field, start_i) where W
    i, c = get_ic(o, field, (0, 0x00), start_i)
    if i == 0
        #throw(KeyError(field))
        v = nothing
        i = start_i
    else
        v = getvalue(W, o.s, i, c)
        i = lastindex_of_value(o.s, i, c)
    end
    return i, v
end

# IOString Wrappers

Base.length(j::JSON.Object{W, IOString{T}}) where W where T =
    pump(() -> div(collection_length(j), 2), j.s)

Base.get(j::JSON.Object{W, IOString{T}}, key, default) where W where T =
    pump(() -> _get(j, key, default), j.s)

Base.iterate(j::JSON.Object{W, IOString{T}}, i = (j.i, 0x00)) where W where T =
    pump(() -> _iterate(j, i), j.s)


text:
import NLsolve

export CollocationResidual, solve_collocation

"""
Return an object that, when given parameters `θ`, construct a
parametric function using that an the family `fam`, and evaluates
`residual(model, f, x)` at the collocation points `x` of
`fam`. The two-argument version puts the residuals in the second
argument (should be a vector), while the first one returns them.

NOTE: Useful for passing to a solver such as `NLsolve.nlsolve`.
"""
immutable CollocationResidual{Tmodel, Tfam <: ParametricFamily, Txs, Tresidual} <: Function
    model::Tmodel
    fam::Tfam
    xs::Txs
    residual::Tresidual
    function CollocationResidual(model, fam, xs, residual)
        @argcheck degf(fam) == length(xs)
        new(model, fam, xs, residual)
    end
end

function CollocationResidual{Tmodel, Tfam, Txs, Tresidual}(model::Tmodel, fam::Tfam,
                                                           xs::Txs, residual::Tresidual)
    CollocationResidual{Tmodel, Tfam, Txs, Tresidual}(model, fam, xs, residual)
end

function CollocationResidual(model, fam, residual)
    CollocationResidual(model, fam, points(fam), residual)
end

@forward CollocationResidual.fam degf

function (cr::CollocationResidual)(θ, r)
    f = ParametricFunction(cr.fam, θ)
    r .= cr.residual.([cr.model], [f], cr.xs)
end

(cr::CollocationResidual)(θ) = cr(θ, Vector{eltype(cr.xs)}(degf(cr)))

function solve_collocation(cr::CollocationResidual, f₀; autodiff = true, options...)
    θ₀ = fit(cr.fam, f₀)
    o = NLsolve.nlsolve(cr, θ₀; autodiff = autodiff, options...)
    NLsolve.converged(o) || warn("optimization did not converge")
    ParametricFunction(cr.fam, o.zero), o
end


text:
using FractionalDiffEq, Plots
tspan = (0, 30); h = 0.01
rightfun(x, y) = 172/125*cos(4/5*x)
prob = MultiTermsFODEProblem([1, 1/16, 4/5, 3/2, 1/25, 6/5], [3, 2.5, 2, 1, 0.5, 0], rightfun, [0, 0, 0, 0, 0, 0], tspan)
sol = solve(prob, h, PIEX())
plot(sol, legend=:bottomright)


text:
# Need to patch up a put! method in order to deal with situation where
# conditions for the channel change while waiting for a lock.

using Base: check_channel_state

function Base.put_unbuffered(c::Channel, v)

    lk = c.cond_take.lock
    lock(lk)

    # relock to support clean exit from receive
    if lk != c.cond_take.lock
        @info "Relocked"

        unlock(lk)
        lk = c.cond_take.lock
        lock(lk)
    end

    taker = try
        while isempty(c.cond_take.waitq)
            check_channel_state(c)
            notify(c.cond_wait)
            wait(c.cond_put)
        end
        # unfair scheduled version of: notify(c.cond_take, v, false, false); yield()
        popfirst!(c.cond_take.waitq)
    finally
        unlock(lk)
    end
    schedule(taker, v)
    yield()  # immediately give taker a chance to run, but don't block the current task
    return v
end


text:
export FactorSize, Factor, VectorFactor, MatrixFactor
export factorsize, factorranks, factorndims
export factor, factormatrix, factorrankselect, block
export factorvcat, factorhcat, factordcat, factorltcat, factorutcat
export factorranktranspose, factormodetranspose
export factorcontract, factormp, factorkp, factorhp
export factorqr!, factorqradd, factorsvd!

const FactorSize = Vector{Int}
const Factor{T,N} = Array{T,N} where {T<:Number,N}
const VectorFactor{T} = Factor{T,3}
const MatrixFactor{T} = Factor{T,4}

function factorsize(U::Factor{T,N}) where {T<:Number,N}
	if length(U) == 0
		throw(ArgumentError("the factor should not be empty"))
	end
	sz = size(U)
	if length(sz) < 3
		throw(ArgumentError("the factor should have two rank dimensions and at least one mode dimension"))
	end
	return FactorSize(collect(sz[2:end-1]))
end

function factorranks(U::Factor{T,N}) where {T<:Number,N}
	if length(U) == 0
		throw(ArgumentError("the factor should not be empty"))
	end
	sz = size(U)
	if length(sz) < 3
		throw(ArgumentError("the factor should have two rank dimensions and at least one mode dimension"))
	end
	return sz[1],sz[end]
end

function factorndims(U::Factor{T,N}) where {T<:Number,N}
	if N < 3
		throw(ArgumentError("the factor should have two rank dimensions and at least one mode dimension"))
	end
	return N-2
end

function factor(U::Array{T,N}) where {T<:Number,N}
	U = reshape(U, (1,size(U)...,1))
	U = Factor{T,N+2}(U)
	return U
end

function factor(U::Matrix{T}, m::Union{Int,NTuple{M,Int},Vector{Int},Vector{Any}}, n::Union{Int,NTuple{N,Int},Vector{Int},Vector{Any}}, π::Union{NTuple{K,Int},Vector{Int}}) where {T<:Number,K,M,N}
	d = length(π)
	if d == 0
		throw(ArgumentError("π is empty"))
	end
	if isa(m, Vector{Any}) && length(m) > 0
		throw(ArgumentError("m should be an integer, a vector or tuple of integers or an empty vector or tuple"))
	end
	isa(m, Int) && (m = [m])
	isa(m, Vector{Any}) && (m = Vector{Int}(undef, 0))
	isa(m, Tuple{}) && (m = Vector{Int}(undef, 0))
	isa(m, Vector{Int}) || (m = collect(m))
	dm = length(m)
	if isa(n, Vector{Any}) && length(n) > 0
		throw(ArgumentError("n should be an integer, a vector or tuple of integers or an empty vector or tuple"))
	end
	isa(n, Int) && (n = [n])
	isa(n, Vector{Any}) && (n = Vector{Int}(undef, 0))
	isa(n, Tuple{}) && (n = Vector{Int}(undef, 0))
	isa(n, Vector{Int}) || (n = collect(n))
	dn = length(n)
	if dm+dn ≠ d || !isperm(π)
		throw(ArgumentError("π is not a valid permutation of 1:length(m)+length(n)"))
	end
	sz = size(U); p = sz[1]÷prod(m); q = sz[2]÷prod(n)
	if sz[1] ≠ prod(m)*p
		throw(DimensionMismatch("the first dimension of U is not divisible by prod(m)"))
	end
	if sz[2] ≠ prod(n)*q
		throw(DimensionMismatch("the second dimension of U is not divisible by prod(n)"))
	end
	U = reshape(U, (m...,p,n...,q))
	U = permutedims(U, (dm+1,π[1:dm]...,(π[dm+1:d].+1)...,d+2))
	return U
end

factor(U::Matrix{T}, m::Union{Int,NTuple{M,Int},Vector{Int},Vector{Any}}, n::Union{Int,NTuple{N,Int},Vector{Int},Vector{Any}}) where {T<:Number,M,N} = factor(U, m, n, collect(1:length(m)+length(n)))
factor(U::Vector{T}, m::Union{Int,NTuple{M,Int},Vector{Int},Vector{Any}}, π::Union{NTuple{M,Int},Vector{Int}}) where {T<:Number,M} = factor(U[:,:], m, (), π)
factor(U::Vector{T}, m::Union{Int,NTuple{M,Int},Vector{Int},Vector{Any}}) where {T<:Number,M} = factor(U[:,:], m, (), collect(1:length(m)))

function factormatrix(U::Factor{T,K}, π::Indices, σ::Indices) where {T<:Number,K}
	d = factorndims(U)
	if d == 0
		throw(ArgumentError("U should have at least one mode dimension"))
	end
	if isa(π, Vector{Any}) && length(π) > 0
		throw(ArgumentError("π should be an integer, a vector or tuple of integers, a colon or an empty vector or tuple"))
	end
	π = indvec(π; min=1, max=d); dπ = length(π)
	if isa(σ, Vector{Any}) && length(σ) > 0
		throw(ArgumentError("σ should be an integer, a vector or tuple of integers, a colon or an empty vector or tuple"))
	end
	σ = indvec(σ; min=1, max=d); dσ = length(σ)
	τ = vcat(π,σ)
	if dπ+dσ ≠ d || !isperm(τ)
		throw(ArgumentError("π and σ do not constitute a valid permutation of the mode dimensions of U"))
	end
	n = factorsize(U); p,q = factorranks(U)
	U = permutedims(U, ((π.+1)...,1,(σ.+1)...,d+2))
	U = reshape(U, p*prod(n[π]), q*prod(n[σ]))
	return U
end

function factorrankselect(U::Factor{T,N}, α::Indices, β::Indices) where {T<:Number,N}
	isa(α, Int) && (α = [α])
	isa(β, Int) && (β = [β])
	p,q = factorranks(U)
	α = indvec(α; min=1, max=p)
	β = indvec(β; min=1, max=q)
	if length(α) == 0
		throw(ArgumentError("the index or range for the first rank is empty"))
	end
	if α ⊈ 1:p
		throw(ArgumentError("the index or range for the first rank is incorrect"))
	end
	if length(β) == 0
		throw(ArgumentError("the index or range for the second rank is empty"))
	end
	if β ⊈ 1:q
		throw(ArgumentError("the index or range for the second rank is incorrect"))
	end
	d = N-2
	V = U
	V = selectdim(V, d+2, β)
	V = selectdim(V, 1, α)
	V = copy(V)
	return V
end

function block(U::Factor{T,N}, α::Int, β::Int) where {T<:Number,N}
	p,q = factorranks(U)
	if α ∉ 1:p
		throw(ArgumentError("the first rank index is out of range"))
	end
	if β ∉ 1:q
		throw(ArgumentError("the second rank index is out of range"))
	end
	d = N-2
	V = U
	V = selectdim(V, d+2, β)
	V = selectdim(V, 1, α)
	V = copy(V)
	return V
end

function factorvcat(U::Factor{T,N}, V::Factor{T,N}, W::Vararg{Factor{T,N},M}) where {T<:Number,N,M}
	m = factorsize(U); _,p = factorranks(U)
	W = (V,W...)
	for V ∈ W
		if factorsize(V) ≠ m
			throw(ArgumentError("the factors are incompatible in mode size"))
		end
		_,q = factorranks(V)
		if q ≠ p
			throw(ArgumentError("the factors are incompatible in the second rank"))
		end
	end
	return cat(U, W...; dims=1)
end

function factorhcat(U::Factor{T,N}, V::Factor{T,N}, W::Vararg{Factor{T,N},M}) where {T<:Number,N,M}
	m = factorsize(U); p,_ = factorranks(U)
	W = (V,W...)
	for V ∈ W
		if factorsize(V) ≠ m
			throw(ArgumentError("the factors are incompatible in mode size"))
		end
		q,_ = factorranks(V)
		if q ≠ p
			throw(ArgumentError("the factors are incompatible in the first rank"))
		end
	end
	d = N-2
	return cat(U, W...; dims=d+2)
end

function factordcat(U::Factor{T,N}, V::Factor{T,N}, W::Vararg{Factor{T,N},M}) where {T<:Number,N,M}
	m = factorsize(U)
	W = (V,W...)
	for V ∈ W
		if factorsize(V) ≠ m
			throw(ArgumentError("the factors are incompatible in mode size"))
		end
	end
	d = N-2
	return cat(U, W...; dims=(1,d+2))
end

function factorutcat(U₁₁::Factor{T,N}, U₁₂::Factor{T,N}, U₂₂::Factor{T,N}) where {T<:Number,N}
	n₁₁ = factorsize(U₁₁); p₁₁,q₁₁ = factorranks(U₁₁)
	n₁₂ = factorsize(U₁₂); p₁₂,q₁₂ = factorranks(U₁₂)
	n₂₂ = factorsize(U₂₂); p₂₂,q₂₂ = factorranks(U₂₂)
	if n₁₂ ≠ n₁₁ || n₂₂ ≠ n₁₁
		throw(ArgumentError("the factors are incompatible in mode size"))
	end
	if p₁₁ ≠ p₁₂
		throw(ArgumentError("U₁₁ and U₁₂ are incompatible in the first rank"))
	end
	if q₁₂ ≠ q₂₂
		throw(ArgumentError("U₁₂ and U₂₂ are incompatible in the second rank"))
	end
	U₂₁ = zeros(T, (p₂₂,n₁₁...,q₁₁))
	d = N-2
	return cat(cat(U₁₁, U₂₁; dims=1), cat(U₁₂, U₂₂; dims=1); dims=d+2)
end

function factorltcat(U₁₁::Factor{T,N}, U₂₁::Factor{T,N}, U₂₂::Factor{T,N}) where {T<:Number,N}
	n₁₁ = factorsize(U₁₁); p₁₁,q₁₁ = factorranks(U₁₁)
	n₂₁ = factorsize(U₂₁); p₂₁,q₂₁ = factorranks(U₂₁)
	n₂₂ = factorsize(U₂₂); p₂₂,q₂₂ = factorranks(U₂₂)
	if n₂₁ ≠ n₁₁ || n₂₂ ≠ n₁₁
		throw(ArgumentError("the factors are incompatible in mode size"))
	end
	if p₂₁ ≠ p₂₂
		throw(ArgumentError("U₂₁ and U₂₂ are incompatible in the first rank"))
	end
	if q₁₁ ≠ q₂₁
		throw(ArgumentError("U₁₁ and U₂₁ are incompatible in the second rank"))
	end
	U₁₂ = zeros(T, (p₁₁,n₁₁...,q₂₂))
	d = N-2
	return cat(cat(U₁₁, U₂₁; dims=1), cat(U₁₂, U₂₂; dims=1); dims=d+2)
end

function factorranktranspose(U::Factor{T,N}) where {T<:Number,N}
	d = N-2
	prm = vcat(d+2,(1:d).+1,1)
	U = permutedims(U, prm)
	return U
end

function factormodetranspose(U::Factor{T,N}, π::Union{NTuple{K,Int},Vector{Int}}) where {T<:Number,N,K}
	d = N-2
	if d == 0
		throw(ArgumentError("the factor should have at least one mode dimension"))
	end
	if length(π) ≠ d || !isperm(π)
		throw(ArgumentError("π is not a valid permutation of the mode dimensions of U"))
	end
	isa(π, Vector{Int}) || (π = collect(π))
	prm = vcat(1,π.+1,d+2)
	U = permutedims(U, prm)
	return U
end

factormodetranspose(U::Factor{T,N}) where {T<:Number,N} = factormodetranspose(U, collect(factorndims(U):-1:1))

function factorcontract(U::Factor{T,N}, V::Factor{T,N}; rev::Bool=false, major::String="last") where {T<:Number,N}
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	rev && ((U,V) = (V,U))
	m = factorsize(U); (p,r) = factorranks(U)
	n = factorsize(V); (s,q) = factorranks(V)
	if r ≠ s
		throw(ArgumentError("U and V have inconsistent ranks"))
	end
	d = length(m)
	if d == 0
		throw(ArgumentError("U has no mode dimensions"))
	end
	U = reshape(U,(p*prod(m),r)); V = reshape(V,(r,prod(n)*q))
	W = U*V; W = reshape(W, (p,m...,n...,q))
	prm = collect(1:2d); prm = reshape(prm, (d,2))
	(major == "first") && (prm = reverse(prm; dims=2))
	prm = prm'; prm = prm[:]
	W = permutedims(W, [1,(prm.+1)...,2*d+2])
	k = [m..., n...]; k = k[prm]; k = reshape(k, (2,d)); k = prod(k; dims=1)
	W = reshape(W, (p,k...,q))
	W = Factor{T}(W)
	return W
end

function factormp(U₁::Factor{T,N₁}, σ₁::Indices, U₂::Factor{T,N₂}, σ₂::Indices) where {T<:Number,N₁,N₂}
	n₁ = factorsize(U₁); d₁ = factorndims(U₁)
	n₂ = factorsize(U₂); d₂ = factorndims(U₂)
	if isa(σ₁, Vector{Any}) && length(σ₁) > 0
		throw(ArgumentError("if σ₁ is passed as a vector, it should be a vector of the type Vector{Int} or an empty vector of the type Vector{Any}"))
	end
	if isa(σ₂, Vector{Any}) && length(σ₂) > 0
		throw(ArgumentError("if σ₂ is passed as a vector, it should be a vector of the type Vector{Int} or an empty vector of the type Vector{Any}"))
	end
	isa(σ₁, Vector{Int}) || (σ₁ = indvec(σ₁; max=d₁))
	isa(σ₂, Vector{Int}) || (σ₂ = indvec(σ₂; max=d₂))
	if length(σ₁) ≠ length(σ₂)
		throw(ArgumentError("the specified sets of modes of σ₁ and σ₂ are inconsistent"))
	end
	if σ₁ ⊈ 1:d₁ || unique(σ₁) ≠ σ₁
		throw(ArgumentError("the set of modes of U₁ is specified incorrectly"))
	end
	if σ₂ ⊈ 1:d₂ || unique(σ₂) ≠ σ₂
		throw(ArgumentError("the set of modes of U₂ is specified incorrectly"))
	end
	if n₁[σ₁] ≠ n₂[σ₂]
		throw(ArgumentError("U₁ and U₂ are inconsistent with respect to the specified modes"))
	end
	τ₁ = setdiff(1:d₁, σ₁)
	τ₂ = setdiff(1:d₂, σ₂)
	(p₁,q₁) = factorranks(U₁); (p₂,q₂) = factorranks(U₂)
	U₁ = permutedims(U₁, (1,d₁+2,(τ₁.+1)...,(σ₁.+1)...))
	U₁ = reshape(U₁, p₁*q₁*prod(n₁[τ₁]), prod(n₁[σ₁]))
	U₂ = permutedims(U₂, [(σ₂.+1)...,(τ₂.+1)...,1,d₂+2])
	U₂ = reshape(U₂, prod(n₂[σ₂]), prod(n₂[τ₂])*p₂*q₂)
	U = U₁*U₂; n = [n₁[τ₁]..., n₂[τ₂]...]; d = length(τ₁) + length(τ₂)
	if d == 0
		n = [1]; d = 1
	end
	U = reshape(U, (p₁,q₁,n...,p₂,q₂))
	U = permutedims(U, (1,d+3,(3:d+2)...,2,d+4))
	U = reshape(U, (p₁*p₂,n...,q₁*q₂))
	return U
end

# function factorkp(U::Factor{T,N}, V::Vararg{Factor{T,N},M}) where {T<:Number,N,M}
# 	d = factorndims(U)
# 	nf = 1+length(V)
# 	n = Matrix{Int}(undef, d, nf)
# 	p = Vector{Int}(undef, nf)
# 	q = Vector{Int}(undef, nf)
# 	n[:,1] = factorsize(U); p[1],q[1] = factorranks(U)
# 	for k ∈ 2:nf
# 		n[:,k] = factorsize(V[k-1]); p[k],q[k] = factorranks(V[k-1])
# 		U = factormp(U, [], V[k-1], [])
# 	end
# 	prm = collect(1:nf*d); prm = reshape(prm, d, nf)
# 	prm = prm'; prm = prm[:]
# 	U = permutedims(U, [1,(prm.+1)...,nf*d+2])
# 	p,q = prod(p),prod(q); n = prod(n; dims=2)
# 	U = reshape(U, p, n..., q)
# 	return U
# end

# function factorkp2(U::Factor{T,N}, V::Factor{T,N}) where {T<:Number,N}
# 	d = factorndims(U)
# 	n = factorsize(U).*factorsize(V)
# 	p,q = factorranks(U).*factorranks(V)
# 	W = factormp(U, [], V, [])
# 	prm = collect(1:2*d); prm = reshape(prm, d, 2)
# 	prm = prm'; prm = prm[:]
# 	W = permutedims(W, [1,(prm.+1)...,2*d+2])
# 	W = reshape(W, p, n..., q)
# 	return W
# end

function factorkp(U::Union{Factor{T,N},Tuple{Factor{T,N},Int}}, V::Vararg{Union{Factor{T,N},Tuple{Factor{T,N},Int}},M}) where {T<:Number,N,M}
	V = (U,V...)
	nf = length(V)
	U = Vector{Factor{T,N}}(undef, nf)
	s = Vector{Int}(undef, nf)
	for k ∈ 1:nf
		W = V[k]
		if isa(W, Tuple)
			if W[2] < 0
				throw(ArgumentError("all the specified exponents should be nonnegative"))
			end
			U[k] = W[1]
			s[k] = W[2]
		else
			U[k] = W
			s[k] = 1
		end
	end
	m = findfirst(s .> 0)
	d = factorndims(U[m])
	W = U[m]
	n = factorsize(U[m])
	p,q = factorranks(U[m])
	for i ∈ 2:s[m]
		W = factormp(W, [], U[m], [])
		n = n.*factorsize(U[m])
		p,q = (p,q).*factorranks(U[m])
	end
	for k ∈ m+1:nf
		for i ∈ 1:s[k]
			W = factormp(W, [], U[k], [])
			n = n.*factorsize(U[k])
			p,q = (p,q).*factorranks(U[k])
		end
	end
	nf = sum(s)
	prm = collect(1:nf*d); prm = reshape(prm, d, nf)
	prm = prm'; prm = prm[:]
	W = permutedims(W, [1,(prm.+1)...,nf*d+2])
	W = reshape(W, p, n..., q)
	return W
end

function factorhp(U::Factor{T,N}, V::Factor{T,N}) where {T<:Number,N}
	m = factorsize(U); (p,q) = factorranks(U)
	n = factorsize(V); (r,s) = factorranks(V)
	if m ≠ n
		throw(ArgumentError("U and V should have the same mode size"))
	end
	m = prod(n); U = reshape(U, (p,m,q)); V = reshape(V, (r,m,s))
	W = Array{T,5}(undef, p, r, m, q, s)
	for i ∈ 1:m
		W[:,:,i,:,:] = reshape(kron(V[:,i,:], U[:,i,:]), (p,r,q,s))
	end
	W = reshape(W, (p*r,n...,q*s))
	return W
end

function factorqr!(U::Factor{T,N}; rev::Bool=false) where {T<:FloatRC{<:AbstractFloat},N}
	n = factorsize(U); (p,q) = factorranks(U); m = ones(Int, length(n))
	if rev
		U = reshape(U, p, prod(n)*q)
		R,U = lq!(U); R,U = Matrix(R),Matrix(U)
		ρ = size(R, 2)
		R,U = reshape(R, (p,m...,ρ)),reshape(U, (ρ,n...,q))
	else
		U = reshape(U, p*prod(n), q)
		U,R = qr!(U); U,R = Matrix(U),Matrix(R)
		ρ = size(R, 1)
		U,R = reshape(U, (p,n...,ρ)),reshape(R, (ρ,m...,q))
	end
	U = Factor{T,N}(U)
	R = Factor{T,N}(R)
	return U,R
end

factorqr!(U::Factor{T,N}, ::Val{false}; rev::Bool=false) where {T<:FloatRC{<:AbstractFloat},N} = factorqr!(U;  rev=rev)

function factorqr!(U::Factor{T,N}, ::Val{true}; rev::Bool=false) where {T<:FloatRC{<:AbstractFloat},N}
	n = factorsize(U); (p,q) = factorranks(U); m = ones(Int, length(n))
	if rev
		U = reshape(U, p, prod(n)*q)
		U = Matrix(U') # reallocation
		f = qr!(U, Val(true))
		U,R = Matrix(f.Q),Matrix(f.R*adjoint(f.P))
		R,U = Matrix(R'),Matrix(U')
		ρ = size(R, 2)
		R,U = reshape(R, (p,m...,ρ)),reshape(U, (ρ,n...,q))
	else
		U = reshape(U, p*prod(n), q)
		f = qr!(U, Val(true))
		U,R = Matrix(f.Q),Matrix(f.R*adjoint(f.P))
		ρ = size(R, 1)
		U,R = reshape(U, (p,n...,ρ)),reshape(R, (ρ,m...,q))
	end
	U = Factor{T,N}(U)
	R = Factor{T,N}(R)
	return U,R
end

function factorqradd(Q::Factor{T,N}, R::Union{Factor{T,N},Nothing}, U::Factor{T,N}; rev::Bool=false) where {T<:FloatRC{<:AbstractFloat},N}
	# assumes that Q is orthogonal w.r.t the first rank if rev==true and w.r.t the second rank if rev==false
	n = factorsize(Q); m = ones(Int, length(n))
	if factorsize(U) ≠ n
		throw(ArgumentError("Q and U are incompatible in mode size"))
	end
	if rev
		r,q = factorranks(Q); s,qq = factorranks(U)
		if isa(R, Nothing)
			R = reshape(Matrix{T}(I, r, r), (r,m...,r)); p = r
		else
			if factorsize(R) ≠ m
				throw(ArgumentError("R should have the same number of mode dimensions as Q and unitary mode sizes"))
			end
			p,rr = factorranks(R)
			if rr ≠ r
				throw(ArgumentError("R and Q are incompatible in rank"))
			end
		end
		if qq ≠ q
			throw(ArgumentError("Q and U are incompatible in the second rank"))
		end
		R = reshape(R, p, r)
		Q = reshape(Q, r, prod(n)*q)
		U = reshape(U, s, prod(n)*q)
		R,Q = lqaddrows(R, Q, U); r = size(R, 2)
		R,Q = reshape(R, (p+s,m...,r)),reshape(Q, (r,n...,q))
	else
		p,r = factorranks(Q); pp,s = factorranks(U)
		if isa(R, Nothing)
			R = reshape(Matrix{S}(I, r, r), (r,m...,r)); q = r
		else
			if factorsize(R) ≠ m
				throw(ArgumentError("R should have the same number of mode dimensions as Q and unitary mode sizes"))
			end
			rr,q = factorranks(R)
			if rr ≠ r
				throw(ArgumentError("Q and R are incompatible in rank"))
			end
		end
		if pp ≠ p
			throw(ArgumentError("Q and U are incompatible in the first rank"))
		end
		R = reshape(R, r, q)
		Q = reshape(Q, p*prod(n), r)
		U = reshape(U, p*prod(n), s)
		Q,R = qraddcols(Q, R, U); r = size(R, 1)
		Q,R = reshape(Q, (p,n...,r)),reshape(R, (r,m...,q+s))
	end
	Q,R = Factor{T,N}(Q),Factor{T,N}(R)
	return Q,R
end

"""
	factorsvd!(W, m, n; aTol=0, rTol=0, rank=0, major="last", rev=false)

produces U and V such that W ≈ U ⋈ V if rev == false and W ≈ V ⋈ U if rev == true
U is orthogonal — with respect to the second or first rank index
                  if rev == false or rev == true respectively
m and n are the mode-size vectors of U and V respectively (at most one of the two arguments may be replaced by ":"
major determines whether the "first" or "last" factor in the SKP carries the major bits
rank=0 means no rank thresholding
"""
function factorsvd!(W::Factor{T,N},
                    m::Union{FactorSize,Colon},
                    n::Union{FactorSize,Colon};
                    aTol::S=convert(S, 0),
                    rTol::S=convert(S, 0),
                    rank::Int=0,
                    major::String="last",
					rev::Bool=false)	where {S<:AbstractFloat,T<:FloatRC{S},N}
	d = factorndims(W)
	k = factorsize(W)
	if isa(m, Colon) && isa(n, Colon)
		throw(ArgumentError("to avoid ambiguity, at lease one of the two mode-size parameters should be specified as a vector"))
	end
	if isa(m, FactorSize)
		if any(m .≤ 0)
			throw(ArgumentError("the elements of the first mode-size parameter, when it is specified, should be positive"))
		end
		if length(m) ≠ d
			throw(DimensionMismatch("the number of entries in the first mode-size vector is inconsistent with the specified factor"))
		end
	else
		m = k.÷n
		if k ≠ m.*n
			throw(DimensionMismatch("not every mode dimension of W is divisible by the corresponding element of the specified mode-size vector"))
		end
	end
	if isa(n, FactorSize)
		if any(n .≤ 0)
			throw(ArgumentError("the elements of the second mode-size vector, when it is specified, should be positive"))
		end
		if length(n) ≠ d
			throw(DimensionMismatch("the number of entries in the second mode-size vector is inconsistent with the specified factor"))
		end
	else
		n = k.÷m
		if k ≠ m.*n
			throw(DimensionMismatch("not every mode dimension of W is divisible by the corresponding element of the specified mode-size vector"))
		end
	end
	if k ≠ m.*n
		throw(DimensionMismatch("the specified mode-size vectors are inconsistent with the specified factor"))
	end
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	#
	if aTol < 0 || !isfinite(aTol)
		throw(ArgumentError("aTol, when specified, should be nonnegative and finite"))
	end
	#
	aTol² = aTol.^2
	if !isfinite(aTol²)
		throw(ErrorException("overflow encountered while squaring aTol, which was passed finite"))
	end
	if aTol > 0 && aTol² == 0
		throw(ErrorException("underflow encountered while squaring aTol, which was passed positive"))
	end
	#
	if rTol < 0 || !isfinite(rTol)
		throw(ArgumentError("rTol, when specified, should be nonnegative and finite"))
	end
	#
	if rank < 0
		throw(ArgumentError("the optional argument rank should be nonnegative"))
	end
	ρ = rank
	#
	(p,q) = factorranks(W); prm = collect(1:d)
	k = rev ? [n; m] : [m; n]
	prm = (major == "last") ? [2*prm.-1; 2*prm] : [2*prm; 2*prm.-1]
	W = reshape(W, (p,k[invperm(prm)]...,q))
	W = permutedims(W, vcat(1, prm.+1, 2d+2))
	W = reshape(W, p*prod(k[1:d]), prod(k[d+1:2d])*q)
	U,Σ,V = svd!(W; full=false); V = V';
	σ = copy(Σ); μ = norm(σ)
	if !isfinite(μ)
		throw(ErrorException("overflow encountered while computing the norm of the decomposition"))
	end
	rTol = min(rTol, 1.0);
	τ = μ*rTol
	if μ > 0 && rTol > 0 && τ == 0
		throw(ErrorException("underflow encountered while computing the absolute accuracy threshold from the relative one"))
	end
	τ² = τ^2
	if μ > 0 && rTol > 0 && τ² == 0
		throw(ErrorException("underflow encountered while computing the squared absolute accuracy threshold from the squared relative one"))
	end
	τ² = [aTol²,τ²]; ind = (τ² .> 0)
	τ² = any(ind) ? minimum(τ²[ind]) : 0.0
	if μ == 0
		ε = 0.0; δ = 0.0; ρ = 1
		U = zeros(T, (p,k[1:d]...,1))
		V = zeros(T, (1,k[d+1:2d]...,q))
	else
		ε²,ρ = threshold(σ.^2, τ², ρ); ε = sqrt(ε²); δ = ε/μ
		ρ = max(ρ,1)
		U = U[:,1:ρ]; Σ = Σ[1:ρ]; V = V[1:ρ,:]
		rev ? U = U*Diagonal(Σ) : V = Diagonal(Σ)*V
		U = reshape(U, (p,k[1:d]...,ρ))
		V = reshape(V, (ρ,k[d+1:2d]...,q))
		rev && ((U,V) = (V,U))
	end
	U = Factor{T,N}(U)
	V = Factor{T,N}(V)
	return U,V,ε,δ,μ,ρ,σ
end


text:
import Base: conj, copy, real, imag
import LinearAlgebra: transpose, transpose!, adjoint!, adjoint
using FillArrays

# IMatrix
for func in (:conj, :real, :transpose, :adjoint, :copy)
    @eval ($func)(M::IMatrix{N, T}) where {N, T} = IMatrix{N, T}()
end
for func in (:adjoint!, :transpose!)
    @eval ($func)(M::IMatrix) = M
end
imag(M::IMatrix{N, T}) where {N, T} = Diagonal(zeros(T,N))

# PermMatrix
for func in (:conj, :real, :imag)
    @eval ($func)(M::PermMatrix) = PermMatrix(M.perm, ($func)(M.vals))
end
copy(M::PermMatrix) = PermMatrix(copy(M.perm), copy(M.vals))

function transpose(M::PermMatrix)
    new_perm = fast_invperm(M.perm)
    return PermMatrix(new_perm, M.vals[new_perm])
end

adjoint(S::PermMatrix{<:Real}) = transpose(S)
adjoint(S::PermMatrix{<:Complex}) = conj(transpose(S))

# scalar
import Base: *, /, ==, +, -, ≈
*(A::IMatrix{N, T}, B::Number) where {N, T} = Diagonal(Fill(promote_type(T, eltype(B))(B), N))
*(B::Number, A::IMatrix{N, T}) where {N, T} = Diagonal(Fill(promote_type(T, eltype(B))(B), N))
/(A::IMatrix{N, T}, B::Number) where {N, T} = Diagonal(Fill(promote_type(T, eltype(B))(1/B), N))

*(A::PermMatrix, B::Number) = PermMatrix(A.perm, A.vals*B)
*(B::Number, A::PermMatrix) = A*B
/(A::PermMatrix, B::Number) = PermMatrix(A.perm, A.vals/B)
#+(A::PermMatrix, B::PermMatrix) = PermMatrix(A.dv+B.dv, A.ev+B.ev)
#-(A::PermMatrix, B::PermMatrix) = PermMatrix(A.dv-B.dv, A.ev-B.ev)

const IDP = Union{Diagonal, PermMatrix, IMatrix}

for op in [:+, :-, :(==), :≈]

    @eval begin
        $op(A::IDP, B::SparseMatrixCSC) = $op(SparseMatrixCSC(A), B)
        $op(B::SparseMatrixCSC, A::IDP) = $op(B, SparseMatrixCSC(A))

        # intra-IDP
        $op(A::PermMatrix, B::IDP) = $op(SparseMatrixCSC(A), SparseMatrixCSC(B))
        $op(A::IDP, B::PermMatrix) = $op(SparseMatrixCSC(A), SparseMatrixCSC(B))
        $op(A::PermMatrix, B::PermMatrix) = $op(SparseMatrixCSC(A), SparseMatrixCSC(B))
    end

    # intra-ID
    if op in [:+, :-]
        @eval begin
            $op(d1::Diagonal, d2::IMatrix) = Diagonal($op(d1.diag, diag(d2)))
            $op(d1::IMatrix, d2::Diagonal) = Diagonal($op(diag(d1), d2.diag))
        end
    else
        @eval begin
            $op(d1::IMatrix, d2::Diagonal) = $op(diag(d1), d2.diag)
            $op(d1::Diagonal, d2::IMatrix) = $op(d1.diag, diag(d2))
            $op(d1::IMatrix{Na}, d2::IMatrix{Nb}) where {Na, Nb} = $op(Na, Nb)
        end
    end

end

# NOTE: promote 2 at least as an integer
+(d1::IMatrix{Na, Ta}, d2::IMatrix{Nb, Tb}) where {Na, Nb, Ta, Tb} = d1==d2 ? Diagonal(Fill(promote_type(Ta, Tb, Int)(2), Na)) : throw(DimensionMismatch())
-(d1::IMatrix{Na, Ta}, d2::IMatrix{Nb, Tb}) where {Na, Ta, Nb, Tb} = d1==d2 ? spzeros(promote_type(Ta, Tb), Na, Na) : throw(DimensionMismatch())


text:
@testset "RegressionBasedDID" begin
    hrs = exampledata("hrs")
    r = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated(11),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=true)
    # Compare estimates with Stata
    @test coef(r, "wave_hosp: 8 & rel: 0") ≈ 2825.5659 atol=1e-4
    @test coef(r, "wave_hosp: 8 & rel: 1") ≈ 825.14585 atol=1e-5
    @test coef(r, "wave_hosp: 8 & rel: 2") ≈ 800.10647 atol=1e-5
    @test coef(r, "wave_hosp: 9 & rel: -2") ≈ 298.97735 atol=1e-5
    @test coef(r, "wave_hosp: 9 & rel: 0") ≈ 3030.8408 atol=1e-4
    @test coef(r, "wave_hosp: 9 & rel: 1") ≈ 106.83785 atol=1e-5
    @test coef(r, "wave_hosp: 10 & rel: -3") ≈ 591.04639 atol=1e-5
    @test coef(r, "wave_hosp: 10 & rel: -2") ≈ 410.58102 atol=1e-5
    @test coef(r, "wave_hosp: 10 & rel: 0") ≈ 3091.5084 atol=1e-4

    @test coef(r) === r.coef
    @test vcov(r) === r.vcov
    @test vce(r) === r.vce
    @test treatment(r) == dynamic(:wave, -1)
    @test nobs(r) == 2624
    @test outcomename(r) == "oop_spend"
    @test coefnames(r) == ["wave_hosp: $w & rel: $r"
        for (w, r) in zip(repeat(8:10, inner=3), [0, 1, 2, -2, 0, 1, -3, -2, 0])]
    @test treatcells(r) == r.treatcells
    @test weights(r) === nothing
    @test ntreatcoef(r) == 9
    @test treatcoef(r) == coef(r)
    @test treatvcov(r) == vcov(r)
    @test treatnames(r) == coefnames(r)
    @test dof_residual(r) == 2610
    @test responsename(r) == "oop_spend"
    @test coefinds(r) == r.coefinds
    @test ncovariate(r) == 0

    @test has_lsweights(r)
    @test r.cellweights == r.cellcounts
    @test r.cellcounts == repeat([252, 176, 163, 65], inner=4)
    @test all(i->r.coef[i] ≈ sum(r.lsweights[:,i].*r.cellymeans), 1:ntreatcoef(r))

    @test has_fe(r)

    @test sprint(show, r) == "Regression-based DID result"
    pv = VERSION < v"1.6.0" ? " <1e-7" : "<1e-07"
    @test sprint(show, MIME("text/plain"), r) == """
        ──────────────────────────────────────────────────────────────────────
        Summary of results: Regression-based DID
        ──────────────────────────────────────────────────────────────────────
        Number of obs:               2624    Degrees of freedom:            14
        F-statistic:                 6.42    p-value:                   $pv
        ──────────────────────────────────────────────────────────────────────
        Cohort-interacted sharp dynamic specification
        ──────────────────────────────────────────────────────────────────────
        Number of cohorts:              3    Interactions within cohorts:    0
        Relative time periods:          5    Excluded periods:              -1
        ──────────────────────────────────────────────────────────────────────
        Fixed effects: fe_hhidpn fe_wave
        ──────────────────────────────────────────────────────────────────────
        Converged:                   true    Singletons dropped:             0
        ──────────────────────────────────────────────────────────────────────"""

    df = DataFrame(hrs)
    df.wave = settime(Date.(hrs.wave), Year(1))
    df.wave_hosp = settime(Date.(hrs.wave_hosp), Year(1), start=Date(7))
    r1 = @did(Reg, data=df, dynamic(:wave, -1), notyettreated(Date(11)),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=true)
    @test coef(r1) ≈ coef(r)
    @test r1.coefnames[1] == "wave_hosp: 0008-01-01 & rel: 0"

    rot = ifelse.(isodd.(hrs.hhidpn), 1, 2)
    df.wave = settime(Date.(hrs.wave), Year(1), rotation=rot)
    df.wave_hosp = settime(Date.(hrs.wave_hosp), Year(1), start=Date(7), rotation=rot)
    e = rotatingtime((1,2), Date(11))
    r2 = @did(Reg, data=df, dynamic(:wave, -1), notyettreated(e),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=true)
    @test length(coef(r2)) == 18
    @test coef(r2)[1] ≈ 3790.7218412450593
    @test r2.coefnames[1] == "wave_hosp: 1_0008-01-01 & rel: 0"

    r = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated([11]),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), cohortinteracted=false, lswtnames=(:wave_hosp, :wave))
    @test all(i->r.coef[i] ≈ sum(r.lsweights[:,i].*r.cellymeans), 1:ntreatcoef(r))

    @test sprint(show, MIME("text/plain"), r) == """
        ──────────────────────────────────────────────────────────────────────
        Summary of results: Regression-based DID
        ──────────────────────────────────────────────────────────────────────
        Number of obs:               2624    Degrees of freedom:             6
        F-statistic:                12.50    p-value:                   <1e-10
        ──────────────────────────────────────────────────────────────────────
        Sharp dynamic specification
        ──────────────────────────────────────────────────────────────────────
        Relative time periods:          5    Excluded periods:              -1
        ──────────────────────────────────────────────────────────────────────
        Fixed effects: none
        ──────────────────────────────────────────────────────────────────────"""

    r0 = @did(Reg, data=hrs, dynamic(:wave, -1), unspecifiedpr(),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave),),
        cohortinteracted=false, solvelsweights=true)
    # Compare estimates with Stata
    # gen rel = wave - wave_hosp
    # gen irel?? = rel==??
    # reghdfe oop_spend irel*, a(wave) cluster(hhidpn)
    @test coef(r0) ≈ [-1029.0482, 245.20926, 188.59266, 3063.2707,
        1060.5317, 1152.3315, 1986.7811] atol=1e-4
    @test diag(vcov(r0)) ≈ [764612.86, 626165.68, 236556.13, 163459.83,
        130471.28, 294368.49, 677821.36] atol=1e0
    pv = VERSION < v"1.6.0" ? " <1e-8" : "<1e-08"
    @test sprint(show, MIME("text/plain"), r0) == """
        ──────────────────────────────────────────────────────────────────────
        Summary of results: Regression-based DID
        ──────────────────────────────────────────────────────────────────────
        Number of obs:               3280    Degrees of freedom:            12
        F-statistic:                 9.12    p-value:                   $pv
        ──────────────────────────────────────────────────────────────────────
        Sharp dynamic specification
        ──────────────────────────────────────────────────────────────────────
        Relative time periods:          7    Excluded periods:              -1
        ──────────────────────────────────────────────────────────────────────
        Fixed effects: fe_wave
        ──────────────────────────────────────────────────────────────────────
        Converged:                   true    Singletons dropped:             0
        ──────────────────────────────────────────────────────────────────────"""

    sr = view(r, 1:3)
    @test coef(sr)[1] == r.coef[1]
    @test vcov(sr)[1] == r.vcov[1]
    @test parent(sr) === r

    tr = rescale(r, fill(2, 5), 1:5)
    @test coef(tr)[1] ≈ 2 * coef(sr)[1]
    @test vcov(tr)[1] ≈ 4 * vcov(sr)[1]
end

@testset "AggregatedRegDIDResult" begin
    hrs = exampledata("hrs")
    r = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated(11),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)))
    a = agg(r)
    @test coef(a) == coef(r)
    @test vcov(a) == vcov(r)
    @test vce(a) == vce(r)
    @test treatment(a) == dynamic(:wave, -1)
    @test nobs(a) == nobs(r)
    @test outcomename(a) == outcomename(r)
    @test coefnames(a) === a.coefnames
    @test treatcells(a) === a.treatcells
    @test weights(a) == weights(r)
    @test ntreatcoef(a) == 9
    @test treatcoef(a) == coef(a)
    @test treatvcov(a) == vcov(a)
    @test treatnames(a) == a.coefnames
    @test parent(a) === r
    @test dof_residual(a) == dof_residual(r)
    @test responsename(a) == "oop_spend"
    @test coefinds(a) === a.coefinds
    @test ncovariate(a) == 0

    @test !has_lsweights(a)

    pv = VERSION < v"1.6.0" ? "<1e-4 " : "<1e-04"
    @test sprint(show, a) === """
        ───────────────────────────────────────────────────────────────────────────────────
                                 Estimate  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
        ───────────────────────────────────────────────────────────────────────────────────
        wave_hosp: 8 & rel: 0    2825.57     1038.18   2.72    0.0065    789.825    4861.31
        wave_hosp: 8 & rel: 1     825.146     912.101  0.90    0.3657   -963.368    2613.66
        wave_hosp: 8 & rel: 2     800.106    1010.81   0.79    0.4287  -1181.97     2782.18
        wave_hosp: 9 & rel: -2    298.977     967.362  0.31    0.7573  -1597.9      2195.85
        wave_hosp: 9 & rel: 0    3030.84      704.631  4.30    $pv   1649.15     4412.53
        wave_hosp: 9 & rel: 1     106.838     652.767  0.16    0.8700  -1173.16     1386.83
        wave_hosp: 10 & rel: -3   591.046    1273.08   0.46    0.6425  -1905.3      3087.39
        wave_hosp: 10 & rel: -2   410.581    1030.4    0.40    0.6903  -1609.9      2431.06
        wave_hosp: 10 & rel: 0   3091.51      998.667  3.10    0.0020   1133.25     5049.77
        ───────────────────────────────────────────────────────────────────────────────────"""

    r = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated(11),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=true)
    # Compare estimates with Stata results from Sun and Abraham (2020)
    a = agg(r, :rel)
    @test coef(a, "rel: -3") ≈ 591.04639 atol=1e-5
    @test coef(a, "rel: -2") ≈ 352.63929 atol=1e-5
    @test coef(a, "rel: 0") ≈ 2960.0448 atol=1e-4
    @test coef(a, "rel: 1") ≈ 529.76686 atol=1e-5
    @test coef(a, "rel: 2") ≈ 800.10647 atol=1e-5

    @test has_lsweights(a)
    @test a.treatweights == a.treatcounts
    @test a.treatcounts == [163, 339, 591, 428, 252]
    @test all(i->a.coef[i] ≈ sum(a.lsweights[:,i].*r.cellymeans), 1:ntreatcoef(a))

    a1 = agg(r, (:rel,), subset=:rel=>isodd)
    @test length(coef(a1)) == 2
    @test coef(a1, "rel: -3") == coef(a, "rel: -3")
    @test coef(a1, "rel: 1") == coef(a, "rel: 1")
    @test sprint(show, a1) === """
        ───────────────────────────────────────────────────────────────────
                 Estimate  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
        ───────────────────────────────────────────────────────────────────
        rel: -3   591.046    1273.08   0.46    0.6425  -1905.3      3087.39
        rel: 1    529.767     586.831  0.90    0.3667   -620.935    1680.47
        ───────────────────────────────────────────────────────────────────"""

    a1 = agg(r, [:rel], bys=:rel=>isodd)
    @test length(coef(a1)) == 2
    @test coef(a1, "rel: false") ≈ sum(coef(a, "rel: $r") for r in ["-2", "0", "2"])
    @test coef(a1, "rel: true") ≈ sum(coef(a, "rel: $r") for r in ["-3", "1"])

    a2 = agg(r, :rel, bys=:rel=>isodd, subset=:rel=>isodd)
    @test coef(a2) == coef(a1, 2:2)

    @test_throws ArgumentError agg(r, subset=:rel=>x->x>10)
end

@testset "@specset" begin
    hrs = exampledata("hrs")
    # The first two specs are identical hence no repetition of steps should occur
    # The third spec should share all the steps until SolveLeastSquares
    # The fourth and fifth specs should not add tasks for MakeFEs and MakeYXCols
    # The sixth spec should not add any task for MakeFEs
    r = @specset [verbose] data=hrs yterm=term(:oop_spend) treatname=:wave_hosp begin
        @did(Reg, dynamic(:wave, -1), notyettreated(11),
            xterms=(fe(:wave)+fe(:hhidpn)))
        @did(Reg, dynamic(:wave, -1), notyettreated(11),
            xterms=[fe(:hhidpn), fe(:wave)])
        @did(Reg, dynamic(:wave, -2:-1), notyettreated(11),
            xterms=[fe(:hhidpn), fe(:wave)])
        @did(Reg, dynamic(:wave, -1), notyettreated(11),
            xterms=[term(:male), fe(:hhidpn), fe(:wave)])
        @did(Reg, dynamic(:wave, -1), notyettreated(11),
            treatintterms=TermSet(:male), xterms=[fe(:hhidpn), fe(:wave)])
        @did(Reg, dynamic(:wave, -1), nevertreated(11),
            xterms=(fe(:wave)+fe(:hhidpn)))
    end
    # Results might differ due to yxterms that include terms from other specs
    @test r[2][4] == didspec(Reg, dynamic(:wave, -1), notyettreated(11), data=hrs,
        yterm=term(:oop_spend), treatname=:wave_hosp, treatintterms=(),
        xterms=TermSet(term(:male), fe(:wave), fe(:hhidpn)))()

    r = @specset data=hrs yterm=term(:oop_spend) treatname=:wave_hosp begin
        @did(Reg, dynamic(:wave, -2:-1), nevertreated(11),
            xterms=(fe(:wave)+fe(:hhidpn)))
        @did(Reg, dynamic(:wave, -1), nevertreated(11),
            xterms=[fe(:hhidpn), fe(:wave)])
    end
    @test r[2][1] == @did(Reg, data=hrs, yterm=term(:oop_spend), treatname=:wave_hosp,
        dynamic(:wave, -2:-1), nevertreated(11), xterms=(fe(:wave)+fe(:hhidpn)))
end

@testset "contrast" begin
    hrs = exampledata("hrs")
    r1 = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated(11),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=true)
    r2 = agg(r1, :rel)
    c1 = contrast(r1)
    @test size(c1) == (16, 10)
    @test parent(c1) == [r1]
    c2 = contrast(r1, r2)
    @test size(c2) == (16, 15)
    @test c2[1] == r1.cellymeans[1]
    @test c2[1,2] == r1.lsweights[1]
    @test IndexStyle(typeof(c2)) == IndexLinear()

    @test c2.r == r1.lsweights.r == r2.lsweights.r
    @test c2.c.name[1] == "cellymeans"
    @test c2.c.name[2] == r1.coefnames[c2.c.icoef[2]]
    @test c2.c.iresult == [0, (1 for i in 1:9)..., (2 for i in 1:5)...]
    @test c2.c.icoef == [0, 1:9..., 1:5...]
    @test parent(c2)[1] === r1
    @test parent(c2)[2] === r2

    @test c1.r !== r1.lsweights.r
    sort!(c1, rev=true)
    @test c1.r == sort(r1.lsweights.r, rev=true)

    @test view(c1, :) == c1
    v = view(c1, :wave=>x->x==10)
    @test v.m == c1.m[c1.r.wave.==10,:]
    v = view(c1, 1:10)
    @test v.m == c1.m[1:10,:]
    @test v.r == view(c1.r, 1:10)
    v = view(c1, [:wave=>x->x==10, :wave_hosp=>x->x==10])
    @test v.m == c1.m[(c1.r.wave.==10).&(c1.r.wave_hosp.==10),:]
    @test v.r == view(c1.r, (c1.r.wave.==10).&(c1.r.wave_hosp.==10))

    c3 = contrast(r1, r2, subset=:wave=>7, coefs=1=>(:wave_hosp,2)=>(x,y)->x+y==8)
    @test size(c3) == (4, 8)
    @test parent(c3)[1] === r1
    @test c3.c.name[2] == "wave_hosp: 8 & rel: 0"
    @test c3.c.name[3] == "wave_hosp: 10 & rel: -2"
    @test c3.c.iresult == [0, (1 for i in 1:2)..., (2 for i in 1:5)...]
    @test c3.c.icoef == [0, 1, 8, 1:5...]
    @test all(c3.r.wave.==7)
    ri = r1.lsweights.r.wave.==7
    @test c3[:,1] == r1.cellymeans[ri]
    @test c3[:,2:end] == c2[ri,[2,9,11:15...]]

    c4 = contrast(r1, r2, subset=1:3, coefs=(1=>1:2, 2=>:))
    @test size(c4) == (3, 8)
    @test c4.c.iresult == [0, (1 for i in 1:2)..., (2 for i in 1:5)...]
    @test c4.c.icoef == [0:2..., 1:5...]
    @test c4.r == view(r1.lsweights.r, 1:3)
    @test c4[:,2:end] == c2[1:3,[2,3,11:15...]]
end

@testset "post!" begin
    hrs = exampledata("hrs")
    r1 = @did(Reg, data=hrs, dynamic(:wave, -1), notyettreated(11),
        vce=Vcov.cluster(:hhidpn), yterm=term(:oop_spend), treatname=:wave_hosp,
        treatintterms=(), xterms=(fe(:wave)+fe(:hhidpn)), solvelsweights=true)
    c1 = contrast(r1)
    gl = Dict{String,Any}()
    gr = Dict{String,Any}()
    gd = Dict{String,Any}()
    post!(gl, gr, gd, StataPostHDF(), c1, model="m")
    @test gl["model"] == "m"
    @test gl["depvar"] == "l_2"
    @test gl["b"] == r1.lsweights[:,1]
    @test gr["depvar"] == "r_3"
    @test gr["b"] == r1.lsweights[:,2]
    @test gd["depvar"] == "d_2_3"
    @test gd["b"] == (r1.lsweights[:,1].-r1.lsweights[:,2]).*r1.cellymeans
    cnames = string.(1:16)
    @test all(g->g["coefnames"]==cnames, (gl, gr, gd))
    @test all(g->!haskey(g, "at"), (gl, gr, gd))

    gl = Dict{String,Any}()
    gr = Dict{String,Any}()
    gd = Dict{String,Any}()
    post!(gl, gr, gd, StataPostHDF(), c1,
        lefttag="x", righttag="y", colnames=16:-1:1, at=1:16)
    @test gl["depvar"] == "l_x"
    @test gl["b"] == r1.lsweights[:,1]
    @test gr["depvar"] == "r_y"
    @test gr["b"] == r1.lsweights[:,2]
    cnames = string.(16:-1:1)
    @test all(g->g["coefnames"]==cnames, (gl, gr, gd))
    @test all(g->g["at"]==1:16, (gl, gr, gd))

    gl = Dict{String,Any}()
    gr = Dict{String,Any}()
    gd = Dict{String,Any}()
    post!(gl, gr, gd, StataPostHDF(), c1, eqnames=1:16)
    cnames = [string(i,":",i) for i in 1:16]
    @test all(g->g["coefnames"]==cnames, (gl, gr, gd))

    gl = Dict{String,Any}()
    gr = Dict{String,Any}()
    gd = Dict{String,Any}()
    post!(gl, gr, gd, StataPostHDF(), c1, eqnames=1:16, colnames=16:-1:1)
    cnames = [string(i,":",17-i) for i in 1:16]
    @test all(g->g["coefnames"]==cnames, (gl, gr, gd))

    gl = Dict{String,Any}()
    gr = Dict{String,Any}()
    gd = Dict{String,Any}()
    @test_throws ArgumentError post!(gl, gr, gd, StataPostHDF(), c1, eqnames=1:2)
    @test_throws ArgumentError post!(gl, gr, gd, StataPostHDF(), c1, colnames=1:2)
    @test_throws ArgumentError post!(gl, gr, gd, StataPostHDF(), c1, at=1:2)
end


text:
#=
    `relspeed` determines the datatype-relative performance of a function
    -  (speedof fn(data::TypeA)) / (speedof fn(data::TypeB))
    -  TypeB(TypeA(data)) == TypeB(data)
    -  TypeA(TypeB(data)) == TypeA(data)
    -  (speedof fn(xdata::TypeA, ydata::T)) / (speedof fn(xdata::TypeB, ydata::T)) where T
    -  TypeB(TypeA(data)) == TypeB(data)
    -  TypeA(TypeB(data)) == TypeA(data)
    -  (speedof fn(xdata::TypeA, ydata::TypeB)) / (speedof fn(xdata::TypeC, ydata::TypeD)) where ..
    -  TypeA(TypeC(xdata)) == TypeA(xdata) && TypeC(TypeA(xdata)) == TypeC(xdata)
    -  TypeB(TypeD(ydata)) == TypeB(ydata) && TypeD(TypeB(ydata)) == TypeD(ydata)
=#

function relspeed(fn, xs, ys; digits::Int=2, broadcasted::Bool=false)
    reltime = if broadcasted
                   relspeed_broadcast_fn(fn, xs, ys)
               else
                   relspeed_apply_fn(fn, xs, ys)
               end
    return round(Float64(reltime), digits=digits)
end

function relspeed(fn, x1s, x2s, y1s, y2s; digits::Int=2, broadcasted::Bool=false)
    reltime = if broadcasted
                  relspeed_broadcast_fn(fn, x1s, x2s, y1s, y2s)
              else
                  relspeed_apply_fn(fn, x1s, x2s, y1s, y2s)
              end
    return round(Float64(reltime), digits=digits)
end

function relspeed_apply_fn(fn, xs, ys)
    xstime = @refd @belapsed ($(fn))($(xs))
    ystime = @refd @belapsed ($(fn))($(ys))
    return Float32(xstime) / Float32(ystime)
end

function relspeed_broadcast_fn(fn, xs, ys)
    xstime = @refd @belapsed ($(fn)).($(xs))
    ystime = @refd @belapsed ($(fn)).($(ys))
    return Float32(xstime) / Float32(ystime)
end

function relspeed_apply_fn(fn, x1s, x2s, y1s, y2s)
    xstime = @refd @belapsed ($(fn))(($(x1s)), ($(x2s)))
    ystime = @refd @belapsed ($(fn))(($(y1s)), ($(y2s)))
    return Float32(xstime) / Float32(ystime)
end

function relspeed_broadcast_fn(fn, x1s, x2s, y1s, y2s)
    xstime = @refd @belapsed ($(fn)).(($(x1s)), ($(x2s)))
    ystime = @refd @belapsed ($(fn)).(($(y1s)), ($(y2s)))
    return Float32(xstime) / Float32(ystime)
end


text:
using DataFrames
using Gadfly

# Input parameters
path = "/home/cns/Downloads/"
abundanceFile = "benchmarks.csv"
plotname = "benchmarks.pdf"

plottheme = Gadfly.Theme(
                        default_point_size=3pt,
                        line_width=2pt,
                        minor_label_font_size=13pt,
                        minor_label_color=colorant"black",
                        major_label_font_size=15pt,
                        major_label_color=colorant"black",
                        key_label_color=colorant"black",
                        key_title_color=colorant"black",
                        key_title_font_size=15pt,
                        key_label_font_size=10pt,
                        grid_color=colorant"black",
                        highlight_width=0pt)

println("Reading input file")
results = readtable(
                  path*abundanceFile,
                  names = [:lang, :test, :time],
                  eltypes = [String, String, Float64],
                  header = false
                 )

p1 = plot(
          results,
           x="time",
           y="lang",
           color="test",
           Geom.point,
#           Geom.contour,
           Guide.xlabel("time"),
           Guide.ylabel("language"),
           Scale.x_log10,
            )
push!(p1,plottheme)

println("Now exporting plot")
draw(PNG("/home/cns/Dropbox/presentations/tuna_talk/benchmarks.png", 6inch, 5inch), p1)
println("Now ending script")


text:
# ---
# title: 938. Range Sum of BST
# id: problem938
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/range-sum-of-bst/description/>
# hidden: true
# ---
#
# Given the `root` node of a binary search tree, return _the sum of values of
# all nodes with a value in the range`[low, high]`_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/bst1.jpg)
#
#
#
#     Input: root = [10,5,15,3,7,null,18], low = 7, high = 15
#     Output: 32
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/05/bst2.jpg)
#
#
#
#     Input: root = [10,5,15,3,7,13,18,1,null,6], low = 6, high = 10
#     Output: 23
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 2 * 104]`.
#   * `1 <= Node.val <= 105`
#   * `1 <= low <= high <= 105`
#   * All `Node.val` are **unique**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file contains the code that allows things to be added to the registry

const ext2sym    = Dict{String, Union{Symbol,Vector{Symbol}}}()
const magic_list = Vector{Pair{Vector{UInt8},Symbol}}()     # sorted, see magic_cmp below
const sym2info   = Dict{Symbol,Tuple{Any,Any}}()   # Symbol=>(magic, extension)
const magic_func = Vector{Pair{Function,Symbol}}() # for formats with complex magic detection
const empty_magic = UInt8[]

## OS:
@enum OS Unix Windows OSX Linux

applies_to_os(oslist) = isempty(oslist) || any(applies_to_os, oslist)
function applies_to_os(os::OS)
    os == Unix && return Sys.isunix()
    os == Windows && return Sys.iswindows()
    os == OSX && return Sys.isapple()
    os == Linux && return Sys.islinux()
    return false
end

## Magic bytes:

# magic_cmp results in magic_list being sorted in order of increasing
# length(magic), then (among sequences with the same length) in
# lexographic order. This ordering has the advantage that you can
# incrementally read bytes from the stream without worrying that
# you'll encounter an EOF yet still have potential matches later in
# the list.
function magic_cmp(a::Vector{UInt8}, b::Vector{UInt8})
    la, lb = length(a), length(b)
    la < lb && return true
    la > lb && return false
    for (ia, ib) in zip(a, b)
        ia < ib && return true
        ia > ib && return false
    end
    return false
end
magic_cmp(p::Pair, m::Vector{UInt8}) = magic_cmp(p.first, m)
magic_cmp(m::Vector{UInt8}, p::Pair) = magic_cmp(m, p.first)

canonicalize_magic(@nospecialize(m::Tuple{Vararg{UInt8}})) = UInt8[m...]
canonicalize_magic(m::AbstractVector{UInt8}) = convert(Vector{UInt8}, m)
canonicalize_magic(m::String) = canonicalize_magic(codeunits(m))

## Load/Save

@enum IOSupport LOAD SAVE

function split_predicates(list)
    os = OS[]
    ls = IOSupport[]
    for item in list
        if isa(item, OS)
            push!(os, item)
        else
            push!(ls, item)
        end
    end
    return os, ls
end

function add_loadsave(format, predicates)
    library = popfirst!(predicates)
    os, loadsave = split_predicates(predicates)
    if applies_to_os(os)
        if isempty(loadsave) || (LOAD in loadsave)
            add_loader(format, library)
        end
        if isempty(loadsave) || (SAVE in loadsave)
            add_saver(format, library)
        end
    end
end

## Add Format:

function add_format(fmt, magic, extension, load_save_libraries...)
    for library in load_save_libraries
        add_loadsave(fmt, library)
    end
    # Add the format after we've validated the packages (to prevent a partially-registered format)
    add_format(fmt, magic, extension)
    fmt
end

"""
`add_format(fmt, magic, extension)` registers a new [`DataFormat`](@ref).
For example:

    add_format(format"TIFF", (UInt8[0x4d,0x4d,0x00,0x2b], UInt8[0x49,0x49,0x2a,0x00]), [".tiff", ".tif"])
    add_format(format"PNG", [0x89,0x50,0x4e,0x47,0x0d,0x0a,0x1a,0x0a], ".png")
    add_format(format"NRRD", "NRRD", [".nrrd",".nhdr"])

Note that extensions, magic numbers, and format-identifiers are case-sensitive.

You can also specify particular packages that support the format with `add_format(fmt, magic, extension, pkgspecifiers...)`,
where example `pkgspecifiers` are:

    add_format(fmt, magic, extension, [:PkgA=>UUID(...)])                     # only PkgA supports the format (load & save)
    add_format(fmt, magic, extension, [:PkgA=>uuidA], [:PkgB=>uuidB])         # try PkgA first, but if it fails try PkgB
    add_format(fmt, magic, extension, [:PkgA=>uuidA, LOAD], [:PkgB=>uuidB])   # try PkgA first for `load`, otherwise use PkgB
    add_format(fmt, magic, extension, [:PkgA=>uuidA, OSX], [:PkgB=>uuidB])    # use PkgA on OSX, and PkgB otherwise

The `uuid`s are all of type `UUID` and can be obtained from the package's `Project.toml` file.

You can combine `LOAD`, `SAVE`, `OSX`, `Unix`, `Windows` and `Linux` arbitrarily to narrow `pkgspecifiers`.
"""
add_format(@nospecialize(fmt::Type), args...) = add_format(formatname(fmt)::Symbol, args...)
add_format(sym::Symbol, magic::Union{Tuple,AbstractVector{UInt8},String}, extension) =
    add_format(sym, canonicalize_magic(magic), extension)
function add_format(sym::Symbol,
                    @nospecialize(magics::Tuple{Vector{UInt8},Vararg{Vector{UInt8}}}), extension)
    add_format(sym, [magics...], extension)
end

function add_format(sym::Symbol, magic::Vector{UInt8}, extension)
    haskey(sym2info, sym) && error("format ", sym, " is already registered")
    rng = searchsorted(magic_list, magic, lt=magic_cmp)
    if !isempty(magic) && !isempty(rng)
        error("magic bytes ", magic, " are already registered")
    end
    insert!(magic_list, first(rng), magic=>sym)
    sym2info[sym] = (magic, extension)
    add_extension(extension, sym)
    nothing
end

# for multiple magic bytes
function add_format(sym::Symbol, magics::Vector{Vector{UInt8}}, extension)
    haskey(sym2info, sym) && error("format ", sym, " is already registered")
    for magic in magics
        rng = searchsorted(magic_list, magic, lt=magic_cmp)
        if !isempty(magic) && !isempty(rng)
            error("magic bytes ", magic, " are already registered")
        end
        insert!(magic_list, first(rng), magic=>sym)
    end
    sym2info[sym] = (sort(magics; lt=magic_cmp), extension)
    add_extension(extension, sym)
    nothing
end

# For when "magic" is supplied as a function (see the HDF5 example in
# registry.jl)
function add_format(sym::Symbol, @nospecialize(magic::Function), extension)
    haskey(sym2info, sym) && error("format ", sym, " is already registered")
    push!(magic_func, Pair(magic,sym))  # magic=>sym in 0.4
    sym2info[sym] = (magic, extension)
    add_extension(extension, sym)
    nothing
end

"""
`del_format(fmt::DataFormat)` deletes `fmt` from the format registry.
"""
del_format(@nospecialize(fmt::Type)) = del_format(formatname(fmt)::Symbol)
function del_format(sym::Symbol)
    magic, extension = sym2info[sym]
    del_magic(magic, sym)
    delete!(sym2info, sym)
    del_extension(extension)
    nothing
end

# # Deletes multiple magic bytes
# del_magic(magic::Tuple, sym) = for m in magic
#     del_magic(m, sym)
# end
# Deletes single magic bytes
del_magic(@nospecialize(magic), sym::Symbol) = del_magic(canonicalize_magic(magic), sym)
function del_magic(magic::Vector{UInt8}, sym::Symbol)
    rng = searchsorted(magic_list, magic, lt=magic_cmp)
    if length(magic) == 0
        fullrng = rng
        found = false
        for idx in fullrng
            if last(magic_list[idx]) == sym
                rng = idx:idx
                found = true
                break
            end
        end
        found || error("format ", sym, " not found")
    end
    @assert length(rng) == 1
    deleteat!(magic_list, first(rng))
    nothing
end
del_magic(magics::Vector{Vector{UInt8}}, sym::Symbol) = foreach(magics) do magic
    del_magic(magic, sym)
end

function del_magic(@nospecialize(magic::Function), sym::Symbol)
    deleteat!(magic_func, something(findfirst(isequal(Pair{Function,Symbol}(magic,sym)), magic_func), 0))
    nothing
end

## File Extensions:

function add_extension(ext::String, sym::Symbol)
    if haskey(ext2sym, ext)
        v = ext2sym[ext]
        if isa(v, Symbol)
            ext2sym[ext] = Symbol[v, sym]
        else
            push!(ext2sym[ext], sym)
        end
        return
    end
    ext2sym[ext] = sym
end
function add_extension(ext::Union{Array,Tuple}, sym)
    for e in ext
        add_extension(e, sym)
    end
end

del_extension(ext::String) = delete!(ext2sym, ext)
function del_extension(ext::Union{Array,Tuple})
    for e in ext
        del_extension(e)
    end
end


text:
module JLLArrays

import ..JLLExceptions: InputException

"""
    move_dim(A, dim)

Permute A so that dimension `dim` is first:

```jldoctest
julia> A = zeros((4,3,2));
julia> A2 = move_dim(A, 2);
julia> size(A2)
(3, 4, 2)
```

    move_dim(A, dim, pos)

Permute A so that dimension `dim` ends up at position `pos`:

```jldoctest
julia> A = zeros((4,3,2));
julia> A2 = move_dim(A, 2, 3);
(4, 2, 3)
```

By default the returned array is a view into `A` (this saves memory but means
changes to values in one will affect the other). This can be changed by setting
the keyword `isview = false`:

```jldoctest
julia> A = zeros((4,3,2));
julia> Aview = move_dim(A, 2, isview=true);
julia> Acopy = move_dim(A, 2, isview=false);
julia> A[1] = 100;
julia> A[1] == Aview[1]
true

julia> A[1] == Acopy[1]
false
```
"""
function move_dim(A, dim, pos=1; isview=true)
    perm_vec = _compute_perm_vec(A, dim, pos);
    if isview
        return PermutedDimsArray(A, perm_vec);
    else
        return permutedims(A, perm_vec);
    end
end

"""
    _compute_perm_vec(A, dim, pos)

Helper function that computes the permutation vector to move `dim` to `pos`.
"""
function _compute_perm_vec(A, dim, pos)
    if dim < 1 || dim > ndims(A)
        throw(InputException("dim must be between 1 and ndims(A)"));
    elseif pos < 1 || pos > ndims(A)
        throw(InputException("pos must be between 1 and ndims(A)"));
    end

    perm_vec = [i for i in 1:ndims(A) if i != dim];
    splice!(perm_vec, pos:(pos-1), dim);
    return perm_vec;
end

"""
    seq_mat(dims...)
    seq_mat(T::DataType, dims...)

Create a matrix (optionally of the type `T`, defaults to Float64) with the
dimensions `dims` where each element's value is its linear index.
"""
function seq_mat(dims...)
    return seq_mat(Float64, dims...);
end

function seq_mat(T::DataType, dims...)
    mat = zeros(T, dims);
    for i = 1:length(mat)
        mat[i] = i;
    end
    return mat;
end

end

max_stars_repo_name: joshua-laughner/Julia-JLLUtils

max_issues_repo_name: joshua-laughner/Julia-JLLUtils

max_forks_repo_name: joshua-laughner/Julia-JLLUtils


text:
using DelimitedFiles, LightGraphs, Combinatorics
include("./misc.jl")
include("./graph.jl")

function construct_network_graph(instance_path::String)
    "In construct_network_g $instance_path"
    data = load_matpower(instance_path)
    println("Okay c passé")
    ## Bus load and shunt information
    i_debut, i_fin = find_numarray(1, data)
    checkfor(data, i_debut - 1, "mpc.bus")
    nb_bus = i_fin - i_debut + 1
    index_bus = Dict( data[i + i_debut - 1,1] => i for i=1:nb_bus)
    if index_bus != Dict( i => i for i=1:nb_bus)
        #println("!!! Specific numerotation of buses !!! \n")
    end
    ## Bus generator information
    i_debut, i_fin = find_numarray(i_fin + 1, data)
    checkfor(data, i_debut - 1, "mpc.gen")
    #initialize network g G
    G = SimpleGraph(nb_bus)
    ## Link information
    i_debut, i_fin = find_numarray(i_fin + 1, data)
    checkfor(data, i_debut - 1, "mpc.branch")
    for i=i_debut:i_fin
        if data[i, 11] == 0
        #@warn("link $(data[i,1])⟶$(data[i,2]) breaker out of service !")
        else
            orig = index_bus[data[i,1]]
            dest = index_bus[data[i,2]]
            LightGraphs.add_edge!(G, orig, dest)
        end
    end
    return G
end

function construct_graph_from_cliques(cliques::AbstractVector)
    nb_vertex = maximum([maximum(c) for c in cliques])
    g = SimpleGraph(nb_vertex)
    for clique in cliques
        make_complete!(g, clique)
    end
    return g
end

function construct_metagraph(g::AbstractGraph)
    mg = MetaGraph(g)
    for v in vertices(g)
        set_indexing_prop!(mg, v, :v, v)
    end
    return mg
end

function construct_graph_from_matrix(L)
    n = size(L,1)
    H = MetaGraph(n)
    for i in 1:n
        set_props!(H, i, Dict(:name => "node$i"))
    end
    for i in 1:n
        for j in 1:i
            if L[i,j] != 0
                MetaGraphs.add_edge!(H,i,j)
            end
        end
    end
    return H
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
"""
    struct MultiFieldFEFunction <: CellField
      # private fields
    end
"""
struct MultiFieldFEFunction{T<:MultiFieldCellField} <: FEFunction
  single_fe_functions::Vector{<:SingleFieldFEFunction}
  free_values::AbstractArray
  fe_space::MultiFieldFESpace
  multi_cell_field::T

  function MultiFieldFEFunction(
    free_values::AbstractVector,
    space::MultiFieldFESpace,
    single_fe_functions::Vector{<:SingleFieldFEFunction})

    multi_cell_field = MultiFieldCellField(map(i->i.cell_field,single_fe_functions))
    T = typeof(multi_cell_field)

    new{T}(
      single_fe_functions,
      free_values,
      space,
      multi_cell_field)
  end
end

CellData.get_data(f::MultiFieldFEFunction) = get_data(f.multi_cell_field)
CellData.get_triangulation(f::MultiFieldFEFunction) = get_triangulation(f.multi_cell_field)
CellData.DomainStyle(::Type{MultiFieldFEFunction{T}}) where T = DomainStyle(T)
FESpaces.get_free_dof_values(f::MultiFieldFEFunction) = f.free_values
FESpaces.get_fe_space(f::MultiFieldFEFunction) = f.fe_space

function FESpaces.get_cell_dof_values(f::MultiFieldFEFunction)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get the DOF values
  on top of the corresponding cells.
  """
  trians = map(get_triangulation,f.fe_space.spaces)
  trian = first(trians)
  @check all(map(t->have_compatible_domains(t,trian),trians)) msg
  get_cell_dof_values(f,trian)
end

function FESpaces.get_cell_dof_values(f::MultiFieldFEFunction,trian::Triangulation)
  function fun(uh)
    trian_i = get_triangulation(uh)
    if have_compatible_domains(trian_i,trian) ||
      have_compatible_domains(trian_i,get_background_triangulation(trian)) ||
      have_compatible_domains(trian_i,get_background_triangulation(get_background_triangulation(trian))) ||
      Geometry.is_included(trian,trian_i)
      cell_dofs = get_cell_dof_values(uh,trian)
    else
      cell_dofs_i = get_cell_dof_values(uh)
      T = eltype(eltype(cell_dofs_i))
      cell_dofs = Fill(T[],num_cells(trian))
    end
    cell_dofs
  end
  uhs = f.single_fe_functions
  blocks = [ fun(uh) for uh in uhs ]
  bids = [ (i,) for i in 1:length(uhs)]
  bsize = (length(uhs),)
  cell_axes_blocks = map(i->lazy_map(axes,i),blocks)
  cell_axes = lazy_map(_multifield_axes_dofs,cell_axes_blocks...)
  lazy_map(BlockArrayCooMap(bsize,bids),cell_axes,blocks...)
end

function FESpaces.get_cell_dof_values(f::MultiFieldFEFunction,trian::SkeletonTriangulation)
  cell_values_plus = get_cell_dof_values(f,trian.plus)
  cell_values_minus = get_cell_dof_values(f,trian.minus)
  cell_axes_plus = lazy_map(axes,cell_values_plus)
  cell_axes_minus = lazy_map(axes,cell_values_minus)
  cell_axes = lazy_map(cell_axes_plus,cell_axes_minus) do axp, axm
    (append_ranges([axp[1],axm[1]]),)
  end
  lazy_map(BlockArrayCooMap((2,),[(1,),(2,)]),cell_axes,cell_values_plus,cell_values_minus)
end

"""
    num_fields(m::MultiFieldFEFunction)
"""
num_fields(m::MultiFieldFEFunction) = length(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction) = iterate(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction,state) = iterate(m.single_fe_functions,state)
Base.getindex(m::MultiFieldFEFunction,field_id::Integer) = m.single_fe_functions[field_id]