Skip to content

Commit

Permalink
Added methods to projection::Jacobian (#93)
Browse files Browse the repository at this point in the history
  • Loading branch information
@odlomax @wdeconinck
odlomax authored and wdeconinck committed Feb 17, 2022
1 parent 09cde23 commit e4b12c7
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 76 deletions.
127 changes: 51 additions & 76 deletions src/atlas/projection/detail/ProjectionImpl.h
View file
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
55 changes: 55 additions & 0 deletions src/tests/projection/test_jacobian.cc
View file
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

0 comments on commit e4b12c7

Please sign in to comment.