juliaapproximation / approxfunbase.jl Goto Github PK
View Code? Open in Web Editor NEWCore functionality of ApproxFun
License: MIT License
Core functionality of ApproxFun
License: MIT License
Are there instances where one needs to modify the fields of a Fun
, as opposed to the coefficients?
using ApproxFun
Dc = Chebyshev(Interval(0,1.))
Df = Fourier(Interval(0,1.))
rx = Fun(t -> cospi(2*t), Df)
ry = Fun(t -> sinpi(2*t), Df)
twoπx = Fun(x -> 2π*x, Dc)
ell = sqrt(differentiate(rx)^2 + differentiate(ry)^2)
L = cumsum(ell)
println(L(1) == twoπx(1)) # expected value: true
println((L - twoπx)(1) ≈ 0) # therefore, expected value: true
println((L - twoπx)(1) == L(1)) # expected value: false
Doing this on a chalkboard we would find L := 2πx and so (L - twoπx) ≡ 0. However, the actual results from this example are: true, false, true.
I want to apply some operators like derivatives to Fun
s and would like to know how many coefficients and points I'll get out. I found that
ApproxFunBase.jl/src/Operators/Operator.jl
Lines 401 to 404 in 25ee448
seems to be where the number of output coefficients is decided. I tried to look into the colstop
function but I can't seem to figure out what it does. The one that gets called in my case is:
ApproxFunBase.jl/src/PDE/KroneckerOperator.jl
Lines 77 to 83 in 25ee448
Could you explain what the two functions shown above do and why? Specifically, I want to understand how to predict what number of output coefficients derivatives on a Chebyshev^2 basis will produce and what number of coefficients I will get when converting from Chebyshev to Ultraspherical basis.
Here's a MWE
using FastTransforms: paduapoints
using ApproxFunOrthogonalPolynomials
np = 1000
ppoints = paduapoints(ceil(Int, (-3 + sqrt(1 + 8 * np)) / 2))
f(v) = cos(v[1]) * sin(v[2])
vals = f.(eachrow(ppoints))
CC = Chebyshev(-1..1) ⊗ Chebyshev(-1..1)
UC = Ultraspherical(1, -1..1) ⊗ Chebyshev(-1..1)
CU = Chebyshev(-1..1) ⊗ Ultraspherical(1, -1..1)
UU = Ultraspherical(1, -1..1) ⊗ Ultraspherical(1, -1..1)
CCtoUU = Conversion(CC, UU)
D1 = Derivative(CC, [1,0])
UCtoUU = Conversion(UC, UU)
coeffs = transform(CC, vals)
pf = Fun(CC, coeffs)
ncoefficients(pf) # 1035
ncoefficients(D1 * pf) # 990: why?
ncoefficients(UCtoUU * (D1 * pf)) # also 990
ncoefficients(CCtoUU * pf) # 1035
This function, at https://github.com/JuliaApproximation/ApproxFunBase.jl/blob/master/src/Operators/general/algebra.jl#L574, doesn't quite work with dual numbers, particularly because dual(0,1) == 0
and dual(1,1) == 1
, which do not behave like other field extensions, e.g. complex.
julia> dual(0,1)*Evaluation(S, 1)
ZeroOperator : Chebyshev() → ConstantSpace(Point(1))
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋯
julia> dual(eps(),1)*Evaluation(S, 1)
ConstantTimesOperator : Chebyshev() → ConstantSpace(Point(1))
2.220446049250313e-16 + 1.0ɛ … 2.220446049250313e-16 + 1.0ɛ ⋯
julia> dual(1,1)*Evaluation(S, 1)
ConcreteEvaluation : Chebyshev() → ConstantSpace(Point(1))
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 ⋯
julia> dual(1+eps(),1)*Evaluation(S, 1)
ConstantTimesOperator : Chebyshev() → ConstantSpace(Point(1))
1.0000000000000002 + 1.0ɛ … 1.0000000000000002 + 1.0ɛ ⋯
Currently,
julia> Derivative(Chebyshev()^2, [1,0])
DerivativeWrapper : Chebyshev() ⊗ Chebyshev() → Ultraspherical(1) ⊗ Chebyshev()
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋯
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋱
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ⋱
⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱
julia> Integral(Chebyshev()^2, [1,0])
ERROR: Implement Integral(Chebyshev() ⊗ Chebyshev(),[1, 0])
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] DefaultIntegral(sp::TensorSpace{Tuple{Chebyshev{ChebyshevInterval{Float64}, Float64}, Chebyshev{ChebyshevInterval{Float64}, Float64}}, DomainSets.FixedIntervalProduct{2, Float64, ChebyshevInterval{Float64}}, Float64}, k::Vector{Int64})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/VD9Rj/src/Operators/banded/CalculusOperator.jl:41
[3] Integral(::TensorSpace{Tuple{Chebyshev{ChebyshevInterval{Float64}, Float64}, Chebyshev{ChebyshevInterval{Float64}, Float64}}, DomainSets.FixedIntervalProduct{2, Float64, ChebyshevInterval{Float64}}, Float64}, ::Vararg{Any})
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/VD9Rj/src/Operators/banded/CalculusOperator.jl:53
[4] top-level scope
@ REPL[94]:1
Integral(Chebyshev()) ⊗ I
KroneckerOperator : Chebyshev() ⊗ ApproxFunBase.UnsetSpace() → Chebyshev() ⊗ ApproxFunBase.UnsetSpace()
It would be nice to have a Syntax for Integral(Chebyshev()^2, [1,0])
for integration over the first variable as well as for DefiniteIntegral
.
Furthermore, currently the following throws an error:
julia> Op = Integral(Chebyshev()) ⊗ I ⊗ I
KroneckerOperator : Chebyshev() ⊗ ApproxFunBase.UnsetSpace() ⊗ ApproxFunBase.UnsetSpace() → Chebyshev() ⊗ ApproxFunBase.UnsetSpace() ⊗ ApproxFunBase.UnsetSpace()
julia> f3 = Fun(Chebyshev()^3, rand(20));
julia> Op * f3
ERROR: Implement Conversion from Chebyshev{ChebyshevInterval{Float64}, Float64} to TensorSpace{Tuple{Chebyshev{ChebyshevInterval{Float64}, Float64}, ApproxFunBase.UnsetSpace}, DomainSets.VcatDomain{2, Float64, (1, 1), Tuple{ChebyshevInterval{Float64}, DomainSets.FullSpace{Float64}}}, Union{}}
Stacktrace:
as well as the same for the Derivative.
Maybe I am misusing the package for something it is not intended to be, but I want to use it for implementing closed form integration, marginalization, and derivatives of SOS polynomials and therefore would need this operations for higher dimensions.
@jishnub Do you have an advice for me where I have to start for implementing this?
At this moment
julia> S = ApproxFunBase.SubSpace(Fourier(), 1:10);
julia> Derivative(S)
DerivativeWrapper : Fourier(【0.0,6.283185307179586❫)|1:10 → SinSpace(【0.0,6.283185307179586❫)⊕CosSpace(【0.0,6.283185307179586❫)
0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -2.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 -3.0 0.0 0.0 0.0
0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -4.0 0.0
0.0 0.0 0.0 0.0 0.0 3.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 0.0 0.0
while the first dimension could be finite.
Where would be the best place to interrupt the ApproxFunBase
flow to create an Operator
that is finite in both dimensions?
ApproxFunBase.jl/src/LinearAlgebra/helper.jl
Lines 57 to 59 in 451f1fe
The comment suggests that ApproxFunBase.real
differs from Base.real
, but it's actually identical as ApproxFunBase
imports Base.real
. This leads to stack overflows if used incorrectly.
Eg.
julia> ApproxFunBase.real === Base.real
true
julia> real(1,2)
ERROR: StackOverflowError:
Stacktrace:
[1] real(::Int64, ::Int64) (repeats 79984 times)
@ ApproxFunBase ~/.julia/packages/ApproxFunBase/gHG1c/src/LinearAlgebra/helper.jl:59
This should ideally be a method error.
julia> f2_1 = Fun(Chebyshev()^2, rand(20));
julia> f2_2 = Fun(Chebyshev()^2, rand(30));
julia> (f2_1*f2_2)(0.2, 0.2) ≈ f2_1(0.2,0.2)*f2_2(0.2,0.2)
true
julia> (f2_1*f2_2)(0.1, 0.8) ≈ f2_1(0.1,0.8)*f2_2(0.1,0.8)
false # should be true
multivariate function multiplication is not working properly
ZeroOperator
maps whatever Fun
to zero on a specific Domain
. Whatever space on that Domain
doesn't matter.
For example, there is an Operator
called C
from space S1
to space S2
. There is a ZeroOperator
called D
from space S1
to space S3
. S2
and S3
are both spaces on [-1,1]
.
Executing C+D
will trigger a conversion between S2
and S3
, which can't succeed in general for random S2
and S3
. But C+D
should be C
AFAIK.
IMO the field rangespace
of ZeroOperator
should be replaced by domain
.
The definition here:
Line 506 in 6637707
changes the definition of iszero
for all numbers (type piracy) and causes a huge number of invalidations
inserting iszero(x::Number) in ApproxFunBase at /home/kc/.julia/packages/ApproxFunBase/r0W9A/src/Fun.jl:506 invalidated:
backedges: 1: superseding iszero(x) in Base at number.jl:40 with MethodInstance for iszero(::UInt32) (1 children)
2: superseding iszero(x) in Base at number.jl:40 with MethodInstance for iszero(::UInt64) (2 children)
3: superseding iszero(x) in Base at number.jl:40 with MethodInstance for iszero(::Int64) (6236 children) <---------------
You can see the effect of this by loading ApproxFunBase and notice that the whole julia session becomes sluggish while it is recompiling everything.
Currently
julia> D1 = Derivative(Chebyshev()) * Derivative(Chebyshev())
ConcreteDerivative : Chebyshev() → Ultraspherical(2)
⋅ ⋅ 4.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 6.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 8.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 12.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 14.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 16.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 18.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
julia> D2 = Derivative(Ultraspherical(1)) * Derivative(Chebyshev())
ConcreteDerivative : Chebyshev() → Ultraspherical(2)
⋅ ⋅ 4.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 6.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 8.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 10.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 12.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 14.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 16.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 18.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋱
julia> D1 == D2 == (Derivative()^2 : Chebyshev())
true
Shouldn't the first one check for compatible domain and range spaces?
I am working with vector valued functions, R -> R^n, and would like to use the Evaluation operator.
If we have f = Fun(x -> [1, 2])
we get the following problems:
Evaluation(0.5)*f
gives an error as the implementation tries to call zero(Vector{Float64})
, which is undefined;Evaluation(-1)*f
and Evaluation(1)*f
still give the same error as before;rangespace(Evaluation(space(f), 0.)
returns ConstantSpace(Point(0.))
instead of ArraySpace(ConstantSpace(Point(0.)), 2)
; andFiniteOperator([1. 2.; 3., 4.]) * Evaluation(space(f), 0.)
errors.I get an stackoverflow error when I try to apply an adjoint of a blockbanded operator on a function. For instance
Op=Derivative(Chebyshev())⊗Derivative(Chebyshev())
f=Fun(sin)⊗Fun(sin)
Op*f # works fine
Op'*f # Stackoverflow
I believe this is due to blockbandwidths(P::AdjointOperator)
not being implemented for adjoint operators (while bandwidths(P::AdjointOperator)
is implemented in OperatorLayout.jl) and thus the operator is not considered blockbanded. Indeed
println(ApproxFunBase.blockbandwidths(Op)) # (-2, 2)
println(ApproxFunBase.blockbandwidths(Op')) # (ℵ₀, ℵ₀)
The following implementation did the job for me:
ApproxFunBase.blockbandwidths(A::AdjointOperator) = reverse(ApproxFunBase.blockbandwidths(A.op))
println(ApproxFunBase.blockbandwidths(Op')) # (2,-2)
and Op'*f
works.
Debug the error code, I found that
almostbanded.jl line81 --->ds=domainspace(io)
almostbanded.jl line106 ---->bl = blocklengths(ds[k])
where, k is the integer, however, In the below case, ds=domainspace(A),is like Laurent(0.1🕒+1.0 + 0.0im)⨄Laurent(0.1🕒+-0.4999999999999998 + 0.8660254037844387im)⨄Laurent(0.1🕒+-0.5000000000000004 - 0.8660254037844385im)
which means ds[1] is no meaning!!!
This headache me for many days, when I try to complete the example code in SingularIntergralEquation Package.
I tested my code here:
using ApproxFun, SingularIntegralEquations,Plots,ApproxFunBase
g1(x,y) = 0.5
N=3
r=0.1
cr = exp.(im*2π*(0:N-1)/N)
crl,crr = (1-2im*r)cr,(1+2im*r)cr
@show dom=Circle(cr[1],r) ∪ Circle(cr[2],r)∪ Circle(cr[3],r)
#@show dom = ∪(Segment.(crl[1:end],crr[1:end])...)
@show sp = Space(dom)
cwsp = CauchyWeight(sp⊗sp,0)
⨍ = DefiniteLineIntegral(dom)
@time G = GreensFun(g1,cwsp;method=:Cholesky)
A=⨍[G]
ApproxFun.qr(A) #error occurs
Here A is union InterlaceOperator such like that: Laurent(0.1🕒+1.0 + 0.0im)⨄Laurent(0.1🕒+-0.4999999999999998 + 0.8660254037844387im)⨄Laurent(0.1🕒+-0.5000000000000004 - 0.8660254037844385im) → Laurent(0.1🕒+1.0 + 0.0im)⨄Laurent(0.1🕒+-0.4999999999999998 + 0.8660254037844387im)⨄Laurent(0.1🕒+-0.5000000000000004 - 0.8660254037844385im)
I see that the line to return such a matrix is commented out below. Why?
ApproxFunBase.jl/src/Operators/Operator.jl
Line 776 in 25ee448
I'd like to preallocate large operator matrices and currently the RaggedMatrix
type seems to be the default. The issue is that RaggedMatrix
uses too much memory. Currently, I'm just manually using BandedBlockBandedMatrix
, but that's not as nice as just being able to use AbstractMatrix(view(SomeOperator, FiniteRange, 1:N)))
, since I don't want to have to care about the exact memory layout myself.
This doesn't feel right...
julia> S = Chebyshev(1..2)⊗Chebyshev(3..4)
Chebyshev(1 .. 2) ⊗ Chebyshev(3 .. 4)
julia> [p ∈ domain(factor(S, 1)) for p in points(S, 10, 10)[1]]
10×10 Matrix{Bool}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
julia> [p ∈ domain(factor(S, 2)) for p in points(S, 10, 10)[2]]
10×10 Matrix{Bool}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
The following works well in Julia 1.9 but fails in Julia 1.10
In Julia 1.9:
julia> Af = Fun(t -> [0 1; -10*cos(t)-1 -24-19*sin(t)],Fourier(0..2pi));
julia> iszero(Af)
false
The old definition of iszero
is
iszero(f::Fun) = all(iszero,f.coefficients)
In Julia 1.10:
julia> Af = Fun(t -> [0 1; -10*cos(t)-1 -24-19*sin(t)],Fourier(0..2pi));
julia> iszero(Af)
ERROR: Override for ApproxFunBase.ArraySpace(ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}[Fourier(【0.0,6.283185307179586❫) Fourier(【0.0,6.283185307179586❫); Fourier(【0.0,6.283185307179586❫) Fourier(【0.0,6.283185307179586❫)])
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] checkcanonicalspace(sp::ApproxFunBase.ArraySpace{ApproxFunBase.SumSpace{…}, 2, PeriodicSegment{…}, Float64, Matrix{…}})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:478
[3] ApproxFunBase.ICanonicalTransformPlan(space::ApproxFunBase.ArraySpace{…}, v::Vector{…}, ip::Val{…})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:493
[4] plan_itransform(sp::ApproxFunBase.ArraySpace{…}, v::Vector{…})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:497
[5] itransform(S::ApproxFunBase.ArraySpace{…}, cfs::Vector{…})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:564
[6] _values
@ C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Fun.jl:428 [inlined]
[7] values
@ C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Fun.jl:426 [inlined]
[8] iszero(f::Fun{ApproxFunBase.ArraySpace{…}, Float64, Vector{…}})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Fun.jl:776
[9] top-level scope
@ REPL[90]:1
Some type information was truncated. Use `show(err)` to see complete types.
The new definition of iszero
is
iszero(f::Fun) = all(iszero, coefficients(f)) || all(iszero, values(f))
so, the failure is caused by
julia> values(Af)
ERROR: Override for ApproxFunBase.ArraySpace(ApproxFunBase.SumSpace{Tuple{CosSpace{PeriodicSegment{Float64}, Float64}, SinSpace{PeriodicSegment{Float64}, Float64}}, PeriodicSegment{Float64}, Float64}[Fourier(【0.0,6.283185307179586❫) Fourier(【0.0,6.283185307179586❫); Fourier(【0.0,6.283185307179586❫) Fourier(【0.0,6.283185307179586❫)])
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] checkcanonicalspace(sp::ApproxFunBase.ArraySpace{ApproxFunBase.SumSpace{…}, 2, PeriodicSegment{…}, Float64, Matrix{…}})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:478
[3] ApproxFunBase.ICanonicalTransformPlan(space::ApproxFunBase.ArraySpace{…}, v::Vector{…}, ip::Val{…})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:493
[4] plan_itransform(sp::ApproxFunBase.ArraySpace{…}, v::Vector{…})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:497
[5] itransform(S::ApproxFunBase.ArraySpace{…}, cfs::Vector{…})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Space.jl:564
[6] _values
@ C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Fun.jl:428 [inlined]
[7] values(::Fun{ApproxFunBase.ArraySpace{…}, Float64, Vector{…}})
@ ApproxFunBase C:\Users\Andreas\.julia\packages\ApproxFunBase\sJBXX\src\Fun.jl:426
[8] top-level scope
@ REPL[91]:1
Some type information was truncated. Use `show(err)` to see complete types.
I wonder why the second test involving all(iszero, values(f))
is necessary!
I would appreciate very much if the fix for this very basic function can be done in a short time. Thanks in advance.
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
If you have a basis of Chebyshev Polynomials, for example, that looks like this:
[T00 T01 T02;
T10 T11 0 ;
T20 0 0 ]
These are stored like so
[T00, T01, T10, T02, T11, T20]
My question is: why? Would it not be more natural to simply keep using the matrix? Or perhaps reshape the matrix to a vector, so columns are concatenated together?
In my case, using the matrix directly means I don't have a differentiation operator with 1e12
entries, but rather one with 1e6
entries, which is the difference between feasible and impossible. (The reduction occurs because in the matrix way of doing things, I can use the 1D operators on each column/row, which amounts to a simple matrix matrix multiplication).
Of course, this comes at the cost of doubling the memory, but it seems more than worth it to double the memory at this stage if you can save many orders of magnitude in memory/time when building operators and applying them.
The following code returns successfully for solving Airy's equation (the extra definitions are work arounds for bugs):
julia> Base.broadcast(::typeof(*), a::IntervalArithmetic.Interval, b::AbstractVector) =
map(x -> a*x, b)
julia> Base.Broadcast.broadcasted(::typeof(*), a::IntervalArithmetic.Interval, b::AbstractVector) =
map(x -> a*x, b)
julia> Base.length(d::UnitRange{<:IntervalArithmetic.Interval}) = Integer(maximum(d.stop.hi)-minimum(d.start.lo)+1)
julia> u = [Dirichlet(); D^2 - x] \ [[1, 0], 0]
Fun(Chebyshev([-1, -1]..[1, 1]),IntervalArithmetic.Interval{Float64}[[0.528647, 0.528648], [-0.522247, -0.522246], [-0.023162, -0.0231619], [0.0223854, 0.0223855], [-0.00557882, -0.00557881], [-0.000122104, -0.000122103], [9.37066e-05, 9.37067e-05], [-1.68122e-05, -1.68121e-05], [-2.4373e-07, -2.43729e-07], [1.63064e-07, 1.63065e-07] … [1.54639e-10, 1.5464e-10], [-1.89644e-11, -1.89643e-11], [-1.63813e-13, -1.63812e-13], [9.21385e-14, 9.21386e-14], [-9.91957e-15, -9.91956e-15], [-7.12466e-17, -7.12465e-17], [3.76658e-17, 3.76659e-17], [-3.63796e-18, -3.63795e-18], [-2.23433e-20, -2.23432e-20], [-1.12162e-20, -1.12161e-20]])
It happens to have worked:
julia> parse(BigFloat,"0.546136459064141770613483785351029091974067802242010192661484687977676021081218637004123196030567166414008915304513669838826570884623515686") in u(0) # true value from mathematica
true
However, this is mostly luck: the convergence criteria doesn't actually control the decay in the tail. To do this properly (and thereby have apriori guaranteed success) we would need to modify the convergence check at
ApproxFunBase.jl/src/Caching/matrix.jl
Line 122 in 7dd9ef0
to properly bound the tail.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.