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

@@ -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"

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++'
     ),

src/pymatlib/core/cpp/CMakeLists.txt

 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"
+)

src/pymatlib/core/cpp/binary_search_interpolation.cpp

+#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);
+}

src/pymatlib/core/cpp/double_lookup_interpolation.cpp

+#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;
+}

src/pymatlib/core/cpp/temperature_from_energy_density_array.h → src/pymatlib/core/cpp/include/binary_search_interpolation.h

 #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);

src/pymatlib/core/cpp/include/double_lookup_interpolation.h

+#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

src/pymatlib/core/cpp/temperature_from_energy_density_array.cpp → src/pymatlib/core/cpp/include/interpolate_binary_search_cpp.h

-#include "temperature_from_energy_density_array.h"
+#pragma once
 #include <cmath>
-
-
-double temperature_from_energy_density_array(
-    const py::array_t<double>& temperature_array,
+#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 py::array_t<double>& energy_density_array) {
+    const Container& 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");
+    // 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 = temp_ptr[0] < temp_ptr[n-1];
+    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_ptr[start_idx] >= energy_ptr[end_idx]) {
+    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_ptr[start_idx]) return temp_ptr[start_idx];
-    if (h_in >= energy_ptr[end_idx]) return temp_ptr[end_idx];
+    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;
@@ -44,10 +48,10 @@ double temperature_from_energy_density_array(
 
     while (left <= right) {
         const size_t mid = (left + right) / 2;
-        const double mid_val = energy_ptr[mid];
+        const double mid_val = energy_density_array[mid];
 
         if (std::abs(mid_val - h_in) < EPSILON) {
-            return temp_ptr[mid];
+            return temperature_array[mid];
         }
 
         if (mid_val > h_in) {
@@ -69,21 +73,11 @@ double temperature_from_energy_density_array(
     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 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);
 }
-
-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"));
-}

src/pymatlib/core/cpp/include/interpolate_double_lookup_cpp.h

+#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);
+}

src/pymatlib/core/cpp/module.cpp

+#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"));
+}

src/pymatlib/core/cpp/test_interpolation.cpp

+#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;
+}

src/pymatlib/core/cpp/test_interpolation.py

+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")

src/pymatlib/core/data_handler.py

+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

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)

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)

src/pymatlib/core/test_interpolate_functs_perf.py

+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")

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

src/pymatlib/data/alloys/SS316L/density_temperature_edited.txt

+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