From 974fe3377d9b8e6b743794d7742a671ccdd67dd7 Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Sun, 17 Nov 2024 09:07:46 +0000 Subject: [PATCH] MIR- eckit::geo ProxyWeightedMethod --- src/mir/method/ProxyWeightedMethod.cc | 116 +++++++++++++++----------- src/mir/method/ProxyWeightedMethod.h | 11 +-- 2 files changed, 71 insertions(+), 56 deletions(-) diff --git a/src/mir/method/ProxyWeightedMethod.cc b/src/mir/method/ProxyWeightedMethod.cc index 680ee5e6d..813f67c4c 100644 --- a/src/mir/method/ProxyWeightedMethod.cc +++ b/src/mir/method/ProxyWeightedMethod.cc @@ -13,7 +13,6 @@ #include "mir/method/ProxyWeightedMethod.h" #include -#include #include "eckit/log/JSON.h" #include "eckit/utils/MD5.h" @@ -21,7 +20,6 @@ #include "mir/param/MIRParametrisation.h" #include "mir/repres/Representation.h" #include "mir/util/Exceptions.h" -#include "mir/util/Log.h" #include "mir/util/allocator/InPlaceAllocator.h" @@ -31,26 +29,72 @@ namespace mir::method { namespace { -struct StructuredBicubic final : public ProxyWeightedMethod { - explicit StructuredBicubic(const param::MIRParametrisation& param) : - ProxyWeightedMethod(param, "structured-bicubic", "serial-halo to serial") {} -}; +#define ATLAS_METHOD(Type, Name) \ + struct Type final : ProxyWeightedMethod { \ + explicit Type(const param::MIRParametrisation& param) : ProxyWeightedMethod(param, Name) {} \ + }; +ATLAS_METHOD(StructuredBilinear, "structured-bilinear"); +ATLAS_METHOD(StructuredBiquasicubic, "structured-biquasicubic"); +ATLAS_METHOD(StructuredBicubic, "structured-bicubic"); +ATLAS_METHOD(FiniteElement, "finite-element"); +ATLAS_METHOD(ConservativeSphericalPolygon, "conservative-spherical-polygon"); +ATLAS_METHOD(GridBoxAverage, "grid-box-average"); +ATLAS_METHOD(GridBoxMaximum, "grid-box-maximum"); -} // namespace +// "nearest-neighbour" (knn) +// "k-nearest-neighbours" (knn) +// "cubedsphere-bilinear" +// "regional-linear-2d" (structured) +// "spherical-vector" + + +#undef ATLAS_INTERPOL -static const MethodFactory* METHODS[]{ - new MethodBuilder("structured-bicubic"), +const MethodFactory* MIR_METHODS[]{ + new MethodBuilder("atlas-structured-bicubic"), + new MethodBuilder("atlas-structured-bilinear"), + new MethodBuilder("atlas-structured-biquasicubic"), + new MethodBuilder("atlas-finite-element"), + new MethodBuilder("atlas-conservative-spherical-polygon"), + new MethodBuilder("atlas-grid-box-average"), + new MethodBuilder("atlas-grid-box-maximum"), }; -ProxyWeightedMethod::ProxyWeightedMethod(const param::MIRParametrisation& param, const std::string& interpolation_type, - const std::string& renumber_type) : +} // namespace + + +ProxyWeightedMethod::ProxyWeightedMethod(const param::MIRParametrisation& param, + const std::string& interpolation_type) : MethodWeighted(param), type_(interpolation_type) { interpol_.set("matrix_free", false); interpol_.set("type", interpolation_type); - renumber_.set("type", renumber_type); +} + + +void ProxyWeightedMethod::foldSourceHalo(const atlas::Interpolation& interpol, size_t Nr, size_t Nc, + WeightMatrix& W) const { + ASSERT(Nr == W.rows()); + ASSERT(Nc < W.cols()); + + const auto& fs = interpol.source(); + const auto global_index{atlas::array::make_view(fs.global_index())}; + ASSERT(global_index.size() == W.cols()); + + auto* a = const_cast(W.data()); + for (auto c = Nc; c < W.cols(); ++c) { + ASSERT(1 <= global_index[c] && global_index[c] < Nc + 1); // (global indexing is 1-based) + a[global_index[c] - 1] += a[c]; + a[c] = 0.; + } + + W.prune(); + WeightMatrix M(new util::allocator::InPlaceAllocator{ + Nr, Nc, W.nonZeros(), const_cast(W.outer()), + const_cast(W.inner()), const_cast(W.data())}); + W.swap(M); } @@ -60,7 +104,7 @@ const char* ProxyWeightedMethod::name() const { int ProxyWeightedMethod::version() const { - return -1; + return 1; } @@ -79,12 +123,12 @@ bool ProxyWeightedMethod::sameAs(const Method& other) const { const auto* o = dynamic_cast(&other); return (o != nullptr) && name() == std::string{o->name()} && digest(interpol_) == digest(o->interpol_) && - digest(renumber_) == digest(o->renumber_) && MethodWeighted::sameAs(*o); + MethodWeighted::sameAs(*o); } void ProxyWeightedMethod::print(std::ostream& out) const { - out << "ProxyWeightedMethod[interpolation=" << interpol_ << ",renumber=" << renumber_ << ","; + out << "ProxyWeightedMethod[interpolation=" << interpol_ << ","; MethodWeighted::print(out); out << "]"; } @@ -92,7 +136,7 @@ void ProxyWeightedMethod::print(std::ostream& out) const { void ProxyWeightedMethod::json(eckit::JSON& j) const { j.startObject(); - j << "options" << interpol_.json(eckit::JSON::Formatting::compact()); + j << "interpolation" << interpol_.json(eckit::JSON::Formatting::compact()); MethodWeighted::json(j); j.endObject(); } @@ -100,43 +144,13 @@ void ProxyWeightedMethod::json(eckit::JSON& j) const { void ProxyWeightedMethod::assemble(util::MIRStatistics&, WeightMatrix& W, const repres::Representation& in, const repres::Representation& out) const { - // interpolation matrix build and assign (with a halo on the source grid) + // build matrix (with a halo on the source grid), move out of cache atlas::Interpolation interpol{interpol_, in.atlasGrid(), out.atlasGrid()}; - { - atlas::interpolation::MatrixCache cache{interpol}; - const auto& M = cache.matrix(); - W.swap(const_cast(M)); - } - - // fold serial+halo into serial - const auto Nr = out.numberOfPoints(); - const auto Nc = in.numberOfPoints(); - ASSERT(Nr == W.rows()); - ASSERT(Nc < W.cols()); - - { - const auto& fs = interpol.source(); - const auto gidx{atlas::array::make_view(fs.global_index())}; - - auto* a = const_cast(W.data()); - for (auto c = Nc; c < W.cols(); ++c) { - ASSERT(1 <= gidx[c] && gidx[c] < Nc + 1); // (global indexing is 1-based) - a[gidx[c] - 1] += a[c]; - a[c] = 0.; - } - W.prune(); - } - - // reshape matrix - WeightMatrix M(new util::allocator::InPlaceAllocator{ - Nr, Nc, W.nonZeros(), const_cast(W.outer()), - const_cast(W.inner()), const_cast(W.data())}); - W.swap(M); -} - + atlas::interpolation::MatrixCache cache{interpol}; + W.swap(const_cast(cache.matrix())); -bool ProxyWeightedMethod::validateMatrixWeights() const { - return true; + // fold source grid halo (from serial + halo into serial) + foldSourceHalo(interpol, out.numberOfPoints(), in.numberOfPoints(), W); } diff --git a/src/mir/method/ProxyWeightedMethod.h b/src/mir/method/ProxyWeightedMethod.h index 19b32cd72..410c5d272 100644 --- a/src/mir/method/ProxyWeightedMethod.h +++ b/src/mir/method/ProxyWeightedMethod.h @@ -24,19 +24,22 @@ class ProxyWeightedMethod : public MethodWeighted { protected: // -- Constructor - ProxyWeightedMethod(const param::MIRParametrisation&, const std::string& interpolation_type, - const std::string& renumber_type = ""); + ProxyWeightedMethod(const param::MIRParametrisation&, const std::string& interpolation_type); // -- Destructor ~ProxyWeightedMethod() override = default; +protected: + // -- Methods + + void foldSourceHalo(const atlas::Interpolation&, size_t Nr, size_t Nc, WeightMatrix&) const; + private: // -- Members const std::string type_; atlas::util::Config interpol_; - atlas::util::Config renumber_; // -- Overridden methods @@ -51,8 +54,6 @@ class ProxyWeightedMethod : public MethodWeighted { void assemble(util::MIRStatistics&, WeightMatrix&, const repres::Representation& in, const repres::Representation& out) const override; - - bool validateMatrixWeights() const override; };