Skip to content

Commit

Permalink
Regular Field tests working
Browse files Browse the repository at this point in the history
  • Loading branch information
shutsch committed Feb 21, 2024
1 parent 9b68ee8 commit 7b76fb1
Show file tree
Hide file tree
Showing 5 changed files with 243 additions and 120 deletions.
88 changes: 57 additions & 31 deletions c_library/test/regular/test_helix.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "RegularModels.h"
#include "../test_helpers.h"


vector model(double x, double y, number ampx, number ampy, number ampz, double rmin, double rmax) {
const double phi = std::atan2(y, x); // azimuthal angle in cylindrical coordinates
const double r = std::sqrt(x * x + y * y); // radius in cylindrical coordinates
Expand All @@ -25,37 +26,38 @@ vector model(double x, double y, number ampx, number ampy, number ampz, double r

#if autodiff_FOUND

Eigen::Matrix3d model_derivative(double x, double y, number ampx, number ampy, number ampz, double rmin, double rmax) {
Eigen::MatrixXd model_derivative(double x, double y, number ampx, number ampy, number ampz, double rmin, double rmax) {
const double phi = std::atan2(y, x); // azimuthal angle in cylindrical coordinates
const double r = std::sqrt(x * x + y * y); // radius in cylindrical coordinates
Eigen::Matrix3d b_deriv;
if ((r > rmin) && (r < rmax))
{
b_deriv(0, 0) = std::cos(phi);
b_deriv(0, 1) = 0.;
b_deriv(0, 2) = 0.;
b_deriv(1, 0) = 0.;
b_deriv(1, 1) = std::sin(phi);
b_deriv(1, 2) = 0.;
b_deriv(2, 0) = 0;
b_deriv(2, 1) = 0.;
b_deriv(2, 2) = 1.;
}
if ((r < rmin) || (r > rmax)) return Eigen::MatrixXd::Zero(3, 3);
Eigen::MatrixXd b_deriv;

b_deriv(0, 0) = std::cos(phi);
b_deriv(0, 1) = 0.;
b_deriv(0, 2) = 0.;

b_deriv(1, 0) = 0.;
b_deriv(1, 1) = std::sin(phi);
b_deriv(1, 2) = 0.;

b_deriv(2, 0) = 0;
b_deriv(2, 1) = 0.;
b_deriv(2, 2) = 1.;

return b_deriv;
}
# endif

TEST_CASE("HelixMagneticField") {

// constructor s
HelixMagneticField helix = HelixMagneticField();
// DEFINITIONS

// define positions in Galactic cartesian coordinates (units are kpc)

std::array<double, 3> origin{{0., 0., 0.}};
std::array<double, 3> smaller_rmin{{.9, -0.87, 0.65}};
std::array<double, 3> smaller_rmin{{.09, -0.47, 0.65}};
std::array<double, 3> bigger_rmin_smaller_rmax{{-17., 13.4, 9.5}};
std::array<double, 3> bigger_rmax{{24., -23.13, 19.5}};

Expand All @@ -78,21 +80,29 @@ TEST_CASE("HelixMagneticField") {
size_t arr_sz = 4*3*2;

// define positions of a irregular grid in Galactic cartesian coordinates (units are kpc)
const std::vector<double> grid_x {{2. , -4. , 0. , 1., 0.4}};
const std::vector<double> grid_y {{4. , 6. , -0.1, 0., 0.2}};
const std::vector<double> grid_z {{-0.2, 0.8, 0.2, 0., 1.}};

const std::vector<double> grid_x {{0., 0. , 0. , 0. , 0. , 0.3, -17.2,
0. , 0. , 0. , 0. , -35.032, 2.2, -6. , .234, -0.00424, -2.2, -2.32 , .234,
0.43 , 8.23 , 3. , 1., -0.001 , -5.342, -5.23, -2.432 }};
const std::vector<double> grid_y {{0., 0. , 0. , 11.5, -7.2, 0. , 0. ,
-0.0032, 23.353, -6.1 , 4. , 0. , 0. , 0. , 0. , -0.00132, 2.2 , -6.42, .234,
6.123, 8.1 , -3. , -1., 0.0033 , 0.45 , -30.2, -1.543 }};
const std::vector<double> grid_z {{0., 3.1, -0.01, 0. , 0. , 0. , 0. ,
-5.32 , 0.012 , 34.02, -8.431, -5.132, 9.3 , -6.1, 4. , 0. , 0. , 0. , 0. ,
10.3 , -0.02 , 3. , -1., 0.000432, -15.43, 0.001, -7.123 }};

const int size_pos = grid_x.size();

std::array<double *, 3> default_irregular_grid;
default_irregular_grid[0] = new double[5];
default_irregular_grid[1] = new double[5];
default_irregular_grid[2] = new double[5];
default_irregular_grid[0] = new double[size_pos];
default_irregular_grid[1] = new double[size_pos];
default_irregular_grid[2] = new double[size_pos];

std::array<double *, 3> updated_irregular_grid;
updated_irregular_grid[0] = new double[5];
updated_irregular_grid[1] = new double[5];
updated_irregular_grid[2] = new double[5];
updated_irregular_grid[0] = new double[size_pos];
updated_irregular_grid[1] = new double[size_pos];
updated_irregular_grid[2] = new double[size_pos];

for (int i=0; i < 5; i++) {
for (int i=0; i < size_pos; i++) {
vector def_mod_at_grid_point = model(grid_x[0], grid_y[0], helix.ampx, helix.ampy, helix.ampz, helix.rmin, helix.rmax);
default_irregular_grid[0][i] = def_mod_at_grid_point[0].val();
default_irregular_grid[1][i] = def_mod_at_grid_point[1].val();
Expand All @@ -102,7 +112,7 @@ TEST_CASE("HelixMagneticField") {
updated_irregular_grid[1][i] = up_mod_at_grid_point[1].val();
updated_irregular_grid[2][i] = up_mod_at_grid_point[2].val();
}

// test at_position -- default parameters
REQUIRE_THAT(helix.at_position(origin[0], origin[1], origin[2]), EqualsVector(zeros));
REQUIRE_THAT(helix.at_position(smaller_rmin[0], smaller_rmin[1], smaller_rmin[2]), EqualsVector(zeros));
Expand All @@ -112,7 +122,19 @@ TEST_CASE("HelixMagneticField") {
// test on_grid -- default parameters
REQUIRE_THROWS_WITH(helix.on_grid(), "The class has not been initialized with a grid, hence on_grid can only be called with a grid provided.");
REQUIRE_NOTHROW(helix.on_grid(shape, refpoint, increment));
REQUIRE_THAT(helix.on_grid(grid_x, grid_y, grid_z), EqualsPointerArray(default_irregular_grid, 5));
REQUIRE_THAT(helix.on_grid(grid_x, grid_y, grid_z), EqualsPointerArray(default_irregular_grid, size_pos));

#if autodiff_FOUND
// test derivative
CAPTURE(helix.derivative(origin[0], origin[1], origin[2]), model_derivative(origin[0], origin[1], helix.ampx, helix.ampy, helix.ampz, helix.rmin, helix.rmax));
REQUIRE_THAT(helix.derivative(origin[0], origin[1], origin[2]), EqualsMatrix(model_derivative(origin[0], origin[1], helix.ampx, helix.ampy, helix.ampz, helix.rmin, helix.rmax)));

REQUIRE_THAT(helix.derivative(smaller_rmin[0], smaller_rmin[1], smaller_rmin[2]), EqualsMatrix(model_derivative(smaller_rmin[0], smaller_rmin[1], helix.ampx, helix.ampy, helix.ampz, helix.rmin, helix.rmax)));
REQUIRE_THAT(helix.derivative(bigger_rmin_smaller_rmax[0], bigger_rmin_smaller_rmax[1], bigger_rmin_smaller_rmax[2]), EqualsMatrix(model_derivative(bigger_rmin_smaller_rmax[0], bigger_rmin_smaller_rmax[1], helix.ampx, helix.ampy, helix.ampz, helix.rmin, helix.rmax)));
REQUIRE_THAT(helix.derivative(bigger_rmax[0], bigger_rmax[1], bigger_rmax[2]), EqualsMatrix(model_derivative(bigger_rmax[0], bigger_rmax[1], helix.ampx, helix.ampy, helix.ampz, helix.rmin, helix.rmax)));
#endif


// parameter update
helix.ampx = ampx;
helix.ampy = ampy;
Expand All @@ -124,6 +146,8 @@ TEST_CASE("HelixMagneticField") {
REQUIRE(helix.ampz == ampz);
REQUIRE(helix.rmin == rmin);
REQUIRE(helix.rmax == rmax);


// test at_position -- default parameters
REQUIRE_THAT(helix.at_position(origin[0], origin[1], origin[2]), EqualsVector(zeros));
REQUIRE_THAT(helix.at_position(smaller_rmin[0], smaller_rmin[1], smaller_rmin[2]), EqualsVector(model(bigger_rmin_smaller_rmax[0], bigger_rmin_smaller_rmax[1], ampx, ampy, ampz, rmin, rmax)));
Expand All @@ -132,7 +156,9 @@ TEST_CASE("HelixMagneticField") {
// test on_grid -- default parameters
REQUIRE_THROWS_WITH(helix.on_grid(), "The class has not been initialized with a grid, hence on_grid can only be called with a grid provided.");
REQUIRE_NOTHROW(helix.on_grid(shape, refpoint, increment));
REQUIRE_THAT(helix.on_grid(grid_x, grid_y, grid_z), EqualsPointerArray(updated_irregular_grid, 5));
REQUIRE_THAT(helix.on_grid(grid_x, grid_y, grid_z), EqualsPointerArray(updated_irregular_grid, size_pos));


#if autodiff_FOUND
// test derivative

Expand Down
108 changes: 74 additions & 34 deletions c_library/test/regular/test_regular.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include <cassert>
#include <iostream>
#include <vector>
#include <cmath>
#include <map>
#include <memory>

Expand All @@ -12,11 +13,22 @@


TEST_CASE("RegularMagneticField") {

const std::vector<double> grid_x {{2., 4., 0., 1., .4}};
const std::vector<double> grid_y {{4., 6., 0.1, 0., .2}};
const std::vector<double> grid_z {{-0.2, 0.8, 0.2, 0., 1.}};

//
const std::vector<double> positions_x {{0., 0. , 0. , 0. , 0. , 0.3, -17.2,
0. , 0. , 0. , 0. , -35.032, 2.2, -6. , .234, -0.00424, -2.2, -2.32 , .234,
0.43 , 8.23 , 3. , 1., -0.001 , -5.342, -5.23, -2.432 }};
const std::vector<double> positions_y {{0., 0. , 0. , 11.5, -7.2, 0. , 0. ,
-0.0032, 23.353, -6.1 , 4. , 0. , 0. , 0. , 0. , -0.00132, 2.2 , -6.42, .234,
6.123, 8.1 , -3. , -1., 0.0033 , 0.45 , -30.2, -1.543 }};
const std::vector<double> positions_z {{0., 3.1, -0.01, 0. , 0. , 0. , 0. ,
-5.32 , 0.012 , 34.02, -8.431, -5.132, 9.3 , -6.1, 4. , 0. , 0. , 0. , 0. ,
10.3 , -0.02 , 3. , -1., 0.000432, -15.43, 0.001, -7.123 }};

int size_pos = positions_x.size();

UNSCOPED_INFO("Test case broken!");
REQUIRE(size_pos == positions_y.size());
REQUIRE(size_pos == positions_z.size());

// Define a regular grid in Galactic cartesian coordinates (units are kpc)
const std::array<int, 3> shape {{4, 3, 2}};
Expand Down Expand Up @@ -56,49 +68,77 @@ TEST_CASE("RegularMagneticField") {
models_w_regular_constructor["WMAP"] = std::shared_ptr<WMAPMagneticField> (new WMAPMagneticField(shape, refpoint, increment));
models_w_regular_constructor["Archimedean spiral"] = std::shared_ptr<ArchimedeanMagneticField> (new ArchimedeanMagneticField(shape, refpoint, increment));

models_w_irregular_constructor["JF12"] = std::shared_ptr<JF12MagneticField> (new JF12MagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Jaffe"] = std::shared_ptr<JaffeMagneticField> (new JaffeMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Sun2007"] = std::shared_ptr<SunMagneticField> (new SunMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Han2018"] = std::shared_ptr<HanMagneticField> (new HanMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Pshirkov"] = std::shared_ptr<PshirkovMagneticField> (new PshirkovMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["HMR"] = std::shared_ptr<HMRMagneticField> (new HMRMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["TT"] = std::shared_ptr<TTMagneticField> (new TTMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["TF"] = std::shared_ptr<TFMagneticField> (new TFMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Fauvet"] = std::shared_ptr<FauvetMagneticField> (new FauvetMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Stanev"] = std::shared_ptr<StanevBSSMagneticField> (new StanevBSSMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["WMAP"] = std::shared_ptr<WMAPMagneticField> (new WMAPMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["Archimedean spiral"] = std::shared_ptr<ArchimedeanMagneticField> (new ArchimedeanMagneticField(grid_x, grid_y, grid_z));
models_w_irregular_constructor["JF12"] = std::shared_ptr<JF12MagneticField> (new JF12MagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Jaffe"] = std::shared_ptr<JaffeMagneticField> (new JaffeMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Sun2007"] = std::shared_ptr<SunMagneticField> (new SunMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Han2018"] = std::shared_ptr<HanMagneticField> (new HanMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Pshirkov"] = std::shared_ptr<PshirkovMagneticField> (new PshirkovMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["HMR"] = std::shared_ptr<HMRMagneticField> (new HMRMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["TT"] = std::shared_ptr<TTMagneticField> (new TTMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["TF"] = std::shared_ptr<TFMagneticField> (new TFMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Fauvet"] = std::shared_ptr<FauvetMagneticField> (new FauvetMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Stanev"] = std::shared_ptr<StanevBSSMagneticField> (new StanevBSSMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["WMAP"] = std::shared_ptr<WMAPMagneticField> (new WMAPMagneticField(positions_x, positions_y, positions_z));
models_w_irregular_constructor["Archimedean spiral"] = std::shared_ptr<ArchimedeanMagneticField> (new ArchimedeanMagneticField(positions_x, positions_y, positions_z));

auto empty_model_iter = models_w_empty_constructor.begin();
auto regular_model_iter = models_w_empty_constructor.begin();
auto irregular_model_iter = models_w_empty_constructor.begin();

while (empty_model_iter != models_w_empty_constructor.end()) {
REQUIRE_THROWS_WITH((*models_w_empty_constructor[empty_model_iter->first]).on_grid(), "The class has not been initialized with a grid, hence on_grid can only be called with a grid provided.");
auto pos_x_iter = positions_x.begin();
auto pos_y_iter = positions_y.begin();
auto pos_z_iter = positions_z.begin();
while (pos_x_iter != positions_x.end()) {
UNSCOPED_INFO("RegularModels testing: at_position call failed for: ");
CAPTURE(irregular_model_iter->first, *pos_x_iter, *pos_y_iter, *pos_z_iter);
REQUIRE_NOTHROW((*models_w_empty_constructor[empty_model_iter->first]).at_position(*pos_x_iter, *pos_y_iter, *pos_z_iter));
UNSCOPED_INFO("RegularModels testing: at_position contained NaN: ");
CAPTURE(irregular_model_iter->first, *pos_x_iter, *pos_y_iter, *pos_z_iter);
auto eval = (*models_w_empty_constructor[empty_model_iter->first]).at_position(*pos_x_iter, *pos_y_iter, *pos_z_iter);
REQUIRE_FALSE(containsNaN(eval));

#if autodiff_FOUND
// test derivative
UNSCOPED_INFO("RegularModels testing: derivative call failed for: ");
CAPTURE(irregular_model_iter->first, pos_x_iter, pos_y_iter, pos_z_iter);
REQUIRE_NOTHROW((*models_w_empty_constructor[empty_model_iter->first]).derivative(*pos_x_iter, *pos_y_iter, *pos_z_iter));
#endif

++pos_x_iter;
++pos_y_iter;
++pos_z_iter;
}

UNSCOPED_INFO("RegularModels testing: on grid failed for: ");
CAPTURE(irregular_model_iter->first);
//std::cout << irregular_model_iter->first << std::endl;

REQUIRE_THROWS_WITH((*models_w_empty_constructor[empty_model_iter->first]).on_grid(), "The class has not been initialized with a grid, hence on_grid can only be called with a grid provided.");
std::array<double*, 3> eval_irregular = (*models_w_irregular_constructor[irregular_model_iter->first]).
on_grid();
std::array<double*, 3> eval_regular = (*models_w_regular_constructor[regular_model_iter->first]).on_grid();
std::array<double*, 3> eval_no_grid_regular = (*models_w_empty_constructor[empty_model_iter->first]).on_grid(shape, refpoint, increment);
std::array<double*, 3> eval_no_grid_irregular = (*models_w_empty_constructor[empty_model_iter->first]).on_grid(grid_x, grid_y, grid_z);

REQUIRE_THAT(eval_regular, EqualsPointerArray(eval_no_grid_regular, 5));
REQUIRE_THAT(eval_irregular, EqualsPointerArray(eval_no_grid_irregular, 5));
std::array<double*, 3> eval_no_grid_irregular = (*models_w_empty_constructor[empty_model_iter->first]).on_grid(positions_x, positions_y, positions_z);

CAPTURE(eval_irregular[0][0], eval_irregular[1][0], eval_irregular[2][0],
eval_no_grid_irregular[0][0], eval_no_grid_irregular[1][0], eval_no_grid_irregular[2][0],
eval_irregular[0][1], eval_irregular[1][1], eval_irregular[2][1],
eval_no_grid_irregular[0][1], eval_no_grid_irregular[1][1], eval_no_grid_irregular[2][1],
eval_irregular[0][2], eval_irregular[1][2], eval_irregular[2][2],
eval_no_grid_irregular[0][2], eval_no_grid_irregular[1][2], eval_no_grid_irregular[2][2],
eval_irregular[0][3], eval_irregular[1][3], eval_irregular[2][3],
eval_no_grid_irregular[0][3], eval_no_grid_irregular[1][3], eval_no_grid_irregular[2][3],
eval_irregular[0][4], eval_irregular[1][4], eval_irregular[2][4],
eval_no_grid_irregular[0][4], eval_no_grid_irregular[1][4], eval_no_grid_irregular[2][4]
);

REQUIRE_THAT(eval_regular, EqualsPointerArray(eval_no_grid_regular, 4*3*2));
REQUIRE_THAT(eval_irregular, EqualsPointerArray(eval_no_grid_irregular, size_pos));

REQUIRE_NOTHROW((*models_w_empty_constructor[empty_model_iter->first]).at_position(1., 2., -3.1));
REQUIRE_NOTHROW((*models_w_irregular_constructor[irregular_model_iter->first]).at_position(0., 0., 0.));
REQUIRE_NOTHROW((*models_w_regular_constructor[regular_model_iter->first]).at_position(-10., 23., -33.1));

#if autodiff_FOUND
// test derivative
REQUIRE_NOTHROW((*models_w_empty_constructor[empty_model_iter->first]).derivative(1., 2., -3.1));
REQUIRE_NOTHROW((*models_w_irregular_constructor[irregular_model_iter->first]).derivative(0., 0., 0.));
REQUIRE_NOTHROW((*models_w_regular_constructor[regular_model_iter->first]).derivative(-10., 23., -33.1));
#endif

++empty_model_iter;
++regular_model_iter;
++irregular_model_iter;
}
}

}
Loading

0 comments on commit 7b76fb1

Please sign in to comment.