Skip to content

Commit

Permalink
Code cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
lucacarniato committed Nov 9, 2023
1 parent 077613a commit b914d45
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 117 deletions.
3 changes: 3 additions & 0 deletions libs/MeshKernel/include/MeshKernel/HessianCalculator.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -43,13 +43,15 @@ namespace meshkernel
/// @brief Intended to be the computation of the Hessian
static void Compute(const std::vector<Sample>& rawSamplePoints,
const Projection projection,
UInt numberOfSmoothingIterations,
const UInt numX,
const UInt numY,
Hessian& hessian);

/// @brief Intended to be the computation of the Hessian
static void Compute(const std::vector<Sample>& rawSamplePoints,
const Projection projection,
UInt numberOfSmoothingIterations,
const UInt numX,
const UInt numY,
std::vector<Sample>& hessianSamples);
Expand Down Expand Up @@ -105,6 +107,7 @@ namespace meshkernel
/// From (prepare_samplehessian.f90)
static void PrepareSampleForHessian(const std::vector<Sample>& samplePoints,
const Projection projection,
UInt numberOfSmoothingIterations,
Hessian& hessian);
};

Expand Down
30 changes: 11 additions & 19 deletions libs/MeshKernel/src/HessianCalculator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -287,8 +287,6 @@ void meshkernel::HessianCalculator::ComputeHessian(const std::vector<Sample>& sa

for (UInt i = 1; i < hessian.size(1) - 1; ++i)
{
// double af = static_cast<double>(i - 1) / static_cast<double>(std::max(hessian.size(1) - 3, 1));

for (UInt j = 1; j < hessian.size(2) - 1; ++j)
{

Expand Down Expand Up @@ -341,56 +339,50 @@ void meshkernel::HessianCalculator::ComputeHessian(const std::vector<Sample>& sa
VV(1, 0) = (gradientiR[1] * SniR[0] - gradientiL[1] * SniL[0] + gradientjR[1] * SnjR[0] - gradientjL[1] * SnjL[0]) * areaInv;
VV(1, 1) = (gradientiR[1] * SniR[1] - gradientiL[1] * SniL[1] + gradientjR[1] * SnjR[1] - gradientjL[1] * SnjL[1]) * areaInv;

// Eigendecompostion
// Eigen decompostion

Check failure on line 342 in libs/MeshKernel/src/HessianCalculator.cpp

View workflow job for this annotation

GitHub Actions / check

decompostion ==> decomposition
Eigen::EigenSolver<Eigen::Matrix2d> eigensolver(VV);

const Eigen::EigenSolver<Eigen::Matrix2d>::EigenvectorsType& eigenvectors = eigensolver.eigenvectors();
Eigen::EigenSolver<Eigen::Matrix2d>::EigenvalueType eigenvalues = eigensolver.eigenvalues();

UInt k = std::abs(eigenvalues[0].real()) > std::abs(eigenvalues[1].real()) ? 0U : 1u;

hessian(1, i, j) = eigenvectors(0, k).real();
hessian(2, i, j) = eigenvectors(1, k).real();
hessian(3, i, j) = eigenvalues[k].real() * area;
hessian(4, i, j) = -(eigenvectors(k, 0).real() * zx - eigenvectors(k, 1).real() * zy) / (eigenvalues[k].real() + 1.0e-8);
const UInt k = std::abs(eigenvalues[0].real()) > std::abs(eigenvalues[1].real()) ? 0U : 1u;
hessian(1, i, j) = eigenvalues[k].real() * area;
}
}
}

void meshkernel::HessianCalculator::PrepareSampleForHessian(const std::vector<Sample>& samplePoints,
const Projection projection,
UInt numberOfSmoothingIterations,
Hessian& hessian)
{

// TODO pass in Compute function
UInt numberOfSmoothingIterations = 0;

SmoothSamples(samplePoints, numberOfSmoothingIterations, hessian);
ComputeHessian(samplePoints, projection, hessian);
}

void meshkernel::HessianCalculator::Compute(const std::vector<Sample>& rawSamplePoints,
const Projection projection,
UInt numberOfSmoothingIterations,
const UInt numX,
const UInt numY,
Hessian& hessian)
{
std::vector<Sample> samplePoints(rawSamplePoints);

hessian.resize(5, numX, numY);
PrepareSampleForHessian(samplePoints, projection, hessian);
hessian.resize(2, numX, numY);
PrepareSampleForHessian(samplePoints, projection, numberOfSmoothingIterations, hessian);
}

void meshkernel::HessianCalculator::Compute(const std::vector<Sample>& rawSamplePoints,
const Projection projection,
UInt numberOfSmoothingIterations,
const UInt numX,
const UInt numY,
std::vector<Sample>& hessianSamples)
{
hessianSamples = rawSamplePoints;

Hessian hessian(5, numY, numX);
PrepareSampleForHessian(hessianSamples, projection, hessian);
Hessian hessian(2, numY, numX);
PrepareSampleForHessian(hessianSamples, projection, numberOfSmoothingIterations, hessian);

UInt count = 0;

Expand All @@ -400,7 +392,7 @@ void meshkernel::HessianCalculator::Compute(const std::vector<Sample>& rawSample
{

count = numY * i + j;
hessianSamples[count].value = hessian(3, j, i);
hessianSamples[count].value = hessian(1, j, i);
}
}
}
107 changes: 9 additions & 98 deletions libs/MeshKernel/tests/src/HessianTests.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ namespace mk = meshkernel;
std::vector<std::vector<mk::Point>> GenerateGridPoints(const mk::UInt rows, const mk::UInt cols)
{

std::vector<std::vector<mk::Point>> meshPoints(cols, std::vector<mk::Point>(rows));
std::vector meshPoints(cols, std::vector<mk::Point>(rows));

mk::Point origin = mk::Point(0.0, 0.0);

Expand All @@ -48,90 +48,6 @@ std::vector<std::vector<mk::Point>> GenerateGridPoints(const mk::UInt rows, cons
return meshPoints;
}

#if 0

TEST(HessianTests, TidyPointsNoDuplicatesTest)
{
// Test that the TidySamples does not change a data set that has no duplicates.

const mk::UInt size = 100;

std::vector<mk::Point> samplePoints(size);
std::vector<mk::Point> generated;
std::vector<double> sampleData(size);

std::random_device rand_dev;
std::mt19937 generator(rand_dev());
std::uniform_real_distribution<double> distribution(0.0, 100.0);

for (mk::UInt i = 0; i < size; ++i)
{
samplePoints[i] = mk::Point(distribution(generator), distribution(generator));
sampleData[i] = distribution(generator);
}

std::vector<mk::Point> expectedSamplePoints = samplePoints;
std::vector<double> expectedSampleData = sampleData;

mk::RidgeRefinement ridgeRefinement;
ridgeRefinement.TidySamples(samplePoints, sampleData);

for (mk::UInt i = 0; i < size; ++i)
{
EXPECT_EQ(expectedSamplePoints[i].x, samplePoints[i].x);
EXPECT_EQ(expectedSamplePoints[i].y, samplePoints[i].y);
EXPECT_EQ(expectedSampleData[i], sampleData[i]);
}
}

TEST(HessianTests, TidyPointsTest)
{

const mk::UInt size = 100;

std::vector<mk::Point> samplePoints(size);
std::vector<mk::Point> expectedSamplePoints;
std::vector<mk::Point> generated;
std::vector<double> sampleData(size);

const double invMax = 1.0 / static_cast<double>(RAND_MAX);

for (mk::UInt i = 0; i < size; ++i)
{
// Use the simple pseudo random number generator to have repeatable values between tests when debugging.
samplePoints[i] = mk::Point(rand() * invMax, rand() * invMax);
sampleData[i] = rand() * invMax;

if (i > 0 && i % 10 == 0)
{
samplePoints[i - 1] = samplePoints[i];
sampleData[i - 1] = sampleData[i];
}
else if (i > 10 && i % 21 == 0)
{
samplePoints[i - 10] = samplePoints[i];
sampleData[i - 10] = sampleData[i];
}
else
{
expectedSamplePoints.push_back(samplePoints[i]);
}
}

mk::RidgeRefinement ridgeRefinement;
ridgeRefinement.TidySamples(samplePoints, sampleData);

EXPECT_EQ(expectedSamplePoints.size(), samplePoints.size());

#if 0
for (const auto& point : expectedSamplePoints)
{
EXPECT_TRUE(std::find(samplePoints.begin(), samplePoints.end(), point) != samplePoints.end());
}
#endif
}

#endif
std::tuple<std::vector<meshkernel::Sample>, std::vector<std::vector<double>>> generateSampleData(meshkernel::UInt nx, meshkernel::UInt ny, double deltaX, double deltaY)
{
meshkernel::UInt start = 0;
Expand All @@ -143,14 +59,8 @@ std::tuple<std::vector<meshkernel::Sample>, std::vector<std::vector<double>>> ge
[[maybe_unused]] double centreX = (static_cast<double>((nx - 1) / 2) * deltaX);
[[maybe_unused]] double centreY = (static_cast<double>((ny - 1) / 2) * deltaY);

std::cout << "centre: " << centreX << " " << centreY << std::endl;

[[maybe_unused]] double xInit = 0.0; // 0.5 * deltaX;
[[maybe_unused]] double yInit = 0.0; // 0.5 * deltaY;

[[maybe_unused]] double scale = (ny / 4.0) * deltaY;
[[maybe_unused]] double x = xInit; // 0.0;
[[maybe_unused]] double y = yInit; // 0.0;

[[maybe_unused]] double r = (nx / 5) * deltaX;
[[maybe_unused]] double maxx = (nx - 1) * deltaX;

Expand All @@ -159,10 +69,10 @@ std::tuple<std::vector<meshkernel::Sample>, std::vector<std::vector<double>>> ge

for (meshkernel::UInt j = start; j < nx; ++j)
{
y = deltaY * i;
x = deltaX * j;
const double y = deltaY * i;
const double x = deltaX * j;

#if 1
#if 0
// Gaussian bump, in the centre of the grid
double centre = (x - centreX) * (x - centreX) + (y - centreY) * (y - centreY);
double sample = 100.0 * std::exp(-0.025 * centre);
Expand Down Expand Up @@ -194,8 +104,8 @@ std::tuple<std::vector<meshkernel::Sample>, std::vector<std::vector<double>>> ge
{
for (int i = ny - 1; i >= 0; --i)
{
y = deltaY * i;
x = deltaX * j;
const double y = deltaY * i;
const double x = deltaX * j;
sampleData[count] = {x, y, sampleDataMatrix[(ny - 1) - i][j]};
count++;
}
Expand Down Expand Up @@ -234,7 +144,8 @@ TEST(HessianTests, CheckHessian)
WriteSampleFileToAsc(sampleDataMatrix, sampleDeltaX, filePath);

std::vector<meshkernel::Sample> samples;
meshkernel::HessianCalculator::Compute(sampleData, mesh->m_projection, sampleNx, sampleNy, samples);
meshkernel::UInt numberOfSmoothingIterations = 0;
meshkernel::HessianCalculator::Compute(sampleData, mesh->m_projection, numberOfSmoothingIterations, sampleNx, sampleNy, samples);

const auto interpolator = std::make_shared<meshkernel::AveragingInterpolation>(*mesh,
samples,
Expand Down

0 comments on commit b914d45

Please sign in to comment.