Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Added multiplication methods to projection::Jacobian. #93

Merged
merged 3 commits into from
Feb 17, 2022
Merged

Added multiplication methods to projection::Jacobian. #93

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d244949
ATLAS-344 Fix PolygonXY::contains, requiring tolerances at edges
wdeconinck Feb 8, 2022
09cde23
Add macro to check eckit version
wdeconinck Feb 17, 2022
e4b12c7
Added methods to projection::Jacobian (#93)
odlomax Feb 17, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 14 additions & 6 deletions src/atlas/grid/detail/partitioner/MatchingMeshPartitionerLonLatPolygon.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,9 +76,13 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti

int min = compute(east - 360.);
constexpr double eps = 1.e-10;
if (min < 0 && east - west > 360. + eps) {
min = compute(west - eps);
}
bool second_try = [&]() {
if (min < 0 && east - west > 360. + eps) {
min = compute(west - eps);
return true;
}
return false;
}();
if (min < 0) {
size_t i = 0;
size_t max_failures = grid.size();
Expand All @@ -95,10 +99,14 @@ void MatchingMeshPartitionerLonLatPolygon::partition(const Grid& grid, int parti
}
size_t nb_failures = failed_index.size();
std::stringstream err;
err.precision(20);
err << "Could not find partition of " << nb_failures
<< " target grid points (source "
"mesh does not contain all target grid points)\n"
"Failed target grid points with global index:\n";
<< " target grid points (source mesh does not contain all target grid points)\n"
<< "Tried first normalizing coordinates with west=" << east - 360.;
if (second_try) {
err << "Tried second time normalizing coordinates with west=" << west - eps << "\n";
}
err << "Failed target grid points with global index:\n";
for (size_t n = 0; n < nb_failures; ++n) {
err << " - " << std::setw(10) << std::left << failed_index[n] + 1 << " {lon,lat} : " << failed_lonlat[n]
<< "\n";
Expand Down
4 changes: 3 additions & 1 deletion src/atlas/library/config.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,9 @@ using idx_t = long;
/// Integer type for unique indices
typedef gidx_t uidx_t;

#if ATLAS_ECKIT_VERSION_INT >= 11900 // eckit >= 1.19
#define ATLAS_ECKIT_VERSION_AT_LEAST(x, y, z) (ATLAS_ECKIT_VERSION_INT >= x * 10000 + y * 100 + z)

#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 19, 0)
#define ATLAS_ECKIT_HAVE_ECKIT_585 1
#else
#define ATLAS_ECKIT_HAVE_ECKIT_585 0
Expand Down
127 changes: 51 additions & 76 deletions src/atlas/projection/detail/ProjectionImpl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,12 @@

#pragma once

#include <initializer_list>
#include <iostream>
#include <memory>
#include <string>

#include "atlas/runtime/Exception.h"
#include "atlas/util/Factory.h"
#include "atlas/util/NormaliseLongitude.h"
#include "atlas/util/Object.h"
Expand Down Expand Up @@ -41,97 +44,68 @@ class ProjectionImpl : public util::Object {
using Spec = atlas::util::Config;
class Jacobian : public std::array<std::array<double, 2>, 2> {
public:
static Jacobian identity() {
Jacobian id;
id[0] = {1.0, 0.0};
id[1] = {0.0, 1.0};
return id;
}
using std::array<std::array<double, 2>, 2>::array;

Jacobian inverse() const {
const Jacobian& jac = *this;
Jacobian inv;
double det = jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0];
inv[0][0] = +jac[1][1] / det;
inv[0][1] = -jac[0][1] / det;
inv[1][0] = -jac[1][0] / det;
inv[1][1] = +jac[0][0] / det;
return inv;
};
static Jacobian identity() { return Jacobian{1., 0., 0., 1.}; }

Jacobian transpose() const {
Jacobian tra = *this;
std::swap(tra[0][1], tra[1][0]);
return tra;
};
Jacobian(double j00, double j01, double j10, double j11):
array{std::array<double, 2>{j00, j01}, std::array<double, 2>{j10, j11}} {}

Jacobian operator-(const Jacobian& jac2) const {
const Jacobian& jac1 = *this;
Jacobian jac;
jac[0][0] = jac1[0][0] - jac2[0][0];
jac[0][1] = jac1[0][1] - jac2[0][1];
jac[1][0] = jac1[1][0] - jac2[1][0];
jac[1][1] = jac1[1][1] - jac2[1][1];
return jac;
}
Jacobian(std::initializer_list<std::initializer_list<double>> list):
Jacobian{*(list.begin()->begin()), *(list.begin()->begin() + 1), *((list.begin() + 1)->begin()),
*((list.begin() + 1)->begin() + 1)} {}

Jacobian operator+(const Jacobian& jac2) const {
const Jacobian& jac1 = *this;
Jacobian jac;
jac[0][0] = jac1[0][0] + jac2[0][0];
jac[0][1] = jac1[0][1] + jac2[0][1];
jac[1][0] = jac1[1][0] + jac2[1][0];
jac[1][1] = jac1[1][1] + jac2[1][1];
return jac;
Jacobian operator-(const Jacobian& jac) const {
return Jacobian{(*this)[0][0] - jac[0][0], (*this)[0][1] - jac[0][1], (*this)[1][0] - jac[1][0],
(*this)[1][1] - jac[1][1]};
}

double norm() const {
const Jacobian& jac = *this;
return sqrt(jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1] + jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1]);
Jacobian operator+(const Jacobian& jac) const {
return Jacobian{(*this)[0][0] + jac[0][0], (*this)[0][1] + jac[0][1], (*this)[1][0] + jac[1][0],
(*this)[1][1] + jac[1][1]};
}

Jacobian operator*(const Jacobian& jac2) const {
const Jacobian& jac1 = *this;
Jacobian jac;
jac[0][0] = jac1[0][0] * jac2[0][0] + jac1[0][1] * jac2[1][0];
jac[0][1] = jac1[0][0] * jac2[0][1] + jac1[0][1] * jac2[1][1];
jac[1][0] = jac1[1][0] * jac2[0][0] + jac1[1][1] * jac2[1][0];
jac[1][1] = jac1[1][0] * jac2[0][1] + jac1[1][1] * jac2[1][1];
return jac;
Jacobian operator*(double a) const {
return Jacobian{(*this)[0][0] * a, (*this)[0][1] * a, (*this)[1][0] * a, (*this)[1][1] * a};
}

double dx_dlon() const {
const Jacobian& jac = *this;
return jac[JDX][JDLON];
}
double dy_dlon() const {
const Jacobian& jac = *this;
return jac[JDY][JDLON];
}
double dx_dlat() const {
const Jacobian& jac = *this;
return jac[JDX][JDLAT];
}
double dy_dlat() const {
const Jacobian& jac = *this;
return jac[JDY][JDLAT];
Jacobian operator*(const Jacobian& jac) const {
return Jacobian{(*this)[0][0] * jac[0][0] + (*this)[0][1] * jac[1][0],
(*this)[0][0] * jac[0][1] + (*this)[0][1] * jac[1][1],
(*this)[1][0] * jac[0][0] + (*this)[1][1] * jac[1][0],
(*this)[1][0] * jac[0][1] + (*this)[1][1] * jac[1][1]};
}

double dlon_dx() const {
const Jacobian& jac = *this;
return jac[JDLON][JDX];
Point2 operator*(const Point2& x) const {
return Point2{(*this)[0][0] * x[0] + (*this)[0][1] * x[1], (*this)[1][0] * x[0] + (*this)[1][1] * x[1]};
}
double dlon_dy() const {
const Jacobian& jac = *this;
return jac[JDLON][JDY];

double norm() const {
return std::sqrt((*this)[0][0] * (*this)[0][0] + (*this)[0][1] * (*this)[0][1] +
(*this)[1][0] * (*this)[1][0] + (*this)[1][1] * (*this)[1][1]);
}
double dlat_dx() const {
const Jacobian& jac = *this;
return jac[JDLAT][JDX];

double determinant() const { return (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0]; }

Jacobian inverse() const {
return Jacobian{(*this)[1][1], -(*this)[0][1], -(*this)[1][0], (*this)[0][0]} * (1. / determinant());
}
double dlat_dy() const {
const Jacobian& jac = *this;
return jac[JDLAT][JDY];

Jacobian transpose() const { return Jacobian{(*this)[1][1], (*this)[0][1], (*this)[1][0], (*this)[0][0]}; }

double dx_dlon() const { return (*this)[JDX][JDLON]; }
double dy_dlon() const { return (*this)[JDY][JDLON]; }
double dx_dlat() const { return (*this)[JDX][JDLAT]; }
double dy_dlat() const { return (*this)[JDY][JDLAT]; }

double dlon_dx() const { return (*this)[JDLON][JDX]; }
double dlon_dy() const { return (*this)[JDLON][JDY]; }
double dlat_dx() const { return (*this)[JDLAT][JDX]; }
double dlat_dy() const { return (*this)[JDLAT][JDY]; }

friend std::ostream& operator<<(std::ostream& os, const Jacobian& jac) {
os << jac[0][0] << " " << jac[0][1] << "\n" << jac[1][0] << " " << jac[1][1];
return os;
}

private:
Expand All @@ -147,6 +121,7 @@ class ProjectionImpl : public util::Object {
};
};


public:
static const ProjectionImpl* create(const eckit::Parametrisation& p);
static const ProjectionImpl* create(const std::string& type, const eckit::Parametrisation& p);
Expand Down
9 changes: 7 additions & 2 deletions src/atlas/util/PolygonXY.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -101,9 +101,10 @@ bool PolygonXY::contains(const Point2& P) const {
return dx * dx + dy * dy;
};

constexpr double eps = 1.e-10;
// check first bounding box
if (coordinatesMax_.y() < P.y() || P.y() < coordinatesMin_.y() || coordinatesMax_.x() < P.x() ||
P.x() < coordinatesMin_.x()) {
if (coordinatesMax_.y() + eps < P.y() || P.y() < coordinatesMin_.y() - eps || coordinatesMax_.x() + eps < P.x() ||
P.x() < coordinatesMin_.x() - eps) {
return false;
}

Expand Down Expand Up @@ -137,6 +138,10 @@ bool PolygonXY::contains(const Point2& P) const {

if (APB != BPA) {
const double side = cross_product_analog(P, A, B);
bool on_edge = std::abs(side) < eps;
if (on_edge) {
return true;
}
if (APB && side > 0) {
++wn;
}
Expand Down
2 changes: 1 addition & 1 deletion src/atlas/util/detail/KDTree.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,7 @@ class KDTree_eckit : public KDTreeBase<PayloadT, PointT> {
}

idx_t size() const override {
#if ATLAS_ECKIT_VERSION_INT >= 11302 // v1.13.2
#if ATLAS_ECKIT_VERSION_AT_LEAST(1, 13, 2)
return static_cast<idx_t>(tree_->size());
#else
// Assume ECKIT-515 not implemented.
Expand Down
55 changes: 55 additions & 0 deletions src/tests/projection/test_jacobian.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,6 +163,61 @@ CASE("test_lonlat") {
doTaylorTest(StructuredGrid("N16"), 1e-9);
}

CASE("test_constructors") {
Projection::Jacobian a = {};
Projection::Jacobian b = {1., 2., 3., 4.};
Projection::Jacobian c = {{5., 6.}, {7., 8.}};

Log::info() << "a =" << std::endl;
Log::info() << a << std::endl;
Log::info() << "b =" << std::endl;
Log::info() << b << std::endl;
Log::info() << "c =" << std::endl;
Log::info() << c << std::endl;

EXPECT_EQ(a[0][0], 0.);
EXPECT_EQ(a[0][1], 0.);
EXPECT_EQ(a[1][0], 0.);
EXPECT_EQ(a[1][1], 0.);

EXPECT_EQ(b[0][0], 1.);
EXPECT_EQ(b[0][1], 2.);
EXPECT_EQ(b[1][0], 3.);
EXPECT_EQ(b[1][1], 4.);

EXPECT_EQ(c[0][0], 5.);
EXPECT_EQ(c[0][1], 6.);
EXPECT_EQ(c[1][0], 7.);
EXPECT_EQ(c[1][1], 8.);
}

CASE("test_multiplication") {
const Projection::Jacobian A = {{0., -1.}, {1., 0.}};
const Point2 x = {2., 1.};

const auto Ax = A * x;
const auto AinvAx = A.inverse() * Ax;
const auto Atimes2 = A * 2.;

Log::info() << "A =" << std::endl;
Log::info() << A << std::endl;
Log::info() << "2A =" << std::endl;
Log::info() << Atimes2 << std::endl;
Log::info() << "x =" << std::endl;
Log::info() << x << std::endl;
Log::info() << "Ax =" << std::endl;
Log::info() << Ax << std::endl;
Log::info() << "A^-1 Ax =" << std::endl;
Log::info() << AinvAx << std::endl;

EXPECT_EQ(Atimes2[0][0], 0.);
EXPECT_EQ(Atimes2[0][1], -2.);
EXPECT_EQ(Atimes2[1][0], 2.);
EXPECT_EQ(Atimes2[1][1], 0.);
EXPECT_EQ(Ax, Point2(-1., 2.));
EXPECT_EQ(AinvAx, x);
}


//-----------------------------------------------------------------------------

Expand Down
2 changes: 1 addition & 1 deletion src/tests/util/test_kdtree.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,7 @@ static const IndexKDTree& search() {
} // namespace
//------------------------------------------------------------------------------------------------

static bool ECKIT_515_implemented = (ATLAS_ECKIT_VERSION_INT >= 11302); // version 1.13.2
static bool ECKIT_515_implemented = ATLAS_ECKIT_VERSION_AT_LEAST(1, 13, 2);
// --> implements eckit::KDTree::size() and eckit::KDTree::empty()

CASE("test kdtree") {
Expand Down