diff --git a/pyproject.toml b/pyproject.toml
index 76a46a4eb4afec8011c531d969d4c8c0845d89c8..ddb7a471426bf0dcf3dda14d8b7a243a08681a7d 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 ed646c5d656344e5ea6ea923100a467229a843c7..594f94fd72e39e0f3f954848372f5b2582c7d8be 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 333ab992a678faa09543bc0cc2793bbd42296b82..26261a13834dbf58386a6d01f5ef6ec414494048 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 0000000000000000000000000000000000000000..97379d9c5bd96801602770acb6944fc54b75c978
--- /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 0000000000000000000000000000000000000000..4f1673e4cd6303a1e981081941e65eac857b58a5
--- /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 0b31321c58c47ef3df497d24d0f73582de6ac81c..2b4dc4cbd5ca06e74f6e008ef7ad5e53daddb534 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 0000000000000000000000000000000000000000..936097b2c832dd3b09863ab9ef1e6a9d961704be
--- /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 0000000000000000000000000000000000000000..ac85c4ebb8de18155de0108e65b8cea4d7b3f477
--- /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 0000000000000000000000000000000000000000..2889a4dfc9b0addc898220a8af66b6e4d7bbfc56
--- /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 0000000000000000000000000000000000000000..3c80f2921a90a69ff2c4c906bcc5d3f75f23d2fb
--- /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 7133bf75198e7afa93457a1f92e539c33477d680..0000000000000000000000000000000000000000
--- 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 0000000000000000000000000000000000000000..fa330f7d89c4899cda470e5a6a17cda8b4badb93
--- /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 0000000000000000000000000000000000000000..17df127f0d586760ec400f37d3a6bf5071747567
--- /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 bda91b600e75be30ccb37039edad23c5c0a4bec7..9e6406324d9987f5aa21df9c671af0d7b8d27a3d 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 6354a88cfe29a56ae88053c129eb50ef2cb0f868..33f259e85fab2b27f9ffe26a33e3ac418e6079ff 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 e7cd9ce8b0b75c5b86a88165c0d3893ba54bf42f..406e0ad318656295aea6e18d13db8b49cfeb40a8 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 0000000000000000000000000000000000000000..bd4b1f355266e3b3880ff388c69c8cbeaf90c9f5
--- /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 19b3ec814ec9554f44612da93c447ec7ae3d1765..8e7ee3433934d9299ff6646bfdbaf2055738712a 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 0000000000000000000000000000000000000000..ac3999bc419ddb58cc72c1062156884067366946
--- /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