From b6e7a8026a0a92854ad259c0439428ee6afda1e6 Mon Sep 17 00:00:00 2001
From: Thomas Knudsen <knudsen.thomas@gmail.com>
Date: Fri, 24 Nov 2023 14:47:27 +0100
Subject: [PATCH] KP: Guess output dimensionality

Make KP guess the required number of output dimensions, by matching
the number of input dimensions
---
 ruminations/003-rumination.md | 78 +++++++++++++++++++----------------
 src/bin/kp.rs                 | 49 ++++++++++++++++++----
 2 files changed, 83 insertions(+), 44 deletions(-)

diff --git a/ruminations/003-rumination.md b/ruminations/003-rumination.md
index adde94d..b685c08 100644
--- a/ruminations/003-rumination.md
+++ b/ruminations/003-rumination.md
@@ -4,13 +4,13 @@
 
 Thomas Knudsen <knudsen.thomas@gmail.com>
 
-2021-08-28. Last [revision](#document-history) 2023-11-20
+2021-08-28. Last [revision](#document-history) 2023-11-24
 
 ### Abstract
 
-```sh
+```console
 $ echo 55 12 | kp "geo:in | utm zone=32"
-> 691875.6321 6098907.8250 0.0000 0.0000
+> 691875.6321 6098907.8250
 ```
 
 ---
@@ -27,20 +27,23 @@ For many years, Poder was in charge of the GI department for computational geode
 
 The basic operation of `kp` is very simple. Any complexity in `kp` usage is related to the description of the operation to carry out, which is the subject of [Rumination 002](/ruminations/002-rumination.md). The `kp` command line syntax is:
 
-```sh
-kp "operation" file1 file2 ...
+```console
+$ kp "operation" file1 file2 ...
+> ...
 ```
 
 or, with input from `stdin`:
 
-```sh
-echo coordinate | kp "operation"
+```console
+$ echo coordinate | kp "operation"
+> ...
 ```
 
 or, with output to the file `result`:
 
-```sh
-kp -o result "operation" file1 file2 ...
+```console
+$ kp -o result "operation" file1 file2 ...
+> ...
 ```
 
 ### Examples
@@ -49,39 +52,40 @@ Convert the coordinate tuple (55 N, 12 E) from geographical coordinates  on the
 
 ```sh
 $ echo 55 12 0 0 | kp "geo:in | utm zone=32"
-> 691875.6321 6098907.8250 0.0000 0.0000
+> 691875.6321 6098907.8250
 ```
 
-While RG coordinates are always 4D, `kp` will provide zero-values for any left-out postfix dimensions:
+While RG coordinates are always 4D, `kp` will provide zero-values for any left-out postfix dimensions, and try to guess a proper number of output dimensions (unless the `-D n` option is given):
 
-```sh
+```console
 $ echo 55 12 | kp "geo:in | utm zone=32"
-> 691875.6321 6098907.8250 0.0000 0.0000
+> 691875.6321 6098907.8250
+
+$ echo 55 12 | kp -D3 "geo:in | utm zone=32"
+> 691875.6321 6098907.8250 0.0000
+
+$ echo 55 | kp "curvature mean"
+> 6385431.75306
+
+$ echo 55 | kp -D4 "curvature mean"
+> 6385431.75306 0.00000 0.00000 NaN
 ```
 
 The `roundtrip` option measures the roundtrip accuracy of a transformation
 (i.e. how close to the origin you end up after a forward+inverse dance). Knud Poder championed this practise with his ingeniously constructed *Poder dual autochecking* method, which was essential at a time where computers were less robust than today (more about that [below](#a-few-more-words-about-knud-poder)).
 
-```sh
+```console
 $ echo 55 12 | kp --roundtrip "geo:in | utm zone=32"
-> 55 12:  d = 0.05 mm
+> -0.0000000000 -0.0000000000
 ```
 
 The `inv` option runs the specified pipeline inversely:
 
-```sh
+```console
 $ echo 691875.6321 6098907.8250 | kp --inv "geo:in | utm zone=32"
 > 54.9999999996 11.9999999994 0.00000 0.00000
 ```
 
-The `inv` and `roundtrip` options are mutually exclusive:
-
-```txt
-$ echo 691875.6321 6098907.8250 | kp --inv --roundtrip "geo:in | utm zone=32"
-> Options `inverse` and `roundtrip` are mutually exclusive
-> error: process didn't exit successfully: ...
-```
-
 ### Options
 
 The `help` option gives the list of options:
@@ -98,18 +102,19 @@ Arguments:
   [ARGS]...    The files to operate on
 
 Options:
-      --inv                  Inverse operation
-  -z, --height <HEIGHT>      Specify a fixed height for all coordinates
-  -t, --time <TIME>          Specify a fixed observation time for all coordinates
-  -d, --decimals <DECIMALS>
-      --debug                Activate debug mode
-  -r, --roundtrip            Report fwd-inv roundtrip deviation
-  -e, --echo                 Echo input to output
-  -v, --verbose...           More output per occurrence
-  -q, --quiet...             Less output per occurrence
-  -o, --output <OUTPUT>      Output file, stdout if not present
-  -h, --help                 Print help
-  -V, --version              Print version
+      --inv                    Inverse operation
+  -z, --height <HEIGHT>        Specify a fixed height for all coordinates
+  -t, --time <TIME>            Specify a fixed observation time for all coordinates
+  -d, --decimals <DECIMALS>    Number of decimals in output
+  -D, --dimension <DIMENSION>  Output dimensionality - default: Estimate from input
+      --debug                  Activate debug mode
+  -r, --roundtrip              Report fwd-inv roundtrip deviation
+  -e, --echo                   Echo input to output
+  -v, --verbose...             More output per occurrence
+  -q, --quiet...               Less output per occurrence
+  -o, --output <OUTPUT>        Output file, stdout if not present
+  -h, --help                   Print help
+  -V, --version                Print version
 ```
 
 ### Operators
@@ -157,3 +162,4 @@ Major revisions and additions:
 - 2022-05-08: Reflect current syntax
 - 2023-08-17: Graphical clean up
 - 2023-11-20: Reflect the current --help text
+- 2023-11-24: Automatic selection of output dimensionality
diff --git a/src/bin/kp.rs b/src/bin/kp.rs
index a2e2683..9e03e20 100644
--- a/src/bin/kp.rs
+++ b/src/bin/kp.rs
@@ -30,9 +30,14 @@ struct Cli {
     #[clap(short = 't', long)]
     time: Option<f64>,
 
+    /// Number of decimals in output
     #[clap(short = 'd', long)]
     decimals: Option<usize>,
 
+    /// Output dimensionality - default: Estimate from input
+    #[clap(short = 'D', long)]
+    dimension: Option<usize>,
+
     /// Activate debug mode
     #[clap(long)]
     debug: bool,
@@ -90,6 +95,7 @@ fn main() -> Result<(), anyhow::Error> {
     // Get ready to read and transform input data
     let mut number_of_operands_read = 0_usize;
     let mut number_of_operands_succesfully_transformed = 0_usize;
+    let mut number_of_dimensions_in_input = 0;
     let mut operands = Vec::new();
     let start = time::Instant::now();
 
@@ -119,6 +125,8 @@ fn main() -> Result<(), anyhow::Error> {
                 continue;
             }
 
+            number_of_dimensions_in_input = number_of_dimensions_in_input.max(n);
+
             // Convert the text representation to a Coor4D
             args.extend(&(["0", "0", "0", "NaN", "0"][args.len()..]));
             let mut b: Vec<f64> = vec![];
@@ -136,15 +144,26 @@ fn main() -> Result<(), anyhow::Error> {
             // on to the transformation factory every time, we have
             // 25000 operands to operate on
             if operands.len() == 25000 {
-                number_of_operands_succesfully_transformed +=
-                    transform(&options, op, &mut operands, &ctx)?;
+                number_of_operands_succesfully_transformed += transform(
+                    &options,
+                    op,
+                    number_of_dimensions_in_input,
+                    &mut operands,
+                    &ctx,
+                )?;
                 operands.truncate(0);
             }
         }
     }
 
     // Transform the remaining coordinates
-    number_of_operands_succesfully_transformed += transform(&options, op, &mut operands, &ctx)?;
+    number_of_operands_succesfully_transformed += transform(
+        &options,
+        op,
+        number_of_dimensions_in_input,
+        &mut operands,
+        &ctx,
+    )?;
 
     let duration = start.elapsed();
     info!("Read {number_of_operands_read} coordinates and succesfully transformed {number_of_operands_succesfully_transformed} in {duration:?}");
@@ -152,13 +171,15 @@ fn main() -> Result<(), anyhow::Error> {
     Ok(())
 }
 
+// Transformation - this is the actual geodetic content
 fn transform(
     options: &Cli,
     op: OpHandle,
+    number_of_dimensions_in_input: usize,
     operands: &mut Vec<Coor4D>,
     ctx: &Plain,
 ) -> Result<usize, geodesy::Error> {
-    // Transformation - this is the actual geodetic content
+    let output_dimension = options.dimension.unwrap_or(number_of_dimensions_in_input);
 
     // When roundtripping, we must keep a copy of the input to be able
     // to compute the roundtrip differences
@@ -208,10 +229,22 @@ fn transform(
 
     // Finally output the transformed coordinates
     for coord in operands {
-        println!(
-            "{1:.0$} {2:.0$} {3:.0$} {4:.0$} ",
-            decimals, coord[0], coord[1], coord[2], coord[3]
-        );
+        match output_dimension {
+            0 | 4 => println!(
+                "{1:.0$} {2:.0$} {3:.0$} {4:.0$} ",
+                decimals, coord[0], coord[1], coord[2], coord[3]
+            ),
+            1 => println!("{1:.0$} ", decimals, coord[0]),
+            2 => println!("{1:.0$} {2:.0$} ", decimals, coord[0], coord[1]),
+            3 => println!(
+                "{1:.0$} {2:.0$} {3:.0$} ",
+                decimals, coord[0], coord[1], coord[2]
+            ),
+            _ => println!(
+                "{1:.0$} {2:.0$} {3:.0$} {4:.0$} ",
+                decimals, coord[0], coord[1], coord[2], coord[3]
+            ),
+        }
     }
     Ok(n)
 }