removed mdspan requirement

.gitignore

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

linear-interpolation/.gitignore

@@ -3,4 +3,5 @@ src/.DS_Store
 src/.cache
 src/build
 .idea
-src/cmake-build-debug
+src/cmake-build-debug
+src/cmake-build-release

linear-interpolation/src/CMakeLists.txt

@@ -4,10 +4,6 @@ project (LinearInterpolate)
 set(CMAKE_CXX_STANDARD 23)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
-# Force use of libc++ for mdspan support
-set(CMAKE_CXX_COMPILER "clang++")
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++ ")
-
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 find_package(netCDF REQUIRED)
@@ -15,14 +11,12 @@ find_package(netCDF REQUIRED)
 add_executable(LinearInterpolate main.cpp
         readdata.cpp
         readdata.h
-        vecutils.cpp
-        vecutils.h
         interpolate.cpp
         interpolate.h
         UVGrid.cpp
         UVGrid.h
-        point.h
+        Vel.h
-        point.cpp)
+        Vel.cpp)
 
 execute_process(
     COMMAND nc-config --includedir

linear-interpolation/src/UVGrid.cpp

@@ -1,18 +1,42 @@
-#include <mdspan>
 #include <ranges>
+#include <print>
 
 #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, sizeU] = readHydrodynamicU();
+    auto us = readHydrodynamicU();
-    auto [vs, sizeV] = readHydrodynamicV();
+    auto vs = readHydrodynamicV();
-    // Assuming sizeU == sizeV
+    if (us.size() != vs.size()) {
-    uvData = views::zip(us, vs) | ranges::to<vector>();
+        throw domain_error(sizeError2);
-    uvMatrix = mdspan(uvData.data(), sizeU);
+    }
+
     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 = views::zip(us, vs)
+             | views::transform([](auto pair) {
+        return Vel(pair);
+    })
+             | ranges::to<vector>();
+}
+
+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 {
@@ -20,9 +44,19 @@ double UVGrid::lonStep() const {
 }
 
 double UVGrid::latStep() const {
-    return lats[1]-lats[0];
+    return lats[1] - lats[0];
 }
 
 int UVGrid::timeStep() const {
-    return times[1]-times[0];
+    return times[1] - times[0];
-}
+}
+
+void UVGrid::printSlice(size_t t) {
+    for (int x = 0; x < latSize; x++) {
+        for (int y = 0; y < lonSize; y++) {
+            auto [u,v] = (*this)[t,x,y];
+            print("({:>7.4f}, {:>7.4f}) ", u, v);
+        }
+        println("");
+    }
+}

linear-interpolation/src/UVGrid.h

@@ -2,29 +2,54 @@
 #define LINEARINTERPOLATE_UVGRID_H
 
 #include <vector>
-#include "vecutils.h"
+#include "Vel.h"
-#include "point.h"
 
 class UVGrid {
 private:
     /**
-     * u == Eastward Current Velocity in the Water Column
+     * 1D data vector of all the us and vs
-     * v == Northward Current Velocity in the Water Column
      */
-    std::vector<point> uvData;
+    std::vector<Vel> uvData;
 public:
     UVGrid();
 
-    // The step functions assume regular spacing
+    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;
 
-    arr3d<point> uvMatrix;
+    /**
+     * The 3D index into the data
+     * @return Velocity at that index
+     */
+    const Vel& operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const;
+
+    // Friend declaration for the stream insertion operator
+    friend std::ostream& operator<<(std::ostream& os, const UVGrid& vel);
+
+    void printSlice(size_t t);
 };
 
+
 #endif //LINEARINTERPOLATE_UVGRID_H

linear-interpolation/src/Vel.cpp

@@ -0,0 +1,33 @@
+#include "Vel.h"
+#include <stdexcept>
+
+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<<(std::ostream& os, const Vel& vel) {
+    os << "(" << vel.u << ", " << vel.v << ")";
+    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

@@ -1,39 +1,37 @@
 #include <cmath>
-
 #include <ranges>
+#include <print>
+
 #include "interpolate.h"
 
 using namespace std;
 
-// Inspired by https://numerical.recipes/book.html Chapter 3.6
+Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon) {
-point bilinearInterpolate(const UVGrid &uvGrid, std::tuple<int, double, double> timeLatLon) {
-    auto [time, lat, lon] = timeLatLon;
-
     double latStep = uvGrid.latStep();
     double lonStep = uvGrid.lonStep();
     int timeStep = uvGrid.timeStep();
 
-    int latIndex = (lat-uvGrid.lats[0])/latStep;
+    int latIndex = (lat - uvGrid.lats[0]) / latStep;
-    int lonIndex = (lon-uvGrid.lons[0])/lonStep;
+    int lonIndex = (lon - uvGrid.lons[0]) / lonStep;
-    int timeIndex = (time-uvGrid.times[0])/timeStep;
+    int timeIndex = (time - uvGrid.times[0]) / timeStep;
 
-    double timeRatio = (static_cast<double>(time)-uvGrid.times[timeIndex])/timeStep;
+    double timeRatio = (static_cast<double>(time) - uvGrid.times[timeIndex]) / timeStep;
-    double latRatio = (lat-uvGrid.lats[latIndex]) / latStep;
+    double latRatio = (lat - uvGrid.lats[latIndex]) / latStep;
-    double lonRatio = (lon-uvGrid.lons[lonIndex]) / lonStep;
+    double lonRatio = (lon - uvGrid.lons[lonIndex]) / lonStep;
 
-    point point = {0, 0};
+    Vel point = {0, 0};
-    for(int time = 0; time <= 1; time++) {
+    for (int timeOffset = 0; timeOffset <= 1; timeOffset++) {
-        for(int lat = 0; lat <= 1; lat++) {
+        for (int latOffset = 0; latOffset <= 1; latOffset++) {
-            for(int lon = 0; lon <= 1; lon++) {
+            for (int lonOffset = 0; lonOffset <= 1; lonOffset++) {
-                auto vertex = uvGrid.uvMatrix[
+                auto vertex = uvGrid[
-                        timeIndex + time < uvGrid.uvMatrix.extent(0) ? timeIndex + 1 : timeIndex,
+                        timeIndex + 1 < uvGrid.timeSize ? timeIndex + timeOffset : timeIndex,
-                        latIndex + lat < uvGrid.uvMatrix.extent(1) ? latIndex + 1 : latIndex,
+                                latIndex + 1 < uvGrid.latSize ? latIndex + latOffset : latIndex,
-                        lonIndex + lon < uvGrid.uvMatrix.extent(2) ? lonIndex + 1 : lonIndex
+                                lonIndex + 1 < uvGrid.lonSize ? lonIndex + lonOffset : lonIndex
                 ];
 
-                double timeRation = (1 - time)*(1-timeRatio) + time*timeRatio;
+                double timeRation = (1 - timeOffset) * (1 - timeRatio) + timeOffset * timeRatio;
-                double latRation = (1 - lat)*(1-latRatio) + lat*latRatio;
+                double latRation = (1 - latOffset) * (1 - latRatio) + latOffset * latRatio;
-                double lonRation = (1 - lon)*(1-lonRatio) + lon*lonRatio;
+                double lonRation = (1 - lonOffset) * (1 - lonRatio) + lonOffset * lonRatio;
                 point += timeRation * latRation * lonRation * vertex;
             }
         }
@@ -41,3 +39,13 @@ point bilinearInterpolate(const UVGrid &uvGrid, std::tuple<int, double, double>
 
     return point;
 }
+
+vector<Vel> bilinearInterpolate(const UVGrid &uvGrid, vector<tuple<int, double, double>> points) {
+    auto results = points
+                   | std::views::transform([&uvGrid](const auto &point) {
+        auto [time, lat, lon] = point;
+        return bilinearInterpolate(uvGrid, time, lat, lon);
+    })
+                   | std::ranges::to<std::vector<Vel>>();
+    return results;
+}

linear-interpolation/src/interpolate.h

@@ -5,8 +5,24 @@
 
 #include "UVGrid.h"
 
-point bilinearInterpolate(const UVGrid &uvGrid, std::tuple<int, double, double> timeLatLong);
+/**
+ * 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);
 
-std::vector<double> bilinearInterpolate();
+/**
+ * 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

@@ -1,13 +1,52 @@
 #include "interpolate.h"
+#include "Vel.h"
+#include <print>
+#include <ranges>
 
 using namespace std;
 
 int main() {
     UVGrid uvGrid;
-    auto p = bilinearInterpolate(uvGrid, {392400, 53, -14.5});
+    uvGrid.printSlice(100);
 
-    println("({}, {})", p.first, p.second);
+    int N = 10000000;  // Number of points
-    p = bilinearInterpolate(uvGrid, {802400, 62, -14.5});
+
+    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;
+
+    // Calculate increments
+    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 = std::chrono::high_resolution_clock::now();
+
+    auto x = bilinearInterpolate(uvGrid, points);
+
+    auto stop = std::chrono::high_resolution_clock::now();
+
+    auto duration = std::chrono::duration_cast<std::chrono::milliseconds >(stop - start);
+
+    println("Time taken for {} points was {} seconds", N, duration.count()/1000.);
+
+    // Do something with result in case of optimisation
+    double sum = 0;
+    for (auto [u,v]: x) {
+        sum += u + v;
+    }
+    println("Sum = {}", sum);
 
     return 0;
 }

linear-interpolation/src/point.cpp

@@ -1,11 +0,0 @@
-#include "point.h"
-
-point operator+(const point& p1, const point& p2) {
-    return {p1.first + p2.first, p1.second + p2.second};
-}
-
-point& operator+=(point& p1, const point& p2) {
-    p1.first += p2.first;
-    p1.second += p2.second;
-    return p1;
-}

linear-interpolation/src/point.h

@@ -1,27 +0,0 @@
-#ifndef LINEARINTERPOLATE_POINT_H
-#define LINEARINTERPOLATE_POINT_H
-
-#include <utility>
-
-using point = std::pair<double, double>; // {u, v}
-
-point operator+(const point& p1, const point& p2);
-
-point& operator+=(point& p1, const point& p2);
-
-template<typename Scalar>
-point operator*(const point& p, Scalar scalar) {
-    return {p.first * scalar, p.second * scalar};
-}
-
-template<typename Scalar>
-point operator*(Scalar scalar, const point& p) {
-    return {p.first * scalar, p.second * scalar};
-}
-
-template<typename Scalar>
-point operator/(const point& p, Scalar scalar) {
-    return {p.first / scalar, p.second / scalar};
-}
-
-#endif //LINEARINTERPOLATE_POINT_H

linear-interpolation/src/readdata.cpp

@@ -22,45 +22,20 @@ vector<T> getVarVector(const NcVar &var) {
     return vec;
 }
 
-/**
+vector<double> readHydrodynamicU() {
- * Read a 3D matrix from a NetCDF variable.
- * Reads data into a contiguous 1D data vector.
- * Returns a pair of the size of the matrix (in the form of an extent) with the data vector.
- *
- * Inteded usage of this function involves using the two returned values
- * to create an mdspan:
- *
- * auto arr = mdspan(vec.data(), size);
- */
-template <typename T>
-pair<vector<T>, std::dextents<std::size_t, 3>> get3DMat(const NcVar &var) {
-    if(var.getDimCount() != 3) {
-        throw invalid_argument("Variable is not 3D");
-    }
-    int timeLength = var.getDim(0).getSize();
-    int latLength = var.getDim(1).getSize();
-    int longLength = var.getDim(2).getSize();
-    vector<T> vec(timeLength*latLength*longLength);
-    var.getVar(vec.data());
-    auto arr = std::mdspan(vec.data(), timeLength, latLength, longLength);
-
-    return {vec, arr.extents()};
-}
-
-pair<vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicU() {
     netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read);
 
     multimap< string, NcVar > vars = data.getVars();
 
-    return get3DMat<double>(vars.find("uo")->second);
+    return getVarVector<double>(vars.find("uo")->second);
 }
 
-pair<vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicV() {
+vector<double> readHydrodynamicV() {
     netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read);
 
     multimap< string, NcVar > vars = data.getVars();
 
-    return get3DMat<double>(vars.find("vo")->second);
+    return getVarVector<double>(vars.find("vo")->second);
 }
 
 tuple<vector<int>, vector<double>, vector<double>> readGrid() {

linear-interpolation/src/readdata.h

@@ -1,19 +1,17 @@
 #ifndef LINEARINTERPOLATE_READDATA_H
 #define LINEARINTERPOLATE_READDATA_H
 
-#include "vecutils.h"
-
 /**
  * reads the file hydrodynamic_U.h5
- * @return a pair of the data vector of the contents and its dimensions to be used with mdspan
+ * @return the data vector of us
  */
-std::pair<std::vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicU();
+std::vector<double> readHydrodynamicU();
 
 /**
  * reads the file hydrodynamic_V.h5
- * @return a pair of the data vector of the contents and its dimensions to be used with mdspan
+ * @return the data vector of vs
  */
-std::pair<std::vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicV();
+std::vector<double> readHydrodynamicV();
 
 /**
  * Reads the file grid.h5

linear-interpolation/src/vecutils.cpp

@@ -1,24 +0,0 @@
-#include <print>
-
-#include "vecutils.h"
-
-using namespace std;
-
-void print3DMatrixSlice(const arr3d<double> &arr, int t) {
-    for (int x = 0; x < arr.extent(1); x++) {
-        for (int y = 0; y < arr.extent(2); y++) {
-            print("{:>10.4f} ", arr[t,x,y]);
-        }
-        println("");
-    }
-}
-
-void print3DMatrixSlice(const arr3d<std::pair<double, double>> &arr, int t) {
-    for (int x = 0; x < arr.extent(1); x++) {
-        for (int y = 0; y < arr.extent(2); y++) {
-            auto [u,v] = arr[t,x,y];
-            print("({:>7.4f}, {:>7.4f}) ", u, v);
-        }
-        println("");
-    }
-}

linear-interpolation/src/vecutils.h

@@ -1,48 +0,0 @@
-#ifndef LINEARINTERPOLATE_VECUTILS_H
-#define LINEARINTERPOLATE_VECUTILS_H
-
-#include <mdspan>
-#include <print>
-
-template <typename T>
-using arr3d = std::mdspan<
-        T,
-        std::dextents<
-                std::size_t,
-                3
-        >
->;
-
-/**
- * Print contents of vector
- * @tparam T The type of the data inside the vector
- * @param vec The vector to be printed
- */
-template <typename T>
-void printContentsOfVec(const std::vector<T>& vec) {
-    for (const auto& element : vec) {
-        std::print("{} ", element);
-
-    }
-    std::println("");
-}
-
-/**
- * Print matrix [x,y] for all values arr[t,x,y]
- * @param arr matrix to be printed
- * @param t value to slice matrix with
- */
-template <typename T>
-void print3DMatrixSlice(const arr3d<T> &arr, int t) {
-    for (int x = 0; x < arr.extent(1); x++) {
-        for (int y = 0; y < arr.extent(2); y++) {
-            std::print("{} ", arr[t,x,y]);
-        }
-        std::println("");
-    }
-}
-
-void print3DMatrixSlice(const arr3d<double> &arr, int t);
-void print3DMatrixSlice(const arr3d<std::pair<double, double>> &arr, int t);
-
-#endif //LINEARINTERPOLATE_VECUTILS_H