Skip to content

Commit

Permalink
MIR- eckit::geo ProxyWeightedMethod
Browse files Browse the repository at this point in the history
  • Loading branch information
pmaciel committed Nov 19, 2024
1 parent 66ecf64 commit 56eb4a0
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 46 deletions.
91 changes: 52 additions & 39 deletions src/mir/method/ProxyWeightedMethod.cc
Original file line number Diff line number Diff line change
Expand Up @@ -13,33 +13,44 @@
#include "mir/method/ProxyWeightedMethod.h"

#include <ostream>
#include <vector>

#include "eckit/log/JSON.h"
#include "eckit/utils/MD5.h"

#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"


namespace mir::method {


namespace {


struct StructuredBicubic final : public ProxyWeightedMethod {
explicit StructuredBicubic(const param::MIRParametrisation& param) :
ProxyWeightedMethod(param, "structured-bicubic") {}
ProxyWeightedMethod(param, "structured-bicubic", "serial-halo to serial") {}
};


} // namespace


static const MethodFactory* METHODS[]{
new MethodBuilder<StructuredBicubic>("structured-bicubic"),
};


ProxyWeightedMethod::ProxyWeightedMethod(const param::MIRParametrisation& param, const std::string& type) :
MethodWeighted(param), type_(type) {
options_ = {"type", type};
options_.set("matrix_free", false);
ProxyWeightedMethod::ProxyWeightedMethod(const param::MIRParametrisation& param, const std::string& interpolation_type,
const std::string& renumber_type) :
MethodWeighted(param), type_(interpolation_type) {
interpol_.set("matrix_free", false);
interpol_.set("type", interpolation_type);
renumber_.set("type", renumber_type);
}


Expand All @@ -55,70 +66,72 @@ int ProxyWeightedMethod::version() const {

void ProxyWeightedMethod::hash(eckit::MD5& h) const {
MethodWeighted::hash(h);
h.add(options_);
h.add(interpol_);
}


bool ProxyWeightedMethod::sameAs(const Method& other) const {
auto digest = [](const auto& options) {
auto digest = [](const auto& config) {
eckit::MD5 h;
h.add(options);
config.hash(h);
return h.digest();
};

const auto* o = dynamic_cast<const ProxyWeightedMethod*>(&other);
return (o != nullptr) && name() == std::string{o->name()} && digest(options_) == digest(o->options_) &&
MethodWeighted::sameAs(*o);
return (o != nullptr) && name() == std::string{o->name()} && digest(interpol_) == digest(o->interpol_) &&
digest(renumber_) == digest(o->renumber_) && MethodWeighted::sameAs(*o);
}


void ProxyWeightedMethod::print(std::ostream& out) const {
out << "ProxyWeightedMethod[options=" << options_ << ",";
out << "ProxyWeightedMethod[interpolation=" << interpol_ << ",renumber=" << renumber_ << ",";
MethodWeighted::print(out);
out << "]";
}


void ProxyWeightedMethod::json(eckit::JSON& j) const {
j.startObject();
j << "options" << options_.json(eckit::JSON::Formatting::compact());
j << "options" << interpol_.json(eckit::JSON::Formatting::compact());
MethodWeighted::json(j);
j.endObject();
}


void ProxyWeightedMethod::assemble(util::MIRStatistics&, WeightMatrix& W, const repres::Representation& in,
const repres::Representation& out) const {
atlas::Interpolation interpol(options_, in.atlasGrid(), out.atlasGrid());

auto cache = interpol.createCache();
ASSERT(cache);

const auto* entry = dynamic_cast<const atlas::interpolation::MatrixCacheEntry*>(cache.get("Matrix"));
ASSERT(entry != nullptr);

W.swap(const_cast<eckit::linalg::SparseMatrix&>(entry->matrix()));

auto& fs = interpol.source();
const auto gi{atlas::array::make_view<atlas::gidx_t, 1>(fs.global_index())};

atlas::gidx_t min = std::numeric_limits<atlas::gidx_t>::max();
atlas::gidx_t max = std::numeric_limits<atlas::gidx_t>::min();
for (size_t i = 0; i < gi.size(); ++i) {
min = std::min(min, gi(i));
max = std::max(max, gi(i));
// interpolation matrix build and assign (with a halo on the source grid)
atlas::Interpolation interpol{interpol_, in.atlasGrid(), out.atlasGrid()};
{
atlas::interpolation::MatrixCache cache{interpol};
const auto& M = cache.matrix();
W.swap(const_cast<eckit::linalg::SparseMatrix&>(M));
}

std::cout << "min: " << min << std::endl;
std::cout << "max: " << max << std::endl;

// 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<atlas::gidx_t, 1>(fs.global_index())};

auto* a = const_cast<eckit::linalg::Scalar*>(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();
}

auto Nr = W.rows();
auto Nc = W.cols();
auto Nin = in.numberOfPoints();
auto Nout = out.numberOfPoints();
ASSERT(Nr == Nout);
ASSERT(Nc == Nin);
// reshape matrix
WeightMatrix M(new util::allocator::InPlaceAllocator{
Nr, Nc, W.nonZeros(), const_cast<eckit::linalg::Index*>(W.outer()),
const_cast<eckit::linalg::Index*>(W.inner()), const_cast<eckit::linalg::Scalar*>(W.data())});
W.swap(M);
}


Expand Down
11 changes: 4 additions & 7 deletions src/mir/method/ProxyWeightedMethod.h
Original file line number Diff line number Diff line change
Expand Up @@ -24,22 +24,19 @@ class ProxyWeightedMethod : public MethodWeighted {
protected:
// -- Constructor

ProxyWeightedMethod(const param::MIRParametrisation&, const std::string& type);
ProxyWeightedMethod(const param::MIRParametrisation&, const std::string& interpolation_type,
const std::string& renumber_type = "");

// -- Destructor

~ProxyWeightedMethod() override = default;

// -- Methods

const auto& options() const { return options_; }
auto& options() { return options_; }

private:
// -- Members

const std::string type_;
atlas::util::Config options_;
atlas::util::Config interpol_;
atlas::util::Config renumber_;

// -- Overridden methods

Expand Down

0 comments on commit 56eb4a0

Please sign in to comment.