Refactor sensor

src/core/sensor/obsel.cpp

@@ -4,8 +4,12 @@
 #include <debug.h>
 
 #include <algorithm>
+#include <memory>
+#include <ranges>
 #include <utility>
 
+#include "matpack_mdspan_helpers_grid_t.h"
+
 bool SensorKey::operator==(const SensorKey& other) const {
   return other.sensor_elem == sensor_elem and
          other.measurement_elem == measurement_elem and other.type == type;
@@ -179,6 +183,10 @@ Obsel::Obsel(const AscendingGrid& fs,
 
 bool Obsel::same_freqs(const Obsel& other) const { return f == other.f; }
 
+bool Obsel::same_freqs(const AscendingGrid& other) const {
+  return f.get() == &other;
+}
+
 bool Obsel::same_freqs(
     const std::shared_ptr<const AscendingGrid>& other) const {
   return f == other;
@@ -193,6 +201,10 @@ bool Obsel::same_poslos(
   return poslos == other;
 }
 
+bool Obsel::same_poslos(const PosLosVector& other) const {
+  return poslos.get() == &other;
+}
+
 void Obsel::set_f_grid_ptr(std::shared_ptr<const AscendingGrid> n) {
   ARTS_USER_ERROR_IF(not n, "Must exist");
   ARTS_USER_ERROR_IF(n->size() != f->size(), "Mismatching size");
@@ -371,12 +383,52 @@ Index Obsel::find(const AscendingGrid& freq_grid) const {
 }
 }  // namespace sensor
 
+std::vector<const AscendingGrid*> unique_frequency_grids(
+    const std::span<const SensorObsel>& obsels) {
+  std::vector<const AscendingGrid*> out;
+
+  for (const auto& obsel : obsels) {
+    const auto& f_grid = obsel.f_grid();
+
+    if (not stdr::contains(out, &f_grid)) {
+      out.push_back(&f_grid);
+    }
+  }
+
+  return out;
+}
+
+std::vector<const SensorPosLosVector*> unique_poslos_grids(
+    const std::span<const SensorObsel>& obsels) {
+  std::vector<const SensorPosLosVector*> out;
+
+  for (const auto& obsel : obsels) {
+    const auto& poslos_grid = obsel.poslos_grid();
+    if (not stdr::contains(out, &poslos_grid)) {
+      out.push_back(&poslos_grid);
+    }
+  }
+
+  return out;
+}
+
 SensorSimulations collect_simulations(
     const std::span<const SensorObsel>& obsels) {
   SensorSimulations out;
 
+  std::unordered_set<const AscendingGrid*> buf{};
+
   for (const auto& obsel : obsels) {
-    out[obsel.f_grid_ptr()].insert(obsel.poslos_grid_ptr());
+    auto& f_grid = obsel.f_grid();
+
+    if (buf.contains(&f_grid)) continue;
+    buf.insert(&f_grid);
+
+    auto& poslos = obsel.poslos_grid();
+
+    for (Size i = 0; i < poslos.size(); ++i) {
+      out.emplace_back(f_grid, poslos, i);
+    }
   }
 
   return out;
@@ -387,21 +439,14 @@ void make_exhaustive(std::span<SensorObsel> obsels) {
 
   // Early return on trivial cases
   if (simuls.size() == 0) return;
-  if (simuls.size() == 1 and simuls.begin()->second.size() <= 1) return;
+  if (simuls.size() == 1 and simuls.front().poslos_grid.size() <= 1) return;
 
   std::vector<Numeric> all_freq;  // Collect all frequencies in the simulations
   std::vector<SensorPosLos> all_geom;  // Collect all poslos in the simulations
 
-  for (auto& [freqs, geoms] : simuls) {
-    all_freq.insert(all_freq.end(), freqs->begin(), freqs->end());
-
-    for (auto& geom : geoms) {
-      for (auto& poslos : *geom) {
-        if (stdr::find(all_geom, poslos) == all_geom.end()) {
-          all_geom.push_back(poslos);
-        }
-      }
-    }
+  for (auto& [freqs, poslos, iposlos] : simuls) {
+    all_freq.insert(all_freq.end(), freqs.begin(), freqs.end());
+    all_geom.emplace_back(poslos[iposlos]);
   }
 
   stdr::sort(all_freq);
@@ -771,47 +816,47 @@ void xml_io_stream<ArrayOfSensorObsel>::write(std::ostream& os,
                                               const ArrayOfSensorObsel& g,
                                               bofstream* pbofs,
                                               std::string_view name) {
-  const auto sen = collect_simulations(g);
-
-  std::vector<std::shared_ptr<const SensorPosLosVector>> plos;
-  plos.reserve(sen.size());
-  for (const auto& i : sen | stdv::values | stdv::join) {
-    if (not stdr::contains(plos, i)) {
-      plos.push_back(i);
-    }
-  }
+  const auto freqs_ptrs = unique_frequency_grids(g);
+  const auto plos_ptrs  = unique_poslos_grids(g);
+
+  const std::vector freqs{
+      std::from_range,
+      freqs_ptrs | stdv::transform([](auto& ptr) { return *ptr; })};
+
+  const std::vector plos{
+      std::from_range,
+      plos_ptrs | stdv::transform([](auto& ptr) { return *ptr; })};
+
+  const std::vector freqs_pos{
+      std::from_range,
+      g | stdv::transform([&freqs_ptrs](const SensorObsel& elem) -> Index {
+        for (Size i = 0; i < freqs_ptrs.size(); ++i) {
+          if (elem.same_freqs(*freqs_ptrs[i])) return i;
+        }
+        return -999;  // Should not happen
+      })};
+
+  const std::vector plos_pos{
+      std::from_range,
+      g | stdv::transform([&plos_ptrs](const SensorObsel& elem) -> Index {
+        for (Size i = 0; i < plos_ptrs.size(); ++i) {
+          if (elem.same_poslos(*plos_ptrs[i])) return i;
+        }
+        return -999;  // Should not happen
+      })};
 
-  std::vector<std::shared_ptr<const AscendingGrid>> freqs;
-  freqs.reserve(sen.size());
-  for (const auto& i : sen | stdv::keys) freqs.push_back(i);
+  const std::vector weights{
+      std::from_range,
+      g | stdv::transform([](auto& elem) { return elem.weight_matrix(); })};
 
   XMLTag tag(type_name, "name", name);
   tag.write_to_stream(os);
 
   xml_write_to_stream(os, freqs, pbofs, "f_grid");
   xml_write_to_stream(os, plos, pbofs, "poslos");
-  xml_write_to_stream(os,
-                      g | stdv::transform([&freqs](auto& elem) {
-                        return std::distance(
-                            freqs.begin(),
-                            stdr::find(freqs, elem.f_grid_ptr()));
-                      }) | stdr::to<std::vector<Index>>(),
-                      pbofs,
-                      "f_grid_pos");
-  xml_write_to_stream(os,
-                      g | stdv::transform([&plos](auto& elem) {
-                        return std::distance(
-                            plos.begin(),
-                            stdr::find(plos, elem.poslos_grid_ptr()));
-                      }) | stdr::to<std::vector<Index>>(),
-                      pbofs,
-                      "poslos_pos");
-  xml_write_to_stream(os,
-                      g | stdv::transform([](auto& elem) {
-                        return elem.weight_matrix();
-                      }) | stdr::to<std::vector<sensor::SparseStokvecMatrix>>(),
-                      pbofs,
-                      "Weights");
+  xml_write_to_stream(os, freqs_pos, pbofs, "f_grid_pos");
+  xml_write_to_stream(os, plos_pos, pbofs, "poslos_pos");
+  xml_write_to_stream(os, weights, pbofs, "Weights");
 
   tag.write_to_end_stream(os);
 }
@@ -825,17 +870,27 @@ void xml_io_stream<ArrayOfSensorObsel>::read(std::istream& is,
 
   g.resize(0);
 
-  std::vector<std::shared_ptr<AscendingGrid>> f_grid;
-  std::vector<std::shared_ptr<SensorPosLosVector>> poslos;
+  std::vector<AscendingGrid> f_grid_vals;
+  std::vector<SensorPosLosVector> poslos_vals;
   std::vector<Index> f_grid_pos, poslos_pos;
   std::vector<sensor::SparseStokvecMatrix> mat;
 
-  xml_read_from_stream(is, f_grid, pbifs);
-  xml_read_from_stream(is, poslos, pbifs);
+  xml_read_from_stream(is, f_grid_vals, pbifs);
+  xml_read_from_stream(is, poslos_vals, pbifs);
   xml_read_from_stream(is, f_grid_pos, pbifs);
   xml_read_from_stream(is, poslos_pos, pbifs);
   xml_read_from_stream(is, mat, pbifs);
 
+  const std::vector f_grid{std::from_range,
+                           f_grid_vals | stdv::transform([](auto& val) {
+                             return std::make_shared<const AscendingGrid>(val);
+                           })};
+
+  const std::vector poslos{
+      std::from_range, poslos_vals | stdv::transform([](auto& val) {
+                         return std::make_shared<const SensorPosLosVector>(val);
+                       })};
+
   const Size N = f_grid_pos.size();
 
   if (N != poslos_pos.size() or N != mat.size()) {

src/core/sensor/obsel.h

@@ -10,6 +10,8 @@
 #include <unordered_map>
 #include <unordered_set>
 
+#include "matpack_mdspan_helpers_grid_t.h"
+
 namespace sensor {
 struct PosLos {
   Vector3 pos;
@@ -111,6 +113,8 @@ class Obsel {
 
   [[nodiscard]] bool same_freqs(const Obsel& other) const;
 
+  [[nodiscard]] bool same_freqs(const AscendingGrid& other) const;
+
   [[nodiscard]] bool same_freqs(
       const std::shared_ptr<const AscendingGrid>& other) const;
 
@@ -119,6 +123,8 @@ class Obsel {
   [[nodiscard]] bool same_poslos(
       const std::shared_ptr<const PosLosVector>& other) const;
 
+  [[nodiscard]] bool same_poslos(const PosLosVector& other) const;
+
   [[nodiscard]] const auto& f_grid_ptr() const { return f; }
   [[nodiscard]] const auto& poslos_grid_ptr() const { return poslos; }
 
@@ -204,9 +210,14 @@ using SensorPosLos       = sensor::PosLos;
 using SensorPosLosVector = sensor::PosLosVector;
 using SensorObsel        = sensor::Obsel;
 using ArrayOfSensorObsel = Array<SensorObsel>;
-using SensorSimulations  = std::unordered_map<
-     std::shared_ptr<const AscendingGrid>,
-     std::unordered_set<std::shared_ptr<const SensorPosLosVector>>>;
+
+struct SensorSimulationsCache {
+  const AscendingGrid& freq_grid;
+  const SensorPosLosVector& poslos_grid;
+  const Size iposlos;
+};
+
+using SensorSimulations = std::vector<SensorSimulationsCache>;
 
 void unflatten(ArrayOfSensorObsel& sensor,
                const ConstVectorView& x,
@@ -235,6 +246,12 @@ void make_exhaustive(std::span<SensorObsel> obsels);
  */
 void make_exclusive(std::span<SensorObsel> obsels);
 
+std::vector<const AscendingGrid*> unique_frequency_grids(
+    const std::span<const SensorObsel>& obsels);
+
+std::vector<const SensorPosLosVector*> unique_poslos_grids(
+    const std::span<const SensorObsel>& obsels);
+
 SensorSimulations collect_simulations(
     const std::span<const SensorObsel>& obsels);
 

src/m_fwd.cc

@@ -11,6 +11,7 @@
 #include <algorithm>
 #include <exception>
 #include <memory>
+#include <ranges>
 
 #include "workspace_class.h"
 
@@ -289,31 +290,26 @@ void measurement_vecFromOperatorPath(
   const SensorSimulations simulations =
       collect_simulations(measurement_vec_sensor);
 
-  for (auto& [f_grid_ptr, poslos_set] : simulations) {
-    for (auto& poslos_gs : poslos_set) {
-      for (Size ip = 0; ip < poslos_gs->size(); ++ip) {
-        ArrayOfPropagationPathPoint ray_path;
-        std::vector<fwd::path> path;
-        StokvecVector spectral_rad;
-
-        const SensorPosLos& poslos = (*poslos_gs)[ip];
-
-        ray_path_observer_agendaExecute(
-            ws, ray_path, poslos.pos, poslos.los, ray_path_observer_agenda);
-        spectral_rad_operator.from_path(path, ray_path);
-        std::transform(f_grid_ptr->begin(),
-                       f_grid_ptr->end(),
-                       spectral_rad.begin(),
-                       [&path, &spectral_rad_operator](Numeric f) {
-                         return spectral_rad_operator(f, path);
-                       });
-
-        for (Size iv = 0; iv < measurement_vec_sensor.size(); ++iv) {
-          const SensorObsel& obsel = measurement_vec_sensor[iv];
-          if (obsel.same_freqs(f_grid_ptr)) {
-            measurement_vec[iv] += obsel.sumup(spectral_rad, ip);
-          }
-        }
+  for (auto& [freq_grid, poslos_grid, iposlos] : simulations) {
+    const auto& poslos = poslos_grid[iposlos];
+
+    ArrayOfPropagationPathPoint ray_path;
+
+    ray_path_observer_agendaExecute(
+        ws, ray_path, poslos.pos, poslos.los, ray_path_observer_agenda);
+
+    const StokvecVector spectral_rad{
+        std::from_range,
+        freq_grid |
+            stdv::transform([path = spectral_rad_operator.from_path(ray_path),
+                             &spectral_rad_operator](Numeric f) {
+              return spectral_rad_operator(f, path);
+            })};
+
+    for (Size iv = 0; iv < measurement_vec_sensor.size(); ++iv) {
+      const SensorObsel& obsel = measurement_vec_sensor[iv];
+      if (obsel.same_freqs(freq_grid)) {
+        measurement_vec[iv] += obsel.sumup(spectral_rad, iposlos);
       }
     }
   }