36#include <pybind11/cast.h>
37#include <pybind11/functional.h>
38#include <pybind11/numpy.h>
43 py::module math_module = root_module.def_submodule(
"math",
"Shamrock math lib");
45 shampylib::init_shamrock_math_AABB<f64_3>(math_module,
"AABB_f64_3");
46 shampylib::init_shamrock_math_Ray<f64_3>(math_module,
"Ray_f64_3");
47 shampylib::init_shamrock_math_RingRay<f64_3>(math_module,
"RingRay_f64_3");
48 shampylib::init_shamrock_math_sfc(math_module);
49 shampylib::init_shamrock_math_sphkernels(math_module);
51 math_module.def(
"derivative_upwind", [](
f64 x,
f64 dx, std::function<
f64(
f64)> fct) {
56 math_module.def(
"derivative_centered", [](
f64 x,
f64 dx, std::function<
f64(
f64)> fct) {
61 math_module.def(
"derivative_3point_forward", [](
f64 x,
f64 dx, std::function<
f64(
f64)> fct) {
66 math_module.def(
"derivative_3point_backward", [](
f64 x,
f64 dx, std::function<
f64(
f64)> fct) {
71 math_module.def(
"derivative_5point_midpoint", [](
f64 x,
f64 dx, std::function<
f64(
f64)> fct) {
82 Estim the correct step to use for a given order when using derivatives
85 py::class_<shammath::paving_function_periodic_3d<f64_3>>(
86 math_module, "paving_function_periodic_3d")
87 .def(py::init([](f64_3 box_size) {
88 return std::make_unique<shammath::paving_function_periodic_3d<f64_3>>(
93 .def(
"f_aabb", &shammath::paving_function_periodic_3d<f64_3>::f_aabb)
94 .def(
"f_aabb_inv", &shammath::paving_function_periodic_3d<f64_3>::f_aabb_inv)
96 "get_paving_index_intersecting",
97 &shammath::paving_function_periodic_3d<f64_3>::get_paving_index_intersecting);
99 py::class_<shammath::paving_function_general_3d<f64_3>>(
100 math_module,
"paving_function_general_3d")
102 py::init([](f64_3 box_size,
106 bool is_z_periodic) {
107 return std::make_unique<shammath::paving_function_general_3d<f64_3>>(
109 .box_size = box_size,
110 .box_center = box_center,
111 .is_x_periodic = is_x_periodic,
112 .is_y_periodic = is_y_periodic,
113 .is_z_periodic = is_z_periodic});
117 .def(
"f_aabb", &shammath::paving_function_general_3d<f64_3>::f_aabb)
118 .def(
"f_aabb_inv", &shammath::paving_function_general_3d<f64_3>::f_aabb_inv)
120 "get_paving_index_intersecting",
121 &shammath::paving_function_general_3d<f64_3>::get_paving_index_intersecting);
123 py::class_<shammath::paving_function_general_3d_shear_x<f64_3>>(
124 math_module,
"paving_function_general_3d_shear_x")
126 py::init([](f64_3 box_size,
132 return std::make_unique<shammath::paving_function_general_3d_shear_x<f64_3>>(
134 .box_size = box_size,
135 .box_center = box_center,
136 .is_x_periodic = is_x_periodic,
137 .is_y_periodic = is_y_periodic,
138 .is_z_periodic = is_z_periodic,
139 .shear_x = shear_x});
143 .def(
"f_aabb", &shammath::paving_function_general_3d_shear_x<f64_3>::f_aabb)
144 .def(
"f_aabb_inv", &shammath::paving_function_general_3d_shear_x<f64_3>::f_aabb_inv)
146 "get_paving_index_intersecting",
147 &shammath::paving_function_general_3d_shear_x<f64_3>::get_paving_index_intersecting);
149 py::class_<f64_4x4>(math_module,
"f64_4x4")
151 return std::make_unique<f64_4x4>();
155 [](
const f64_4x4 &m, std::pair<int, int> idx) ->
double {
156 return m(idx.first, idx.second);
160 [](
f64_4x4 &m, std::pair<int, int> idx,
double value) {
161 m(idx.first, idx.second) = value;
166 std::ostringstream oss;
168 for (
size_t i = 0; i < 4; ++i) {
170 for (
size_t j = 0; j < 4; ++j) {
190 .def(
"to_pyarray", [](
const f64_4x4 &self) {
191 py::array_t<f64> ret({4, 4});
192 for (
u32 i = 0; i < 4; i++) {
193 for (
u32 j = 0; j < 4; j++) {
194 ret.mutable_at(i, j) = self(i, j);
201 math_module.def(
"get_identity_f64_4x4", []() ->
f64_4x4 {
211 math_module.def(
"mat_set_identity", [](
f64_4x4 &a) {
216 py::class_<shammath::SymTensor3d_1<f64>>(math_module,
"SymTensor3d_1_f64")
217 .def(py::init<f64, f64, f64>(), py::arg(
"v_0"), py::arg(
"v_1"), py::arg(
"v_2"))
219 .def_readwrite(
"v_0", &shammath::SymTensor3d_1<f64>::v_0)
220 .def_readwrite(
"v_1", &shammath::SymTensor3d_1<f64>::v_1)
221 .def_readwrite(
"v_2", &shammath::SymTensor3d_1<f64>::v_2)
225 &shammath::SymTensor3d_1<f64>::inner, py::const_),
226 "Inner product with another SymTensor3d_1")
229 py::overload_cast<const f64>(&shammath::SymTensor3d_1<f64>::inner, py::const_),
230 "Scalar multiplication")
237 return fmt::format(
"SymTensor3d_1(v_0={}, v_1={}, v_2={})", t.v_0, t.v_1, t.v_2);
241 py::class_<shammath::SymTensor3d_2<f64>>(math_module,
"SymTensor3d_2_f64")
243 py::init<f64, f64, f64, f64, f64, f64>(),
251 .def_readwrite(
"v_00", &shammath::SymTensor3d_2<f64>::v_00)
252 .def_readwrite(
"v_01", &shammath::SymTensor3d_2<f64>::v_01)
253 .def_readwrite(
"v_02", &shammath::SymTensor3d_2<f64>::v_02)
254 .def_readwrite(
"v_11", &shammath::SymTensor3d_2<f64>::v_11)
255 .def_readwrite(
"v_12", &shammath::SymTensor3d_2<f64>::v_12)
256 .def_readwrite(
"v_22", &shammath::SymTensor3d_2<f64>::v_22)
260 &shammath::SymTensor3d_2<f64>::inner, py::const_),
261 "Inner product with another SymTensor3d_2")
265 &shammath::SymTensor3d_2<f64>::inner, py::const_),
266 "Inner product with SymTensor3d_1")
269 py::overload_cast<const f64>(&shammath::SymTensor3d_2<f64>::inner, py::const_),
270 "Scalar multiplication")
278 "SymTensor3d_2(v_00={}, v_01={}, v_02={}, v_11={}, v_12={}, v_22={})",
288 py::class_<shammath::SymTensor3d_3<f64>>(math_module,
"SymTensor3d_3_f64")
290 py::init<f64, f64, f64, f64, f64, f64, f64, f64, f64, f64>(),
302 .def_readwrite(
"v_000", &shammath::SymTensor3d_3<f64>::v_000)
303 .def_readwrite(
"v_001", &shammath::SymTensor3d_3<f64>::v_001)
304 .def_readwrite(
"v_002", &shammath::SymTensor3d_3<f64>::v_002)
305 .def_readwrite(
"v_011", &shammath::SymTensor3d_3<f64>::v_011)
306 .def_readwrite(
"v_012", &shammath::SymTensor3d_3<f64>::v_012)
307 .def_readwrite(
"v_022", &shammath::SymTensor3d_3<f64>::v_022)
308 .def_readwrite(
"v_111", &shammath::SymTensor3d_3<f64>::v_111)
309 .def_readwrite(
"v_112", &shammath::SymTensor3d_3<f64>::v_112)
310 .def_readwrite(
"v_122", &shammath::SymTensor3d_3<f64>::v_122)
311 .def_readwrite(
"v_222", &shammath::SymTensor3d_3<f64>::v_222)
315 &shammath::SymTensor3d_3<f64>::inner, py::const_),
316 "Inner product with another SymTensor3d_3")
320 &shammath::SymTensor3d_3<f64>::inner, py::const_),
321 "Inner product with SymTensor3d_2")
325 &shammath::SymTensor3d_3<f64>::inner, py::const_),
326 "Inner product with SymTensor3d_1")
329 py::overload_cast<const f64>(&shammath::SymTensor3d_3<f64>::inner, py::const_),
330 "Scalar multiplication")
338 "SymTensor3d_3(v_000={}, v_001={}, v_002={}, v_011={}, v_012={}, v_022={}, "
339 "v_111={}, v_112={}, v_122={}, v_222={})",
353 py::class_<shammath::SymTensor3d_4<f64>>(math_module,
"SymTensor3d_4_f64")
355 py::init<f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64>(),
372 .def_readwrite(
"v_0000", &shammath::SymTensor3d_4<f64>::v_0000)
373 .def_readwrite(
"v_0001", &shammath::SymTensor3d_4<f64>::v_0001)
374 .def_readwrite(
"v_0002", &shammath::SymTensor3d_4<f64>::v_0002)
375 .def_readwrite(
"v_0011", &shammath::SymTensor3d_4<f64>::v_0011)
376 .def_readwrite(
"v_0012", &shammath::SymTensor3d_4<f64>::v_0012)
377 .def_readwrite(
"v_0022", &shammath::SymTensor3d_4<f64>::v_0022)
378 .def_readwrite(
"v_0111", &shammath::SymTensor3d_4<f64>::v_0111)
379 .def_readwrite(
"v_0112", &shammath::SymTensor3d_4<f64>::v_0112)
380 .def_readwrite(
"v_0122", &shammath::SymTensor3d_4<f64>::v_0122)
381 .def_readwrite(
"v_0222", &shammath::SymTensor3d_4<f64>::v_0222)
382 .def_readwrite(
"v_1111", &shammath::SymTensor3d_4<f64>::v_1111)
383 .def_readwrite(
"v_1112", &shammath::SymTensor3d_4<f64>::v_1112)
384 .def_readwrite(
"v_1122", &shammath::SymTensor3d_4<f64>::v_1122)
385 .def_readwrite(
"v_1222", &shammath::SymTensor3d_4<f64>::v_1222)
386 .def_readwrite(
"v_2222", &shammath::SymTensor3d_4<f64>::v_2222)
390 &shammath::SymTensor3d_4<f64>::inner, py::const_),
391 "Inner product with another SymTensor3d_4")
395 &shammath::SymTensor3d_4<f64>::inner, py::const_),
396 "Inner product with SymTensor3d_3")
400 &shammath::SymTensor3d_4<f64>::inner, py::const_),
401 "Inner product with SymTensor3d_2")
405 &shammath::SymTensor3d_4<f64>::inner, py::const_),
406 "Inner product with SymTensor3d_1")
409 py::overload_cast<const f64>(&shammath::SymTensor3d_4<f64>::inner, py::const_),
410 "Scalar multiplication")
418 "SymTensor3d_4(v_0000={}, v_0001={}, v_0002={}, v_0011={}, v_0012={}, v_0022={}, "
419 "v_0111={}, v_0112={}, v_0122={}, v_0222={}, v_1111={}, v_1112={}, v_1122={}, "
420 "v_1222={}, v_2222={})",
439 py::class_<shammath::SymTensor3d_5<f64>>(math_module,
"SymTensor3d_5_f64")
485 .def_readwrite(
"v_00000", &shammath::SymTensor3d_5<f64>::v_00000)
486 .def_readwrite(
"v_00001", &shammath::SymTensor3d_5<f64>::v_00001)
487 .def_readwrite(
"v_00002", &shammath::SymTensor3d_5<f64>::v_00002)
488 .def_readwrite(
"v_00011", &shammath::SymTensor3d_5<f64>::v_00011)
489 .def_readwrite(
"v_00012", &shammath::SymTensor3d_5<f64>::v_00012)
490 .def_readwrite(
"v_00022", &shammath::SymTensor3d_5<f64>::v_00022)
491 .def_readwrite(
"v_00111", &shammath::SymTensor3d_5<f64>::v_00111)
492 .def_readwrite(
"v_00112", &shammath::SymTensor3d_5<f64>::v_00112)
493 .def_readwrite(
"v_00122", &shammath::SymTensor3d_5<f64>::v_00122)
494 .def_readwrite(
"v_00222", &shammath::SymTensor3d_5<f64>::v_00222)
495 .def_readwrite(
"v_01111", &shammath::SymTensor3d_5<f64>::v_01111)
496 .def_readwrite(
"v_01112", &shammath::SymTensor3d_5<f64>::v_01112)
497 .def_readwrite(
"v_01122", &shammath::SymTensor3d_5<f64>::v_01122)
498 .def_readwrite(
"v_01222", &shammath::SymTensor3d_5<f64>::v_01222)
499 .def_readwrite(
"v_02222", &shammath::SymTensor3d_5<f64>::v_02222)
500 .def_readwrite(
"v_11111", &shammath::SymTensor3d_5<f64>::v_11111)
501 .def_readwrite(
"v_11112", &shammath::SymTensor3d_5<f64>::v_11112)
502 .def_readwrite(
"v_11122", &shammath::SymTensor3d_5<f64>::v_11122)
503 .def_readwrite(
"v_11222", &shammath::SymTensor3d_5<f64>::v_11222)
504 .def_readwrite(
"v_12222", &shammath::SymTensor3d_5<f64>::v_12222)
505 .def_readwrite(
"v_22222", &shammath::SymTensor3d_5<f64>::v_22222)
509 &shammath::SymTensor3d_5<f64>::inner, py::const_),
510 "Inner product with another SymTensor3d_5")
514 &shammath::SymTensor3d_5<f64>::inner, py::const_),
515 "Inner product with SymTensor3d_4")
519 &shammath::SymTensor3d_5<f64>::inner, py::const_),
520 "Inner product with SymTensor3d_3")
524 &shammath::SymTensor3d_5<f64>::inner, py::const_),
525 "Inner product with SymTensor3d_2")
529 &shammath::SymTensor3d_5<f64>::inner, py::const_),
530 "Inner product with SymTensor3d_1")
533 py::overload_cast<const f64>(&shammath::SymTensor3d_5<f64>::inner, py::const_),
534 "Scalar multiplication")
542 "SymTensor3d_5(v_00000={}, v_00001={}, v_00002={}, v_00011={}, v_00012={}, "
543 "v_00022={}, v_00111={}, v_00112={}, v_00122={}, v_00222={}, v_01111={}, "
544 "v_01112={}, v_01122={}, v_01222={}, v_02222={}, v_11111={}, v_11112={}, "
545 "v_11122={}, v_11222={}, v_12222={}, v_22222={})",
571 py::class_<shammath::SymTensorCollection<f64, 0, 5>>(math_module,
"SymTensorCollection_f64_0_5")
600 "SymTensorCollection_f64_0_5(\n t0={},\n t1={},\n t2={},\n t3={},\n t4={},\n "
603 py::str(py::cast(c.t1)).cast<std::string>(),
604 py::str(py::cast(c.t2)).cast<std::string>(),
605 py::str(py::cast(c.t3)).cast<std::string>(),
606 py::str(py::cast(c.t4)).cast<std::string>(),
607 py::str(py::cast(c.t5)).cast<std::string>());
611 py::class_<shammath::SymTensorCollection<f64, 0, 4>>(math_module,
"SymTensorCollection_f64_0_4")
637 "SymTensorCollection_f64_0_4(\n t0={},\n t1={},\n t2={},\n t3={},\n t4={}\n)",
639 py::str(py::cast(c.t1)).cast<std::string>(),
640 py::str(py::cast(c.t2)).cast<std::string>(),
641 py::str(py::cast(c.t3)).cast<std::string>(),
642 py::str(py::cast(c.t4)).cast<std::string>());
646 py::class_<shammath::SymTensorCollection<f64, 0, 3>>(math_module,
"SymTensorCollection_f64_0_3")
669 "SymTensorCollection_f64_0_3(\n t0={},\n t1={},\n t2={},\n t3={}\n)",
671 py::str(py::cast(c.t1)).cast<std::string>(),
672 py::str(py::cast(c.t2)).cast<std::string>(),
673 py::str(py::cast(c.t3)).cast<std::string>());
677 py::class_<shammath::SymTensorCollection<f64, 0, 2>>(math_module,
"SymTensorCollection_f64_0_2")
694 "SymTensorCollection_f64_0_2(\n t0={},\n t1={},\n t2={}\n)",
696 py::str(py::cast(c.t1)).cast<std::string>(),
697 py::str(py::cast(c.t2)).cast<std::string>());
701 py::class_<shammath::SymTensorCollection<f64, 0, 1>>(math_module,
"SymTensorCollection_f64_0_1")
713 "SymTensorCollection_f64_0_1(\n t0={},\n t1={}\n)",
715 py::str(py::cast(c.t1)).cast<std::string>());
719 py::class_<shammath::SymTensorCollection<f64, 0, 0>>(math_module,
"SymTensorCollection_f64_0_0")
721 .def(py::init<f64>(), py::arg(
"t0"))
729 return fmt::format(
"SymTensorCollection_f64_0_0(t0={})", c.t0);
733 py::class_<shammath::SymTensorCollection<f64, 1, 5>>(math_module,
"SymTensorCollection_f64_1_5")
759 "SymTensorCollection_f64_1_5(\n t1={},\n t2={},\n t3={},\n t4={},\n t5={}\n)",
760 py::str(py::cast(c.t1)).cast<std::string>(),
761 py::str(py::cast(c.t2)).cast<std::string>(),
762 py::str(py::cast(c.t3)).cast<std::string>(),
763 py::str(py::cast(c.t4)).cast<std::string>(),
764 py::str(py::cast(c.t5)).cast<std::string>());
768 py::class_<shammath::SymTensorCollection<f64, 1, 4>>(math_module,
"SymTensorCollection_f64_1_4")
791 "SymTensorCollection_f64_1_4(\n t1={},\n t2={},\n t3={},\n t4={}\n)",
792 py::str(py::cast(c.t1)).cast<std::string>(),
793 py::str(py::cast(c.t2)).cast<std::string>(),
794 py::str(py::cast(c.t3)).cast<std::string>(),
795 py::str(py::cast(c.t4)).cast<std::string>());
799 py::class_<shammath::SymTensorCollection<f64, 1, 3>>(math_module,
"SymTensorCollection_f64_1_3")
819 "SymTensorCollection_f64_1_3(\n t1={},\n t2={},\n t3={}\n)",
820 py::str(py::cast(c.t1)).cast<std::string>(),
821 py::str(py::cast(c.t2)).cast<std::string>(),
822 py::str(py::cast(c.t3)).cast<std::string>());
826 py::class_<shammath::SymTensorCollection<f64, 1, 2>>(math_module,
"SymTensorCollection_f64_1_2")
841 "SymTensorCollection_f64_1_2(\n t1={},\n t2={}\n)",
842 py::str(py::cast(c.t1)).cast<std::string>(),
843 py::str(py::cast(c.t2)).cast<std::string>());
847 py::class_<shammath::SymTensorCollection<f64, 1, 1>>(math_module,
"SymTensorCollection_f64_1_1")
858 "SymTensorCollection_f64_1_1(\n t1={}\n)",
859 py::str(py::cast(c.t1)).cast<std::string>());
875 Solve ODE with Euler method
876 start : Lower bound of integration
877 end : Higher bound of integration
878 step : Step of integration
880 x0 : Initial coordinate
886 [](
const std::function<
f64(
const std::vector<f64> &,
f64)> &func,
887 const std::vector<f64> &x_data,
888 const std::vector<f64> &y_data,
889 const std::vector<f64> &p0) {
898 Fit data with a given function by least squares method
899 f: Function (1d values)
902 p0: Initial parameters estimated
905 math_module.def("get_ideal_hcp_box", [](
f64 dr, f64_3 box_min, f64_3 box_max) {
906 return shammath::LatticeHCP<f64_3>::get_ideal_hcp_box(dr, {box_min, box_max});
910 "get_periodic_hcp_box",
911 [](
f64 dr, std::array<i32, 3> box_min, std::array<i32, 3> box_max) {
913 return std::tuple<f64_3, f64_3>{ret.lower, ret.upper};
919 Get the periodic box corresponding to integer lattice coordinates
920 this function will throw if the coordinates asked cannot make a periodic lattice
923 dr: the particle spacing in the lattice
924 box_min: integer triplet for the minimal coordinates on the lattice
925 box_max: integer triplet for the maximal coordinates on the lattice
double f64
Alias for double.
std::uint32_t u32
32 bit unsigned integer
static constexpr CoordRange< Tvec > get_periodic_box(Tscal dr, std::array< i32, dim > coord_min, std::array< i32, dim > coord_max)
Get the periodic box corresponding to integer lattice coordinates this function will throw if the coo...
constexpr auto get_mdspan()
Get the matrix data as a mdspan.
shammath::mat< f64, 4, 4 > f64_4x4
Alias for 4x4 double matrix.
constexpr mat< T, n, n > mat_identity()
Returns the identity matrix of size n.
T estim_deriv_step(u32 order)
Estimate the best step size for numerical differentiation of given order.
T derivative_3point_forward(T x, T dx, std::function< T(T)> &&fct)
Compute the derivative of a function at x using a 3-point forward finite difference.
constexpr std::pair< std::vector< T >, std::vector< T > > euler_ode(T start, T end, T step, Lambda &&ode, T x0, T u0)
Euler solving of ODE The ode has the form.
T derivative_3point_backward(T x, T dx, std::function< T(T)> &&fct)
Compute the derivative of a function at x using a 3-point backward finite difference.
std::pair< std::vector< T >, T > least_squares(const Lambda &f, const std::vector< T > &X, const std::vector< T > &Y, const std::vector< T > &p0, int maxits=1000, T tolerance=1e-9)
This function determines the best fit parameters for a given function with least squares.
T derivative_5point_midpoint(T x, T dx, std::function< T(T)> &&fct)
Compute the derivative of a function at x using a 5-point centered finite difference.
T derivative_upwind(T x, T dx, std::function< T(T)> &&fct)
Compute the derivative of a function at x using the upwind method.
void mat_set_identity(const std::mdspan< T, Extents, Layout, Accessor > &input1)
Set the content of a matrix to the identity matrix.
void mat_prod(const std::mdspan< Ta, Extents1, Layout1, Accessor1 > &input1, const std::mdspan< Ta, Extents2, Layout2, Accessor2 > &input2, const std::mdspan< Tb, Extents3, Layout3, Accessor3 > &output)
Compute the product of two matrices.
T derivative_centered(T x, T dx, std::function< T(T)> &&fct)
Compute the derivative of a function at x using the centered difference method.
Pybind11 include and definitions.
#define ON_PYTHON_INIT
Register a Python module init function using static initialization.
A structure for 3D paving functions with shearing along the x-axis and general boundary conditions.
Tvec f_inv(Tvec x, int i, int j, int k) const
Applies the inverse of the paving function with shearing and boundary conditions.
Tvec f(Tvec x, int i, int j, int k) const
Applies the paving function with shearing and boundary conditions.
A structure for 3D paving functions with general boundary conditions (periodic or reflective per dire...
Tvec f_inv(Tvec x, int i, int j, int k) const
Applies the inverse of the paving function.
Tvec f(Tvec x, int i, int j, int k) const
Applies the paving function with periodic or reflective boundary conditions.
A structure for 3D paving functions with periodic boundary conditions.
Tvec f_inv(Tvec x, int i, int j, int k) const
Applies the inverse of the paving function.
Tvec f(Tvec x, int i, int j, int k) const
Applies the paving function with periodic boundary conditions.