Add standard chop (#94)

* Use standard chop to see if tail is zero

* add standardchop

* weaker zero test

* tests pass

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.6"
+version = "0.6.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -28,7 +28,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ArrayLayouts = "0.8"
 BandedMatrices = "0.17"
 BlockArrays = "0.16.9"
-BlockBandedMatrices = "0.11"
+BlockBandedMatrices = "0.11.6"
 ContinuumArrays = "0.10"
 DomainSets = "0.5.6"
 FFTW = "1.1"
@@ -37,20 +37,21 @@ FastTransforms = "0.13"
 FillArrays = "0.13"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.3"
-InfiniteLinearAlgebra = "0.6.5"
+InfiniteLinearAlgebra = "0.6.7"
 IntervalSets = "0.5, 0.6"
 LazyArrays = "0.22"
-LazyBandedMatrices = "0.7"
+LazyBandedMatrices = "0.7.14"
 QuasiArrays = "0.9"
 SpecialFunctions = "1.0, 2"
 julia = "1.7"
 
 [extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays"]
+test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays", "Random"]

src/ClassicalOrthogonalPolynomials.jl

@@ -32,7 +32,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     AbstractQuasiFill, _dot, _equals, QuasiArrayLayout, PolynomialLayout
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
-import InfiniteLinearAlgebra: chop!, chop
+import InfiniteLinearAlgebra: chop!, chop, choplength, compatible_resize!
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     inbounds_getindex, grid, plotgrid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
     AffineQuasiVector, AffineMap, WeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
@@ -66,6 +66,7 @@ cardinality(::EuclideanDomain) = ℵ₁
 cardinality(::Union{DomainSets.RealNumbers,DomainSets.ComplexNumbers}) = ℵ₁
 cardinality(::Union{DomainSets.Integers,DomainSets.Rationals,DomainSets.NaturalNumbers}) = ℵ₀
 
+include("standardchop.jl")
 include("adaptivetransform.jl")
 
 const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}

src/adaptivetransform.jl

@@ -11,6 +11,10 @@ padchop(cfs, tol, ax...) = pad(chop(cfs, tol), ax...)
 padchop!(cfs::PseudoBlockVector, tol, ax...) = padchop!(cfs.blocks, tol, ax...)
 
 
+padresize!(cfs, m, ax...) = pad(compatible_resize!(cfs, m), ax...)
+padresize!(cfs::PseudoBlockVector, m, ax...) = padresize!(cfs.blocks, m, ax...)
+
+
 increasingtruncations(::OneToInf) = oneto.(2 .^ (4:∞))
 increasingtruncations(::BlockedUnitRange) = broadcast(n -> Block.(oneto(n)), (2 .^ (4:∞)))
 function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiVector{V}) where {U,V}
@@ -27,16 +31,18 @@ function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiVector
         An = A[:,jr]
         cfs = An \ f
         maxabsc = maximum(abs, cfs)
-        if maxabsc == 0 && maxabsfr == 0
+        if maxabsc ≤ tol && maxabsfr ≤ tol # probably zero
             return pad(similar(cfs,0), ax)
         end
 
         m = length(cfs)
-        if maximum(abs,@views(cfs[max(1,m-2):m])) < 10tol*maxabsc
-            n = length(cfs)
-            c = padchop!(cfs, tol, ax)
+        m̃ = standardchoplength(cfs, tol)
+
+        
+        if m̃ < m-1 # coefficient tail is "zero" based on standard chop
+            c = padresize!(cfs, m̃, ax)
             un = A * c
-            if all(norm.(un[r] - fr, 1) .< tol * n * maxabsfr*1000)
+            if all(norm.(un[r] - fr, 1) .<  m*tol*1000*max(maxabsfr, 1))
                 return c
             end
         end

src/standardchop.jl

@@ -0,0 +1,97 @@
+#  Contains code that is based in part on Chebfun v5's chebfun/standardChop.m,
+# which is distributed with the following license:
+
+# Copyright (c) 2017, The Chancellor, Masters and Scholars of the University
+# of Oxford, and the Chebfun Developers. All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#     * Redistributions of source code must retain the above copyright
+#       notice, this list of conditions and the following disclaimer.
+#     * Redistributions in binary form must reproduce the above copyright
+#       notice, this list of conditions and the following disclaimer in the
+#       documentation and/or other materials provided with the distribution.
+#     * Neither the name of the University of Oxford nor the names of its
+#       contributors may be used to endorse or promote products derived from
+#       this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+# Jared Aurentz and Nick Trefethen, July 2015.
+#
+# Copyright 2017 by The University of Oxford and The Chebfun Developers.
+# See http://www.chebfun.org/ for Chebfun information.
+
+function standardchoplength(coeffs, tol)
+    # Check magnitude of TOL:
+    if tol ≥ 1
+        throw(ArgumentError("tolerance must be less than 1"))
+    end
+
+    # Make sure COEFFS has length at least 17:
+    n = length(coeffs)
+
+    if  n < 17
+        # resort to naive chop
+        mx = maximum(abs,coeffs)
+        return choplength(coeffs, tol*mx)
+    end
+
+    # Step 1: Convert COEFFS to a new monotonically nonincreasing
+    #         vector ENVELOPE normalized to begin with the value 1.
+
+    b = convert(Vector, abs.(coeffs))
+    for j = n-1:-1:1
+        b[j] = max(b[j], b[j+1])
+    end
+    iszero(b[1]) && return 1
+
+    b .= b ./ b[1]
+
+
+    plateauPoint = 0
+
+    for j = 2:n
+        j2 = round(Int,1.25*j + 5);
+        if  j2 > n
+            # there is no plateau: exit
+            return n
+        end
+        e1 = b[j]
+        e2 = b[j2]
+        r = 3*(1 - log(e1)/log(tol))
+        plateau = (e1 == 0) || (e2/e1 > r)
+        if plateau
+            # a plateau has been found: go to Step 3
+            plateauPoint = j - 1
+            break
+        end
+    end
+
+    # Step 3: fix CUTOFF at a point where ENVELOPE, plus a linear function
+    # included to bias the result towards the left end, is minimal.
+    #
+
+    if  plateauPoint ≠ 0 && b[plateauPoint] == 0
+        return plateauPoint
+    end
+
+    j3 = sum(b .≥ tol^(7/6))
+    if j3 < j2
+        j2 = j3 + 1
+        b[j2] = tol^(7/6)
+    end
+    cc = log10.(b[1:j2])
+    cc .+= range(0, stop=(-1/3)*log10(tol), length=j2)
+    d = argmin(cc)
+    return max(d - 1, 1)
+end

test/runtests.jl

@@ -1,6 +1,6 @@
 using Base, ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, FillArrays,
         LazyArrays, BandedMatrices, LinearAlgebra, FastTransforms, IntervalSets,
-        InfiniteLinearAlgebra, Test
+        InfiniteLinearAlgebra, Random, Test
 using ForwardDiff, SemiseparableMatrices, SpecialFunctions, LazyBandedMatrices
 import ContinuumArrays: BasisLayout, MappedBasisLayout
 import ClassicalOrthogonalPolynomials: jacobimatrix, ∞, ChebyshevInterval, LegendreWeight,
@@ -12,6 +12,8 @@ import ClassicalOrthogonalPolynomials: oneto
 import InfiniteLinearAlgebra: KronTrav, Block
 import FastTransforms: clenshaw!
 
+Random.seed!(0)
+
 @testset "singularities" begin
     x = Inclusion(ChebyshevInterval())
     @test singularities(x) == singularities(exp.(x)) == singularities(x.^2) == 

test/test_interlace.jl

@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, FillArrays, ContinuumArrays, StaticArrays, Test
-import ClassicalOrthogonalPolynomials: PiecewiseInterlace, SetindexInterlace, plotgrid
+import ClassicalOrthogonalPolynomials: PiecewiseInterlace, SetindexInterlace, plotgrid, BroadcastQuasiVector
 
 @testset "Interlace" begin
     @testset "Piecewise" begin
@@ -152,5 +152,17 @@ import ClassicalOrthogonalPolynomials: PiecewiseInterlace, SetindexInterlace, pl
             U = V / V \ broadcast(x -> cos.((1:10) .* x), x)
             @test U[0.1] ≈ cos.((1:10) .* 0.1)
         end
+
+        @testset "zero" begin
+            Bˣ = [-1/4 1/4 1/2 1/2; 1/4 -1/4 1/2 1/2; 0 0 -1/4 1/4; 0 0 1/4 -1/4]
+            Bʸ = [1/4 1/4 -1/2 1/2; 1/4 1/4 1/2 -1/2; 0 0 1/4 1/4; 0 0 1/4 1/4]
+            X = z -> Bˣ/z + Bˣ'*z
+            Y = z -> Bʸ/z + Bʸ'*z
+            F = Fourier{ComplexF64}()
+            S = SetindexInterlace(SMatrix{4,4,ComplexF64},fill(F,4^2)...)
+            θ = axes(S,1)
+            XY = S \ BroadcastQuasiVector{eltype(S)}(θ -> (I-X(exp(im*θ))^2)*(I-Y(exp(im*θ))^2), θ)
+            @test iszero(norm(XY))
+        end
     end
 end

test/test_legendre.jl

@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, LazyArrays, QuasiArrays, BandedMatrices, ContinuumArrays, ForwardDiff, Test
-import ClassicalOrthogonalPolynomials: recurrencecoefficients, jacobimatrix, Clenshaw, weighted
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, jacobimatrix, Clenshaw, weighted, oneto
 import QuasiArrays: MulQuasiArray
 
 @testset "Legendre" begin

test/test_stieltjes.jl

@@ -310,8 +310,8 @@ end
 
         @testset "inversion" begin
             H̃ = BlockHcat(Eye((axes(H,1),))[:,Block(1)], H)
-            @test BlockArrays.blockcolsupport(H̃,1) == Block.(1:1)
-            @test last(BlockArrays.blockcolsupport(H̃,2)) ≤ Block(30)
+            @test BlockArrays.blockcolsupport(H̃,Block(1)) == Block.(1:1)
+            @test last(BlockArrays.blockcolsupport(H̃,Block(2))) ≤ Block(30)
 
             UT = U \ T
             D = U \ Derivative(x) * T