Skip to content

Conversation

@J-Hizzle
Copy link
Contributor

Solves #47, albeit with linear interpolation as a first, crude way. We might eventually implement a better interpolation scheme (e.g. cubic spline interpolation) here.

@J-Hizzle
Copy link
Contributor Author

I suspect that the error arising in the SPARTIAN algorithm has something to do with the fact that the interpolation alorithm produces additional points of the thermodynamic force right-hand of the center that are not there on the left-hand side.
This might be remedied by setting all points of the thermodynamic force to zero outside the application region, which would anyways make more sense, as the thermodynamic force seen by the particles IS ZERO outside the application region.

@J-Hizzle
Copy link
Contributor Author

I suspect that the error arising in the SPARTIAN algorithm has something to do with the fact that the interpolation alorithm produces additional points of the thermodynamic force right-hand of the center that are not there on the left-hand side. This might be remedied by setting all points of the thermodynamic force to zero outside the application region, which would anyways make more sense, as the thermodynamic force seen by the particles IS ZERO outside the application region.

Different from what I first thought, the problem seems to be with CUDA, so I assume that I'm using Kokkos wrong somewhere.

@J-Hizzle
Copy link
Contributor Author

I suspect that the error arising in the SPARTIAN algorithm has something to do with the fact that the interpolation alorithm produces additional points of the thermodynamic force right-hand of the center that are not there on the left-hand side. This might be remedied by setting all points of the thermodynamic force to zero outside the application region, which would anyways make more sense, as the thermodynamic force seen by the particles IS ZERO outside the application region.

Different from what I first thought, the problem seems to be with CUDA, so I assume that I'm using Kokkos wrong somewhere.

This problem arose due to the out-of-bounds access to the "grid" and "gradient" views, which happened because the result of the getRightBin function for the zeroth entry of the input grid taken -1 results in a value of -1, which is then taken to access the -1th element of either of the two views and hence produces the error.

I fixed it by manually setting the interpolated output data points to zero when they are beyond the input grid points.

One could also implement periodic boundary conditions here, but since the density profile also seems to be somewhat artificial close to the periodic boundaries, it's best practice not to use those values in the thermodynamic force anyways.

@J-Hizzle J-Hizzle force-pushed the implement-smooth-interpolation-of-thermodynamic-force branch from ef5a5d0 to bb5e16c Compare August 29, 2025 14:17
@J-Hizzle
Copy link
Contributor Author

J-Hizzle commented Sep 4, 2025

As indicated by the commit messages, I switched the interpolation strategies back and forth between a straight-forward linear interpolation and an "interpolation" based on overlapping density bins.
The problem with the overlapping bins seems to be that there are substantial oscillations in the measured density profiles at neighboring bins that inevitably affect the stability of the thermodynamic force calculation. As those oscillations stem from the parts of the bins that do NOT overlap (the parts that do overlap yield the same density after all), this creates correlations of the measured densities with statistics at positions far away from the bin's center. I'm not entirely sure whether this is even physically reasonable, it sure sounds shaky.
Hence, I reverted back to linear interpolation.

@J-Hizzle
Copy link
Contributor Author

Should I resolve the merge conflicts with the new main?

@XzzX
Copy link
Owner

XzzX commented Nov 27, 2025

Should I resolve the merge conflicts with the new main?

Yes and no. I prefer having smaller PRs with clear scope. Let's start by putting the changes to ApplicationRegion and interpolation into their own PRs. They should be independent of the rest. We will then continue to split this PR into pieces.

Copy link
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR implements smooth interpolation of thermodynamic forces to address issue #47. The implementation uses linear interpolation as an initial approach, with the intention to potentially upgrade to cubic spline interpolation in the future. The changes introduce the ability to use different grid spacings for density sampling and force application, allowing for more flexible thermodynamic force calculations.

Key Changes:

  • Added new interpolation utilities (lerp, findRightBin, interpolate, constrainToApplicationRegion) to enable smooth force field interpolation between different grids
  • Extended MultiHistogram to support custom bin sizes independent of the histogram range
  • Modified the density profile calculation to support overlapping bins with configurable bin widths for smoother sampling

Reviewed changes

Copilot reviewed 19 out of 19 changed files in this pull request and generated 11 comments.

Show a summary per file
File Description
mrmd/util/interpolation.hpp New utility header defining linear interpolation functions and grid search functionality
mrmd/util/interpolation.cpp Implementation of histogram interpolation and application region constraining
mrmd/data/MultiHistogram.hpp Added constructor overload accepting explicit binSize parameter and device grid creation method
mrmd/data/MultiHistogram.cpp Implemented device-side grid creation for use in Kokkos kernels
mrmd/data/MultiHistogram.test.cpp Added test for consistency between bin-to-position conversions
mrmd/action/ThermodynamicForce.hpp Extended with new constructors distinguishing grid spacing from bin width; added ApplicationRegion-aware update method
mrmd/action/ThermodynamicForce.cpp Implemented interpolation in update methods and added multiple constructor overloads
mrmd/analysis/AxialDensityProfile.hpp Added binWidth parameter to support overlapping bins
mrmd/analysis/AxialDensityProfile.cpp Reimplemented binning algorithm to support variable bin widths with overlapping regions
mrmd/analysis/AxialDensityProfile.test.cpp Added test for overlapping bin functionality
mrmd/util/ApplicationRegion.hpp Changed boundary condition from inclusive to exclusive; added getter methods
mrmd/util/ApplicationRegion.test.cpp Updated tests to match new exclusive boundary behavior
mrmd/io/RestoreThermoForce.hpp Added optional requestedDensityBinWidth parameter
mrmd/io/RestoreThermoForce.cpp Refactored to support variable density bin widths and renamed variables for clarity
mrmd/io/ThermoForce.test.cpp Updated test to use new gridSpacing terminology
pyMRMD/analysis.cpp Updated Python binding to include new binWidth parameter
pyMRMD/action.cpp Updated Python binding for ThermodynamicForce update method signature
examples/DensityFluctuations/DensityFluctuations.cpp Added densityBinWidth parameter to density profile calculation
examples/BinaryLennardJones/SPARTIAN.cpp Updated to pass binSize to density profile calculation
Comments suppressed due to low confidence (1)

mrmd/data/MultiHistogram.hpp:68

  • The MultiHistogram copy constructor doesn't preserve the custom binSize. The copy constructor (line 63-68) delegates to the constructor that recalculates binSize as (max - min) / numBins. This means if the original histogram was created with a custom binSize (using the 6-parameter constructor), the copy will not preserve this custom bin size. Consider updating the copy constructor to use the 6-parameter constructor:
MultiHistogram(const std::string& label, const MultiHistogram& histogram)
    : MultiHistogram(
          label, histogram.min, histogram.max, histogram.numBins, histogram.binSize, histogram.numHistograms)
{
    Kokkos::deep_copy(data, histogram.data);
}
    MultiHistogram(const std::string& label, const MultiHistogram& histogram)
        : MultiHistogram(
              label, histogram.min, histogram.max, histogram.numBins, histogram.numHistograms)
    {
        Kokkos::deep_copy(data, histogram.data);
    }

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

You can also share your feedback on Copilot code review for a chance to win a $100 gift card. Take the survey.

data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid)
{
data::MultiHistogram output(
"interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms);
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The output MultiHistogram created on line 23-24 doesn't preserve the input's binSize. Similar to the copy constructor issue, this creates an output histogram that recalculates binSize as (max - min) / numBins, which may not match the original grid spacing when the input has a custom bin size. Since the output is meant to match a specific target grid, this might be intentional, but it should be verified that the output's binSize correctly matches the spacing of the provided grid parameter.

Suggested change
"interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms);
"interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms);
// Set output binSize to match grid spacing
if (grid.extent(0) > 1)
{
output.binSize = (grid(grid.extent(0) - 1) - grid(0)) / (grid.extent(0) - 1);
output.inverseBinSize = 1.0_r / output.binSize;
}
else
{
output.binSize = 0.0_r;
output.inverseBinSize = 0.0_r;
}

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting point, I'll have to look into it.

Comment on lines +45 to +48
output.data(binIdx, histogramIdx) =
lerp(inputDataLeft,
inputDataRight,
(grid(binIdx) - inputGrid(leftBin)) * input.inverseBinSize);
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The interpolation factor calculation may be incorrect. On line 48, the factor is calculated as (grid(binIdx) - inputGrid(leftBin)) * input.inverseBinSize. This assumes the distance between leftBin and rightBin equals input.binSize. However, findRightBin returns the first bin where grid(idx) >= value, so the distance between inputGrid(leftBin) and inputGrid(rightBin) is actually input.binSize, making the formula correct. However, this should be verified through testing. Consider adding a check or assertion that rightBin == leftBin + 1 to ensure the bins are adjacent.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will have to look into it.

Comment on lines +35 to +56
int binNumForce = 0;
int histNum = 0;
real_t grid0 = 0;
real_t grid1 = 0;

std::ifstream fileThermoForce(filename);
std::getline(fileThermoForce, line);
std::stringstream gridLineStream(line);
while (gridLineStream >> word)
{
if (binNum == 0)
{
grid0 = std::stod(word);
}
if (binNum == 1)
{
grid1 = std::stod(word);
}
binNum++;
binNumForce++;
}
MRMD_HOST_ASSERT_GREATER(binNum, 1);
real_t binWidth = grid1 - grid0;
MRMD_HOST_ASSERT_GREATER(binNumForce, 1);

MultiView::HostMirror h_forcesRead("h_forcesRead", binNum, maxNumForces);
MultiView::HostMirror h_forcesRead("h_forcesRead", binNumForce, maxNumForces);

while (std::getline(fileThermoForce, line))
{
binNum = 0;
binNumForce = 0;
std::stringstream forceLineStream(line);
while (forceLineStream >> word)
{
h_forcesRead(binNum, histNum) = std::stod(word);
binNum++;
h_forcesRead(binNumForce, histNum) = std::stod(word);
binNumForce++;
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The variable name binNumForce is ambiguous. It's unclear whether this represents the number of force bins, a bin number for forces, or something else. The name was changed from binNum which was also not very descriptive. Consider using a more descriptive name like numBinsInForceFile or forceBinCount to clarify its purpose.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So far, my understanding was that numBins would signify the total number of bins and binNum the number of a specific bin. I was hoping to disambiguate a little here, but if this didn't improve things enough, I'd be happy to change it again.

Comment on lines +75 to 81
ThermodynamicForce(const std::vector<real_t>& targetDensities,
const data::Subdomain& subdomain,
const real_t& requestedDensityGridSpacing,
const real_t& requestedDensityBinWidth,
const std::vector<real_t>& thermodynamicForceModulation,
const std::vector<real_t>& thermodynamicForceModulations,
const bool enforceSymmetry = false,
const bool usePeriodicity = false);
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing documentation to distinguish requestedDensityGridSpacing from requestedDensityBinWidth. The ThermodynamicForce constructors use both parameters but their distinct purposes are not documented. Based on the implementation, gridSpacing appears to define the spacing between bin centers in the force grid, while binWidth defines the width of bins used during density profile sampling. This distinction should be documented in comments for the constructor parameters to avoid confusion.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, could be a naming inconsistency that slipped my attention.

Comment on lines +55 to +61
/**
* Interpolate data values of MultiHistogram to a new grid.
* @param input input MultiHistogram
* @param grid grid to interpolate to
* @return MultiHistogram with interpolated values
*/
data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid);
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing documentation for the input parameter's expected binSize. The documentation for the interpolate function doesn't specify whether the input MultiHistogram is expected to have uniform or non-uniform bin spacing. Additionally, the implementation assumes uniform spacing when calculating the interpolation factor on line 48 using input.inverseBinSize. This assumption should be documented.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point.

Comment on lines +21 to +81
data::MultiHistogram interpolate(const data::MultiHistogram& input, const ScalarView& grid)
{
data::MultiHistogram output(
"interpolated-profile", input.min, input.max, idx_c(grid.extent(0)), input.numHistograms);
const ScalarView& inputGrid = input.createGrid_d();

auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<2>>(
{idx_t(0), idx_t(0)}, {idx_c(grid.extent(0)), input.numHistograms});
auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx)
{
// find the two enclosing bins in the input histogram
auto rightBin = findRightBin(inputGrid, grid(binIdx));
auto leftBin = rightBin - 1;

// handle boundaries
if (leftBin < 0 || rightBin >= input.numBins)
{
output.data(binIdx, histogramIdx) = 0.0_r; // out of bounds, set to zero
return;
}

auto inputDataLeft = input.data(leftBin, histogramIdx);
auto inputDataRight = input.data(rightBin, histogramIdx);

output.data(binIdx, histogramIdx) =
lerp(inputDataLeft,
inputDataRight,
(grid(binIdx) - inputGrid(leftBin)) * input.inverseBinSize);
};
Kokkos::parallel_for("MultiHistogram::interpolate", policy, kernel);
Kokkos::fence();

return output;
}

data::MultiHistogram constrainToApplicationRegion(const data::MultiHistogram& input,
const util::ApplicationRegion& applicationRegion)
{
data::MultiHistogram constrainedProfile(
"constrained-profile", input.min, input.max, input.numBins, input.numHistograms);

auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<2>>({idx_t(0), idx_t(0)},
{input.numBins, input.numHistograms});
auto kernel = KOKKOS_LAMBDA(const idx_t binIdx, const idx_t histogramIdx)
{
if (!applicationRegion.isInApplicationRegion(input.getBinPosition(binIdx), 0_r, 0_r))
{
constrainedProfile.data(binIdx, histogramIdx) = 0_r;
}
else
{
constrainedProfile.data(binIdx, histogramIdx) = input.data(binIdx, histogramIdx);
}
};
Kokkos::parallel_for("MultiHistogram::constrainToApplicationRegion", policy, kernel);
Kokkos::fence();

return constrainedProfile;
}
} // namespace util
} // namespace mrmd No newline at end of file
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing test coverage for the new interpolation functionality. The new interpolate and constrainToApplicationRegion functions in mrmd/util/interpolation.cpp lack dedicated unit tests. Given that other utility functions in this directory have corresponding .test.cpp files (e.g., Random.test.cpp, math.test.cpp), consider adding interpolation.test.cpp to verify the correctness of the linear interpolation, boundary handling, and application region constraining logic.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, unit tests are necessary. I will implement them.

Comment on lines +42 to 54
auto policy = Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {numAtoms, numBins});
auto kernel = KOKKOS_LAMBDA(const idx_t atomIdx, const idx_t binIdx)
{
MRMD_DEVICE_ASSERT_GREATEREQUAL(types(idx), 0);
MRMD_DEVICE_ASSERT_LESS(types(idx), numTypes);
auto bin = histogram.getBin(positions(idx, axis));
if (bin == -1) return;
MRMD_DEVICE_ASSERT_GREATEREQUAL(types(atomIdx), 0);
MRMD_DEVICE_ASSERT_LESS(types(atomIdx), numTypes);
auto binPos = histogram.getBinPosition(binIdx);
auto atomPos = positions(atomIdx, axis);
auto access = scatter.access();
access(bin, types(idx)) += 1_r;
if (atomPos >= binPos - binWidth / 2 && atomPos <= binPos + binWidth / 2)
{
access(binIdx, types(atomIdx)) += 1_r;
}
};
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new implementation has significantly worse computational complexity. The algorithm changed from O(N) to O(N*M) where N is the number of atoms and M is the number of bins. The old implementation (commented out below) iterated only over atoms and directly computed which bin each atom belonged to. The new implementation uses a 2D iteration over all (atom, bin) pairs and checks if each atom falls within each bin's range. For large systems, this could cause significant performance degradation. Consider optimizing by iterating only over atoms and identifying which bins each atom should contribute to based on the binWidth parameter.

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a great point. We do indeed observe a significant loss in performance for this implementation. Should be changed accordingly.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What was the problem with the initial implementation? I do not get what your changes do differently? Has it to do with the out-of-bounds access you mentioned?

.def("sample", &action::ThermodynamicForce::sample)
.def("update", &action::ThermodynamicForce::update)
.def("update",
static_cast<void (action::ThermodynamicForce::*)(const real_t& smoothingSigma, const real_t& smoothingIntensity)>(&action::ThermodynamicForce::update))
Copy link

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The new update method overload with ApplicationRegion parameter is not exposed to Python. The C++ code adds a new overload of the update method that takes an ApplicationRegion parameter (line 66-68 in ThermodynamicForce.hpp), but this overload is not bound in the Python interface. If Python users need access to application region constrained updates, consider adding:

.def("update",
     static_cast<void (action::ThermodynamicForce::*)(const real_t&, const real_t&, const util::ApplicationRegion&)>(&action::ThermodynamicForce::update))
Suggested change
static_cast<void (action::ThermodynamicForce::*)(const real_t& smoothingSigma, const real_t& smoothingIntensity)>(&action::ThermodynamicForce::update))
static_cast<void (action::ThermodynamicForce::*)(const real_t& smoothingSigma, const real_t& smoothingIntensity)>(&action::ThermodynamicForce::update))
.def("update",
static_cast<void (action::ThermodynamicForce::*)(const real_t&, const real_t&, const util::ApplicationRegion&)>(&action::ThermodynamicForce::update))

Copilot uses AI. Check for mistakes.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will adapt the python interface accordingly.

@J-Hizzle
Copy link
Contributor Author

J-Hizzle commented Dec 4, 2025

Yes and no. I prefer having smaller PRs with clear scope. Let's start by putting the changes to ApplicationRegion and interpolation into their own PRs. They should be independent of the rest. We will then continue to split this PR into pieces.

Should I open a new PR for the changes in ApplicationRegion or can we somehow decouple the suggestions here into their own PR?

XzzX and others added 2 commits December 5, 2025 09:07
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
const real_t minArg,
const real_t maxArg,
idx_t numBinsArg,
real_t binSizeArg,
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The additional bin size argument is redundant and might be conflicting with min/max/num. Why do you need it?

const real_t min_ = min;
const real_t binSize_ = binSize;
auto policy = Kokkos::RangePolicy<>(0, numBins);
auto kernel = KOKKOS_LAMBDA(const idx_t idx)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not use getBinPosition?

Copy link
Contributor Author

@J-Hizzle J-Hizzle Dec 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did it to avoid capturing the this pointer therefore I opted for this redundant but inline implementation. How can I use getBinPosition here without capturing the this pointer?

{
grid[idx] = min_ + (real_c(idx) + 0.5_r) * binSize_;
};
Kokkos::parallel_for("MultiHistogram::scale", policy, kernel);
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wrong label

@J-Hizzle J-Hizzle mentioned this pull request Jan 20, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants