From ef72d145f2accbda3f14ebb77a14ebde8d1c60f1 Mon Sep 17 00:00:00 2001 From: Rahil Doshi <rahil.doshi@fau.de> Date: Mon, 3 Feb 2025 08:08:09 +0100 Subject: [PATCH] Merge development into release/v0.1.0: Add binary search and double lookup interpolation - Implement new interpolation methods in C++ - Update core functionality for improved performance - Add test files for interpolation methods - Update build configuration - Clean code for release --- pyproject.toml | 1 + setup.py | 9 +- src/pymatlib/core/cpp/CMakeLists.txt | 50 ++- .../core/cpp/binary_search_interpolation.cpp | 311 ++++++++++++++++++ .../core/cpp/double_lookup_interpolation.cpp | 213 ++++++++++++ .../binary_search_interpolation.h} | 3 +- .../cpp/include/double_lookup_interpolation.h | 17 + .../include/interpolate_binary_search_cpp.h | 83 +++++ .../include/interpolate_double_lookup_cpp.h | 43 +++ src/pymatlib/core/cpp/module.cpp | 26 ++ .../temperature_from_energy_density_array.cpp | 89 ----- src/pymatlib/core/cpp/test_interpolation.cpp | 57 ++++ src/pymatlib/core/cpp/test_interpolation.py | 30 ++ src/pymatlib/core/data_handler.py | 79 +++++ src/pymatlib/core/interpolators.py | 187 ++++++++++- src/pymatlib/core/models.py | 84 ----- .../core/test_interpolate_functs_perf.py | 222 +++++++++++++ src/pymatlib/data/alloys/SS316L/SS316L.py | 10 +- .../SS316L/density_temperature_edited.txt | 301 +++++++++++++++++ 19 files changed, 1610 insertions(+), 205 deletions(-) create mode 100644 src/pymatlib/core/cpp/binary_search_interpolation.cpp create mode 100644 src/pymatlib/core/cpp/double_lookup_interpolation.cpp rename src/pymatlib/core/cpp/{temperature_from_energy_density_array.h => include/binary_search_interpolation.h} (89%) create mode 100644 src/pymatlib/core/cpp/include/double_lookup_interpolation.h create mode 100644 src/pymatlib/core/cpp/include/interpolate_binary_search_cpp.h create mode 100644 src/pymatlib/core/cpp/include/interpolate_double_lookup_cpp.h create mode 100644 src/pymatlib/core/cpp/module.cpp delete mode 100644 src/pymatlib/core/cpp/temperature_from_energy_density_array.cpp create mode 100644 src/pymatlib/core/cpp/test_interpolation.cpp create mode 100644 src/pymatlib/core/cpp/test_interpolation.py create mode 100644 src/pymatlib/core/test_interpolate_functs_perf.py create mode 100644 src/pymatlib/data/alloys/SS316L/density_temperature_edited.txt diff --git a/pyproject.toml b/pyproject.toml index 76a46a4..ddb7a47 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,6 +6,7 @@ authors = [ {name = "Rahil Doshi", email = "rahil.doshi@fau.de"}, ] dependencies = [ + "numpy>=1.18.0", "pystencils@git+https://i10git.cs.fau.de/pycodegen/pystencils.git@v2.0-dev" ] requires-python = ">=3.10" diff --git a/setup.py b/setup.py index ed646c5..594f94f 100644 --- a/setup.py +++ b/setup.py @@ -6,8 +6,13 @@ import pybind11 ext_modules = [ Extension( "pymatlib.core.cpp.fast_interpolation", # Module name in Python - ["src/pymatlib/core/cpp/temperature_from_energy_density_array.cpp"], - include_dirs=[pybind11.get_include()], + [ + "src/pymatlib/core/cpp/module.cpp", + "src/pymatlib/core/cpp/binary_search_interpolation.cpp", + "src/pymatlib/core/cpp/double_lookup_interpolation.cpp", + ], + include_dirs=[pybind11.get_include(), + "src/pymatlib/core/cpp/include"], extra_compile_args=['-O3', '-std=c++11'], # Enable high optimization and C++11 language='c++' ), diff --git a/src/pymatlib/core/cpp/CMakeLists.txt b/src/pymatlib/core/cpp/CMakeLists.txt index 333ab99..26261a1 100644 --- a/src/pymatlib/core/cpp/CMakeLists.txt +++ b/src/pymatlib/core/cpp/CMakeLists.txt @@ -1,13 +1,51 @@ cmake_minimum_required(VERSION 3.10) project(fast_interpolation) -find_package(Python COMPONENTS Development NumPy REQUIRED) +# Find required packages +find_package(Python COMPONENTS Interpreter Development NumPy REQUIRED) +execute_process( + COMMAND ${Python_EXECUTABLE} -m pybind11 --cmakedir + OUTPUT_VARIABLE PYBIND11_CMAKE_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE +) +find_package(pybind11 PATHS ${PYBIND11_CMAKE_DIR} NO_DEFAULT_PATH REQUIRED) -execute_process(COMMAND ${Python_EXECUTABLE} -m pybind11 --cmakedir OUTPUT_VARIABLE PYBIND11_CMAKE_PATH) -string(STRIP "${PYBIND11_CMAKE_PATH}" PYBIND11_CMAKE_PATH) -find_package(pybind11 PATHS "${PYBIND11_CMAKE_PATH}" NO_DEFAULT_PATH REQUIRED) +# Create interface library for header-only template implementation +add_library(interpolation_template INTERFACE) +target_include_directories(interpolation_template INTERFACE + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/include +) -pybind11_add_module(fast_interpolation temperature_from_energy_density_array.cpp) +# Add the Pybind11 module +pybind11_add_module(fast_interpolation + module.cpp + binary_search_interpolation.cpp + double_lookup_interpolation.cpp +) +target_include_directories(fast_interpolation PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/include +) target_compile_features(fast_interpolation PRIVATE cxx_std_11) target_compile_options(fast_interpolation PRIVATE -O3) -target_link_libraries(fast_interpolation PRIVATE Python::NumPy) +target_link_libraries(fast_interpolation PRIVATE + interpolation_template + pybind11::module + Python::NumPy +) + +# Add C++ test executable +add_executable(fast_interpolation_test_cpp test_interpolation.cpp) +target_compile_features(fast_interpolation_test_cpp PRIVATE cxx_std_11) +target_compile_options(fast_interpolation_test_cpp PRIVATE -O3) +target_link_libraries(fast_interpolation_test_cpp PRIVATE + interpolation_template +) + +# Add custom target for running Python tests +add_custom_target(fast_interpolation_test_py + COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_interpolation.py + DEPENDS fast_interpolation + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMENT "Running Python tests" +) diff --git a/src/pymatlib/core/cpp/binary_search_interpolation.cpp b/src/pymatlib/core/cpp/binary_search_interpolation.cpp new file mode 100644 index 0000000..97379d9 --- /dev/null +++ b/src/pymatlib/core/cpp/binary_search_interpolation.cpp @@ -0,0 +1,311 @@ +#include "include/binary_search_interpolation.h" +#include <iostream> +#include <cmath> + + +double interpolate_binary_search( + const py::array_t<double>& temperature_array, + double h_in, + const py::array_t<double>& energy_density_array) { + + static constexpr double EPSILON = 1e-6; + + // Get array views + auto temp_arr = temperature_array.unchecked<1>(); + auto energy_arr = energy_density_array.unchecked<1>(); + const size_t n = temp_arr.shape(0); + + // Input validation + if (temperature_array.size() != energy_density_array.size() || n < 2) { + throw std::runtime_error("Invalid array sizes"); + } + + // Determine array order and cache indices + const bool is_ascending = temp_arr(0) < temp_arr(n-1); + + const size_t start_idx = is_ascending ? 0 : n-1; + const size_t end_idx = is_ascending ? n-1 : 0; + + // Boundary checks + if (h_in <= energy_arr(start_idx)) return temp_arr(start_idx); + if (h_in >= energy_arr(end_idx)) return temp_arr(end_idx); + + // Binary search with optimized memory access + size_t left = 0; + size_t right = n - 1; + + while (left <= right) { + const size_t mid = (left + right) / 2; + const double mid_val = energy_arr(mid); + + if (std::abs(mid_val - h_in) < EPSILON) { + return temp_arr(mid); + } + + const bool go_left = (mid_val > h_in) == is_ascending; + if (go_left) { + right = mid - 1; + } else { + left = mid + 1; + } + } + + + // Linear interpolation + const double x0 = energy_arr(right); + const double x1 = energy_arr(left); + const double y0 = temp_arr(right); + const double y1 = temp_arr(left); + + return y0 + (y1 - y0) * (h_in - x0) / (x1 - x0); +} + +/* +double interpolate_binary_search( + const py::array_t<double>& temperature_array, + double h_in, + const py::array_t<double>& energy_density_array) { + + static constexpr double EPSILON = 1e-6; + + // Cache buffer requests and get pointers + const auto temp_buf = temperature_array.request(); + const auto energy_buf = energy_density_array.request(); + + // Use __restrict__ to inform compiler about non-aliasing + const double* __restrict__ temp_ptr = static_cast<double*>(temp_buf.ptr); + const double* __restrict__ energy_ptr = static_cast<double*>(energy_buf.ptr); + const size_t n = temp_buf.size; + + // Input validation + if (temp_buf.size != energy_buf.size || n < 2) { + throw std::runtime_error("Invalid array sizes"); + } + + // Determine array order and cache indices + const bool is_ascending = temp_ptr[0] < temp_ptr[n-1]; + const size_t start_idx = is_ascending ? 0 : n-1; + const size_t end_idx = is_ascending ? n-1 : 0; + std::cout << "C++ - is_ascending: " << is_ascending << std::endl; + std::cout << "C++ - h_in: " << h_in << std::endl; + std::cout << "C++ - First few T values: "; + for(size_t i = 0; i < 5; i++) std::cout << temp_ptr[i] << " "; + std::cout << std::endl; + std::cout << "C++ - First few E values: "; + for(size_t i = 0; i < 5; i++) std::cout << energy_ptr[i] << " "; + std::cout << std::endl; + + if (energy_ptr[start_idx] >= energy_ptr[end_idx]) { + throw std::runtime_error("Energy density must increase with temperature"); + } + + // Boundary checks + if (h_in <= energy_ptr[start_idx]) return temp_ptr[start_idx]; + if (h_in >= energy_ptr[end_idx]) return temp_ptr[end_idx]; + + // Binary search with optimized memory access + size_t left = 0; + size_t right = n - 1; + + while (left <= right) { + const size_t mid = (left + right) / 2; + const double mid_val = energy_ptr[mid]; + + if (std::abs(mid_val - h_in) < EPSILON) { + return temp_ptr[mid]; + } + + const bool go_left = (mid_val > h_in) == is_ascending; + if (go_left) { + right = mid - 1; + } else { + left = mid + 1; + } + } + + // Linear interpolation + const size_t idx1 = is_ascending ? right : left; + const size_t idx2 = is_ascending ? left : right; + const double x0 = energy_ptr[idx1]; + const double dx = energy_ptr[idx2] - x0; + const double y0 = temp_ptr[idx1]; + + return y0 + (temp_ptr[idx2] - y0) * (h_in - x0) / dx; +} +*/ +/* +double interpolate_binary_search( + const py::array_t<double>& temperature_array, + double h_in, + const py::array_t<double>& energy_density_array) { + + static constexpr double EPSILON = 1e-6; + + const auto temp_buf = temperature_array.request(); + const auto energy_buf = energy_density_array.request(); + + if (temp_buf.size != energy_buf.size) { + throw std::runtime_error("Arrays must have same length"); + } + + const auto temp_ptr = static_cast<double*>(temp_buf.ptr); + const auto energy_ptr = static_cast<double*>(energy_buf.ptr); + const size_t n = temp_buf.size; + + if (n < 2) { + throw std::runtime_error("Arrays must have at least 2 elements"); + } + + // Determine array order + const bool is_ascending = temp_ptr[0] < temp_ptr[n-1]; + const size_t start_idx = is_ascending ? 0 : n-1; + const size_t end_idx = is_ascending ? n-1 : 0; + + // Validate energy density increases with temperature + if (energy_ptr[start_idx] >= energy_ptr[end_idx]) { + throw std::runtime_error("Energy density must increase with temperature"); + } + + // Quick boundary checks + if (h_in <= energy_ptr[start_idx]) return temp_ptr[start_idx]; + if (h_in >= energy_ptr[end_idx]) return temp_ptr[end_idx]; + + // Binary search + size_t left = 0; + size_t right = n - 1; + + while (left <= right) { + const size_t mid = (left + right) / 2; + const double mid_val = energy_ptr[mid]; + + if (std::abs(mid_val - h_in) < EPSILON) { + return temp_ptr[mid]; + } + + if (mid_val > h_in) { + if (is_ascending) { + right = mid - 1; + } else { + left = mid + 1; + } + } else { + if (is_ascending) { + left = mid + 1; + } else { + right = mid - 1; + } + } + } + + // Linear interpolation + const size_t idx1 = is_ascending ? right : left; + const size_t idx2 = is_ascending ? left : right; + + const double x0 = energy_ptr[idx1]; + const double x1 = energy_ptr[idx2]; + const double y0 = temp_ptr[idx1]; + const double y1 = temp_ptr[idx2]; + + const double slope = (y1 - y0) / (x1 - x0); + return y0 + slope * (h_in - x0); +} +*/ +/* +namespace detail { + static constexpr double EPSILON = 1e-6; + + inline size_t interpolation_guess(const size_t left, const size_t right, + const double* energy_ptr, const double h_in) { + const double left_val = energy_ptr[left]; + const double right_val = energy_ptr[right]; + + // Interpolation guess + const double pos = static_cast<double>(left) + + ((h_in - left_val) * static_cast<double>(right - left)) / + (right_val - left_val); + + return static_cast<size_t>(std::round(pos)); + } +} + +double interpolate_binary_search( + const py::array_t<double>& temperature_array, + double h_in, + const py::array_t<double>& energy_density_array) { + + const auto temp_buf = temperature_array.request(); + const auto energy_buf = energy_density_array.request(); + + if (temp_buf.size != energy_buf.size) { + throw std::runtime_error("Arrays must have same length"); + } + + const double* temp_ptr = static_cast<double*>(temp_buf.ptr); + const double* energy_ptr = static_cast<double*>(energy_buf.ptr); + const size_t n = temp_buf.size; + + if (n < 2) { + throw std::runtime_error("Arrays must have at least 2 elements"); + } + + // Determine array order + const bool is_ascending = temp_ptr[0] < temp_ptr[n-1]; + const size_t start_idx = is_ascending ? 0 : n-1; + const size_t end_idx = is_ascending ? n-1 : 0; + + // Validate energy density increases with temperature + if (energy_ptr[start_idx] >= energy_ptr[end_idx]) { + throw std::runtime_error("Energy density must increase with temperature"); + } + + // Quick boundary checks + if (h_in <= energy_ptr[start_idx]) return temp_ptr[start_idx]; + if (h_in >= energy_ptr[end_idx]) return temp_ptr[end_idx]; + + // Search with interpolation and binary fallback + size_t left = 0; + size_t right = n - 1; + + while (left <= right) { + // Try interpolation search first + size_t mid = detail::interpolation_guess(left, right, energy_ptr, h_in); + + // Fallback to binary search if interpolation guess is out of bounds + if (mid < left || mid > right) { + mid = left + (right - left) / 2; + } + + const double mid_val = energy_ptr[mid]; + + if (std::abs(mid_val - h_in) < detail::EPSILON) { + return temp_ptr[mid]; + } + + if (mid_val > h_in) { + if (is_ascending) { + right = mid - 1; + } else { + left = mid + 1; + } + } else { + if (is_ascending) { + left = mid + 1; + } else { + right = mid - 1; + } + } + } + + // Linear interpolation + const size_t idx1 = is_ascending ? right : left; + const size_t idx2 = is_ascending ? left : right; + + const double x0 = energy_ptr[idx1]; + const double x1 = energy_ptr[idx2]; + const double y0 = temp_ptr[idx1]; + const double y1 = temp_ptr[idx2]; + + const double slope = (y1 - y0) / (x1 - x0); + return y0 + slope * (h_in - x0); +} +*/ \ No newline at end of file diff --git a/src/pymatlib/core/cpp/double_lookup_interpolation.cpp b/src/pymatlib/core/cpp/double_lookup_interpolation.cpp new file mode 100644 index 0000000..4f1673e --- /dev/null +++ b/src/pymatlib/core/cpp/double_lookup_interpolation.cpp @@ -0,0 +1,213 @@ +#include "include/double_lookup_interpolation.h" +#include <iostream> +#include <cmath> + + +double interpolate_double_lookup( + double E_target, + const py::array_t<double>& T_eq, + const py::array_t<double>& E_neq, + const py::array_t<double>& E_eq, + const py::array_t<int32_t>& idx_map) { + + // Get array views with bounds checking disabled for performance + auto T_eq_arr = T_eq.unchecked<1>(); + auto E_neq_arr = E_neq.unchecked<1>(); + auto E_eq_arr = E_eq.unchecked<1>(); + auto idx_map_arr = idx_map.unchecked<1>(); + + // Cache frequently accessed values + const size_t last_idx = E_neq_arr.shape(0) - 1; + const double first_e = E_neq_arr(0); + const double last_e = E_neq_arr(last_idx); + + // Quick boundary checks with cached values + if (E_target <= first_e) { + return T_eq_arr(0); + } + + if (E_target >= last_e) { + return T_eq_arr(last_idx); + } + + // Precompute and cache interpolation parameters + const double E_eq_start = E_eq_arr(0); + const double delta_E = E_eq_arr(1) - E_eq_start; + const double inv_delta_E = 1.0 / delta_E; + + const int idx_E_eq = std::min( + static_cast<int>((E_target - E_eq_start) * inv_delta_E), + static_cast<int>(idx_map_arr.shape(0) - 1) + ); + //std::cout << "idx_E_eq: " << idx_E_eq << std::endl; + + int idx_E_neq = idx_map_arr(idx_E_eq); + + if (E_neq_arr(idx_E_neq + 1) < E_target) { + ++idx_E_neq; + } + + // Get interpolation index + const double E1 = E_neq_arr(idx_E_neq); + const double E2 = E_neq_arr(idx_E_neq + 1); + const double T1 = T_eq_arr(idx_E_neq); + const double T2 = T_eq_arr(idx_E_neq + 1); + + // Optimized linear interpolation + const double inv_dE = 1.0 / (E2 - E1); + + return T1 + (T2 - T1) * (E_target - E1) * inv_dE; +} + +/* +double interpolate_double_lookup( + double E_target, + const py::array_t<double>& T_eq, + const py::array_t<double>& E_neq, + const py::array_t<double>& E_eq, + const py::array_t<int32_t>& idx_map) { + + // Cache all buffer requests at start + const auto T_eq_buf = T_eq.request(); + const auto E_neq_buf = E_neq.request(); + const auto E_eq_buf = E_eq.request(); + const auto idx_map_buf = idx_map.request(); + + // Get all pointers + const double* __restrict__ T_eq_ptr = static_cast<double*>(T_eq_buf.ptr); + const double* __restrict__ E_neq_ptr = static_cast<double*>(E_neq_buf.ptr); + const double* __restrict__ E_eq_ptr = static_cast<double*>(E_eq_buf.ptr); + const int32_t* __restrict__ idx_map_ptr = static_cast<int32_t*>(idx_map_buf.ptr); + std::cout << "C++ - First few T_eq_ptr values: "; + for(size_t i = 0; i < 5; i++) std::cout << T_eq_ptr[i] << " "; + std::cout << std::endl; + std::cout << "C++ - First few E_neq_ptr values: "; + for(size_t i = 0; i < 5; i++) std::cout << E_neq_ptr[i] << " "; + std::cout << std::endl; + + // Quick boundary checks + if (E_target <= E_neq_ptr[0]) { + return T_eq_ptr[0]; + } + + const size_t last_idx = E_neq_buf.size - 1; + if (E_target >= E_neq_ptr[last_idx]) { + return T_eq_ptr[last_idx]; + } + + // Cache delta_E to avoid recomputation + const double delta_E = E_eq_ptr[1] - E_eq_ptr[0]; + const double E_eq_start = E_eq_ptr[0]; + + // Compute indices + const int idx_E_eq = std::min( + static_cast<int>((E_target - E_eq_start) / delta_E), + static_cast<int>(idx_map_buf.size - 1) + ); + + int idx_E_neq = idx_map_ptr[idx_E_eq]; + if (E_neq_ptr[idx_E_neq + 1] < E_target) { + ++idx_E_neq; + } + + // Linear interpolation with minimal temporaries + const double E1 = E_neq_ptr[idx_E_neq]; + const double dE = E_neq_ptr[idx_E_neq + 1] - E1; + const double T1 = T_eq_ptr[idx_E_neq]; + + return T1 + (T_eq_ptr[idx_E_neq + 1] - T1) * (E_target - E1) / dE; +} +*/ +/* +double interpolate_double_lookup( + double E_target, + const py::array_t<double>& T_eq, + const py::array_t<double>& E_neq, + const py::array_t<double>& E_eq, + const py::array_t<int32_t>& idx_map) { + + auto T_eq_buf = T_eq.request(); + auto E_neq_buf = E_neq.request(); + auto E_eq_buf = E_eq.request(); + auto idx_map_buf = idx_map.request(); + + const double* T_eq_ptr = static_cast<double*>(T_eq_buf.ptr); + const double* E_neq_ptr = static_cast<double*>(E_neq_buf.ptr); + const double* E_eq_ptr = static_cast<double*>(E_eq_buf.ptr); + const int32_t* idx_map_ptr = static_cast<int32_t*>(idx_map_buf.ptr); + + if (E_target <= E_neq_ptr[0]) { + return T_eq_ptr[0]; + } + if (E_target >= E_neq_ptr[E_neq_buf.size - 1]) { + return T_eq_ptr[T_eq_buf.size - 1]; + } + + int idx_E_eq = static_cast<int>((E_target - E_eq_ptr[0]) / (E_eq_ptr[1] - E_eq_ptr[0])); + idx_E_eq = std::min(idx_E_eq, static_cast<int>(idx_map_buf.size - 1)); + + int idx_E_neq = idx_map_ptr[idx_E_eq]; + + if (E_neq_ptr[idx_E_neq + 1] < E_target) { + ++idx_E_neq; + } + + double E1 = E_neq_ptr[idx_E_neq]; + double E2 = E_neq_ptr[idx_E_neq + 1]; + double T1 = T_eq_ptr[idx_E_neq]; + double T2 = T_eq_ptr[idx_E_neq + 1]; + + return T1 + (T2 - T1) * (E_target - E1) / (E2 - E1); +} +*/ +/* +double interpolate_double_lookup( + double E_target, + const py::array_t<double>& T_eq, + const py::array_t<double>& E_neq, + const py::array_t<double>& E_eq, + const py::array_t<int32_t>& idx_map) { + + // Get buffer info once and cache sizes + const auto& E_neq_buf = E_neq.request(); + const auto E_neq_size = E_neq_buf.size; + const double* const E_neq_ptr = static_cast<double*>(E_neq_buf.ptr); + + // Quick boundary checks + if (E_target <= E_neq_ptr[0]) { + return static_cast<double*>(T_eq.request().ptr)[0]; + } + if (E_target >= E_neq_ptr[E_neq_size - 1]) { + return static_cast<double*>(T_eq.request().ptr)[E_neq_size - 1]; + } + + // Cache frequently used values + const double* const E_eq_ptr = static_cast<double*>(E_eq.request().ptr); + const double delta_E = E_eq_ptr[1] - E_eq_ptr[0]; + + // Calculate index in equidistant grid + const int idx_E_eq = std::min( + static_cast<int>((E_target - E_eq_ptr[0]) / delta_E), + static_cast<int>(idx_map.request().size - 1) + ); + + // Get mapped index and data pointers + int idx_E_neq = static_cast<int32_t*>(idx_map.request().ptr)[idx_E_eq]; + const double* const T_eq_ptr = static_cast<double*>(T_eq.request().ptr); + + // Adjust index if needed + if (E_neq_ptr[idx_E_neq + 1] < E_target) { + ++idx_E_neq; + } + + // Linear interpolation using direct array access + const double E1 = E_neq_ptr[idx_E_neq]; + const double E2 = E_neq_ptr[idx_E_neq + 1]; + const double T1 = T_eq_ptr[idx_E_neq]; + const double T2 = T_eq_ptr[idx_E_neq + 1]; + + return T1 + ((T2 - T1) * (E_target - E1)) / (E2 - E1); +} +*/ + + diff --git a/src/pymatlib/core/cpp/temperature_from_energy_density_array.h b/src/pymatlib/core/cpp/include/binary_search_interpolation.h similarity index 89% rename from src/pymatlib/core/cpp/temperature_from_energy_density_array.h rename to src/pymatlib/core/cpp/include/binary_search_interpolation.h index 0b31321..2b4dc4c 100644 --- a/src/pymatlib/core/cpp/temperature_from_energy_density_array.h +++ b/src/pymatlib/core/cpp/include/binary_search_interpolation.h @@ -1,6 +1,7 @@ #pragma once #include <pybind11/pybind11.h> #include <pybind11/numpy.h> +#include "interpolate_binary_search_cpp.h" namespace py = pybind11; @@ -13,7 +14,7 @@ namespace py = pybind11; * @return Interpolated temperature value * @throws std::runtime_error if arrays have different lengths or wrong monotonicity */ -double temperature_from_energy_density_array( +double interpolate_binary_search( const py::array_t<double>& temperature_array, double h_in, const py::array_t<double>& energy_density_array); diff --git a/src/pymatlib/core/cpp/include/double_lookup_interpolation.h b/src/pymatlib/core/cpp/include/double_lookup_interpolation.h new file mode 100644 index 0000000..936097b --- /dev/null +++ b/src/pymatlib/core/cpp/include/double_lookup_interpolation.h @@ -0,0 +1,17 @@ +#ifndef DOUBLE_LOOKUP_INTERPOLATION_H +#define DOUBLE_LOOKUP_INTERPOLATION_H + +#include <pybind11/pybind11.h> +#include <pybind11/numpy.h> +#include "interpolate_double_lookup_cpp.h" + +namespace py = pybind11; + +double interpolate_double_lookup( + double E_target, + const py::array_t<double>& T_eq, + const py::array_t<double>& E_neq, + const py::array_t<double>& E_eq, + const py::array_t<int32_t>& idx_map); + +#endif // DOUBLE_LOOKUP_INTERPOLATION_H diff --git a/src/pymatlib/core/cpp/include/interpolate_binary_search_cpp.h b/src/pymatlib/core/cpp/include/interpolate_binary_search_cpp.h new file mode 100644 index 0000000..ac85c4e --- /dev/null +++ b/src/pymatlib/core/cpp/include/interpolate_binary_search_cpp.h @@ -0,0 +1,83 @@ +#pragma once +#include <cmath> +#include <vector> +#include <stdexcept> + +/** + * Fast temperature interpolation using C++ containers + * + * @param temperature_array Container of temperature values + * @param h_in Energy density value to interpolate + * @param energy_density_array Container of energy density values + * @return Interpolated temperature value + * @throws std::runtime_error if arrays have different lengths or wrong monotonicity + */ + + +template<typename Container> +double interpolate_binary_search_cpp( + const Container& temperature_array, + double h_in, + const Container& energy_density_array) { + + static constexpr double EPSILON = 1e-6; + + // Input validation + const size_t n = temperature_array.size(); + if (n != energy_density_array.size() || n < 2) { + throw std::runtime_error("Invalid array sizes"); + } + + // Determine array order + const bool is_ascending = temperature_array[0] < temperature_array[n-1]; + const size_t start_idx = is_ascending ? 0 : n-1; + const size_t end_idx = is_ascending ? n-1 : 0; + + // Validate energy density increases with temperature + if (energy_density_array[start_idx] >= energy_density_array[end_idx]) { + throw std::runtime_error("Energy density must increase with temperature"); + } + + // Quick boundary checks + if (h_in <= energy_density_array[start_idx]) return temperature_array[start_idx]; + if (h_in >= energy_density_array[end_idx]) return temperature_array[end_idx]; + + // Binary search + size_t left = 0; + size_t right = n - 1; + + while (left <= right) { + const size_t mid = (left + right) / 2; + const double mid_val = energy_density_array[mid]; + + if (std::abs(mid_val - h_in) < EPSILON) { + return temperature_array[mid]; + } + + if (mid_val > h_in) { + if (is_ascending) { + right = mid - 1; + } else { + left = mid + 1; + } + } else { + if (is_ascending) { + left = mid + 1; + } else { + right = mid - 1; + } + } + } + + // Linear interpolation + const size_t idx1 = is_ascending ? right : left; + const size_t idx2 = is_ascending ? left : right; + + const double x0 = energy_density_array[idx1]; + const double x1 = energy_density_array[idx2]; + const double y0 = temperature_array[idx1]; + const double y1 = temperature_array[idx2]; + + const double slope = (y1 - y0) / (x1 - x0); + return y0 + slope * (h_in - x0); +} diff --git a/src/pymatlib/core/cpp/include/interpolate_double_lookup_cpp.h b/src/pymatlib/core/cpp/include/interpolate_double_lookup_cpp.h new file mode 100644 index 0000000..2889a4d --- /dev/null +++ b/src/pymatlib/core/cpp/include/interpolate_double_lookup_cpp.h @@ -0,0 +1,43 @@ +#pragma once +#include <cmath> +#include <vector> +#include <stdexcept> + + +template<typename Container> +double interpolate_double_lookup_cpp( + double E_target, + const Container& T_eq, + const Container& E_neq, + const Container& E_eq, + const Container& idx_map) { + + // Handle boundary cases + if (E_target <= E_neq[0]) { + return T_eq[0]; + } + if (E_target >= E_neq.back()) { + return T_eq.back(); + } + + // Calculate index in equidistant grid + int idx_E_eq = static_cast<int>((E_target - E_eq[0]) / (E_eq[1] - E_eq[0])); + idx_E_eq = std::min(idx_E_eq, static_cast<int>(idx_map.size() - 1)); + + // Get index from mapping + int idx_E_neq = idx_map[idx_E_eq]; + + // Adjust index if needed + if (E_neq[idx_E_neq + 1] < E_target) { + ++idx_E_neq; + } + + // Get interpolation points + const double E1 = E_neq[idx_E_neq]; + const double E2 = E_neq[idx_E_neq + 1]; + const double T1 = T_eq[idx_E_neq]; + const double T2 = T_eq[idx_E_neq + 1]; + + // Linear interpolation + return T1 + (T2 - T1) * (E_target - E1) / (E2 - E1); +} diff --git a/src/pymatlib/core/cpp/module.cpp b/src/pymatlib/core/cpp/module.cpp new file mode 100644 index 0000000..3c80f29 --- /dev/null +++ b/src/pymatlib/core/cpp/module.cpp @@ -0,0 +1,26 @@ +#include <pybind11/pybind11.h> +#include <pybind11/numpy.h> +#include "include/binary_search_interpolation.h" +#include "include/double_lookup_interpolation.h" + +namespace py = pybind11; + +PYBIND11_MODULE(fast_interpolation, m) { + m.doc() = "Fast interpolation methods implementation"; + + m.def("interpolate_binary_search", + &interpolate_binary_search, + "Find temperature using binary search and linear interpolation", + py::arg("temperature_array"), + py::arg("h_in"), + py::arg("energy_density_array")); + + m.def("interpolate_double_lookup", + &interpolate_double_lookup, + "Fast interpolation using double lookup method and linear interpolation", + py::arg("E_target"), + py::arg("T_eq"), + py::arg("E_neq"), + py::arg("E_eq"), + py::arg("idx_map")); +} diff --git a/src/pymatlib/core/cpp/temperature_from_energy_density_array.cpp b/src/pymatlib/core/cpp/temperature_from_energy_density_array.cpp deleted file mode 100644 index 7133bf7..0000000 --- a/src/pymatlib/core/cpp/temperature_from_energy_density_array.cpp +++ /dev/null @@ -1,89 +0,0 @@ -#include "temperature_from_energy_density_array.h" -#include <cmath> - - -double temperature_from_energy_density_array( - const py::array_t<double>& temperature_array, - double h_in, - const py::array_t<double>& energy_density_array) { - - static constexpr double EPSILON = 1e-6; - - const auto temp_buf = temperature_array.request(); - const auto energy_buf = energy_density_array.request(); - - if (temp_buf.size != energy_buf.size) { - throw std::runtime_error("Arrays must have same length"); - } - - const auto temp_ptr = static_cast<double*>(temp_buf.ptr); - const auto energy_ptr = static_cast<double*>(energy_buf.ptr); - const size_t n = temp_buf.size; - - if (n < 2) { - throw std::runtime_error("Arrays must have at least 2 elements"); - } - - // Determine array order - const bool is_ascending = temp_ptr[0] < temp_ptr[n-1]; - const size_t start_idx = is_ascending ? 0 : n-1; - const size_t end_idx = is_ascending ? n-1 : 0; - - // Validate energy density increases with temperature - /*if (energy_ptr[start_idx] >= energy_ptr[end_idx]) { - throw std::runtime_error("Energy density must increase with temperature"); - }*/ - - // Quick boundary checks - if (h_in <= energy_ptr[start_idx]) return temp_ptr[start_idx]; - if (h_in >= energy_ptr[end_idx]) return temp_ptr[end_idx]; - - // Binary search - size_t left = 0; - size_t right = n - 1; - - while (left <= right) { - const size_t mid = (left + right) / 2; - const double mid_val = energy_ptr[mid]; - - if (std::abs(mid_val - h_in) < EPSILON) { - return temp_ptr[mid]; - } - - if (mid_val > h_in) { - if (is_ascending) { - right = mid - 1; - } else { - left = mid + 1; - } - } else { - if (is_ascending) { - left = mid + 1; - } else { - right = mid - 1; - } - } - } - - // Linear interpolation - const size_t idx1 = is_ascending ? right : left; - const size_t idx2 = is_ascending ? left : right; - - const double x0 = energy_ptr[idx1]; - const double x1 = energy_ptr[idx2]; - const double y0 = temp_ptr[idx1]; - const double y1 = temp_ptr[idx2]; - - const double slope = (y1 - y0) / (x1 - x0); - return y0 + slope * (h_in - x0); -} - -PYBIND11_MODULE(fast_interpolation, m) { - m.doc() = "Fast temperature interpolation using binary search"; - m.def("temperature_from_energy_density_array", - &temperature_from_energy_density_array, - "Find temperature using binary search and linear interpolation", - py::arg("temperature_array"), - py::arg("h_in"), - py::arg("energy_density_array")); -} diff --git a/src/pymatlib/core/cpp/test_interpolation.cpp b/src/pymatlib/core/cpp/test_interpolation.cpp new file mode 100644 index 0000000..fa330f7 --- /dev/null +++ b/src/pymatlib/core/cpp/test_interpolation.cpp @@ -0,0 +1,57 @@ +#include "interpolate_binary_search_cpp.h" +#include <vector> +#include <array> +#include <iostream> +#include <cassert> + +void run_cpp_tests() { + // Test with std::vector + std::vector<double> temperatures = {3273.15, 3263.15, 3253.15, 3243.15}; + std::vector<double> energy_densities = {1.71e10, 1.70e10, 1.69e10, 1.68e10}; + + + try { + // Test normal interpolation + double h_in = 1.695e10; + double result = interpolate_binary_search_cpp( + temperatures, h_in, energy_densities); + std::cout << "Vector test temperature: " << result << std::endl; + assert(result > 3243.15 && result < 3273.15); + + // Test with std::array + std::array<double, 4> temp_arr = {3273.15, 3263.15, 3253.15, 3243.15}; + std::array<double, 4> energy_arr = {1.71e10, 1.70e10, 1.69e10, 1.68e10}; + result = interpolate_binary_search_cpp(temp_arr, h_in, energy_arr); + std::cout << "Array test temperature: " << result << std::endl; + assert(result > 3243.15 && result < 3273.15); + + // Test boundary values + result = interpolate_binary_search_cpp( + temperatures, 1.67e10, energy_densities); + assert(std::abs(result - 3243.15) < 1e-6); + + result = interpolate_binary_search_cpp( + temperatures, 1.72e10, energy_densities); + assert(std::abs(result - 3273.15) < 1e-6); + + // Test error cases + std::vector<double> wrong_size = {1.0, 2.0}; + try { + interpolate_binary_search_cpp(temperatures, h_in, wrong_size); + std::cerr << "Expected size mismatch error" << std::endl; + assert(false); + } catch (const std::runtime_error&) { + std::cout << "Size mismatch test passed" << std::endl; + } + + } catch (const std::exception& e) { + std::cerr << "Test error: " << e.what() << std::endl; + assert(false); + } +} + +int main() { + run_cpp_tests(); + std::cout << "All C++ tests completed successfully" << std::endl; + return 0; +} diff --git a/src/pymatlib/core/cpp/test_interpolation.py b/src/pymatlib/core/cpp/test_interpolation.py new file mode 100644 index 0000000..17df127 --- /dev/null +++ b/src/pymatlib/core/cpp/test_interpolation.py @@ -0,0 +1,30 @@ +import numpy as np +import pytest +from pymatlib.core.cpp.fast_interpolation import interpolate_binary_search + +def test_numpy_implementation(): + # Test data + temperatures = np.array([3273.15, 3263.15, 3253.15, 3243.15]) + energy_densities = np.array([1.71e10, 1.70e10, 1.69e10, 1.68e10]) + h_in = 1.695e10 + + # Test normal interpolation + result = interpolate_binary_search(temperatures, h_in, energy_densities) + assert 3243.15 < result < 3273.15 + print(f"NumPy test temperature: {result}") + + # Test boundary values + result_min = interpolate_binary_search(temperatures, 1.67e10, energy_densities) + assert abs(result_min - 3243.15) < 1e-6 + + result_max = interpolate_binary_search(temperatures, 1.72e10, energy_densities) + assert abs(result_max - 3273.15) < 1e-6 + + # Test error cases + with pytest.raises(RuntimeError): + wrong_size = np.array([1.0, 2.0]) + interpolate_binary_search(temperatures, h_in, wrong_size) + +if __name__ == "__main__": + test_numpy_implementation() + print("All Python tests completed successfully") diff --git a/src/pymatlib/core/data_handler.py b/src/pymatlib/core/data_handler.py index bda91b6..9e64063 100644 --- a/src/pymatlib/core/data_handler.py +++ b/src/pymatlib/core/data_handler.py @@ -1,6 +1,8 @@ +import os import numpy as np from pathlib import Path from typing import Union, Tuple +from matplotlib import pyplot as plt def print_results(file_path: str, temperatures: np.ndarray, material_property: np.ndarray) -> None: @@ -100,6 +102,58 @@ def thousand_times(q: Union[float, np.ndarray]) -> Union[float, np.ndarray]: return q * 1000 +# Moved from interpolators.py to data_handler.py +def check_equidistant(temp: np.ndarray) -> float: + """ + Tests if the temperature values are equidistant. + + :param temp: Array of temperature values. + :return: The common difference if equidistant, otherwise 0. + """ + if len(temp) < 2: + raise ValueError(f"{temp} array has length < 2") + + temperature_diffs = np.diff(temp) + if np.allclose(temperature_diffs, temperature_diffs[0], atol=1e-10): + return float(temperature_diffs[0]) + return 0.0 + + +def check_strictly_increasing(arr, name="Array", threshold=0.1): + """ + Check if array is strictly monotonically increasing and raise ValueError if not. + + Args: + arr: numpy array to check + name: name of array for reporting + threshold: minimum required difference between consecutive elements + + Raises: + ValueError: If array is not strictly increasing, with detailed information + about where the violation occurs + """ + for i in range(1, len(arr)): + diff = arr[i] - arr[i-1] + if diff <= threshold: + # Prepare error message with context + start_idx = max(0, i-2) + end_idx = min(len(arr), i+3) + context = "\nSurrounding values:\n" + for j in range(start_idx, end_idx): + context += f"Index {j}: {arr[j]:.10e}\n" + + error_msg = ( + f"{name} is not strictly increasing at index {i}:\n" + f"Previous value ({i-1}): {arr[i-1]:.10e}\n" + f"Current value ({i}): {arr[i]:.10e}\n" + f"Difference: {diff:.10e}\n" + f"{context}" + ) + raise ValueError(error_msg) + print(f"{name} is strictly monotonically increasing") + return True + + def find_min_max_temperature(temperatures_input) -> tuple: """ Find the minimum and maximum temperature from either a text file or a NumPy array. @@ -145,6 +199,31 @@ def find_min_max_temperature(temperatures_input) -> tuple: except Exception as e: raise ValueError(f"An error occurred while processing the input: {e}") + +def plot_arrays(x_arr: np.ndarray, y_arr: np.ndarray, x_label: str = None, y_label: str = None) -> None: + # Set labels and titles + x_label = x_label or "x-axis" + y_label = y_label or "y-axis" + + # Create the plot + plt.figure(figsize=(10, 6)) + plt.plot(x_arr, y_arr, 'b-', linewidth=1) + plt.xlabel(x_label) + plt.ylabel(y_label) + plt.title(f'{y_label} vs {x_label}') + plt.grid(True) + + # Define filename and directory + filename = f"{y_label.replace('/', '_')}_vs_{x_label.replace('/', '_')}.png" + directory = "plots" + os.makedirs(directory, exist_ok=True) # Ensure the directory exists + + filepath = os.path.join(directory, filename) + plt.savefig(filepath, dpi=300, bbox_inches="tight") + plt.show() + print(f"Plot saved as {filepath}") + + if __name__ == '__main__': # Example usage: # 1. Using a file path diff --git a/src/pymatlib/core/interpolators.py b/src/pymatlib/core/interpolators.py index 6354a88..33f259e 100644 --- a/src/pymatlib/core/interpolators.py +++ b/src/pymatlib/core/interpolators.py @@ -3,6 +3,7 @@ import sympy as sp from typing import Union, List, Tuple from pymatlib.core.models import wrapper, material_property_wrapper from pymatlib.core.typedefs import Assignment, ArrayTypes, MaterialProperty +from pymatlib.core.data_handler import check_equidistant COUNT = 0 @@ -129,22 +130,6 @@ def interpolate_lookup( raise ValueError(f"Invalid input for T: {type(T)}") -def check_equidistant(temp: np.ndarray) -> float: - """ - Tests if the temperature values are equidistant. - - :param temp: Array of temperature values. - :return: The common difference if equidistant, otherwise 0. - """ - if len(temp) <= 1: - return 0.0 - - temperature_diffs = np.diff(temp) - if np.allclose(temperature_diffs, temperature_diffs[0], atol=1e-10): - return float(temperature_diffs[0]) - return 0 - - def interpolate_property( T: Union[float, sp.Symbol], temp_array: ArrayTypes, @@ -192,5 +177,171 @@ def interpolate_property( if force_lookup or incr == 0 or len(temp_array) < temp_array_limit: return interpolate_lookup(T, temp_array, prop_array) else: - return interpolate_equidistant( - T, float(temp_array[0]), incr, prop_array) + return interpolate_equidistant(T, float(temp_array[0]), incr, prop_array) + + +# Moved from models.py to interpolators.py +def temperature_from_energy_density( + T: sp.Expr, + temperature_array: np.ndarray, + h_in: float, + energy_density: MaterialProperty) -> float: + """ + Compute the temperature for energy density using linear interpolation. + Args: + T: Symbol for temperature in material property. + temperature_array: Array of temperature values. + h_in: Target property value. + energy_density: Material property for energy_density. + Returns: + Corresponding temperature(s) for the input energy density value(s). + """ + T_min, T_max = np.min(temperature_array), np.max(temperature_array) + h_min, h_max = energy_density.evalf(T, T_min), energy_density.evalf(T, T_max) + + if h_in < h_min or h_in > h_max: + raise ValueError(f"The input energy density value of {h_in} is outside the computed range {h_min, h_max}.") + + tolerance: float = 5e-6 + max_iterations: int = 5000 + + iteration_count = 0 + + for _ in range(max_iterations): + iteration_count += 1 + # Linear interpolation to find T_1 + T_1 = T_min + (h_in - h_min) * (T_max - T_min) / (h_max - h_min) + # Evaluate h_1 at T_1 + h_1 = energy_density.evalf(T, T_1) + + if abs(h_1 - h_in) < tolerance: + # print(f"Linear interpolation converged in {iteration_count} iterations.") + return T_1 + + if h_1 < h_in: + T_min, h_min = T_1, h_1 + else: + T_max, h_max = T_1, h_1 + + raise RuntimeError(f"Linear interpolation did not converge within {max_iterations} iterations.") + + +# Moved from models.py to interpolators.py +def interpolate_binary_search( + temperature_array: np.ndarray, + h_in: float, + energy_density_array: np.ndarray, + epsilon: float = 1e-6) -> float: + + # Input validation + if len(temperature_array) != len(energy_density_array): + raise ValueError("temperature_array and energy_density_array must have the same length.") + + n = len(temperature_array) + is_ascending = temperature_array[0] < temperature_array[-1] + + # Critical boundary indices + start_idx = 0 if is_ascending else n - 1 + end_idx = n - 1 if is_ascending else 0 + + # Boundary checks + if h_in <= energy_density_array[start_idx]: + return float(temperature_array[start_idx]) + if h_in >= energy_density_array[end_idx]: + return float(temperature_array[end_idx]) + + # Binary search + left, right = 0, n - 1 + while left <= right: + mid = (left + right) // 2 + mid_val = energy_density_array[mid] + + if abs(mid_val - h_in) < epsilon: + return float(temperature_array[mid]) + + if (mid_val > h_in) == is_ascending: + right = mid - 1 + else: + left = mid + 1 + + # Linear interpolation + x0 = energy_density_array[right] + x1 = energy_density_array[left] + y0 = temperature_array[right] + y1 = temperature_array[left] + + return float(y0 + (y1 - y0) * (h_in - x0) / (x1 - x0)) + + +def E_eq_from_E_neq(E_neq: np.ndarray) -> np.ndarray: + # delta_E_neq = np.diff(E_neq) + delta_min: float = np.min(np.diff(E_neq)) + if delta_min < 1.: + raise ValueError(f"Energy density array points are very closely spaced, delta = {delta_min}") + delta_E_eq = max(np.floor(delta_min * 0.95), 1.) + E_eq = np.arange(E_neq[0], E_neq[-1] + delta_E_eq, delta_E_eq, dtype=np.float64) + + return E_eq + + +def create_idx_mapping(E_neq: np.ndarray, E_eq: np.ndarray) -> np.ndarray: + """idx_map = np.zeros(len(E_eq), dtype=int) + for i, e in enumerate(E_eq): + idx = int(np.searchsorted(E_neq, e) - 1) + idx = max(0, min(idx, len(E_neq) - 2)) # Bound check + idx_map[i] = idx""" + # idx_map = np.searchsorted(E_neq, E_eq) - 1 + idx_map = np.searchsorted(E_neq, E_eq, side='right') - 1 + idx_map = np.clip(idx_map, 0, len(E_neq) - 2) + + return idx_map.astype(np.int32) + + +def prepare_interpolation_arrays(temperature_array: np.ndarray, energy_density_array: np.ndarray)\ + -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + + # Input validation + if len(temperature_array) != len(energy_density_array): + raise ValueError("temperature_array and energy_density_array must have the same length.") + T_incr = check_equidistant(temperature_array) + + if T_incr == 0.0: + raise ValueError("Temperature array must be equidistant") + + # Convert to numpy arrays if not already + T_eq = np.asarray(temperature_array) + E_neq = np.asarray(energy_density_array) + + # Flip arrays if temperature increment is negative + if T_incr < 0.0: + T_eq = np.flip(T_eq) + E_neq = np.flip(E_neq) + + if E_neq[0] >= E_neq[-1]: + raise ValueError("Energy density must increase with temperature") + + # Create equidistant energy array and index mapping + E_eq = E_eq_from_E_neq(E_neq) + idx_mapping = create_idx_mapping(E_neq, E_eq) + + return T_eq, E_neq, E_eq, idx_mapping + + +def interpolate_double_lookup(E_target, T_eq, E_neq, E_eq, idx_map) -> float: + + if E_target <= E_neq[0]: + return T_eq[0] + if E_target >= E_neq[-1]: + return T_eq[-1] + + idx_E_eq = int((E_target-E_eq[0]) / (E_eq[1] - E_eq[0])) + idx_E_eq = min(idx_E_eq, len(idx_map) - 1) + + idx_E_neq = idx_map[idx_E_eq] + if E_neq[idx_E_neq + 1] < E_target: + idx_E_neq += 1 + + E1, E2 = E_neq[idx_E_neq], E_neq[idx_E_neq + 1] + T1, T2 = T_eq[idx_E_neq], T_eq[idx_E_neq + 1] + + return T1 + (T2 - T1) * (E_target - E1) / (E2 - E1) diff --git a/src/pymatlib/core/models.py b/src/pymatlib/core/models.py index e7cd9ce..406e0ad 100644 --- a/src/pymatlib/core/models.py +++ b/src/pymatlib/core/models.py @@ -166,87 +166,3 @@ def energy_density( _energy_density = density_expr * (temperature * heat_capacity_expr + latent_heat_expr) return MaterialProperty(_energy_density, sub_assignments) - - -def temperature_from_energy_density( - T: sp.Expr, - temperature_array: np.ndarray, - h_in: float, - energy_density: MaterialProperty - ) -> float: - """ - Compute the temperature for energy density using linear interpolation. - Args: - T: Symbol for temperature in material property. - temperature_array: Array of temperature values. - h_in: Target property value. - energy_density: Material property for energy_density. - Returns: - Corresponding temperature(s) for the input energy density value(s). - """ - T_min, T_max = np.min(temperature_array), np.max(temperature_array) - h_min, h_max = energy_density.evalf(T, T_min), energy_density.evalf(T, T_max) - - if h_in < h_min or h_in > h_max: - raise ValueError(f"The input energy density value of {h_in} is outside the computed range {h_min, h_max}.") - - tolerance: float = 5e-6 - max_iterations: int = 5000 - - iteration_count = 0 - - for _ in range(max_iterations): - iteration_count += 1 - # Linear interpolation to find T_1 - T_1 = T_min + (h_in - h_min) * (T_max - T_min) / (h_max - h_min) - # Evaluate h_1 at T_1 - h_1 = energy_density.evalf(T, T_1) - - if abs(h_1 - h_in) < tolerance: - return T_1 - - if h_1 < h_in: - T_min, h_min = T_1, h_1 - else: - T_max, h_max = T_1, h_1 - - raise RuntimeError(f"Linear interpolation did not converge within {max_iterations} iterations.") - - -def temperature_from_energy_density_array( - temperature_array: np.ndarray, - h_in: float, - energy_density_array: np.ndarray - ) -> float: - # Input validation - if len(temperature_array) != len(energy_density_array): - raise ValueError("temperature_array and energy_density_array must have the same length.") - - # Binary search to find the interval - left, right = 0, len(energy_density_array) - 1 - - while left <= right: - mid = (left + right) // 2 - if energy_density_array[mid] == h_in: - return float(temperature_array[mid]) - elif energy_density_array[mid] > h_in: - left = mid + 1 - else: - right = mid - 1 - - # After binary search, 'right' points to the upper bound - # and 'left' points to the lower bound of our interval - if left >= len(energy_density_array): - return float(temperature_array[-1]) - if right < 0: - return float(temperature_array[0]) - - # Linear interpolation within the found interval - x0, x1 = energy_density_array[right], energy_density_array[left] - y0, y1 = temperature_array[right], temperature_array[left] - - # Use high-precision arithmetic for interpolation - slope = (y1 - y0) / (x1 - x0) - temperature = y0 + slope * (h_in - x0) - - return float(temperature) diff --git a/src/pymatlib/core/test_interpolate_functs_perf.py b/src/pymatlib/core/test_interpolate_functs_perf.py new file mode 100644 index 0000000..bd4b1f3 --- /dev/null +++ b/src/pymatlib/core/test_interpolate_functs_perf.py @@ -0,0 +1,222 @@ +import time +import numpy as np +import sympy as sp +from pathlib import Path +from typing import Union +from pymatlib.core.alloy import Alloy +from pymatlib.data.element_data import Fe, Cr, Ni, Mo, Mn +from pymatlib.core.models import thermal_diffusivity_by_heat_conductivity, density_by_thermal_expansion, energy_density +from pymatlib.core.data_handler import read_data, celsius_to_kelvin, thousand_times, check_equidistant, check_strictly_increasing, plot_arrays +from pymatlib.core.interpolators import (interpolate_property, temperature_from_energy_density, + E_eq_from_E_neq, create_idx_mapping, prepare_interpolation_arrays) +#from pymatlib.core.interpolators import interpolate_binary_search, interpolate_double_lookup +from pymatlib.core.cpp.fast_interpolation import interpolate_binary_search, interpolate_double_lookup + + +def generate_target_points(E_min: float, E_max: float, num_points: int) -> np.ndarray: + """Generate target energy points with specified distribution.""" + below_count = int(num_points * 0.225) # 22.5% below minimum + above_count = int(num_points * 0.225) # 22.5% above maximum + boundary_count = int(num_points * 0.05) # 5% boundary points + inside_count = num_points - (below_count + above_count + boundary_count) + + below_points = np.random.uniform(E_min - 5_000_000, E_min, size=below_count) + above_points = np.random.uniform(E_max, E_max + 5_000_000, size=above_count) + boundary_points = np.array([E_min, E_max] * (boundary_count // 2), dtype=np.float64) + inside_points = np.random.uniform(E_min, E_max, size=inside_count) + + points = np.concatenate([ + np.float64(inside_points), + np.float64(below_points), + np.float64(above_points), + boundary_points + ]) + np.random.shuffle(points) + return points + + +def compare_interpolation_methods(T_eq: np.ndarray, E_neq: np.ndarray, E_target: np.ndarray, label: str = "") -> None: + """Compare binary search and double lookup interpolation methods.""" + E_eq = E_eq_from_E_neq(E_neq) + idx_mapping = create_idx_mapping(E_neq, E_eq) + + # Time binary search method + start_time = time.perf_counter() + T_binary = [interpolate_binary_search(T_eq, float(E), E_neq) for E in E_target] + binary_time = time.perf_counter() - start_time + + # Time double lookup method + start_time = time.perf_counter() + T_double = [interpolate_double_lookup(float(E), T_eq, E_neq, E_eq, idx_mapping) + for E in E_target] + double_time = time.perf_counter() - start_time + + print(f"\nResults for {label}:") + print(f"Binary search time: {binary_time:.8f} s") + print(f"Double lookup time: {double_time:.8f} s") + + # Check for mismatches + mismatches = [(E, T1, T2) for E, T1, T2 in zip(E_target, T_binary, T_double) + if abs(T1 - T2) > 1e-8] + + if mismatches: + print("Mismatches found (E_target, T_binary, T_double):") + for E, T1, T2 in mismatches[:5]: # Show only first 5 mismatches + print(f"E={E:.8f}, T1={T1:.8f}, T2={T2:.8f}, diff={abs(T1-T2):.8f}") + else: + print("No mismatches found between methods") + + +def create_alloy(T: Union[float, sp.Symbol]) -> Alloy: + alloy = Alloy( + elements=[Fe, Cr, Ni, Mo, Mn], + composition=[0.675, 0.17, 0.12, 0.025, 0.01], # Fe: 67.5%, Cr: 17%, Ni: 12%, Mo: 2.5%, Mn: 1% + temperature_solidus=1653.15, # Solidus temperature in Kelvin (test at 1653.15 K = 1380 C) + temperature_liquidus=1723.15, # Liquidus temperature in Kelvin (test at 1723.15 K = 1450 C) + thermal_expansion_coefficient=16.3e-6 # in 1/K + ) + + # Determine the base directory + base_dir = Path(__file__).parent + + # Paths to data files using relative paths + density_data_file_path = str(base_dir/'..'/'data'/'alloys'/'SS316L'/'density_temperature_edited.txt') + heat_capacity_data_file_path = str(base_dir/'..'/'data'/'alloys'/'SS316L'/'heat_capacity_temperature_edited.txt') + + # Read temperature and material property data from the files + density_temp_array, density_array = read_data(density_data_file_path) + heat_capacity_temp_array, heat_capacity_array = read_data(heat_capacity_data_file_path) + + # Ensure the data was loaded correctly + if any(arr.size < 2 for arr in [density_temp_array, density_array, heat_capacity_temp_array, heat_capacity_array]): + raise ValueError("Failed to load temperature or material data from the file.") + + # Conversion: temperature to K and/or other material properties to SI units if necessary + density_temp_array = celsius_to_kelvin(density_temp_array) + density_array = thousand_times(density_array) # Density in kg/m³ # gm/cm³ -> kg/m³ + # plot_arrays(density_temp_array, density_array, "Temperature (K)", "Density (Kg/m^3)") + heat_capacity_temp_array = celsius_to_kelvin(heat_capacity_temp_array) + heat_capacity_array = thousand_times(heat_capacity_array) # Specific heat capacity in J/(kg·K) # J/g-K -> J/kg-K + # plot_arrays(heat_capacity_temp_array, heat_capacity_array, "Temperature (K)", "Heat Capacity (J/Kg-K)") + + alloy.density = interpolate_property(T, density_temp_array, density_array) + alloy.heat_capacity = interpolate_property(T, heat_capacity_temp_array, heat_capacity_array) + alloy.latent_heat_of_fusion = interpolate_property(T, alloy.solidification_interval(), np.array([0.0, 260000.0])) + alloy.energy_density = energy_density(T, alloy.density, alloy.heat_capacity, alloy.latent_heat_of_fusion) + + # Populate temperature_array and energy_density_array + alloy.temperature_array = density_temp_array + + alloy.energy_density_array = np.array([ + alloy.energy_density.evalf(T, temp) for temp in density_temp_array + ]) + # plot_arrays(alloy.temperature_array, alloy.energy_density_array, "Temperature (K)", "Energy Density (J/m^3)") + + T_eq, E_neq, E_eq, idx_map = prepare_interpolation_arrays( + alloy.temperature_array, + alloy.energy_density_array + ) + + check_strictly_increasing(T_eq, "T_eq") + check_strictly_increasing(E_neq, "E_neq") + + # Online execution + E = alloy.energy_density.evalf(T, alloy.temperature_liquidus) + + start_time1 = time.perf_counter() + T_interpolate1 = interpolate_binary_search(alloy.temperature_array, E, alloy.energy_density_array) + execution_time1 = time.perf_counter() - start_time1 + print(f"Interpolated temperature: {T_interpolate1}") + print(f"Execution time: {execution_time1:.8f} seconds") + + start_time2 = time.perf_counter() + T_interpolate2 = interpolate_binary_search(T_eq, E, E_neq) + execution_time2 = time.perf_counter() - start_time2 + print(f"Interpolated temperature: {T_interpolate2}") + print(f"Execution time: {execution_time2:.8f} seconds") + + start_time3 = time.perf_counter() + T_interpolate3 = interpolate_double_lookup(E, T_eq, E_neq, E_eq, idx_map) + execution_time3 = time.perf_counter() - start_time3 + print(f"Interpolated temperature: {T_interpolate3}") + print(f"Execution time: {execution_time3:.8f} seconds") + + if not (T_interpolate1 == T_interpolate2 == T_interpolate3): + raise ValueError(f"Mismatch value. Temperature value should be {alloy.temperature_liquidus}") + + + def measure_performance(iterations=1): + all_execution_times = np.zeros(iterations) + + for i in range(iterations): + start_measure_performance = time.perf_counter() + + results = [interpolate_binary_search(T_eq, E, E_neq) for E in E_neq] + # results = [interpolate_double_lookup(E, T_eq, E_neq, E_eq, idx_map) for E in E_neq] + + all_execution_times[i] = time.perf_counter() - start_measure_performance + + # Calculate statistics using numpy + total_time = np.sum(all_execution_times) + avg_time = np.mean(all_execution_times) + avg_per_iteration = avg_time / len(E_neq) + + print(f"Total execution time ({iterations} runs): {total_time:.8f} seconds") + print(f"Average execution time per run: {avg_time:.8f} seconds") + print(f"Average execution time per iteration: {avg_per_iteration:.8f} seconds") + + # Run the performance test + measure_performance(iterations=10_000) + + E_target = generate_target_points(float(E_neq[0]), float(E_neq[-1]), 100) + compare_interpolation_methods(T_eq, E_neq, E_target, "Alloy Dataset") + + return alloy + + +if __name__ == '__main__': + Temp = sp.Symbol('T') + alloy = create_alloy(Temp) + + # Test arrays + T_eq_small = np.array([100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.], dtype=np.float64) + E_neq_small = np.array([1000., 1220., 1650., 2020., 2260., 2609., 3050., 3623., 3960., 4210.], dtype=np.float64) + # 0 1 2 3 4 5 6 7 8 9 + + # Equidistant energy array with delta_E_eq = 200 (smaller than min_delta = 220) + E_eq_small = E_eq_from_E_neq(E_neq_small) + # E_eq = np.array([1000., 1209., 1418., 1627., 1836., 2045., 2254., 2463., 2672., + # 2881., 3090., 3299., 3508., 3717., 3926., 4135., 4344.]) + + # Index mapping array + idx_mapping = create_idx_mapping(E_neq_small, E_eq_small) + # idx_map = np.array([0, 0, 1, 1, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8]) + + # Test target energy values + E_target = 1.*np.array([1000, 1585, 2688, 3960, 4210]) + + for target in E_target: + T_star = interpolate_binary_search(T_eq_small, target, E_neq_small) + T_interpolate_double_lookup = interpolate_double_lookup(target, T_eq_small, E_neq_small, E_eq_small, idx_mapping) + if T_star != T_interpolate_double_lookup: + raise ValueError(f"Value Mismatch. {T_star} != {T_interpolate_double_lookup}") + + E_target_small = generate_target_points(float(E_neq_small[0]), float(E_neq_small[-1]), 100) + + compare_interpolation_methods(T_eq_small, E_neq_small, E_target_small, "Small Dataset") + + + # Create larger test arrays + size = 5_000_000 + T_eq_large = np.linspace(0.0, 5_000_000.0, size, dtype=np.float64) # Equidistant temperature values + # Generate non-equidistant energy density array + # Using cumsum of random values ensures monotonically increasing values + E_neq_large = np.cumsum(np.random.uniform(1, 10, size)) + 5_000_000.0 + + E_eq_large = E_eq_from_E_neq(E_neq_large) + + idx_mapping_large = create_idx_mapping(E_neq_large, E_eq_large) + + E_target_large = generate_target_points(float(E_neq_large[0]), float(E_neq_large[-1]), 1_000_000) + + compare_interpolation_methods(T_eq_small, E_neq_small, E_target_small, "Large Dataset") diff --git a/src/pymatlib/data/alloys/SS316L/SS316L.py b/src/pymatlib/data/alloys/SS316L/SS316L.py index 19b3ec8..8e7ee34 100644 --- a/src/pymatlib/data/alloys/SS316L/SS316L.py +++ b/src/pymatlib/data/alloys/SS316L/SS316L.py @@ -6,10 +6,10 @@ from typing import Union from matplotlib import pyplot as plt from pymatlib.core.alloy import Alloy from pymatlib.data.element_data import Fe, Cr, Ni, Mo, Mn -from pymatlib.core.models import thermal_diffusivity_by_heat_conductivity, density_by_thermal_expansion, energy_density, temperature_from_energy_density +from pymatlib.core.models import thermal_diffusivity_by_heat_conductivity, density_by_thermal_expansion, energy_density from pymatlib.core.data_handler import read_data, celsius_to_kelvin, thousand_times -from pymatlib.core.interpolators import interpolate_property -from pymatlib.core.cpp.fast_interpolation import temperature_from_energy_density_array +from pymatlib.core.interpolators import interpolate_property, prepare_interpolation_arrays#, interpolate_binary_search +from pymatlib.core.cpp.fast_interpolation import interpolate_binary_search, interpolate_double_lookup def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy: @@ -62,7 +62,7 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy: base_dir = Path(__file__).parent # Directory of the current file # Paths to data files using relative paths - density_data_file_path = str(base_dir / 'density_temperature.txt') + density_data_file_path = str(base_dir / 'density_temperature_edited.txt') heat_capacity_data_file_path = str(base_dir / 'heat_capacity_temperature_edited.txt') heat_conductivity_data_file_path = str(base_dir / '..' / 'SS316L' / 'heat_conductivity_temperature.txt') @@ -105,7 +105,7 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy: SS316L.energy_density_liquidus, SS316L.energy_density_array) - T_star = temperature_from_energy_density_array(*args) + T_star = interpolate_binary_search(*args) print(f"T_star: {T_star}") return SS316L diff --git a/src/pymatlib/data/alloys/SS316L/density_temperature_edited.txt b/src/pymatlib/data/alloys/SS316L/density_temperature_edited.txt new file mode 100644 index 0000000..ac3999b --- /dev/null +++ b/src/pymatlib/data/alloys/SS316L/density_temperature_edited.txt @@ -0,0 +1,301 @@ +T(C) 'D-1.0(C/s)' +3000.0 5.667564404 +2990.0 5.676952519 +2980.0 5.686342789 +2970.0 5.695735083 +2960.0 5.705129266 +2950.0 5.714525204 +2940.0 5.72392276 +2930.0 5.733321798 +2920.0 5.74272218 +2910.0 5.752123768 +2900.0 5.761526422 +2890.0 5.770930002 +2880.0 5.780334366 +2870.0 5.789739372 +2860.0 5.799144876 +2850.0 5.808550735 +2840.0 5.817956803 +2830.0 5.827362933 +2820.0 5.83676898 +2810.0 5.846174794 +2800.0 5.855580226 +2790.0 5.864985128 +2780.0 5.874389347 +2770.0 5.883792732 +2760.0 5.893195129 +2750.0 5.902596387 +2740.0 5.911996349 +2730.0 5.92139486 +2720.0 5.930791763 +2710.0 5.940186901 +2700.0 5.949580115 +2690.0 5.958971247 +2680.0 5.968360134 +2670.0 5.977746618 +2660.0 5.987130534 +2650.0 5.99651172 +2640.0 6.005890011 +2630.0 6.015265243 +2620.0 6.02463725 +2610.0 6.034005865 +2600.0 6.043370919 +2590.0 6.052732245 +2580.0 6.062089672 +2570.0 6.07144303 +2560.0 6.080792147 +2550.0 6.090136851 +2540.0 6.099476968 +2530.0 6.108812324 +2520.0 6.118142744 +2510.0 6.127468052 +2500.0 6.136788071 +2490.0 6.146102622 +2480.0 6.155411527 +2470.0 6.164714607 +2460.0 6.17401168 +2450.0 6.183302565 +2440.0 6.19258708 +2430.0 6.201865041 +2420.0 6.211136264 +2410.0 6.220400565 +2400.0 6.229657757 +2390.0 6.238907653 +2380.0 6.248150066 +2370.0 6.257384807 +2360.0 6.266611687 +2350.0 6.275830516 +2340.0 6.285041102 +2330.0 6.294243254 +2320.0 6.303436779 +2310.0 6.312621483 +2300.0 6.321797171 +2290.0 6.33096365 +2280.0 6.340120722 +2270.0 6.34926819 +2260.0 6.358405857 +2250.0 6.367533525 +2240.0 6.376650993 +2230.0 6.385758063 +2220.0 6.394854532 +2210.0 6.4039402 +2200.0 6.413014864 +2190.0 6.422078321 +2180.0 6.431130366 +2170.0 6.440170795 +2160.0 6.449199403 +2150.0 6.458215983 +2140.0 6.467220328 +2130.0 6.476212231 +2120.0 6.485191483 +2110.0 6.494157875 +2100.0 6.503111197 +2090.0 6.512051239 +2080.0 6.520977789 +2070.0 6.529890636 +2060.0 6.538789566 +2050.0 6.547674367 +2040.0 6.556544824 +2030.0 6.565400724 +2020.0 6.57424185 +2010.0 6.583067987 +2000.0 6.591878918 +1990.0 6.600674426 +1980.0 6.609454294 +1970.0 6.618218303 +1960.0 6.626966234 +1950.0 6.635697868 +1940.0 6.644412984 +1930.0 6.653111361 +1920.0 6.661792779 +1910.0 6.670457016 +1900.0 6.67910385 +1890.0 6.687733056 +1880.0 6.696344413 +1870.0 6.704937696 +1860.0 6.713512681 +1850.0 6.722069143 +1840.0 6.730606857 +1830.0 6.739125596 +1820.0 6.747625135 +1810.0 6.756105246 +1800.0 6.764565703 +1790.0 6.773006278 +1780.0 6.781426743 +1770.0 6.78982687 +1760.0 6.79820643 +1750.0 6.806565193 +1740.0 6.81490293 +1730.0 6.823219411 +1720.0 6.831514406 +1710.0 6.839787684 +1700.0 6.848039015 +1690.0 6.856268167 +1680.0 6.864474909 +1670.0 6.872659008 +1660.0 6.880820234 +1650.0 6.888958353 +1640.0 6.897073133 +1630.0 6.905164341 +1620.0 6.913231745 +1610.0 6.92127511 +1600.0 6.929294204 +1590.0 6.937288793 +1580.0 6.945258643 +1570.0 6.953203521 +1560.0 6.961123192 +1550.0 6.969017421 +1540.0 6.976885976 +1530.0 6.984728621 +1520.0 6.992545122 +1510.0 7.000335245 +1500.0 7.008098755 +1490.0 7.015835416 +1480.0 7.023544996 +1470.0 7.031227258 +#1460.0 7.038881969 +#heavily manipulated +1460.0 7.065976992 +1450.0 7.101066042 +1440.0 7.134387853 +1430.0 7.177458788 +1420.0 7.211193546 +1410.0 7.228319303 +1400.0 7.239681101 +1390.0 7.249455888 +1380.0 7.258134874 +1370.0 7.266781373 +1360.0 7.275378823 +1350.0 7.283907087 +1340.0 7.292372432 +1330.0 7.300806759 +1320.0 7.309120646 +1310.0 7.317361161 +1300.0 7.325511618 +1290.0 7.333644738 +1280.0 7.341683473 +1270.0 7.349597336 +1260.0 7.357432453 +1250.0 7.365195996 +1240.0 7.372887415 +1230.0 7.380522323 +1220.0 7.388029457 +1210.0 7.395402654 +1200.0 7.401625175 +1190.0 7.406659562 +1180.0 7.411685346 +1170.0 7.416702547 +1160.0 7.42171118 +1150.0 7.426711264 +1140.0 7.431702814 +1130.0 7.436685849 +1120.0 7.441660386 +1110.0 7.446626444 +1100.0 7.451584039 +1090.0 7.456533191 +1080.0 7.461473918 +1070.0 7.466406239 +1060.0 7.471330172 +1050.0 7.476245738 +1040.0 7.481152955 +1030.0 7.486051844 +1020.0 7.490942424 +1010.0 7.495824716 +1000.0 7.500698741 +990.0 7.505564519 +980.0 7.510422072 +970.0 7.515271421 +960.0 7.520112589 +950.0 7.524945596 +940.0 7.529770467 +930.0 7.534587223 +920.0 7.539395888 +910.0 7.544196486 +900.0 7.548989038 +890.0 7.553773571 +880.0 7.558550107 +870.0 7.563318673 +860.0 7.568079292 +850.0 7.57283199 +840.0 7.577576792 +830.0 7.582313725 +820.0 7.587042814 +810.0 7.591764087 +800.0 7.59647757 +790.0 7.601183289 +780.0 7.605881273 +770.0 7.61057155 +760.0 7.615254146 +750.0 7.619929092 +740.0 7.624596415 +730.0 7.629256145 +720.0 7.633908311 +710.0 7.638552943 +700.0 7.64319007 +690.0 7.647819723 +680.0 7.652441932 +670.0 7.657056729 +660.0 7.661664145 +650.0 7.666264211 +640.0 7.67085696 +630.0 7.675442422 +620.0 7.680020632 +610.0 7.684591621 +600.0 7.689155423 +590.0 7.693712072 +580.0 7.6982616 +570.0 7.702804042 +560.0 7.707339433 +550.0 7.711867807 +540.0 7.7163892 +530.0 7.720903645 +520.0 7.72541118 +510.0 7.72991184 +500.0 7.734405661 +490.0 7.73889268 +480.0 7.743372934 +470.0 7.74784646 +460.0 7.752313295 +450.0 7.756773477 +440.0 7.761227045 +430.0 7.765674036 +420.0 7.77011449 +410.0 7.774548445 +400.0 7.778975941 +390.0 7.783397018 +380.0 7.787811714 +370.0 7.792220072 +360.0 7.796622131 +350.0 7.801017932 +340.0 7.805407516 +330.0 7.809790925 +320.0 7.8141682 +310.0 7.818539384 +300.0 7.82290452 +290.0 7.827263649 +280.0 7.831616814 +270.0 7.83596406 +260.0 7.84030543 +250.0 7.844640967 +240.0 7.848970716 +230.0 7.853294721 +220.0 7.857613028 +210.0 7.861925681 +200.0 7.866232725 +190.0 7.870534208 +180.0 7.874830173 +170.0 7.879120669 +160.0 7.883405741 +150.0 7.887685437 +140.0 7.891959803 +130.0 7.896228887 +120.0 7.900492738 +110.0 7.904751403 +100.0 7.90900493 +90.0 7.913253369 +80.0 7.917496768 +70.0 7.921735177 +60.0 7.925968645 +50.0 7.930197222 +40.0 7.93442096 +30.0 7.938639907 -- GitLab