git subrepo pull src/vgrid

subrepo:
  subdir:   "src/vgrid"
  merged:   "aff34f4"
upstream:
  origin:   "[email protected]:jreniel/schismrs-vgrid.git"
  branch:   "main"
  commit:   "aff34f4"
git-subrepo:
  version:  "0.4.5"
  origin:   "https://github.com/ingydotnet/git-subrepo"
  commit:   "dbb99be"

Author: jreniel

src/vgrid/.gitrepo

@@ -6,7 +6,7 @@
 [subrepo]
 	remote = [email protected]:jreniel/schismrs-vgrid.git
 	branch = main
-	commit = 6c8ef6b58fb355fb8e7b6b4ab070be57e4adf18b
-	parent = 2969634fb366a60e1c695458be74af57634dceb1
+	commit = aff34f42690ee15e2e1770a89fce82743bbf0e2b
+	parent = a2c13d1087c0e22b7a6d63bacefc68b5f53aa259
 	method = merge
 	cmdver = 0.4.5

src/vgrid/build.rs

@@ -1,4 +1,3 @@
-
 use std::error::Error;
 use vergen::EmitBuilder;
 

src/vgrid/src/bin/gen_sz.rs

@@ -11,34 +11,72 @@ const VERSION: &'static str = concat! {
     env! {"VERGEN_GIT_DESCRIBE"}
 };
 
+fn greater_than_two(s: &str) -> Result<usize, String> {
+    let value: usize = s
+        .parse()
+        .map_err(|_| format!("`{}` isn't a valid positive integer", s))?;
+
+    if value >= 2 {
+        Ok(value)
+    } else {
+        Err(format!(
+            "The value must be greater or equal than 2, got {}",
+            value
+        ))
+    }
+}
+
 #[derive(Parser, Debug)]
 #[command(author, about, long_about = None)]
 #[command(version = VERSION)]
 struct Cli {
     hgrid_path: PathBuf,
     #[clap(short, long)]
     output_filepath: Option<PathBuf>,
-    #[clap(long)]
-    slevels: usize,
-    #[clap(long, value_delimiter = ' ', num_args = 1..)]
+    #[clap(
+        long,
+        help = "Number of sigma-levels. Must be an integer greater or equal than 0.",
+        value_parser=clap::builder::ValueParser::new(greater_than_two),
+        default_value = "2",
+        )
+    ]
+    slevels: Option<usize>,
+    #[clap(long, value_delimiter = ' ', num_args = 1..,
+        help = "Space delimited list of depths for each z-level. \
+                Optional. Defaults to pure sigma grid."
+        )]
     zlevels: Option<Vec<f64>>,
+
     #[clap(
         long,
+        default_value = "0.1",
         help = "Range is (0., 20.]. Values closer to 0. make the transformation \
                 more similar to traditional sigma. Larger values will increase \
                 resolution at the top and bottom."
     )]
-    theta_f: f64,
+    theta_f: Option<f64>,
+
     #[clap(
         long,
+        default_value = "0.",
         help = "Range is [0., 1.]. For values closer to 0. the surface is \
                 resolved. For values closer to 1., both the surface and bottom \
                 are resolved."
     )]
-    theta_b: f64,
-    #[clap(long, help = "Critical layer depth. Value must be > 5.", alias = "hc")]
+    theta_b: Option<f64>,
+    #[clap(
+        long,
+        help = "Critical layer depth. Value must be > 5.",
+        alias = "hc",
+        default_value = "5."
+    )]
     critical_depth: f64,
-    #[clap(short, long, default_value = "0.")]
+    #[clap(
+        short,
+        long,
+        default_value = "0.",
+        help = "Water level offset. Not typically needed."
+    )]
     etal: Option<f64>,
     #[clap(long, action)]
     show_plot: bool,
@@ -52,9 +90,9 @@ fn entrypoint() -> Result<(), Box<dyn Error>> {
     let hgrid = Hgrid::try_from(&cli.hgrid_path)?;
     let mut builder = SZBuilder::default();
     builder.hgrid(&hgrid);
-    builder.slevels(&cli.slevels);
-    builder.theta_f(&cli.theta_f);
-    builder.theta_b(&cli.theta_b);
+    builder.slevels(&cli.slevels.as_ref().unwrap());
+    builder.theta_f(&cli.theta_f.as_ref().unwrap());
+    builder.theta_b(&cli.theta_b.as_ref().unwrap());
     builder.critical_depth(&cli.critical_depth);
     builder.etal(&cli.etal.as_ref().unwrap());
     if cli.zlevels.is_some() {
@@ -63,6 +101,8 @@ fn entrypoint() -> Result<(), Box<dyn Error>> {
     let sz = builder.build()?;
     if cli.output_filepath.is_some() {
         sz.write_to_file(&cli.output_filepath.as_ref().unwrap())?;
+    } else {
+        println!("{}", sz)
     };
 
     if cli.show_plot || cli.save_plot.is_some() {

src/vgrid/src/bin/gen_vqs.rs

@@ -78,12 +78,21 @@ enum Modes {
 
 #[derive(Args, Debug)]
 struct KmeansCliOpts {
-    #[clap(short, long)]
+    #[clap(short, long, help = "Number of clusters. Must be an interger >= 1")]
     clusters: usize,
-    #[clap(short, long, default_value = "2")]
+    #[clap(
+        short,
+        long,
+        default_value = "2",
+        help = "Controls the initial number of layers. Must be an integer >= 2."
+    )]
     shallow_levels: Option<usize>,
-    #[clap(long)]
-    max_levels: usize,
+    #[clap(
+        long,
+        help = "Controls the maximum number of layers in the clustering hierarchy. \
+                Defaults to shallow_levels + clusters - 1"
+    )]
+    max_levels: Option<usize>,
 }
 
 #[derive(Args, Debug)]
@@ -96,18 +105,27 @@ struct HsmCliOpts {
 
 #[derive(Args, Debug)]
 struct AutoCliOpts {
-    #[clap(long)]
+    #[clap(long, help = "Number of master grids to generate. Must be an int >= 1")]
     ngrids: usize,
     #[clap(
         long,
         default_value = "1.",
         help = "This is the first depth below etal. This input is positive down."
     )]
     initial_depth: Option<f64>,
-    #[clap(long, default_value = "2")]
+    #[clap(
+        short,
+        long,
+        default_value = "2",
+        help = "Controls the initial number of layers. Must be an integer >= 2."
+    )]
     shallow_levels: Option<usize>,
-    #[clap(long)]
-    max_levels: usize,
+    #[clap(
+        long,
+        help = "Controls the maximum number of layers in the clustering hierarchy. \
+                Defaults to shallow_levels + clusters - 1"
+    )]
+    max_levels: Option<usize>,
 }
 
 fn entrypoint() -> Result<(), Box<dyn Error>> {
@@ -151,7 +169,9 @@ fn entrypoint() -> Result<(), Box<dyn Error>> {
             if let Some(shallow_levels) = &opts.shallow_levels {
                 builder.shallow_levels(shallow_levels);
             }
-            builder.max_levels(&opts.max_levels);
+            if let Some(max_levels) = &opts.max_levels {
+                builder.max_levels(max_levels);
+            }
             builder.build()?
         }
         Modes::Auto(opts) => {
@@ -162,7 +182,9 @@ fn entrypoint() -> Result<(), Box<dyn Error>> {
             builder.dz_bottom_min(&cli.dz_bottom_min);
             builder.initial_depth(&opts.initial_depth.as_ref().unwrap());
             builder.shallow_levels(&opts.shallow_levels.as_ref().unwrap());
-            builder.max_levels(&opts.max_levels);
+            if let Some(max_levels) = &opts.max_levels {
+                builder.max_levels(max_levels);
+            }
             builder.build()?
         }
     };

src/vgrid/src/sz.rs

@@ -100,13 +100,11 @@ pub struct SZBuilder<'a> {
 
 impl<'a> SZBuilder<'a> {
     pub fn build(&self) -> Result<SZ, SZBuilderError> {
-        let hgrid = self
-            .hgrid
-            .ok_or_else(|| SZBuilderError::UninitializedFieldError("hgrid".to_string()))?;
         let slevels = self
             .slevels
             .as_ref()
             .ok_or_else(|| SZBuilderError::UninitializedFieldError("slevels".to_string()))?;
+        Self::validate_s_levels(slevels)?;
         let theta_f = self
             .theta_f
             .as_ref()
@@ -126,6 +124,9 @@ impl<'a> SZBuilder<'a> {
             .etal
             .as_ref()
             .ok_or_else(|| SZBuilderError::UninitializedFieldError("etal".to_string()))?;
+        let hgrid = self
+            .hgrid
+            .ok_or_else(|| SZBuilderError::UninitializedFieldError("hgrid".to_string()))?;
         let depths = hgrid.depths();
         let deepest_point = depths.min()?;
         let z_array: Array1<f64> = match &self.zlevels {
@@ -164,6 +165,12 @@ impl<'a> SZBuilder<'a> {
         };
         Ok(())
     }
+    fn validate_s_levels(s_levels: &usize) -> Result<(), SZBuilderError> {
+        if *s_levels < 2 {
+            return Err(SZBuilderError::InvalidSLevels);
+        };
+        Ok(())
+    }
     fn validate_z_levels(deepest_point: &f64, zlevels: &Vec<f64>) -> Result<(), SZBuilderError> {
         if !zlevels.iter().all(|&val| val <= 0.0) {
             return Err(SZBuilderError::InvalidZLevels);
@@ -225,6 +232,8 @@ pub enum SZBuilderError {
     MinMaxError(#[from] ndarray_stats::errors::MinMaxError),
     #[error("zlevels must be all negative and increasing")]
     InvalidZLevels,
+    #[error("slevels must be >= 2")]
+    InvalidSLevels,
     #[error("The first point of zlevels must be smaller or equal to the deepest point in the mesh ({0}) but got {1}")]
     InvalidZLevelsValues(f64, f64),
     #[error("theta_b must be in [0., 1.], but got {0}")]

src/vgrid/src/vqs.rs

@@ -360,17 +360,22 @@ impl<'a> VQSKMeansBuilder<'a> {
             VQSKMeansBuilderError::UninitializedFieldError("shallow_levels".to_string())
         })?;
         Self::validate_shallow_levels(shallow_levels)?;
-        let max_levels = self.max_levels.ok_or_else(|| {
-            VQSKMeansBuilderError::UninitializedFieldError("max_levels".to_string())
-        })?;
-        Self::validate_max_levels(shallow_levels, max_levels)?;
+        // let max_levels = self.max_levels.ok_or_else(|| {
+        //     VQSKMeansBuilderError::UninitializedFieldError("max_levels".to_string())
+        // })?;
+        let max_levels = match self.max_levels {
+            Some(max_levels) => *max_levels,
+            None => Self::calculate_max_levels(shallow_levels, nclusters),
+        };
+        Self::validate_max_levels(shallow_levels, &max_levels)?;
+
         let dz_bottom_min = self.dz_bottom_min.ok_or_else(|| {
             VQSKMeansBuilderError::UninitializedFieldError("dz_bottom_min".to_string())
         })?;
         let mut hsm = kmeans_hsm(hgrid, nclusters, etal)?;
         hsm.iter_mut().for_each(|depth| *depth = depth.abs());
         let mut nlevels = Vec::<usize>::with_capacity(*nclusters);
-        let levels = Array::linspace(*shallow_levels as f64, *max_levels as f64, *nclusters);
+        let levels = Array::linspace(*shallow_levels as f64, max_levels as f64, *nclusters);
         for level in levels.iter() {
             let mut level = level.round() as usize;
             if level < *shallow_levels {
@@ -421,9 +426,13 @@ impl<'a> VQSKMeansBuilder<'a> {
         }
         Ok(())
     }
+
+    fn calculate_max_levels(shallow_levels: &usize, clusters: &usize) -> usize {
+        shallow_levels + clusters - 1
+    }
     fn validate_max_levels(
-        shallow_levels: &'a usize,
-        max_levels: &'a usize,
+        shallow_levels: &usize,
+        max_levels: &usize,
     ) -> Result<(), VQSKMeansBuilderError> {
         if *shallow_levels > *max_levels {
             return Err(VQSKMeansBuilderError::InvalidMaxLevels(
@@ -484,12 +493,17 @@ impl<'a> VQSAutoBuilder<'a> {
             VQSAutoBuilderError::UninitializedFieldError("shallow_levels".to_string())
         })?;
         Self::validate_shallow_levels(shallow_levels)?;
-        let max_levels = self.max_levels.ok_or_else(|| {
-            VQSAutoBuilderError::UninitializedFieldError("max_levels".to_string())
-        })?;
-        Self::validate_max_levels(shallow_levels, max_levels)?;
+        let max_levels = match self.max_levels {
+            Some(max_levels) => *max_levels,
+            None => Self::calculate_max_levels(shallow_levels, ngrids),
+        };
+        Self::validate_max_levels(shallow_levels, &max_levels)?;
+        // let max_levels = self.max_levels.ok_or_else(|| {
+        //     VQSAutoBuilderError::UninitializedFieldError("max_levels".to_string())
+        // })?;
+        // Self::validate_max_levels(shallow_levels, max_levels)?;
         let (hsm, nlevels) =
-            Self::build_hsm_and_nlevels(hgrid, ngrids, initial_depth, shallow_levels, max_levels)?;
+            Self::build_hsm_and_nlevels(hgrid, ngrids, initial_depth, shallow_levels, &max_levels)?;
         Ok(VQSBuilder::default()
             .hgrid(&hgrid)
             .depths(&hsm)
@@ -525,6 +539,10 @@ impl<'a> VQSAutoBuilder<'a> {
         Ok(())
     }
 
+    fn calculate_max_levels(shallow_levels: &usize, clusters: &usize) -> usize {
+        shallow_levels + clusters - 1
+    }
+
     fn exponential_samples(start: f64, end: f64, steps: usize) -> Vec<f64> {
         let mut samples = Vec::with_capacity(steps);
         let scale = (end / start).powf(1.0 / (steps as f64 - 1.0));