Bessel zero finding cleanup #596

chris-langfield · 2022-03-15T17:03:17Z

The functions num_besselj_zeros, besselj_zeros, check_besselj_zeros are almost exact ports from the MATLAB code, but are lacking clear comments.

These in turn are modifications of this script by Jonas Lundgren:

https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/26639/versions/5/previews/zerobess.m/index.html

When attempting to work out what exactly these functions are doing to estimate the first few zeros, I came across several tricks that I haven't been able to find a single definitive source for. e.g.,

ASPIRE-Python/src/aspire/basis/basis_utils.py

Lines 165 to 174 in 6cb389e

    
           # Guess first zeros using powers of nu 
        
           c0 = np.array( 
        
               [ 
        
                   [0.1701, -0.6563, 1.0355, 1.8558], 
        
                   [0.1608, -1.0189, 3.1348, 3.2447], 
        
                   [-0.2005, -1.2542, 5.7249, 4.3817], 
        
               ] 
        
           ) 
        
           z0 = nu + c0 @ ((nu + 1) ** np.array([[-1, -2 / 3, -1 / 3, 1 / 3]]).T)

Roughly, this computes initial guesses based on powers of nu, the real number order of the Bessel function. The closest literature I could find for this is in the following paper:

Elbert, Árpád. “Some recent results on the zeros of Bessel functions and orthogonal polynomials.” Journal of Computational and Applied Mathematics 133 (2001): 65-83.

(Section 1.3 starting around eqn 1.7)

The powers as well as the hardcoded values in our code also appear here.

See also: https://www.boost.org/doc/libs/1_63_0/libs/math/doc/html/math_toolkit/bessel/bessel_root.html

It may be a lesser concern to be able to find the theoretical justification here, but it would probably be ideal if we could. Regardless, I'm going to use this PR to comment and clean up these functions as much as possible so that it is easier to understand in the future

codecov · 2022-03-15T17:16:19Z

Codecov Report

Merging #596 (1a6a273) into develop (b7a3398) will increase coverage by 0.89%.
The diff coverage is n/a.

@@             Coverage Diff             @@
##           develop     #596      +/-   ##
===========================================
+ Coverage    86.48%   87.38%   +0.89%     
===========================================
  Files          108      108              
  Lines         8273     8942     +669     
===========================================
+ Hits          7155     7814     +659     
- Misses        1118     1128      +10

Impacted Files	Coverage Δ
src/aspire/basis/basis.py	`87.83% <ø> (ø)`
src/aspire/basis/basis_utils.py	`98.17% <ø> (ø)`
src/aspire/noise/noise.py	`97.50% <0.00%> (-0.84%)`	⬇️
src/aspire/utils/__init__.py	`100.00% <0.00%> (ø)`
src/aspire/covariance/covar.py	`88.88% <0.00%> (+0.09%)`	⬆️
src/aspire/reconstruction/estimator.py	`92.66% <0.00%> (+0.16%)`	⬆️
src/aspire/source/simulation.py	`99.61% <0.00%> (+0.19%)`	⬆️
src/aspire/source/image.py	`97.15% <0.00%> (+0.40%)`	⬆️
src/aspire/ctf/ctf_estimator.py	`99.25% <0.00%> (+0.67%)`	⬆️
src/aspire/volume/__init__.py	`96.48% <0.00%> (+1.74%)`	⬆️
... and 1 more

📣 Codecov can now indicate which changes are the most critical in Pull Requests. Learn more

chris-langfield · 2022-03-15T18:19:59Z

It would be good to have another pair of eyes on this line:

ASPIRE-Python/src/aspire/basis/basis_utils.py

Lines 27 to 31 in 6cb389e

    
           if nu >= 0.5: 
        
               result = result and all(ddz < 16 * np.spacing(z[1:-1])) 
        
           else: 
        
               result = result and all(ddz > -16 * np.spacing(z[1:-1]))

np.spacing returns the distance to the nearest machine-representable float. My understanding of this check is that ddz, the second order differences between the zeros, should be zero up to 16x machine precision.

chris-langfield · 2022-03-16T19:28:11Z

Created #597 for continuing conversation, as we can merge this PR with just the docstrings maybe

j-c-c

Looks good. Much clearer!

src/aspire/basis/basis_utils.py

janden · 2022-03-31T13:31:22Z

It would be good to have another pair of eyes on this line:

ASPIRE-Python/src/aspire/basis/basis_utils.py

Lines 27 to 31 in 6cb389e

if nu >= 0.5:

result = result and all(ddz < 16 * np.spacing(z[1:-1]))

else:

result = result and all(ddz > -16 * np.spacing(z[1:-1]))

np.spacing returns the distance to the nearest machine-representable float. My understanding of this check is that ddz, the second order differences between the zeros, should be zero up to 16x machine precision.

I think it's rather checking that the second-order differences are zero or negative. If they were all zero, we'd have a three-term recurrence formula for the zeros and wouldn't need to find them this way (although maybe that's true for large zeros).

janden · 2022-03-31T13:31:50Z

(although maybe that's true for large zeros)

Sorry I should say asymptotically. It might be a decent approximation.

janden · 2022-03-31T13:40:31Z

Right. Of course. From Zhao & Singer (2013):

So at large q the approximate spacing between zeros is π. In fact, we use this in

ASPIRE-Python/src/aspire/basis/basis_utils.py

Line 193 in 6cb389e

dz = pi * np.ones((j, 1))

to initialize the search for zeros by using the last zero and adding multiples of π when the number of zeros is large enough.

chris-langfield · 2022-03-31T14:23:31Z

Added your comments to #597 and made the docstring changes to this PR. Thanks!

janden

LGTM!

* docstring for besselj_zeros * rest of docstrings * remove pdb import * second order differences comment * comments and docstring for num_besselj_zeros * formatting * tiny fixes; * docstring fixes * NumPy * NumPy again * sanity check * parentheses inside sentence

* Fix `k` bug in BFSReddyChetterji * restrict corr to a disc * Change `self.L` and `self.n` to `self.src.L` and `self.src.n` for estimators (#550) * changed self.L and self.n to self.src.L and self.src.n for estimators * extraneous comma * Add `aspire.utils.coor_trans.crop_2d` and move `coor_trans` objects to second-level API access (#565) * utils.coor_trans methods accessible from second level * crop_2d method * crop2d test * clarify flooring behavior for centering and add apdding * tests finished * isort shenanigans * fix even to odd n->n+1 padding bug * slightly better comment * better behavior for dtype * remove mistaken characters * fixes after merge * last ensure removed * rect tests * padding tests and fill value test * better comments * even better comments * np.diag * crop_pad_2d * clarifying comment when padding * Logical rather than hardcoded downsample tests (#567) * downsample tests * cleanup and deleting test downsample from preprocess_pipeline * remove hardcoded ds / whiten test * flake8 * 3d volume downsample test * remove redundant test from tests/test_preprocess.py * remove unused import * remove numbered tests for uniformity * typo * skip one of the ds tests * use gaussian_3d to generate test volume * isort * checkCenterPoint * max_resolution -> L_ds * parametrize ds tests * Bessel zero finding cleanup (#596) * docstring for besselj_zeros * rest of docstrings * remove pdb import * second order differences comment * comments and docstring for num_besselj_zeros * formatting * tiny fixes; * docstring fixes * NumPy * NumPy again * sanity check * parentheses inside sentence * `Basis` size fixes (#598) * fix size checking bug * fb2d * fb3d and pswf2d * dirac * fpswf * polar * format * sz -> size * test init with int * docstring * typo * the same typo * Authors page grammar fix. * cmap=gray * Gallery verb tense. Change cmap to gray. * Update gallery/experiments/README.rst header language. * Encode size in 2D FB test * Add test for 2D FB Gaussian expansion * Add eval–expand test for 2D FB * Add adjoint test for 2D FB * Add test for isotropic and modulated in 2D FB * Add test for indices * Add test for specific FB 2D basis element * Make 2D FB tests FFB friendly Cast `Images` to arrays when needed. * Move tests not specific to FB into a mixin * Add mixin to FFB 2D tests * Remove hardcoding from complex conversion tests * First attempt at testElement for FFB2D * Remove hardcoded (F)FB2D tests * Parameterize tests * Pass dtype to gaussian_2d * Expand on comment regarding factor of 1/2 * Fix rtoc bug * Add shape checks To make sure we don't get the same inadvertent broadcasting bug again. * Loop over different indices in 2D (F)FB tests * Fix #604 (#607) * Bump version: 0.9.0 → 0.9.1 * Fourier-Bessel Mixin class (#599) Fourier Bessel Mixin Class: fb.py * initial changes to code and tests * FFB bases * dirac * tentative fixes to mean estimator tests * import * covar3d mat_evaluate fix * unnecessary check * comments * black * fb2d check evaluate * fb3d * ffb3d * fb3d comments * dirac * fb2d * ffb2d * importvolume * pswf * fpswf * class checking logic moved to Basis * asnumpy in evaluate_t fb 3d bases * evaluate_t fixes * f string * fixes * fb2d image check consistent with fb3d * import * ISinstance * an Image * fix Dirac basis evaluate_t * dirac basis comments * remove volume check * Revert "remove volume check" This reverts commit 86d920e. * repr for Volume class * right place to squeeze in estimate() * no pdb * adjusted comments * remove new space clutter * clarified kernel_convolve squeeze * some changes * FB evaluate_t inputs are always Images/Volumes, just take transpose of asnumpy * remove dtype warnings * no pdb * diract evaluate_t take images * move repetitive function tests up to basis_util * cleanups * no self.L pswf * mat_evaluate_t * remove dirac * remove dirac tests * pswf * int size test in _basis_util * polar2d * remove image assert * remove type checks * space * spaces * final format * fix self.size typo * comment in mat_evaluate_t * replace DiracBasis with PolarBasis2D in unit test Co-authored-by: Garrett Wright <47759732+garrettwrong@users.noreply.github.com> Co-authored-by: Garrett Wright <garrettwrong@gmail.com> Co-authored-by: Joakim Andén <janden@kth.se> Co-authored-by: Josh Carmichael <carmichael@princeton.edu> Co-authored-by: Josh Carmichael <66811911+j-c-c@users.noreply.github.com>

* Fix `k` bug in BFSReddyChetterji * restrict corr to a disc * Change `self.L` and `self.n` to `self.src.L` and `self.src.n` for estimators (#550) * changed self.L and self.n to self.src.L and self.src.n for estimators * extraneous comma * Add `aspire.utils.coor_trans.crop_2d` and move `coor_trans` objects to second-level API access (#565) * utils.coor_trans methods accessible from second level * crop_2d method * crop2d test * clarify flooring behavior for centering and add apdding * tests finished * isort shenanigans * fix even to odd n->n+1 padding bug * slightly better comment * better behavior for dtype * remove mistaken characters * fixes after merge * last ensure removed * rect tests * padding tests and fill value test * better comments * even better comments * np.diag * crop_pad_2d * clarifying comment when padding * Logical rather than hardcoded downsample tests (#567) * downsample tests * cleanup and deleting test downsample from preprocess_pipeline * remove hardcoded ds / whiten test * flake8 * 3d volume downsample test * remove redundant test from tests/test_preprocess.py * remove unused import * remove numbered tests for uniformity * typo * skip one of the ds tests * use gaussian_3d to generate test volume * isort * checkCenterPoint * max_resolution -> L_ds * parametrize ds tests * Bessel zero finding cleanup (#596) * docstring for besselj_zeros * rest of docstrings * remove pdb import * second order differences comment * comments and docstring for num_besselj_zeros * formatting * tiny fixes; * docstring fixes * NumPy * NumPy again * sanity check * parentheses inside sentence * `Basis` size fixes (#598) * fix size checking bug * fb2d * fb3d and pswf2d * dirac * fpswf * polar * format * sz -> size * test init with int * docstring * typo * the same typo * Authors page grammar fix. * cmap=gray * Gallery verb tense. Change cmap to gray. * Update gallery/experiments/README.rst header language. * Encode size in 2D FB test * Add test for 2D FB Gaussian expansion * Add eval–expand test for 2D FB * Add adjoint test for 2D FB * Add test for isotropic and modulated in 2D FB * Add test for indices * Add test for specific FB 2D basis element * Make 2D FB tests FFB friendly Cast `Images` to arrays when needed. * Move tests not specific to FB into a mixin * Add mixin to FFB 2D tests * Remove hardcoding from complex conversion tests * First attempt at testElement for FFB2D * Remove hardcoded (F)FB2D tests * Parameterize tests * Pass dtype to gaussian_2d * Expand on comment regarding factor of 1/2 * Fix rtoc bug * Add shape checks To make sure we don't get the same inadvertent broadcasting bug again. * Loop over different indices in 2D (F)FB tests * Fix #604 (#607) * Bump version: 0.9.0 → 0.9.1 * Fourier-Bessel Mixin class (#599) Fourier Bessel Mixin Class: fb.py * New downsample method (#606) * dropped in new ds method * adjust tests temporarily * some imports * change to centered_fft2 and allow gridpoint checks * format * remove debug show() * comment * simplify centered_fft2 and centered_ifft2 * dict of dicts * adjusted tests * Zenodo file with authors and affiliations (#615) * add zenodo file with just creators * name * add zenodo file to manifest * json errors * chris orcid * updated affiliations * update tox.ini to include json.tool check * josh orcid * author ordering * review changes * Remove offending sk import * Explicitly use svd_solver='full' for repro unit test * Fix PIL deprecation warning * sensible defaults * optional saving mrc * populate_unique_filters and read_ctf_file * body of populate_unique_filters * loop var * tests * work around ImageSource.set_metadata * update test * final test * rationalize sample CTF values * readme space * correct metadata test * write_all_star * remove write_one_star * method name * reading from Relion micrograph_ctf.star * format * import_ctf() instead of __init__ option * added tests and clarified import_ctf method * condense ctf args coordinate source * .values() instead of items() * add save image flags to increase test converage * docstrings * rest of docstrings * rename fun to populate_ctf * first pass refactor * refactor 2 * docstring * refactor 3 Co-authored-by: Garrett Wright <47759732+garrettwrong@users.noreply.github.com> Co-authored-by: Garrett Wright <garrettwrong@gmail.com> Co-authored-by: Joakim Andén <janden@kth.se> Co-authored-by: Josh Carmichael <carmichael@princeton.edu> Co-authored-by: Josh Carmichael <66811911+j-c-c@users.noreply.github.com>

docstring for besselj_zeros

db7b619

chris-langfield added the cleanup label Mar 15, 2022

chris-langfield self-assigned this Mar 15, 2022

chris-langfield added 2 commits March 15, 2022 14:16

rest of docstrings

3ebaf7c

remove pdb import

69c7cea

chris-langfield added 4 commits March 15, 2022 14:21

second order differences comment

497fd8f

comments and docstring for num_besselj_zeros

d48797c

formatting

6c65fb9

tiny fixes;

936b07d

chris-langfield mentioned this pull request Mar 16, 2022

Theoretical basis of Bessel root-finding #597

Closed

chris-langfield requested a review from j-c-c March 16, 2022 19:32

j-c-c previously approved these changes Mar 18, 2022

View reviewed changes

chris-langfield marked this pull request as ready for review March 18, 2022 15:02

chris-langfield requested a review from janden as a code owner March 18, 2022 15:02

janden requested changes Mar 31, 2022

View reviewed changes

src/aspire/basis/basis_utils.py Outdated Show resolved Hide resolved

src/aspire/basis/basis_utils.py Outdated Show resolved Hide resolved

src/aspire/basis/basis_utils.py Outdated Show resolved Hide resolved

src/aspire/basis/basis_utils.py Outdated Show resolved Hide resolved

chris-langfield added 2 commits March 31, 2022 09:33

docstring fixes

7fd28f8

NumPy

95152b3

chris-langfield dismissed j-c-c’s stale review via 95152b3 March 31, 2022 13:47

chris-langfield added 3 commits March 31, 2022 09:48

NumPy again

6828481

sanity check

a3ac8b1

parentheses inside sentence

1a6a273

chris-langfield requested a review from janden March 31, 2022 14:23

janden approved these changes Mar 31, 2022

View reviewed changes

chris-langfield merged commit 01204af into develop Mar 31, 2022

chris-langfield deleted the basis_utils_cleanup branch March 31, 2022 17:39

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Bessel zero finding cleanup #596

Bessel zero finding cleanup #596

Uh oh!

chris-langfield commented Mar 15, 2022 •

edited

Loading

Uh oh!

codecov bot commented Mar 15, 2022 •

edited

Loading

Uh oh!

chris-langfield commented Mar 15, 2022

Uh oh!

chris-langfield commented Mar 16, 2022

Uh oh!

j-c-c left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

janden commented Mar 31, 2022

Uh oh!

janden commented Mar 31, 2022

Uh oh!

janden commented Mar 31, 2022

Uh oh!

chris-langfield commented Mar 31, 2022

Uh oh!

janden left a comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

	# Guess first zeros using powers of nu
	c0 = np.array(
	[
	[0.1701, -0.6563, 1.0355, 1.8558],
	[0.1608, -1.0189, 3.1348, 3.2447],
	[-0.2005, -1.2542, 5.7249, 4.3817],
	]
	)
	z0 = nu + c0 @ ((nu + 1) ** np.array([[-1, -2 / 3, -1 / 3, 1 / 3]]).T)

Bessel zero finding cleanup #596

Bessel zero finding cleanup #596

Uh oh!

Conversation

chris-langfield commented Mar 15, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov bot commented Mar 15, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

chris-langfield commented Mar 15, 2022

Uh oh!

chris-langfield commented Mar 16, 2022

Uh oh!

j-c-c left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

janden commented Mar 31, 2022

Uh oh!

janden commented Mar 31, 2022

Uh oh!

janden commented Mar 31, 2022

Uh oh!

chris-langfield commented Mar 31, 2022

Uh oh!

janden left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

chris-langfield commented Mar 15, 2022 •

edited

Loading

codecov bot commented Mar 15, 2022 •

edited

Loading