[Math] Add least squares method

src/shammath/include/shammath/integrator.hpp

@@ -11,6 +11,7 @@
 
 /**
  * @file integrator.hpp
+ * @author David Fang (david.fang@ikmail.com)
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
  * @brief
  *
@@ -30,4 +31,47 @@ namespace shammath {
         return acc;
     }
 
+    /**
+     * @brief Euler solving of ODE
+     * The ode has the form
+     * \f{eqnarray*}{
+     * u'(x) &=& f(u,x) \\
+     * u(x_0) &=& u_0
+     * \f}
+     * and will be solved between start and end with step \f$ \mathrm{d}t \f$.
+     *
+     * @param start Lower bound of integration
+     * @param end   Higher bound of integration
+     * @param step  Step of integration \f$ \mathrm{d}t \f$
+     * @param ode   Ode function \f$ f \f$
+     * @param x0    Initial coordinate \f$ x_0 \f$
+     * @param u0    Initial value \f$ u_0 \f$
+     */
+    template<class T, class Lambda>
+    inline constexpr std::pair<std::vector<T>, std::vector<T>> euler_ode(
+        T start, T end, T step, Lambda &&ode, T x0, T u0) {
+        std::vector<T> U = {u0};
+        std::vector<T> X = {x0};
+
+        T u_prev = u0;
+        T u      = u0;
+        for (T x = x0 + step; x < end; x += step) {
+            u = u_prev + ode(u_prev, x) * step;
+            X.push_back(x);
+            U.push_back(u);
+            u_prev = u;
+        };
+        u_prev = u0;
+        std::vector<T> X_backward, U_backward;
+        for (T x = x0 - step; x > start; x -= step) {
+            u = u_prev - ode(u_prev, x) * step;
+            X_backward.push_back(x);
+            U_backward.push_back(u);
+            u_prev = u;
+        }
+        X.insert(X.begin(), X_backward.rbegin(), X_backward.rend());
+        U.insert(U.begin(), U_backward.rbegin(), U_backward.rend());
+        return {X, U};
+    }
+
 } // namespace shammath

src/shammath/include/shammath/matrix.hpp

@@ -11,6 +11,7 @@
 
 /**
  * @file matrix.hpp
+ * @author David Fang (david.fang@ikmail.com)
  * @author Léodasce Sewanou (leodasce.sewanou@ens-lyon.fr)
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
  * @author Yann Bernard (yann.bernard@univ-grenoble-alpes.fr)
@@ -57,6 +58,7 @@ namespace shammath {
         /// Check if this matrix is equal to another one
         bool operator==(const mat<T, m, n> &other) const { return data == other.data; }
 
+        // Addition operator for matrices
         inline mat &operator+=(const mat &other) {
 #pragma unroll
             for (size_t i = 0; i < m * n; i++) {
@@ -117,9 +119,105 @@ namespace shammath {
         inline constexpr T &operator[](int i) { return get_mdspan()(i); }
 
         /// Check if this vector is equal to another one
-        bool operator==(const vec<T, n> &other) { return data == other.data; }
+        bool operator==(const vec<T, n> &other) const { return data == other.data; }
     };
 
+    /**
+     * @brief Matrix class with runtime size based on std::vector storage and mdspan
+     * @tparam T the type of the matrix entries
+     */
+    template<class T>
+    class mat_d {
+        public:
+        /// The matrix data
+        std::vector<T> data;
+        /// Number of rows
+        int rows;
+        /// Number of columns
+        int columns;
+
+        /// Constructor
+        mat_d(int rows, int columns) : rows(rows), columns(columns), data(rows * columns) {}
+
+        /// Get the matrix data as a mdspan
+        inline constexpr auto get_mdspan() {
+            return std::mdspan<T, std::dextents<size_t, 2>>(data.data(), rows, columns);
+        }
+
+        /// const overload
+        inline constexpr auto get_mdspan() const {
+            return std::mdspan<const T, std::dextents<size_t, 2>>(data.data(), rows, columns);
+        }
+
+        /// Access the matrix entry at position (i, j)
+        inline constexpr T &operator()(int i, int j) { return get_mdspan()(i, j); }
+
+        /// const overload
+        inline constexpr const T &operator()(int i, int j) const { return get_mdspan()(i, j); }
+
+        /// Check if this matrix is equal to another one
+        bool operator==(const mat_d<T> &other) const { return data == other.data; }
+
+        /// Addition operator for matrices
+        inline mat_d &operator+=(const mat_d &other) {
+            for (size_t i = 0; i < get_mdspan().extent(0) * get_mdspan().extent(1); i++) {
+                data[i] += other.data[i];
+            }
+            return *this;
+        }
+
+        /// check if this matrix is equal to another one at a given precison
+        bool equal_at_precision(const mat_d<T> &other, const T precision) const {
+            bool res = true;
+            for (auto i = 0; i < rows; i++) {
+                for (auto j = 0; j < columns; j++) {
+                    if (sham::abs(data[i * columns + j] - other.data[i * columns + j])
+                        >= precision) {
+                        res = false;
+                    }
+                }
+            }
+            return res;
+        }
+    };
+
+    /**
+     * @brief Vector class with runtime size based on std::vector storage and mdspan
+     * @tparam T the type of the vector entries
+     * @tparam n the number of entries
+     */
+    template<class T>
+    class vec_d {
+        public:
+        /// The vector data
+        std::vector<T> data;
+        /// The vector size
+        int size;
+
+        /// Constructor
+        vec_d(int size) : size(size), data(size) {}
+
+        /// Get the vector data as a mdspan
+        inline constexpr auto get_mdspan() {
+            return std::mdspan<T, std::dextents<size_t, 1>>(data.data(), size);
+        }
+
+        /// Get the vector data as a mdspan of a matrix with one column
+        inline constexpr auto get_mdspan_mat_col() {
+            return std::mdspan<T, std::dextents<size_t, 2>>(data.data(), size, 1);
+        }
+
+        /// Get the vector data as a mdspan of a matrix with one row
+        inline constexpr auto get_mdspan_mat_row() {
+            return std::mdspan<T, std::dextents<size_t, 2>>(data.data(), 1, size);
+        }
+
+        /// Access the vector entry at position i
+        inline constexpr T &operator[](int i) { return get_mdspan()(i); }
+
+        /// Check if this vector is equal to another one
+        bool operator==(const vec_d<T> &other) const { return data == other.data; }
+    };
 } // namespace shammath
 
 template<class T, int m, int n>

src/shammath/include/shammath/matrix_op.hpp

@@ -11,6 +11,7 @@
 
 /**
  * @file matrix_op.hpp
+ * @author David Fang (david.fang@ikmail.com)
  * @author Léodasce Sewanou (leodasce.sewanou@ens-lyon.fr)
  * @author Timothée David--Cléris (tim.shamrock@proton.me)
  * @author Yann Bernard (yann.bernard@univ-grenoble-alpes.fr)
@@ -52,6 +53,27 @@ namespace shammath {
         }
     }
 
+    /**
+     * @brief Set the elements of a vector according to a user-provided function
+     *
+     * @param input The vector to update the elements of
+     * @param func The function to use to update the elements of the matrix. The
+     * function must take one argument being the index.
+     *
+     * @details The function `func` is called for each element of the vector, and
+     * the value returned by the function is used to set the corresponding
+     * element of the vector.
+     */
+    template<class T, class Extents, class Layout, class Accessor, class Func>
+    inline void vec_set_vals(const std::mdspan<T, Extents, Layout, Accessor> &input, Func &&func) {
+
+        shambase::check_functor_signature<T, int>(func);
+
+        for (int i = 0; i < input.extent(0); i++) {
+            input(i) = func(i);
+        }
+    }
+
     /**
      * @brief Update the elements of a matrix according to a user-provided function
      *
@@ -653,4 +675,127 @@ namespace shammath {
         }
     }
 
+    /**
+     * @brief This function transposes a matrix
+     * @param input matrix to transpose
+     * @param output matrix to store the transposed matrix
+     */
+    template<
+        class T,
+        class Extents1,
+        class Extents2,
+        class Layout1,
+        class Layout2,
+        class Accessor1,
+        class Accessor2>
+    inline void mat_transpose(
+        const std::mdspan<T, Extents1, Layout1, Accessor1> &input,
+        const std::mdspan<T, Extents2, Layout2, Accessor2> &output) {
+
+        SHAM_ASSERT(input.extent(0) == output.extent(1));
+        SHAM_ASSERT(input.extent(1) == output.extent(0));
+
+        for (int i = 0; i < output.extent(0); i++) {
+            for (int j = 0; j < output.extent(1); j++) {
+                output(i, j) = input(j, i);
+            }
+        }
+    }
+
+    /**
+     * @brief This function performs Cholesky decomposition. From a (real) symmetric,
+     definite-positive square matrix \f$ M \f$, return a lower triangular matrix \f$ L \f$ such that
+     \f[
+        M = L L^T
+     \f]
+     * @param M a square symmetric, definite-positive matrix
+     * @param L the output matrix to store the lower triangular matrix obtained by Cholesky
+     decomposition
+     */
+    template<
+        class T,
+        class Extents1,
+        class Extents2,
+        class Layout1,
+        class Layout2,
+        class Accessor1,
+        class Accessor2>
+    inline void Cholesky_decomp(
+        const std::mdspan<T, Extents1, Layout1, Accessor1> &M,
+        const std::mdspan<T, Extents2, Layout2, Accessor2> &L) {
+
+        SHAM_ASSERT(M.extent(1) == M.extent(0));
+        SHAM_ASSERT(M.extent(0) == L.extent(0));
+        SHAM_ASSERT(L.extent(1) == L.extent(0));
+
+        for (int i = 0; i < M.extent(0); i++) {
+            T sum_ik = 0.0;
+            for (int k = 0; k < i; k++) {
+                sum_ik += L(i, k) * L(i, k);
+            }
+            L(i, i) = sycl::sqrt(M(i, i) - sum_ik);
+            for (int j = i + 1; j < M.extent(1); j++) {
+                T sum_ikjk = 0.0;
+                for (int k = 0; k < i; k++) {
+                    sum_ikjk += L(i, k) * L(j, k);
+                }
+                L(j, i) = (M(i, j) - sum_ikjk) / L(i, i);
+                L(i, j) = 0.0;
+            }
+        }
+    }
+
+    /**
+    * @brief This function solves a system of linear equations with Cholesky decomposition. The
+    system must have the form
+    \f[
+        Mx = y
+    \f]
+    where \f$ M \f$ is a (real) symmetric, definite-positive square matrix.
+    * @param M a square symmetric, definite-positive matrix
+    * @param y a vector, right hand side of the system
+    * @param x the ouput vector to store the solution of the system
+    */
+    template<
+        class T,
+        class Extents1,
+        class Extents2,
+        class Extents3,
+        class Layout1,
+        class Layout2,
+        class Layout3,
+        class Accessor1,
+        class Accessor2,
+        class Accessor3>
+    inline void Cholesky_solve(
+        const std::mdspan<T, Extents1, Layout1, Accessor1> &M,
+        const std::mdspan<T, Extents2, Layout2, Accessor2> &y,
+        const std::mdspan<T, Extents3, Layout3, Accessor3> &x) {
+
+        SHAM_ASSERT(M.extent(1) == M.extent(0));
+        SHAM_ASSERT(M.extent(1) == x.extent(0));
+        SHAM_ASSERT(M.extent(0) == y.extent(0));
+
+        std::vector<T> a(M.extent(0));
+        std::vector<T> L_storage(M.extent(0) * M.extent(1));
+
+        std::mdspan<T, Extents1> L{L_storage.data(), M.extent(0), M.extent(1)};
+        Cholesky_decomp(M, L);
+
+        for (int i = 0; i < M.extent(0); i++) {
+            T sum = 0.0;
+            for (int k = 0; k < i; k++) {
+                sum += L(i, k) * a[k];
+            }
+            a[i] = (y(i) - sum) / L(i, i);
+        }
+        for (int i = M.extent(0) - 1; i >= 0; i--) {
+            T sum = 0.0;
+            for (int k = i + 1; k < M.extent(0); k++) {
+                sum += L(k, i) * x(k);
+            }
+            x(i) = (a[i] - sum) / L(i, i);
+        }
+    }
+
 } // namespace shammath