From 1c93024ff24309eb930a2519a2984c8d62201913 Mon Sep 17 00:00:00 2001
From: Gitea CI <ci@bausager.org>
Date: Wed, 8 Oct 2025 16:08:04 +0000
Subject: [PATCH] Sync public subset from Flux

---
 .../neural_networks/datasets/vertical.h       |  40 +++++
 .../modules/neural_networks/neural_networks.h |   1 +
 include/numerics/matadd.h                     | 160 ++++--------------
 include/numerics/matargmax.h                  |  62 +++++++
 include/numerics/matscalar.h                  |  32 ++++
 include/numerics/numerics.h                   |  12 +-
 include/numerics/vecargmax.h                  |  26 +++
 include/numerics/vecmean.h                    |  38 +++++
 include/utils/matrix.h                        |   2 +-
 include/utils/vector.h                        |  20 ++-
 10 files changed, 260 insertions(+), 133 deletions(-)
 create mode 100644 include/modules/neural_networks/datasets/vertical.h
 create mode 100644 include/numerics/matargmax.h
 create mode 100644 include/numerics/matscalar.h
 create mode 100644 include/numerics/vecargmax.h

diff --git a/include/modules/neural_networks/datasets/vertical.h b/include/modules/neural_networks/datasets/vertical.h
new file mode 100644
index 0000000..a5c1dc8
--- /dev/null
+++ b/include/modules/neural_networks/datasets/vertical.h
@@ -0,0 +1,40 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "./utils/matrix.h"
+#include "./utils/vector.h"
+#include "./utils/random.h"
+
+//#include <math.h>
+
+
+namespace neural_networks{
+
+		template <typename TX, typename Ty>
+		void create_vertical_data(const uint64_t samples, const uint64_t classes, utils::Matrix<TX>& X, utils::Matrix<Ty>& y) {
+
+			const uint64_t rows = samples*classes;
+			uint64_t row_idx;
+
+
+			if ((rows != X.rows()) || (X.cols() != 2)){
+				X.resize(samples*classes, 2);
+			}
+			if (rows != y.rows()){
+				y.resize(rows, 1);
+			}
+
+			for (uint64_t i = 0; i < classes; ++i){
+				for (uint64_t j = 0; j < samples; ++j){
+
+					row_idx = (i*samples) + j;
+
+					X(row_idx, 0) = static_cast<TX>(i)/static_cast<TX>(classes) + utils::random(TX{-0.1}, TX{0.1});
+					X(row_idx, 1) = TX{0.5} + utils::random(TX{-0.5}, TX{0.5});
+					y(row_idx, 0) = static_cast<Ty>(i);
+				}
+			}
+		}
+
+
+} // end namesoace NN
\ No newline at end of file
diff --git a/include/modules/neural_networks/neural_networks.h b/include/modules/neural_networks/neural_networks.h
index 3db2b2b..85523ed 100644
--- a/include/modules/neural_networks/neural_networks.h
+++ b/include/modules/neural_networks/neural_networks.h
@@ -2,6 +2,7 @@
 #pragma once
 
 #include "datasets/spiral.h"
+#include "datasets/vertical.h"
 
 #include "layers/Dense_Layer.h"
 
diff --git a/include/numerics/matadd.h b/include/numerics/matadd.h
index d3890c8..624e3c3 100644
--- a/include/numerics/matadd.h
+++ b/include/numerics/matadd.h
@@ -7,6 +7,39 @@
 
 namespace numerics{
 	
+	template <typename T>
+	void inplace_matadd_mat(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_matadd: dimension mismatch");
+	    }
+
+	    for (uint64_t i = 0; i < cols; ++i) {
+	        for (uint64_t j = 0; j < rows; ++j) {
+	            A(j, i) +=  B(j, i);
+	        }
+	    }
+	}
+
+	template <typename T>
+	utils::Matrix<T> matadd_mat(const utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    utils::Matrix<T> C = A;
+
+	    if (rows != B.rows() || cols != B.cols()) {
+	        throw std::runtime_error("inplace_matadd: dimension mismatch");
+	    }
+
+	    inplace_matadd_mat(C, B);
+	    return C;
+	}
+
 	template <typename T>
 	void inplace_matadd_colvec(utils::Matrix<T>& A, const utils::Vector<T>& x) {
 	    
@@ -94,133 +127,6 @@ namespace numerics{
 		}
     }
 
-
-
-
-	/*
-	// -------------- Collapse(2) OpenMP ----------------
-	template <typename T>
-	utils::Vector<T> matvec_omp(const utils::Matrix<T>& A, const utils::Vector<T>& x) {
-	    if (A.cols() != x.size()) {
-	        throw std::runtime_error("matvec: dimension mismatch");
-	    }
-
-	    const uint64_t m = A.rows();
-	    const uint64_t n = A.cols();
-
-	    utils::Vector<T> y(m, T{0}); // <-- y has length m (rows)
-
-
-        const T* xptr = x.data();
-        const T* Aptr = A.data();	// row-major: A(i,j) == Aptr[i*n + j]
-
-        // Each row i is an independent dot product: y[i] = dot(A[i,*], x)
-	    #pragma omp parallel for schedule(static)
-	    for (uint64_t i = 0; i < m; ++i) {
-            const T* row = Aptr + i * n;     // contiguous row i
-            T acc = T{0};
-            #pragma omp simd reduction(+:acc)
-	        for (uint64_t j = 0; j < n; ++j) {
-	            acc += row[j] * xptr[j];
-	        }
-	        y[i] = acc;
-		}
-
-	    return y;
-	}
-
-	// -------------- Auto OpenMP ----------------
-	template <typename T>
-	utils::Vector<T> matvec_auto(const utils::Matrix<T>& A,
-	                             const utils::Vector<T>& x) {
-
-	    
-	    uint64_t work = A.rows() * A.cols();
-
-	    bool can_parallel = omp_config::omp_parallel_allowed();
-		#ifdef _OPENMP
-		    int threads = omp_get_max_threads();
-		#else
-		    int threads = 1;
-		#endif
-
-	    if (can_parallel || work > static_cast<uint64_t>(threads) * 4ull) {
-	        return matvec_omp(A,x);
-	    }
-	    else{
-	    	// Safe fallback
-	    	return matvec(A,x);
-	    }
-	    
-	}
-
-// =================================================
-//   y = x * A    (Vector–Matrix product)
-// =================================================
-	template <typename T>
-	utils::Vector<T> vecmat(const utils::Vector<T>& x, const utils::Matrix<T>& A) {
-	    if (x.size() != A.rows()) {
-	        throw std::runtime_error("vecmat: dimension mismatch");
-	    }
-	    const uint64_t m = A.rows();
-	    const uint64_t n = A.cols();
-
-	    utils::Vector<T> y(n, T{0});
-
-	    for (uint64_t j = 0; j < n; ++j) {
-	        for (uint64_t i = 0; i < m; ++i) {
-	            y[j] += x[i] * A(i, j);
-	        }
-	    }
-	    return y;
-	}
-
-	// -------------- Collapse(2) OpenMP ----------------
-	template <typename T>
-	utils::Vector<T> vecmat_omp(const utils::Vector<T>& x, const utils::Matrix<T>& A) {
-	    if (x.size() != A.rows()) {
-	        throw std::runtime_error("vecmat: dimension mismatch");
-	    }
-	    const uint64_t m = A.rows();
-	    const uint64_t n = A.cols();
-
-	    utils::Vector<T> y(n, T{0});
-	    #pragma omp parallel for schedule(static)
-	    for (uint64_t j = 0; j < n; ++j) {
-	        T acc = T{0};
-	        for (uint64_t i = 0; i < m; ++i) {
-	            acc += x[i] * A(i, j);
-	        }
-	        y[j] = acc;
-	    }
-	    return y;
-	}
-
-	// -------------- Auto OpenMP ----------------
-	template <typename T>
-	utils::Vector<T> vecmat_auto(const utils::Vector<T>& x,
-								 const utils::Matrix<T>& A) {
-	    
-	    uint64_t work = A.rows() * A.cols();
-
-	    bool can_parallel = omp_config::omp_parallel_allowed();
-		#ifdef _OPENMP
-		    int threads = omp_get_max_threads();
-		#else
-		    int threads = 1;
-		#endif
-
-	    if (can_parallel || work > static_cast<uint64_t>(threads) * 4ull) {
-	        return vecmat_omp(x,A);
-	    }
-	    else{
-	    	// Safe fallback
-	    	return vecmat(x,A);
-	    }
-	    
-	}
-*/
-
 } // namespace numerics
 
 #endif // _matadd_n_
\ No newline at end of file
diff --git a/include/numerics/matargmax.h b/include/numerics/matargmax.h
new file mode 100644
index 0000000..7fd411f
--- /dev/null
+++ b/include/numerics/matargmax.h
@@ -0,0 +1,62 @@
+#pragma once
+
+#include "./utils/matrix.h"
+
+namespace numerics{
+	
+	template <typename Ti, typename Td>
+	void inplace_matargmax_row(const utils::Matrix<Td>& A, utils::Vector<Ti>& b){
+
+		if (b.size() != A.rows()){
+			b.resize(A.rows(), Ti{0});
+		}
+		Td value;
+		
+		for (uint64_t i = 0; i < A.rows(); ++i){
+			value = Td{0};
+			for (uint64_t j = 0; j < A.cols(); ++j){
+				if (value < A(i,j)){
+					value = A(i,j);
+					b[i] = j;
+				}
+			}
+		}
+	}
+
+	template <typename Ti, typename Td>
+	void inplace_matargmax_col(const utils::Matrix<Td>& A, utils::Vector<Ti>& b){
+
+		if (b.size() != A.cols()){
+			b(A.cols(), Ti{0});
+		}
+		Td value;
+		
+		for (uint64_t j = 0; j < A.cols(); ++j){
+			value = Td{0};
+			for (uint64_t i = 0; i < A.cols(); ++i){
+				if (value < A(i,j)){
+					value = A(i,j);
+					b[j] = i;
+				}
+			}
+		}
+	}
+
+	template <typename Ti, typename Td>
+	utils::Vector<Ti> matargmax_row(const utils::Matrix<Td>& A){
+		utils::Vector<Ti> b(A.rows(), Ti{0});
+		inplace_matargmax_row(A, b);
+		return b;
+
+	}
+
+	template <typename Ti, typename Td>
+	utils::Vector<Ti> matargmax_col(const utils::Matrix<Td>& A){
+		utils::Vector<Ti> b(A.rows(), Ti{0});
+		inplace_matargmax_col(A, b);
+		return b;
+
+	}
+
+} // namespace numerics
+
diff --git a/include/numerics/matscalar.h b/include/numerics/matscalar.h
new file mode 100644
index 0000000..a4c2dfb
--- /dev/null
+++ b/include/numerics/matscalar.h
@@ -0,0 +1,32 @@
+#pragma once
+
+#include "./utils/matrix.h"
+
+namespace numerics{
+	
+	template <typename T>
+	void inplace_matscalar(utils::Matrix<T>& A, const T scalar){
+
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				A(i,j) *= scalar;
+			}
+		}
+	}
+
+	template <typename T>
+	utils::Matrix<T> matscalar(const utils::Matrix<T>& A, T scalar){
+
+		utils::Matrix<T> B = A;
+		inplace_matscalar(B, scalar);
+		return B;
+
+	}
+
+
+
+} // namespace numerics
+
diff --git a/include/numerics/numerics.h b/include/numerics/numerics.h
index d048dab..27f2a10 100644
--- a/include/numerics/numerics.h
+++ b/include/numerics/numerics.h
@@ -3,11 +3,17 @@
 #include "./numerics/max.h"
 #include "./numerics/exp.h"
 #include "./numerics/log.h"
+#include "./numerics/vecclip.h"
+#include "./numerics/vecexp.h"
+#include "./numerics/vecmax.h"
+#include "./numerics/veclog.h"
+#include "./numerics/vecargmax.h"
 #include "./numerics/initializers/eye.h"
 #include "./numerics/matequal.h"
 #include "./numerics/transpose.h"
 #include "./numerics/inverse.h"
 #include "./numerics/matmul.h"
+#include "./numerics/matscalar.h"
 #include "./numerics/matmax.h"
 #include "./numerics/matdiv.h"
 #include "./numerics/matvec.h"
@@ -18,12 +24,10 @@
 #include "./numerics/matexp.h"
 #include "./numerics/matlog.h"
 #include "./numerics/matdot.h"
+#include "./numerics/matargmax.h"
 #include "./numerics/min.h"
 #include "./numerics/abs.h"
-#include "./numerics/vecclip.h"
-#include "./numerics/vecexp.h"
-#include "./numerics/vecmax.h"
-#include "./numerics/veclog.h"
+
 
 
 #include "./numerics/interpolation1d.h"                       		// base
diff --git a/include/numerics/vecargmax.h b/include/numerics/vecargmax.h
new file mode 100644
index 0000000..f61b99d
--- /dev/null
+++ b/include/numerics/vecargmax.h
@@ -0,0 +1,26 @@
+#pragma once
+
+#include "./utils/vector.h"
+
+namespace numerics{
+	
+	template <typename T>
+	T vecarmax(const utils::Vector<T>& a){
+
+		const uint64_t N = a.size();
+		uint64_t idx = 0;
+		T value;
+
+		for (uint64_t i = 0; i < N; ++i){
+			if (value > a[i]){
+				value = a[i];
+				idx = i;
+
+			}
+		}
+		return idx;
+	}
+
+
+} // namespace numerics
+
diff --git a/include/numerics/vecmean.h b/include/numerics/vecmean.h
index 268b5aa..65b0b7a 100644
--- a/include/numerics/vecmean.h
+++ b/include/numerics/vecmean.h
@@ -22,6 +22,44 @@ namespace numerics{
 	    return mean;
 	}
 
+	template <typename T>
+	T vecmean_isequal(utils::Vector<T>& a, utils::Vector<T>& b, T tol=1e-12) {
+
+		uint64_t count = T{0};
+
+		const uint64_t N = a.size();
+
+		if (a.size() != b.size()){
+			throw std::runtime_error("vecmean_equal: dimension mismatch");
+		}
+
+
+	    for (uint64_t i = 0; i < N; ++i) {
+	    	if ((a[i]-b[i]) < tol){
+	    		count += 1;
+	    	}
+	    }
+	    return static_cast<T>(count)/static_cast<T>(N);
+	}
+	template <typename Td, typename Ti>
+	Td vecmean_equal(utils::Vector<Ti>& a, utils::Vector<Ti>& b) {
+
+		Ti count = Ti{0};
+
+		const uint64_t N = a.size();
+
+		if (a.size() != b.size()){
+			throw std::runtime_error("vecmean_equal: dimension mismatch");
+		}
+
+
+	    for (uint64_t i = 0; i < N; ++i) {
+	    	if (a[i]==b[i]){
+	    		count += Ti{1};
+	    	}
+	    }
+	    return static_cast<Td>(count)/static_cast<Td>(N);
+	}
 
 } // namespace numerics
 
diff --git a/include/utils/matrix.h b/include/utils/matrix.h
index c13cb3d..8929e72 100644
--- a/include/utils/matrix.h
+++ b/include/utils/matrix.h
@@ -23,7 +23,7 @@ namespace utils{
 	class Matrix{
 public:
 	    Matrix() : rows_(0), cols_(0), data_() {} // Default constructor
-
+#include <random>
 	   // Constructor to initialize matrix with rows × cols and a fill value
 	    Matrix(uint64_t rows, uint64_t cols, const T& value = T()) 
 	    : rows_(rows), cols_(cols), data_(rows * cols, value) {}
diff --git a/include/utils/vector.h b/include/utils/vector.h
index 5780b19..03c03ce 100644
--- a/include/utils/vector.h
+++ b/include/utils/vector.h
@@ -7,7 +7,7 @@
 
 #include <initializer_list>  
 
-
+#include "./utils/random.h"
 #include <cstdint>
 #include <type_traits>
 #include <stdexcept>
@@ -38,6 +38,24 @@ public:
 
 
 
+
+
+		void random(const uint64_t size, const T& lower, const T& higher){
+
+			v.resize(size, T{0});
+
+		    // Copy data row by row
+		    for (uint64_t i = 0; i < size; ++i) {
+	            v[i] = utils::random(lower, higher);
+		    }
+
+		}
+
+
+
+
+
+
 	    //##########################################################
 		//#            VECTOR: --- basic properties ---            #
 		//##########################################################