[Math] Add least squares method by DavidFang03 · Pull Request #1591 · Shamrock-code/Shamrock

DavidFang03 · 2026-01-25T13:34:28Z

This PR implements

Least squares method (Levenberg-Marquardt)
Cholesky decomposition (required for least squares)
Some matrices handling: transpose and runtime dimensions (required for least squares)
Euler ODE solving

This will be useful for #1445 as I would like to numerically solve $\frac{du}{d \rho} = p(u, \rho)/\rho^2$ and fit with a non-linear function

I wasn't very familiar with std::mdspan and The Levenberg-Marquardt has room for improvement but I believe it works well enough for #1145

github-actions · 2026-01-25T13:34:36Z

Thanks @DavidFang03 for opening this PR!

You can do multiple things directly here:
1 - Comment pre-commit.ci run to run pre-commit checks.
2 - Comment pre-commit.ci autofix to apply fixes.
3 - Add label autofix.ci to fix authorship & pre-commit for every commit made.
4 - Add label light-ci to only trigger a reduced & faster version of the CI (need the full one before merge).
5 - Add label trigger-ci to create an empty commit to trigger the CI.

Once the workflow completes a message will appear displaying informations related to the run.

Also the PR gets automatically reviewed by gemini, you can:
1 - Comment /gemini review to trigger a review
2 - Comment /gemini summary for a summary
3 - Tag it using @gemini-code-assist either in the PR or in review comments on files

gemini-code-assist · 2026-01-25T13:34:48Z

Summary of Changes

Hello @DavidFang03, I'm Gemini Code Assist¹! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!

This pull request significantly expands the shammath library's capabilities for numerical analysis and data fitting. It introduces methods for solving ordinary differential equations and performing non-linear curve fitting, which are foundational for various scientific and engineering applications. The changes also include essential linear algebra components like dynamic matrix types and Cholesky decomposition, providing a more flexible and powerful mathematical toolkit.

Highlights

Least Squares Method: Implemented the Levenberg-Marquardt algorithm for non-linear least squares fitting, enabling the determination of best-fit parameters for given functions.
Cholesky Decomposition and Solver: Added Cholesky decomposition and a corresponding solver for linear systems, which are essential components for the least squares method and other numerical tasks.
Dynamic Matrix and Vector Handling: Introduced new dynamic-sized matrix (mat_d) and vector (vec_d) classes using std::vector and std::mdspan, along with a matrix transpose function, to provide more flexible and runtime-adaptable linear algebra operations.
Euler ODE Solver: Included an Euler method for numerically solving ordinary differential equations, expanding the library's capabilities for integration problems.
Python Bindings: Exposed the newly added Euler ODE solver and least squares method to Python, allowing these functionalities to be used from Python applications.

Using Gemini Code Assist

The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.

Invoking Gemini

You can request assistance from Gemini at any point by creating a comment using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands on the current page.

Feature	Command	Description
Code Review	`/gemini review`	Performs a code review for the current pull request in its current state.
Pull Request Summary	`/gemini summary`	Provides a summary of the current pull request in its current state.
Comment	@gemini-code-assist	Responds in comments when explicitly tagged, both in pull request comments and review comments.
Help	`/gemini help`	Displays a list of available commands.

Customization

To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.

Limitations & Feedback

Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here.

You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.

Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution. ↩

gemini-code-assist

Code Review

This PR introduces a significant amount of new numerical functionality, including a Levenberg-Marquardt solver, Cholesky decomposition, and an ODE integrator. The implementation is a good start, but there are several critical bugs and performance issues that need to be addressed. I've found incorrect assertions, indexing errors in matrix operations, and significant inefficiencies in the core optimization loop. The Python bindings also have performance problems due to unnecessary data copying. Additionally, some of the new tests contain bugs that could prevent them from correctly verifying the code's behavior, which aligns with the guideline to investigate the root cause in tests rather than allowing them to pass trivially. My review includes detailed suggestions to fix these issues and improve the overall quality and performance of the new code.

src/shammath/include/shammath/matrix.hpp

src/shammath/include/shammath/matrix_op.hpp

src/shammath/include/shammath/solve.hpp

src/shampylib/src/math/pyshammath.cpp

src/shammath/include/shammath/integrator.hpp

src/shammath/include/shammath/matrix.hpp

src/shammath/include/shammath/solve.hpp

src/tests/shammath/matrixTests.cpp

DavidFang03 · 2026-01-25T15:14:44Z

/gemini review

gemini-code-assist

Code Review

This pull request introduces several new mathematical functionalities, including an Euler ODE solver, dynamic-sized matrix/vector classes (mat_d, vec_d), matrix transpose, Cholesky decomposition, Cholesky solver, and a Levenberg-Marquardt least squares method. The changes also include necessary pybind11 bindings to expose these functionalities to Python. Overall, the new features are valuable additions to the library. The provided feedback addresses various areas for improvement regarding const-correctness, potential runtime issues, code clarity, and adherence to best practices for numerical stability and portability.

src/shammath/include/shammath/solve.hpp

src/shammath/include/shammath/matrix_op.hpp

src/shammath/include/shammath/matrix.hpp

gemini-code-assist · 2026-01-25T15:19:46Z

src/shammath/include/shammath/integrator.hpp

+
+        T u_prev = u0;
+        T u      = u0;
+        for (T x = x0 + step; x < end; x += step) {


Using floating-point numbers directly in loop conditions like x < end can lead to precision issues, potentially causing the loop to terminate one step early or run one step too long due to accumulated floating-point errors. It's generally safer to iterate a fixed number of times or use a small epsilon for comparison, e.g., x < end - std::numeric_limits<T>::epsilon() * std::abs(end) * N where N is the number of steps.

src/shammath/include/shammath/solve.hpp

gemini-code-assist · 2026-01-25T15:19:46Z

src/shammath/include/shammath/integrator.hpp

+        };
+        u_prev = u0;
+        std::vector<T> X_backward, U_backward;
+        for (T x = x0 - step; x > start; x -= step) {


Similar to the forward loop, using floating-point numbers directly in loop conditions like x > start can lead to precision issues. Consider using a fixed number of iterations or a small epsilon for comparison to ensure robust behavior.

…ameters

github-actions · 2026-01-27T14:48:31Z

Workflow report

workflow report corresponding to commit d2cbdd1
Commiter email is david.fang@ikmail.com

Light CI is enabled. This will only run the basic tests and not the full tests.
Merging a PR require the job "on PR / all" to pass which is disabled in this case.

Pre-commit check report

Pre-commit check: ✅

trim trailing whitespace.................................................Passed
fix end of files.........................................................Passed
check for merge conflicts................................................Passed
check that executables have shebangs.....................................Passed
check that scripts with shebangs are executable..........................Passed
check for added large files..............................................Passed
check for case conflicts.................................................Passed
check for broken symlinks................................................Passed
check yaml...............................................................Passed
detect private key.......................................................Passed
No-tabs checker..........................................................Passed
Tabs remover.............................................................Passed
Validate GitHub Workflows................................................Passed
clang-format.............................................................Passed
ruff check...............................................................Passed
ruff format..............................................................Passed
Check doxygen headers....................................................Passed
Check license headers....................................................Passed
Check #pragma once.......................................................Passed
Check SYCL #include......................................................Passed
No ssh in git submodules remote..........................................Passed
No UTF-8 in files (except for authors)...................................Passed

Test pipeline can run.

Doxygen diff with `main`

Removed warnings : 20
New warnings : 24
Warnings count : 8049 → 8053 (0.0%)

Detailed changes :

- src/shammath/include/shammath/integrator.hpp:24: warning: Member integ_riemann_sum(T start, T end, T step, Lambda &&fct) (function) of namespace shammath is not documented.
- src/shammath/include/shammath/integrator.hpp:24: warning: Member integ_riemann_sum(T start, T end, T step, Lambda &&fct) (function) of namespace shammath is not documented.
+ src/shammath/include/shammath/integrator.hpp:25: warning: Member integ_riemann_sum(T start, T end, T step, Lambda &&fct) (function) of namespace shammath is not documented.
+ src/shammath/include/shammath/integrator.hpp:25: warning: Member integ_riemann_sum(T start, T end, T step, Lambda &&fct) (function) of namespace shammath is not documented.
- src/shammath/include/shammath/matrix.hpp:126: warning: Compound sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:127: warning: Member component_type (typedef) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:128: warning: Member dimension (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:130: warning: Member is_float_based (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:132: warning: Member is_uint_based (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:135: warning: Member is_int_based (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:138: warning: Member has_info (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:140: warning: Member get_min() (function) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:144: warning: Member get_max() (function) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:148: warning: Member get_zero() (function) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:225: warning: Compound sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:226: warning: Member component_type (typedef) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:227: warning: Member dimension (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:229: warning: Member is_float_based (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:231: warning: Member is_uint_based (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:234: warning: Member is_int_based (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:237: warning: Member has_info (variable) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:239: warning: Member get_min() (function) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:243: warning: Member get_max() (function) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
+ src/shammath/include/shammath/matrix.hpp:247: warning: Member get_zero() (function) of struct sham::VectorProperties< shammath::mat< T, m, n > > is not documented.
- src/shammath/include/shammath/matrix.hpp:60: warning: Member operator+=(const mat &other) (function) of class shammath::mat is not documented.
+ src/shammath/include/shammath/matrix.hpp:62: warning: Member operator+=(const mat &other) (function) of class shammath::mat is not documented.
- src/shammath/include/shammath/solve.hpp:25: warning: Member newton_rhaphson(std::function< T(T)> &&f, std::function< T(T)> &&df, T epsilon_c, T x_0) (function) of namespace shammath is not documented.
- src/shammath/include/shammath/solve.hpp:25: warning: Member newton_rhaphson(std::function< T(T)> &&f, std::function< T(T)> &&df, T epsilon_c, T x_0) (function) of namespace shammath is not documented.
+ src/shammath/include/shammath/solve.hpp:28: warning: Member newton_rhaphson(std::function< T(T)> &&f, std::function< T(T)> &&df, T epsilon_c, T x_0) (function) of namespace shammath is not documented.
+ src/shammath/include/shammath/solve.hpp:28: warning: Member newton_rhaphson(std::function< T(T)> &&f, std::function< T(T)> &&df, T epsilon_c, T x_0) (function) of namespace shammath is not documented.
- src/shampylib/src/math/pyshammath.cpp:37: warning: Member Register_pymod(pysham_mathinit) (function) of file pyshammath.cpp is not documented.
+ src/shampylib/src/math/pyshammath.cpp:40: warning: Member Register_pymod(pysham_mathinit) (function) of file pyshammath.cpp is not documented.

DavidFang03 added 13 commits January 22, 2026 01:32

add Euler method

47bc414

add Cholesky decomposition

db367cf

add Cholesky solving

c9dd6fd

oops

04d385c

typo

c13398a

add matrix tranpose

979b3e5

add vec_update_vals

e7468f3

whoops

594076f

reshape

9487349

add dynamic matrices when dimensions are known at runtime

39d0094

add least squares

a10b2d7

authorship

b6a165a

Merge remote-tracking branch 'upstream/main' into math_for_tillotson

5e31c24

gemini-code-assist bot reviewed Jan 25, 2026

View reviewed changes

DavidFang03 added 6 commits January 25, 2026 15:46

fixes

1c42e54

fix

9ce02cd

add tolerance on sse for convergence criteria in least_squares

3e6e72f

fix

247e17d

least_squares now also returns R^2

33629bd

one last thing

002ea61

gemini-code-assist bot reviewed Jan 25, 2026

View reviewed changes

typos

ee1090f

tdavidcl added need-review light-ci labels Jan 25, 2026

DavidFang03 added 3 commits January 26, 2026 16:16

the number of observations needs to be greater than the number of par…

5e6c1ca

…ameters

doxygen

2d7614b

remove magic numbers

90d5570

DavidFang03 added 7 commits January 26, 2026 17:18

== overload constant

e4d071f

std::sqrt instead of sycl::sqrt

0ebcc10

Merge remote-tracking branch 'upstream/main' into math_for_tillotson

cd3fd51

last assert

49fdfb6

revert useless change

ef2f577

remove last doxygen warnings

087764a

Merge remote-tracking branch 'upstream/main' into math_for_tillotson

d2cbdd1

Conversation

DavidFang03 commented Jan 25, 2026

Uh oh!

github-actions bot commented Jan 25, 2026

Uh oh!

gemini-code-assist bot commented Jan 25, 2026

Summary of Changes

Highlights

Footnotes

Uh oh!

gemini-code-assist bot left a comment

Choose a reason for hiding this comment

Code Review

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

DavidFang03 commented Jan 25, 2026

Uh oh!

gemini-code-assist bot left a comment

Choose a reason for hiding this comment

Code Review

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

gemini-code-assist bot Jan 25, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

gemini-code-assist bot Jan 25, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot commented Jan 27, 2026

Workflow report

Pre-commit check report

Doxygen diff with main

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Doxygen diff with `main`