Use ContinuumArrays.jl

Project.toml

@@ -1,13 +1,14 @@
 name = "ApproxFunBase"
 uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
-version = "0.6.6"
+version = "0.7"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+ContinuumArrays = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
@@ -18,13 +19,13 @@ InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
@@ -34,6 +35,7 @@ BandedMatrices = "0.16, 0.17"
 BlockArrays = "0.14, 0.15, 0.16"
 BlockBandedMatrices = "0.10, 0.11"
 Calculus = "0.5"
+ContinuumArrays = "0.10"
 DSP = "0.6, 0.7"
 DomainSets = "0.5"
 DualNumbers = "0.6.2"
@@ -44,15 +46,17 @@ InfiniteArrays = "0.11, 0.12"
 InfiniteLinearAlgebra = "0.5, 0.6"
 IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.20, 0.21, 0.22"
-LowRankApprox = "0.2, 0.3, 0.4, 0.5"
+LowRankApprox = "0.4, 0.5"
 SpecialFunctions = "0.10, 1.0, 2"
 StaticArrays = "0.12, 1.0"
 ToeplitzMatrices = "0.6, 0.7"
-julia = "1.5"
+julia = "1.7"
 
 [extras]
+ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 ApproxFunOrthogonalPolynomials = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ApproxFunOrthogonalPolynomials", "Random"]
+test = ["Random", "ApproxFunOrthogonalPolynomials", "ClassicalOrthogonalPolynomials", "Test"]

src/ApproxFunBase.jl

@@ -2,7 +2,8 @@ module ApproxFunBase
     using Base: AnyDict
 using Base, BlockArrays, BandedMatrices, BlockBandedMatrices, DomainSets, IntervalSets,
             SpecialFunctions, AbstractFFTs, FFTW, SpecialFunctions, DSP, DualNumbers,
-            LinearAlgebra, SparseArrays, LowRankApprox, FillArrays, InfiniteArrays, InfiniteLinearAlgebra #, Arpack
+            LinearAlgebra, SparseArrays, LowRankApprox, FillArrays, InfiniteArrays,
+            InfiniteLinearAlgebra, ContinuumArrays
 import StaticArrays, Calculus
 
 import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
@@ -26,7 +27,7 @@ import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !,
                 getproperty, findfirst, unsafe_getindex, fld, cld, div, imag,
                 @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                 Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
-                split, iszero, permutedims
+                split, iszero, permutedims, vec
 
 import Base.Broadcast: BroadcastStyle, Broadcasted, AbstractArrayStyle, broadcastable,
                         DefaultArrayStyle, broadcasted
@@ -40,6 +41,8 @@ import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, cross
 
 import SparseArrays: blockdiag
 
+
+import ContinuumArrays: AbstractQuasiMatrix, AbstractQuasiVector, AbstractQuasiArray, arguments, checkpoints
 # import Arpack: eigs
 
 # we need to import all special functions to use Calculus.symbolic_derivatives_1arg
@@ -78,6 +81,7 @@ import InfiniteArrays: PosInfinity, InfRanges, AbstractInfUnitRange, OneToInf, I
 
 # convenience for 1-d block ranges
 const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
+const Space{T} = AbstractQuasiMatrix{T}
 
 import Base: view
 
@@ -95,6 +99,7 @@ export pad!, pad, chop!, sample,
 
 export .., Interval, ChebyshevInterval, leftendpoint, rightendpoint, endpoints, cache
 
+export bilinearform, linebilinearform, innerproduct, lineinnerproduct
 
 if VERSION < v"1.6-"
 	oneto(n) = Base.OneTo(n)
@@ -103,17 +108,6 @@ else
 end
 
 
-include("LinearAlgebra/LinearAlgebra.jl")
 include("Fun.jl")
-include("Domains/Domains.jl")
-include("Multivariate/Multivariate.jl")
-include("Operators/Operator.jl")
-include("Caching/caching.jl")
-include("PDE/PDE.jl")
-include("Spaces/Spaces.jl")
-include("hacks.jl")
-include("testing.jl")
-include("specialfunctions.jl")
-include("show.jl")
 
 end #module

src/Fun.jl

@@ -2,97 +2,83 @@ export Fun, evaluate, values, points, extrapolate, setdomain
 export coefficients, ncoefficients, coefficient
 export integrate, differentiate, domain, space, linesum, linenorm
 
-include("Domain.jl")
-include("Space.jl")
-
 ##  Constructors
 
 
 mutable struct Fun{S,T,VT} <: Function
     space::S
     coefficients::VT
-    function Fun{S,T,VT}(sp::S,coeff::VT) where {S,T,VT}
-        @assert length(coeff) ≤ dimension(sp)
+    function Fun{S,T,VT}(sp::S, coeff::VT) where {S,T,VT}
+        axes(sp,2) == axes(coeff,1) || throw(DimensionMismatch(""))
         new{S,T,VT}(sp,coeff)
     end
 end
 
 const VFun{S,T} = Fun{S,T,Vector{T}}
 
-Fun(sp::Space,coeff::AbstractVector) = Fun{typeof(sp),eltype(coeff),typeof(coeff)}(sp,coeff)
+_pad(c, _, n) = [c; zeros(eltype(c), n-length(c))]
+_pad(c, ::PosInfinity, ::PosInfinity) = c
+pad(c, n::Integer) = _pad(c, length(c), n)
+function Fun(sp::Space, c::AbstractVector)
+    coeff = pad(c, size(sp,2))
+    Fun{typeof(sp),eltype(coeff),typeof(coeff)}(sp, coeff)
+end
 Fun() = Fun(identity)
-Fun(d::Domain) = Fun(identity,d)
-Fun(d::Space) = Fun(identity,d)
+Fun(d) = Fun(identity, d)
+
+Fun(v::AbstractQuasiVector) = Fun(arguments(v)...)
+function Fun(f::Function, S::Space)
+    if !applicable(f, checkpoints(S)[1])
+        Fun(x -> f(x...), S)
+    else
+        Fun(S, S \ f.(axes(S,1)))
+    end
+end
 
+Fun(f::Function, d) = Fun(f, Space(d))
 
-function Fun(sp::Space,v::AbstractVector{Any})
+
+function Fun(sp::Space, v::AbstractVector{Any})
     if isempty(v)
-        Fun(sp,Float64[])
+        Fun(sp, Float64[])
     elseif all(x->isa(x,Number),v)
-        Fun(sp,Vector{mapreduce(typeof,promote_type,v)}(v))
+        Fun(sp, Vector{mapreduce(typeof,promote_type,v)}(v))
     else
         error("Cannot construct Fun with coefficients $v and space $sp")
     end
 end
 
 
-hasnumargs(f::Fun,k) = k == 1 || domaindimension(f) == k  # all funs take a single argument as a Vec
+hasnumargs(f::Fun, k) = k == 1 || domaindimension(f) == k  # all funs take a single argument as a Vec
 
 ##Coefficient routines
 #TODO: domainscompatible?
 
-function coefficients(f::Fun,msp::Space)
-    #zero can always be converted
-    fc = f.coefficients
-    if ncoefficients(f) == 0 || (ncoefficients(f) == 1 && fc[1] == 0)
-        fc
-    else
-        coefficients(fc, space(f), msp)
-    end
-end
-coefficients(f::Fun,::Type{T}) where {T<:Space} = coefficients(f,T(domain(f)))
-coefficients(f::Fun) = f.coefficients
-coefficients(c::Number,sp::Space) = Fun(c,sp).coefficients
-
-function coefficient(f::Fun,k::Integer)
-    if k > dimension(space(f)) || k < 1
-        throw(BoundsError())
-    elseif k > ncoefficients(f)
-        zero(cfstype(f))
-    else
-        f.coefficients[k]
-    end
-end
 
-function coefficient(f::Fun,kr::AbstractRange)
-    b = maximum(kr)
-    f.coefficients[first(kr):min(b, end)]
-end
+coefficients(f::Fun) = f.coefficients
 
-coefficient(f::Fun,K::Block) = coefficient(f,blockrange(space(f),K.n[1]))
-coefficient(f::Fun,::Colon) = coefficient(f,1:dimension(space(f)))
+vec(f::Fun) = f.space * coefficients(f)
+coefficients(f::Fun, msp::Space) = msp \ vec(f)
+coefficients(c::Number,sp::Space) = coefficients(Fun(c,sp))
+Fun(f::Fun, S::Space) = Fun(S, coefficients(f, S))
 
 ##Convert routines
 
 
-convert(::Type{Fun{S,T,VT}},f::Fun{S}) where {T,S,VT} =
-    Fun{S,T,VT}(f.space, convert(VT,f.coefficients)::VT)
-function convert(::Type{Fun{S,T,VT}},f::Fun) where {T,S,VT}
-    g = Fun(Fun(f.space, convert(VT,f.coefficients)::VT), convert(S,space(f))::S)
-    Fun{S,T,VT}(g.space, g.coefficients)
-end
+convert(::Type{Fun{S,T,VT}}, f::Fun{S}) where {T,S,VT} =
+    Fun(f.space,convert(VT,f.coefficients))
+convert(::Type{Fun{S,T,VT}}, f::Fun) where {T,S,VT} =
+    Fun(Fun(f.space,convert(VT,f.coefficients)),convert(S,space(f)))
 
-function convert(::Type{Fun{S,T}},f::Fun{S}) where {T,S}
-    coeff = convert(AbstractVector{T},f.coefficients)::AbstractVector{T}
-    Fun{S, T, typeof(coeff)}(f.space, coeff)
-end
+convert(::Type{Fun{S,T}}, f::Fun{S}) where {T,S} =
+    Fun(f.space,convert(AbstractVector{T},f.coefficients))
 
 
-convert(::Type{VFun{S,T}},x::Number) where {T,S} =
-    (x==0 ? zeros(T,S(AnyDomain())) : x*ones(T,S(AnyDomain())))::VFun{S,T}
-convert(::Type{Fun{S}},x::Number) where {S} =
-    (x==0 ? zeros(S(AnyDomain())) : x*ones(S(AnyDomain())))::Fun{S}
-convert(::Type{IF},x::Number) where {IF<:Fun} = convert(IF,Fun(x))::IF
+convert(::Type{VFun{S,T}}, x::Number) where {T,S} =
+    x==0 ? zeros(T,S(AnyDomain())) : x*ones(T,S(AnyDomain()))
+convert(::Type{Fun{S}}, x::Number) where {S} =
+    x==0 ? zeros(S(AnyDomain())) : x*ones(S(AnyDomain()))
+convert(::Type{IF}, x::Number) where {IF<:Fun} = convert(IF,Fun(x))
 
 Fun{S,T,VT}(f::Fun) where {S,T,VT} = convert(Fun{S,T,VT}, f)
 Fun{S,T}(f::Fun) where {S,T} = convert(Fun{S,T}, f)
@@ -104,7 +90,7 @@ Base.promote_rule(::Type{Fun{S,T,VT1}},::Type{Fun{S,V,VT2}}) where {T,V,S,VT1,VT
 
 
 # TODO: Never assume!
-Base.promote_op(::typeof(*),::Type{F1},::Type{F2}) where {F1<:Fun,F2<:Fun} =
+Base.promote_op(::typeof(*), ::Type{F1}, ::Type{F2}) where {F1<:Fun,F2<:Fun} =
     promote_type(F1,F2) # assume multiplication is defined between same types
 
 # we know multiplication by numbers preserves types
@@ -143,14 +129,15 @@ coefficients(f::AbstractArray) = f
 
 
 #supports broadcasting and scalar iterator
-const ScalarFun = Fun{S} where S<:Space{D,R} where {D,R<:Number}
-const ArrayFun = Fun{S} where {S<:Space{D,R}} where {D,R<:AbstractArray}
-const MatrixFun = Fun{S} where {S<:Space{D,R}} where {D,R<:AbstractMatrix}
-const VectorFun = Fun{S} where {S<:Space{D,R}} where {D,R<:AbstractVector}
-
-size(f::Fun,k...) = size(space(f),k...)
-length(f::Fun) = length(space(f))
+const ScalarSpace = Space{<:Number}
+const RealSpace = Space{<:Real}
+const ScalarFun = Fun{<:ScalarSpace}
+const ArrayFun = Fun{<:Space{<:AbstractArray}}
+const MatrixFun = Fun{<:Space{<:AbstractMatrix}}
+const VectorFun = Fun{<:Space{<:AbstractVector}}
 
+size(f::Fun, k...) = size(first(axes(space(f),1)), k...)
+length(f::Fun) = length(first(axes(space(f),1)))
 
 getindex(f::ScalarFun, ::CartesianIndex{0}) = f
 getindex(f::ScalarFun, k::Integer) = k == 1 ? f : throw(BoundsError())
@@ -163,8 +150,7 @@ iterate(A::ArrayFun, i=1) = (@_inline_meta; (i % UInt) - 1 < length(A) ? (@inbou
 
 in(x::ScalarFun, y::ScalarFun) = x == y
 
-setspace(v::AbstractVector,s::Space) = Fun(s,v)
-setspace(f::Fun,s::Space) = Fun(s,f.coefficients)
+setspace(f::Fun, s::Space) = Fun(s, coefficients(f))
 
 
 ## domain
@@ -173,61 +159,33 @@ setspace(f::Fun,s::Space) = Fun(s,f.coefficients)
 ## General routines
 
 
-domain(f::Fun) = domain(f.space)
-domain(v::AbstractMatrix{T}) where {T<:Fun} = map(domain,v)
+domain(f::Fun) = axes(f.space, 1)
+domain(v::AbstractMatrix{T}) where {T<:Fun} = map(domain, v)
 domaindimension(f::Fun) = domaindimension(f.space)
 
 setdomain(f::Fun,d::Domain) = Fun(setdomain(space(f),d),f.coefficients)
 
-for op in (:tocanonical,:tocanonicalD,:fromcanonical,:fromcanonicalD,:invfromcanonicalD)
-    @eval $op(f::Fun,x...) = $op(space(f),x...)
-end
-
-for op in (:tocanonical,:tocanonicalD)
-    @eval $op(d::Domain) = $op(d,Fun(identity,d))
-end
-for op in (:fromcanonical,:fromcanonicalD,:invfromcanonicalD)
-    @eval $op(d::Domain) = $op(d,Fun(identity,canonicaldomain(d)))
-end
-
-
 space(f::Fun) = f.space
-spacescompatible(f::Fun,g::Fun) = spacescompatible(space(f),space(g))
-pointscompatible(f::Fun,g::Fun) = pointscompatible(space(f),space(g))
+spacescompatible(f::Fun, g::Fun) = spacescompatible(space(f),space(g))
+pointscompatible(f::Fun, g::Fun) = pointscompatible(space(f),space(g))
 canonicalspace(f::Fun) = canonicalspace(space(f))
 canonicaldomain(f::Fun) = canonicaldomain(space(f))
 
 
 ##Evaluation
 
-function evaluate(f::AbstractVector,S::Space,x...)
-    csp=canonicalspace(S)
-    if spacescompatible(csp,S)
-        error("Override evaluate for " * string(typeof(csp)))
-    else
-        evaluate(coefficients(f,S,csp),csp,x...)
-    end
-end
-
-evaluate(f::Fun,x) = evaluate(f.coefficients,f.space,x)
-evaluate(f::Fun,x,y,z...) = evaluate(f.coefficients,f.space,Vec(x,y,z...))
-
-
-(f::Fun)(x...) = evaluate(f,x...)
+evaluate(f::Fun, x) = vec(f)[x]
+evaluate(f::Fun, x, y, z...) = evaluate(f, Vec(x,y,z...))
+(f::Fun)(x...) = evaluate(f, x...)
 
 dynamic(f::Fun) = f # Fun's are already dynamic in that they compile by type
 
-for (op,dop) in ((:first,:leftendpoint),(:last,:rightendpoint))
-    @eval $op(f::Fun{S,T}) where {S,T} = f($dop(domain(f)))
-end
-
-
 
 ## Extrapolation
 
 
 # Default extrapolation is evaluation. Override this function for extrapolation enabled spaces.
-extrapolate(f::AbstractVector,S::Space,x...) = evaluate(f,S,x...)
+extrapolate(f::AbstractVector, S::Space, x...) = evaluate(f, S, x...)
 
 # Do not override these
 extrapolate(f::Fun,x) = extrapolate(f.coefficients,f.space,x)
@@ -241,38 +199,19 @@ values(f::Fun,dat...) = _values(f.space, f.coefficients, dat...)
 _values(sp, v, dat...) = itransform(sp, v, dat...)
 _values(sp, v::Vector{T}, dat...) where {T} = itransform(sp, v, dat...)::Vector{float(T)}
 points(f::Fun) = points(f.space,ncoefficients(f))
-ncoefficients(f::Fun)::Int = length(f.coefficients)
-blocksize(f::Fun) = (block(space(f),ncoefficients(f)).n[1],)
-
-function stride(f::Fun)
-    # Check only for stride 2 at the moment
-    # as higher stride is very rare anyways
-    M=maximum(abs,f.coefficients)
-    for k=2:2:ncoefficients(f)
-        if abs(f.coefficients[k])>40*M*eps()
-            return 1
-        end
-    end
-
-    2
-end
-
-
+ncoefficients(f::Fun) = nzeros(f.coefficients)
 
 ## Manipulate length
 
-pad!(f::Fun,n::Integer) = (pad!(f.coefficients,n);f)
-pad(f::Fun,n::Integer) = Fun(f.space,pad(f.coefficients,n))
-
 
-function chop!(sp::UnivariateSpace,cfs,tol::Real)
-    n=standardchoplength(cfs,tol)
+function chop!(sp::ScalarSpace, cfs, tol::Real)
+    n = standardchoplength(cfs, tol)
     resize!(cfs,n)
     cfs
 end
 
-chop!(sp::Space,cfs,tol::Real) = chop!(cfs,maximum(abs,cfs)*tol)
-chop!(sp::Space,cfs) = chop!(sp,cfs,10eps())
+chop!(sp::Space, cfs, tol::Real) = chop!(cfs,maximum(abs,cfs)*tol)
+chop!(sp::Space, cfs) = chop!(sp,cfs,10eps())
 
 function chop!(f::Fun,tol...)
     chop!(space(f),f.coefficients,tol...)
@@ -286,32 +225,31 @@ copy(f::Fun) = Fun(space(f),copy(f.coefficients))
 
 ## Addition and multiplication
 
+
+
 for op in (:+,:-)
     @eval begin
-        function $op(f::Fun,g::Fun)
-            if spacescompatible(f,g)
-                n = max(ncoefficients(f),ncoefficients(g))
-                f2 = pad(f,n); g2 = pad(g,n)
-
-                Fun(isambiguous(domain(f)) ? g.space : f.space,($op)(f2.coefficients,g2.coefficients))
+        function $op(f::Fun, g::Fun)
+            if space(f) == space(g)
+                Fun(space(f), ($op)(coefficients(f), coefficients(g)))
             else
-                m=union(f.space,g.space)
-                if isa(m,NoSpace)
+                m = broadcastspace($op, f.space, g.space)
+                if m isa NoSpace
                     error("Cannot "*string($op)*" because no space is the union of "*string(typeof(f.space))*" and "*string(typeof(g.space)))
                 end
-                $op(Fun(f,m),Fun(g,m)) # convert to same space
+                $op(Fun(f,m), Fun(g,m)) # convert to same space
             end
         end
         $op(f::Fun{S,T},c::T) where {S,T<:Number} = c==0 ? f : $op(f,Fun(c))
-        $op(f::Fun,c::Number) = $op(f,Fun(c))
+        $op(f::Fun,c::Number) = $op(f,Fun(c,space(f)))
         $op(f::Fun,c::UniformScaling) = $op(f,c.λ)
         $op(c::UniformScaling,f::Fun) = $op(c.λ,f)
     end
 end
 
 
 # equivalent to Y+=a*X
-axpy!(a,X::Fun,Y::Fun)=axpy!(a,coefficients(X,space(Y)),Y)
+axpy!(a, X::Fun, Y::Fun)=axpy!(a,coefficients(X,space(Y)),Y)
 function axpy!(a,xcfs::AbstractVector,Y::Fun)
     if a!=0
         n=ncoefficients(Y); m=length(xcfs)
@@ -335,8 +273,9 @@ function axpy!(a,xcfs::AbstractVector,Y::Fun)
 end
 
 
+
 +(a::Fun) = copy(a)
--(f::Fun) = Fun(f.space,-f.coefficients)
+-(f::Fun) = Fun(f.space, -f.coefficients)
 -(c::Number,f::Fun) = -(f-c)
 
 
@@ -373,7 +312,6 @@ inv(f::Fun) = 1/f
 
 # Integrals over two Funs, which are fast with the orthogonal weight.
 
-export bilinearform, linebilinearform, innerproduct, lineinnerproduct
 
 # Having fallbacks allow for the fast implementations.
 
@@ -406,7 +344,7 @@ for (OP,SUM) in ((:(norm),:(sum)),(:linenorm,:linesum))
     @eval begin
         $OP(f::Fun) = $OP(f,2)
 
-        function $OP(f::Fun{S},p::Real) where S<:Space{D,R} where {D,R<:Number}
+        function $OP(f::ScalarFun, p::Real)
             if p < 1
                 return error("p should be 1 ≤ p ≤ ∞")
             elseif 1 ≤ p < Inf
@@ -416,7 +354,7 @@ for (OP,SUM) in ((:(norm),:(sum)),(:linenorm,:linesum))
             end
         end
 
-        function $OP(f::Fun{S},p::Int) where S<:Space{D,R} where {D,R<:Number}
+        function $OP(f::ScalarFun,p::Int)
             if 1 ≤ p < Inf
                 return iseven(p) ? abs($SUM(abs2(f)^(p÷2)))^(1/p) : abs($SUM(abs2(f)^(p/2)))^(1/p)
             else
@@ -459,6 +397,11 @@ function differentiate(f::Fun,k::Integer)
 end
 
 # use conj(transpose(f)) for ArraySpace
+function differentiate(f)
+    v = vec(f)
+    D = Derivative(axes(v,1))
+    Fun(D*v)
+end
 adjoint(f::Fun) = differentiate(f)
 
 
@@ -601,5 +544,3 @@ function copyto!(dest::Fun, bc::Broadcasted{FunStyle})
     dest.coefficients[:] = cfs
     dest
 end
-
-include("constructors.jl")