From 2443d96733c671906450617e82529fcd6d72b510 Mon Sep 17 00:00:00 2001 From: robin Date: Thu, 18 Apr 2024 13:21:53 +0200 Subject: [PATCH 01/15] feat: implemented reading hdf5 --- .gitignore | 2 + opening-hdf5/.gitignore | 4 ++ opening-hdf5/src/CMakeLists.txt | 28 ++++++++++ opening-hdf5/src/main.cpp | 99 +++++++++++++++++++++++++++++++++ 4 files changed, 133 insertions(+) create mode 100644 .gitignore create mode 100644 opening-hdf5/.gitignore create mode 100644 opening-hdf5/src/CMakeLists.txt create mode 100644 opening-hdf5/src/main.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5ca0973 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +.DS_Store + diff --git a/opening-hdf5/.gitignore b/opening-hdf5/.gitignore new file mode 100644 index 0000000..13c6291 --- /dev/null +++ b/opening-hdf5/.gitignore @@ -0,0 +1,4 @@ +.DS_Store +src/.DS_Store +src/.cache +src/build diff --git a/opening-hdf5/src/CMakeLists.txt b/opening-hdf5/src/CMakeLists.txt new file mode 100644 index 0000000..07d1489 --- /dev/null +++ b/opening-hdf5/src/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required (VERSION 3.28) +project (ReadHDF5) + +set(CMAKE_CXX_STANDARD 23) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +find_package(netCDF REQUIRED) + +add_executable(ReadHDF5 main.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(ReadHDF5 PUBLIC ${netCDF_INCLUDE_DIR}) + +find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH) +target_link_libraries(ReadHDF5 ${NETCDF_LIB}) diff --git a/opening-hdf5/src/main.cpp b/opening-hdf5/src/main.cpp new file mode 100644 index 0000000..1acc55d --- /dev/null +++ b/opening-hdf5/src/main.cpp @@ -0,0 +1,99 @@ +#include +#include +#include +#include + +#include + +using namespace std; +using namespace netCDF; + +template +void printContentsOfVec(const std::vector& vec) { + for (const auto& element : vec) { + cout << element << " "; + } + cout << std::endl; +} + +template +vector getVarVector(NcVar var) { + int length = 1; + for (NcDim dim : var.getDims()) { + length *= dim.getSize(); + } + + vector vec(length); + + var.getVar(vec.data()); + + return vec; +} + +template +using arr3d = std::mdspan< + T, + std::dextents< + std::size_t, + 3 + > + >; + + +template +void print3DMatrixSlice(arr3d arr, int t) { + for (int x = 0; x < arr.extent(1); x++) { + for (int y = 0; y < arr.extent(2); y++) { + print("{} ", arr[t,x,y]); + } + println(""); + } +} + +void print3DMatrixSlice(arr3d 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(""); + } +} + +/** + * 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 +pair, std::dextents> get3DMat(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 vec(timeLength*latLength*longLength); + var.getVar(vec.data()); + auto arr = std::mdspan(vec.data(), timeLength, latLength, longLength); + + return {vec, arr.extents()}; +} + +int main () { + netCDF::NcFile data("../../../hdf5/hydrodynamic_U.h5", netCDF::NcFile::read); + + multimap< string, NcVar > vars = data.getVars(); + + auto [vec, size] = get3DMat(vars.find("uo")->second); + + auto arr = std::mdspan(vec.data(), size); + + print3DMatrixSlice(arr, 100); + + return 0; +} From 4c2dfa8e2b9fac0fc84143eda22710ba2d9b59f2 Mon Sep 17 00:00:00 2001 From: robin Date: Thu, 18 Apr 2024 13:42:31 +0200 Subject: [PATCH 02/15] readme and path --- opening-hdf5/README.md | 21 +++++++++++++++++++++ opening-hdf5/src/main.cpp | 2 +- 2 files changed, 22 insertions(+), 1 deletion(-) create mode 100644 opening-hdf5/README.md diff --git a/opening-hdf5/README.md b/opening-hdf5/README.md new file mode 100644 index 0000000..301bc25 --- /dev/null +++ b/opening-hdf5/README.md @@ -0,0 +1,21 @@ +## 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 +Makes use of `mdspan` which is not supported by GCC libstdc++ at time of writing. See [compiler support](https://en.cppreference.com/w/cpp/compiler_support/23) for `mdspan`. +Let the current directory be the `src` directory. Run: +```shell +mkdir build +cd build +cmake .. +make +``` diff --git a/opening-hdf5/src/main.cpp b/opening-hdf5/src/main.cpp index 1acc55d..691e160 100644 --- a/opening-hdf5/src/main.cpp +++ b/opening-hdf5/src/main.cpp @@ -85,7 +85,7 @@ pair, std::dextents> get3DMat(NcVar var) { } int main () { - netCDF::NcFile data("../../../hdf5/hydrodynamic_U.h5", netCDF::NcFile::read); + netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); multimap< string, NcVar > vars = data.getVars(); From ce91a9384d0682af20e22a4dc18603cd03db573f Mon Sep 17 00:00:00 2001 From: robin Date: Thu, 18 Apr 2024 14:44:48 +0200 Subject: [PATCH 03/15] readme update --- opening-hdf5/README.md | 5 ++++- opening-hdf5/src/CMakeLists.txt | 2 +- opening-hdf5/src/main.cpp | 1 + 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/opening-hdf5/README.md b/opening-hdf5/README.md index 301bc25..44aca1b 100644 --- a/opening-hdf5/README.md +++ b/opening-hdf5/README.md @@ -11,7 +11,6 @@ interactive-track-and-trace/ ``` ## Compiling -Makes use of `mdspan` which is not supported by GCC libstdc++ at time of writing. See [compiler support](https://en.cppreference.com/w/cpp/compiler_support/23) for `mdspan`. Let the current directory be the `src` directory. Run: ```shell mkdir build @@ -19,3 +18,7 @@ cd build cmake .. make ``` + +### Building with Linux +Makes use of `mdspan` which is not supported by GCC libstdc++ at time of writing. See [compiler support](https://en.cppreference.com/w/cpp/compiler_support/23) for `mdspan`. Probably you will need to install Clang libc++. + diff --git a/opening-hdf5/src/CMakeLists.txt b/opening-hdf5/src/CMakeLists.txt index 07d1489..aad0646 100644 --- a/opening-hdf5/src/CMakeLists.txt +++ b/opening-hdf5/src/CMakeLists.txt @@ -6,7 +6,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -find_package(netCDF REQUIRED) +# find_package(netCDF REQUIRED) add_executable(ReadHDF5 main.cpp) diff --git a/opening-hdf5/src/main.cpp b/opening-hdf5/src/main.cpp index 691e160..68e0ee8 100644 --- a/opening-hdf5/src/main.cpp +++ b/opening-hdf5/src/main.cpp @@ -91,6 +91,7 @@ int main () { auto [vec, size] = get3DMat(vars.find("uo")->second); + auto arr = std::mdspan(vec.data(), size); print3DMatrixSlice(arr, 100); From f492ebd91e42157115566471e19b11b544a49f15 Mon Sep 17 00:00:00 2001 From: robin Date: Thu, 18 Apr 2024 14:45:37 +0200 Subject: [PATCH 04/15] fix: cmakelists.txt --- opening-hdf5/src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opening-hdf5/src/CMakeLists.txt b/opening-hdf5/src/CMakeLists.txt index aad0646..07d1489 100644 --- a/opening-hdf5/src/CMakeLists.txt +++ b/opening-hdf5/src/CMakeLists.txt @@ -6,7 +6,7 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -# find_package(netCDF REQUIRED) +find_package(netCDF REQUIRED) add_executable(ReadHDF5 main.cpp) From ffdfd0a40be2bb81c5368ad4250c70b21c16cc36 Mon Sep 17 00:00:00 2001 From: djairoh Date: Thu, 18 Apr 2024 15:18:08 +0200 Subject: [PATCH 05/15] fix: compiling for linux --- opening-hdf5/src/CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opening-hdf5/src/CMakeLists.txt b/opening-hdf5/src/CMakeLists.txt index 07d1489..e1a4067 100644 --- a/opening-hdf5/src/CMakeLists.txt +++ b/opening-hdf5/src/CMakeLists.txt @@ -4,6 +4,10 @@ project (ReadHDF5) 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) From f8bcf177307b29c6d9dab373299203c46737f237 Mon Sep 17 00:00:00 2001 From: djairoh Date: Thu, 18 Apr 2024 15:51:02 +0200 Subject: [PATCH 06/15] fix: upload README update --- opening-hdf5/README.md | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/opening-hdf5/README.md b/opening-hdf5/README.md index 44aca1b..890a39d 100644 --- a/opening-hdf5/README.md +++ b/opening-hdf5/README.md @@ -20,5 +20,16 @@ make ``` ### Building with Linux -Makes use of `mdspan` which is not supported by GCC libstdc++ at time of writing. See [compiler support](https://en.cppreference.com/w/cpp/compiler_support/23) for `mdspan`. Probably you will need to install Clang libc++. +Makes use of `mdspan` which is not supported by GCC libstdc++ 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 derivatves 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. Build the source code with the following flags: + ```sh + mkdir build && cd build + cmake .. -DCMAKE_CXX_COMPILER=/usr/bin/clang++ + make + ctest + make install + ``` + 3. Now the code should compile through the standard steps described in the Compiling section. From 1eac15304dd4427a7f46807360c557b180f8a6e5 Mon Sep 17 00:00:00 2001 From: djairoh Date: Fri, 19 Apr 2024 10:20:32 +0200 Subject: [PATCH 07/15] fix: forgot a line under compiling for linux --- opening-hdf5/README.md | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/opening-hdf5/README.md b/opening-hdf5/README.md index 890a39d..ec2a268 100644 --- a/opening-hdf5/README.md +++ b/opening-hdf5/README.md @@ -20,16 +20,20 @@ make ``` ### Building with Linux -Makes use of `mdspan` which is not supported by GCC libstdc++ 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 derivatves 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. +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. Build the source code with the following flags: + 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 - make install + sudo make install ``` 3. Now the code should compile through the standard steps described in the Compiling section. From 4cd545c78833df0d4cf362c8189fbc073ccf3dea Mon Sep 17 00:00:00 2001 From: robin Date: Sun, 21 Apr 2024 12:35:46 +0200 Subject: [PATCH 08/15] updated gitignore --- opening-hdf5/.gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opening-hdf5/.gitignore b/opening-hdf5/.gitignore index 13c6291..95abd3d 100644 --- a/opening-hdf5/.gitignore +++ b/opening-hdf5/.gitignore @@ -2,3 +2,5 @@ src/.DS_Store src/.cache src/build +.idea +src/cmake-build-debug \ No newline at end of file From a123071fec7743c250629c61f9e1ebae2aa83968 Mon Sep 17 00:00:00 2001 From: robin Date: Sun, 21 Apr 2024 14:20:23 +0200 Subject: [PATCH 09/15] feat: modules --- linear-interpolation/.gitignore | 6 +++ linear-interpolation/README.md | 39 +++++++++++++++ linear-interpolation/src/CMakeLists.txt | 36 ++++++++++++++ linear-interpolation/src/main.cpp | 18 +++++++ linear-interpolation/src/readdata.cpp | 66 +++++++++++++++++++++++++ linear-interpolation/src/readdata.h | 18 +++++++ linear-interpolation/src/vecutils.cpp | 14 ++++++ linear-interpolation/src/vecutils.h | 47 ++++++++++++++++++ 8 files changed, 244 insertions(+) create mode 100644 linear-interpolation/.gitignore create mode 100644 linear-interpolation/README.md create mode 100644 linear-interpolation/src/CMakeLists.txt create mode 100644 linear-interpolation/src/main.cpp create mode 100644 linear-interpolation/src/readdata.cpp create mode 100644 linear-interpolation/src/readdata.h create mode 100644 linear-interpolation/src/vecutils.cpp create mode 100644 linear-interpolation/src/vecutils.h diff --git a/linear-interpolation/.gitignore b/linear-interpolation/.gitignore new file mode 100644 index 0000000..95abd3d --- /dev/null +++ b/linear-interpolation/.gitignore @@ -0,0 +1,6 @@ +.DS_Store +src/.DS_Store +src/.cache +src/build +.idea +src/cmake-build-debug \ No newline at end of file diff --git a/linear-interpolation/README.md b/linear-interpolation/README.md new file mode 100644 index 0000000..ec2a268 --- /dev/null +++ b/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. diff --git a/linear-interpolation/src/CMakeLists.txt b/linear-interpolation/src/CMakeLists.txt new file mode 100644 index 0000000..b8b6abb --- /dev/null +++ b/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) + +# 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) + +add_executable(LinearInterpolate main.cpp + readdata.cpp + readdata.h + vecutils.cpp + vecutils.h) + +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}) diff --git a/linear-interpolation/src/main.cpp b/linear-interpolation/src/main.cpp new file mode 100644 index 0000000..8792393 --- /dev/null +++ b/linear-interpolation/src/main.cpp @@ -0,0 +1,18 @@ +#include "readdata.h" + +#include + +using namespace std; + +int main() { + auto [vec, size] = readHydrodynamicU(); + + auto arr = std::mdspan(vec.data(), size); + + print3DMatrixSlice(arr, 100); + + auto [times, lats, longs] = readGrid(); + printContentsOfVec(lats); + + return 0; +} diff --git a/linear-interpolation/src/readdata.cpp b/linear-interpolation/src/readdata.cpp new file mode 100644 index 0000000..6cc7726 --- /dev/null +++ b/linear-interpolation/src/readdata.cpp @@ -0,0 +1,66 @@ +#include +#include + +#include + +#include "readdata.h" + +using namespace std; +using namespace netCDF; + +template +vector getVarVector(const NcVar &var) { + int length = 1; + for (NcDim dim : var.getDims()) { + length *= dim.getSize(); + } + + vector vec(length); + + var.getVar(vec.data()); + + return vec; +} + +/** + * 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 +pair, std::dextents> 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 vec(timeLength*latLength*longLength); + var.getVar(vec.data()); + auto arr = std::mdspan(vec.data(), timeLength, latLength, longLength); + + return {vec, arr.extents()}; +} + +pair, std::dextents> readHydrodynamicU() { + netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); + + multimap< string, NcVar > vars = data.getVars(); + + return get3DMat(vars.find("uo")->second); +} + +tuple, vector, vector> readGrid() { + netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read); + multimap< string, NcVar > vars = data.getVars(); + vector time = getVarVector(vars.find("times")->second); + vector longitude = getVarVector(vars.find("longitude")->second); + vector latitude = getVarVector(vars.find("latitude")->second); + + return {time, latitude, longitude}; +} \ No newline at end of file diff --git a/linear-interpolation/src/readdata.h b/linear-interpolation/src/readdata.h new file mode 100644 index 0000000..c5effb5 --- /dev/null +++ b/linear-interpolation/src/readdata.h @@ -0,0 +1,18 @@ +#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 + */ +std::pair, std::dextents> readHydrodynamicU(); + +/** + * Reads the file grid.h5 + * @return a tuple of (times, latitude, longitude) + */ +std::tuple, std::vector, std::vector> readGrid(); + +#endif //LINEARINTERPOLATE_READDATA_H diff --git a/linear-interpolation/src/vecutils.cpp b/linear-interpolation/src/vecutils.cpp new file mode 100644 index 0000000..f367b3d --- /dev/null +++ b/linear-interpolation/src/vecutils.cpp @@ -0,0 +1,14 @@ +#include + +#include "vecutils.h" + +using namespace std; + +void print3DMatrixSlice(const arr3d &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(""); + } +} diff --git a/linear-interpolation/src/vecutils.h b/linear-interpolation/src/vecutils.h new file mode 100644 index 0000000..5ddeec6 --- /dev/null +++ b/linear-interpolation/src/vecutils.h @@ -0,0 +1,47 @@ +#ifndef LINEARINTERPOLATE_VECUTILS_H +#define LINEARINTERPOLATE_VECUTILS_H + +#include +#include + +template +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 +void printContentsOfVec(const std::vector& 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 +void print3DMatrixSlice(const arr3d &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 &arr, int t); + +#endif //LINEARINTERPOLATE_VECUTILS_H From 51fe302920978757c42243e872c381b65c5a38ed Mon Sep 17 00:00:00 2001 From: robin Date: Mon, 22 Apr 2024 00:21:11 +0200 Subject: [PATCH 10/15] Implement bilinear interpolation --- linear-interpolation/src/CMakeLists.txt | 8 ++++- linear-interpolation/src/UVGrid.cpp | 28 +++++++++++++++ linear-interpolation/src/UVGrid.h | 30 +++++++++++++++++ linear-interpolation/src/interpolate.cpp | 43 ++++++++++++++++++++++++ linear-interpolation/src/interpolate.h | 12 +++++++ linear-interpolation/src/main.cpp | 15 +++------ linear-interpolation/src/point.cpp | 11 ++++++ linear-interpolation/src/point.h | 27 +++++++++++++++ linear-interpolation/src/readdata.cpp | 12 +++++-- linear-interpolation/src/readdata.h | 10 ++++-- linear-interpolation/src/vecutils.cpp | 10 ++++++ linear-interpolation/src/vecutils.h | 1 + 12 files changed, 192 insertions(+), 15 deletions(-) create mode 100644 linear-interpolation/src/UVGrid.cpp create mode 100644 linear-interpolation/src/UVGrid.h create mode 100644 linear-interpolation/src/interpolate.cpp create mode 100644 linear-interpolation/src/interpolate.h create mode 100644 linear-interpolation/src/point.cpp create mode 100644 linear-interpolation/src/point.h diff --git a/linear-interpolation/src/CMakeLists.txt b/linear-interpolation/src/CMakeLists.txt index b8b6abb..e293e85 100644 --- a/linear-interpolation/src/CMakeLists.txt +++ b/linear-interpolation/src/CMakeLists.txt @@ -16,7 +16,13 @@ add_executable(LinearInterpolate main.cpp readdata.cpp readdata.h vecutils.cpp - vecutils.h) + vecutils.h + interpolate.cpp + interpolate.h + UVGrid.cpp + UVGrid.h + point.h + point.cpp) execute_process( COMMAND nc-config --includedir diff --git a/linear-interpolation/src/UVGrid.cpp b/linear-interpolation/src/UVGrid.cpp new file mode 100644 index 0000000..aa493eb --- /dev/null +++ b/linear-interpolation/src/UVGrid.cpp @@ -0,0 +1,28 @@ +#include +#include + +#include "UVGrid.h" +#include "readdata.h" + +using namespace std; + +UVGrid::UVGrid() { + auto [us, sizeU] = readHydrodynamicU(); + auto [vs, sizeV] = readHydrodynamicV(); + // Assuming sizeU == sizeV + uvData = views::zip(us, vs) | ranges::to(); + uvMatrix = mdspan(uvData.data(), sizeU); + tie(times, lats, lons) = readGrid(); +} + +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]; +} \ No newline at end of file diff --git a/linear-interpolation/src/UVGrid.h b/linear-interpolation/src/UVGrid.h new file mode 100644 index 0000000..23642f3 --- /dev/null +++ b/linear-interpolation/src/UVGrid.h @@ -0,0 +1,30 @@ +#ifndef LINEARINTERPOLATE_UVGRID_H +#define LINEARINTERPOLATE_UVGRID_H + +#include +#include "vecutils.h" +#include "point.h" + +class UVGrid { +private: + /** + * u == Eastward Current Velocity in the Water Column + * v == Northward Current Velocity in the Water Column + */ + std::vector uvData; +public: + UVGrid(); + + // The step functions assume regular spacing + double lonStep() const; + double latStep() const; + int timeStep() const; + + std::vector times; + std::vector lats; + std::vector lons; + + arr3d uvMatrix; +}; + +#endif //LINEARINTERPOLATE_UVGRID_H diff --git a/linear-interpolation/src/interpolate.cpp b/linear-interpolation/src/interpolate.cpp new file mode 100644 index 0000000..f39e801 --- /dev/null +++ b/linear-interpolation/src/interpolate.cpp @@ -0,0 +1,43 @@ +#include + +#include +#include "interpolate.h" + +using namespace std; + +// Inspired by https://numerical.recipes/book.html Chapter 3.6 +point bilinearInterpolate(const UVGrid &uvGrid, std::tuple 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 lonIndex = (lon-uvGrid.lons[0])/lonStep; + int timeIndex = (time-uvGrid.times[0])/timeStep; + + double timeRatio = (static_cast(time)-uvGrid.times[timeIndex])/timeStep; + double latRatio = (lat-uvGrid.lats[latIndex]) / latStep; + double lonRatio = (lon-uvGrid.lons[lonIndex]) / lonStep; + + point point = {0, 0}; + for(int time = 0; time <= 1; time++) { + for(int lat = 0; lat <= 1; lat++) { + for(int lon = 0; lon <= 1; lon++) { + auto vertex = uvGrid.uvMatrix[ + timeIndex + time < uvGrid.uvMatrix.extent(0) ? timeIndex + 1 : timeIndex, + latIndex + lat < uvGrid.uvMatrix.extent(1) ? latIndex + 1 : latIndex, + lonIndex + lon < uvGrid.uvMatrix.extent(2) ? lonIndex + 1 : lonIndex + ]; + + double timeRation = (1 - time)*(1-timeRatio) + time*timeRatio; + double latRation = (1 - lat)*(1-latRatio) + lat*latRatio; + double lonRation = (1 - lon)*(1-lonRatio) + lon*lonRatio; + point += timeRation * latRation * lonRation * vertex; + } + } + } + + return point; +} diff --git a/linear-interpolation/src/interpolate.h b/linear-interpolation/src/interpolate.h new file mode 100644 index 0000000..131b747 --- /dev/null +++ b/linear-interpolation/src/interpolate.h @@ -0,0 +1,12 @@ +#ifndef LINEARINTERPOLATE_INTERPOLATE_H +#define LINEARINTERPOLATE_INTERPOLATE_H + +#include + +#include "UVGrid.h" + +point bilinearInterpolate(const UVGrid &uvGrid, std::tuple timeLatLong); + +std::vector bilinearInterpolate(); + +#endif //LINEARINTERPOLATE_INTERPOLATE_H diff --git a/linear-interpolation/src/main.cpp b/linear-interpolation/src/main.cpp index 8792393..dc36897 100644 --- a/linear-interpolation/src/main.cpp +++ b/linear-interpolation/src/main.cpp @@ -1,18 +1,13 @@ -#include "readdata.h" - -#include +#include "interpolate.h" using namespace std; int main() { - auto [vec, size] = readHydrodynamicU(); + UVGrid uvGrid; + auto p = bilinearInterpolate(uvGrid, {392400, 53, -14.5}); - auto arr = std::mdspan(vec.data(), size); - - print3DMatrixSlice(arr, 100); - - auto [times, lats, longs] = readGrid(); - printContentsOfVec(lats); + println("({}, {})", p.first, p.second); + p = bilinearInterpolate(uvGrid, {802400, 62, -14.5}); return 0; } diff --git a/linear-interpolation/src/point.cpp b/linear-interpolation/src/point.cpp new file mode 100644 index 0000000..a84800f --- /dev/null +++ b/linear-interpolation/src/point.cpp @@ -0,0 +1,11 @@ +#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; +} \ No newline at end of file diff --git a/linear-interpolation/src/point.h b/linear-interpolation/src/point.h new file mode 100644 index 0000000..4a22f35 --- /dev/null +++ b/linear-interpolation/src/point.h @@ -0,0 +1,27 @@ +#ifndef LINEARINTERPOLATE_POINT_H +#define LINEARINTERPOLATE_POINT_H + +#include + +using point = std::pair; // {u, v} + +point operator+(const point& p1, const point& p2); + +point& operator+=(point& p1, const point& p2); + +template +point operator*(const point& p, Scalar scalar) { + return {p.first * scalar, p.second * scalar}; +} + +template +point operator*(Scalar scalar, const point& p) { + return {p.first * scalar, p.second * scalar}; +} + +template +point operator/(const point& p, Scalar scalar) { + return {p.first / scalar, p.second / scalar}; +} + +#endif //LINEARINTERPOLATE_POINT_H diff --git a/linear-interpolation/src/readdata.cpp b/linear-interpolation/src/readdata.cpp index 6cc7726..bc4f364 100644 --- a/linear-interpolation/src/readdata.cpp +++ b/linear-interpolation/src/readdata.cpp @@ -55,10 +55,18 @@ pair, std::dextents> readHydrodynamicU() { return get3DMat(vars.find("uo")->second); } -tuple, vector, vector> readGrid() { +pair, std::dextents> readHydrodynamicV() { + netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read); + + multimap< string, NcVar > vars = data.getVars(); + + return get3DMat(vars.find("vo")->second); +} + +tuple, vector, vector> readGrid() { netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read); multimap< string, NcVar > vars = data.getVars(); - vector time = getVarVector(vars.find("times")->second); + vector time = getVarVector(vars.find("times")->second); vector longitude = getVarVector(vars.find("longitude")->second); vector latitude = getVarVector(vars.find("latitude")->second); diff --git a/linear-interpolation/src/readdata.h b/linear-interpolation/src/readdata.h index c5effb5..3008924 100644 --- a/linear-interpolation/src/readdata.h +++ b/linear-interpolation/src/readdata.h @@ -4,15 +4,21 @@ #include "vecutils.h" /** - * Reads the file hydrodynamic_U.h5 + * reads the file hydrodynamic_U.h5 * @return a pair of the data vector of the contents and its dimensions to be used with mdspan */ std::pair, std::dextents> 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 + */ +std::pair, std::dextents> readHydrodynamicV(); + /** * Reads the file grid.h5 * @return a tuple of (times, latitude, longitude) */ -std::tuple, std::vector, std::vector> readGrid(); +std::tuple, std::vector, std::vector> readGrid(); #endif //LINEARINTERPOLATE_READDATA_H diff --git a/linear-interpolation/src/vecutils.cpp b/linear-interpolation/src/vecutils.cpp index f367b3d..1b02df9 100644 --- a/linear-interpolation/src/vecutils.cpp +++ b/linear-interpolation/src/vecutils.cpp @@ -12,3 +12,13 @@ void print3DMatrixSlice(const arr3d &arr, int t) { println(""); } } + +void print3DMatrixSlice(const arr3d> &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(""); + } +} diff --git a/linear-interpolation/src/vecutils.h b/linear-interpolation/src/vecutils.h index 5ddeec6..8bf3691 100644 --- a/linear-interpolation/src/vecutils.h +++ b/linear-interpolation/src/vecutils.h @@ -43,5 +43,6 @@ void print3DMatrixSlice(const arr3d &arr, int t) { } void print3DMatrixSlice(const arr3d &arr, int t); +void print3DMatrixSlice(const arr3d> &arr, int t); #endif //LINEARINTERPOLATE_VECUTILS_H From b2415d8c4d361ceb5cb44145b9237ea3c9231953 Mon Sep 17 00:00:00 2001 From: robin Date: Tue, 23 Apr 2024 12:21:08 +0200 Subject: [PATCH 11/15] removed mdspan requirement --- .gitignore | 1 + linear-interpolation/.gitignore | 3 +- linear-interpolation/src/CMakeLists.txt | 10 +---- linear-interpolation/src/UVGrid.cpp | 52 ++++++++++++++++++++---- linear-interpolation/src/UVGrid.h | 39 ++++++++++++++---- linear-interpolation/src/Vel.cpp | 33 +++++++++++++++ linear-interpolation/src/Vel.h | 44 ++++++++++++++++++++ linear-interpolation/src/interpolate.cpp | 52 ++++++++++++++---------- linear-interpolation/src/interpolate.h | 20 ++++++++- linear-interpolation/src/main.cpp | 45 ++++++++++++++++++-- linear-interpolation/src/point.cpp | 11 ----- linear-interpolation/src/point.h | 27 ------------ linear-interpolation/src/readdata.cpp | 33 ++------------- linear-interpolation/src/readdata.h | 10 ++--- linear-interpolation/src/vecutils.cpp | 24 ----------- linear-interpolation/src/vecutils.h | 48 ---------------------- 16 files changed, 255 insertions(+), 197 deletions(-) create mode 100644 linear-interpolation/src/Vel.cpp create mode 100644 linear-interpolation/src/Vel.h delete mode 100644 linear-interpolation/src/point.cpp delete mode 100644 linear-interpolation/src/point.h delete mode 100644 linear-interpolation/src/vecutils.cpp delete mode 100644 linear-interpolation/src/vecutils.h diff --git a/.gitignore b/.gitignore index 5ca0973..08302ef 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ .DS_Store +.idea diff --git a/linear-interpolation/.gitignore b/linear-interpolation/.gitignore index 95abd3d..34ddd4b 100644 --- a/linear-interpolation/.gitignore +++ b/linear-interpolation/.gitignore @@ -3,4 +3,5 @@ src/.DS_Store src/.cache src/build .idea -src/cmake-build-debug \ No newline at end of file +src/cmake-build-debug +src/cmake-build-release diff --git a/linear-interpolation/src/CMakeLists.txt b/linear-interpolation/src/CMakeLists.txt index e293e85..9428e41 100644 --- a/linear-interpolation/src/CMakeLists.txt +++ b/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 - point.cpp) + Vel.h + Vel.cpp) execute_process( COMMAND nc-config --includedir diff --git a/linear-interpolation/src/UVGrid.cpp b/linear-interpolation/src/UVGrid.cpp index aa493eb..16eabc6 100644 --- a/linear-interpolation/src/UVGrid.cpp +++ b/linear-interpolation/src/UVGrid.cpp @@ -1,18 +1,42 @@ -#include #include +#include #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 [vs, sizeV] = readHydrodynamicV(); - // Assuming sizeU == sizeV - uvData = views::zip(us, vs) | ranges::to(); - uvMatrix = mdspan(uvData.data(), sizeU); + 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 = views::zip(us, vs) + | views::transform([](auto pair) { + return Vel(pair); + }) + | ranges::to(); +} + +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]; -} \ No newline at end of file + 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(""); + } +} diff --git a/linear-interpolation/src/UVGrid.h b/linear-interpolation/src/UVGrid.h index 23642f3..613bd77 100644 --- a/linear-interpolation/src/UVGrid.h +++ b/linear-interpolation/src/UVGrid.h @@ -2,29 +2,54 @@ #define LINEARINTERPOLATE_UVGRID_H #include -#include "vecutils.h" -#include "point.h" +#include "Vel.h" class UVGrid { private: /** - * u == Eastward Current Velocity in the Water Column - * v == Northward Current Velocity in the Water Column + * 1D data vector of all the us and vs */ - std::vector uvData; + std::vector 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 times; std::vector lats; std::vector lons; - arr3d 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 diff --git a/linear-interpolation/src/Vel.cpp b/linear-interpolation/src/Vel.cpp new file mode 100644 index 0000000..803e8f2 --- /dev/null +++ b/linear-interpolation/src/Vel.cpp @@ -0,0 +1,33 @@ +#include "Vel.h" +#include + +Vel::Vel(double u, double v) : u(u), v(v) {} + +Vel::Vel(const std::pair& p) : u(p.first), v(p.second) {} + +Vel& Vel::operator=(const std::pair& 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 +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; +} diff --git a/linear-interpolation/src/Vel.h b/linear-interpolation/src/Vel.h new file mode 100644 index 0000000..ec7ce52 --- /dev/null +++ b/linear-interpolation/src/Vel.h @@ -0,0 +1,44 @@ +#ifndef LINEARINTERPOLATE_VEL_H +#define LINEARINTERPOLATE_VEL_H + +#include +#include +#include +#include + +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& p); // Conversion constructor + Vel& operator=(const std::pair& 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 + Vel operator*(Scalar scalar) const { + return Vel(u * scalar, v * scalar); + } + + // Operator / to divide Vel by a scalar, defined as a member template + template + 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 +Vel operator*(Scalar scalar, const Vel& p) { + return Vel(p.u * scalar, p.v * scalar); +} + +#endif //LINEARINTERPOLATE_VEL_H diff --git a/linear-interpolation/src/interpolate.cpp b/linear-interpolation/src/interpolate.cpp index f39e801..579e241 100644 --- a/linear-interpolation/src/interpolate.cpp +++ b/linear-interpolation/src/interpolate.cpp @@ -1,39 +1,37 @@ #include - #include +#include + #include "interpolate.h" using namespace std; -// Inspired by https://numerical.recipes/book.html Chapter 3.6 -point bilinearInterpolate(const UVGrid &uvGrid, std::tuple timeLatLon) { - auto [time, lat, lon] = timeLatLon; - +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; + 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(time)-uvGrid.times[timeIndex])/timeStep; - double latRatio = (lat-uvGrid.lats[latIndex]) / latStep; - double lonRatio = (lon-uvGrid.lons[lonIndex]) / lonStep; + double timeRatio = (static_cast(time) - uvGrid.times[timeIndex]) / timeStep; + double latRatio = (lat - uvGrid.lats[latIndex]) / latStep; + double lonRatio = (lon - uvGrid.lons[lonIndex]) / lonStep; - point point = {0, 0}; - for(int time = 0; time <= 1; time++) { - for(int lat = 0; lat <= 1; lat++) { - for(int lon = 0; lon <= 1; lon++) { - auto vertex = uvGrid.uvMatrix[ - timeIndex + time < uvGrid.uvMatrix.extent(0) ? timeIndex + 1 : timeIndex, - latIndex + lat < uvGrid.uvMatrix.extent(1) ? latIndex + 1 : latIndex, - lonIndex + lon < uvGrid.uvMatrix.extent(2) ? lonIndex + 1 : lonIndex + 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 - time)*(1-timeRatio) + time*timeRatio; - double latRation = (1 - lat)*(1-latRatio) + lat*latRatio; - double lonRation = (1 - lon)*(1-lonRatio) + lon*lonRatio; + 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; } } @@ -41,3 +39,13 @@ point bilinearInterpolate(const UVGrid &uvGrid, std::tuple return point; } + +vector bilinearInterpolate(const UVGrid &uvGrid, vector> 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>(); + return results; +} diff --git a/linear-interpolation/src/interpolate.h b/linear-interpolation/src/interpolate.h index 131b747..68b58f7 100644 --- a/linear-interpolation/src/interpolate.h +++ b/linear-interpolation/src/interpolate.h @@ -5,8 +5,24 @@ #include "UVGrid.h" -point bilinearInterpolate(const UVGrid &uvGrid, std::tuple 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 bilinearInterpolate(); +/** + * Helper function for bilnearly interpolating a vector of points + * @param uvGrid velocity grid + * @param points vector of points + * @return interpolated velocities + */ +std::vector bilinearInterpolate(const UVGrid &uvGrid, std::vector> points); #endif //LINEARINTERPOLATE_INTERPOLATE_H diff --git a/linear-interpolation/src/main.cpp b/linear-interpolation/src/main.cpp index dc36897..8d2c16e 100644 --- a/linear-interpolation/src/main.cpp +++ b/linear-interpolation/src/main.cpp @@ -1,13 +1,52 @@ #include "interpolate.h" +#include "Vel.h" +#include +#include using namespace std; int main() { UVGrid uvGrid; - auto p = bilinearInterpolate(uvGrid, {392400, 53, -14.5}); + uvGrid.printSlice(100); - println("({}, {})", p.first, p.second); - p = bilinearInterpolate(uvGrid, {802400, 62, -14.5}); + 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; + + // 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> 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(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; } diff --git a/linear-interpolation/src/point.cpp b/linear-interpolation/src/point.cpp deleted file mode 100644 index a84800f..0000000 --- a/linear-interpolation/src/point.cpp +++ /dev/null @@ -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; -} \ No newline at end of file diff --git a/linear-interpolation/src/point.h b/linear-interpolation/src/point.h deleted file mode 100644 index 4a22f35..0000000 --- a/linear-interpolation/src/point.h +++ /dev/null @@ -1,27 +0,0 @@ -#ifndef LINEARINTERPOLATE_POINT_H -#define LINEARINTERPOLATE_POINT_H - -#include - -using point = std::pair; // {u, v} - -point operator+(const point& p1, const point& p2); - -point& operator+=(point& p1, const point& p2); - -template -point operator*(const point& p, Scalar scalar) { - return {p.first * scalar, p.second * scalar}; -} - -template -point operator*(Scalar scalar, const point& p) { - return {p.first * scalar, p.second * scalar}; -} - -template -point operator/(const point& p, Scalar scalar) { - return {p.first / scalar, p.second / scalar}; -} - -#endif //LINEARINTERPOLATE_POINT_H diff --git a/linear-interpolation/src/readdata.cpp b/linear-interpolation/src/readdata.cpp index bc4f364..d556c6f 100644 --- a/linear-interpolation/src/readdata.cpp +++ b/linear-interpolation/src/readdata.cpp @@ -22,45 +22,20 @@ vector getVarVector(const NcVar &var) { return vec; } -/** - * 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 -pair, std::dextents> 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 vec(timeLength*latLength*longLength); - var.getVar(vec.data()); - auto arr = std::mdspan(vec.data(), timeLength, latLength, longLength); - - return {vec, arr.extents()}; -} - -pair, std::dextents> readHydrodynamicU() { +vector readHydrodynamicU() { netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); multimap< string, NcVar > vars = data.getVars(); - return get3DMat(vars.find("uo")->second); + return getVarVector(vars.find("uo")->second); } -pair, std::dextents> readHydrodynamicV() { +vector readHydrodynamicV() { netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read); multimap< string, NcVar > vars = data.getVars(); - return get3DMat(vars.find("vo")->second); + return getVarVector(vars.find("vo")->second); } tuple, vector, vector> readGrid() { diff --git a/linear-interpolation/src/readdata.h b/linear-interpolation/src/readdata.h index 3008924..6e4c2c9 100644 --- a/linear-interpolation/src/readdata.h +++ b/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::dextents> readHydrodynamicU(); +std::vector 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::dextents> readHydrodynamicV(); +std::vector readHydrodynamicV(); /** * Reads the file grid.h5 diff --git a/linear-interpolation/src/vecutils.cpp b/linear-interpolation/src/vecutils.cpp deleted file mode 100644 index 1b02df9..0000000 --- a/linear-interpolation/src/vecutils.cpp +++ /dev/null @@ -1,24 +0,0 @@ -#include - -#include "vecutils.h" - -using namespace std; - -void print3DMatrixSlice(const arr3d &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> &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(""); - } -} diff --git a/linear-interpolation/src/vecutils.h b/linear-interpolation/src/vecutils.h deleted file mode 100644 index 8bf3691..0000000 --- a/linear-interpolation/src/vecutils.h +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef LINEARINTERPOLATE_VECUTILS_H -#define LINEARINTERPOLATE_VECUTILS_H - -#include -#include - -template -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 -void printContentsOfVec(const std::vector& 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 -void print3DMatrixSlice(const arr3d &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 &arr, int t); -void print3DMatrixSlice(const arr3d> &arr, int t); - -#endif //LINEARINTERPOLATE_VECUTILS_H From d06f287a91fa5b21ea3a201f21feefdb530a2161 Mon Sep 17 00:00:00 2001 From: robin Date: Tue, 23 Apr 2024 12:40:40 +0200 Subject: [PATCH 12/15] removed opening-hdf5 directory --- opening-hdf5/.gitignore | 6 -- opening-hdf5/README.md | 39 ------------- opening-hdf5/src/CMakeLists.txt | 32 ---------- opening-hdf5/src/main.cpp | 100 -------------------------------- 4 files changed, 177 deletions(-) delete mode 100644 opening-hdf5/.gitignore delete mode 100644 opening-hdf5/README.md delete mode 100644 opening-hdf5/src/CMakeLists.txt delete mode 100644 opening-hdf5/src/main.cpp diff --git a/opening-hdf5/.gitignore b/opening-hdf5/.gitignore deleted file mode 100644 index 95abd3d..0000000 --- a/opening-hdf5/.gitignore +++ /dev/null @@ -1,6 +0,0 @@ -.DS_Store -src/.DS_Store -src/.cache -src/build -.idea -src/cmake-build-debug \ No newline at end of file diff --git a/opening-hdf5/README.md b/opening-hdf5/README.md deleted file mode 100644 index ec2a268..0000000 --- a/opening-hdf5/README.md +++ /dev/null @@ -1,39 +0,0 @@ -## 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. diff --git a/opening-hdf5/src/CMakeLists.txt b/opening-hdf5/src/CMakeLists.txt deleted file mode 100644 index e1a4067..0000000 --- a/opening-hdf5/src/CMakeLists.txt +++ /dev/null @@ -1,32 +0,0 @@ -cmake_minimum_required (VERSION 3.28) -project (ReadHDF5) - -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) - -add_executable(ReadHDF5 main.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(ReadHDF5 PUBLIC ${netCDF_INCLUDE_DIR}) - -find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH) -target_link_libraries(ReadHDF5 ${NETCDF_LIB}) diff --git a/opening-hdf5/src/main.cpp b/opening-hdf5/src/main.cpp deleted file mode 100644 index 68e0ee8..0000000 --- a/opening-hdf5/src/main.cpp +++ /dev/null @@ -1,100 +0,0 @@ -#include -#include -#include -#include - -#include - -using namespace std; -using namespace netCDF; - -template -void printContentsOfVec(const std::vector& vec) { - for (const auto& element : vec) { - cout << element << " "; - } - cout << std::endl; -} - -template -vector getVarVector(NcVar var) { - int length = 1; - for (NcDim dim : var.getDims()) { - length *= dim.getSize(); - } - - vector vec(length); - - var.getVar(vec.data()); - - return vec; -} - -template -using arr3d = std::mdspan< - T, - std::dextents< - std::size_t, - 3 - > - >; - - -template -void print3DMatrixSlice(arr3d arr, int t) { - for (int x = 0; x < arr.extent(1); x++) { - for (int y = 0; y < arr.extent(2); y++) { - print("{} ", arr[t,x,y]); - } - println(""); - } -} - -void print3DMatrixSlice(arr3d 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(""); - } -} - -/** - * 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 -pair, std::dextents> get3DMat(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 vec(timeLength*latLength*longLength); - var.getVar(vec.data()); - auto arr = std::mdspan(vec.data(), timeLength, latLength, longLength); - - return {vec, arr.extents()}; -} - -int main () { - netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read); - - multimap< string, NcVar > vars = data.getVars(); - - auto [vec, size] = get3DMat(vars.find("uo")->second); - - - auto arr = std::mdspan(vec.data(), size); - - print3DMatrixSlice(arr, 100); - - return 0; -} From 855ac4d65493091c4e6940d800412227ab94426b Mon Sep 17 00:00:00 2001 From: djairoh Date: Tue, 23 Apr 2024 15:02:06 +0200 Subject: [PATCH 13/15] fix: manual implementation of ranges::to as this is not supported by gcc 13 --- linear-interpolation/src/CMakeLists.txt | 5 +++++ linear-interpolation/src/UVGrid.cpp | 17 +++++++++-------- linear-interpolation/src/interpolate.cpp | 10 ++++------ linear-interpolation/src/main.cpp | 6 +++--- linear-interpolation/src/readdata.cpp | 3 +-- 5 files changed, 22 insertions(+), 19 deletions(-) diff --git a/linear-interpolation/src/CMakeLists.txt b/linear-interpolation/src/CMakeLists.txt index 9428e41..faeeaa7 100644 --- a/linear-interpolation/src/CMakeLists.txt +++ b/linear-interpolation/src/CMakeLists.txt @@ -6,6 +6,10 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +# Make flags for release compilation +# set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -DNDEBUG") + + find_package(netCDF REQUIRED) add_executable(LinearInterpolate main.cpp @@ -16,6 +20,7 @@ add_executable(LinearInterpolate main.cpp UVGrid.cpp UVGrid.h Vel.h + to_vector.h Vel.cpp) execute_process( diff --git a/linear-interpolation/src/UVGrid.cpp b/linear-interpolation/src/UVGrid.cpp index 16eabc6..f209a12 100644 --- a/linear-interpolation/src/UVGrid.cpp +++ b/linear-interpolation/src/UVGrid.cpp @@ -1,8 +1,8 @@ #include -#include #include "UVGrid.h" #include "readdata.h" +#include "to_vector.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" @@ -27,11 +27,12 @@ UVGrid::UVGrid() { throw domain_error(sizeError); } - uvData = views::zip(us, vs) - | views::transform([](auto pair) { - return Vel(pair); - }) - | ranges::to(); + uvData = to_vector(views::transform(views::zip(us, vs), [](auto pair){return Vel(pair);})); +// uvData = views::zip(us, vs) +// | views::transform([](auto pair) { +// return Vel(pair); +// }) +// | ranges::to(); } const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const { @@ -55,8 +56,8 @@ 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); + printf("%7.4f, %7.4f", u, v); } - println(""); + cout << endl; } } diff --git a/linear-interpolation/src/interpolate.cpp b/linear-interpolation/src/interpolate.cpp index 579e241..26c46b3 100644 --- a/linear-interpolation/src/interpolate.cpp +++ b/linear-interpolation/src/interpolate.cpp @@ -1,8 +1,8 @@ #include #include -#include #include "interpolate.h" +#include "to_vector.h" using namespace std; @@ -41,11 +41,9 @@ Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon) } vector bilinearInterpolate(const UVGrid &uvGrid, vector> points) { - auto results = points - | std::views::transform([&uvGrid](const auto &point) { + auto results = points | std::views::transform([&uvGrid](const auto &point) { auto [time, lat, lon] = point; return bilinearInterpolate(uvGrid, time, lat, lon); - }) - | std::ranges::to>(); - return results; + }); + return to_vector(results); } diff --git a/linear-interpolation/src/main.cpp b/linear-interpolation/src/main.cpp index 8d2c16e..ebcd487 100644 --- a/linear-interpolation/src/main.cpp +++ b/linear-interpolation/src/main.cpp @@ -1,7 +1,7 @@ #include "interpolate.h" #include "Vel.h" -#include #include +#include using namespace std; @@ -39,14 +39,14 @@ int main() { auto duration = std::chrono::duration_cast(stop - start); - println("Time taken for {} points was {} seconds", N, duration.count()/1000.); + printf("Time taken for %d points was %lf seconds\n", 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); + printf("Sum = %lf\n", sum); return 0; } diff --git a/linear-interpolation/src/readdata.cpp b/linear-interpolation/src/readdata.cpp index d556c6f..85a56e8 100644 --- a/linear-interpolation/src/readdata.cpp +++ b/linear-interpolation/src/readdata.cpp @@ -1,4 +1,3 @@ -#include #include #include @@ -46,4 +45,4 @@ tuple, vector, vector> readGrid() { vector latitude = getVarVector(vars.find("latitude")->second); return {time, latitude, longitude}; -} \ No newline at end of file +} From 6288a43da7f951c6e631dbf3524f22b8864f9c1d Mon Sep 17 00:00:00 2001 From: djairoh Date: Tue, 23 Apr 2024 15:02:21 +0200 Subject: [PATCH 14/15] fix: forgot to include a file --- linear-interpolation/src/to_vector.h | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 linear-interpolation/src/to_vector.h diff --git a/linear-interpolation/src/to_vector.h b/linear-interpolation/src/to_vector.h new file mode 100644 index 0000000..989085b --- /dev/null +++ b/linear-interpolation/src/to_vector.h @@ -0,0 +1,24 @@ +#ifndef TO_VECTOR_H +#define TO_VECTOR_H + +#include +#include + +template +auto to_vector(R&& r) { + std::vector> v; + + // if we can get a size, reserve that much + if constexpr (requires { std::ranges::size(r); }) { + v.reserve(std::ranges::size(r)); + } + + // push all the elements + for (auto&& e : r) { + v.push_back(static_cast(e)); + } + + return v; +} + +#endif From ea23ed9d68210e2d0132f3dad5e4b8b1d93404bf Mon Sep 17 00:00:00 2001 From: robin Date: Tue, 23 Apr 2024 17:41:02 +0200 Subject: [PATCH 15/15] removed some functional stuff --- linear-interpolation/src/CMakeLists.txt | 5 ----- linear-interpolation/src/UVGrid.cpp | 20 +++++++++----------- linear-interpolation/src/UVGrid.h | 6 +----- linear-interpolation/src/Vel.cpp | 13 ++++++++++--- linear-interpolation/src/interpolate.cpp | 16 +++++++--------- linear-interpolation/src/main.cpp | 13 ++++++------- linear-interpolation/src/to_vector.h | 24 ------------------------ 7 files changed, 33 insertions(+), 64 deletions(-) delete mode 100644 linear-interpolation/src/to_vector.h diff --git a/linear-interpolation/src/CMakeLists.txt b/linear-interpolation/src/CMakeLists.txt index faeeaa7..9428e41 100644 --- a/linear-interpolation/src/CMakeLists.txt +++ b/linear-interpolation/src/CMakeLists.txt @@ -6,10 +6,6 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) -# Make flags for release compilation -# set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -DNDEBUG") - - find_package(netCDF REQUIRED) add_executable(LinearInterpolate main.cpp @@ -20,7 +16,6 @@ add_executable(LinearInterpolate main.cpp UVGrid.cpp UVGrid.h Vel.h - to_vector.h Vel.cpp) execute_process( diff --git a/linear-interpolation/src/UVGrid.cpp b/linear-interpolation/src/UVGrid.cpp index f209a12..6673237 100644 --- a/linear-interpolation/src/UVGrid.cpp +++ b/linear-interpolation/src/UVGrid.cpp @@ -2,7 +2,6 @@ #include "UVGrid.h" #include "readdata.h" -#include "to_vector.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" @@ -27,12 +26,11 @@ UVGrid::UVGrid() { throw domain_error(sizeError); } - uvData = to_vector(views::transform(views::zip(us, vs), [](auto pair){return Vel(pair);})); -// uvData = views::zip(us, vs) -// | views::transform([](auto pair) { -// return Vel(pair); -// }) -// | ranges::to(); + 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 { @@ -52,12 +50,12 @@ int UVGrid::timeStep() const { return times[1] - times[0]; } -void UVGrid::printSlice(size_t t) { +void UVGrid::streamSlice(ostream &os, size_t t) { for (int x = 0; x < latSize; x++) { for (int y = 0; y < lonSize; y++) { - auto [u,v] = (*this)[t,x,y]; - printf("%7.4f, %7.4f", u, v); + auto vel = (*this)[t,x,y]; + os << vel << " "; } - cout << endl; + os << endl; } } diff --git a/linear-interpolation/src/UVGrid.h b/linear-interpolation/src/UVGrid.h index 613bd77..6c3f8f6 100644 --- a/linear-interpolation/src/UVGrid.h +++ b/linear-interpolation/src/UVGrid.h @@ -45,11 +45,7 @@ public: */ 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); + void streamSlice(std::ostream &os, size_t t); }; - #endif //LINEARINTERPOLATE_UVGRID_H diff --git a/linear-interpolation/src/Vel.cpp b/linear-interpolation/src/Vel.cpp index 803e8f2..6bb5990 100644 --- a/linear-interpolation/src/Vel.cpp +++ b/linear-interpolation/src/Vel.cpp @@ -1,5 +1,8 @@ #include "Vel.h" #include +#include + +using namespace std; Vel::Vel(double u, double v) : u(u), v(v) {} @@ -27,7 +30,11 @@ Vel Vel::operator/(Scalar scalar) const { return Vel(u / scalar, v / scalar); } -std::ostream& operator<<(std::ostream& os, const Vel& vel) { - os << "(" << vel.u << ", " << vel.v << ")"; +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; -} +} \ No newline at end of file diff --git a/linear-interpolation/src/interpolate.cpp b/linear-interpolation/src/interpolate.cpp index 26c46b3..98fe42d 100644 --- a/linear-interpolation/src/interpolate.cpp +++ b/linear-interpolation/src/interpolate.cpp @@ -1,8 +1,4 @@ -#include -#include - #include "interpolate.h" -#include "to_vector.h" using namespace std; @@ -41,9 +37,11 @@ Vel bilinearInterpolate(const UVGrid &uvGrid, int time, double lat, double lon) } vector bilinearInterpolate(const UVGrid &uvGrid, vector> points) { - auto results = points | std::views::transform([&uvGrid](const auto &point) { - auto [time, lat, lon] = point; - return bilinearInterpolate(uvGrid, time, lat, lon); - }); - return to_vector(results); + vector result; + result.reserve(points.size()); + for (auto [time, lat, lon]: points) { + result.push_back(bilinearInterpolate(uvGrid, time, lat, lon)); + } + + return result; } diff --git a/linear-interpolation/src/main.cpp b/linear-interpolation/src/main.cpp index ebcd487..77f68ff 100644 --- a/linear-interpolation/src/main.cpp +++ b/linear-interpolation/src/main.cpp @@ -7,7 +7,7 @@ using namespace std; int main() { UVGrid uvGrid; - uvGrid.printSlice(100); + uvGrid.streamSlice(cout, 100); int N = 10000000; // Number of points @@ -20,7 +20,6 @@ int main() { 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); @@ -31,22 +30,22 @@ int main() { 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 start = chrono::high_resolution_clock::now(); auto x = bilinearInterpolate(uvGrid, points); - auto stop = std::chrono::high_resolution_clock::now(); + auto stop = chrono::high_resolution_clock::now(); - auto duration = std::chrono::duration_cast(stop - start); + auto duration = chrono::duration_cast(stop - start); - printf("Time taken for %d points was %lf seconds\n", N, duration.count()/1000.); + 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; } - printf("Sum = %lf\n", sum); + cout << "Sum = " << sum << endl; return 0; } diff --git a/linear-interpolation/src/to_vector.h b/linear-interpolation/src/to_vector.h deleted file mode 100644 index 989085b..0000000 --- a/linear-interpolation/src/to_vector.h +++ /dev/null @@ -1,24 +0,0 @@ -#ifndef TO_VECTOR_H -#define TO_VECTOR_H - -#include -#include - -template -auto to_vector(R&& r) { - std::vector> v; - - // if we can get a size, reserve that much - if constexpr (requires { std::ranges::size(r); }) { - v.reserve(std::ranges::size(r)); - } - - // push all the elements - for (auto&& e : r) { - v.push_back(static_cast(e)); - } - - return v; -} - -#endif