Repository: UXARRAY/uxarray

Issue:
Incorporating GMPY2 for multiple precision calculation

## Overview
<!--  Please provide a few bullet points summarizing the changes in this PR. This should include
      points on any bug fixes, new functions, or other changes that have been made. -->
This feature introduces multi-precision calculation support in Project UXArray. It offers advantages such as detecting unique coordinates and exact calculation of intersection points, which are crucial in computational geometry. With gmpy2 library, we can achieve higher precision and accuracy in computations, making our project reliable and robust. 

For more information about the multi-precision and exact computation paradigm, please refers to the following materials
[GMPY2](https://gmpy2.readthedocs.io/en/latest/index.html)
[The Exact Computation Paradigm](https://www.cgal.org/exact.html)

## Expected Usage
<!--  If you are adding a feature into the Internal API, please produce a short example of it in action.
      You may ignore this step if it is not applicable (comment out this section). -->
```Python
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
```
## Proposed Progress
**New project structure that incorporates the multi-precision mode**
- [ ] Figure out the API for turning on the multi-precision modes #339 

**Data Read in**
- [x] Support multi-precision when initializing through  #`Grid.__from_vert__()`  #321
- [x] #353

**Grid Pre-processing**
- [x] Precise spherical and cartesian coordinates conversion #326 
- [ ] Build exact Latitude Longitude bound
    - [ ] The customized intervaltree #113 



**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.

```python

# 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](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](https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
- [Precision and robustness in geometric computations](https://link.springer.com/chapter/10.1007/3-540-63818-0_9)

### 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.


Discussion:

Patch:
diff --git a/uxarray/core/grid.py b/uxarray/core/grid.py
index f341d7e9..e3897045 100644
--- a/uxarray/core/grid.py
+++ b/uxarray/core/grid.py
@@ -374,7 +374,10 @@ class Grid:
                 })
 
     # load mesh from a file
-    def __from_ds__(self, dataset):
+    def __from_ds__(self,
+                    dataset,
+                    multi_precision=False,
+                    precision=FLOAT_PRECISION_BITS):
         """Loads a mesh dataset."""
         # call reader as per mesh_type
         if self.mesh_type == "exo":
@@ -395,8 +398,39 @@ class Grid:
         else:
             raise RuntimeError("unknown mesh type")
 
+        # set precision flags
+        self._multi_precision = multi_precision
+        self._precision = precision
+
+        # convert coordinates to gmpy2.mpfr format
+        if self._multi_precision:
+            if "Mesh2_node_y" in self._ds and "Mesh2_node_x" in self._ds:
+                converted_y_cords = convert_to_multiprecision(
+                    self._ds["Mesh2_node_y"].values,
+                    str_mode=False,
+                    precision=self._precision)
+                converted_x_cords = convert_to_multiprecision(
+                    self._ds["Mesh2_node_x"].values,
+                    str_mode=False,
+                    precision=self._precision)
+
+                self._ds["Mesh2_node_y"] = converted_y_cords
+                self._ds["Mesh2_node_x"] = converted_x_cords
+            if not self.islatlon:
+                if "Mesh2_node_z" in self._ds:
+                    converted_z_cords = convert_to_multiprecision(
+                        self._ds["Mesh2_node_z"].values,
+                        str_mode=False,
+                        precision=self._precision)
+                    self._ds["Mesh2_node_z"] = converted_z_cords
+
+            self._ds["Mesh2_node_y"] = converted_y_cords
+            self._ds["Mesh2_node_x"] = converted_x_cords
+
         dataset.close()
 
+        return self._ds
+
     def __repr__(self):
         """Constructs a string representation of the contents of a ``Grid``."""
 
diff --git a/uxarray/utils/constants.py b/uxarray/utils/constants.py
index 0be6cd24..1f337d91 100644
--- a/uxarray/utils/constants.py
+++ b/uxarray/utils/constants.py
@@ -10,4 +10,3 @@ FLOAT_PRECISION_BITS = 53
 
 ENABLE_JIT_CACHE = True
 ENABLE_JIT = True
-


Test patch:
diff --git a/test/test_grid.py b/test/test_grid.py
index 3913f59f..aed910c7 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -209,6 +209,35 @@ class TestGrid(TestCase):
         assert (vgrid.nMesh2_face == 4)
         assert (vgrid.nMesh2_node == 4)
 
+    def test_from_ds(self):
+        exodus_filepath = current_path / "meshfiles" / "exodus" / "outCSne8" / "outCSne8.g"
+
+        xrds_exodus = xr.open_dataset(exodus_filepath)
+
+        # Create the grid with multi-precision off
+        gridOff = ux.Grid(xrds_exodus)
+
+        # Create the grid with multi-precision on
+        gridOn = ux.Grid(xrds_exodus)
+        resultsPrecisionOn = gridOn.__from_ds__(xrds_exodus,
+                                                multi_precision=True,
+                                                precision=128)
+
+        # Get the node values
+        firstNodePrecisionOn = resultsPrecisionOn["Mesh2_node_y"].values[0]
+        firstNodePrecisionOff = gridOff._ds["Mesh2_node_y"].values[0]
+
+        # Cast the floating-point value to mpfr
+        firstNodePrecisionOff_mpfr = gmpy2.mpfr(str(firstNodePrecisionOff))
+
+        # Assert that the values are not equal
+        self.assertTrue(
+            gmpy2.cmp(firstNodePrecisionOff_mpfr, firstNodePrecisionOn) != 0)
+
+        # Check that the flags are set properly
+        self.assertTrue(gridOn._multi_precision)
+        self.assertEqual(gridOn._precision, 128)
+
     def test_init_verts_multi_precision_fill_values(self):
         # Set the desired precision