Grid provider for Plain (#77)

* Grid provider for Plain

The GridCollection struct provides access to Arc<dyn Grid>, such that
any grid needed by any Operator in any instantiation of Plain, always
refer to the same heap allocation.

* Better grid handling for the Plain Context

Plain now uses a global BTreeMap for holding
Arc<dyn Grid>, and hands out Arc::clone()s to
Ops needing access. This means that all uses
of a grid refers to the same heap allocation.

Unfortunately, we need to keep the "original"
instantiation in the list, so the Arcs will
never be fully released. To mitigate this, we
empty the BTreeMap every time a new Plain is
instantiated. While not ideal, it limits the
risk of hanging on to needless heap allocations

src/context/plain.rs

@@ -1,7 +1,10 @@
 #[cfg(feature = "with_plain")]
 use crate::authoring::*;
 use crate::grid::ntv2::Ntv2Grid;
-use std::{path::PathBuf, sync::Arc};
+use std::{
+    path::PathBuf,
+    sync::{Arc, Mutex},
+};
 
 // ----- T H E   P L A I N   C O N T E X T ---------------------------------------------
 
@@ -18,6 +21,51 @@ pub struct Plain {
     paths: Vec<std::path::PathBuf>,
 }
 
+// Helper for Plain: Provide grid access for all `Op`s
+// in all instantiations of `Plain` by handing out
+// reference counted clones to a single heap allocation
+static mut GRIDS: Mutex<GridCollection> =
+    Mutex::new(GridCollection(BTreeMap::<String, Arc<dyn Grid>>::new()));
+
+struct GridCollection(BTreeMap<String, Arc<dyn Grid>>);
+impl GridCollection {
+    fn get_grid(&mut self, name: &str, paths: &[PathBuf]) -> Result<Arc<dyn Grid>, Error> {
+        // If the grid is already there, just return a reference clone
+        if let Some(grid) = self.0.get(name) {
+            return Ok(grid.clone());
+        }
+
+        // Otherwise, we must look for it in the data path
+        let n = PathBuf::from(name);
+        let ext = n
+            .extension()
+            .unwrap_or_default()
+            .to_str()
+            .unwrap_or_default();
+
+        for path in paths {
+            let mut path = path.clone();
+            path.push(ext);
+            path.push(name);
+            let Ok(grid) = std::fs::read(path) else {
+                continue;
+            };
+
+            if ext == "gsb" {
+                self.0
+                    .insert(name.to_string(), Arc::new(Ntv2Grid::new(&grid)?));
+            } else {
+                self.0
+                    .insert(name.to_string(), Arc::new(BaseGrid::gravsoft(&grid)?));
+            }
+            if let Some(grid) = self.0.get(name) {
+                return Ok(grid.clone());
+            }
+        }
+        Err(Error::NotFound(name.to_string(), ": Grid".to_string()))
+    }
+}
+
 const BAD_ID_MESSAGE: Error = Error::General("Plain: Unknown operator id");
 
 impl Default for Plain {
@@ -27,6 +75,10 @@ impl Default for Plain {
         let operators = BTreeMap::new();
         let mut paths = Vec::new();
 
+        // To avoid having GRIDS growing through the roof, we clear it
+        // out every time a new Plain context is instantiated
+        unsafe { GRIDS.lock().unwrap().0.clear() }
+
         let localpath: PathBuf = [".", "geodesy"].iter().collect();
         paths.push(localpath);
 
@@ -199,25 +251,10 @@ impl Context for Plain {
 
     /// Access grid resources by identifier
     fn get_grid(&self, name: &str) -> Result<Arc<dyn Grid>, Error> {
-        let n = PathBuf::from(name);
-        let ext = n
-            .extension()
-            .unwrap_or_default()
-            .to_str()
-            .unwrap_or_default();
-        for path in &self.paths {
-            let mut path = path.clone();
-            path.push(ext);
-            path.push(name);
-            let Ok(grid) = std::fs::read(path) else {
-                continue;
-            };
-            if ext == "gsb" {
-                return Ok(Arc::new(Ntv2Grid::new(&grid)?));
-            }
-            return Ok(Arc::new(BaseGrid::gravsoft(&grid)?));
-        }
-        Err(Error::NotFound(name.to_string(), ": Blob".to_string()))
+        // The GridCollection does all the hard work here, but accessing GRIDS,
+        // which is a mutable static is (mis-)diagnosed as unsafe by the compiler,
+        // even though the mutable static is behind a Mutex guard
+        unsafe { GRIDS.lock().unwrap().get_grid(name, &self.paths) }
     }
 }
 
@@ -303,4 +340,18 @@ mod tests {
 
         Ok(())
     }
+
+    #[test]
+    fn grids() -> Result<(), Error> {
+        let mut ctx = Plain::new();
+
+        // Here, we only invoke reference counting in the GridCollection. The tests in
+        // gridshift and deformation makes sure that the correct grids are actually
+        // provided by GridCollection::get_grid()
+        let _op1 = ctx.op("gridshift grids=5458.gsb, 5458_with_subgrid.gsb")?;
+        let _op2 = ctx.op("gridshift grids=5458.gsb, 5458_with_subgrid.gsb")?;
+        let _op3 = ctx.op("gridshift grids=test.geoid")?;
+        assert!(ctx.op("gridshift grids=non.existing").is_err());
+        Ok(())
+    }
 }

src/inner_op/gridshift.rs

@@ -270,19 +270,34 @@ mod tests {
         let mut ctx = Plain::default();
         let op = ctx.op("gridshift grids=@test_subset.datum, @missing.gsb, test.datum")?;
 
+        // Check that the Arc reference counting and the clear-GRIDS-by-instantiation
+        // works together correctly (i.e. the Arcs used by op are kept alive when
+        // GRIDS is cleared on the instantiation of ctx2)
+        let mut ctx2 = Plain::default();
+        let op2 = ctx2.op("gridshift grids=@test_subset.datum, @missing.gsb, test.datum")?;
+
         // Copenhagen is outside of the (optional, but present, subset grid)
         let cph = Coor4D::geo(55., 12., 0., 0.);
         let mut data = [cph];
+        let mut data2 = [cph];
 
         ctx.apply(op, Fwd, &mut data)?;
         let res = data[0].to_geo();
         assert!((res[0] - 55.015278).abs() < 1e-6);
         assert!((res[1] - 12.003333).abs() < 1e-6);
 
+        // Same procedure on the other context gives same result
+        ctx2.apply(op2, Fwd, &mut data2)?;
+        assert_eq!(data, data2);
+
         ctx.apply(op, Inv, &mut data)?;
         assert!((data[0][0] - cph[0]).abs() < 1e-10);
         assert!((data[0][1] - cph[1]).abs() < 1e-10);
 
+        // Same procedure on the other context gives same result
+        ctx2.apply(op2, Inv, &mut data2)?;
+        assert_eq!(data, data2);
+
         // Havnebyen (a small town with a large geodetic installation) is inside the subset grid
         let haby = Coor4D::geo(55.97, 11.33, 0., 0.);
         let mut data = [haby];