Skip to content

Commit

Permalink
Added methods to projection::Jacobian.
Browse files Browse the repository at this point in the history
  • Loading branch information
@odlomax
odlomax committed Feb 17, 2022
1 parent 7afd801 commit 39477a6
Show file tree
Hide file tree
Showing 2 changed files with 126 additions and 61 deletions.
126 changes: 65 additions & 61 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,97 @@ class ProjectionImpl : public util::Object {
using Spec = atlas::util::Config;
class Jacobian : public std::array<std::array<double, 2>, 2> {
public:

using std::array<std::array<double, 2>, 2>::array;

static Jacobian identity() {
Jacobian id;
id[0] = {1.0, 0.0};
id[1] = {0.0, 1.0};
return id;
return Jacobian{1., 0., 0., 1.};
}

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;
};
Jacobian(double j00, double j01, double j10, double j11) :
array{std::array<double, 2>{j00, j01},
std::array<double, 2>{j10, j11}} {}

Jacobian transpose() const {
Jacobian tra = *this;
std::swap(tra[0][1], tra[1][0]);
return tra;
};
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]};
}

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]};
}

Jacobian operator*(double a) const {
return Jacobian{(*this)[0][0] * a, (*this)[0][1] * a,
(*this)[1][0] * a, (*this)[1][1] * a};
}

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]};
}

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 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]);
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 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());
}

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 transpose() const {
return Jacobian{(*this)[1][1], (*this)[0][1], (*this)[1][0], (*this)[0][0]};
}

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

double dlon_dx() const {
const Jacobian& jac = *this;
return jac[JDLON][JDX];
return (*this)[JDLON][JDX];
}
double dlon_dy() const {
const Jacobian& jac = *this;
return jac[JDLON][JDY];
return (*this)[JDLON][JDY];
}
double dlat_dx() const {
const Jacobian& jac = *this;
return jac[JDLAT][JDX];
return (*this)[JDLAT][JDX];
}
double dlat_dy() const {
const Jacobian& jac = *this;
return jac[JDLAT][JDY];
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 +150,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
61 changes: 61 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,67 @@ 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 39477a6

Please sign in to comment.