Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
pyshammath.cpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
17
25#include "shammath/matrix.hpp"
28#include "shammath/solve.hpp"
35#include <fmt/core.h>
36#include <pybind11/cast.h>
37#include <pybind11/functional.h>
38#include <pybind11/numpy.h>
39#include <functional>
40
42
43 py::module math_module = root_module.def_submodule("math", "Shamrock math lib");
44
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);
50
51 math_module.def("derivative_upwind", [](f64 x, f64 dx, std::function<f64(f64)> fct) {
52 return shammath::derivative_upwind<f64>(x, dx, [&](f64 x) {
53 return fct(x);
54 });
55 });
56 math_module.def("derivative_centered", [](f64 x, f64 dx, std::function<f64(f64)> fct) {
57 return shammath::derivative_centered<f64>(x, dx, [&](f64 x) {
58 return fct(x);
59 });
60 });
61 math_module.def("derivative_3point_forward", [](f64 x, f64 dx, std::function<f64(f64)> fct) {
63 return fct(x);
64 });
65 });
66 math_module.def("derivative_3point_backward", [](f64 x, f64 dx, std::function<f64(f64)> fct) {
68 return fct(x);
69 });
70 });
71 math_module.def("derivative_5point_midpoint", [](f64 x, f64 dx, std::function<f64(f64)> fct) {
73 return fct(x);
74 });
75 });
76 math_module.def(
77 "estim_deriv_step",
78 [](u32 order) {
80 },
81 R"pbdoc(
82 Estim the correct step to use for a given order when using derivatives
83 )pbdoc");
84
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>>(
90 }))
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)
95 .def(
96 "get_paving_index_intersecting",
97 &shammath::paving_function_periodic_3d<f64_3>::get_paving_index_intersecting);
98
99 py::class_<shammath::paving_function_general_3d<f64_3>>(
100 math_module, "paving_function_general_3d")
101 .def(
102 py::init([](f64_3 box_size,
103 f64_3 box_center,
104 bool is_x_periodic,
105 bool is_y_periodic,
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});
114 }))
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)
119 .def(
120 "get_paving_index_intersecting",
121 &shammath::paving_function_general_3d<f64_3>::get_paving_index_intersecting);
122
123 py::class_<shammath::paving_function_general_3d_shear_x<f64_3>>(
124 math_module, "paving_function_general_3d_shear_x")
125 .def(
126 py::init([](f64_3 box_size,
127 f64_3 box_center,
128 bool is_x_periodic,
129 bool is_y_periodic,
130 bool is_z_periodic,
131 f64 shear_x) {
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});
140 }))
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)
145 .def(
146 "get_paving_index_intersecting",
147 &shammath::paving_function_general_3d_shear_x<f64_3>::get_paving_index_intersecting);
148
149 py::class_<f64_4x4>(math_module, "f64_4x4")
150 .def(py::init([]() {
151 return std::make_unique<f64_4x4>();
152 }))
153 .def(
154 "__getitem__",
155 [](const f64_4x4 &m, std::pair<int, int> idx) -> double {
156 return m(idx.first, idx.second);
157 })
158 .def(
159 "__setitem__",
160 [](f64_4x4 &m, std::pair<int, int> idx, double value) {
161 m(idx.first, idx.second) = value;
162 })
163 .def(
164 "__repr__",
165 [](const f64_4x4 &m) {
166 std::ostringstream oss;
167 oss << "[";
168 for (size_t i = 0; i < 4; ++i) {
169 oss << "[";
170 for (size_t j = 0; j < 4; ++j) {
171 oss << m(i, j);
172 if (j + 1 < 4)
173 oss << ", ";
174 }
175 oss << "]";
176 if (i + 1 < 4)
177 oss << ",\n ";
178 }
179 oss << "]";
180 return oss.str();
181 })
182 .def(
183 "__matmul__",
184 [](const f64_4x4 &a, const f64_4x4 &b) {
185 f64_4x4 ret;
187 return ret;
188 },
189 py::is_operator())
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);
195 }
196 }
197
198 return ret;
199 });
200
201 math_module.def("get_identity_f64_4x4", []() -> f64_4x4 {
203 });
204
205 math_module.def("mat_mul", [](const f64_4x4 &a, const f64_4x4 &b) -> f64_4x4 {
206 f64_4x4 ret;
208 return ret;
209 });
210
211 math_module.def("mat_set_identity", [](f64_4x4 &a) {
213 });
214
215 // SymTensor3d_1 bindings
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"))
218 .def(py::init<>())
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)
222 .def(
223 "inner",
224 py::overload_cast<const shammath::SymTensor3d_1<f64> &>(
225 &shammath::SymTensor3d_1<f64>::inner, py::const_),
226 "Inner product with another SymTensor3d_1")
227 .def(
228 "inner",
229 py::overload_cast<const f64>(&shammath::SymTensor3d_1<f64>::inner, py::const_),
230 "Scalar multiplication")
231 .def("__mul__", &shammath::SymTensor3d_1<f64>::operator*, "Multiply by scalar")
232 .def("__imul__", &shammath::SymTensor3d_1<f64>::operator*=, "In-place multiply by scalar")
233 .def("__add__", &shammath::SymTensor3d_1<f64>::operator+, "Add two tensors")
234 .def("__iadd__", &shammath::SymTensor3d_1<f64>::operator+=, "In-place add")
235 .def("__sub__", &shammath::SymTensor3d_1<f64>::operator-, "Subtract two tensors")
236 .def("__repr__", [](const shammath::SymTensor3d_1<f64> &t) {
237 return fmt::format("SymTensor3d_1(v_0={}, v_1={}, v_2={})", t.v_0, t.v_1, t.v_2);
238 });
239
240 // SymTensor3d_2 bindings
241 py::class_<shammath::SymTensor3d_2<f64>>(math_module, "SymTensor3d_2_f64")
242 .def(
243 py::init<f64, f64, f64, f64, f64, f64>(),
244 py::arg("v_00"),
245 py::arg("v_01"),
246 py::arg("v_02"),
247 py::arg("v_11"),
248 py::arg("v_12"),
249 py::arg("v_22"))
250 .def(py::init<>())
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)
257 .def(
258 "inner",
259 py::overload_cast<const shammath::SymTensor3d_2<f64> &>(
260 &shammath::SymTensor3d_2<f64>::inner, py::const_),
261 "Inner product with another SymTensor3d_2")
262 .def(
263 "inner",
264 py::overload_cast<const shammath::SymTensor3d_1<f64> &>(
265 &shammath::SymTensor3d_2<f64>::inner, py::const_),
266 "Inner product with SymTensor3d_1")
267 .def(
268 "inner",
269 py::overload_cast<const f64>(&shammath::SymTensor3d_2<f64>::inner, py::const_),
270 "Scalar multiplication")
271 .def("__mul__", &shammath::SymTensor3d_2<f64>::operator*, "Multiply by scalar")
272 .def("__imul__", &shammath::SymTensor3d_2<f64>::operator*=, "In-place multiply by scalar")
273 .def("__add__", &shammath::SymTensor3d_2<f64>::operator+, "Add two tensors")
274 .def("__iadd__", &shammath::SymTensor3d_2<f64>::operator+=, "In-place add")
275 .def("__sub__", &shammath::SymTensor3d_2<f64>::operator-, "Subtract two tensors")
276 .def("__repr__", [](const shammath::SymTensor3d_2<f64> &t) {
277 return fmt::format(
278 "SymTensor3d_2(v_00={}, v_01={}, v_02={}, v_11={}, v_12={}, v_22={})",
279 t.v_00,
280 t.v_01,
281 t.v_02,
282 t.v_11,
283 t.v_12,
284 t.v_22);
285 });
286
287 // SymTensor3d_3 bindings
288 py::class_<shammath::SymTensor3d_3<f64>>(math_module, "SymTensor3d_3_f64")
289 .def(
290 py::init<f64, f64, f64, f64, f64, f64, f64, f64, f64, f64>(),
291 py::arg("v_000"),
292 py::arg("v_001"),
293 py::arg("v_002"),
294 py::arg("v_011"),
295 py::arg("v_012"),
296 py::arg("v_022"),
297 py::arg("v_111"),
298 py::arg("v_112"),
299 py::arg("v_122"),
300 py::arg("v_222"))
301 .def(py::init<>())
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)
312 .def(
313 "inner",
314 py::overload_cast<const shammath::SymTensor3d_3<f64> &>(
315 &shammath::SymTensor3d_3<f64>::inner, py::const_),
316 "Inner product with another SymTensor3d_3")
317 .def(
318 "inner",
319 py::overload_cast<const shammath::SymTensor3d_2<f64> &>(
320 &shammath::SymTensor3d_3<f64>::inner, py::const_),
321 "Inner product with SymTensor3d_2")
322 .def(
323 "inner",
324 py::overload_cast<const shammath::SymTensor3d_1<f64> &>(
325 &shammath::SymTensor3d_3<f64>::inner, py::const_),
326 "Inner product with SymTensor3d_1")
327 .def(
328 "inner",
329 py::overload_cast<const f64>(&shammath::SymTensor3d_3<f64>::inner, py::const_),
330 "Scalar multiplication")
331 .def("__mul__", &shammath::SymTensor3d_3<f64>::operator*, "Multiply by scalar")
332 .def("__imul__", &shammath::SymTensor3d_3<f64>::operator*=, "In-place multiply by scalar")
333 .def("__add__", &shammath::SymTensor3d_3<f64>::operator+, "Add two tensors")
334 .def("__iadd__", &shammath::SymTensor3d_3<f64>::operator+=, "In-place add")
335 .def("__sub__", &shammath::SymTensor3d_3<f64>::operator-, "Subtract two tensors")
336 .def("__repr__", [](const shammath::SymTensor3d_3<f64> &t) {
337 return fmt::format(
338 "SymTensor3d_3(v_000={}, v_001={}, v_002={}, v_011={}, v_012={}, v_022={}, "
339 "v_111={}, v_112={}, v_122={}, v_222={})",
340 t.v_000,
341 t.v_001,
342 t.v_002,
343 t.v_011,
344 t.v_012,
345 t.v_022,
346 t.v_111,
347 t.v_112,
348 t.v_122,
349 t.v_222);
350 });
351
352 // SymTensor3d_4 bindings
353 py::class_<shammath::SymTensor3d_4<f64>>(math_module, "SymTensor3d_4_f64")
354 .def(
355 py::init<f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64, f64>(),
356 py::arg("v_0000"),
357 py::arg("v_0001"),
358 py::arg("v_0002"),
359 py::arg("v_0011"),
360 py::arg("v_0012"),
361 py::arg("v_0022"),
362 py::arg("v_0111"),
363 py::arg("v_0112"),
364 py::arg("v_0122"),
365 py::arg("v_0222"),
366 py::arg("v_1111"),
367 py::arg("v_1112"),
368 py::arg("v_1122"),
369 py::arg("v_1222"),
370 py::arg("v_2222"))
371 .def(py::init<>())
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)
387 .def(
388 "inner",
389 py::overload_cast<const shammath::SymTensor3d_4<f64> &>(
390 &shammath::SymTensor3d_4<f64>::inner, py::const_),
391 "Inner product with another SymTensor3d_4")
392 .def(
393 "inner",
394 py::overload_cast<const shammath::SymTensor3d_3<f64> &>(
395 &shammath::SymTensor3d_4<f64>::inner, py::const_),
396 "Inner product with SymTensor3d_3")
397 .def(
398 "inner",
399 py::overload_cast<const shammath::SymTensor3d_2<f64> &>(
400 &shammath::SymTensor3d_4<f64>::inner, py::const_),
401 "Inner product with SymTensor3d_2")
402 .def(
403 "inner",
404 py::overload_cast<const shammath::SymTensor3d_1<f64> &>(
405 &shammath::SymTensor3d_4<f64>::inner, py::const_),
406 "Inner product with SymTensor3d_1")
407 .def(
408 "inner",
409 py::overload_cast<const f64>(&shammath::SymTensor3d_4<f64>::inner, py::const_),
410 "Scalar multiplication")
411 .def("__mul__", &shammath::SymTensor3d_4<f64>::operator*, "Multiply by scalar")
412 .def("__imul__", &shammath::SymTensor3d_4<f64>::operator*=, "In-place multiply by scalar")
413 .def("__add__", &shammath::SymTensor3d_4<f64>::operator+, "Add two tensors")
414 .def("__iadd__", &shammath::SymTensor3d_4<f64>::operator+=, "In-place add")
415 .def("__sub__", &shammath::SymTensor3d_4<f64>::operator-, "Subtract two tensors")
416 .def("__repr__", [](const shammath::SymTensor3d_4<f64> &t) {
417 return fmt::format(
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={})",
421 t.v_0000,
422 t.v_0001,
423 t.v_0002,
424 t.v_0011,
425 t.v_0012,
426 t.v_0022,
427 t.v_0111,
428 t.v_0112,
429 t.v_0122,
430 t.v_0222,
431 t.v_1111,
432 t.v_1112,
433 t.v_1122,
434 t.v_1222,
435 t.v_2222);
436 });
437
438 // SymTensor3d_5 bindings
439 py::class_<shammath::SymTensor3d_5<f64>>(math_module, "SymTensor3d_5_f64")
440 .def(
441 py::init<
442 f64,
443 f64,
444 f64,
445 f64,
446 f64,
447 f64,
448 f64,
449 f64,
450 f64,
451 f64,
452 f64,
453 f64,
454 f64,
455 f64,
456 f64,
457 f64,
458 f64,
459 f64,
460 f64,
461 f64,
462 f64>(),
463 py::arg("v_00000"),
464 py::arg("v_00001"),
465 py::arg("v_00002"),
466 py::arg("v_00011"),
467 py::arg("v_00012"),
468 py::arg("v_00022"),
469 py::arg("v_00111"),
470 py::arg("v_00112"),
471 py::arg("v_00122"),
472 py::arg("v_00222"),
473 py::arg("v_01111"),
474 py::arg("v_01112"),
475 py::arg("v_01122"),
476 py::arg("v_01222"),
477 py::arg("v_02222"),
478 py::arg("v_11111"),
479 py::arg("v_11112"),
480 py::arg("v_11122"),
481 py::arg("v_11222"),
482 py::arg("v_12222"),
483 py::arg("v_22222"))
484 .def(py::init<>())
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)
506 .def(
507 "inner",
508 py::overload_cast<const shammath::SymTensor3d_5<f64> &>(
509 &shammath::SymTensor3d_5<f64>::inner, py::const_),
510 "Inner product with another SymTensor3d_5")
511 .def(
512 "inner",
513 py::overload_cast<const shammath::SymTensor3d_4<f64> &>(
514 &shammath::SymTensor3d_5<f64>::inner, py::const_),
515 "Inner product with SymTensor3d_4")
516 .def(
517 "inner",
518 py::overload_cast<const shammath::SymTensor3d_3<f64> &>(
519 &shammath::SymTensor3d_5<f64>::inner, py::const_),
520 "Inner product with SymTensor3d_3")
521 .def(
522 "inner",
523 py::overload_cast<const shammath::SymTensor3d_2<f64> &>(
524 &shammath::SymTensor3d_5<f64>::inner, py::const_),
525 "Inner product with SymTensor3d_2")
526 .def(
527 "inner",
528 py::overload_cast<const shammath::SymTensor3d_1<f64> &>(
529 &shammath::SymTensor3d_5<f64>::inner, py::const_),
530 "Inner product with SymTensor3d_1")
531 .def(
532 "inner",
533 py::overload_cast<const f64>(&shammath::SymTensor3d_5<f64>::inner, py::const_),
534 "Scalar multiplication")
535 .def("__mul__", &shammath::SymTensor3d_5<f64>::operator*, "Multiply by scalar")
536 .def("__imul__", &shammath::SymTensor3d_5<f64>::operator*=, "In-place multiply by scalar")
537 .def("__add__", &shammath::SymTensor3d_5<f64>::operator+, "Add two tensors")
538 .def("__iadd__", &shammath::SymTensor3d_5<f64>::operator+=, "In-place add")
539 .def("__sub__", &shammath::SymTensor3d_5<f64>::operator-, "Subtract two tensors")
540 .def("__repr__", [](const shammath::SymTensor3d_5<f64> &t) {
541 return fmt::format(
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={})",
546 t.v_00000,
547 t.v_00001,
548 t.v_00002,
549 t.v_00011,
550 t.v_00012,
551 t.v_00022,
552 t.v_00111,
553 t.v_00112,
554 t.v_00122,
555 t.v_00222,
556 t.v_01111,
557 t.v_01112,
558 t.v_01122,
559 t.v_01222,
560 t.v_02222,
561 t.v_11111,
562 t.v_11112,
563 t.v_11122,
564 t.v_11222,
565 t.v_12222,
566 t.v_22222);
567 });
568
569 // SymTensorCollection bindings
570 // SymTensorCollection<f64, 0, 5>
571 py::class_<shammath::SymTensorCollection<f64, 0, 5>>(math_module, "SymTensorCollection_f64_0_5")
572 .def(py::init<>())
573 .def(
574 py::init<
575 f64,
581 py::arg("t0"),
582 py::arg("t1"),
583 py::arg("t2"),
584 py::arg("t3"),
585 py::arg("t4"),
586 py::arg("t5"))
594 .def_static("from_vec", &shammath::SymTensorCollection<f64, 0, 5>::from_vec, py::arg("v"))
598 .def("__repr__", [](const shammath::SymTensorCollection<f64, 0, 5> &c) {
599 return fmt::format(
600 "SymTensorCollection_f64_0_5(\n t0={},\n t1={},\n t2={},\n t3={},\n t4={},\n "
601 "t5={}\n)",
602 c.t0,
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>());
608 });
609
610 // SymTensorCollection<f64, 0, 4>
611 py::class_<shammath::SymTensorCollection<f64, 0, 4>>(math_module, "SymTensorCollection_f64_0_4")
612 .def(py::init<>())
613 .def(
614 py::init<
615 f64,
620 py::arg("t0"),
621 py::arg("t1"),
622 py::arg("t2"),
623 py::arg("t3"),
624 py::arg("t4"))
631 .def_static("from_vec", &shammath::SymTensorCollection<f64, 0, 4>::from_vec, py::arg("v"))
635 .def("__repr__", [](const shammath::SymTensorCollection<f64, 0, 4> &c) {
636 return fmt::format(
637 "SymTensorCollection_f64_0_4(\n t0={},\n t1={},\n t2={},\n t3={},\n t4={}\n)",
638 c.t0,
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>());
643 });
644
645 // SymTensorCollection<f64, 0, 3>
646 py::class_<shammath::SymTensorCollection<f64, 0, 3>>(math_module, "SymTensorCollection_f64_0_3")
647 .def(py::init<>())
648 .def(
649 py::init<
650 f64,
654 py::arg("t0"),
655 py::arg("t1"),
656 py::arg("t2"),
657 py::arg("t3"))
663 .def_static("from_vec", &shammath::SymTensorCollection<f64, 0, 3>::from_vec, py::arg("v"))
667 .def("__repr__", [](const shammath::SymTensorCollection<f64, 0, 3> &c) {
668 return fmt::format(
669 "SymTensorCollection_f64_0_3(\n t0={},\n t1={},\n t2={},\n t3={}\n)",
670 c.t0,
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>());
674 });
675
676 // SymTensorCollection<f64, 0, 2>
677 py::class_<shammath::SymTensorCollection<f64, 0, 2>>(math_module, "SymTensorCollection_f64_0_2")
678 .def(py::init<>())
679 .def(
681 py::arg("t0"),
682 py::arg("t1"),
683 py::arg("t2"))
688 .def_static("from_vec", &shammath::SymTensorCollection<f64, 0, 2>::from_vec, py::arg("v"))
692 .def("__repr__", [](const shammath::SymTensorCollection<f64, 0, 2> &c) {
693 return fmt::format(
694 "SymTensorCollection_f64_0_2(\n t0={},\n t1={},\n t2={}\n)",
695 c.t0,
696 py::str(py::cast(c.t1)).cast<std::string>(),
697 py::str(py::cast(c.t2)).cast<std::string>());
698 });
699
700 // SymTensorCollection<f64, 0, 1>
701 py::class_<shammath::SymTensorCollection<f64, 0, 1>>(math_module, "SymTensorCollection_f64_0_1")
702 .def(py::init<>())
703 .def(py::init<f64, shammath::SymTensor3d_1<f64>>(), py::arg("t0"), py::arg("t1"))
707 .def_static("from_vec", &shammath::SymTensorCollection<f64, 0, 1>::from_vec, py::arg("v"))
711 .def("__repr__", [](const shammath::SymTensorCollection<f64, 0, 1> &c) {
712 return fmt::format(
713 "SymTensorCollection_f64_0_1(\n t0={},\n t1={}\n)",
714 c.t0,
715 py::str(py::cast(c.t1)).cast<std::string>());
716 });
717
718 // SymTensorCollection<f64, 0, 0>
719 py::class_<shammath::SymTensorCollection<f64, 0, 0>>(math_module, "SymTensorCollection_f64_0_0")
720 .def(py::init<>())
721 .def(py::init<f64>(), py::arg("t0"))
724 .def_static("from_vec", &shammath::SymTensorCollection<f64, 0, 0>::from_vec, py::arg("v"))
728 .def("__repr__", [](const shammath::SymTensorCollection<f64, 0, 0> &c) {
729 return fmt::format("SymTensorCollection_f64_0_0(t0={})", c.t0);
730 });
731
732 // SymTensorCollection<f64, 1, 5>
733 py::class_<shammath::SymTensorCollection<f64, 1, 5>>(math_module, "SymTensorCollection_f64_1_5")
734 .def(py::init<>())
735 .def(
736 py::init<
742 py::arg("t1"),
743 py::arg("t2"),
744 py::arg("t3"),
745 py::arg("t4"),
746 py::arg("t5"))
753 .def_static("from_vec", &shammath::SymTensorCollection<f64, 1, 5>::from_vec, py::arg("v"))
757 .def("__repr__", [](const shammath::SymTensorCollection<f64, 1, 5> &c) {
758 return fmt::format(
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>());
765 });
766
767 // SymTensorCollection<f64, 1, 4>
768 py::class_<shammath::SymTensorCollection<f64, 1, 4>>(math_module, "SymTensorCollection_f64_1_4")
769 .def(py::init<>())
770 .def(
771 py::init<
776 py::arg("t1"),
777 py::arg("t2"),
778 py::arg("t3"),
779 py::arg("t4"))
785 .def_static("from_vec", &shammath::SymTensorCollection<f64, 1, 4>::from_vec, py::arg("v"))
789 .def("__repr__", [](const shammath::SymTensorCollection<f64, 1, 4> &c) {
790 return fmt::format(
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>());
796 });
797
798 // SymTensorCollection<f64, 1, 3>
799 py::class_<shammath::SymTensorCollection<f64, 1, 3>>(math_module, "SymTensorCollection_f64_1_3")
800 .def(py::init<>())
801 .def(
802 py::init<
806 py::arg("t1"),
807 py::arg("t2"),
808 py::arg("t3"))
813 .def_static("from_vec", &shammath::SymTensorCollection<f64, 1, 3>::from_vec, py::arg("v"))
817 .def("__repr__", [](const shammath::SymTensorCollection<f64, 1, 3> &c) {
818 return fmt::format(
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>());
823 });
824
825 // SymTensorCollection<f64, 1, 2>
826 py::class_<shammath::SymTensorCollection<f64, 1, 2>>(math_module, "SymTensorCollection_f64_1_2")
827 .def(py::init<>())
828 .def(
830 py::arg("t1"),
831 py::arg("t2"))
835 .def_static("from_vec", &shammath::SymTensorCollection<f64, 1, 2>::from_vec, py::arg("v"))
839 .def("__repr__", [](const shammath::SymTensorCollection<f64, 1, 2> &c) {
840 return fmt::format(
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>());
844 });
845
846 // SymTensorCollection<f64, 1, 1>
847 py::class_<shammath::SymTensorCollection<f64, 1, 1>>(math_module, "SymTensorCollection_f64_1_1")
848 .def(py::init<>())
849 .def(py::init<shammath::SymTensor3d_1<f64>>(), py::arg("t1"))
852 .def_static("from_vec", &shammath::SymTensorCollection<f64, 1, 1>::from_vec, py::arg("v"))
856 .def("__repr__", [](const shammath::SymTensorCollection<f64, 1, 1> &c) {
857 return fmt::format(
858 "SymTensorCollection_f64_1_1(\n t1={}\n)",
859 py::str(py::cast(c.t1)).cast<std::string>());
860 });
861
862 math_module.def(
863 "euler_ode",
864 [](f64 start, f64 end, f64 step, std::function<f64(f64, f64)> &&ode, f64 x0, f64 u0) {
865 return shammath::euler_ode<f64>(start, end, step, ode, x0, u0);
866 },
867 py::kw_only(),
868 py::arg("start"),
869 py::arg("end"),
870 py::arg("step"),
871 py::arg("ode"),
872 py::arg("x0"),
873 py::arg("u0"),
874 R"pbdoc(
875 Solve ODE with Euler method
876 start : Lower bound of integration
877 end : Higher bound of integration
878 step : Step of integration
879 ode : Ode function
880 x0 : Initial coordinate
881 u0 : Initial value
882 )pbdoc");
883
884 math_module.def(
885 "least_squares",
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) {
890 return shammath::least_squares(func, x_data, y_data, p0);
891 },
892 py::kw_only(),
893 py::arg("func"),
894 py::arg("x_data"),
895 py::arg("y_data"),
896 py::arg("p0"),
897 R"pbdoc(
898 Fit data with a given function by least squares method
899 f: Function (1d values)
900 X: $x$ Data to fit
901 Y: $y$ Data to fit
902 p0: Initial parameters estimated
903 )pbdoc");
904
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});
907 });
908
909 math_module.def(
910 "get_periodic_hcp_box",
911 [](f64 dr, std::array<i32, 3> box_min, std::array<i32, 3> box_max) {
912 auto ret = shammath::LatticeHCP<f64_3>::get_periodic_box(dr, box_min, box_max);
913 return std::tuple<f64_3, f64_3>{ret.lower, ret.upper};
914 },
915 py::arg("dr"),
916 py::arg("box_min"),
917 py::arg("box_max"),
918 R"pbdoc(
919 Get the periodic box corresponding to integer lattice coordinates
920 this function will throw if the coordinates asked cannot make a periodic lattice
921
922 Args:
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
926 )pbdoc");
927}
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.
Definition matrix.hpp:43
shammath::mat< f64, 4, 4 > f64_4x4
Alias for 4x4 double matrix.
Definition matrix.hpp:255
constexpr mat< T, n, n > mat_identity()
Returns the identity matrix of size n.
Definition matrix.hpp:86
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.
Definition solve.hpp:70
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.