diff --git a/src/core/sensor/obsel.cpp b/src/core/sensor/obsel.cpp index f42df16314..b158e0d8e3 100644 --- a/src/core/sensor/obsel.cpp +++ b/src/core/sensor/obsel.cpp @@ -4,8 +4,12 @@ #include #include +#include +#include #include +#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& 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 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 unique_frequency_grids( + const std::span& obsels) { + std::vector 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 unique_poslos_grids( + const std::span& obsels) { + std::vector 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& obsels) { SensorSimulations out; + std::unordered_set 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 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 all_freq; // Collect all frequencies in the simulations std::vector 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::write(std::ostream& os, const ArrayOfSensorObsel& g, bofstream* pbofs, std::string_view name) { - const auto sen = collect_simulations(g); - - std::vector> 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> 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>(), - 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>(), - pbofs, - "poslos_pos"); - xml_write_to_stream(os, - g | stdv::transform([](auto& elem) { - return elem.weight_matrix(); - }) | stdr::to>(), - 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::read(std::istream& is, g.resize(0); - std::vector> f_grid; - std::vector> poslos; + std::vector f_grid_vals; + std::vector poslos_vals; std::vector f_grid_pos, poslos_pos; std::vector 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(val); + })}; + + const std::vector poslos{ + std::from_range, poslos_vals | stdv::transform([](auto& val) { + return std::make_shared(val); + })}; + const Size N = f_grid_pos.size(); if (N != poslos_pos.size() or N != mat.size()) { diff --git a/src/core/sensor/obsel.h b/src/core/sensor/obsel.h index 7f5d5d8c9d..4517ace946 100644 --- a/src/core/sensor/obsel.h +++ b/src/core/sensor/obsel.h @@ -10,6 +10,8 @@ #include #include +#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& other) const; @@ -119,6 +123,8 @@ class Obsel { [[nodiscard]] bool same_poslos( const std::shared_ptr& 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; -using SensorSimulations = std::unordered_map< - std::shared_ptr, - std::unordered_set>>; + +struct SensorSimulationsCache { + const AscendingGrid& freq_grid; + const SensorPosLosVector& poslos_grid; + const Size iposlos; +}; + +using SensorSimulations = std::vector; void unflatten(ArrayOfSensorObsel& sensor, const ConstVectorView& x, @@ -235,6 +246,12 @@ void make_exhaustive(std::span obsels); */ void make_exclusive(std::span obsels); +std::vector unique_frequency_grids( + const std::span& obsels); + +std::vector unique_poslos_grids( + const std::span& obsels); + SensorSimulations collect_simulations( const std::span& obsels); diff --git a/src/m_fwd.cc b/src/m_fwd.cc index 9f4a1048b6..3c36e71181 100644 --- a/src/m_fwd.cc +++ b/src/m_fwd.cc @@ -11,6 +11,7 @@ #include #include #include +#include #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 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); } } } diff --git a/src/m_rad.cc b/src/m_rad.cc index f4eed6a73f..d0be247c20 100644 --- a/src/m_rad.cc +++ b/src/m_rad.cc @@ -320,50 +320,7 @@ void measurement_vecFromSensor( for (auto &obsel : measurement_sensor) obsel.check(); const SensorSimulations simulations = collect_simulations(measurement_sensor); - - const auto flat_size = [](const SensorSimulations &simulations) { - Size size = 0; - for (const auto &[f_grid_ptr, poslos_set] : simulations) { - for (const auto &poslos_gs : poslos_set) { - size += poslos_gs->size(); - } - } - - return size; - }; - - const Size N = flat_size(simulations); - - struct unflatten_data { - const std::shared_ptr *f_grid_ptr{nullptr}; - const std::shared_ptr *poslos_ptr{nullptr}; - Size ip{std::numeric_limits::max()}; - - unflatten_data(const SensorSimulations &simulations, Size n) { - for (auto &[f_grid_ptr, poslos_set] : simulations) { - for (auto &poslos_gs : poslos_set) { - const Size np = poslos_gs->size(); - if (n < np) { - this->f_grid_ptr = &f_grid_ptr; - this->poslos_ptr = &poslos_gs; - this->ip = n; - goto end; - } - n -= np; - } - } - - end: - ARTS_USER_ERROR_IF( - f_grid_ptr == nullptr or poslos_ptr == nullptr, - "Failed to find f_grid_ptr and poslos_ptr for index {} in simulations", - n) - ARTS_USER_ERROR_IF(ip >= (**poslos_ptr).size(), - "Index {} out of bounds for poslos_ptr with size {}", - ip, - (**poslos_ptr).size()); - } - }; + const Size N = simulations.size(); std::string error{}; @@ -371,11 +328,9 @@ void measurement_vecFromSensor( #pragma omp parallel for schedule(dynamic) if (arts_omp_parallel(-1, N > 1)) for (Size i = 0; i < N; i++) { try { - const unflatten_data unflat(simulations, i); - - const Size ip = unflat.ip; - auto &f_grid_ptr = *unflat.f_grid_ptr; - auto &poslos = (**unflat.poslos_ptr)[ip]; + const Size ip = simulations[i].iposlos; + const auto &freq_grid = simulations[i].freq_grid; + const auto &poslos = simulations[i].poslos_grid[ip]; StokvecVector spectral_rad; StokvecMatrix spectral_rad_jac; @@ -385,7 +340,7 @@ void measurement_vecFromSensor( spectral_rad, spectral_rad_jac, ray_path, - *f_grid_ptr, + freq_grid, jac_targets, poslos.pos, poslos.los, @@ -396,12 +351,12 @@ void measurement_vecFromSensor( ARTS_USER_ERROR_IF(ray_path.empty(), "No ray path found"); spectral_rad_transform_operator( - spectral_rad, spectral_rad_jac, *f_grid_ptr, ray_path.front()); + spectral_rad, spectral_rad_jac, freq_grid, ray_path.front()); #pragma omp critical for (Size iv = 0; iv < measurement_sensor.size(); ++iv) { const SensorObsel &obsel = measurement_sensor[iv]; - if (obsel.same_freqs(f_grid_ptr)) { + if (obsel.same_freqs(freq_grid)) { measurement_vec[iv] += obsel.sumup(spectral_rad, ip); obsel.sumup(measurement_jac[iv], spectral_rad_jac, ip); diff --git a/src/python_interface/py_sensor.cpp b/src/python_interface/py_sensor.cpp index cb4c469b6f..a5b986e64b 100644 --- a/src/python_interface/py_sensor.cpp +++ b/src/python_interface/py_sensor.cpp @@ -270,25 +270,18 @@ Numeric, Vector, or Matrix a1.def( "unique_freq_grids", [](const ArrayOfSensorObsel& x) { - const SensorSimulations simuls = collect_simulations(x); - std::vector out; - out.reserve(simuls.size()); - for (auto& k : simuls | stdv::keys) out.push_back(*k); - return out; + return std::vector{std::from_range, + unique_frequency_grids(x) | + stdv::transform([](auto& ptr) { return *ptr; })}; }, "List of the unique frequency grids"); a1.def( "unique_poslos_grids", [](const ArrayOfSensorObsel& x) { - const SensorSimulations simuls = collect_simulations(x); - std::unordered_set> vecs; - for (auto& k : simuls | stdv::values) vecs.insert(k.begin(), k.end()); - - std::vector out; - out.reserve(vecs.size()); - for (auto& v : vecs) out.push_back(*v); - return out; + return std::vector{std::from_range, + unique_poslos_grids(x) | + stdv::transform([](auto& ptr) { return *ptr; })}; }, "List of the unique poslos grids");