Arbitrary-Precision Arithmetic in UXarray

[hongyuchen1030]

Overview

UXarray allows for the reading of grid files and provides functionality that lets users calculate and modify grids. Multi-precision allows this same functionality, but instead of having a fixed precision of 53 bits which is the default in Python, it lets us choose the precision we want. The proposal is to overload the algorithms with multi-precision. Then the user can choose whether to use it or completely ignore it and stick with the default 53-bit precision. It would be as simple as

ux.Grid(face_nodes_connectivity, multi_precision=True, precision= 65)

if users want it on, or

ux.Grid(face_nodes_connectivity) 

if the user would like to stick with the default value.

More resources for information can be found here:
GMPY2
The Exact Computation Paradigm
Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates
Precision and robustness in geometric computations

Expected Usage

import uxarray as ux

# Provide nodes that are extremely closed to each other
face_nodes_connectivity = [['1.000000000000000000', '1.000000000000000000001', '1.0000000000000000000012'],
 ['1.000000000000000001', '1.000000000000000000002', '1.0000000000000000000013']]

uxgrid =   ux.Grid(face_nodes_connectivity, multi_precision=True, precision= 65,
                        vertices=True,
                        islatlon=False,
                        concave=False)
print(uxgrid.nMesh2_face )# It will still detect two faces without degenration
>>> 2

# Perform any provided grid manipulation and then the results is precise up the the 65 bits

Why Arbitrary-Precision

Using standard default (53-bit) precision works in most cases for grids that aren’t too complex. However, when we have grids that include points that are extremely close to one another, we run into problems with the default precision. Take the three points [‘1.0000000000000’, ‘1.0000000000001’, ‘1.0000000000012’]. In default precision, these points would all be read as equal. This causes problems with precise grids, however. A singular case of this is not a big deal, however as this problem repeats over and over, it builds up to make inaccurate calculations and grids. Thus, having the ability to specify the exact number of bit precision in a given grid instead of just using the default value would allow users to control the precision of their grids when they have points requiring such use.

Roadmap / Goals

• Data Read in
  • Support for multi-precision through Grid.__from_vert__()
    - Converts user input coordinates to any specified precision using gmpy2.mpfr
  • Support for multi-precision through Grid.__from_ds__()
• Grid Pre-processing
  • Precise spherical and cartesian coordinates conversion
  • Build exact Latitude Longitude bound
• Grid Operation
  • Exact intersection between two great circle arcs (GCR)
  • Exact intersection between a GCR and a constant latitude
  • Exact Non-conservative Zonal Average
• Grid Construction
  • Exact conservative zonal regrinding

Q&A

1. Q: From 53 to 60 bit precision - is that necessary to use another library for this?

A: The Python float type is limited to 53 bits of precision, which means we can't go beyond that using standard libraries. By integrating gmpy2, designed for multi-precision arithmetic, we can surpass the 53-bit limit and achieve higher precision.

While the PR mentions a precision increase to 60 bits as an example, multi-precision support allows for arbitrary precision settings. We can choose any desired value, such as 100 bits or 200 bits, depending on project requirements. The 60-bit example is just illustrative.

Enabling multi-precision support facilitates not only higher precision calculations but also exact computation. The exact computation is crucial in computational geometry to ensure precise and accurate results, especially when rounding errors and numerical instability could impact the outcome.

This feature allows for higher precision calculations and provides the capability for exact computation. Exact computation is an important paradigm in computational geometry as it ensures precise and accurate results, particularly in cases where rounding errors and numerical instability can affect the outcome.

Incorporating multi-precision support isn't just about adding another package. It means embracing the exact computation paradigm, providing robust analysis for both the geoscience and computational geometry communities.

2. Q: Can you say a little more about how this will be used in practice? Normally everything in Uxarray will be in default python precision but sometimes a calculation must be done in increased precision. How do we decide and what has to happen to "convert" to using this library?

A: So by default, our UXArray project uses the python float for all calculations. But if users want to achieve the arbitrary-precision, they can set the multi-precision flag as True and the desired precision when initializing the Grid object. Then all operations in the UXarray will be done with the desired precision.

# The default floating point Grid object
float_grid = ux.Grid(verts) 

# The multiplrecision Grid object that's set to the 100 precision
multi_precision_grid = ux.Grid(verts, multiprecision=True, precision=100)

In other words, we never mix the usage of multi-precision and the python float, if the grid object is set to default float operation, all the operations will be in floating points. If the grid object is set to multi-precision, then all operations will be in arbitrary-precision. And we can always round the multi-precision number back to the float precision with different rounding options

With this feature, we can achieve

• The unique coordinate determination: We can distinguish two coordinates even though they're extremely close to each other
• The exact intersection points between GCRs or GCR and constant latitude
• We can also achieve exact computations of conservative zonal re-gridding

3. Q: Other than initializing a grid from vertices, what other applications would there be for this? Is there good references/papers about this ideas?

A: The application of the exact computation of the https://www.cgal.org/exact.html will be not only useful in determining unique coordinates but also in

• Determining is a point really on that GCR? (What if it just 1.0e-20 magnitude away from the GCR?). A single point doesn't seems important here, but what if there'are 200000 points that each one of it is 1.0e-20 from the last one? The error will accumulate

When dealing with large datasets, even small errors can accumulate and lead to significant discrepancies. Using exact computation can help ensure that each point is accurately placed on the GCR.

• Get the exact intersection point between two GCRs: Are you sure such an intersection point is exactly on both GCRs. not just 1.0e-20 away from the 1st GCR? And if I calculate this point from the 1st GCR, it will yield the exact same result as when I calculate it from the 2nd GCR?

Without exact computation, it's possible that the intersection point may be off by a small margin, which can cause issues in certain applications. By using exact computation, you can be sure that the intersection point is exactly on both GCRs and that the same result will be obtained regardless of which GCR is used to calculate it.

If you're interested, you can also refer to the following:

• Adaptive Precision Floating-Point Arithmetic
  and Fast Robust Geometric Predicates
• Precision and robustness in geometric computations

4. Q: What happens if you want to do an operation between two grids, like an intersection, and one is default precision and one is arbitrary? As for the performance. Do operations on a grid declared with arbitrary precision take a lot more time?

A: I have tested the performance difference on the coordinates conversion so far, and here're the results summary:

Using the gmpy2 library, precision can be increased by 2 to 3 orders of magnitude compared to the standard float64 type. This significant improvement in precision comes at the cost of only 1 to 2 orders of magnitude increase in runtime.