From 952ff64350d60026e22813f41569e927527d181d Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Mon, 18 Nov 2024 13:09:27 +0100 Subject: [PATCH 1/7] Rearrange some cmake --- src/core/CMakeLists.txt | 3 +- src/core/absorption/CMakeLists.txt | 2 +- src/core/atm/CMakeLists.txt | 2 +- src/core/fwd/CMakeLists.txt | 2 +- src/core/jacobian/CMakeLists.txt | 6 ++ src/core/{ => jacobian}/jacobian.cc | 15 ++-- src/core/{ => jacobian}/jacobian.h | 83 +++++++++++++--------- src/core/lbl/CMakeLists.txt | 4 +- src/core/options/arts_options.cc | 25 +++++++ src/core/path/CMakeLists.txt | 2 +- src/core/physics/CMakeLists.txt | 2 +- src/core/quantum/CMakeLists.txt | 2 +- src/core/rtepack/CMakeLists.txt | 2 +- src/core/sensor/CMakeLists.txt | 2 +- src/core/sensor/obsel.cpp | 104 ++++++++++++++++++++++++++++ src/core/sensor/obsel.h | 49 +++++++++++-- src/m_jactargets.cc | 86 ++++++++++++++++------- 17 files changed, 310 insertions(+), 81 deletions(-) create mode 100644 src/core/jacobian/CMakeLists.txt rename src/core/{ => jacobian}/jacobian.cc (97%) rename src/core/{ => jacobian}/jacobian.h (89%) diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 3b61ec740d..0bf2877ecd 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -21,6 +21,7 @@ add_subdirectory(matpack) add_subdirectory(rtepack) add_subdirectory(predefined) add_subdirectory(operators) +add_subdirectory(jacobian) add_subdirectory(fwd) add_subdirectory(lbl) add_subdirectory(lookup) @@ -48,7 +49,6 @@ add_library(artscore STATIC geodetic.cc igrf13.cc interpolation.cc - jacobian.cc jpl_species.cc lineshapemodel.cc mc_antenna.cc @@ -93,6 +93,7 @@ target_link_libraries(artscore PUBLIC util disort-cpp legendre + jacobian ) target_include_directories(artscore PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") diff --git a/src/core/absorption/CMakeLists.txt b/src/core/absorption/CMakeLists.txt index ce70fed72c..1997e570f5 100644 --- a/src/core/absorption/CMakeLists.txt +++ b/src/core/absorption/CMakeLists.txt @@ -5,4 +5,4 @@ add_library(absorption STATIC predefined_absorption_models.cc ) target_link_libraries(absorption PUBLIC path predef physics lbl) -target_include_directories(absorption PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(absorption PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/atm/CMakeLists.txt b/src/core/atm/CMakeLists.txt index 0ea76716b7..9a520923b9 100644 --- a/src/core/atm/CMakeLists.txt +++ b/src/core/atm/CMakeLists.txt @@ -1,4 +1,4 @@ add_library(atm STATIC atm.cpp) target_link_libraries(atm PUBLIC matpack operators quantum species physics) -target_include_directories(atm PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) \ No newline at end of file +target_include_directories(atm PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) \ No newline at end of file diff --git a/src/core/fwd/CMakeLists.txt b/src/core/fwd/CMakeLists.txt index fbb1f9c9bf..cb09dc445d 100644 --- a/src/core/fwd/CMakeLists.txt +++ b/src/core/fwd/CMakeLists.txt @@ -7,4 +7,4 @@ add_library(fwd STATIC fwd_spectral_radiance.cpp ) target_link_libraries(fwd PUBLIC path absorption) -target_include_directories(fwd PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(fwd PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/jacobian/CMakeLists.txt b/src/core/jacobian/CMakeLists.txt new file mode 100644 index 0000000000..7987df4e2f --- /dev/null +++ b/src/core/jacobian/CMakeLists.txt @@ -0,0 +1,6 @@ +add_library(jacobian STATIC +jacobian.cc +) + +target_link_libraries(jacobian PUBLIC atm sensor scattering lbl) +target_include_directories(jacobian PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/jacobian.cc b/src/core/jacobian/jacobian.cc similarity index 97% rename from src/core/jacobian.cc rename to src/core/jacobian/jacobian.cc index 23bc435098..5deb87ea85 100644 --- a/src/core/jacobian.cc +++ b/src/core/jacobian/jacobian.cc @@ -166,7 +166,7 @@ void default_sensor_x_set(ExhaustiveVectorView x, using enum SensorKeyType; switch (key.model) { - case SensorModelType::None: { + case SensorJacobianModelType::None: { switch (key.type) { case Frequency: ARTS_USER_ERROR_IF(x.size() != v.f_grid().size(), @@ -243,7 +243,7 @@ void default_sensor_x_set(ExhaustiveVectorView x, } break; } } break; - case SensorModelType::PolynomialOffset: { + case SensorJacobianModelType::PolynomialOffset: { const Vector& o = key.original_grid; switch (key.type) { case Frequency: polyfit(x, o, rem_frq(v, o)); break; @@ -287,7 +287,7 @@ void set_frq(const SensorObsel& v, const auto xs = std::make_shared( x.begin(), x.end(), [](auto& x) { return x; }); - // Must copy, as we may change the shared_ptr later, regrdless of clang-tidy + // Must copy, as we may change the shared_ptr later const auto fs = v.f_grid_ptr(); for (auto& elem : sensor) { @@ -323,7 +323,7 @@ void set_poslos(const SensorObsel& v, const auto xs = std::make_shared(std::move(xsv)); - // Must copy, as we may change the shared_ptr later, regrdless of clang-tidy + // Must copy, as we may change the shared_ptr later const auto ps = v.poslos_grid_ptr(); for (auto& elem : sensor) { @@ -370,7 +370,7 @@ void default_x_sensor_set(ArrayOfSensorObsel& sensor, using enum SensorKeyType; switch (key.model) { - case SensorModelType::None: { + case SensorJacobianModelType::None: { switch (key.type) { case Frequency: set_frq(v, sensor, x); break; case PointingZenith: set_zag(v, sensor, x); break; @@ -401,7 +401,7 @@ void default_x_sensor_set(ArrayOfSensorObsel& sensor, } break; } } break; - case SensorModelType::PolynomialOffset: { + case SensorJacobianModelType::PolynomialOffset: { auto& o = key.original_grid; Vector r = polynomial_offset_evaluate(x, o); @@ -481,7 +481,8 @@ std::vector& Targets::line() { return target(); } void Targets::finalize(const AtmField& atmospheric_field, const SurfaceField& surface_field, - const AbsorptionBands&) { + const AbsorptionBands&, + const ArrayOfSensorObsel& measurement_sensor) { zero_out_x(); const Size natm = atm().size(); diff --git a/src/core/jacobian.h b/src/core/jacobian/jacobian.h similarity index 89% rename from src/core/jacobian.h rename to src/core/jacobian/jacobian.h index 81328091ca..58c09eb879 100644 --- a/src/core/jacobian.h +++ b/src/core/jacobian/jacobian.h @@ -3,7 +3,7 @@ #include #include #include -#include +#include #include #include @@ -16,31 +16,6 @@ #include "operators.h" -enum class SensorKeyType { - Frequency, - PointingZenith, - PointingAzimuth, - PointingAltitude, - PointingLatitude, - PointingLongitude, - Weights, -}; - -enum class SensorModelType { - None, - PolynomialOffset, -}; - -struct SensorKey { - SensorKeyType type; - - Index elem; - - SensorModelType model{SensorModelType::None}; - - Vector original_grid{}; -}; - using AtmTargetSetState = CustomOperator(), - [&](auto& a) { + [xsize, t_size](auto& a) { ARTS_USER_ERROR_IF((a.x_start + a.x_size) > xsize, "The target {}" " is out of bounds of the x-vector. (xsize: {})", @@ -300,29 +275,33 @@ struct targets_t { } }; -struct Targets final : targets_t { +struct Targets final : targets_t { [[nodiscard]] const std::vector& atm() const; [[nodiscard]] const std::vector& surf() const; [[nodiscard]] const std::vector& line() const; + [[nodiscard]] const std::vector& sensor() const; [[nodiscard]] std::vector& atm(); [[nodiscard]] std::vector& surf(); [[nodiscard]] std::vector& line(); + [[nodiscard]] std::vector& sensor(); //! Sets the sizes and x-positions of the targets. void finalize(const AtmField& atmospheric_field, const SurfaceField& surface_field, - const AbsorptionBands& absorption_bands); + const AbsorptionBands& absorption_bands, + const ArrayOfSensorObsel& measurement_sensor); }; struct TargetType { - using variant_t = std::variant; + using variant_t = std::variant; variant_t target; - template + template [[nodiscard]] constexpr auto apply(const AtmKeyValFunc& ifatm, const SurfaceKeyValFunc& ifsurf, - const LblLineKeyFunc& ifline) const { + const LblLineKeyFunc& ifline, + const SensorKeyFunc& ifsensor) const { return std::visit( [&](const auto& t) { using T = std::decay_t; @@ -332,6 +311,8 @@ struct TargetType { return ifsurf(t); } else if constexpr (std::same_as) { return ifline(t); + } else if constexpr (std::same_as) { + return ifsensor(t); } }, target); @@ -342,7 +323,8 @@ struct TargetType { [[nodiscard]] std::string type() const { return apply([](auto&) { return "AtmKeyVal"s; }, [](auto&) { return "SurfaceKeyVal"s; }, - [](auto&) { return "LblLineKey"s; }); + [](auto&) { return "LblLineKey"s; }, + [](auto&) { return "SensorKey"s; }); } }; } // namespace Jacobian @@ -390,6 +372,9 @@ struct std::formatter { }, [this, &ctx](const LblLineKey& key) { tags.format(ctx, "LblLineKey::"sv, key); + }, + [this, &ctx](const SensorKey& key) { + tags.format(ctx, "SensorKey::"sv, key); }); return ctx.out(); @@ -492,6 +477,38 @@ struct std::formatter { } }; +template <> +struct std::formatter { + format_tags tags; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const Jacobian::SensorTarget& v, + FmtContext& ctx) const { + const std::string_view sep = tags.sep(); + return tags.format(ctx, + v.type, + ": "sv, + v.d, + sep, + "target_pos: "sv, + v.target_pos, + sep, + "x_start: "sv, + v.x_start, + sep, + "x_size: "sv, + v.x_size); + } +}; + template <> struct std::formatter { format_tags tags; diff --git a/src/core/lbl/CMakeLists.txt b/src/core/lbl/CMakeLists.txt index 1ba6d26019..fcf852ea0b 100644 --- a/src/core/lbl/CMakeLists.txt +++ b/src/core/lbl/CMakeLists.txt @@ -15,5 +15,5 @@ add_library(lbl STATIC lbl_zeeman.cpp ) -target_link_libraries(lbl PUBLIC matpack atm Faddeeva predef quantum surface wigner util arts_options) -target_include_directories(lbl PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_link_libraries(lbl PUBLIC matpack atm jacobian Faddeeva predef quantum surface wigner util arts_options) +target_include_directories(lbl PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/options/arts_options.cc b/src/core/options/arts_options.cc index 52722e67b4..fececec53c 100644 --- a/src/core/options/arts_options.cc +++ b/src/core/options/arts_options.cc @@ -14,6 +14,31 @@ using Value = std::vector; std::vector internal_options_create() { std::vector opts; + opts.emplace_back(EnumeratedOption{ + .name = "SensorKeyType", + .desc = + R"(A key for identifying a sensor property +)", + .values_and_desc = { + Value{"Frequency", "Frequency of the sensor"}, + Value{"PointingZenith", "Pointing Zenith of the sensor"}, + Value{"PointingAzimuth", "Pointing Azimuth of the sensor"}, + Value{"PointingAltitude", "Pointing Altitude of the sensor"}, + Value{"PointingLatitude", "Pointing Latitude of the sensor"}, + Value{"PointingLongitude", "Pointing Longitude of the sensor"}, + Value{"Weights", "Weights of the sensor"}, + }}); + + opts.emplace_back(EnumeratedOption{ + .name = "SensorJacobianModelType", + .desc = + R"(How to model the sensor Jacobian model target. +)", + .values_and_desc = { + Value{"None", "No model, work purely on sensor data"}, + Value{"PolynomialOffset", "The sensor Jacobian is modeled as a polynomial offset"}, + }}); + opts.emplace_back(EnumeratedOption{ .name = "HitranLineStrengthOption", .desc = diff --git a/src/core/path/CMakeLists.txt b/src/core/path/CMakeLists.txt index cca53a99e2..cd7a516207 100644 --- a/src/core/path/CMakeLists.txt +++ b/src/core/path/CMakeLists.txt @@ -1,4 +1,4 @@ add_library(path STATIC path_point.cpp atm_path.cpp) target_link_libraries(path PUBLIC matpack atm surface) -target_include_directories(path PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(path PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/physics/CMakeLists.txt b/src/core/physics/CMakeLists.txt index b8330c3066..e8cb2e6684 100644 --- a/src/core/physics/CMakeLists.txt +++ b/src/core/physics/CMakeLists.txt @@ -3,4 +3,4 @@ physics_funcs.cc wigner_functions.cc ) target_link_libraries(physics PUBLIC matpack wigner) -target_include_directories(physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/quantum/CMakeLists.txt b/src/core/quantum/CMakeLists.txt index 493b057db5..9842467d0c 100644 --- a/src/core/quantum/CMakeLists.txt +++ b/src/core/quantum/CMakeLists.txt @@ -4,4 +4,4 @@ add_library(quantum STATIC quantum_term_symbol.cc ) target_link_libraries(quantum PUBLIC species matpack binio) -target_include_directories(quantum PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(quantum PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/core/rtepack/CMakeLists.txt b/src/core/rtepack/CMakeLists.txt index 863c82d615..fd796446c6 100644 --- a/src/core/rtepack/CMakeLists.txt +++ b/src/core/rtepack/CMakeLists.txt @@ -9,5 +9,5 @@ add_library(rtepack STATIC rtepack_transmission.cc rtepack.cc ) -target_include_directories(rtepack PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/..) +target_include_directories(rtepack PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) target_link_libraries(rtepack PUBLIC matpack physics) diff --git a/src/core/sensor/CMakeLists.txt b/src/core/sensor/CMakeLists.txt index b2fe8837b6..b99c25e1e6 100644 --- a/src/core/sensor/CMakeLists.txt +++ b/src/core/sensor/CMakeLists.txt @@ -1,3 +1,3 @@ add_library(sensor STATIC obsel.cpp) +target_link_libraries(sensor PUBLIC matpack rtepack arts_enum_options) target_include_directories(sensor PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) -target_link_libraries(sensor PUBLIC matpack rtepack) diff --git a/src/core/sensor/obsel.cpp b/src/core/sensor/obsel.cpp index 417e2d6c36..76c6ba94ab 100644 --- a/src/core/sensor/obsel.cpp +++ b/src/core/sensor/obsel.cpp @@ -1,10 +1,15 @@ #include "obsel.h" #include +#include #include "compare.h" #include "debug.h" +bool SensorKey::operator==(const SensorKey& other) const { + return other.elem == elem and other.model == model and other.type == type; +} + namespace sensor { std::ostream& operator<<(std::ostream& os, const PosLos& obsel) { return os << "[pos: [" << obsel.pos << "], los: [" << obsel.los << "]]"; @@ -70,6 +75,105 @@ void Obsel::sumup(VectorView out, const StokvecMatrixView& j, Index ip) const { out[ij] += std::transform_reduce(jac.begin(), jac.end(), ws.begin(), 0.0); } } + +Size Obsel::flat_size(const SensorKeyType& key) const { + switch (key) { + using enum SensorKeyType; + case Frequency: return f->size(); + case PointingZenith: return poslos->size(); + case PointingAzimuth: return poslos->size(); + case PointingAltitude: return poslos->size(); + case PointingLatitude: return poslos->size(); + case PointingLongitude: return poslos->size(); + case Weights: return f->size() * poslos->size() * 4; + } + + std::unreachable(); +} + +void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { + switch (key) { + using enum SensorKeyType; + case Frequency: + ARTS_USER_ERROR_IF(x.size() != f_grid().size(), + "Bad size. x.size(): {}, f_grid().size(): {}", + x.size(), + f_grid().size()) + x = f_grid(); + break; + case PointingZenith: + ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), + "Bad size. x.size(): {}, poslos_grid().size(): {}", + x.size(), + poslos_grid().size()) + std::transform(poslos_grid().begin(), + poslos_grid().end(), + x.begin(), + [](auto& poslos) { return poslos.los[0]; }); + break; + case PointingAzimuth: + ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), + "Bad size. x.size(): {}, poslos_grid().size(): {}", + x.size(), + poslos_grid().size()) + std::transform(poslos_grid().begin(), + poslos_grid().end(), + x.begin(), + [](auto& poslos) { return poslos.los[1]; }); + break; + case PointingAltitude: + ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), + "Bad size. x.size(): {}, poslos_grid().size(): {}", + x.size(), + poslos_grid().size()) + std::transform(poslos_grid().begin(), + poslos_grid().end(), + x.begin(), + [](auto& poslos) { return poslos.pos[0]; }); + break; + case PointingLatitude: + ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), + "Bad size. x.size(): {}, poslos_grid().size(): {}", + x.size(), + poslos_grid().size()) + std::transform(poslos_grid().begin(), + poslos_grid().end(), + x.begin(), + [](auto& poslos) { return poslos.pos[1]; }); + break; + case PointingLongitude: + ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), + "Bad size. x.size(): {}, poslos_grid().size(): {}", + x.size(), + poslos_grid().size()) + std::transform(poslos_grid().begin(), + poslos_grid().end(), + x.begin(), + [](auto& poslos) { return poslos.pos[2]; }); + break; + case Weights: { + ARTS_USER_ERROR_IF(w.size() * 4 != x.size(), + "Bad size. x.size(): {} for w.size()*4: {}", + x.size(), + w.size() * 4) + auto X = x.reshape_as(w.nrows(), w.ncols(), 4); + for (Index i = 0; i < X.npages(); ++i) { + for (Index j = 0; j < X.nrows(); ++j) { + X(i, j, 0) = w(i, j).I(); + X(i, j, 1) = w(i, j).Q(); + X(i, j, 2) = w(i, j).U(); + X(i, j, 3) = w(i, j).V(); + } + } + } break; + } +} + +Vector Obsel::flat(const SensorKeyType& key) const { + Vector out(flat_size(key)); + flat(out, key); + return out; +} } // namespace sensor SensorSimulations collect_simulations(const ArrayOfSensorObsel& obsels) { diff --git a/src/core/sensor/obsel.h b/src/core/sensor/obsel.h index 440d89a3c2..805bd5e640 100644 --- a/src/core/sensor/obsel.h +++ b/src/core/sensor/obsel.h @@ -9,10 +9,8 @@ #include #include -#include "debug.h" -#include "format_tags.h" -#include "matpack_constexpr.h" -#include "sorted_grid.h" +#include +#include namespace sensor { struct PosLos { @@ -136,11 +134,54 @@ class Obsel { [[nodiscard]] Numeric sumup(const StokvecVectorView& i, Index ip) const; void sumup(VectorView out, const StokvecMatrixView& j, Index ip) const; + + [[nodiscard]] Size flat_size(const SensorKeyType& key) const; + void flat(ExhaustiveVectorView x, const SensorKeyType& key) const; + [[nodiscard]] Vector flat(const SensorKeyType& key) const; }; std::ostream& operator<<(std::ostream& os, const Array& obsel); } // namespace sensor +struct SensorKey { + SensorKeyType type; + + Index elem; + + SensorJacobianModelType model; + + Index polyorder{-1}; + + Vector original_grid{}; + + bool operator==(const SensorKey& other) const; +}; + +template<> +struct std::hash { + std::size_t operator()(const SensorKey& g) const { + return std::hash{}(g.type) ^ (std::hash{}(g.model) << 10) ^ (std::hash{}(g.elem) << 20); + } +}; + +template <> +struct std::formatter { + format_tags tags{}; + + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } + + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } + + template + FmtContext::iterator format(const SensorKey& v, FmtContext& ctx) const { + return tags.format(ctx, v.type, tags.sep(), v.elem, tags.sep(), v.model); + } +}; + using SensorPosLos = sensor::PosLos; using SensorPosLosVector = sensor::PosLosVector; using SensorObsel = sensor::Obsel; diff --git a/src/m_jactargets.cc b/src/m_jactargets.cc index 2dc7ddd6c1..67cdfade33 100644 --- a/src/m_jactargets.cc +++ b/src/m_jactargets.cc @@ -1,5 +1,5 @@ -#include #include +#include #include #include @@ -19,67 +19,61 @@ void jacobian_targetsFinalize(JacobianTargets& jacobian_targets, const AtmField& atmospheric_field, const SurfaceField& surface_field, const AbsorptionBands& absorption_bands) { - jacobian_targets.finalize(atmospheric_field, surface_field, absorption_bands); + jacobian_targets.finalize( + atmospheric_field, surface_field, absorption_bands, {}); } void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, const SurfaceKey& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.surf().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, const SurfaceTypeTag& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.surf().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, const SurfacePropertyTag& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.surf().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const AtmKey& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const SpeciesEnum& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const SpeciesIsotope& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const QuantumIdentifier& key, const Numeric& d) { - jacobian_targets.target().emplace_back( - key, d, jacobian_targets.target_count()); + jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); } void jacobian_targetsAddTemperature(JacobianTargets& jacobian_targets, const Numeric& d) { - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::t, d, jacobian_targets.target_count()); } void jacobian_targetsAddPressure(JacobianTargets& jacobian_targets, const Numeric& d) { - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::p, d, jacobian_targets.target_count()); } @@ -89,15 +83,15 @@ void jacobian_targetsAddMagneticField(JacobianTargets& jacobian_targets, using enum FieldComponent; switch (to(component)) { case u: - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::mag_u, d, jacobian_targets.target_count()); break; case v: - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::mag_v, d, jacobian_targets.target_count()); break; case w: - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::mag_w, d, jacobian_targets.target_count()); break; } @@ -109,15 +103,15 @@ void jacobian_targetsAddWindField(JacobianTargets& jacobian_targets, using enum FieldComponent; switch (to(component)) { case u: - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::wind_u, d, jacobian_targets.target_count()); break; case v: - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::wind_v, d, jacobian_targets.target_count()); break; case w: - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( AtmKey::wind_w, d, jacobian_targets.target_count()); break; } @@ -126,7 +120,7 @@ void jacobian_targetsAddWindField(JacobianTargets& jacobian_targets, void jacobian_targetsAddSpeciesVMR(JacobianTargets& jacobian_targets, const SpeciesEnum& species, const Numeric& d) { - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( species, d, jacobian_targets.target_count()); } @@ -139,6 +133,46 @@ void jacobian_targetsAddSpeciesIsotopologueRatio( "Unknown isotopologue: \"{}\"", species.FullName()); - jacobian_targets.target().emplace_back( + jacobian_targets.atm().emplace_back( species, d, jacobian_targets.target_count()); } + +void jacobian_targetsAddSensor(JacobianTargets& jacobian_targets, + const ArrayOfSensorObsel& measurement_sensor, + const SensorKeyType& key, + const Numeric& d, + const Index& elem) { + ARTS_USER_ERROR_IF(measurement_sensor.size() <= static_cast(elem), + "Sensor element out of bounds: {}", + elem) + + jacobian_targets.sensor().push_back({ + .type = {.type = key, .elem = elem, .model = SensorJacobianModelType::None}, + .d = d, + .target_pos = jacobian_targets.target_count(), + }); +} + +void jacobian_targetsAddSensorPolyFit( + JacobianTargets& jacobian_targets, + const ArrayOfSensorObsel& measurement_sensor, + const SensorKeyType& key, + const Numeric& d, + const Index& elem, + const Index& polyorder) { + ARTS_USER_ERROR_IF(measurement_sensor.size() <= static_cast(elem), + "Sensor element out of bounds: {}", + elem) + ARTS_USER_ERROR_IF( + polyorder < 0, "Polyorder must be non-negative: {}", polyorder) + + Jacobian::SensorTarget target{ + .type = {.type = key, + .elem = elem, + .model = SensorJacobianModelType::PolynomialOffset, + .polyorder = polyorder}, + // .original_grid = sensor_grid(measurement_sensor[elem], key)}, + .d = d, + .target_pos = jacobian_targets.target_count(), + }; +} From 6aa5a594935def961d6ef2eef74c84d87441ca8c Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Tue, 19 Nov 2024 17:28:56 +0100 Subject: [PATCH 2/7] Simplified sensor jacobian setups --- .../2.zeeman-sensor.py | 10 +- .../3-disort/1.clearsky-radiance.py | 3 +- .../3-disort/2.clearsky-flux.py | 2 +- python/src/pyarts/plots/ppvar_atm.py | 18 +- python/test/workspace/test_variables.py | 6 + src/CMakeLists.txt | 1 + src/core/jacobian/jacobian.cc | 269 ++++-------------- src/core/options/arts_options.cc | 85 +++--- src/core/sensor/obsel.cpp | 158 ++++++++-- src/core/sensor/obsel.h | 19 +- src/m_covmat.cc | 53 +++- src/m_jactargets.cc | 19 +- src/m_rad.cc | 229 ++++++++++++--- src/python_interface/gen_auto_py_options.cpp | 4 +- src/python_interface/py_global.cpp | 26 +- src/python_interface/py_lbl.cpp | 2 +- src/python_interface/py_path.cpp | 4 +- src/python_interface/py_quantum.cpp | 2 +- src/python_interface/py_species.cpp | 2 +- src/python_interface/py_spectroscopy.cpp | 12 +- src/workspace_agenda_creator.cpp | 11 +- src/workspace_agendas.cpp | 10 +- src/workspace_groups.cpp | 54 +--- src/workspace_meta_methods.cpp | 33 ++- src/workspace_methods.cpp | 64 ++++- src/workspace_variables.cpp | 12 +- src/xml_io_compound_types.cc | 41 +++ tests/core/jac/full_arts_emission.py | 1 + tests/core/jac/full_optimal_estimation.py | 15 +- .../jac/full_optimal_estimation_multiparam.py | 15 +- tests/core/jac/model_state_atm.py | 2 + tests/core/sensor/jac.py | 4 + tests/core/wind/propmat_jac.py | 2 +- .../wind/spectral_radiance_jacobian_wind.py | 10 +- tests/core/zeeman/propmat_jac.py | 2 +- ...ectral_radiance_jacobian_magnetic_field.py | 11 +- 36 files changed, 707 insertions(+), 504 deletions(-) create mode 100644 tests/core/sensor/jac.py diff --git a/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py b/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py index e681547dbd..c1e793da22 100644 --- a/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py +++ b/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py @@ -35,7 +35,7 @@ ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground") ws.spectral_radiance_surface_agendaSet(option="Blackbody") ws.ray_path_observer_agendaSet(option="Geometric") -ws.spectral_radiance_observer_agendaSet(option="EmissionUnits") +ws.spectral_radiance_observer_agendaSet(option="Emission") # %% Set up a sensor with Gaussian channel widths on individual frequency ranges @@ -45,13 +45,11 @@ # %% Core calculations -result = pyarts.arts.Vector() -result_jac = pyarts.arts.Matrix() -ws.measurement_vectorFromSensor(result, result_jac) +ws.measurement_vectorFromSensor() # %% Show results -plt.plot((ws.frequency_grid - line_f0) / 1e6, result) +plt.plot((ws.frequency_grid - line_f0) / 1e6, ws.measurement_vector) plt.xlabel("Frequency offset [MHz]") plt.ylabel("Spectral radiance [K]") plt.title( @@ -61,7 +59,7 @@ # %% Test assert np.allclose( - result[::100], + ws.measurement_vector[::100], np.array( [ 227.78646795, diff --git a/examples/getting-started/3-disort/1.clearsky-radiance.py b/examples/getting-started/3-disort/1.clearsky-radiance.py index eee6f8459d..ebfbe3a506 100644 --- a/examples/getting-started/3-disort/1.clearsky-radiance.py +++ b/examples/getting-started/3-disort/1.clearsky-radiance.py @@ -32,7 +32,6 @@ toa=toa, basename="planets/Earth/afgl/tropical/", missing_is_zero=1 ) ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37") -# ws.atmospheric_field[pyarts.arts.AtmKey.t] = 300.0 # %% Checks and settings @@ -42,7 +41,7 @@ # %% Core Disort calculations -ws.disort_spectral_radiance_fieldClearsky( +ws.disort_spectral_radiance_fieldSunlessClearsky( longitude=lon, latitude=lat, disort_quadrature_dimension=NQuad, diff --git a/examples/getting-started/3-disort/2.clearsky-flux.py b/examples/getting-started/3-disort/2.clearsky-flux.py index eb52955465..2b8512d87e 100644 --- a/examples/getting-started/3-disort/2.clearsky-flux.py +++ b/examples/getting-started/3-disort/2.clearsky-flux.py @@ -41,7 +41,7 @@ # %% Core Disort calculations -ws.disort_settings_agendaSet(option="Clearsky") +ws.disort_settings_agendaSet(option="SunlessClearsky") ws.ray_pathGeometricDownlooking(longitude=lon, latitude=lat, max_step=40_000) ws.ray_path_atmospheric_pointFromPath() diff --git a/python/src/pyarts/plots/ppvar_atm.py b/python/src/pyarts/plots/ppvar_atm.py index 937ad1349e..bbd691018f 100644 --- a/python/src/pyarts/plots/ppvar_atm.py +++ b/python/src/pyarts/plots/ppvar_atm.py @@ -41,22 +41,22 @@ def plt_info(atm, keys=[]): for key in keys: if isinstance(key, pyarts.arts.ArrayOfSpeciesTag): out[key] = Info("linear", f"VMR {key} [-]", 1) - elif isinstance(key, pyarts.arts.options.AtmKey): - if key == pyarts.arts.options.AtmKey.t: + elif isinstance(key, pyarts.arts.AtmKey): + if key == pyarts.arts.AtmKey.t: out[key] = Info("linear", "Temperature [K]", 1) - elif key == pyarts.arts.options.AtmKey.p: + elif key == pyarts.arts.AtmKey.p: out[key] = Info("log", "Pressure [Pa]", 1) - elif key == pyarts.arts.options.AtmKey.mag_u: + elif key == pyarts.arts.AtmKey.mag_u: out[key] = Info("linear", "Magnetic u-component [µT]", 1e6) - elif key == pyarts.arts.options.AtmKey.mag_v: + elif key == pyarts.arts.AtmKey.mag_v: out[key] = Info("linear", "Magnetic v-component [µT]", 1e6) - elif key == pyarts.arts.options.AtmKey.mag_w: + elif key == pyarts.arts.AtmKey.mag_w: out[key] = Info("linear", "Magnetic w-component [µT]", 1e6) - elif key == pyarts.arts.options.AtmKey.wind_u: + elif key == pyarts.arts.AtmKey.wind_u: out[key] = Info("linear", "Wind u-component [m/s]", 1) - elif key == pyarts.arts.options.AtmKey.wind_v: + elif key == pyarts.arts.AtmKey.wind_v: out[key] = Info("linear", "Wind v-component [m/s]", 1) - elif key == pyarts.arts.options.AtmKey.wind_w: + elif key == pyarts.arts.AtmKey.wind_w: out[key] = Info("linear", "Wind w-component [m/s]", 1) else: assert False, "Unknown key type" diff --git a/python/test/workspace/test_variables.py b/python/test/workspace/test_variables.py index 282f8367fc..a1115e49d8 100644 --- a/python/test/workspace/test_variables.py +++ b/python/test/workspace/test_variables.py @@ -169,6 +169,7 @@ def test_method_agenda_variables_io(self): vars = list(allvars.keys()) mets = pyarts.arts.globals.workspace_methods() ages = pyarts.arts.globals.workspace_agendas() + opts = pyarts.arts.globals.option_groups() res = {} for var in vars: @@ -189,6 +190,10 @@ def test_method_agenda_variables_io(self): if var in ages[ag].input or var in ages[ag].output: ok_anyways = True break + + # Allow options to be pure input + if nout == 0 and allvars[var].type in opts: + ok_anyways = True if nin != 0 and allvars[var].type == "Agenda": ok_anyways = True @@ -196,6 +201,7 @@ def test_method_agenda_variables_io(self): if not ok_anyways: bad_vars.append(f"{var}") + bad_vars.sort() assert len(bad_vars) == 0, ( "Should have no pure input/output variables,\n" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 32a7b9ef7d..b73a35b734 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -22,6 +22,7 @@ add_custom_target(UtilityHeadersArts SOURCES ${HEADERFILES}) add_library(workspace_group STATIC workspace_groups.cpp ) +target_link_libraries(workspace_group PRIVATE arts_options) add_executable(make_auto_wsg make_auto_wsg.cpp) add_custom_command(OUTPUT diff --git a/src/core/jacobian/jacobian.cc b/src/core/jacobian/jacobian.cc index 5deb87ea85..2a53bffb43 100644 --- a/src/core/jacobian/jacobian.cc +++ b/src/core/jacobian/jacobian.cc @@ -163,98 +163,21 @@ void default_sensor_x_set(ExhaustiveVectorView x, const ArrayOfSensorObsel& sensor, const SensorKey& key) { const SensorObsel& v = sensor.at(key.elem); + const Vector& o = key.original_grid; using enum SensorKeyType; switch (key.model) { - case SensorJacobianModelType::None: { + case SensorJacobianModelType::None: v.flat(x, key.type); break; + case SensorJacobianModelType::PolynomialOffset: switch (key.type) { - case Frequency: - ARTS_USER_ERROR_IF(x.size() != v.f_grid().size(), - "Bad size. x.size(): {}, f_grid().size(): {}", - x.size(), - v.f_grid().size()) - x = v.f_grid(); - break; - case PointingZenith: - ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), - "Bad size. x.size(): {}, poslos_grid().size(): {}", - x.size(), - v.poslos_grid().size()) - std::transform(v.poslos_grid().begin(), - v.poslos_grid().end(), - x.begin(), - [](auto& poslos) { return poslos.los[0]; }); - break; - case PointingAzimuth: - ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), - "Bad size. x.size(): {}, poslos_grid().size(): {}", - x.size(), - v.poslos_grid().size()) - std::transform(v.poslos_grid().begin(), - v.poslos_grid().end(), - x.begin(), - [](auto& poslos) { return poslos.los[1]; }); - break; - case PointingAltitude: - ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), - "Bad size. x.size(): {}, poslos_grid().size(): {}", - x.size(), - v.poslos_grid().size()) - std::transform(v.poslos_grid().begin(), - v.poslos_grid().end(), - x.begin(), - [](auto& poslos) { return poslos.pos[0]; }); - break; - case PointingLatitude: - ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), - "Bad size. x.size(): {}, poslos_grid().size(): {}", - x.size(), - v.poslos_grid().size()) - std::transform(v.poslos_grid().begin(), - v.poslos_grid().end(), - x.begin(), - [](auto& poslos) { return poslos.pos[1]; }); - break; - case PointingLongitude: - ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), - "Bad size. x.size(): {}, poslos_grid().size(): {}", - x.size(), - v.poslos_grid().size()) - std::transform(v.poslos_grid().begin(), - v.poslos_grid().end(), - x.begin(), - [](auto& poslos) { return poslos.pos[2]; }); - break; - case Weights: { - auto& w = v.weight_matrix(); - ARTS_USER_ERROR_IF(w.size() * 4 != x.size(), - "Bad size. x.size(): {} for w.size()*4: {}", - x.size(), - w.size() * 4) - auto X = x.reshape_as(w.nrows(), w.ncols(), 4); - for (Index i = 0; i < X.npages(); ++i) { - for (Index j = 0; j < X.nrows(); ++j) { - X(i, j, 0) = w(i, j).I(); - X(i, j, 1) = w(i, j).Q(); - X(i, j, 2) = w(i, j).U(); - X(i, j, 3) = w(i, j).V(); - } - } - } break; + case f: polyfit(x, o, rem_frq(v, o)); break; + case za: polyfit(x, o, rem_zag(v, o)); break; + case aa: polyfit(x, o, rem_aag(v, o)); break; + case alt: polyfit(x, o, rem_alt(v, o)); break; + case lat: polyfit(x, o, rem_lat(v, o)); break; + case lon: polyfit(x, o, rem_lon(v, o)); break; } - } break; - case SensorJacobianModelType::PolynomialOffset: { - const Vector& o = key.original_grid; - switch (key.type) { - case Frequency: polyfit(x, o, rem_frq(v, o)); break; - case PointingZenith: polyfit(x, o, rem_zag(v, o)); break; - case PointingAzimuth: polyfit(x, o, rem_aag(v, o)); break; - case PointingAltitude: polyfit(x, o, rem_alt(v, o)); break; - case PointingLatitude: polyfit(x, o, rem_lat(v, o)); break; - case PointingLongitude: polyfit(x, o, rem_lon(v, o)); break; - case Weights: ARTS_USER_ERROR("Not available for polynomial model."); - } - } break; + break; } } @@ -276,93 +199,6 @@ Vector polynomial_offset_evaluate(const ExhaustiveConstVectorView x, return out; } -void set_frq(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - ARTS_USER_ERROR_IF(x.size() != v.f_grid().size(), - "Bad size. x.size(): {}, f_grid().size(): {}", - x.size(), - v.f_grid().size()) - - const auto xs = std::make_shared( - x.begin(), x.end(), [](auto& x) { return x; }); - - // Must copy, as we may change the shared_ptr later - const auto fs = v.f_grid_ptr(); - - for (auto& elem : sensor) { - if (elem.f_grid_ptr() == fs) { - elem.set_f_grid_ptr(xs); // may change here - } - } -} - -template -void set_poslos(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), - "Bad size. x.size(): {}, poslos_grid().size(): {}", - x.size(), - v.poslos_grid().size()) - - SensorPosLosVector xsv = v.poslos_grid(); - - std::transform(xsv.begin(), - xsv.end(), - x.begin(), - xsv.begin(), - [](auto poslos, Numeric val) { - if constexpr (pos) { - poslos.pos[k] = val; - } else { - poslos.los[k] = val; - } - return poslos; - }); - - const auto xs = std::make_shared(std::move(xsv)); - - // Must copy, as we may change the shared_ptr later - const auto ps = v.poslos_grid_ptr(); - - for (auto& elem : sensor) { - if (elem.poslos_grid_ptr() == ps) { - elem.set_poslos_grid_ptr(ps); - } - } -} - -void set_alt(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - set_poslos(v, sensor, x); -} - -void set_lat(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - set_poslos(v, sensor, x); -} - -void set_lon(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - set_poslos(v, sensor, x); -} - -void set_zag(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - set_poslos(v, sensor, x); -} - -void set_aag(const SensorObsel& v, - ArrayOfSensorObsel& sensor, - const ExhaustiveConstVectorView x) { - set_poslos(v, sensor, x); -} - void default_x_sensor_set(ArrayOfSensorObsel& sensor, const SensorKey& key, const ExhaustiveConstVectorView x) { @@ -370,50 +206,13 @@ void default_x_sensor_set(ArrayOfSensorObsel& sensor, using enum SensorKeyType; switch (key.model) { - case SensorJacobianModelType::None: { - switch (key.type) { - case Frequency: set_frq(v, sensor, x); break; - case PointingZenith: set_zag(v, sensor, x); break; - case PointingAzimuth: set_aag(v, sensor, x); break; - case PointingAltitude: set_alt(v, sensor, x); break; - case PointingLatitude: set_lat(v, sensor, x); break; - case PointingLongitude: set_lon(v, sensor, x); break; - case Weights: { - auto [r, c] = v.weight_matrix().shape(); - - ARTS_USER_ERROR_IF(r * c * 4 != x.size(), - "Bad size. x.size(): {} for w.size()*4: {}", - x.size(), - r * c * 4) - - auto X = x.reshape_as(r, c, 4); - StokvecMatrix ws{r, c}; - for (Index i = 0; i < r; ++i) { - for (Index j = 0; j < c; ++j) { - ws(i, j).I() = X(i, j, 0); - ws(i, j).Q() = X(i, j, 1); - ws(i, j).U() = X(i, j, 2); - ws(i, j).V() = X(i, j, 3); - } - } - - v.set_weight_matrix(std::move(ws)); - } break; - } - } break; + case SensorJacobianModelType::None: + unflatten(sensor, x, v, key.type); + break; case SensorJacobianModelType::PolynomialOffset: { - auto& o = key.original_grid; - Vector r = polynomial_offset_evaluate(x, o); - - switch (key.type) { - case Frequency: set_frq(v, sensor, r); break; - case PointingZenith: set_zag(v, sensor, r); break; - case PointingAzimuth: set_aag(v, sensor, r); break; - case PointingAltitude: set_alt(v, sensor, r); break; - case PointingLatitude: set_lat(v, sensor, r); break; - case PointingLongitude: set_lon(v, sensor, r); break; - case Weights: ARTS_USER_ERROR("Not available for polynomial model."); - } + auto& o = key.original_grid; + const Vector r = polynomial_offset_evaluate(x, o); + unflatten(sensor, r, v, key.type); } break; } } @@ -473,22 +272,29 @@ const std::vector& Targets::line() const { return target(); } +const std::vector& Targets::sensor() const { + return target(); +} + std::vector& Targets::atm() { return target(); } std::vector& Targets::surf() { return target(); } std::vector& Targets::line() { return target(); } +std::vector& Targets::sensor() { return target(); } + void Targets::finalize(const AtmField& atmospheric_field, const SurfaceField& surface_field, const AbsorptionBands&, const ArrayOfSensorObsel& measurement_sensor) { zero_out_x(); - const Size natm = atm().size(); - const Size nsurf = surf().size(); - const Size nline = line().size(); - ARTS_ASSERT(target_count() == (natm + nsurf + nline)); + const Size natm = atm().size(); + const Size nsurf = surf().size(); + const Size nline = line().size(); + const Size nsensor = sensor().size(); + ARTS_ASSERT(target_count() == (natm + nsurf + nline + nsensor)); Size last_size = 0; @@ -528,6 +334,27 @@ void Targets::finalize(const AtmField& atmospheric_field, last_size += t.x_size; } + for (Size i = 0; i < nsensor; i++) { + SensorTarget& t = sensor()[i]; + ARTS_USER_ERROR_IF(std::ranges::any_of(sensor() | std::views::drop(i + 1), + Cmp::eq(t.type), + &SensorTarget::type), + "Multiple targets of the same type: {}", + t.type) + t.x_start = last_size; + switch (t.type.model) { + case SensorJacobianModelType::None: + t.x_size = measurement_sensor.at(t.type.elem).flat_size(t.type.type); + break; + case SensorJacobianModelType::PolynomialOffset: + ARTS_USER_ERROR_IF(t.type.polyorder < 0, + "Must have a polynomial order.") + t.x_size = t.type.polyorder + 1; + break; + } + last_size += t.x_size; + } + finalized = true; throwing_check(last_size); diff --git a/src/core/options/arts_options.cc b/src/core/options/arts_options.cc index fececec53c..d2ff90e8d0 100644 --- a/src/core/options/arts_options.cc +++ b/src/core/options/arts_options.cc @@ -20,13 +20,12 @@ std::vector internal_options_create() { R"(A key for identifying a sensor property )", .values_and_desc = { - Value{"Frequency", "Frequency of the sensor"}, - Value{"PointingZenith", "Pointing Zenith of the sensor"}, - Value{"PointingAzimuth", "Pointing Azimuth of the sensor"}, - Value{"PointingAltitude", "Pointing Altitude of the sensor"}, - Value{"PointingLatitude", "Pointing Latitude of the sensor"}, - Value{"PointingLongitude", "Pointing Longitude of the sensor"}, - Value{"Weights", "Weights of the sensor"}, + Value{"f", "Frequency", "Frequency of the sensor"}, + Value{"za", "PointingZenith", "Pointing Zenith of the sensor"}, + Value{"aa", "PointingAzimuth", "Pointing Azimuth of the sensor"}, + Value{"alt", "PointingAltitude", "Pointing Altitude of the sensor"}, + Value{"lat", "PointingLatitude", "Pointing Latitude of the sensor"}, + Value{"lon", "PointingLongitude", "Pointing Longitude of the sensor"}, }}); opts.emplace_back(EnumeratedOption{ @@ -36,7 +35,8 @@ std::vector internal_options_create() { )", .values_and_desc = { Value{"None", "No model, work purely on sensor data"}, - Value{"PolynomialOffset", "The sensor Jacobian is modeled as a polynomial offset"}, + Value{"PolynomialOffset", + "The sensor Jacobian is modeled as a polynomial offset"}, }}); opts.emplace_back(EnumeratedOption{ @@ -523,12 +523,6 @@ if good cases, so we have provide this selection mechanism to make them match. - *ray_path_observer_agendaExecute* - *spectral_radianceClearskyEmission* -)"}, - Value{"EmissionUnits", R"( - - - *ray_path_observer_agendaExecute* - - *spectral_radianceClearskyEmission* - - *spectral_radianceApplyUnitFromSpectralRadiance* )"}, }, }); @@ -539,10 +533,9 @@ if good cases, so we have provide this selection mechanism to make them match. )", .values_and_desc = { - Value{"Clearsky", R"( + Value{"SunlessClearsky", R"( - - *jacobian_targetsInit* - - *jacobian_targetsFinalize* + - Set *jacobian_targets* empty - *ray_path_atmospheric_pointFromPath* - *ray_path_frequency_gridFromPath* - *ray_path_propagation_matrixFromPath* @@ -710,53 +703,49 @@ which is [W / m :math:`^{2}` Hz sr]. For the description below, :math:`c` is the speed of light, :math:`k` is the Boltzmann constant, :math:`f` is the frequency, :math:`h` is the Planck constant, :math:`F(x)` is the conversion -function, and :math:`[I, Q, U, V]` is the internal representation of *spectral_radiance*. +function, and :math:`[I,\; Q,\; U,\; V]` is the internal representation of *spectral_radiance*. For Rayleigh-Jeans brightness temperature the conversion is: .. math:: - [I, Q, U, V] \rightarrow [F(I), F(Q), F(U), F(V)] - -with + [I,\; Q,\; U,\; V] \rightarrow [F(I),\; F(Q),\; F(U),\; F(V)] .. math:: - F(x)=\frac{c^2 x}{2 k f^2}. + F(x) = \frac{c^2 x}{2 k f^2}. For Planck brightness temperature the conversion is: .. math:: - [I, Q, U, V] \rightarrow [ - F(I), - F( \frac{1}{2}\left( I + Q \right) ) - F( \frac{1}{2}\left( I - Q \right) ), - F( \frac{1}{2}\left( I + U \right) ) - F( \frac{1}{2}\left( I - U \right) ), - F( \frac{1}{2}\left( I + V \right) ) - F( \frac{1}{2}\left( I - V \right) )] - -with + [I,\; Q,\; U,\; V] \rightarrow + \left[ + F(I),\; + F\left( \frac{I + Q}{2} \right) - F\left( \frac{I - Q}{2} \right),\; + F\left( \frac{I + U}{2} \right) - F\left( \frac{I - U}{2} \right),\; + F\left( \frac{I + V}{2} \right) - F\left( \frac{I - V}{2} \right) + \right] .. math:: F(x) = \frac{h f}{k}\left[\log\left(1 + \frac{2hf^3}{c^2x}\right)\right]^{-1}. -The spectral radiance per wavelength [W / m :math:`^{2}` m sr`] is: +The spectral radiance per wavelength [W / m :math:`^{2}` m sr] is: .. math:: - [I, Q, U, V] \rightarrow [F(I), F(Q), F(U), F(V)] - -with + [I,\; Q,\; U,\; V] \rightarrow [F(I),\; F(Q),\; F(U),\; F(V)] .. math:: F(x) = \frac{f^2 x}{c}. -The spectral radiance per wavenumber [W / m :math:`^{2}` m :math:`^{-1}` sr`] is: +The spectral radiance per wavenumber [W / m :math:`^{2}` m :math:`^{-1}` sr] is: .. math:: - [I, Q, U, V] \rightarrow [F(I), F(Q), F(U), F(V)] + [I,\; Q,\; U,\; V] \rightarrow [F(I),\; F(Q),\; F(U),\; F(V)] .. math:: @@ -766,7 +755,7 @@ Lastly, the unit option of course just retains the current state [W / m :math:`^ .. math:: - [I, Q, U, V] \rightarrow [I, Q, U, V] + [I,\; Q,\; U,\; V] \rightarrow [I,\; Q,\; U,\; V]. )", .values_and_desc = { @@ -780,11 +769,11 @@ Lastly, the unit option of course just retains the current state [W / m :math:`^ )"}, Value{"W_m2_m_sr", "W/(m^2 m sr)", - "Spectral radiance wavelength [W / m :math:`^{2}` m sr`]"}, + "Spectral radiance wavelength [W / m :math:`^{2}` m sr]"}, Value{ "W_m2_m1_sr", "W/(m^2 m-1 sr)", - "Spectral radiance wavenumber [W / m :math:`^{2}` m :math:`^{-1}` sr`]"}, + "Spectral radiance wavenumber [W / m :math:`^{2}` m :math:`^{-1}` sr]"}, Value{"unit", "1", "Unit spectral radiance [W / m :math:`^{2}` Hz sr]"}, @@ -1114,18 +1103,14 @@ const std::vector& internal_options() { } std::string_view EnumeratedOption::sz() const { - if (values_and_desc.size() < std::numeric_limits::max()) { - return "char"; - } - if (values_and_desc.size() < std::numeric_limits::max()) { - return "short"; - } - if (values_and_desc.size() < std::numeric_limits::max()) { - return "int"; - } - if (values_and_desc.size() < std::numeric_limits::max()) { - return "long"; - } + const auto n = values_and_desc.size(); + + if (n < std::numeric_limits::max()) return "bool"; + if (n < std::numeric_limits::max()) return "char"; + if (n < std::numeric_limits::max()) return "short"; + if (n < std::numeric_limits::max()) return "int"; + if (n < std::numeric_limits::max()) return "long"; + if (n < std::numeric_limits::max()) return "long long"; throw std::runtime_error("Too many values for enum class"); } diff --git a/src/core/sensor/obsel.cpp b/src/core/sensor/obsel.cpp index 76c6ba94ab..8bc3c65afe 100644 --- a/src/core/sensor/obsel.cpp +++ b/src/core/sensor/obsel.cpp @@ -79,13 +79,12 @@ void Obsel::sumup(VectorView out, const StokvecMatrixView& j, Index ip) const { Size Obsel::flat_size(const SensorKeyType& key) const { switch (key) { using enum SensorKeyType; - case Frequency: return f->size(); - case PointingZenith: return poslos->size(); - case PointingAzimuth: return poslos->size(); - case PointingAltitude: return poslos->size(); - case PointingLatitude: return poslos->size(); - case PointingLongitude: return poslos->size(); - case Weights: return f->size() * poslos->size() * 4; + case f: return this->f->size(); + case za: return poslos->size(); + case aa: return poslos->size(); + case alt: return poslos->size(); + case lat: return poslos->size(); + case lon: return poslos->size(); } std::unreachable(); @@ -94,14 +93,14 @@ Size Obsel::flat_size(const SensorKeyType& key) const { void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { switch (key) { using enum SensorKeyType; - case Frequency: + case f: ARTS_USER_ERROR_IF(x.size() != f_grid().size(), "Bad size. x.size(): {}, f_grid().size(): {}", x.size(), f_grid().size()) x = f_grid(); break; - case PointingZenith: + case za: ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), "Bad size. x.size(): {}, poslos_grid().size(): {}", x.size(), @@ -111,7 +110,7 @@ void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { x.begin(), [](auto& poslos) { return poslos.los[0]; }); break; - case PointingAzimuth: + case aa: ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), "Bad size. x.size(): {}, poslos_grid().size(): {}", x.size(), @@ -121,7 +120,7 @@ void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { x.begin(), [](auto& poslos) { return poslos.los[1]; }); break; - case PointingAltitude: + case alt: ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), "Bad size. x.size(): {}, poslos_grid().size(): {}", x.size(), @@ -131,7 +130,7 @@ void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { x.begin(), [](auto& poslos) { return poslos.pos[0]; }); break; - case PointingLatitude: + case lat: ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), "Bad size. x.size(): {}, poslos_grid().size(): {}", x.size(), @@ -141,7 +140,7 @@ void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { x.begin(), [](auto& poslos) { return poslos.pos[1]; }); break; - case PointingLongitude: + case lon: ARTS_USER_ERROR_IF(x.size() != poslos_grid().size(), "Bad size. x.size(): {}, poslos_grid().size(): {}", x.size(), @@ -151,21 +150,6 @@ void Obsel::flat(ExhaustiveVectorView x, const SensorKeyType& key) const { x.begin(), [](auto& poslos) { return poslos.pos[2]; }); break; - case Weights: { - ARTS_USER_ERROR_IF(w.size() * 4 != x.size(), - "Bad size. x.size(): {} for w.size()*4: {}", - x.size(), - w.size() * 4) - auto X = x.reshape_as(w.nrows(), w.ncols(), 4); - for (Index i = 0; i < X.npages(); ++i) { - for (Index j = 0; j < X.nrows(); ++j) { - X(i, j, 0) = w(i, j).I(); - X(i, j, 1) = w(i, j).Q(); - X(i, j, 2) = w(i, j).U(); - X(i, j, 3) = w(i, j).V(); - } - } - } break; } } @@ -174,6 +158,22 @@ Vector Obsel::flat(const SensorKeyType& key) const { flat(out, key); return out; } + +Index Obsel::find(const Vector3& pos, const Vector2& los) const { + const auto first = poslos->cbegin(); + const auto last = poslos->cend(); + const auto same_poslos = [&](const auto& p) { + return pos == p.pos and los == p.los; + }; + + const auto it = std::find_if(first, last, same_poslos); + + return (it == last) ? dont_have : std::distance(first, it); +} + +Index Obsel::find(const AscendingGrid& frequency_grid) const { + return (*f != frequency_grid) ? dont_have : 0; +} } // namespace sensor SensorSimulations collect_simulations(const ArrayOfSensorObsel& obsels) { @@ -251,3 +251,105 @@ void make_exhaustive(ArrayOfSensorObsel& obsels) { obsel = SensorObsel(f_grid_ptr, poslos_grid_ptr, std::move(weights)); } } + +void set_frq(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + ARTS_USER_ERROR_IF(x.size() != v.f_grid().size(), + "Bad size. x.size(): {}, f_grid().size(): {}", + x.size(), + v.f_grid().size()) + + const auto xs = std::make_shared( + x.begin(), x.end(), [](auto& x) { return x; }); + + // Must copy, as we may change the shared_ptr later + const auto fs = v.f_grid_ptr(); + + for (auto& elem : sensor) { + if (elem.f_grid_ptr() == fs) { + elem.set_f_grid_ptr(xs); // may change here + } + } +} + +template +void set_poslos(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + ARTS_USER_ERROR_IF(x.size() != v.poslos_grid().size(), + "Bad size. x.size(): {}, poslos_grid().size(): {}", + x.size(), + v.poslos_grid().size()) + + SensorPosLosVector xsv = v.poslos_grid(); + + std::transform(xsv.begin(), + xsv.end(), + x.begin(), + xsv.begin(), + [](auto poslos, Numeric val) { + if constexpr (pos) { + poslos.pos[k] = val; + } else { + poslos.los[k] = val; + } + return poslos; + }); + + const auto xs = std::make_shared(std::move(xsv)); + + // Must copy, as we may change the shared_ptr later + const auto ps = v.poslos_grid_ptr(); + + for (auto& elem : sensor) { + if (elem.poslos_grid_ptr() == ps) { + elem.set_poslos_grid_ptr(ps); + } + } +} + +void set_alt(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + set_poslos(v, sensor, x); +} + +void set_lat(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + set_poslos(v, sensor, x); +} + +void set_lon(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + set_poslos(v, sensor, x); +} + +void set_zag(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + set_poslos(v, sensor, x); +} + +void set_aag(const SensorObsel& v, + ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView x) { + set_poslos(v, sensor, x); +} + +void unflatten(ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView& x, + const SensorObsel& v, + const SensorKeyType& key) { + switch (key) { + using enum SensorKeyType; + case f: set_frq(v, sensor, x); break; + case za: set_zag(v, sensor, x); break; + case aa: set_aag(v, sensor, x); break; + case alt: set_alt(v, sensor, x); break; + case lat: set_lat(v, sensor, x); break; + case lon: set_lon(v, sensor, x); break; + } +} \ No newline at end of file diff --git a/src/core/sensor/obsel.h b/src/core/sensor/obsel.h index 805bd5e640..6d8e87d746 100644 --- a/src/core/sensor/obsel.h +++ b/src/core/sensor/obsel.h @@ -1,5 +1,7 @@ #pragma once +#include +#include #include #include @@ -9,9 +11,6 @@ #include #include -#include -#include - namespace sensor { struct PosLos { Vector3 pos; @@ -138,6 +137,9 @@ class Obsel { [[nodiscard]] Size flat_size(const SensorKeyType& key) const; void flat(ExhaustiveVectorView x, const SensorKeyType& key) const; [[nodiscard]] Vector flat(const SensorKeyType& key) const; + + [[nodiscard]] Index find(const Vector3&, const Vector2&) const; + [[nodiscard]] Index find(const AscendingGrid& frequency_grid) const; }; std::ostream& operator<<(std::ostream& os, const Array& obsel); @@ -157,10 +159,12 @@ struct SensorKey { bool operator==(const SensorKey& other) const; }; -template<> +template <> struct std::hash { std::size_t operator()(const SensorKey& g) const { - return std::hash{}(g.type) ^ (std::hash{}(g.model) << 10) ^ (std::hash{}(g.elem) << 20); + return std::hash{}(g.type) ^ + (std::hash{}(g.model) << 10) ^ + (std::hash{}(g.elem) << 20); } }; @@ -190,6 +194,11 @@ using SensorSimulations = std::unordered_map< std::shared_ptr, std::unordered_set>>; +void unflatten(ArrayOfSensorObsel& sensor, + const ExhaustiveConstVectorView& x, + const SensorObsel& v, + const SensorKeyType& key); + void make_exhaustive(ArrayOfSensorObsel& obsels); SensorSimulations collect_simulations(const ArrayOfSensorObsel& obsels); diff --git a/src/m_covmat.cc b/src/m_covmat.cc index aad910c7d4..44f3ac3406 100644 --- a/src/m_covmat.cc +++ b/src/m_covmat.cc @@ -19,15 +19,15 @@ void add_diagonal_covmat(CovarianceMatrix& covmat, ARTS_USER_ERROR_IF( matrix.ncols() != colrow.extent or matrix.nrows() != colrow.extent, - R"(The matrix must be square. It must also have the same size as the target. + R"(The matrix must be square. It must also have the same size as the target. shape(matrix) = {:B,}, shape(target) = [{}, {}] Target: {} )", - matrix.shape(), - colrow.extent, - colrow.extent, - target.type); + matrix.shape(), + colrow.extent, + colrow.extent, + target.type); covmat.add_correlation({colrow, colrow, @@ -37,15 +37,15 @@ Target: {} if (inverse.not_null()) { ARTS_USER_ERROR_IF( inverse.ncols() != colrow.extent or inverse.nrows() != colrow.extent, - R"(The inverse matrix must be square. It must also have the same size as the target. + R"(The inverse matrix must be square. It must also have the same size as the target. shape(matrix) = {:B,}, shape(target) = [{}, {}] Target: {} )", - inverse.shape(), - colrow.extent, - colrow.extent, - target.type); + inverse.shape(), + colrow.extent, + colrow.extent, + target.type); covmat.add_correlation({colrow, colrow, @@ -128,6 +128,29 @@ void model_state_covariance_matrixAdd( not found, "No target found for surface target : {}", new_target); } +void model_state_covariance_matrixAdd( + CovarianceMatrix& model_state_covariance_matrix, + const JacobianTargets& jacobian_targets, + const SensorKey& new_target, + const BlockMatrix& matrix, + const BlockMatrix& inverse) { + ARTS_USER_ERROR_IF(not jacobian_targets.finalized, + "Jacobian targets not finalized."); + + bool found = false; + + for (const auto& target : jacobian_targets.sensor()) { + if (target.type == new_target) { + found = true; + add_diagonal_covmat( + model_state_covariance_matrix, target, matrix, inverse); + } + } + + ARTS_USER_ERROR_IF( + not found, "No target found for sensor target : {}", new_target); +} + void model_state_covariance_matrixAddAtmosphere( CovarianceMatrix& model_state_covariance_matrix, const JacobianTargets& jacobian_targets, @@ -285,9 +308,13 @@ void RetrievalFinalizeDiagonal(CovarianceMatrix& model_state_covariance_matrix, covariance_matrix_diagonal_blocks, const AtmField& atmospheric_field, const SurfaceField& surface_field, - const AbsorptionBands& absorption_bands) { - jacobian_targetsFinalize( - jacobian_targets, atmospheric_field, surface_field, absorption_bands); + const AbsorptionBands& absorption_bands, + const ArrayOfSensorObsel& measurement_sensor) { + jacobian_targetsFinalize(jacobian_targets, + atmospheric_field, + surface_field, + absorption_bands, + measurement_sensor); for (auto& key_data : covariance_matrix_diagonal_blocks) { std::visit( diff --git a/src/m_jactargets.cc b/src/m_jactargets.cc index 67cdfade33..f8c04c9654 100644 --- a/src/m_jactargets.cc +++ b/src/m_jactargets.cc @@ -18,9 +18,10 @@ void jacobian_targetsInit(JacobianTargets& jacobian_targets) { void jacobian_targetsFinalize(JacobianTargets& jacobian_targets, const AtmField& atmospheric_field, const SurfaceField& surface_field, - const AbsorptionBands& absorption_bands) { + const AbsorptionBands& absorption_bands, + const ArrayOfSensorObsel& measurement_sensor) { jacobian_targets.finalize( - atmospheric_field, surface_field, absorption_bands, {}); + atmospheric_field, surface_field, absorption_bands, measurement_sensor); } void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, @@ -147,7 +148,9 @@ void jacobian_targetsAddSensor(JacobianTargets& jacobian_targets, elem) jacobian_targets.sensor().push_back({ - .type = {.type = key, .elem = elem, .model = SensorJacobianModelType::None}, + .type = {.type = key, + .elem = elem, + .model = SensorJacobianModelType::None}, .d = d, .target_pos = jacobian_targets.target_count(), }); @@ -167,11 +170,11 @@ void jacobian_targetsAddSensorPolyFit( polyorder < 0, "Polyorder must be non-negative: {}", polyorder) Jacobian::SensorTarget target{ - .type = {.type = key, - .elem = elem, - .model = SensorJacobianModelType::PolynomialOffset, - .polyorder = polyorder}, - // .original_grid = sensor_grid(measurement_sensor[elem], key)}, + .type = {.type = key, + .elem = elem, + .model = SensorJacobianModelType::PolynomialOffset, + .polyorder = polyorder}, + // .original_grid = sensor_grid(measurement_sensor[elem], key)}, .d = d, .target_pos = jacobian_targets.target_count(), }; diff --git a/src/m_rad.cc b/src/m_rad.cc index ae6b8f5ab8..3a6f337d1a 100644 --- a/src/m_rad.cc +++ b/src/m_rad.cc @@ -8,9 +8,11 @@ #include #include +#include #include "arts_omp.h" #include "debug.h" +#include "enumsSensorJacobianModelType.h" #include "fwd.h" #include "rtepack_stokes_vector.h" #include "sorted_grid.h" @@ -135,7 +137,7 @@ void spectral_radiance_jacobianApplyUnit( const StokvecVector &spectral_radiance, const AscendingGrid &frequency_grid, const PropagationPathPoint &ray_path_point, - const String &spectral_radiance_unit) try { + const SpectralRadianceUnitType &spectral_radiance_unit) try { ARTS_USER_ERROR_IF(spectral_radiance.size() != frequency_grid.size(), "spectral_radiance must have same size as frequency_grid") ARTS_USER_ERROR_IF( @@ -143,9 +145,8 @@ void spectral_radiance_jacobianApplyUnit( spectral_radiance_jacobian.ncols() != frequency_grid.size(), "spectral_radiance must have same size as frequency_grid") - const auto dF = rtepack::dunit_converter( - to(spectral_radiance_unit), - ray_path_point.nreal); + const auto dF = + rtepack::dunit_converter(spectral_radiance_unit, ray_path_point.nreal); //! Must apply the unit to the spectral radiance jacobian first for (Index i = 0; i < spectral_radiance_jacobian.nrows(); i++) { @@ -158,15 +159,15 @@ void spectral_radiance_jacobianApplyUnit( } ARTS_METHOD_ERROR_CATCH -void spectral_radianceApplyUnit(StokvecVector &spectral_radiance, - const AscendingGrid &frequency_grid, - const PropagationPathPoint &ray_path_point, - const String &spectral_radiance_unit) try { +void spectral_radianceApplyUnit( + StokvecVector &spectral_radiance, + const AscendingGrid &frequency_grid, + const PropagationPathPoint &ray_path_point, + const SpectralRadianceUnitType &spectral_radiance_unit) try { ARTS_USER_ERROR_IF(spectral_radiance.size() != frequency_grid.size(), "spectral_radiance must have same size as frequency_grid") - const auto F = rtepack::unit_converter( - to(spectral_radiance_unit), - ray_path_point.nreal); + const auto F = + rtepack::unit_converter(spectral_radiance_unit, ray_path_point.nreal); std::transform(spectral_radiance.begin(), spectral_radiance.end(), @@ -176,36 +177,193 @@ void spectral_radianceApplyUnit(StokvecVector &spectral_radiance, } ARTS_METHOD_ERROR_CATCH +void spectral_radiance_jacobianAddSensorJacobianPerturbations( + const Workspace &ws, + StokvecMatrix &spectral_radiance_jacobian, + const StokvecVector &spectral_radiance, + const ArrayOfSensorObsel &measurement_sensor, + const AscendingGrid &frequency_grid, + const JacobianTargets &jacobian_targets, + const Vector3 &pos, + const Vector2 &los, + const AtmField &atmospheric_field, + const SurfaceField &surface_field, + const Agenda &spectral_radiance_observer_agenda) try { + /* + + This method likely calls itself "recursively" for sensor parameters + + However, the flag for bailing and stopping this recursion + is that there are no more jacobian targets. So it calls + itself always with an empty jacobian_targets. This is how + it bails out of the recursion. + + */ + if (jacobian_targets.sensor().empty()) return; + + ARTS_USER_ERROR_IF( + spectral_radiance.size() != frequency_grid.size(), + R"(spectral_radiance must have same size as element frequency grid + +spectral_radiance.size() = {}, +frequency_grid.size() = {} +)", + spectral_radiance.size(), + frequency_grid.size()) + + ARTS_USER_ERROR_IF( + (spectral_radiance_jacobian.shape() != + std::array{static_cast(jacobian_targets.x_size()), + frequency_grid.size()}), + R"(spectral_radiance_jacobian must be x-grid times frequency grid + +spectral_radiance_jacobian.shape() = {:B,}, +jacobian_targets.x_size() = {}, +frequency_grid.size() = {} +)", + spectral_radiance_jacobian.shape(), + jacobian_targets.x_size(), + frequency_grid.size()) + + const JacobianTargets jacobian_targets_empty{}; + StokvecMatrix spectral_radiance_jacobian_empty{}; + ArrayOfPropagationPathPoint ray_path{}; + + StokvecVector dsrad; + auto call = [&](const AscendingGrid &frequency_grid_2, + const Vector3 &pos2, + const Vector2 &los2, + const Numeric d) { + spectral_radiance_observer_agendaExecute(ws, + dsrad, + spectral_radiance_jacobian_empty, + ray_path, + frequency_grid_2, + jacobian_targets_empty, + pos2, + los2, + atmospheric_field, + surface_field, + spectral_radiance_observer_agenda); + + ARTS_USER_ERROR_IF(dsrad.size() != spectral_radiance.size(), + R"(Wrong size of perturbed spectral radiance: + + dsrad.size() = {}, + spectral_radiance.size() = {} +)", + dsrad.size(), + spectral_radiance.size()) + + // Convert to perturbed Jacobian + dsrad -= spectral_radiance; + for (auto &v : dsrad) v /= d; + return dsrad; + }; + + const auto &x = frequency_grid; + const auto b = x.begin(); + const auto e = x.end(); + + for (auto &target : jacobian_targets.sensor()) { + ARTS_USER_ERROR_IF( + measurement_sensor.size() <= static_cast(target.type.elem), + "Sensor element out of bounds"); + + auto m = spectral_radiance_jacobian.slice(target.x_start, target.x_size); + const Numeric d = target.d; + + // Check that the Jacobian targets are represented by this frequency grid and this pos-los pair + const Index iposlos = measurement_sensor[target.type.elem].find(pos, los); + if (iposlos == SensorObsel::dont_have) continue; + if (measurement_sensor[target.type.elem].find(frequency_grid) == + SensorObsel::dont_have) + continue; + + using enum SensorKeyType; + switch (target.type.type) { + case f: call({b, e, [d](auto x) { return x + d; }}, pos, los, d); + case za: call(x, pos, {los[0] + d, los[1]}, d); break; + case aa: call(x, pos, {los[0], los[1] + d}, d); break; + case alt: call(x, {pos[0] + d, pos[1], pos[2]}, los, d); break; + case lat: call(x, {pos[0], pos[1] + d, pos[2]}, los, d); break; + case lon: call(x, {pos[0], pos[1], pos[2] + d}, los, d); break; + } + + switch (target.type.model) { + using enum SensorJacobianModelType; + case PolynomialOffset: { + const auto &o = target.type.original_grid; + + if (target.type.type == f) { + ARTS_USER_ERROR_IF( + x.size() != o.size(), + "Expects the frequency grid to be the same size as the original grid") + + for (Size i = 0; i < target.x_size; i++) { + for (Index iv = 0; iv < x.size(); iv++) { + m(i, iv) += dsrad[iv] * std::pow(o[iv], i); + } + } + } else { + ARTS_USER_ERROR_IF( + static_cast(target.x_size) != o.size(), + "Expects original grid to be the same as the required x-parameters in the jacobian target") + + for (Size i = 0; i < target.x_size; i++) { + const Numeric g = std::pow(o[i], i); + for (Index iv = 0; iv < x.size(); iv++) { + m(i, iv) += dsrad[iv] * g; + } + } + } + } break; + case None: { + if (target.type.type == f) { + ARTS_USER_ERROR_IF(m.ncols() != m.nrows(), + "Expects square matrix for frequency derivative") + m.diagonal() += dsrad; + } else { + ARTS_USER_ERROR_IF(static_cast(iposlos) >= target.x_size, + "Bad pos-los index"); + m[iposlos] += dsrad; + } + } break; + } + } +} +ARTS_METHOD_ERROR_CATCH + void measurement_vectorFromSensor( const Workspace &ws, Vector &measurement_vector, - Matrix &measurement_vector_jacobian, - const ArrayOfSensorObsel &measurement_vector_sensor, + Matrix &measurement_jacobian, + const ArrayOfSensorObsel &measurement_sensor, const JacobianTargets &jacobian_targets, const AtmField &atmospheric_field, const SurfaceField &surface_field, - const String &spectral_radiance_unit, + const SpectralRadianceUnitType &spectral_radiance_unit, const Agenda &spectral_radiance_observer_agenda) try { - measurement_vector.resize(measurement_vector_sensor.size()); + measurement_vector.resize(measurement_sensor.size()); measurement_vector = 0.0; - measurement_vector_jacobian.resize(measurement_vector_sensor.size(), - jacobian_targets.x_size()); - measurement_vector_jacobian = 0.0; + measurement_jacobian.resize(measurement_sensor.size(), + jacobian_targets.x_size()); + measurement_jacobian = 0.0; - if (measurement_vector_sensor.empty()) return; + if (measurement_sensor.empty()) return; //! Check the observational elements that their dimensions are correct - for (auto &obsel : measurement_vector_sensor) obsel.check(); + for (auto &obsel : measurement_sensor) obsel.check(); - const SensorSimulations simulations = - collect_simulations(measurement_vector_sensor); + const SensorSimulations simulations = collect_simulations(measurement_sensor); for (auto &[f_grid_ptr, poslos_set] : simulations) { for (auto &poslos_gs : poslos_set) { for (Index ip = 0; ip < poslos_gs->size(); ++ip) { StokvecVector spectral_radiance; StokvecMatrix spectral_radiance_jacobian; + ArrayOfPropagationPathPoint ray_path; const SensorPosLos &poslos = (*poslos_gs)[ip]; @@ -213,13 +371,13 @@ void measurement_vectorFromSensor( ws, spectral_radiance, spectral_radiance_jacobian, + ray_path, *f_grid_ptr, jacobian_targets, poslos.pos, poslos.los, atmospheric_field, surface_field, - spectral_radiance_unit, spectral_radiance_observer_agenda); ARTS_USER_ERROR_IF( @@ -232,27 +390,32 @@ f_grid_ptr->size() = {} spectral_radiance.size(), f_grid_ptr->size()) ARTS_USER_ERROR_IF( - spectral_radiance_jacobian.nrows() != - measurement_vector_jacobian.ncols() or - spectral_radiance_jacobian.ncols() != f_grid_ptr->size(), + (spectral_radiance_jacobian.shape() != + std::array{measurement_jacobian.ncols(), f_grid_ptr->size()}), R"(spectral_radiance_jacobian must be targets x frequency grid size spectral_radiance_jacobian.shape() = {:B,}, f_grid_ptr->size() = {}, -measurement_vector_jacobian.ncols() = {} +measurement_jacobian.ncols() = {} )", spectral_radiance_jacobian.shape(), f_grid_ptr->size(), - measurement_vector_jacobian.ncols()) + measurement_jacobian.ncols()) + + spectral_radianceApplyUnitFromSpectralRadiance( + spectral_radiance, + spectral_radiance_jacobian, + *f_grid_ptr, + ray_path, + spectral_radiance_unit); - for (Size iv = 0; iv < measurement_vector_sensor.size(); ++iv) { - const SensorObsel &obsel = measurement_vector_sensor[iv]; + for (Size iv = 0; iv < measurement_sensor.size(); ++iv) { + const SensorObsel &obsel = measurement_sensor[iv]; if (obsel.same_freqs(f_grid_ptr)) { measurement_vector[iv] += obsel.sumup(spectral_radiance, ip); - obsel.sumup(measurement_vector_jacobian[iv], - spectral_radiance_jacobian, - ip); + obsel.sumup( + measurement_jacobian[iv], spectral_radiance_jacobian, ip); } } } diff --git a/src/python_interface/gen_auto_py_options.cpp b/src/python_interface/gen_auto_py_options.cpp index 9aa7a8256d..3ca5fe5172 100644 --- a/src/python_interface/gen_auto_py_options.cpp +++ b/src/python_interface/gen_auto_py_options.cpp @@ -14,7 +14,7 @@ os << std::format(R"-x-( void enum_{0}(py::module_& m) {{ py::class_<{0}> _g{0}(m, "{0}"); - _g{0}.doc() = R"-ENUMDOC-({1})-ENUMDOC-"; + _g{0}.doc() = PythonWorkspaceGroupInfo<{0}>::desc; xml_interface(_g{0}); @@ -55,7 +55,7 @@ void enum_{0}(py::module_& m) {{ _g{0}.def_static("get_options_as_strings", [](){{return enumstrs::{0}Names<>;}}, "Get a list of all options as strings"); -)-x-", wso.name, unwrap_stars(wso.docs())); +)-x-", wso.name); static std::array pykeywords{"None", "any", "all", "print"}; constexpr std::string_view ignore_str = diff --git a/src/python_interface/py_global.cpp b/src/python_interface/py_global.cpp index dd3d35effa..08a3548c01 100644 --- a/src/python_interface/py_global.cpp +++ b/src/python_interface/py_global.cpp @@ -1,5 +1,7 @@ #include +#include +#include #include #include #include @@ -7,7 +9,6 @@ #include #include -#include "nanobind/nanobind.h" #include "python_interface.h" struct global_data { @@ -127,6 +128,18 @@ Return "Return\n------\n:class:`dict`" "\n Map of agendas"); + global.def( + "option_groups", + []() -> std::vector { + std::vector out(internal_options().size()); + for (Size i = 0; i < out.size(); i++) { + out[i] = internal_options()[i].name; + } + return out; + }, + "Get a copy of all named options\n\n" + "Return\n------\n:class:`list`"); + global.def( "all_isotopologues", []() { return Species::Isotopologues; }, @@ -164,13 +177,16 @@ Parameters .def_ro_static("is_lgpl", &global_data::is_lgpl, "Whether the ARTS library is licensed under the LGPL") - .def("__repr__", [](const global_data&) { - return std::format(R"(Global state of ARTS: + .def("__repr__", + [](const global_data&) { + return std::format(R"(Global state of ARTS: is_lgpl: {} )", - global_data::is_lgpl); - }).doc() = "A set of global data that we might need from ARTS inside pyarts"; + global_data::is_lgpl); + }) + .doc() = + "A set of global data that we might need from ARTS inside pyarts"; } catch (std::exception& e) { throw std::runtime_error( var_string("DEV ERROR:\nCannot initialize global\n", e.what())); diff --git a/src/python_interface/py_lbl.cpp b/src/python_interface/py_lbl.cpp index 2d283af50a..63739997cb 100644 --- a/src/python_interface/py_lbl.cpp +++ b/src/python_interface/py_lbl.cpp @@ -49,7 +49,7 @@ void py_lbl(py::module_& m) try { [](lbl::temperature::data& self, LineShapeModelType x) { self = lbl::temperature::data{x, self.X()}; }, - ":class:`~pyarts.arts.options.TemperatureModelType` The type of the model") + ":class:`~pyarts.arts.TemperatureModelType` The type of the model") .def_prop_rw( "data", [](lbl::temperature::data& self) { return self.X(); }, diff --git a/src/python_interface/py_path.cpp b/src/python_interface/py_path.cpp index 2235842a7e..0d6cb9be96 100644 --- a/src/python_interface/py_path.cpp +++ b/src/python_interface/py_path.cpp @@ -13,11 +13,11 @@ void py_path(py::module_& m) try { pppp.def_rw( "pos_type", &PropagationPathPoint::pos_type, - ":class:`~pyarts.arts.options.PathPositionType` Path position type") + ":class:`~pyarts.arts.PathPositionType` Path position type") .def_rw( "los_type", &PropagationPathPoint::los_type, - ":class:`~pyarts.arts.options.PathPositionType` Path line-of-sight type") + ":class:`~pyarts.arts.PathPositionType` Path line-of-sight type") .def_rw("pos", &PropagationPathPoint::pos, ":class:`~pyarts.arts.Vector3` Path position") diff --git a/src/python_interface/py_quantum.cpp b/src/python_interface/py_quantum.cpp index 2650413c90..8eeecd6976 100644 --- a/src/python_interface/py_quantum.cpp +++ b/src/python_interface/py_quantum.cpp @@ -22,7 +22,7 @@ void py_quantum(py::module_& m) try { .PythonInterfaceBasicRepresentation(QuantumNumberValue) .def_rw("type", &QuantumNumberValue::type, - ":class:`~pyarts.arts.options.QuantumNumberType` Type of number") + ":class:`~pyarts.arts.QuantumNumberType` Type of number") .def_prop_rw( "str_upp", &QuantumNumberValue::str_upp, diff --git a/src/python_interface/py_species.cpp b/src/python_interface/py_species.cpp index f9442ea57b..52778b2292 100644 --- a/src/python_interface/py_species.cpp +++ b/src/python_interface/py_species.cpp @@ -149,7 +149,7 @@ void py_species(py::module_& m) try { stag.def_rw("spec_ind", &SpeciesTag::spec_ind, ":class:`int` Species index") .def_rw("type", &SpeciesTag::type, - ":class:`~pyarts.arts.options.SpeciesTagType` Type of tag") + ":class:`~pyarts.arts.SpeciesTagType` Type of tag") .def_rw("cia_2nd_species", &SpeciesTag::cia_2nd_species, ":class:`~pyarts.arts.Species` CIA species") diff --git a/src/python_interface/py_spectroscopy.cpp b/src/python_interface/py_spectroscopy.cpp index d757d8d1fd..a2f7f5e711 100644 --- a/src/python_interface/py_spectroscopy.cpp +++ b/src/python_interface/py_spectroscopy.cpp @@ -45,7 +45,7 @@ void py_spectroscopy(py::module_& m) try { .def_rw( "type", &LineShapeModelParameters::type, - ":class:`~pyarts.arts.options.LineShapeTemperatureModel` The temperature model") + ":class:`~pyarts.arts.LineShapeTemperatureModel` The temperature model") .def_rw( "X0", &LineShapeModelParameters::X0, ":class:`float` 1st coefficient") .def_rw( @@ -566,22 +566,22 @@ void py_spectroscopy(py::module_& m) try { .def_rw( "cutoff", &AbsorptionLines::cutoff, - ":class:`~pyarts.arts.options.AbsorptionCutoffTypeOld` Cutoff type") + ":class:`~pyarts.arts.AbsorptionCutoffTypeOld` Cutoff type") .def_rw( "mirroring", &AbsorptionLines::mirroring, - ":class:`~pyarts.arts.options.AbsorptionMirroringype` Mirroring type") + ":class:`~pyarts.arts.AbsorptionMirroringype` Mirroring type") .def_rw( "population", &AbsorptionLines::population, - ":class:`~pyarts.arts.options.AbsorptionPopulationTypeOld` Line population distribution") + ":class:`~pyarts.arts.AbsorptionPopulationTypeOld` Line population distribution") .def_rw( "normalization", &AbsorptionLines::normalization, - ":class:`~pyarts.arts.options.AbsorptionNormalizationTypeOld` Normalization type") + ":class:`~pyarts.arts.AbsorptionNormalizationTypeOld` Normalization type") .def_rw("lineshapetype", &AbsorptionLines::lineshapetype, - ":class:`~pyarts.arts.options.LineShapeTypeOld` Line shape type") + ":class:`~pyarts.arts.LineShapeTypeOld` Line shape type") .def_rw( "T0", &AbsorptionLines::T0, diff --git a/src/workspace_agenda_creator.cpp b/src/workspace_agenda_creator.cpp index aaf8b8afaa..081ec5fa4b 100644 --- a/src/workspace_agenda_creator.cpp +++ b/src/workspace_agenda_creator.cpp @@ -89,11 +89,7 @@ Agenda get_spectral_radiance_observer_agenda(const std::string& option) { case Emission: agenda.add("ray_path_observer_agendaExecute"); agenda.add("spectral_radianceClearskyEmission"); - break; - case EmissionUnits: - agenda.add("ray_path_observer_agendaExecute"); - agenda.add("spectral_radianceClearskyEmission"); - agenda.add("spectral_radianceApplyUnitFromSpectralRadiance"); + agenda.add("spectral_radiance_jacobianAddSensorJacobianPerturbations"); break; } @@ -158,9 +154,8 @@ Agenda get_disort_settings_agenda(const std::string& option) { using enum disort_settings_agendaPredefined; switch (to(option)) { - case Clearsky: - agenda.add("jacobian_targetsInit"); - agenda.add("jacobian_targetsFinalize"); + case SunlessClearsky: + agenda.set("jacobian_targets", JacobianTargets{}); agenda.add("ray_path_atmospheric_pointFromPath"); agenda.add("ray_path_frequency_gridFromPath"); agenda.add("ray_path_propagation_matrixFromPath"); diff --git a/src/workspace_agendas.cpp b/src/workspace_agendas.cpp index 88c87336c4..eeeea0394f 100644 --- a/src/workspace_agendas.cpp +++ b/src/workspace_agendas.cpp @@ -46,19 +46,23 @@ position and line of sight. The intent of this agenda is to provide a spectral radiance as seen from the observer position and line of sight. +It also outputs the *ray_path* as seen from the observer position and line of sight. +This is useful in-case a call to the destructive *spectral_radianceApplyUnitFromSpectralRadiance* +is warranted + The output must be sized as: - *spectral_radiance* : (*frequency_grid*) - *spectral_radiance_jacobian* : (*jacobian_targets*, *frequency_grid*) +- *ray_path* : (Unknown) )--", - .output = {"spectral_radiance", "spectral_radiance_jacobian"}, + .output = {"spectral_radiance", "spectral_radiance_jacobian", "ray_path"}, .input = {"frequency_grid", "jacobian_targets", "spectral_radiance_observer_position", "spectral_radiance_observer_line_of_sight", "atmospheric_field", - "surface_field", - "spectral_radiance_unit"}, + "surface_field"}, }; wsa_data["spectral_radiance_space_agenda"] = { diff --git a/src/workspace_groups.cpp b/src/workspace_groups.cpp index 940a2eeeaa..1e1fed0206 100644 --- a/src/workspace_groups.cpp +++ b/src/workspace_groups.cpp @@ -1,5 +1,7 @@ #include "workspace_groups.h" +#include + std::unordered_map internal_workspace_groups_creator() { std::unordered_map wsg_data; @@ -182,7 +184,7 @@ internal_workspace_groups_creator() { wsg_data["ScatteringSpeciesProperty"] = { .file = "scattering/properties.h", .desc = "Meta data for scattering spefcies.", -}; + }; wsg_data["ArrayOfScatteringMetaData"] = { .file = "optproperties.h", @@ -209,12 +211,6 @@ about the isotopologue, the absorption scheme, and the frequency limits )--", }; - wsg_data["SurfaceKey"] = { - .file = "enumsSurfaceKey.h", - .desc = R"--(A surface key -)--", - }; - wsg_data["SurfaceTypeTag"] = { .file = "surf.h", .desc = R"--(A surface type @@ -233,30 +229,6 @@ about the isotopologue, the absorption scheme, and the frequency limits )--", }; - wsg_data["SpeciesEnum"] = { - .file = "enumsSpeciesEnum.h", - .desc = R"--(An atmospheric species -)--", - }; - - wsg_data["AtmKey"] = { - .file = "enumsAtmKey.h", - .desc = R"--(An atmospheric key -)--", - }; - - wsg_data["LineShapeModelVariable"] = { - .file = "enumsLineShapeModelVariable.h", - .desc = R"--(A line shape model parameter -)--", - }; - - wsg_data["LineByLineVariable"] = { - .file = "enumsLineByLineVariable.h", - .desc = R"--(An line-by-line variable parameter -)--", - }; - wsg_data["ArrayOfSparse"] = { .file = "matpack_sparse.h", .desc = "A list of *Sparse*\n", @@ -338,12 +310,6 @@ and interpolate the cross-section to other temperatures and pressures )--", }; - wsg_data["AtmKey"] = { - .file = "enumsAtmKey.h", - .desc = R"--(A key for atmospheric data -)--", - }; - wsg_data["AtmField"] = { .file = "atm.h", .desc = R"--(An atmospheric field @@ -1010,11 +976,21 @@ well as the sampling device's polarization response. }; wsg_data["AbsorptionLookupTables"] = { - .file = "lookup_map.h", - .desc = "A map of tables of of lookup calculations.\n", + .file = "lookup_map.h", + .desc = "A map of tables of of lookup calculations.\n", .map_type = true, }; + for (auto& g : internal_options()) { + if (wsg_data.find(g.name) != wsg_data.end()) + throw std::runtime_error("Duplicate workspace group name: " + g.name); + + wsg_data[g.name] = { + .file = "enums.h", + .desc = g.docs(), + }; + } + return wsg_data; } diff --git a/src/workspace_meta_methods.cpp b/src/workspace_meta_methods.cpp index 239fd41758..303a6eb663 100644 --- a/src/workspace_meta_methods.cpp +++ b/src/workspace_meta_methods.cpp @@ -22,17 +22,6 @@ std::vector internal_meta_methods_creator() { .out = {"disort_spectral_flux_field"}, }); - out.push_back(WorkspaceMethodInternalMetaRecord{ - .name = "disort_spectral_flux_fieldClearsky", - .desc = "Use Disort for clearsky calculations of spectral flux field", - .author = {"Richard Larsson"}, - .methods = {"ray_pathGeometricUplooking", - "disort_settings_agendaSet", - "disort_spectral_flux_fieldFromAgenda"}, - .out = {"disort_spectral_flux_field"}, - .preset_gin = {"option"}, - .preset_gin_value = {String{"Clearsky"}}}); - out.push_back(WorkspaceMethodInternalMetaRecord{ .name = "disort_spectral_radiance_fieldFromAgenda", .desc = "Use the disort settings agenda to calculate spectral radiance", @@ -45,7 +34,18 @@ std::vector internal_meta_methods_creator() { }); out.push_back(WorkspaceMethodInternalMetaRecord{ - .name = "disort_spectral_radiance_fieldClearsky", + .name = "disort_spectral_flux_fieldSunlessClearsky", + .desc = "Use Disort for clearsky calculations of spectral flux field", + .author = {"Richard Larsson"}, + .methods = {"ray_pathGeometricUplooking", + "disort_settings_agendaSet", + "disort_spectral_flux_fieldFromAgenda"}, + .out = {"disort_spectral_flux_field"}, + .preset_gin = {"option"}, + .preset_gin_value = {String{"SunlessClearsky"}}}); + + out.push_back(WorkspaceMethodInternalMetaRecord{ + .name = "disort_spectral_radiance_fieldSunlessClearsky", .desc = "Use Disort for clearsky calculations of spectral flux field", .author = {"Richard Larsson"}, .methods = {"ray_pathGeometricDownlooking", @@ -55,11 +55,16 @@ std::vector internal_meta_methods_creator() { "disort_quadrature_angles", "disort_quadrature_weights"}, .preset_gin = {"option"}, - .preset_gin_value = {String{"Clearsky"}}}); + .preset_gin_value = {String{"SunlessClearsky"}}}); out.push_back(WorkspaceMethodInternalMetaRecord{ .name = "spectral_radianceApplyUnitFromSpectralRadiance", - .desc = "Apply unit changes to spectral radiance and its Jacobian", + .desc = R"(Apply unit changes to spectral radiance and its Jacobian + +.. warning:: + This is a destructive method. Any use of it means that it is undefined behavior + to use *spectral_radiance* or *spectral_radiance_jacobian* in future methods. +)", .author = {"Richard Larsson"}, .methods = {"ray_path_pointForeground", "spectral_radiance_jacobianApplyUnit", diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index e9b6bde75d..9495a81551 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -589,9 +589,10 @@ common form of a predefined model. .gin = {"basename", "name_missing", "ignore_missing"}, .gin_type = {"String", "Index", "Index"}, .gin_value = {std::nullopt, Index{1}, Index{0}}, - .gin_desc = {R"--(The path to the split catalog files)--", - R"--(Flag to name models that are missing)--", - R"--(Flag to otherwise (if not name_missing is true) ignore missing models)--"}, + .gin_desc = + {R"--(The path to the split catalog files)--", + R"--(Flag to name models that are missing)--", + R"--(Flag to otherwise (if not name_missing is true) ignore missing models)--"}, }; wsm_data["absorption_bandsFromAbsorbtionLines"] = { @@ -1787,6 +1788,32 @@ The Jacobian variable is all 0s, the background is [1 0 0 0] everywhere "ray_path_point"}, }; + wsm_data["spectral_radiance_jacobianAddSensorJacobianPerturbations"] = { + .desc = R"--(Adds sensor properties to the *spectral_radiance_jacobian*. + +This is done via perturbation based on the input delta values to the sensor +Jacobian targets and a callback to *spectral_radiance_observer_agenda* with +a modified *jacobian_targets*, making it safe to use this method inside +*spectral_radiance_observer_agenda*. +)--", + .author = {"Richard Larsson"}, + .out = {"spectral_radiance_jacobian"}, + .in = + { + "spectral_radiance_jacobian", + "spectral_radiance", + "measurement_sensor", + "frequency_grid", + "jacobian_targets", + "spectral_radiance_observer_position", + "spectral_radiance_observer_line_of_sight", + "atmospheric_field", + "surface_field", + "spectral_radiance_observer_agenda", + }, + .pass_workspace = true, + }; + wsm_data["spectral_radiance_jacobianEmpty"] = { .desc = R"--(Set the cosmic background radiation derivative to empty. @@ -1833,10 +1860,11 @@ Size : (*jacobian_targets*, *frequency_grid*) wsm_data["spectral_radiance_jacobianApplyUnit"] = { .desc = R"(Applies a unit to *spectral_radiance*, returning a new field -See *SpectralRadianceUnitType* and *spectral_radiance_unit* for valid use cases -and limitations. - Also be aware that *spectral_radiance_jacobianApplyUnit* must be called before *spectral_radianceApplyUnit*. + +.. warning:: + This is a destructive method. Any use of it means that it is undefined behavior + to use *spectral_radiance* or *spectral_radiance_jacobian* in future methods. )", .author = {"Richard Larsson"}, .out = {"spectral_radiance_jacobian"}, @@ -1854,6 +1882,10 @@ See *SpectralRadianceUnitType* and *spectral_radiance_unit* for valid use cases and limitations. Also be aware that *spectral_radiance_jacobianApplyUnit* must be called before *spectral_radianceApplyUnit*. + +.. warning:: + This is a destructive method. Any use of it means that it is undefined behavior + to use *spectral_radiance* or *spectral_radiance_jacobian* in future methods. )", .author = {"Richard Larsson"}, .out = {"spectral_radiance"}, @@ -1946,14 +1978,19 @@ matrix to be calculated will work. }; wsm_data["jacobian_targetsFinalize"] = { - .desc = R"--(Finalize *jacobian_targets* for use in RT methods + .desc = R"--(Finalize *jacobian_targets*. + +The finalization computes the size of the required *model_state_vector*. +It is thus necessary if any OEM or other functionality that requires the +building of an actual Jacobian matrix. )--", .author = {"Richard Larsson"}, .out = {"jacobian_targets"}, .in = {"jacobian_targets", "atmospheric_field", "surface_field", - "absorption_bands"}, + "absorption_bands", + "measurement_sensor"}, }; wsm_data["jacobian_targetsAddTemperature"] = { @@ -3064,7 +3101,9 @@ The core calculations happens inside the *spectral_radiance_operator*. The core calculations happens inside the *spectral_radiance_observer_agenda*. -User choices of *spectral_radiance_unit* does not adversely affect this method. +User choices of *spectral_radiance_unit* does not adversely affect this method +unless the *measurement_vector* or *measurement_jacobian* are further modified +before consumption by, e.g., *OEM* )--", .author = {"Richard Larsson"}, .out = {"measurement_vector", "measurement_jacobian"}, @@ -4014,7 +4053,9 @@ Note that you must have set the optical thickness before calling this. }; wsm_data["RetrievalFinalizeDiagonal"] = { - .desc = R"(Add species VMR to the retrieval setup. + .desc = R"(Finalize the retrieval setup. + +See *jacobian_targetsFinalize* for more information. )", .author = {"Richard Larsson"}, .out = {"model_state_covariance_matrix", "jacobian_targets"}, @@ -4022,7 +4063,8 @@ Note that you must have set the optical thickness before calling this. "covariance_matrix_diagonal_blocks", "atmospheric_field", "surface_field", - "absorption_bands"}, + "absorption_bands", + "measurement_sensor"}, }; /* diff --git a/src/workspace_variables.cpp b/src/workspace_variables.cpp index bdef1aee6e..8e4e9bbe49 100644 --- a/src/workspace_variables.cpp +++ b/src/workspace_variables.cpp @@ -11,7 +11,8 @@ internal_workspace_variables() { }; wsv_data["absorption_lookup_table"] = { - .desc = R"--(Absorption lookup table for scalar gas absorption coefficients. + .desc = + R"--(Absorption lookup table for scalar gas absorption coefficients. Precomputing this table replaces the need for the calculation of scalar gas line-by-line absorption. @@ -155,7 +156,7 @@ wind shift Jacobian calculations. The order is ``[df_du, df_dv, df_fw]`` )--", - .type = "Vector3", + .type = "Vector3", .default_value = Vector3{0.0, 0.0, 0.0}, }; @@ -455,9 +456,6 @@ gravity : Numeric wsv_data["spectral_radiance_unit"] = { .desc = R"--(The spectral radiance unit after conversion. -If a unit conversion is desired, the user has to set this variable to one -of the valid options in *SpectralRadianceUnitType*. - Internally, it is always assumed that this is set to "1" and that no unit conversion are taking place. @@ -466,8 +464,8 @@ Unless a method or variable explicitly mention that a unit conversion is support it is called, the use of *spectral_radiance_unit* with a different unit than "1" may lead to undesired results. )--", - .type = "String", - .default_value = String{"1"}, + .type = "SpectralRadianceUnitType", + .default_value = SpectralRadianceUnitType::unit, }; wsv_data["spectral_radiance"] = { diff --git a/src/xml_io_compound_types.cc b/src/xml_io_compound_types.cc index 9bd41d76d7..69b77de726 100644 --- a/src/xml_io_compound_types.cc +++ b/src/xml_io_compound_types.cc @@ -3064,6 +3064,45 @@ void xml_write_to_stream(std::ostream& os_xml, close_tag.write_to_stream(os_xml); } +//=== SensorKey ========================================================= + +//! Reads SensorKey from XML input stream +/*! + \param is_xml XML Input stream + \param key SensorKey + \param pbifs Pointer to binary file stream. NULL for ASCII output. +*/ +void xml_read_from_stream(std::istream& is_xml, + SensorKey& key, + bifstream* pbifs) { + ArtsXMLTag tag; + tag.read_from_stream(is_xml); + tag.check_name("SensorKey"); + + tag.read_from_stream(is_xml); + tag.check_name("/SensorKey"); +} + +//! Write SensorKey to XML output stream +/*! + \param os_xml XML output stream + \param key SensorKey + \param pbofs Pointer to binary file stream. NULL for ASCII output. + \param name Unused +*/ +void xml_write_to_stream(std::ostream& os_xml, + const SensorKey& key, + bofstream* pbofs, + const String& name [[maybe_unused]]) { + ArtsXMLTag open_tag; + ArtsXMLTag close_tag; + + open_tag.set_name("SensorKey"); + + close_tag.set_name("/SensorKey"); + close_tag.write_to_stream(os_xml); +} + //=== JacobianTargetType ========================================================= //! Reads JacobianTargetType from XML input stream @@ -3087,6 +3126,8 @@ void xml_read_from_stream(std::istream& is_xml, jtt.target = SurfaceKeyVal{}; } else if (type == "LblLineKey"s) { jtt.target = LblLineKey{}; + } else if (type == "SensorKey"s) { + jtt.target = SensorKey{}; } else { ARTS_USER_ERROR(R"(Cannot understand type: "{}")", type); } diff --git a/tests/core/jac/full_arts_emission.py b/tests/core/jac/full_arts_emission.py index ec965dde08..4e37bc9044 100644 --- a/tests/core/jac/full_arts_emission.py +++ b/tests/core/jac/full_arts_emission.py @@ -46,6 +46,7 @@ # %% Jacobian +ws.measurement_sensor = [] ws.jacobian_targetsInit() ws.jacobian_targetsAddSpeciesVMR(species="O2") ws.jacobian_targetsFinalize() diff --git a/tests/core/jac/full_optimal_estimation.py b/tests/core/jac/full_optimal_estimation.py index 110debe020..76f71a71e2 100644 --- a/tests/core/jac/full_optimal_estimation.py +++ b/tests/core/jac/full_optimal_estimation.py @@ -39,7 +39,7 @@ # %% Checks and settings ws.spectral_radiance_unit = "Tb" -ws.spectral_radiance_observer_agendaSet(option="EmissionUnits") +ws.spectral_radiance_observer_agendaSet(option="Emission") ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground") ws.spectral_radiance_surface_agendaSet(option="Blackbody") ws.ray_path_observer_agendaSet(option="Geometric") @@ -59,22 +59,23 @@ ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].lon_low = "Nearest" ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].lon_upp = "Nearest" + +# %% Set up sensor + +pos = [100e3, 0, 0] +los = [180.0, 0.0] +ws.measurement_sensorSimple(pos=pos, los=los) + # %% Jacobian ws.RetrievalInit() ws.RetrievalAddSpeciesVMR(species="O2", matrix=np.diag(np.ones((3)) * 5)) ws.RetrievalFinalizeDiagonal() - # %% Core calculations -pos = [100e3, 0, 0] -los = [180.0, 0.0] -ws.ray_pathGeometric(pos=pos, los=los, max_step=1000.0) - for i in range(LIMIT): ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].data = grid - ws.measurement_sensorSimple(pos=pos, los=los) ws.measurement_vectorFromSensor() ws.measurement_vector_fitted = [] diff --git a/tests/core/jac/full_optimal_estimation_multiparam.py b/tests/core/jac/full_optimal_estimation_multiparam.py index b1c2be5977..943c0049c6 100644 --- a/tests/core/jac/full_optimal_estimation_multiparam.py +++ b/tests/core/jac/full_optimal_estimation_multiparam.py @@ -41,7 +41,7 @@ # %% Checks and settings ws.spectral_radiance_unit = "Tb" -ws.spectral_radiance_observer_agendaSet(option="EmissionUnits") +ws.spectral_radiance_observer_agendaSet(option="Emission") ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground") ws.spectral_radiance_surface_agendaSet(option="Blackbody") ws.ray_path_observer_agendaSet(option="Geometric") @@ -57,10 +57,6 @@ ) ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2] = vmr_grid -ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].lat_low = "Nearest" -ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].lat_upp = "Nearest" -ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].lon_low = "Nearest" -ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].lon_upp = "Nearest" # %% Artificial Temperature @@ -78,6 +74,12 @@ ws.atmospheric_field[pyarts.arts.AtmKey.temperature].lon_low = "Nearest" ws.atmospheric_field[pyarts.arts.AtmKey.temperature].lon_upp = "Nearest" +# %% Set up sensor + +pos = [100e3, 0, 0] +los = [180.0, 0.0] +ws.measurement_sensorSimple(pos=pos, los=los) + # %% Jacobian ws.RetrievalInit() @@ -114,14 +116,11 @@ def condition(x, xas): # %% Core calculations -pos = [100e3, 0, 0] -los = [180.0, 0.0] works = False for i in range(LIMIT): ws.atmospheric_field[pyarts.arts.SpeciesEnum.O2].data = vmr_grid ws.atmospheric_field[pyarts.arts.AtmKey.t].data = temp_grid - ws.measurement_sensorSimple(pos=pos, los=los) ws.measurement_vectorFromSensor() ws.measurement_vector_fitted = [] diff --git a/tests/core/jac/model_state_atm.py b/tests/core/jac/model_state_atm.py index 3a67acee84..5346e7876a 100644 --- a/tests/core/jac/model_state_atm.py +++ b/tests/core/jac/model_state_atm.py @@ -47,6 +47,8 @@ def ops(ws, op): ws = pyarts.Workspace() +ws.measurement_sensor = [] + ws.jacobian_targetsInit() ws.atmospheric_fieldInit(toa=toa) diff --git a/tests/core/sensor/jac.py b/tests/core/sensor/jac.py new file mode 100644 index 0000000000..03144efd69 --- /dev/null +++ b/tests/core/sensor/jac.py @@ -0,0 +1,4 @@ +import pyarts +import numpy as np +import matplotlib.pyplot as plt + diff --git a/tests/core/wind/propmat_jac.py b/tests/core/wind/propmat_jac.py index ec3b47a924..58edb57a14 100644 --- a/tests/core/wind/propmat_jac.py +++ b/tests/core/wind/propmat_jac.py @@ -47,7 +47,7 @@ # Set up the wind field Jacobian ws.jacobian_targetsInit() ws.jacobian_targetsAddWindField(component=keys[i]) - ws.jacobian_targetsFinalize() + ws.jacobian_targetsFinalize(measurement_sensor=[]) # Reset ws.atmospheric_point = ws.atmospheric_field.at(*ws.ray_path_point.pos) diff --git a/tests/core/wind/spectral_radiance_jacobian_wind.py b/tests/core/wind/spectral_radiance_jacobian_wind.py index bfcd3ce52c..40c40b301d 100644 --- a/tests/core/wind/spectral_radiance_jacobian_wind.py +++ b/tests/core/wind/spectral_radiance_jacobian_wind.py @@ -40,7 +40,7 @@ # %% Checks and settings ws.spectral_radiance_unit = "Tb" -ws.spectral_radiance_observer_agendaSet(option="EmissionUnits") +ws.spectral_radiance_observer_agendaSet(option="Emission") ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground") ws.spectral_radiance_surface_agendaSet(option="Blackbody") ws.ray_path_observer_agendaSet(option="Geometric") @@ -68,20 +68,20 @@ def inversion_iterate_agenda(ws): ws.measurement_vector_fittedFromMeasurement() +pos = [110e3, 0, 0] +los = [140.0, 30.0] +ws.measurement_sensorSimple(pos=pos, los=los) + for fc in [uf, vf, wf]: ws.RetrievalInit() ws.RetrievalAddWindField(component=str(fc), matrix=np.diag(np.ones((1)) * 100)) ws.RetrievalFinalizeDiagonal() - - pos = [110e3, 0, 0] - los = [140.0, 30.0] fail = True for i in range(LIMIT): ws.atmospheric_field["wind_u"] = wind[fc] ws.atmospheric_field["wind_v"] = wind[fc] ws.atmospheric_field["wind_w"] = wind[fc] - ws.measurement_sensorSimple(pos=pos, los=los) ws.measurement_vectorFromSensor() ws.measurement_vector_fitted = [] diff --git a/tests/core/zeeman/propmat_jac.py b/tests/core/zeeman/propmat_jac.py index 8878baa427..2a7c42ce2d 100644 --- a/tests/core/zeeman/propmat_jac.py +++ b/tests/core/zeeman/propmat_jac.py @@ -30,7 +30,7 @@ ws.atmospheric_field["mag_w"] = B[2] ws.jacobian_targetsInit() ws.jacobian_targetsAddMagneticField(component="w") -ws.jacobian_targetsFinalize() +ws.jacobian_targetsFinalize(measurement_sensor=[]) ws.atmospheric_point.temperature = 250 ws.atmospheric_point.pressure = 1.0 diff --git a/tests/core/zeeman/spectral_radiance_jacobian_magnetic_field.py b/tests/core/zeeman/spectral_radiance_jacobian_magnetic_field.py index 969d5dfdce..2ffc367d92 100644 --- a/tests/core/zeeman/spectral_radiance_jacobian_magnetic_field.py +++ b/tests/core/zeeman/spectral_radiance_jacobian_magnetic_field.py @@ -4,7 +4,7 @@ PLOT = False # Plot for debug LIMIT = 50 # Rerun limit for finding a fit -ATOL = 50 +ATOL = 100 NF = 1001 noise = 0.5 @@ -41,7 +41,7 @@ # %% Checks and settings ws.spectral_radiance_unit = "Tb" -ws.spectral_radiance_observer_agendaSet(option="EmissionUnits") +ws.spectral_radiance_observer_agendaSet(option="Emission") ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground") ws.spectral_radiance_surface_agendaSet(option="Blackbody") ws.ray_path_observer_agendaSet(option="Geometric") @@ -61,13 +61,15 @@ # %% Retrieval agenda - @pyarts.workspace.arts_agenda(ws=ws, fix=True) def inversion_iterate_agenda(ws): ws.UpdateModelStates() ws.measurement_vectorFromSensor() ws.measurement_vector_fittedFromMeasurement() +pos = [110e3, 0, 0] +los = [160.0, 0.0] +ws.measurement_sensorSimple(pos=pos, los=los) for fc in [uf, vf, wf]: ws.RetrievalInit() @@ -76,13 +78,10 @@ def inversion_iterate_agenda(ws): ) ws.RetrievalFinalizeDiagonal() - pos = [110e3, 0, 0] - los = [160.0, 0.0] fail = True for i in range(LIMIT): ws.atmospheric_field["mag_" + str(fc)] = mag[fc] - ws.measurement_sensorSimple(pos=pos, los=los) ws.measurement_vectorFromSensor() ws.measurement_vector_fitted = [] From 59eefb8f9c51677ff0fb2ac9592217b650b1e81d Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Wed, 20 Nov 2024 09:08:57 +0100 Subject: [PATCH 3/7] Make init-add split --- python/test/workspace/test_variables.py | 1 + src/m_obsel.cc | 76 ++++++++++++++----------- src/workspace_meta_methods.cpp | 27 +++++++++ src/workspace_methods.cpp | 48 +++++++++++++--- 4 files changed, 110 insertions(+), 42 deletions(-) diff --git a/python/test/workspace/test_variables.py b/python/test/workspace/test_variables.py index a1115e49d8..8c72eb2059 100644 --- a/python/test/workspace/test_variables.py +++ b/python/test/workspace/test_variables.py @@ -195,6 +195,7 @@ def test_method_agenda_variables_io(self): if nout == 0 and allvars[var].type in opts: ok_anyways = True + # Allow agendas to be pure input if nin != 0 and allvars[var].type == "Agenda": ok_anyways = True diff --git a/src/m_obsel.cc b/src/m_obsel.cc index a008dc8d2b..de11e04b51 100644 --- a/src/m_obsel.cc +++ b/src/m_obsel.cc @@ -9,22 +9,28 @@ #include "debug.h" #include "matpack_data.h" -void measurement_sensorSimple(ArrayOfSensorObsel& measurement_sensor, - const AscendingGrid& frequency_grid, - const Vector3& pos, - const Vector2& los, - const Stokvec& pol) try { - measurement_sensor.resize(frequency_grid.size()); +void measurement_sensorInit(ArrayOfSensorObsel& measurement_sensor) { + measurement_sensor = {}; +} + +void measurement_sensorAddSimple(ArrayOfSensorObsel& measurement_sensor, + const AscendingGrid& frequency_grid, + const Vector3& pos, + const Vector2& los, + const Stokvec& pol) try { + const Index n = frequency_grid.size(); + const Size sz = measurement_sensor.size(); + + measurement_sensor.resize(sz + n); auto f = std::make_shared(frequency_grid); auto p = std::make_shared( SensorPosLosVector{{.pos = pos, .los = los}}); - for (Size i = 0; i < measurement_sensor.size(); i++) { - StokvecMatrix w(1, frequency_grid.size(), 0); - - w(0, i) = pol; - measurement_sensor[i] = SensorObsel(f, p, std::move(w)); + for (Index i = 0; i < n; i++) { + StokvecMatrix w(1, n, 0); + w(0, i) = pol; + measurement_sensor[i + sz] = {f, p, std::move(w)}; } } ARTS_METHOD_ERROR_CATCH @@ -33,17 +39,19 @@ Numeric gauss(Numeric f0, Numeric f, Numeric fwhm) { return std::exp(-4 * Constant::ln_2 * Math::pow2((f - f0) / fwhm)); } -void measurement_sensorSimpleGaussian(ArrayOfSensorObsel& measurement_sensor, - const AscendingGrid& frequency_grid, - const Vector& fwhm, - const Vector3& pos, - const Vector2& los, - const Stokvec& pol) try { +void measurement_sensorAddVectorGaussian(ArrayOfSensorObsel& measurement_sensor, + const AscendingGrid& frequency_grid, + const Vector& fwhm, + const Vector3& pos, + const Vector2& los, + const Stokvec& pol) try { const Index n = frequency_grid.size(); + const Size sz = measurement_sensor.size(); + ARTS_USER_ERROR_IF(n < 2, "Must have a frequency grid") ARTS_USER_ERROR_IF(n != fwhm.size(), "Must have a frequency grid") - measurement_sensor.resize(n); + measurement_sensor.resize(sz + n); auto f = std::make_shared(frequency_grid); auto p = std::make_shared( @@ -51,16 +59,16 @@ void measurement_sensorSimpleGaussian(ArrayOfSensorObsel& measurement_sensor, String error; #pragma omp parallel for if (not arts_omp_in_parallel()) - for (Size i = 0; i < measurement_sensor.size(); i++) { + for (Index i = 0; i < n; i++) { try { - StokvecMatrix w(1, frequency_grid.size(), pol); + StokvecMatrix w(1, n, pol); for (Index j = 0; j < n; j++) { w(0, j) *= gauss(frequency_grid[i], frequency_grid[j], fwhm[i]); } - measurement_sensor[i] = SensorObsel(f, p, std::move(w)); - measurement_sensor[i].normalize(pol, sum(pol)); + measurement_sensor[i + sz] = {f, p, std::move(w)}; + measurement_sensor[i + sz].normalize(pol, sum(pol)); } catch (std::runtime_error& e) { #pragma omp critical if (error.empty()) error = e.what(); @@ -71,16 +79,16 @@ void measurement_sensorSimpleGaussian(ArrayOfSensorObsel& measurement_sensor, } ARTS_METHOD_ERROR_CATCH -void measurement_sensorSimpleGaussian(ArrayOfSensorObsel& measurement_sensor, - const AscendingGrid& frequency_grid, - const Numeric& fwhm, - const Vector3& pos, - const Vector2& los, - const Stokvec& pol) { - measurement_sensorSimpleGaussian(measurement_sensor, - frequency_grid, - Vector(frequency_grid.size(), fwhm), - pos, - los, - pol); +void measurement_sensorAddSimpleGaussian(ArrayOfSensorObsel& measurement_sensor, + const AscendingGrid& frequency_grid, + const Numeric& fwhm, + const Vector3& pos, + const Vector2& los, + const Stokvec& pol) { + measurement_sensorAddVectorGaussian(measurement_sensor, + frequency_grid, + Vector(frequency_grid.size(), fwhm), + pos, + los, + pol); } diff --git a/src/workspace_meta_methods.cpp b/src/workspace_meta_methods.cpp index 303a6eb663..5c672a789c 100644 --- a/src/workspace_meta_methods.cpp +++ b/src/workspace_meta_methods.cpp @@ -13,6 +13,33 @@ std::vector internal_meta_methods_creator() { std::vector out; + out.push_back(WorkspaceMethodInternalMetaRecord{ + .name = "measurement_sensorSimple", + .desc = "Wrapper for a single simple dirac-opening sensor", + .author = {"Richard Larsson"}, + .methods = {"measurement_sensorInit", + "measurement_sensorAddSimple"}, + .out = {"measurement_sensor"}, + }); + + out.push_back(WorkspaceMethodInternalMetaRecord{ + .name = "measurement_sensorSimpleGaussian", + .desc = "Wrapper for a single simple Gaussian-opening sensor", + .author = {"Richard Larsson"}, + .methods = {"measurement_sensorInit", + "measurement_sensorAddSimpleGaussian"}, + .out = {"measurement_sensor"}, + }); + + out.push_back(WorkspaceMethodInternalMetaRecord{ + .name = "measurement_sensorVectorGaussian", + .desc = "Wrapper for a single simple Gaussian-opening sensor", + .author = {"Richard Larsson"}, + .methods = {"measurement_sensorInit", + "measurement_sensorAddVectorGaussian"}, + .out = {"measurement_sensor"}, + }); + out.push_back(WorkspaceMethodInternalMetaRecord{ .name = "disort_spectral_flux_fieldFromAgenda", .desc = "Use Disort for clearsky calculations of spectral flux field", diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index 9495a81551..d58b1faa1a 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -3116,15 +3116,23 @@ before consumption by, e.g., *OEM* .pass_workspace = true, }; - wsm_data["measurement_sensorSimple"] = { + wsm_data["measurement_sensorInit"] = { .desc = - R"--(Sets a sensor with a Gaussian channel opening around the frequency grid. + R"--(Initialize *measurement_sensor* to empty. +)--", + .author = {"Richard Larsson"}, + .out = {"measurement_sensor"}, + }; + + wsm_data["measurement_sensorAddSimple"] = { + .desc = + R"--(Adds a sensor with a dirac channel opening around the frequency grid. All elements share position, line-of-sight, and frequency grid. )--", .author = {"Richard Larsson"}, .out = {"measurement_sensor"}, - .in = {"frequency_grid"}, + .in = {"measurement_sensor", "frequency_grid"}, .gin = {"pos", "los", "pol"}, .gin_type = {"Vector3", "Vector2", "Stokvec"}, .gin_value = {std::nullopt, @@ -3136,9 +3144,9 @@ All elements share position, line-of-sight, and frequency grid. "The polarization whos dot-product with the spectral radiance becomes the measurement"}, }; - wsm_data["measurement_sensorSimpleGaussian"] = { + wsm_data["measurement_sensorAddSimpleGaussian"] = { .desc = - R"--(Sets a sensor with a Gaussian channel opening around the frequency grid. + R"--(Adds a sensor with a Gaussian channel opening around the frequency grid. All elements share position, line-of-sight, and frequency grid. @@ -3146,15 +3154,39 @@ Note that this means you only get "half" a Gaussian channel for the outermost ch )--", .author = {"Richard Larsson"}, .out = {"measurement_sensor"}, - .in = {"frequency_grid"}, + .in = {"measurement_sensor", "frequency_grid"}, + .gin = {"fwhm", "pos", "los", "pol"}, + .gin_type = {"Numeric", "Vector3", "Vector2", "Stokvec"}, + .gin_value = {std::nullopt, + std::nullopt, + std::nullopt, + rtepack::to_stokvec(PolarizationChoice::I)}, + .gin_desc = + {"The full-width half-maximum of the channel", + "A position [alt, lat, lon]", + "A line of sight [zenith, azimuth]", + "The polarization whos dot-product with the spectral radiance becomes the measurement"}, + }; + + wsm_data["measurement_sensorAddVectorGaussian"] = { + .desc = + R"--(Adds a sensor with a Gaussian channel opening around the frequency grid. + +All elements share position, line-of-sight, and frequency grid. + +Note that this means you only get "half" a Gaussian channel for the outermost channels. +)--", + .author = {"Richard Larsson"}, + .out = {"measurement_sensor"}, + .in = {"measurement_sensor", "frequency_grid"}, .gin = {"fwhm", "pos", "los", "pol"}, - .gin_type = {"Vector, Numeric", "Vector3", "Vector2", "Stokvec"}, + .gin_type = {"Vector", "Vector3", "Vector2", "Stokvec"}, .gin_value = {std::nullopt, std::nullopt, std::nullopt, rtepack::to_stokvec(PolarizationChoice::I)}, .gin_desc = - {"The full-width half-maximum of the channel(s)", + {"The full-width half-maximum of the channels", "A position [alt, lat, lon]", "A line of sight [zenith, azimuth]", "The polarization whos dot-product with the spectral radiance becomes the measurement"}, From f8c8913a17eed25b457b5c624e64137b61d20c3e Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Thu, 21 Nov 2024 09:44:56 +0100 Subject: [PATCH 4/7] tmp --- .../2.zeeman-sensor.py | 6 +- src/core/CMakeLists.txt | 1 - src/core/matpack/matpack_math_extra.h | 25 +- src/core/retrieval_target.cc | 29 --- src/core/retrieval_target.h | 55 ++--- src/core/rtepack/rtepack_stokes_vector.h | 39 ++- src/core/sensor/obsel.cpp | 29 ++- src/core/sensor/obsel.h | 15 +- src/m_jactargets.cc | 13 +- src/m_obsel.cc | 35 +-- src/m_rad.cc | 2 +- src/m_retrieval.cc | 130 +++++----- src/python_interface/py_retrieval.cpp | 13 +- src/workspace_groups.cpp | 6 + src/workspace_methods.cpp | 224 ++++++------------ src/xml_io_compound_types.cc | 77 +++++- tests/core/sensor/jac.py | 91 +++++++ 17 files changed, 421 insertions(+), 369 deletions(-) delete mode 100644 src/core/retrieval_target.cc diff --git a/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py b/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py index c1e793da22..34b4631177 100644 --- a/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py +++ b/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py @@ -37,11 +37,13 @@ ws.ray_path_observer_agendaSet(option="Geometric") ws.spectral_radiance_observer_agendaSet(option="Emission") -# %% Set up a sensor with Gaussian channel widths on individual frequency ranges +# %% Set up a sensor with Gaussian FWHM channel widths on individual frequency ranges pos = [100e3, 0, 0] los = [180.0, 0.0] -ws.measurement_sensorSimpleGaussian(fwhm=1e5, pos=pos, los=los, pol="RC") +ws.measurement_sensorSimpleGaussian( + std=1e5 / (2 * np.sqrt(2 * np.log(2))), pos=pos, los=los, pol="RC" +) # %% Core calculations diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 0bf2877ecd..0566d75da8 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -60,7 +60,6 @@ add_library(artscore STATIC radiation_field.cc raw.cc refraction.cc - retrieval_target.cc sensor.cc special_interp.cc sun.cc diff --git a/src/core/matpack/matpack_math_extra.h b/src/core/matpack/matpack_math_extra.h index 142881f51e..e8ddb323b0 100644 --- a/src/core/matpack/matpack_math_extra.h +++ b/src/core/matpack/matpack_math_extra.h @@ -9,28 +9,31 @@ all functionality should mimic the matpack functionality, but not all matpack functionality */ -#include "matpack_concepts.h" - #include +#include "matpack_concepts.h" + template -concept standard_iterable = not matpack::any_matpack_type and matpack::rank() == 1 and requires (T a) { - { a.begin() } -> std::random_access_iterator; - { a.end() } -> std::random_access_iterator; -}; +concept standard_iterable = not matpack::any_matpack_type and + matpack::rank() == 1 and requires(T a) { + { a.begin() } -> std::random_access_iterator; + { a.end() } -> std::random_access_iterator; + }; /** Find minimum of x by reduction */ -constexpr auto min(standard_iterable auto &&x) { +constexpr auto min(standard_iterable auto&& x) { return std::reduce( - x.begin(), x.end(), + x.begin(), + x.end(), std::numeric_limits>::max(), [](auto a, auto b) { return a < b ? a : b; }); } /** Find miximum of x by reduction */ -constexpr auto max(standard_iterable auto &&x) { +constexpr auto max(standard_iterable auto&& x) { return std::reduce( - x.begin(), x.end(), + x.begin(), + x.end(), std::numeric_limits>::lowest(), [](auto a, auto b) { return a > b ? a : b; }); } @@ -43,4 +46,4 @@ constexpr auto minmax(standard_iterable auto&& x) { /** Sum up the values */ constexpr auto sum(standard_iterable auto&& x) { return std::reduce(x.begin(), x.end()); -} +} \ No newline at end of file diff --git a/src/core/retrieval_target.cc b/src/core/retrieval_target.cc deleted file mode 100644 index 5907666bd3..0000000000 --- a/src/core/retrieval_target.cc +++ /dev/null @@ -1,29 +0,0 @@ -#include "retrieval_target.h" - -void JacobianTargetsDiagonalCovarianceMatrixMap::set( - const JacobianTargetType& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse) { - map[type] = {matrix, inverse}; -} - -void JacobianTargetsDiagonalCovarianceMatrixMap::set( - const AtmKeyVal& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse) { - set(JacobianTargetType{type}, matrix, inverse); -} - -void JacobianTargetsDiagonalCovarianceMatrixMap::set( - const SurfaceKeyVal& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse) { - set(JacobianTargetType{type}, matrix, inverse); -} - -void JacobianTargetsDiagonalCovarianceMatrixMap::set( - const LblLineKey& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse) { - set(JacobianTargetType{type}, matrix, inverse); -} diff --git a/src/core/retrieval_target.h b/src/core/retrieval_target.h index 28e04e74ee..b49a0f3a85 100644 --- a/src/core/retrieval_target.h +++ b/src/core/retrieval_target.h @@ -2,47 +2,32 @@ #include -#include "atm.h" #include "covariance_matrix.h" -#include "enums.h" -#include "jacobian.h" -#include "lbl_data.h" -#include "surf.h" +#include -struct JacobianTargetsDiagonalCovarianceMatrixMap { - std::unordered_map> - map; +struct PairOfBlockMatrix { + BlockMatrix first, second; +}; + +using JacobianTargetsDiagonalCovarianceMatrixMap = + std::unordered_map; - [[nodiscard]] auto begin() { return map.begin(); } - [[nodiscard]] auto end() { return map.end(); } - [[nodiscard]] auto begin() const { return map.begin(); } - [[nodiscard]] auto end() const { return map.end(); } - [[nodiscard]] auto cbegin() const { return map.cbegin(); } - [[nodiscard]] auto cend() const { return map.cend(); } - [[nodiscard]] auto size() const { return map.size(); } - [[nodiscard]] auto empty() const { return map.empty(); } - [[nodiscard]] auto clear() { return map.clear(); } - [[nodiscard]] auto find(const JacobianTargetType& x) const { return map.find(x); } +template <> +struct std::formatter { + format_tags tags; - void set(const JacobianTargetType& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse = BlockMatrix{}); + [[nodiscard]] constexpr auto& inner_fmt() { return *this; } + [[nodiscard]] constexpr auto& inner_fmt() const { return *this; } - void set(const AtmKeyVal& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse = BlockMatrix{}); - - void set(const SurfaceKeyVal& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse = BlockMatrix{}); - - void set(const LblLineKey& type, - const BlockMatrix& matrix, - const BlockMatrix& inverse = BlockMatrix{}); + constexpr std::format_parse_context::iterator parse( + std::format_parse_context& ctx) { + return parse_format_tags(tags, ctx); + } - friend std::ostream& operator<<( - std::ostream& os, const JacobianTargetsDiagonalCovarianceMatrixMap&) { - return os << "NotImplementedYet"; + template + FmtContext::iterator format(const PairOfBlockMatrix& x, + FmtContext& ctx) const { + return tags.format(ctx, x.first, tags.sep(), x.second); } }; diff --git a/src/core/rtepack/rtepack_stokes_vector.h b/src/core/rtepack/rtepack_stokes_vector.h index 672913a345..e93f57d9b2 100644 --- a/src/core/rtepack/rtepack_stokes_vector.h +++ b/src/core/rtepack/rtepack_stokes_vector.h @@ -51,31 +51,30 @@ struct stokvec final : vec4 { [[nodiscard]] constexpr bool is_zero() const { return I() == 0.0 && Q() == 0.0 && U() == 0.0 && V() == 0.0; } + + [[nodiscard]] constexpr bool is_polarized() const { + return Q() != 0.0 or U() != 0.0 or V() != 0.0; + } + + [[nodiscard]] constexpr Size nonzero_components() const { + return static_cast(I() != 0.0) + static_cast(Q() != 0.0) + + static_cast(U() != 0.0) + static_cast(V() != 0.0); + } }; constexpr stokvec to_stokvec(PolarizationChoice p) { using enum PolarizationChoice; switch (p) { - case I: - return {1.0, 0.0, 0.0, 0.0}; - case Q: - return {0.0, 1.0, 0.0, 0.0}; - case U: - return {0.0, 0.0, 1.0, 0.0}; - case V: - return {0.0, 0.0, 0.0, 1.0}; - case Iv: - return {1.0, 1.0, 0.0, 0.0}; - case Ih: - return {1.0, -1.0, 0.0, 0.0}; - case Ip45: - return {1.0, 0.0, 1.0, 0.0}; - case Im45: - return {1.0, 0.0, -1.0, 0.0}; - case Ilhc: - return {1.0, 0.0, 0.0, -1.0}; - case Irhc: - return {1.0, 0.0, 0.0, 1.0}; + case I: return {1.0, 0.0, 0.0, 0.0}; + case Q: return {0.0, 1.0, 0.0, 0.0}; + case U: return {0.0, 0.0, 1.0, 0.0}; + case V: return {0.0, 0.0, 0.0, 1.0}; + case Iv: return {1.0, 1.0, 0.0, 0.0}; + case Ih: return {1.0, -1.0, 0.0, 0.0}; + case Ip45: return {1.0, 0.0, 1.0, 0.0}; + case Im45: return {1.0, 0.0, -1.0, 0.0}; + case Ilhc: return {1.0, 0.0, 0.0, -1.0}; + case Irhc: return {1.0, 0.0, 0.0, 1.0}; } std::unreachable(); } diff --git a/src/core/sensor/obsel.cpp b/src/core/sensor/obsel.cpp index 8bc3c65afe..063a6ea9f7 100644 --- a/src/core/sensor/obsel.cpp +++ b/src/core/sensor/obsel.cpp @@ -43,17 +43,22 @@ void Obsel::check() const { f->size()); } -void Obsel::normalize(Stokvec pol, Numeric new_value) { - const Numeric sum = std::transform_reduce( - w.elem_begin(), - w.elem_end(), - 0.0, - std::plus<>(), - [pol](const Stokvec& ws) -> Numeric { return pol * ws; }); - - ARTS_USER_ERROR_IF(sum == 0.0, "Cannot normalize, sum is zero"); +void Obsel::normalize(Stokvec pol) { + const Stokvec x = sum(w); + + if (x.I() != 0.0) pol.I() = std::abs(pol.I() / x.I()); + if (x.is_polarized()) { + const Numeric hyp = std::hypot(x.Q(), x.U(), x.V()); + pol.Q() = std::abs(pol.Q() / hyp); + pol.U() = std::abs(pol.U() / hyp); + pol.V() = std::abs(pol.V() / hyp); + } - w *= new_value / sum; + std::transform( + w.elem_begin(), w.elem_end(), w.elem_begin(), [pol](auto& e) -> Stokvec { + return { + e.I() * pol.I(), e.Q() * pol.Q(), e.U() * pol.U(), e.V() * pol.V()}; + }); } Numeric Obsel::sumup(const StokvecVectorView& i, Index ip) const { @@ -163,7 +168,7 @@ Index Obsel::find(const Vector3& pos, const Vector2& los) const { const auto first = poslos->cbegin(); const auto last = poslos->cend(); const auto same_poslos = [&](const auto& p) { - return pos == p.pos and los == p.los; + return &pos == &p.pos and & los == &p.los; }; const auto it = std::find_if(first, last, same_poslos); @@ -172,7 +177,7 @@ Index Obsel::find(const Vector3& pos, const Vector2& los) const { } Index Obsel::find(const AscendingGrid& frequency_grid) const { - return (*f != frequency_grid) ? dont_have : 0; + return (f.get() != &frequency_grid) ? dont_have : 0; } } // namespace sensor diff --git a/src/core/sensor/obsel.h b/src/core/sensor/obsel.h index 6d8e87d746..c1a65e49cb 100644 --- a/src/core/sensor/obsel.h +++ b/src/core/sensor/obsel.h @@ -127,9 +127,8 @@ class Obsel { * polarization will be equal to the new value. * * @param pol The polarization that is sampled for the sum. The default is [1, 0, 0, 0]. - * @param new_value The new value for the normalization. The default is 1.0. */ - void normalize(Stokvec pol = {1., 0., 0., 0.}, Numeric new_value = 1.0); + void normalize(Stokvec pol = {1., 0., 0., 0.}); [[nodiscard]] Numeric sumup(const StokvecVectorView& i, Index ip) const; void sumup(VectorView out, const StokvecMatrixView& j, Index ip) const; @@ -146,11 +145,11 @@ std::ostream& operator<<(std::ostream& os, const Array& obsel); } // namespace sensor struct SensorKey { - SensorKeyType type; + SensorKeyType type{}; - Index elem; + Index elem{}; - SensorJacobianModelType model; + SensorJacobianModelType model{}; Index polyorder{-1}; @@ -163,8 +162,10 @@ template <> struct std::hash { std::size_t operator()(const SensorKey& g) const { return std::hash{}(g.type) ^ - (std::hash{}(g.model) << 10) ^ - (std::hash{}(g.elem) << 20); + (std::hash{}(g.model) + << (8 * sizeof(SensorKeyType))) ^ + (std::hash{}(g.elem) + << (8 * (sizeof(SensorKeyType) + sizeof(SensorJacobianModelType)))); } }; diff --git a/src/m_jactargets.cc b/src/m_jactargets.cc index f8c04c9654..4a45995f27 100644 --- a/src/m_jactargets.cc +++ b/src/m_jactargets.cc @@ -156,10 +156,9 @@ void jacobian_targetsAddSensor(JacobianTargets& jacobian_targets, }); } -void jacobian_targetsAddSensorPolyFit( +void jacobian_targetsAddSensorFrequencyPolyFit( JacobianTargets& jacobian_targets, const ArrayOfSensorObsel& measurement_sensor, - const SensorKeyType& key, const Numeric& d, const Index& elem, const Index& polyorder) { @@ -170,11 +169,11 @@ void jacobian_targetsAddSensorPolyFit( polyorder < 0, "Polyorder must be non-negative: {}", polyorder) Jacobian::SensorTarget target{ - .type = {.type = key, - .elem = elem, - .model = SensorJacobianModelType::PolynomialOffset, - .polyorder = polyorder}, - // .original_grid = sensor_grid(measurement_sensor[elem], key)}, + .type = {.type = SensorKeyType::f, + .elem = elem, + .model = SensorJacobianModelType::PolynomialOffset, + .polyorder = polyorder, + .original_grid = measurement_sensor[elem].flat(SensorKeyType::f)}, .d = d, .target_pos = jacobian_targets.target_count(), }; diff --git a/src/m_obsel.cc b/src/m_obsel.cc index de11e04b51..0cb986e54f 100644 --- a/src/m_obsel.cc +++ b/src/m_obsel.cc @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -35,21 +36,23 @@ void measurement_sensorAddSimple(ArrayOfSensorObsel& measurement_sensor, } ARTS_METHOD_ERROR_CATCH -Numeric gauss(Numeric f0, Numeric f, Numeric fwhm) { - return std::exp(-4 * Constant::ln_2 * Math::pow2((f - f0) / fwhm)); -} - void measurement_sensorAddVectorGaussian(ArrayOfSensorObsel& measurement_sensor, const AscendingGrid& frequency_grid, - const Vector& fwhm, + const Vector& stds, const Vector3& pos, const Vector2& los, const Stokvec& pol) try { - const Index n = frequency_grid.size(); - const Size sz = measurement_sensor.size(); + using gauss = boost::math::normal_distribution; + using boost::math::pdf; + + const Index n = frequency_grid.size(); + const Size sz = measurement_sensor.size(); + const Size nonzero = pol.nonzero_components(); ARTS_USER_ERROR_IF(n < 2, "Must have a frequency grid") - ARTS_USER_ERROR_IF(n != fwhm.size(), "Must have a frequency grid") + ARTS_USER_ERROR_IF(n != stds.size(), + "Must have a standard deviation for each frequency point") + ARTS_USER_ERROR_IF(nonzero == 0, "pol is 0") measurement_sensor.resize(sz + n); @@ -63,12 +66,13 @@ void measurement_sensorAddVectorGaussian(ArrayOfSensorObsel& measurement_sensor, try { StokvecMatrix w(1, n, pol); + const gauss dist(frequency_grid[i], stds[i]); for (Index j = 0; j < n; j++) { - w(0, j) *= gauss(frequency_grid[i], frequency_grid[j], fwhm[i]); + w(0, j) *= pdf(dist, frequency_grid[j]); } measurement_sensor[i + sz] = {f, p, std::move(w)}; - measurement_sensor[i + sz].normalize(pol, sum(pol)); + measurement_sensor[i + sz].normalize(pol); } catch (std::runtime_error& e) { #pragma omp critical if (error.empty()) error = e.what(); @@ -81,14 +85,11 @@ ARTS_METHOD_ERROR_CATCH void measurement_sensorAddSimpleGaussian(ArrayOfSensorObsel& measurement_sensor, const AscendingGrid& frequency_grid, - const Numeric& fwhm, + const Numeric& std, const Vector3& pos, const Vector2& los, const Stokvec& pol) { - measurement_sensorAddVectorGaussian(measurement_sensor, - frequency_grid, - Vector(frequency_grid.size(), fwhm), - pos, - los, - pol); + const Vector stds(frequency_grid.size(), std); + measurement_sensorAddVectorGaussian( + measurement_sensor, frequency_grid, stds, pos, los, pol); } diff --git a/src/m_rad.cc b/src/m_rad.cc index 3a6f337d1a..e848f7055e 100644 --- a/src/m_rad.cc +++ b/src/m_rad.cc @@ -396,7 +396,7 @@ f_grid_ptr->size() = {} spectral_radiance_jacobian.shape() = {:B,}, f_grid_ptr->size() = {}, -measurement_jacobian.ncols() = {} +measurement_jacobian.ncols() = {} )", spectral_radiance_jacobian.shape(), f_grid_ptr->size(), diff --git a/src/m_retrieval.cc b/src/m_retrieval.cc index 11046d2070..3377821efd 100644 --- a/src/m_retrieval.cc +++ b/src/m_retrieval.cc @@ -20,97 +20,96 @@ void RetrievalAddSurface(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const SurfaceKey& key, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddSurface(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = SurfaceKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.surf().back().type}] = {matrix, inverse}; } void RetrievalAddSurface(JacobianTargets& jacobian_targets, - JacobianTargetsDiagonalCovarianceMatrixMap& - covariance_matrix_diagonal_blocks, - const SurfaceTypeTag& key, - const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + JacobianTargetsDiagonalCovarianceMatrixMap& + covariance_matrix_diagonal_blocks, + const SurfaceTypeTag& key, + const Numeric& d, + const BlockMatrix& matrix, + const BlockMatrix& inverse) { jacobian_targetsAddSurface(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = SurfaceKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.surf().back().type}] = {matrix, inverse}; } void RetrievalAddSurface(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const SurfacePropertyTag& key, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddSurface(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = SurfaceKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.surf().back().type}] = {matrix, inverse}; } void RetrievalAddAtmosphere(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const AtmKey& key, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddAtmosphere(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddAtmosphere(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const SpeciesEnum& key, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddAtmosphere(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddAtmosphere(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const SpeciesIsotope& key, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddAtmosphere(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddAtmosphere(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const QuantumIdentifier& key, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddAtmosphere(jacobian_targets, key, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{key}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } - void RetrievalAddSpeciesVMR(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const SpeciesEnum& species, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddSpeciesVMR(jacobian_targets, species, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{species}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddSpeciesIsotopologueRatio( @@ -118,56 +117,73 @@ void RetrievalAddSpeciesIsotopologueRatio( JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const SpeciesIsotope& species, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddSpeciesIsotopologueRatio(jacobian_targets, species, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{species}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddMagneticField(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const String& component, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddMagneticField(jacobian_targets, component, d); - covariance_matrix_diagonal_blocks.set( - {.target = to_mag(component)}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddWindField(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, const String& component, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddWindField(jacobian_targets, component, d); - covariance_matrix_diagonal_blocks.set( - {.target = to_wind(component)}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddTemperature(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddTemperature(jacobian_targets, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{AtmKey::t}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; } void RetrievalAddPressure(JacobianTargets& jacobian_targets, JacobianTargetsDiagonalCovarianceMatrixMap& covariance_matrix_diagonal_blocks, + const Numeric& d, const BlockMatrix& matrix, - const BlockMatrix& inverse, - const Numeric& d) { + const BlockMatrix& inverse) { jacobian_targetsAddPressure(jacobian_targets, d); - covariance_matrix_diagonal_blocks.set( - {.target = AtmKeyVal{AtmKey::p}}, matrix, inverse); + covariance_matrix_diagonal_blocks[JacobianTargetType{ + jacobian_targets.atm().back().type}] = {matrix, inverse}; +} + +void RetrievalAddSensorFrequencyPolyFit( + JacobianTargets& jacobian_targets, + JacobianTargetsDiagonalCovarianceMatrixMap& + covariance_matrix_diagonal_blocks, + const ArrayOfSensorObsel& measurement_sensor, + const Numeric& d, + const Index& elem, + const Index& polyorder, + const BlockMatrix& matrix, + const BlockMatrix& inverse) { + jacobian_targetsAddSensorFrequencyPolyFit( + jacobian_targets, measurement_sensor, d, elem, polyorder); + auto keyk = JacobianTargetType{ + jacobian_targets.sensor().back().type}; + covariance_matrix_diagonal_blocks[keyk] = {matrix, inverse}; } diff --git a/src/python_interface/py_retrieval.cpp b/src/python_interface/py_retrieval.cpp index b8e0613848..aa2fef8eb7 100644 --- a/src/python_interface/py_retrieval.cpp +++ b/src/python_interface/py_retrieval.cpp @@ -1,12 +1,21 @@ +#include +#include +#include + #include "hpy_arts.h" #include "python_interface.h" -#include namespace Python { void py_retrieval(py::module_& m) try { - py::class_ jtdcmm( + auto jtdcmm = py::bind_map( m, "JacobianTargetsDiagonalCovarianceMatrixMap"); workspace_group_interface(jtdcmm); + + py::class_ pobm(m, "PairOfBlockMatrix"); + workspace_group_interface(pobm); + pobm.def_rw("first", &PairOfBlockMatrix::first, "Matrix"); + pobm.def_rw("second", &PairOfBlockMatrix::second, "Inverse of Matrix"); py::class_ jtt(m, "JacobianTargetType"); jtt.def_rw("value", &JacobianTargetType::target, "Target"); diff --git a/src/workspace_groups.cpp b/src/workspace_groups.cpp index 1e1fed0206..09ff1aaf58 100644 --- a/src/workspace_groups.cpp +++ b/src/workspace_groups.cpp @@ -847,6 +847,11 @@ when computing the Jacobian matrix or partial derivatives. )--", }; + wsg_data["PairOfBlockMatrix"] = { + .file = "retrieval_target.h", + .desc = R"--(A pair of *BlockMatrix* objects)--", + }; + wsg_data["JacobianTargetsDiagonalCovarianceMatrixMap"] = { .file = "retrieval_target.h", .desc = @@ -855,6 +860,7 @@ when computing the Jacobian matrix or partial derivatives. The intended use of this type is to store required *BlockMatrix* objects so that the user-interface for setting up retrieval targets can be simplified. )--", + .map_type = true, }; wsg_data["PropagationPathPoint"] = { diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index d58b1faa1a..2e470f94bb 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -1993,6 +1993,56 @@ building of an actual Jacobian matrix. "measurement_sensor"}, }; + const auto jac2ret = [&wsm_data](const std::string& name) { + auto v = wsm_data[name]; + v.desc += std::format(R"( +This method wraps *{}* together with adding the covariance matrices, +to the *covariance_matrix_diagonal_blocks*, which are required to perform *OEM*. + +The input covariance matrices must fit the size of the later computed model state +represented by the *jacobian_targets*. The covariance matrix inverse +)", name); + v.out.insert(v.out.begin() + 1, "covariance_matrix_diagonal_blocks"); + v.in.insert(v.in.begin() + 1, "covariance_matrix_diagonal_blocks"); + v.gin.insert(v.gin.end(), "matrix"); + v.gin.insert(v.gin.end(), "inverse"); + v.gin_type.insert(v.gin_type.end(), "BlockMatrix"); + v.gin_type.insert(v.gin_type.end(), "BlockMatrix"); + v.gin_value.insert(v.gin_value.end(), std::nullopt); + v.gin_value.insert(v.gin_value.end(), BlockMatrix{}); + v.gin_desc.insert(v.gin_desc.end(), "The covariance diagonal block matrix"); + v.gin_desc.insert(v.gin_desc.end(), + "The inverse covariance diagonal block matrix"); + return v; + }; + + wsm_data["jacobian_targetsAddSensorFrequencyPolyFit"] = { + .desc = + R"--(Set sensor frequency derivative to use polynomial fitting offset + +Order 0 means constant: f := f0 + a +Order 1 means linear: f := f0 + a + b * f0 +and so on. The derivatives that are added to the *model_state_vector* are +those with regards to a, b, etc.. + +Note that sensor elements that share the same frequency grid will all +be affected by this method. In fact, finalizing the *jacobian_targets* +with several elements of the frequency grids is a runtime error. +)--", + .author = {"Richard Larsson"}, + .out = {"jacobian_targets"}, + .in = {"jacobian_targets", "measurement_sensor"}, + .gin = {"d", "elem", "polyorder"}, + .gin_type = {"Numeric", "Index", "Index"}, + .gin_value = {Numeric{0.1}, std::nullopt, Index{0}}, + .gin_desc = + {"The perturbation used in methods that cannot compute derivatives analytically", + "The sensor element whose frequency grid to use", + "The order of the polynomial fit"}, + }; + wsm_data["RetrievalAddSensorFrequencyPolyFit"] = + jac2ret("jacobian_targetsAddSensorFrequencyPolyFit"); + wsm_data["jacobian_targetsAddTemperature"] = { .desc = R"--(Set temperature derivative )--", @@ -2005,6 +2055,8 @@ building of an actual Jacobian matrix. .gin_desc = {"The perturbation used in methods that cannot compute derivatives analytically"}, }; + wsm_data["RetrievalAddTemperature"] = + jac2ret("jacobian_targetsAddTemperature"); wsm_data["jacobian_targetsAddPressure"] = { .desc = R"--(Set pressure derivative @@ -2018,6 +2070,7 @@ building of an actual Jacobian matrix. .gin_desc = {"The perturbation used in methods that cannot compute derivatives analytically"}, }; + wsm_data["RetrievalAddPressure"] = jac2ret("jacobian_targetsAddPressure"); wsm_data["jacobian_targetsAddMagneticField"] = { .desc = R"--(Set magnetic field derivative @@ -2034,6 +2087,7 @@ See *FieldComponent* for valid ``component`` {"The component to use [u, v, w]", "The perturbation used in methods that cannot compute derivatives analytically"}, }; + wsm_data["RetrievalAddMagneticField"] = jac2ret("jacobian_targetsAddMagneticField"); wsm_data["jacobian_targetsAddWindField"] = { .desc = R"--(Set wind field derivative @@ -2053,6 +2107,7 @@ See *FieldComponent* for valid ``component`` {"The component to use [u, v, w]", "The perturbation used in methods that cannot compute derivatives analytically"}, }; +wsm_data["RetrievalAddWindField"] = jac2ret("jacobian_targetsAddWindField"); wsm_data["jacobian_targetsAddSpeciesVMR"] = { .desc = R"--(Set volume mixing ratio derivative @@ -2069,6 +2124,7 @@ See *SpeciesEnum* for valid ``species`` {"The species of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; +wsm_data["RetrievalAddSpeciesVMR"] = jac2ret("jacobian_targetsAddSpeciesVMR"); wsm_data["jacobian_targetsAddAtmosphere"] = { .desc = R"--(Sets an atmospheric target @@ -2084,6 +2140,7 @@ See *SpeciesEnum* for valid ``species`` {"The target of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; +wsm_data["RetrievalAddAtmosphere"] = jac2ret("jacobian_targetsAddAtmosphere"); wsm_data["jacobian_targetsAddSurface"] = { .desc = R"--(Sets a surface target @@ -2098,6 +2155,7 @@ See *SpeciesEnum* for valid ``species`` {"The target of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; +wsm_data["RetrievalAddSurface"] = jac2ret("jacobian_targetsAddSurface"); wsm_data["jacobian_targetsAddSpeciesIsotopologueRatio"] = { .desc = R"--(Set isotopologue ratio derivative @@ -2114,6 +2172,7 @@ See *SpeciesIsotope* for valid ``species`` {"The species isotopologue of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; +wsm_data["RetrievalAddSpeciesIsotopologueRatio"] = jac2ret("jacobian_targetsAddSpeciesIsotopologueRatio"); wsm_data["absorption_bandsReadHITRAN"] = { .desc = @@ -3151,46 +3210,29 @@ All elements share position, line-of-sight, and frequency grid. All elements share position, line-of-sight, and frequency grid. Note that this means you only get "half" a Gaussian channel for the outermost channels. + +The I component's distribution is normalized to 1 or 0 by itself, while +the Q, U, and V components' hypotenuse are normalized to 1 or 0 together. )--", .author = {"Richard Larsson"}, .out = {"measurement_sensor"}, .in = {"measurement_sensor", "frequency_grid"}, - .gin = {"fwhm", "pos", "los", "pol"}, + .gin = {"std", "pos", "los", "pol"}, .gin_type = {"Numeric", "Vector3", "Vector2", "Stokvec"}, .gin_value = {std::nullopt, std::nullopt, std::nullopt, rtepack::to_stokvec(PolarizationChoice::I)}, .gin_desc = - {"The full-width half-maximum of the channel", + {"The standard deviations of the channels", "A position [alt, lat, lon]", "A line of sight [zenith, azimuth]", "The polarization whos dot-product with the spectral radiance becomes the measurement"}, }; - wsm_data["measurement_sensorAddVectorGaussian"] = { - .desc = - R"--(Adds a sensor with a Gaussian channel opening around the frequency grid. - -All elements share position, line-of-sight, and frequency grid. - -Note that this means you only get "half" a Gaussian channel for the outermost channels. -)--", - .author = {"Richard Larsson"}, - .out = {"measurement_sensor"}, - .in = {"measurement_sensor", "frequency_grid"}, - .gin = {"fwhm", "pos", "los", "pol"}, - .gin_type = {"Vector", "Vector3", "Vector2", "Stokvec"}, - .gin_value = {std::nullopt, - std::nullopt, - std::nullopt, - rtepack::to_stokvec(PolarizationChoice::I)}, - .gin_desc = - {"The full-width half-maximum of the channels", - "A position [alt, lat, lon]", - "A line of sight [zenith, azimuth]", - "The polarization whos dot-product with the spectral radiance becomes the measurement"}, - }; + wsm_data["measurement_sensorAddVectorGaussian"] = + wsm_data["measurement_sensorAddSimpleGaussian"]; + wsm_data["measurement_sensorAddVectorGaussian"].gin_type[0] = "Vector"; wsm_data["sun_pathFromObserverAgenda"] = { .desc = @@ -3952,138 +3994,6 @@ Note that you must have set the optical thickness before calling this. "covariance_matrix_diagonal_blocks"}, }; - wsm_data["RetrievalAddSurface"] = { - .desc = R"(Add surface property the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"species", "matrix", "inverse", "d"}, - .gin_type = {"SurfaceKey,SurfaceTypeTag,SurfacePropertyTag", - "BlockMatrix", - "BlockMatrix", - "Numeric"}, - .gin_value = {std::nullopt, std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The property added to the retrieval system", - "The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddAtmosphere"] = { - .desc = R"(Add atmospheric property the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"species", "matrix", "inverse", "d"}, - .gin_type = {"AtmKey,SpeciesEnum,SpeciesIsotope,QuantumIdentifier", - "BlockMatrix", - "BlockMatrix", - "Numeric"}, - .gin_value = {std::nullopt, std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The property added to the retrieval system", - "The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddSpeciesVMR"] = { - .desc = R"(Add species VMR to the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"species", "matrix", "inverse", "d"}, - .gin_type = {"SpeciesEnum", "BlockMatrix", "BlockMatrix", "Numeric"}, - .gin_value = {std::nullopt, std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The species added to the retrieval system", - "The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddSpeciesIsotopologueRatio"] = { - .desc = R"(Add species isotopologue ratio to the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"species", "matrix", "inverse", "d"}, - .gin_type = {"SpeciesIsotope", "BlockMatrix", "BlockMatrix", "Numeric"}, - .gin_value = {std::nullopt, std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The isotopologue ratio added to the retrieval system", - "The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddMagneticField"] = { - .desc = R"(Add magnetic field component to the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"component", "matrix", "inverse", "d"}, - .gin_type = {"String", "BlockMatrix", "BlockMatrix", "Numeric"}, - .gin_value = {std::nullopt, std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The magnetic field component added to the retrieval system", - "The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddWindField"] = { - .desc = R"(Add species VMR to the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"component", "matrix", "inverse", "d"}, - .gin_type = {"String", "BlockMatrix", "BlockMatrix", "Numeric"}, - .gin_value = {std::nullopt, std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The magnetic field component added to the retrieval system", - "The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddTemperature"] = { - .desc = R"(Add atmospheric temperature to the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"matrix", "inverse", "d"}, - .gin_type = {"BlockMatrix", "BlockMatrix", "Numeric"}, - .gin_value = {std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - - wsm_data["RetrievalAddPressure"] = { - .desc = R"(Add atmospheric pressure to the retrieval setup. -)", - .author = {"Richard Larsson"}, - .out = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .in = {"jacobian_targets", "covariance_matrix_diagonal_blocks"}, - .gin = {"matrix", "inverse", "d"}, - .gin_type = {"BlockMatrix", "BlockMatrix", "Numeric"}, - .gin_value = {std::nullopt, BlockMatrix{}, Numeric{0.1}}, - .gin_desc = - {"The covariance diagonal block matrix", - "The inverse covariance diagonal block matrix", - "The delta of the value incase manual perturbation is needed"}, - }; - wsm_data["RetrievalFinalizeDiagonal"] = { .desc = R"(Finalize the retrieval setup. diff --git a/src/xml_io_compound_types.cc b/src/xml_io_compound_types.cc index 69b77de726..2bf834dfea 100644 --- a/src/xml_io_compound_types.cc +++ b/src/xml_io_compound_types.cc @@ -3078,7 +3078,13 @@ void xml_read_from_stream(std::istream& is_xml, ArtsXMLTag tag; tag.read_from_stream(is_xml); tag.check_name("SensorKey"); - + + xml_read_from_stream(is_xml, key.type, pbifs); + xml_read_from_stream(is_xml, key.elem, pbifs); + xml_read_from_stream(is_xml, key.model, pbifs); + xml_read_from_stream(is_xml, key.polyorder, pbifs); + xml_read_from_stream(is_xml, key.original_grid, pbifs); + tag.read_from_stream(is_xml); tag.check_name("/SensorKey"); } @@ -3098,6 +3104,14 @@ void xml_write_to_stream(std::ostream& os_xml, ArtsXMLTag close_tag; open_tag.set_name("SensorKey"); + open_tag.write_to_stream(os_xml); + os_xml << '\n'; + + xml_write_to_stream(os_xml, key.type, pbofs, "type"); + xml_write_to_stream(os_xml, key.elem, pbofs, "elem"); + xml_write_to_stream(os_xml, key.model, pbofs, "model"); + xml_write_to_stream(os_xml, key.polyorder, pbofs, "polyorder"); + xml_write_to_stream(os_xml, key.original_grid, pbofs, "original_grid"); close_tag.set_name("/SensorKey"); close_tag.write_to_stream(os_xml); @@ -3169,6 +3183,53 @@ void xml_write_to_stream(std::ostream& os_xml, close_tag.write_to_stream(os_xml); } +//=== PairOfBlockMatrix ========================================================= + +//! Reads PairOfBlockMatrix from XML input stream +/*! + \param is_xml XML Input stream + \param x PairOfBlockMatrix + \param pbifs Pointer to binary file stream. NULL for ASCII output. +*/ +void xml_read_from_stream(std::istream& is_xml, + PairOfBlockMatrix& x, + bifstream* pbifs) { + ArtsXMLTag tag; + tag.read_from_stream(is_xml); + tag.check_name("PairOfBlockMatrix"); + + xml_read_from_stream(is_xml, x.first, pbifs); + xml_read_from_stream(is_xml, x.second, pbifs); + + tag.read_from_stream(is_xml); + tag.check_name("/PairOfBlockMatrix"); +} + +//! Write PairOfBlockMatrix to XML output stream +/*! + \param os_xml XML output stream + \param x PairOfBlockMatrix + \param pbofs Pointer to binary file stream. NULL for ASCII output. + \param name +*/ +void xml_write_to_stream(std::ostream& os_xml, + const PairOfBlockMatrix& x, + bofstream* pbofs, + const String&) { + ArtsXMLTag open_tag; + ArtsXMLTag close_tag; + + open_tag.set_name("PairOfBlockMatrix"); + open_tag.write_to_stream(os_xml); + os_xml << '\n'; + + xml_write_to_stream(os_xml, x.first, pbofs, "first"); + xml_write_to_stream(os_xml, x.second, pbofs, "second"); + + close_tag.set_name("/PairOfBlockMatrix"); + close_tag.write_to_stream(os_xml); +} + //=== JacobianTargetsDiagonalCovarianceMatrixMap ========================================================= //! Reads JacobianTargetsDiagonalCovarianceMatrixMap from XML input stream @@ -3180,7 +3241,7 @@ void xml_write_to_stream(std::ostream& os_xml, void xml_read_from_stream(std::istream& is_xml, JacobianTargetsDiagonalCovarianceMatrixMap& jtmap, bifstream* pbifs) { - jtmap.map.clear(); + jtmap.clear(); ArtsXMLTag tag; tag.read_from_stream(is_xml); @@ -3191,14 +3252,9 @@ void xml_read_from_stream(std::istream& is_xml, for (Index i = 0; i < size; i++) { JacobianTargetType key; - BlockMatrix matrix; - BlockMatrix inverse; xml_read_from_stream(is_xml, key, pbifs); - xml_read_from_stream(is_xml, matrix, pbifs); - xml_read_from_stream(is_xml, inverse, pbifs); - - jtmap.set(key, matrix, inverse); + xml_read_from_stream(is_xml, jtmap[key], pbifs); } tag.read_from_stream(is_xml); @@ -3221,14 +3277,13 @@ void xml_write_to_stream( ArtsXMLTag close_tag; open_tag.set_name("JacobianTargetsDiagonalCovarianceMatrixMap"); - open_tag.add_attribute("size", static_cast(jtmap.map.size())); + open_tag.add_attribute("size", static_cast(jtmap.size())); open_tag.write_to_stream(os_xml); os_xml << '\n'; for (const auto& [key, value] : jtmap) { xml_write_to_stream(os_xml, key, pbofs, "key"); - xml_write_to_stream(os_xml, value.first, pbofs, "matrix"); - xml_write_to_stream(os_xml, value.second, pbofs, "inverse"); + xml_write_to_stream(os_xml, value, pbofs, "matrix-pair"); } close_tag.set_name("/JacobianTargetsDiagonalCovarianceMatrixMap"); diff --git a/tests/core/sensor/jac.py b/tests/core/sensor/jac.py index 03144efd69..7abb537673 100644 --- a/tests/core/sensor/jac.py +++ b/tests/core/sensor/jac.py @@ -2,3 +2,94 @@ import numpy as np import matplotlib.pyplot as plt +ws = pyarts.workspace.Workspace() + +# %% Sampled frequency range + +line_f0 = 118750348044.712 +f = np.linspace(-5e6, 5e6, 1001) + line_f0 + +# %% Species and line absorption + +ws.absorption_speciesSet(species=["O2-66"]) +ws.ReadCatalogData() +ws.absorption_bandsSelectFrequency(fmin=40e9, fmax=120e9, by_line=1) +ws.absorption_bandsSetZeeman(species="O2-66", fmin=118e9, fmax=119e9) +ws.WignerInit() + +# %% Use the automatic agenda setter for propagation matrix calculations +ws.propagation_matrix_agendaAuto() + +# %% Grids and planet + +ws.surface_fieldSetPlanetEllipsoid(option="Earth") +ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0 +ws.atmospheric_fieldRead( + toa=100e3, basename="planets/Earth/afgl/tropical/", missing_is_zero=1 +) +ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37") + +# %% Checks and settings + +ws.spectral_radiance_unit = "Tb" +ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground") +ws.spectral_radiance_surface_agendaSet(option="Blackbody") +ws.ray_path_observer_agendaSet(option="Geometric") +ws.spectral_radiance_observer_agendaSet(option="Emission") + +# %% Set up a sensor with Gaussian channel widths on individual frequency ranges + +pos = [100e3, 0, 0] +los = [180.0, 0.0] + +ws.measurement_sensorSimpleGaussian( + frequency_grid=f, std=1e5, pos=pos, los=los, pol="Ih" +) +ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f, std=1e5, pos=pos, los=los, pol="Iv" +) +ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f, std=1e5, pos=pos, los=los, pol="RC" +) +ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f, std=1e5, pos=pos, los=los, pol="LC" +) + +# %% Original calculations + +ws.measurement_vectorFromSensor() +orig = ws.measurement_vector * 1.0 + +# %% Modify the sensor + +DF1 = 1e6 +DF2 = 2e6 +DF3 = 3e6 +DF4 = -4e6 + +ws.measurement_sensorSimpleGaussian( + frequency_grid=f + DF1, std=1e5, pos=pos, los=los, pol="Ih" +) +ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f + DF2, std=1e5, pos=pos, los=los, pol="Iv" +) +ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f + DF3, std=1e5, pos=pos, los=los, pol="RC" +) +ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f + DF4, std=1e5, pos=pos, los=los, pol="LC" +) + +ws.measurement_vectorFromSensor() +mod = ws.measurement_vector * 1.0 + +# %% Set up the retrieval + +ws.RetrievalInit() +ws.RetrievalAddSensorFrequencyPolyFit(elem=0, d=1e3, matrix=np.diag(np.ones((1)) * 1e6)) +ws.RetrievalFinalizeDiagonal() + +ws.measurement_vector_fitted = [] +ws.model_state_vector = [] +ws.measurement_jacobian = [[]] +# %% From 207b148c2cf011f074fb3c65f6fd404a298476b6 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Thu, 21 Nov 2024 15:06:11 +0100 Subject: [PATCH 5/7] Retrieve multiple sensor's frequency grids --- .../2.zeeman-sensor.py | 28 +++-- src/core/jacobian/jacobian.cc | 54 ++++++++- src/core/jacobian/jacobian.h | 25 ++++- src/core/matpack/matpack_view.h | 2 +- src/core/minimize.cc | 72 ++++++------ src/core/minimize.h | 14 +-- src/m_jactargets.cc | 106 +++++------------- src/m_model_state.cc | 16 +++ src/m_rad.cc | 29 ++++- src/workspace_groups.cpp | 5 +- src/workspace_meta_methods.cpp | 22 ++-- src/workspace_methods.cpp | 36 ++++-- tests/core/sensor/jac.py | 99 +++++++++++----- 13 files changed, 310 insertions(+), 198 deletions(-) diff --git a/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py b/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py index 34b4631177..c1a0759172 100644 --- a/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py +++ b/examples/getting-started/2-clearsky-radiative-transfer/2.zeeman-sensor.py @@ -37,13 +37,11 @@ ws.ray_path_observer_agendaSet(option="Geometric") ws.spectral_radiance_observer_agendaSet(option="Emission") -# %% Set up a sensor with Gaussian FWHM channel widths on individual frequency ranges +# %% Set up a sensor with Gaussian standard deviation channel widths on individual frequency ranges pos = [100e3, 0, 0] los = [180.0, 0.0] -ws.measurement_sensorSimpleGaussian( - std=1e5 / (2 * np.sqrt(2 * np.log(2))), pos=pos, los=los, pol="RC" -) +ws.measurement_sensorSimpleGaussian(std=1e5, pos=pos, los=los, pol="RC") # %% Core calculations @@ -64,17 +62,17 @@ ws.measurement_vector[::100], np.array( [ - 227.78646795, - 230.8638575, - 234.80652899, - 240.36081974, - 249.78247057, - 207.62113428, - 249.78190355, - 240.35972683, - 234.80495168, - 230.86180615, - 227.78395156, + 227.85626271, + 230.93430882, + 234.89998118, + 240.50100578, + 250.05272664, + 209.9140708, + 249.51258095, + 240.21976958, + 234.71155853, + 230.79135297, + 227.73970752, ] ), ) diff --git a/src/core/jacobian/jacobian.cc b/src/core/jacobian/jacobian.cc index 2a53bffb43..0fdf228be1 100644 --- a/src/core/jacobian/jacobian.cc +++ b/src/core/jacobian/jacobian.cc @@ -150,7 +150,15 @@ void polyfit(ExhaustiveVectorView param, const Vector& x, const Vector& y) { using namespace Minimize; auto [success, fit] = curve_fit(x, y, param.size() - 1); - ARTS_USER_ERROR_IF(not success, "Could not fit polynomial.") + ARTS_USER_ERROR_IF(not success, + R"(Could not fit polynomial: + x: {:B,}, + y: {:B,} + fit: {:B,} +)", + x, + y, + Vector{fit}) ARTS_USER_ERROR_IF(fit.size() != param.size(), "Bad size. fit.size(): {}, param.size(): {}", fit.size(), @@ -234,6 +242,18 @@ bool AtmTarget::is_wind() const { type == AtmKey::wind_w; } +void SensorTarget::update(ArrayOfSensorObsel& sens, const Vector& x) const { + const auto sz = static_cast(x.size()); + ARTS_USER_ERROR_IF(sz < (x_start + x_size), "Got too small vector.") + set_model(sens, type, x.slice(x_start, x_size)); +} + +void SensorTarget::update(Vector& x, const ArrayOfSensorObsel& sens) const { + const auto sz = static_cast(x.size()); + ARTS_USER_ERROR_IF(sz < (x_start + x_size), "Got too small vector.") + set_state(x.slice(x_start, x_size), sens, type); +} + void SurfaceTarget::update(SurfaceField& surf, const Vector& x) const { const auto sz = static_cast(x.size()); ARTS_USER_ERROR_IF(sz < (x_start + x_size), "Got too small vector.") @@ -359,4 +379,36 @@ void Targets::finalize(const AtmField& atmospheric_field, throwing_check(last_size); } + +AtmTarget& Targets::emplace_back(AtmKeyVal&& t, Numeric d) { + return atm().emplace_back(std::move(t), d, target_count()); +} + +SurfaceTarget& Targets::emplace_back(SurfaceKeyVal&& t, Numeric d) { + return surf().emplace_back(std::move(t), d, target_count()); +} + +LineTarget& Targets::emplace_back(LblLineKey&& t, Numeric d) { + return line().emplace_back(std::move(t), d, target_count()); +} + +SensorTarget& Targets::emplace_back(SensorKey&& t, Numeric d) { + return sensor().emplace_back(std::move(t), d, target_count()); +} + +AtmTarget& Targets::emplace_back(const AtmKeyVal& t, Numeric d) { + return atm().emplace_back(t, d, target_count()); +} + +SurfaceTarget& Targets::emplace_back(const SurfaceKeyVal& t, Numeric d) { + return surf().emplace_back(t, d, target_count()); +} + +LineTarget& Targets::emplace_back(const LblLineKey& t, Numeric d) { + return line().emplace_back(t, d, target_count()); +} + +SensorTarget& Targets::emplace_back(const SensorKey& t, Numeric d) { + return sensor().emplace_back(t, d, target_count()); +} } // namespace Jacobian diff --git a/src/core/jacobian/jacobian.h b/src/core/jacobian/jacobian.h index 58c09eb879..6c54774e4f 100644 --- a/src/core/jacobian/jacobian.h +++ b/src/core/jacobian/jacobian.h @@ -275,7 +275,8 @@ struct targets_t { } }; -struct Targets final : targets_t { +struct Targets final + : targets_t { [[nodiscard]] const std::vector& atm() const; [[nodiscard]] const std::vector& surf() const; [[nodiscard]] const std::vector& line() const; @@ -291,13 +292,26 @@ struct Targets final : targets_t; + using variant_t = + std::variant; variant_t target; - template + template [[nodiscard]] constexpr auto apply(const AtmKeyValFunc& ifatm, const SurfaceKeyValFunc& ifsurf, const LblLineKeyFunc& ifline, @@ -534,7 +548,10 @@ struct std::formatter { v.surf(), sep, R"("line": )"sv, - v.line()); + v.line(), + sep, + R"("sensor": )"sv, + v.sensor()); tags.add_if_bracket(ctx, '}'); return ctx.out(); diff --git a/src/core/matpack/matpack_view.h b/src/core/matpack/matpack_view.h index 5a1f5afae8..0bd36713a5 100644 --- a/src/core/matpack/matpack_view.h +++ b/src/core/matpack/matpack_view.h @@ -1556,7 +1556,7 @@ struct std::formatter> { if (tags.short_str and n > 8) { for (auto&& a : v | take(3) | drop(1)) tags.format(ctx, sep, a); - tags.format(ctx, sep, "..."); + tags.format(ctx, sep, "..."sv); for (auto&& a : v | drop(n - 3)) tags.format(ctx, sep, a); } else { for (auto&& a : v | drop(1)) tags.format(ctx, sep, a); diff --git a/src/core/minimize.cc b/src/core/minimize.cc index 4026227399..61aa30f89a 100644 --- a/src/core/minimize.cc +++ b/src/core/minimize.cc @@ -1,13 +1,14 @@ #include "minimize.h" + #include "matpack_math.h" namespace Minimize { int Polynom::operator()(const T4::InputType& p, T4::ValueType& f) const { - for (Index i=0; i #include -#include -#include - -#include "configtypes.h" -#include "debug.h" -#include "isotopologues.h" -#include "lbl_data.h" -#include "lbl_lineshape_model.h" -#include "quantum_numbers.h" - void jacobian_targetsInit(JacobianTargets& jacobian_targets) { jacobian_targets.clear(); } @@ -27,55 +17,53 @@ void jacobian_targetsFinalize(JacobianTargets& jacobian_targets, void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, const SurfaceKey& key, const Numeric& d) { - jacobian_targets.surf().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(SurfaceKeyVal{key}, d); } void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, const SurfaceTypeTag& key, const Numeric& d) { - jacobian_targets.surf().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(SurfaceKeyVal{key}, d); } void jacobian_targetsAddSurface(JacobianTargets& jacobian_targets, const SurfacePropertyTag& key, const Numeric& d) { - jacobian_targets.surf().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(SurfaceKeyVal{key}, d); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const AtmKey& key, const Numeric& d) { - jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{key}, d); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const SpeciesEnum& key, const Numeric& d) { - jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{key}, d); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const SpeciesIsotope& key, const Numeric& d) { - jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{key}, d); } void jacobian_targetsAddAtmosphere(JacobianTargets& jacobian_targets, const QuantumIdentifier& key, const Numeric& d) { - jacobian_targets.atm().emplace_back(key, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{key}, d); } void jacobian_targetsAddTemperature(JacobianTargets& jacobian_targets, const Numeric& d) { - jacobian_targets.atm().emplace_back( - AtmKey::t, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{AtmKey::t}, d); } void jacobian_targetsAddPressure(JacobianTargets& jacobian_targets, const Numeric& d) { - jacobian_targets.atm().emplace_back( - AtmKey::p, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{AtmKey::p}, d); } void jacobian_targetsAddMagneticField(JacobianTargets& jacobian_targets, @@ -83,18 +71,9 @@ void jacobian_targetsAddMagneticField(JacobianTargets& jacobian_targets, const Numeric& d) { using enum FieldComponent; switch (to(component)) { - case u: - jacobian_targets.atm().emplace_back( - AtmKey::mag_u, d, jacobian_targets.target_count()); - break; - case v: - jacobian_targets.atm().emplace_back( - AtmKey::mag_v, d, jacobian_targets.target_count()); - break; - case w: - jacobian_targets.atm().emplace_back( - AtmKey::mag_w, d, jacobian_targets.target_count()); - break; + case u: jacobian_targets.emplace_back(AtmKeyVal{AtmKey::mag_u}, d); break; + case v: jacobian_targets.emplace_back(AtmKeyVal{AtmKey::mag_v}, d); break; + case w: jacobian_targets.emplace_back(AtmKeyVal{AtmKey::mag_w}, d); break; } } @@ -103,57 +82,23 @@ void jacobian_targetsAddWindField(JacobianTargets& jacobian_targets, const Numeric& d) { using enum FieldComponent; switch (to(component)) { - case u: - jacobian_targets.atm().emplace_back( - AtmKey::wind_u, d, jacobian_targets.target_count()); - break; - case v: - jacobian_targets.atm().emplace_back( - AtmKey::wind_v, d, jacobian_targets.target_count()); - break; - case w: - jacobian_targets.atm().emplace_back( - AtmKey::wind_w, d, jacobian_targets.target_count()); - break; + case u: jacobian_targets.emplace_back(AtmKeyVal{AtmKey::wind_u}, d); break; + case v: jacobian_targets.emplace_back(AtmKeyVal{AtmKey::wind_v}, d); break; + case w: jacobian_targets.emplace_back(AtmKeyVal{AtmKey::wind_w}, d); break; } } void jacobian_targetsAddSpeciesVMR(JacobianTargets& jacobian_targets, const SpeciesEnum& species, const Numeric& d) { - jacobian_targets.atm().emplace_back( - species, d, jacobian_targets.target_count()); + jacobian_targets.emplace_back(AtmKeyVal{species}, d); } void jacobian_targetsAddSpeciesIsotopologueRatio( JacobianTargets& jacobian_targets, const SpeciesIsotope& species, const Numeric& d) { - ARTS_USER_ERROR_IF( - std::ranges::none_of(Species::Isotopologues, Cmp::eq(species)), - "Unknown isotopologue: \"{}\"", - species.FullName()); - - jacobian_targets.atm().emplace_back( - species, d, jacobian_targets.target_count()); -} - -void jacobian_targetsAddSensor(JacobianTargets& jacobian_targets, - const ArrayOfSensorObsel& measurement_sensor, - const SensorKeyType& key, - const Numeric& d, - const Index& elem) { - ARTS_USER_ERROR_IF(measurement_sensor.size() <= static_cast(elem), - "Sensor element out of bounds: {}", - elem) - - jacobian_targets.sensor().push_back({ - .type = {.type = key, - .elem = elem, - .model = SensorJacobianModelType::None}, - .d = d, - .target_pos = jacobian_targets.target_count(), - }); + jacobian_targets.emplace_back(AtmKeyVal{species}, d); } void jacobian_targetsAddSensorFrequencyPolyFit( @@ -168,13 +113,12 @@ void jacobian_targetsAddSensorFrequencyPolyFit( ARTS_USER_ERROR_IF( polyorder < 0, "Polyorder must be non-negative: {}", polyorder) - Jacobian::SensorTarget target{ - .type = {.type = SensorKeyType::f, - .elem = elem, - .model = SensorJacobianModelType::PolynomialOffset, - .polyorder = polyorder, - .original_grid = measurement_sensor[elem].flat(SensorKeyType::f)}, - .d = d, - .target_pos = jacobian_targets.target_count(), - }; + jacobian_targets.emplace_back( + SensorKey{ + .type = SensorKeyType::f, + .elem = elem, + .model = SensorJacobianModelType::PolynomialOffset, + .polyorder = polyorder, + .original_grid = measurement_sensor[elem].flat(SensorKeyType::f)}, + d); } diff --git a/src/m_model_state.cc b/src/m_model_state.cc index 95b7b4c0fb..a9d316099a 100644 --- a/src/m_model_state.cc +++ b/src/m_model_state.cc @@ -56,3 +56,19 @@ void model_state_vectorFromBands(Vector& model_state_vector, target.update(model_state_vector, absorption_bands); } } + +void measurement_sensorFromModelState(ArrayOfSensorObsel& measurement_sensor, + const Vector& model_state_vector, + const JacobianTargets& jacobian_targets) { + for (auto& target : jacobian_targets.sensor()) { + target.update(measurement_sensor, model_state_vector); + } +} + +void model_state_vectorFromSensor(Vector& model_state_vector, + const ArrayOfSensorObsel& measurement_sensor, + const JacobianTargets& jacobian_targets) { + for (auto& target : jacobian_targets.sensor()) { + target.update(model_state_vector, measurement_sensor); + } +} diff --git a/src/m_rad.cc b/src/m_rad.cc index e848f7055e..658901fc11 100644 --- a/src/m_rad.cc +++ b/src/m_rad.cc @@ -265,24 +265,26 @@ frequency_grid.size() = {} const auto b = x.begin(); const auto e = x.end(); + bool find_any = false; for (auto &target : jacobian_targets.sensor()) { ARTS_USER_ERROR_IF( measurement_sensor.size() <= static_cast(target.type.elem), "Sensor element out of bounds"); + auto &elem = measurement_sensor[target.type.elem]; auto m = spectral_radiance_jacobian.slice(target.x_start, target.x_size); const Numeric d = target.d; // Check that the Jacobian targets are represented by this frequency grid and this pos-los pair - const Index iposlos = measurement_sensor[target.type.elem].find(pos, los); + const Index iposlos = elem.find(pos, los); if (iposlos == SensorObsel::dont_have) continue; - if (measurement_sensor[target.type.elem].find(frequency_grid) == - SensorObsel::dont_have) - continue; + if (elem.find(frequency_grid) == SensorObsel::dont_have) continue; + + find_any = true; using enum SensorKeyType; switch (target.type.type) { - case f: call({b, e, [d](auto x) { return x + d; }}, pos, los, d); + case f: call({b, e, [d](auto x) { return x + d; }}, pos, los, d); break; case za: call(x, pos, {los[0] + d, los[1]}, d); break; case aa: call(x, pos, {los[0], los[1] + d}, d); break; case alt: call(x, {pos[0] + d, pos[1], pos[2]}, los, d); break; @@ -331,6 +333,22 @@ frequency_grid.size() = {} } break; } } + + ARTS_USER_ERROR_IF(not find_any, + R"(No sensor element found for pos-los/frequency grid pair + + frequency_grid: {:Bs,} + pos: {:B,} + los: {:B,} + +Note: It is not allowed to change the frequency grid or the pos-los pair in an agenda +that calls this function. This is because the actual memory address is used to identify +a sensor element. Modifying pos, los or frequency_grid will copy the data to a new memory +location and the sensor element will not be found. +)", + frequency_grid, + pos, + los); } ARTS_METHOD_ERROR_CATCH @@ -389,6 +407,7 @@ f_grid_ptr->size() = {} )", spectral_radiance.size(), f_grid_ptr->size()) + ARTS_USER_ERROR_IF( (spectral_radiance_jacobian.shape() != std::array{measurement_jacobian.ncols(), f_grid_ptr->size()}), diff --git a/src/workspace_groups.cpp b/src/workspace_groups.cpp index 09ff1aaf58..52a5ef27db 100644 --- a/src/workspace_groups.cpp +++ b/src/workspace_groups.cpp @@ -2,6 +2,8 @@ #include +#include + std::unordered_map internal_workspace_groups_creator() { std::unordered_map wsg_data; @@ -640,7 +642,6 @@ if it is contained in the Atlas and NAN otherwise. .desc = "A map of vibrational energy levels for NLTE calculations\n", }; - // rtepack types wsg_data["Propmat"] = { .file = "rtepack.h", .desc = R"--(A single propagation matrix. @@ -989,7 +990,7 @@ well as the sampling device's polarization response. for (auto& g : internal_options()) { if (wsg_data.find(g.name) != wsg_data.end()) - throw std::runtime_error("Duplicate workspace group name: " + g.name); + throw std::runtime_error("Duplicate workspace group name (name is reserved as options-group): " + g.name); wsg_data[g.name] = { .file = "enums.h", diff --git a/src/workspace_meta_methods.cpp b/src/workspace_meta_methods.cpp index 5c672a789c..0c969ac5ce 100644 --- a/src/workspace_meta_methods.cpp +++ b/src/workspace_meta_methods.cpp @@ -17,8 +17,7 @@ std::vector internal_meta_methods_creator() { .name = "measurement_sensorSimple", .desc = "Wrapper for a single simple dirac-opening sensor", .author = {"Richard Larsson"}, - .methods = {"measurement_sensorInit", - "measurement_sensorAddSimple"}, + .methods = {"measurement_sensorInit", "measurement_sensorAddSimple"}, .out = {"measurement_sensor"}, }); @@ -195,8 +194,12 @@ std::vector internal_meta_methods_creator() { .author = {"Richard Larsson"}, .methods = {"absorption_bandsFromModelState", "surface_fieldFromModelState", - "atmospheric_fieldFromModelState"}, - .out = {"absorption_bands", "surface_field", "atmospheric_field"}, + "atmospheric_fieldFromModelState", + "measurement_sensorFromModelState"}, + .out = {"absorption_bands", + "surface_field", + "atmospheric_field", + "measurement_sensor"}, }); out.push_back(WorkspaceMethodInternalMetaRecord{ @@ -207,7 +210,8 @@ std::vector internal_meta_methods_creator() { "model_state_vectorZero", "model_state_vectorFromAtmosphere", "model_state_vectorFromSurface", - "model_state_vectorFromBands"}, + "model_state_vectorFromBands", + "model_state_vectorFromSensor"}, .out = {"model_state_vector"}, }); @@ -266,15 +270,15 @@ WorkspaceMethodInternalRecord WorkspaceMethodInternalMetaRecord::create( const auto ptr = wsms.find(m); if (ptr == wsms.end()) { - throw std::runtime_error("Method " + m + " not found"); + throw std::runtime_error(std::format(R"(Method "{}" not found)", m)); } const auto& wm = ptr->second; if (wm.has_any() or wm.has_overloads()) { - throw std::runtime_error( - "Method " + m + - " has overloads and does not work with meta-functions"); + throw std::runtime_error(std::format( + R"(Method "{}" has overloads and does not work with meta-functions)", + m)); } wsm.author.insert(wsm.author.end(), wm.author.begin(), wm.author.end()); diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index 2e470f94bb..c8a533c041 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -1994,14 +1994,15 @@ building of an actual Jacobian matrix. }; const auto jac2ret = [&wsm_data](const std::string& name) { - auto v = wsm_data[name]; + auto v = wsm_data[name]; v.desc += std::format(R"( This method wraps *{}* together with adding the covariance matrices, to the *covariance_matrix_diagonal_blocks*, which are required to perform *OEM*. The input covariance matrices must fit the size of the later computed model state represented by the *jacobian_targets*. The covariance matrix inverse -)", name); +)", + name); v.out.insert(v.out.begin() + 1, "covariance_matrix_diagonal_blocks"); v.in.insert(v.in.begin() + 1, "covariance_matrix_diagonal_blocks"); v.gin.insert(v.gin.end(), "matrix"); @@ -2087,7 +2088,8 @@ See *FieldComponent* for valid ``component`` {"The component to use [u, v, w]", "The perturbation used in methods that cannot compute derivatives analytically"}, }; - wsm_data["RetrievalAddMagneticField"] = jac2ret("jacobian_targetsAddMagneticField"); + wsm_data["RetrievalAddMagneticField"] = + jac2ret("jacobian_targetsAddMagneticField"); wsm_data["jacobian_targetsAddWindField"] = { .desc = R"--(Set wind field derivative @@ -2107,7 +2109,7 @@ See *FieldComponent* for valid ``component`` {"The component to use [u, v, w]", "The perturbation used in methods that cannot compute derivatives analytically"}, }; -wsm_data["RetrievalAddWindField"] = jac2ret("jacobian_targetsAddWindField"); + wsm_data["RetrievalAddWindField"] = jac2ret("jacobian_targetsAddWindField"); wsm_data["jacobian_targetsAddSpeciesVMR"] = { .desc = R"--(Set volume mixing ratio derivative @@ -2124,7 +2126,7 @@ See *SpeciesEnum* for valid ``species`` {"The species of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; -wsm_data["RetrievalAddSpeciesVMR"] = jac2ret("jacobian_targetsAddSpeciesVMR"); + wsm_data["RetrievalAddSpeciesVMR"] = jac2ret("jacobian_targetsAddSpeciesVMR"); wsm_data["jacobian_targetsAddAtmosphere"] = { .desc = R"--(Sets an atmospheric target @@ -2140,7 +2142,7 @@ wsm_data["RetrievalAddSpeciesVMR"] = jac2ret("jacobian_targetsAddSpeciesVMR"); {"The target of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; -wsm_data["RetrievalAddAtmosphere"] = jac2ret("jacobian_targetsAddAtmosphere"); + wsm_data["RetrievalAddAtmosphere"] = jac2ret("jacobian_targetsAddAtmosphere"); wsm_data["jacobian_targetsAddSurface"] = { .desc = R"--(Sets a surface target @@ -2155,7 +2157,7 @@ wsm_data["RetrievalAddAtmosphere"] = jac2ret("jacobian_targetsAddAtmosphere"); {"The target of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; -wsm_data["RetrievalAddSurface"] = jac2ret("jacobian_targetsAddSurface"); + wsm_data["RetrievalAddSurface"] = jac2ret("jacobian_targetsAddSurface"); wsm_data["jacobian_targetsAddSpeciesIsotopologueRatio"] = { .desc = R"--(Set isotopologue ratio derivative @@ -2172,7 +2174,8 @@ See *SpeciesIsotope* for valid ``species`` {"The species isotopologue of interest", "The perturbation used in methods that cannot compute derivatives analytically"}, }; -wsm_data["RetrievalAddSpeciesIsotopologueRatio"] = jac2ret("jacobian_targetsAddSpeciesIsotopologueRatio"); + wsm_data["RetrievalAddSpeciesIsotopologueRatio"] = + jac2ret("jacobian_targetsAddSpeciesIsotopologueRatio"); wsm_data["absorption_bandsReadHITRAN"] = { .desc = @@ -3175,6 +3178,15 @@ before consumption by, e.g., *OEM* .pass_workspace = true, }; + wsm_data["measurement_sensorFromModelState"] = { + .desc = + R"--(Update *measurement_sensor* from *model_state_vector*. +)--", + .author = {"Richard Larsson"}, + .out = {"measurement_sensor"}, + .in = {"measurement_sensor", "model_state_vector", "jacobian_targets"}, + }; + wsm_data["measurement_sensorInit"] = { .desc = R"--(Initialize *measurement_sensor* to empty. @@ -3464,6 +3476,14 @@ Hence, a temperature of 0 means 0s the edges of the *frequency_grid*. .in = {"model_state_vector"}, }; + wsm_data["model_state_vectorFromSensor"] = { + .desc = R"--(Sets *model_state_vector*'s sensor part. +)--", + .author = {"Richard Larsson"}, + .out = {"model_state_vector"}, + .in = {"model_state_vector", "measurement_sensor", "jacobian_targets"}, + }; + wsm_data["model_state_vectorFromAtmosphere"] = { .desc = R"--(Sets *model_state_vector*'s atmospheric part. )--", diff --git a/tests/core/sensor/jac.py b/tests/core/sensor/jac.py index 7abb537673..832a74e9a2 100644 --- a/tests/core/sensor/jac.py +++ b/tests/core/sensor/jac.py @@ -2,12 +2,20 @@ import numpy as np import matplotlib.pyplot as plt +LIMIT = 50 +noise = 0.5 +NF = 1001 +std = 1 +pol1 = "RC" +pol2 = "Ih" +ATOL = 1000 + ws = pyarts.workspace.Workspace() # %% Sampled frequency range line_f0 = 118750348044.712 -f = np.linspace(-5e6, 5e6, 1001) + line_f0 +f = np.linspace(-5e6, 5e6, NF) + line_f0 # %% Species and line absorption @@ -43,16 +51,10 @@ los = [180.0, 0.0] ws.measurement_sensorSimpleGaussian( - frequency_grid=f, std=1e5, pos=pos, los=los, pol="Ih" -) -ws.measurement_sensorAddSimpleGaussian( - frequency_grid=f, std=1e5, pos=pos, los=los, pol="Iv" -) -ws.measurement_sensorAddSimpleGaussian( - frequency_grid=f, std=1e5, pos=pos, los=los, pol="RC" + frequency_grid=f, std=std, pos=pos, los=los, pol=pol1 ) ws.measurement_sensorAddSimpleGaussian( - frequency_grid=f, std=1e5, pos=pos, los=los, pol="LC" + frequency_grid=f, std=std, pos=pos, los=los, pol=pol2 ) # %% Original calculations @@ -62,34 +64,79 @@ # %% Modify the sensor -DF1 = 1e6 -DF2 = 2e6 -DF3 = 3e6 -DF4 = -4e6 +DF1 = 1e5 ws.measurement_sensorSimpleGaussian( - frequency_grid=f + DF1, std=1e5, pos=pos, los=los, pol="Ih" -) -ws.measurement_sensorAddSimpleGaussian( - frequency_grid=f + DF2, std=1e5, pos=pos, los=los, pol="Iv" + frequency_grid=f + DF1, std=std, pos=pos, los=los, pol=pol1 ) ws.measurement_sensorAddSimpleGaussian( - frequency_grid=f + DF3, std=1e5, pos=pos, los=los, pol="RC" -) -ws.measurement_sensorAddSimpleGaussian( - frequency_grid=f + DF4, std=1e5, pos=pos, los=los, pol="LC" + frequency_grid=f + 2 * DF1, std=std, pos=pos, los=los, pol=pol2 ) ws.measurement_vectorFromSensor() mod = ws.measurement_vector * 1.0 +# %% Retrieval agenda + + +@pyarts.workspace.arts_agenda(ws=ws, fix=True) +def inversion_iterate_agenda(ws): + ws.UpdateModelStates() + ws.measurement_vectorFromSensor() + ws.measurement_vector_fittedFromMeasurement() + + # %% Set up the retrieval ws.RetrievalInit() -ws.RetrievalAddSensorFrequencyPolyFit(elem=0, d=1e3, matrix=np.diag(np.ones((1)) * 1e6)) +ws.RetrievalAddSensorFrequencyPolyFit( + elem=0, d=1e3, matrix=np.diag(np.ones((1)) * 1e10), polyorder=0 +) +ws.RetrievalAddSensorFrequencyPolyFit( + elem=1001, d=1e3, matrix=np.diag(np.ones((1)) * 1e10), polyorder=0 +) ws.RetrievalFinalizeDiagonal() -ws.measurement_vector_fitted = [] -ws.model_state_vector = [] -ws.measurement_jacobian = [[]] -# %% +# %% Perform OEM retrieval of frequency grid + +fail = True +for i in range(LIMIT): + ws.measurement_vector_fitted = [] + ws.model_state_vector = [] + ws.measurement_jacobian = [[]] + + ws.measurement_sensorSimpleGaussian( + frequency_grid=f, std=std, pos=pos, los=los, pol=pol1 + ) + ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f, std=std, pos=pos, los=los, pol=pol2 + ) + ws.measurement_vectorFromSensor() + ws.measurement_vector += np.random.normal(0, noise, 2 * NF) + + ws.measurement_sensorSimpleGaussian( + frequency_grid=f + DF1, std=std, pos=pos, los=los, pol=pol1 + ) + ws.measurement_sensorAddSimpleGaussian( + frequency_grid=f + 2 * DF1, std=std, pos=pos, los=los, pol=pol2 + ) + + ws.measurement_vector_fitted = [] + ws.model_state_vector = [] + ws.measurement_jacobian = [[]] + + ws.model_state_vector_aprioriFromData() + + ws.measurement_vector_error_covariance_matrixConstant(value=noise**2) + + ws.OEM(method="gn") + + print(f"got: [{ws.model_state_vector[0]}, {ws.model_state_vector[1]}]") + print(f"expected: [{-DF1}, {-2*DF1}]") + if np.allclose(ws.model_state_vector, [-DF1, -2 * DF1], atol=ATOL): + print(f"Within {ATOL} Hz. Success!") + fail = False + break + print(f"AbsDiff not less than {ATOL} HZ, rerunning with new random noise") + +assert not fail, "Failed to retrieve frequency grid" From 00f6251b3ff1e5db86343e0a2627c75eeea9c46b Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Thu, 21 Nov 2024 15:07:15 +0100 Subject: [PATCH 6/7] Rename --- ...ac.py => spectral_radiance_jacobian_multiple_sensor_f_grid.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/core/sensor/{jac.py => spectral_radiance_jacobian_multiple_sensor_f_grid.py} (100%) diff --git a/tests/core/sensor/jac.py b/tests/core/sensor/spectral_radiance_jacobian_multiple_sensor_f_grid.py similarity index 100% rename from tests/core/sensor/jac.py rename to tests/core/sensor/spectral_radiance_jacobian_multiple_sensor_f_grid.py From bb0a7481ce1b4481aba4c3cf4b6f75a9186633c1 Mon Sep 17 00:00:00 2001 From: Richard Larsson Date: Thu, 21 Nov 2024 15:46:12 +0100 Subject: [PATCH 7/7] Use correct include --- src/workspace_groups.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/workspace_groups.cpp b/src/workspace_groups.cpp index 52a5ef27db..05a799efad 100644 --- a/src/workspace_groups.cpp +++ b/src/workspace_groups.cpp @@ -2,7 +2,7 @@ #include -#include +#include std::unordered_map internal_workspace_groups_creator() {