Fast approximate Barnes interpolation: illustrated by Python-Numba implementation fast-barnes-py v1.0

Zürcher, Bruno K.

doi:https://doi.org/10.5194/gmd-16-1697-2023

Articles | Volume 16, issue 6

https://doi.org/10.5194/gmd-16-1697-2023

© Author(s) 2023. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/gmd-16-1697-2023

© Author(s) 2023. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 16, issue 6

Development and technical paper

|

27 Mar 2023

Development and technical paper |

| 27 Mar 2023

Fast approximate Barnes interpolation: illustrated by Python-Numba implementation fast-barnes-py v1.0

Bruno K. Zürcher

Download

Final revised paper (published on 27 Mar 2023)
Preprint (discussion started on 19 Jul 2022)

Interactive discussion

Status: closed

RC1:
'Comment on gmd-2022-116', Anonymous Referee #1, 31 Aug 2022

It is a good paper. The idea of approximating Guassian weights using the central limit theorem (CLT) is very good. The Approximation 3 is a happy result mainly because of cancellation of the common factors in the numerator and the denominator. The algorithms for implementation are also nice and thorough. It is surprising to read how the computing time is reduced from a naive Barnes interpolation. Nevertheless, I have some comments and questions, basically for a better presentation.

1. It uses a special kind of the uniform random variable on an interval with a threshold. Hope the author describe why this threshold (sqrt(3/n) sigma) was chosen.

2. page 3, line 54: hope the author provide the definition of the n-fold convolution of p(x) with itself.

3. p5 line 92: hope the author provide more concretely what *^nx and *^yn are, because some readers may know it after reading the Algorithms 2, 3, and 4.

4. In Section 5, the author applied 4-fold convolution. It is known that a summation of uniform random variables converge slowly to the CLT compared to other uni-modal random variables, even though the convergence in distribution and the convergence in pdf are different. Thus I think only 4-fold is too small. I would recommend to apply at least 10-fold convolutions in real applications to ensure a good approximation.

5. The result in this paper is only a nice approximation to Barnes interpolation. I think the title can be changed to 'Fast approximation to Barnes interpolation'.

6. It uses a special kind of the uniform random variable and a numerical integration for convolution in Algorithm 2. It is still good, but I guess one may consider other random variables which may accelerate the convergene to the CLT. Then numerical integration may be a little more complex than that in this paper. This may increase computing time in Algorithm 2, but may be alright with a smaller n.

7. In Algorithms 2 and 3, I was a bit confused with the notations g, r_T, and F[i,j]. Hope the author kindly give short descriptions on these notations to improve reader's understanding.

8. Finally, I think that the idea of approximating Guassian weights using the CLT is very good and so it can be applied to the other area of numerical computations which use Gaussian weights. Hope the author try to search the literature if any body has already published or applied this idea.

Citation: https://doi.org/10.5194/gmd-2022-116-RC1
- AC1: 'Reply on RC1', Bruno Zürcher, 08 Oct 2022
  
  Dear Referee,
  Thank you very much for reading the manuscript carefully and for your valuable comments. Below I give a first response to these comments, point after point.
  1. page 3, line 65: The calculation, why the interval has the length (sqrt(3/n) sigma) is actually not obvious. It is basically chosen in this way, such that the variance of the resulting uniform distribution u_n(x) is (sigma^2/n).
  
  I will add in the manuscript a short comment about that and a calculation that verifies this.
  2. page 3, line 54: I propose to provide the definition of the n-fold convolution of p(x) with itself in an appendix to the manuscript.
  3. page 5, line 92: I agree, the explanation of the operators *^nx and *^yn in the text is minimalistic. I propose to give a thorough derivation of them in the appendix.
  4. Section 5, rate of convergence: Maybe one has to differentiate from apllication to application. When Barnes interpolation is used to visualize data by means of isolines - as done in the manuscript - I noticed that the resulting isolines found more or less their stable form already after applying a 3-fold convolution. For other applications, which require a higher precision, a 10-fold convolution or even more might make sense.
  
  I extend the corresponding remark on p. 18 line 295 accordingly.
  5. I can change the title to something like "Fast Approximative Barnes Interpolation". I did not use the term "approximative" in the title, because other _feasible_ methods to compute Barnes interpolation are approximations as well, but do not emphasize this in an obvious way.
  6. Yes, there are many other PDFs that have a faster CLT-convergence to a Gaussian than the uniform PDF - in the extreme case one could even take a normal distribution itself for which we would have n = 1.
  
  But when using a non-trivial PDF, it is clear that putting and moving the weight window (as described on p. 9 line 161) over the data to be convolved requires in general 2T+1 multiplications and 2T additions per element. Thus, the simplicity of algorithm 2 is lost to a far part and the algorithmic complexity of it grows to O(L*T) instead of O(L). Overall we then could not claim a complexity of O(N+W*H) anymore, this would be rather O(N+W*H*T).
  
  Your comment leads me to consider to add a remark on p. 6 line 128 that Approximation 3 is valid in a much more general context, i.e. also for other (non-uniform) PDFs. In the case of a normal distribution even equality holds.
  7. I try to write this more clear.
  8. This is a good question. I know that the CLT is used by some applications to quickly generate the weights of a normal distribution (as above, of course only in an approximative way). Further I assume that some image processing programs like Photoshop employ this technique in order to apply a Gaussian blur filter on an image. This process can be regarded as a special case of the method presented in this manuscript, where the observation points in question are given by the image matrix. In this case the points are located on a regular grid and the denominator of equation (8) simplifies to 1.
  
  I try to find references for such applications and will mention them in the text of the conclusions section.
  Best regards and thank you once more
  
  Bruno
  
  Citation: https://doi.org/10.5194/gmd-2022-116-AC1
AC2: 'Update of Preprint of 19.07.2022', Bruno Zürcher, 11 Dec 2022

I attach to this comment a revised version of the manuscript that takes the comments of the first reviewer into account.

In particular, the manuscript now includes an appendix that provides explaining details of the used mathematics, which might make the work of further reviewers somewhat easier.

Citation: https://doi.org/10.5194/gmd-2022-116-AC2
AC3: 'Corresponding Difference PDF', Bruno Zürcher, 11 Dec 2022

The difference PDF showing the deltas of the updated preprint version - mentioned in the previous comment - compared to the preprint version of 19.07.2022 can be found in the attachment of this comment.

Citation: https://doi.org/10.5194/gmd-2022-116-AC3
RC2:
'Comment on gmd-2022-116', Anonymous Referee #2, 04 Jan 2023

This is a nice paper. It presents an efficient method to perform a Gaussian smoothing, which takes advantage of the CLT by using the repeated convolution of a computationally cheap (uniform) smoothing to approximate the desired Gaussian result (and further that the 2D result is the product of two 1D calculations).

I see that the author already improved the first version of the manuscript in response to previous comments, and I found it very readable. I found that the questions that arose on my first reading were answered later on. I have a minor suggestion about the presentation of the computational complexity of the main result, which is that the wording could be more precise. For example, in the abstract "When implemented naively, the computational complexity of Barnes interpolation depends directly on both the number of sample points and the number of grid points." This initially confused me as all algorithms (including the proposed one as well as the naive approach) must surely depend directly on these numbers, simply from the time taken to read and write them. It would be better expressed as varying with the *product* of the number of sample points and grid points. Or else (additionally?) write it in O() notation so there can be no confusion. The Abstract should also explicitly present the complexity of the new result (and perhaps both old and new complexities could be placed again in the Conclusions).

Also one piece of ungrammatical English at l133: "are sufficiently good contained" could be "sufficiently well contained" perhaps?

I suggest that the Conclusion contain a repetition of the author's suggestion that 4 convolutions seems to work well in most practical applications.

Citation: https://doi.org/10.5194/gmd-2022-116-RC2
- AC4: 'Reply on RC2', Bruno Zürcher, 07 Jan 2023
  
  Thank you very much for your review, I am very glad to hear that you found it very readable. And also thank you for keeping an eye on my English, I really appreciate that. Below I briefly repeat your suggestions before I respond to them.
  - wording of computational complexity:
  This can actually be written more clearly. As you suggest, I will adapt the corresponding text passages in the Abstract and the Conclusions sections and also explicitly use the O() notation.
  - "well" instead "good":
  Many thanks for proofreading!
  - extend Conclusions with hint to 4 convolutions to achieve good results for most applications:
  As suggested by you, I will include this in the Conclusions.
  
  Citation: https://doi.org/10.5194/gmd-2022-116-AC4

Peer review completion

AR: Author's response | RR: Referee report | ED: Editor decision | EF: Editorial file upload

AR by Bruno Zürcher on behalf of the Authors (14 Jan 2023) Author's response Author's tracked changes Manuscript

ED: Publish subject to technical corrections (16 Jan 2023) by Sylwester Arabas

Dear Author,

Thank you for addressing the reviewers' comments in the revised manuscript and for providing the point-by-point reply. Let me also underline that I am sorry that the review process took longer than expected. Congratulations for excellent reviewers' scoring: in six out of eight received marks the paper was ranked as "Excellent".

Below, I'm listing technical remarks I'm hereby asking to address before moving on with acceptance, typesetting and publication:

1. Please attach a version number for the software release, e.g. v1.0 and ensure that the same version number and license metadata are featured on Zenodo. As of now, the GitHub repository has a release without an explicit number, while the only tag indicates a pre-release: v0.8-alpha.

2. While it is in not required for publication, I highly recommend to consider disseminating the software as a Python package (by including a setup.py, __init__.py, pyproject.toml files). This way, users could easily import the code with proper version identification. This would also enable dissemination of the code on PyPI.org or conda forge package sites.

3. According to the GMD guidelines (https://www.geoscientific-model-development.net/about/manuscript_types.html#item2), the following applies to papers submitted as "Development and technical paper": "If the main intention of an article is to make a general (i.e. model independent) statement about the usefulness of a new development, but the usefulness is shown with the help of one specific model, the model name and version number must be stated in the title." Please thus append to the title a phrase or a subtitle featuring the version number, e.g.: "Fast approximate Barnes interpolation and its Python/Numba implementation fast-barnes-py v1.0" (or alike).

4. Language/typesetting remarks (disclaimer: I'm not a native speaker; note also that as any other accepted GMD manuscript, the paper will go through English copy-editing before galley proofs):
- p1/l15: change "commonly" to "jointly"
- p4/l85: rephrase around "writes" (?), e.g. with "can be expressed as" or "reduces to"
- p5/l88: remove "better"
- p6/l97: comma after "definition"
- p7/l140: comma after "distribution", rephrase around "even equality holds"
- p9/l191: comma before "the outer loop"
- p10/l193: change "After," into "Subsequently, "
- p10/l201: write "NaN" without italics
- p11/l207: change "in the further discourse" to "hereinafter"
- p11/l207: change "in an experimental setup" to "on a dataset" (otherwise unclear if experimental refers to testing or measurements)
- p12/l220: change "reduce the unsteadiness of purely interpreted" to "achieve compiled-code performance using ordinary"
- p12/l225: change "Barnes'" into "Barnes" (as used elsewhere)?
- p13/l236: comma after "In this sense"
- p13/l247: comma after "For smaller grid sizes"
- p16/l283: change "Summarized can be stated that" to "In summary, " or "In summary, there are ... that determine"
- p17/l309: comma after "Therefore"
- p17/l317: change "This finds also its visual correspondence" to "A corresponding feature is visible" (?)
- p18/Table 4: rephrase "Key numbers"
- p18/Table 4: change "in dependency" to "as a function of" or alike
- p19/l327: change "data, which is" into "data, which are" or "dataset which is"
- p19/l333: change "effortful" into "costly" or alike
- p22/l371: remove "intelligent"
- p22/l373: rephrase to "Another error reduction technique which is effective..."
- p23/l389: change "which promise" into "which offers"
- p24/l414: change "earlier" to "in the paper"
- p24/l109: move "Getreuer" out of the parentheses

5. Reading through the final sentences of the conclusions, I found it a bit puzzling to discern which sentences refer to the paper and which refer to outlined future endeavors. Perhaps naming the section "Summary and outlook", and separating the summary from the outlook remarks in a clearer manner would improve flow?

6. Algorithms: please add a number for the algorithm on page 2 and use it when referring to it (all other algorithms are numbered).

7. Figures: please use vector graphics format instead of raster graphics (e.g., by saving to pdf, postscript or svg formats).

8. "Python" remarks:
- p15/l259 (computational cost increases in Python): if the code is compiled by Numba/LLVM, this would not be a Python behavior anymore - please clarify.
- p10/l195 & p12/l221 (potential for concurrent computations): just a comment for future developments (outside of the scope of the paper) - it seems likely that Numba `prange` constructs could facilitate introducing concurrency in the code (see https://numba.pydata.org/numba-doc/latest/user/performance-tips.html#parallel-true) by just changing selected `range` statements into `numba.prange` (akin to OpenMP directives).
- Also, please consider either specifying or allowing users of the package to specify performance-oriented `njit` parameters such as `error_model='numpy'` (instead of the default 'python') and `fastmath=True` (assuming the code does not rely on NaN propagation, etc) - again, this is just a comment, outside of scope of the paper.

9. Appendices: please consider changing the nested numbering in appendices (A: A1, A2, A3) into Appendix A, Appendix B and Appendix C.

10. In the bibliography:
- for Bentley 75, please replace the url with https://doi.org/10.1145/361002.361007
- for Bergthórsson et al., please correct the DOI to: https://doi.org/10.3402/tellusa.v7i3.8902
- for Daley, the DOI refers to a review of the book, and not the book itself
- for Gwosdek et al., please correct DOI to https://doi.org/10.1007/978-3-642-24785-9_38
- for Koppert, please consider adding an URL: https://ams.confex.com/ams/84Annual/techprogram/paper_71789.htm
- for Wells 86, the "dx." in the URL can be removed

Thank you, best regards,
Sylwester Arabas

Hide

AR by Bruno Zürcher on behalf of the Authors (18 Feb 2023) Author's response Manuscript

Short summary

We present a novel algorithm to efficiently compute Barnes interpolation, which is a method for transforming data values recorded at irregularly spaced points into a corresponding regular grid. In contrast to naive implementations with an algorithmic complexity that depends on the product of the number of sample points and the number of grid points, our approach reduces this dependency to their sum.