|
| 1 | +# Contains code that is based in part on Chebfun v5's chebfun/standardChop.m, |
| 2 | +# which is distributed with the following license: |
| 3 | + |
| 4 | +# Copyright (c) 2017, The Chancellor, Masters and Scholars of the University |
| 5 | +# of Oxford, and the Chebfun Developers. All rights reserved. |
| 6 | +# |
| 7 | +# Redistribution and use in source and binary forms, with or without |
| 8 | +# modification, are permitted provided that the following conditions are met: |
| 9 | +# * Redistributions of source code must retain the above copyright |
| 10 | +# notice, this list of conditions and the following disclaimer. |
| 11 | +# * Redistributions in binary form must reproduce the above copyright |
| 12 | +# notice, this list of conditions and the following disclaimer in the |
| 13 | +# documentation and/or other materials provided with the distribution. |
| 14 | +# * Neither the name of the University of Oxford nor the names of its |
| 15 | +# contributors may be used to endorse or promote products derived from |
| 16 | +# this software without specific prior written permission. |
| 17 | +# |
| 18 | +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND |
| 19 | +# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED |
| 20 | +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
| 21 | +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR |
| 22 | +# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
| 23 | +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; |
| 24 | +# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND |
| 25 | +# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
| 26 | +# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| 27 | +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 28 | + |
| 29 | +# Jared Aurentz and Nick Trefethen, July 2015. |
| 30 | +# |
| 31 | +# Copyright 2017 by The University of Oxford and The Chebfun Developers. |
| 32 | +# See http://www.chebfun.org/ for Chebfun information. |
| 33 | + |
| 34 | +function standardchoplength(coeffs, tol) |
| 35 | + # Check magnitude of TOL: |
| 36 | + if tol ≥ 1 |
| 37 | + throw(ArgumentError("tolerance must be less than 1")) |
| 38 | + end |
| 39 | + |
| 40 | + # Make sure COEFFS has length at least 17: |
| 41 | + n = length(coeffs) |
| 42 | + |
| 43 | + if n < 17 |
| 44 | + # resort to naive chop |
| 45 | + mx = maximum(abs,coeffs) |
| 46 | + return choplength(coeffs, tol*mx) |
| 47 | + end |
| 48 | + |
| 49 | + # Step 1: Convert COEFFS to a new monotonically nonincreasing |
| 50 | + # vector ENVELOPE normalized to begin with the value 1. |
| 51 | + |
| 52 | + b = convert(Vector, abs.(coeffs)) |
| 53 | + for j = n-1:-1:1 |
| 54 | + b[j] = max(b[j], b[j+1]) |
| 55 | + end |
| 56 | + iszero(b[1]) && return 1 |
| 57 | + |
| 58 | + b .= b ./ b[1] |
| 59 | + |
| 60 | + |
| 61 | + plateauPoint = 0 |
| 62 | + |
| 63 | + for j = 2:n |
| 64 | + j2 = round(Int,1.25*j + 5); |
| 65 | + if j2 > n |
| 66 | + # there is no plateau: exit |
| 67 | + return n |
| 68 | + end |
| 69 | + e1 = b[j] |
| 70 | + e2 = b[j2] |
| 71 | + r = 3*(1 - log(e1)/log(tol)) |
| 72 | + plateau = (e1 == 0) || (e2/e1 > r) |
| 73 | + if plateau |
| 74 | + # a plateau has been found: go to Step 3 |
| 75 | + plateauPoint = j - 1 |
| 76 | + break |
| 77 | + end |
| 78 | + end |
| 79 | + |
| 80 | + # Step 3: fix CUTOFF at a point where ENVELOPE, plus a linear function |
| 81 | + # included to bias the result towards the left end, is minimal. |
| 82 | + # |
| 83 | + |
| 84 | + if plateauPoint ≠ 0 && b[plateauPoint] == 0 |
| 85 | + return plateauPoint |
| 86 | + end |
| 87 | + |
| 88 | + j3 = sum(b .≥ tol^(7/6)) |
| 89 | + if j3 < j2 |
| 90 | + j2 = j3 + 1 |
| 91 | + b[j2] = tol^(7/6) |
| 92 | + end |
| 93 | + cc = log10.(b[1:j2]) |
| 94 | + cc .+= range(0, stop=(-1/3)*log10(tol), length=j2) |
| 95 | + d = argmin(cc) |
| 96 | + return max(d - 1, 1) |
| 97 | +end |
0 commit comments