Merge pull request #7 from MakeNEnjoy/robin-test-linear-interpolation

Implemented bilinear interpolation

.gitignore

@@ -0,0 +1,3 @@
+.DS_Store
+.idea
+

linear-interpolation/.gitignore

@@ -0,0 +1,7 @@
+.DS_Store
+src/.DS_Store
+src/.cache
+src/build
+.idea
+src/cmake-build-debug
+src/cmake-build-release

linear-interpolation/README.md

@@ -0,0 +1,39 @@
+## Location of data
+The data path is hardcoded such that the following tree structure is assumed:
+```
+data/
+  grid.h5
+  hydrodynamic_U.h5
+  hydrodynamic_V.h5
+interactive-track-and-trace/
+  opening-hdf5/
+     ...
+```
+
+## Compiling
+Let the current directory be the `src` directory. Run:
+```shell
+mkdir build
+cd build
+cmake ..
+make
+```
+
+### Building with Linux
+Makes use of `mdspan` which is not supported by glibc++ at time of writing. See [compiler support](https://en.cppreference.com/w/cpp/compiler_support/23) for `mdspan`. The solution to this is to use Clang and libc++; this is configured in our CMake setup, however the default installation of the `netcdf-cxx` package on at least Arch linux (and suspectedly Debian derivatives as well) specifically builds for the glibc implementation. To get the netcdf C++ bindings functional with the libc++ implementation, one needs to build from source. On Linux, this requires a few changes to the CMake file included with the netcdf-cxx source code, which are detailed below.
+
+Step-by-step to build the program using clang++ and libc++ on linux:
+ 1. Download the source code of netcdf-cxx, found at 'https://github.com/Unidata/netcdf-cxx4/releases/tag/v4.3.1' (make sure to download the release source code, as the master branch contains non-compilable code).
+ 2. Edit the CMakeLists.txt file, by appending '-stdlib=libc++' to the `CMAKE_CXX_FLAGS` variable in line 430. This means line 430 should read: 
+    ```cmake 
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -Wall -Wno-unused-variable -Wno-unused-parameter -stdlib=libc++")
+    ```
+ 2. Build the source code with the following: 
+  ```sh 
+    mkdir build && cd build
+    cmake .. -DCMAKE_CXX_COMPILER=/usr/bin/clang++
+    make
+    ctest
+    sudo make install
+  ```
+ 3. Now the code should compile through the standard steps described in the Compiling section.

linear-interpolation/src/CMakeLists.txt

@@ -0,0 +1,36 @@
+cmake_minimum_required (VERSION 3.28)
+project (LinearInterpolate)
+
+set(CMAKE_CXX_STANDARD 23)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+find_package(netCDF REQUIRED)
+
+add_executable(LinearInterpolate main.cpp
+        readdata.cpp
+        readdata.h
+        interpolate.cpp
+        interpolate.h
+        UVGrid.cpp
+        UVGrid.h
+        Vel.h
+        Vel.cpp)
+
+execute_process(
+    COMMAND nc-config --includedir
+    OUTPUT_VARIABLE NETCDF_INCLUDE_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+execute_process(
+    COMMAND ncxx4-config --libdir
+    OUTPUT_VARIABLE NETCDFCXX_LIB_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+target_include_directories(LinearInterpolate PUBLIC ${netCDF_INCLUDE_DIR})
+
+find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH)
+target_link_libraries(LinearInterpolate ${NETCDF_LIB})

linear-interpolation/src/UVGrid.cpp

@@ -0,0 +1,61 @@
+#include <ranges>
+
+#include "UVGrid.h"
+#include "readdata.h"
+
+#define sizeError2 "The sizes of the hydrodynamic data files are different"
+#define sizeError "The sizes of the hydrodynamicU or -V files does not correspond with the sizes of the grid file"
+
+using namespace std;
+
+UVGrid::UVGrid() {
+    auto us = readHydrodynamicU();
+    auto vs = readHydrodynamicV();
+    if (us.size() != vs.size()) {
+        throw domain_error(sizeError2);
+    }
+
+    tie(times, lats, lons) = readGrid();
+
+    timeSize = times.size();
+    latSize = lats.size();
+    lonSize = lons.size();
+
+    size_t gridSize = timeSize * latSize * lonSize;
+    if (gridSize != us.size()) {
+        throw domain_error(sizeError);
+    }
+
+    uvData.reserve(gridSize);
+
+    for (auto vel: views::zip(us, vs)) {
+        uvData.push_back(Vel(vel));
+    }
+}
+
+const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const {
+    size_t index = timeIndex * (latSize * lonSize) + latIndex * lonIndex + lonIndex;
+    return uvData[index];
+}
+
+double UVGrid::lonStep() const {
+    return lons[1] - lons[0];
+}
+
+double UVGrid::latStep() const {
+    return lats[1] - lats[0];
+}
+
+int UVGrid::timeStep() const {
+    return times[1] - times[0];
+}
+
+void UVGrid::streamSlice(ostream &os, size_t t) {
+    for (int x = 0; x < latSize; x++) {
+        for (int y = 0; y < lonSize; y++) {
+            auto vel = (*this)[t,x,y];
+            os << vel << " ";
+        }
+      os << endl;
+    }
+}

linear-interpolation/src/UVGrid.h

@@ -0,0 +1,51 @@
+#ifndef LINEARINTERPOLATE_UVGRID_H
+#define LINEARINTERPOLATE_UVGRID_H
+
+#include <vector>
+#include "Vel.h"
+
+class UVGrid {
+private:
+    /**
+     * 1D data vector of all the us and vs
+     */
+    std::vector<Vel> uvData;
+public:
+    UVGrid();
+
+    size_t timeSize;
+    size_t latSize;
+    size_t lonSize;
+
+    /**
+     * Assuming grid is a regular grid, gives the longitudinal spacing of grid.
+     * @return longitudinal spacing
+     */
+    double lonStep() const;
+
+    /**
+     * Assuming grid is a regular grid, gives the latitudinal spacing of grid.
+     * @return latitudinal spacing
+     */
+    double latStep() const;
+
+    /**
+     * Assuming grid is a regular grid, gives the time spacing of grid.
+     * @return time spacing
+     */
+    int timeStep() const;
+
+    std::vector<int> times;
+    std::vector<double> lats;
+    std::vector<double> lons;
+
+    /**
+     * The 3D index into the data
+     * @return Velocity at that index
+     */
+    const Vel& operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const;
+
+    void streamSlice(std::ostream &os, size_t t);
+};
+
+#endif //LINEARINTERPOLATE_UVGRID_H

linear-interpolation/src/Vel.cpp

@@ -0,0 +1,40 @@
+#include "Vel.h"
+#include <stdexcept>
+#include <iomanip>
+
+using namespace std;
+
+Vel::Vel(double u, double v) : u(u), v(v) {}
+
+Vel::Vel(const std::pair<double, double>& p) : u(p.first), v(p.second) {}
+
+Vel& Vel::operator=(const std::pair<double, double>& p) {
+    u = p.first;
+    v = p.second;
+    return *this;
+}
+
+Vel Vel::operator+(const Vel& other) const {
+    return Vel(u + other.u, v + other.v);
+}
+
+Vel& Vel::operator+=(const Vel& other) {
+    u += other.u;
+    v += other.v;
+    return *this;
+}
+
+template<typename Scalar>
+Vel Vel::operator/(Scalar scalar) const {
+    if (scalar == 0) throw std::runtime_error("Division by zero");
+    return Vel(u / scalar, v / scalar);
+}
+
+std::ostream& operator<<(ostream& os, const Vel& vel) {
+    os << "(";
+    os << fixed << setprecision(2) << setw(5) << vel.u;
+    os << ", ";
+    os << fixed << setprecision(2) << setw(5) << vel.v;
+    os << ")";
+    return os;
+}

linear-interpolation/src/Vel.h

@@ -0,0 +1,44 @@
+#ifndef LINEARINTERPOLATE_VEL_H
+#define LINEARINTERPOLATE_VEL_H
+
+#include <utility>
+#include <stdexcept>
+#include <iostream>
+#include <format>
+
+class Vel {
+public:
+    double u; // Eastward Current Velocity in the Water Column
+    double v; // Northward Current Velocity in the Water Column
+
+    Vel(double u, double v);
+    Vel(const std::pair<double, double>& p);  // Conversion constructor
+    Vel& operator=(const std::pair<double, double>& p);
+
+    // Operator + to add two Vel objects
+    Vel operator+(const Vel& other) const;
+
+    // Operator += to add another Vel object to this object
+    Vel& operator+=(const Vel& other);
+
+    // Operator * to multiply Vel by a scalar, defined as a member template
+    template<typename Scalar>
+    Vel operator*(Scalar scalar) const {
+        return Vel(u * scalar, v * scalar);
+    }
+
+    // Operator / to divide Vel by a scalar, defined as a member template
+    template<typename Scalar>
+    Vel operator/(Scalar scalar) const;
+
+    // Friend declaration for the stream insertion operator
+    friend std::ostream& operator<<(std::ostream& os, const Vel& vel);
+};
+
+// Non-member function for scalar multiplication on the left
+template<typename Scalar>
+Vel operator*(Scalar scalar, const Vel& p) {
+    return Vel(p.u * scalar, p.v * scalar);
+}
+
+#endif //LINEARINTERPOLATE_VEL_H

linear-interpolation/src/interpolate.cpp

@@ -0,0 +1,47 @@
+#include "interpolate.h"
+
+using namespace std;
+
+Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon) {
+    double latStep = uvGrid.latStep();
+    double lonStep = uvGrid.lonStep();
+    int timeStep = uvGrid.timeStep();
+
+    int latIndex = (lat - uvGrid.lats[0]) / latStep;
+    int lonIndex = (lon - uvGrid.lons[0]) / lonStep;
+    int timeIndex = (time - uvGrid.times[0]) / timeStep;
+
+    double timeRatio = (static_cast<double>(time) - uvGrid.times[timeIndex]) / timeStep;
+    double latRatio = (lat - uvGrid.lats[latIndex]) / latStep;
+    double lonRatio = (lon - uvGrid.lons[lonIndex]) / lonStep;
+
+    Vel point = {0, 0};
+    for (int timeOffset = 0; timeOffset <= 1; timeOffset++) {
+        for (int latOffset = 0; latOffset <= 1; latOffset++) {
+            for (int lonOffset = 0; lonOffset <= 1; lonOffset++) {
+                auto vertex = uvGrid[
+                        timeIndex + 1 < uvGrid.timeSize ? timeIndex + timeOffset : timeIndex,
+                                latIndex + 1 < uvGrid.latSize ? latIndex + latOffset : latIndex,
+                                lonIndex + 1 < uvGrid.lonSize ? lonIndex + lonOffset : lonIndex
+                ];
+
+                double timeRation = (1 - timeOffset) * (1 - timeRatio) + timeOffset * timeRatio;
+                double latRation = (1 - latOffset) * (1 - latRatio) + latOffset * latRatio;
+                double lonRation = (1 - lonOffset) * (1 - lonRatio) + lonOffset * lonRatio;
+                point += timeRation * latRation * lonRation * vertex;
+            }
+        }
+    }
+
+    return point;
+}
+
+vector<Vel> bilinearInterpolate(const UVGrid &uvGrid, vector<tuple<int, double, double>> points) {
+    vector<Vel> result;
+    result.reserve(points.size());
+    for (auto [time, lat, lon]: points) {
+        result.push_back(bilinearInterpolate(uvGrid, time, lat, lon));
+    }
+
+    return result;
+}

linear-interpolation/src/interpolate.h

@@ -0,0 +1,28 @@
+#ifndef LINEARINTERPOLATE_INTERPOLATE_H
+#define LINEARINTERPOLATE_INTERPOLATE_H
+
+#include <vector>
+
+#include "UVGrid.h"
+
+/**
+ * Bilinearly interpolate the point (time, lat, lon) to produce the interpolated velocity.
+ * Since it is in 3D, this means that it interpolates against 8 points (excluding edges).
+ * As described in https://numerical.recipes/book.html Chapter 3.6
+ * @param uvGrid velocity grid
+ * @param time time of point
+ * @param lat latitude of point
+ * @param lon longitude of point
+ * @return interpolated velocity
+ */
+Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon);
+
+/**
+ * Helper function for bilnearly interpolating a vector of points
+ * @param uvGrid velocity grid
+ * @param points vector of points
+ * @return interpolated velocities
+ */
+std::vector<Vel> bilinearInterpolate(const UVGrid &uvGrid, std::vector<std::tuple<int, double, double>> points);
+
+#endif //LINEARINTERPOLATE_INTERPOLATE_H

linear-interpolation/src/main.cpp

@@ -0,0 +1,51 @@
+#include "interpolate.h"
+#include "Vel.h"
+#include <ranges>
+#include <chrono>
+
+using namespace std;
+
+int main() {
+    UVGrid uvGrid;
+    uvGrid.streamSlice(cout, 100);
+
+    int N = 10000000;  // Number of points
+
+    int time_start = 0;
+    int time_end = 31536000;
+
+    double lat_start = 46.125;
+    double lat_end = 62.625;
+
+    double lon_start = -15.875;
+    double lon_end = 12.875;
+
+    double lat_step = (lat_end - lat_start) / (N - 1);
+    double lon_step = (lon_end - lon_start) / (N - 1);
+    int time_step = (time_end - time_start) / (N - 1);
+
+    vector<tuple<int, double, double>> points;
+
+    for(int i = 0; i < N; i++) {
+        points.push_back({time_start+time_step*i, lat_start+lat_step*i, lon_start+lon_step*i});
+    }
+
+    auto start = chrono::high_resolution_clock::now();
+
+    auto x = bilinearInterpolate(uvGrid, points);
+
+    auto stop = chrono::high_resolution_clock::now();
+
+    auto duration = chrono::duration_cast<std::chrono::milliseconds >(stop - start);
+
+    cout << "Time taken for " << N << " points was " << duration.count()/1000. << " seconds\n";
+
+    // Do something with result in case of optimisation
+    double sum = 0;
+    for (auto [u,v]: x) {
+        sum += u + v;
+    }
+    cout << "Sum = " << sum << endl;
+
+    return 0;
+}

linear-interpolation/src/readdata.cpp

@@ -0,0 +1,48 @@
+#include <stdexcept>
+
+#include <netcdf>
+
+#include "readdata.h"
+
+using namespace std;
+using namespace netCDF;
+
+template <typename T>
+vector<T> getVarVector(const NcVar &var) {
+    int length = 1;
+    for (NcDim dim : var.getDims()) {
+        length *= dim.getSize();
+    }
+
+    vector<T> vec(length);
+
+    var.getVar(vec.data());
+
+    return vec;
+}
+
+vector<double> readHydrodynamicU() {
+    netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return getVarVector<double>(vars.find("uo")->second);
+}
+
+vector<double> readHydrodynamicV() {
+    netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return getVarVector<double>(vars.find("vo")->second);
+}
+
+tuple<vector<int>, vector<double>, vector<double>> readGrid() {
+    netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
+    multimap< string, NcVar > vars = data.getVars();
+    vector<int> time = getVarVector<int>(vars.find("times")->second);
+    vector<double> longitude = getVarVector<double>(vars.find("longitude")->second);
+    vector<double> latitude = getVarVector<double>(vars.find("latitude")->second);
+
+    return {time, latitude, longitude};
+}

linear-interpolation/src/readdata.h

@@ -0,0 +1,22 @@
+#ifndef LINEARINTERPOLATE_READDATA_H
+#define LINEARINTERPOLATE_READDATA_H
+
+/**
+ * reads the file hydrodynamic_U.h5
+ * @return the data vector of us
+ */
+std::vector<double> readHydrodynamicU();
+
+/**
+ * reads the file hydrodynamic_V.h5
+ * @return the data vector of vs
+ */
+std::vector<double> readHydrodynamicV();
+
+/**
+ * Reads the file grid.h5
+ * @return a tuple of (times, latitude, longitude)
+ */
+std::tuple<std::vector<int>, std::vector<double>, std::vector<double>> readGrid();
+
+#endif //LINEARINTERPOLATE_READDATA_H