Starting on the model.h, but need to make layer structs and structs for loss and optimizers

include/modules/neural_networks/activation_functions/Activation_Linear.h

@@ -0,0 +1,33 @@
+#pragma once
+
+#include "core/omp_config.h"
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+#include "modules/neural_networks/layers/Layer.h"
+
+
+
+namespace neural_networks{
+
+	template <typename T>
+	struct Activation_Linear : Layer<T>{
+
+		utils::Matrix<T> _inputs;
+		utils::Matrix<T> outputs;
+
+		utils::Matrix<T> dinputs;
+
+        void forward(const utils::Matrix<T>& inputs){
+        	_inputs = inputs;
+		    outputs = inputs;
+        }
+        void backward(const utils::Matrix<T>& dvalues){
+        	
+        	dinputs = dvalues;
+        }
+	};
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/activation_functions/Activation_ReLU.h

@@ -4,12 +4,12 @@
 
 #include "utils/vector.h"
 #include "utils/matrix.h"
-
+#include "modules/neural_networks/layers/Layer.h"
 
 namespace neural_networks{
 
 	template <typename T>
-	struct Activation_ReLU{
+	struct Activation_ReLU : Layer<T>{
 
 		utils::Matrix<T> _inputs;
 		utils::Matrix<T> outputs;

include/modules/neural_networks/activation_functions/Activation_Sigmoid.h

@@ -4,6 +4,7 @@
 
 #include "utils/vector.h"
 #include "utils/matrix.h"
+#include "modules/neural_networks/layers/Layer.h"
 
 #include "numerics/neg.h"
 #include "numerics/exp.h"
@@ -15,7 +16,7 @@
 namespace neural_networks{
 
 	template <typename T>
-	struct Activation_Sigmoid{
+	struct Activation_Sigmoid : Layer<T>{
 
 		utils::Matrix<T> _inputs;
 		utils::Matrix<T> outputs;

include/modules/neural_networks/activation_functions/Activation_Softmax.h

@@ -4,6 +4,7 @@
 
 #include "utils/vector.h"
 #include "utils/matrix.h"
+#include "modules/neural_networks/layers/Layer.h"
 
 #include "numerics/max.h"
 #include "numerics/sub.h"
@@ -15,7 +16,7 @@
 namespace neural_networks{
 
 	template <typename T>
-	struct Activation_Softmax{
+	struct Activation_Softmax : Layer<T>{
 
 			//utils::Matrix<T> exp_values;
 			//utils::Matrix<T> probabilities;

include/modules/neural_networks/activation_functions/Activation_Softmax_Loss_CategoricalCrossentropy.h

@@ -4,6 +4,7 @@
 
 #include "utils/vector.h"
 #include "utils/matrix.h"
+#include "modules/neural_networks/layers/Layer.h"
 
 #include "numerics/max.h"
 #include "numerics/sub.h"
@@ -17,7 +18,7 @@
 namespace neural_networks{
 
 	template <typename Td, typename Ti>
-	struct Activation_Softmax_Loss_CategoricalCrossentropy{
+	struct Activation_Softmax_Loss_CategoricalCrossentropy : Layer<Td>{
 
 			neural_networks::Activation_Softmax<Td> activation;
 			neural_networks::Loss_CategoricalCrossentrophy<Td, Ti> loss;

include/modules/neural_networks/datasets/sine.h

@@ -0,0 +1,37 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "./utils/matrix.h"
+#include "./utils/vector.h"
+#include "utils/generators/linspace.h"
+
+#include <math.h>       /* sin */
+
+//#include "./utils/random.h"
+
+//#include <math.h>
+
+
+namespace neural_networks{
+
+		template <typename T>
+		void create_sine_data(const uint64_t samples, const T length, utils::Matrix<T>& X, utils::Matrix<T>& y) {
+
+			X.resize(samples, 1);
+			X.set_col(0, utils::linspace(
+				T{0},								// start
+				length,								// stop
+				samples,							// N
+				true								// endpoint
+				));
+
+			y.resize(samples, 1);
+
+
+			for (uint64_t i = 0; i < samples; ++i){
+				y(i,0) = std::sin(X(i,0));
+			}
+		}
+
+
+} // end namesoace NN

include/modules/neural_networks/layers/Dense_Layer.h

@@ -4,13 +4,17 @@
 
 #include "utils/vector.h"
 #include "utils/matrix.h"
+#include "Layer.h"
+
+
 #include "utils/random.h"
+#include "random/random.h"
 
 
 namespace neural_networks{
 
 	template <typename T>
-	struct Dense_Layer{
+	struct Dense_Layer : Layer<T>{
 
 		T weight_regularizer_l1 = {0};
 		T weight_regularizer_l2 = {0};
@@ -48,7 +52,7 @@ namespace neural_networks{
         	this->bias_regularizer_l1 = bias_regularizer_l1;
         	this->bias_regularizer_l2 = bias_regularizer_l2;
 
-        	weights.random(n_inputs, n_neurons, -1, 1);
+        	weights = numerics::mul( rng::normal(n_inputs, n_neurons, 0.0f, 1.0f), 0.1f);
         	//weights = numerics::random_matrix(n_inputs, n_neurons, -1, 1);
         	biases.resize(n_neurons, T{0});
         	

include/modules/neural_networks/layers/Dropout_Layer.h

@@ -4,13 +4,15 @@
 
 #include "utils/vector.h"
 #include "utils/matrix.h"
+#include "modules/neural_networks/layers/Layer.h"
+
 #include "random/random.h"
 
 
 namespace neural_networks{
 
 	template <typename T>
-	struct Dropout_Layer{
+	struct Dropout_Layer : Layer<T>{
 
 		// Store rate, we invert it as for example for dropout
 		// of 0.1 we need a success rate of 0.9

include/modules/neural_networks/layers/Layer.h

@@ -0,0 +1,21 @@
+#pragma once
+
+#include "core/omp_config.h"
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+
+namespace neural_networks{
+
+
+template <typename T>
+	struct Layer {
+	    utils::Matrix<T> outputs;
+	    utils::Matrix<T> dinputs;
+
+	    virtual void forward(const utils::Matrix<T>& inputs) = 0;
+	    virtual void backward(const utils::Matrix<T>& dvalues) = 0;
+
+	    virtual ~Layer() = default;
+	};
+}

include/modules/neural_networks/loss/Loss.h

@@ -10,7 +10,7 @@
 
 namespace neural_networks{
 
-	template <typename Td, typename Ti>
+	template <typename Td, typename Ti = Td>
 	struct Loss{
 
 		utils::Vector<Td> sample_losses;

include/modules/neural_networks/loss/Loss_MeanAbsoluteError.h

@@ -0,0 +1,58 @@
+#pragma once
+
+#include "core/omp_config.h"
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+#include "utils/matcast.h"
+
+#include "numerics/clip.h"
+#include "numerics/log.h"
+#include "numerics/sub.h"
+
+#include "Loss.h"
+
+
+namespace neural_networks{
+
+	template <typename T>
+	struct Loss_MeanAbsoluteError : Loss<T>  {
+
+		utils::Matrix<T> dinputs;
+		utils::Matrix<T> y_true;
+
+		utils::Vector<T> sample_losses;
+
+	
+		utils::Vector<T> forward(const utils::Matrix<T>& y_pred, const utils::Matrix<T>& y_true) override{
+			
+			
+			// Calculate loss
+			sample_losses = numerics::mean_rowwise(numerics::abs(numerics::sub(y_true, y_pred)));
+
+
+			// Return losses
+			return sample_losses;
+		}
+		
+		void backward(const utils::Matrix<T>& dvalues, const utils::Matrix<T>& y_true) override{
+
+
+			// Number of samples
+			const T samples = static_cast<T> (dvalues.rows());
+			// Number of outputs in every sample
+			const T outputs = static_cast<T> (dvalues.cols());
+
+
+			// Gradient values
+			dinputs = numerics::div(numerics::sign(numerics::sub(y_true, dvalues)), outputs);
+			// Normalise gradient
+			dinputs = numerics::div(dinputs, samples);
+
+
+		}
+	};
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/loss/Loss_MeanSquaredError.h

@@ -0,0 +1,58 @@
+#pragma once
+
+#include "core/omp_config.h"
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+#include "utils/matcast.h"
+
+#include "numerics/clip.h"
+#include "numerics/log.h"
+#include "numerics/sub.h"
+
+#include "Loss.h"
+
+
+namespace neural_networks{
+
+	template <typename T>
+	struct Loss_MeanSquaredError : Loss<T>  {
+
+		utils::Matrix<T> dinputs;
+		utils::Matrix<T> y_true;
+
+		utils::Vector<T> sample_losses;
+
+	
+		utils::Vector<T> forward(const utils::Matrix<T>& y_pred, const utils::Matrix<T>& y_true) override{
+			
+			
+			// Calculate loss
+			sample_losses = numerics::mean_rowwise(numerics::pow(numerics::sub(y_true, y_pred), T{2}));
+
+
+			// Return losses
+			return sample_losses;
+		}
+		
+		void backward(const utils::Matrix<T>& dvalues, const utils::Matrix<T>& y_true) override{
+
+
+			// Number of samples
+			const T samples = static_cast<T> (dvalues.rows());
+			// Number of outputs in every sample
+			const T outputs = static_cast<T> (dvalues.cols());
+
+
+			// Gradient values
+			dinputs = numerics::mul(numerics::div(numerics::sub(y_true, dvalues), outputs), T{-2});
+			// Normalise gradient
+			dinputs = numerics::div(dinputs, samples);
+
+
+		}
+	};
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/model/model.h

@@ -0,0 +1,30 @@
+#pragma once
+
+#include "core/omp_config.h"
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+
+#include "modules/neural_networks/layers/Layer.h"
+
+
+namespace neural_networks {
+
+    template <typename T>
+    struct Model {
+
+        std::vector<Layer<T>*> layers;
+
+        Model() = default;
+
+        void add(Layer<T>& layer) {
+            layers.push_back(&layer);
+        }
+
+    };
+
+ 
+
+
+
+} // end namespace neural_networks

include/modules/neural_networks/neural_networks.h

@@ -3,8 +3,9 @@
 
 #include "datasets/spiral.h"
 #include "datasets/vertical.h"
+#include "datasets/sine.h"
 
-
+#include "layers/Layer.h"
 #include "layers/Dense_Layer.h"
 #include "layers/Dropout_Layer.h"
 
@@ -13,14 +14,21 @@
 #include "activation_functions/Activation_Softmax.h"
 #include "activation_functions/Activation_Softmax_Loss_CategoricalCrossentropy.h"
 #include "activation_functions/Activation_Sigmoid.h"
+#include "activation_functions/Activation_Linear.h"
 
 
 #include "loss/Loss.h" // Base
 #include "loss/Loss_CategoricalCrossentrophy.h"
 #include "loss/Loss_BinaryCrossentropy.h"
+#include "loss/Loss_MeanSquaredError.h"
+
+
 
 
 #include "optimizers/Optimizer_SGD.h"
 #include "optimizers/Optimizer_Adagrad.h"
 #include "optimizers/Optimizer_RMSprop.h"
-#include "optimizers/Optimizer_Adam.h"
+#include "optimizers/Optimizer_Adam.h"
+
+
+#include "model/model.h"

include/numerics/detail/greater_or_less_serial.h

@@ -0,0 +1,622 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+
+
+namespace numerics::detail{
+	
+	
+	//
+	// ---------------- Greater - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_greater_serial(utils::Vector<T>& v, const T c) {
+		const uint64_t n = v.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_greater_serial: empty vector");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] < c){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		const uint64_t n = v.size();
+		const uint64_t m = p.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_greater_serial: empty vector");
+		}
+
+		if (n != m){
+			throw std::runtime_error("inplace_greater_serial: vector size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] < p[i]){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	//
+	// ---------------- Greater - Matrices ----------------
+	//
+	template <typename T>
+	inline void inplace_greater_serial(utils::Matrix<T>& A, const T c) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_serial: empty matrix");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) < c){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_rowwise_serial: empty matrix");
+		}
+
+		if (rows != v.size()){
+			throw std::runtime_error("inplace_greater_rowwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) < v[i]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_colwise_serial: empty matrix");
+		}
+
+		if (cols != v.size()){
+			throw std::runtime_error("inplace_greater_colwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) < v[j]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_serial: empty matrix");
+		}
+
+		if (rows != B.rows() || cols != B.cols()){
+			throw std::runtime_error("inplace_greater_serial: Matrix size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) < B(i,j)){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+
+
+
+
+
+
+
+	//
+	// ---------------- Greater or Equal - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_greater_equal_serial(utils::Vector<T>& v, const T c) {
+		const uint64_t n = v.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_greater_equal_serial: empty vector");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] <= c){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_equal_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		const uint64_t n = v.size();
+		const uint64_t m = p.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_greater_equal_serial: empty vector");
+		}
+
+		if (n != m){
+			throw std::runtime_error("inplace_greater_equal_serial: vector size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] <= p[i]){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	//
+	// ---------------- Greater - Matrices ----------------
+	//
+	template <typename T>
+	inline void inplace_greater_equal_serial(utils::Matrix<T>& A, const T c) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_equal_serial: empty matrix");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) <= c){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_equal_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_equal_rowwise_serial: empty matrix");
+		}
+
+		if (rows != v.size()){
+			throw std::runtime_error("inplace_greater_equal_rowwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) <= v[i]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_equal_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_equal_colwise_serial: empty matrix");
+		}
+
+		if (cols != v.size()){
+			throw std::runtime_error("inplace_greater_equal_colwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) <= v[j]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_greater_equal_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_greater_equal_serial: empty matrix");
+		}
+
+		if (rows != B.rows() || cols != B.cols()){
+			throw std::runtime_error("inplace_greater_equal_serial: Matrix size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) <= B(i,j)){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+	//
+	// ---------------- Less - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_less_serial(utils::Vector<T>& v, const T c) {
+		const uint64_t n = v.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_less_serial: empty vector");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] > c){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		const uint64_t n = v.size();
+		const uint64_t m = p.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_less_serial: empty vector");
+		}
+
+		if (n != m){
+			throw std::runtime_error("inplace_less_serial: vector size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] > p[i]){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	//
+	// ---------------- Greater - Matrices ----------------
+	//
+	template <typename T>
+	inline void inplace_less_serial(utils::Matrix<T>& A, const T c) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_serial: empty matrix");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) > c){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_rowwise_serial: empty matrix");
+		}
+
+		if (rows != v.size()){
+			throw std::runtime_error("inplace_less_rowwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) > v[i]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_colwise_serial: empty matrix");
+		}
+
+		if (cols != v.size()){
+			throw std::runtime_error("inplace_less_colwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) > v[j]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_serial: empty matrix");
+		}
+
+		if (rows != B.rows() || cols != B.cols()){
+			throw std::runtime_error("inplace_less_serial: Matrix size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) > B(i,j)){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+
+
+
+
+
+
+
+	//
+	// ---------------- Less or Equal - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_less_equal_serial(utils::Vector<T>& v, const T c) {
+		const uint64_t n = v.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_less_equal_serial: empty vector");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] >= c){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_equal_serial(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		const uint64_t n = v.size();
+		const uint64_t m = p.size();
+
+		if (n == 0){
+			throw std::runtime_error("inplace_less_equal_serial: empty vector");
+		}
+
+		if (n != m){
+			throw std::runtime_error("inplace_less_equal_serial: vector size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			if (v[i] >= p[i]){
+				v[i] = T{0};
+			}
+			else{
+				v[i] = T{1};
+			}
+		}
+	}
+	//
+	// ---------------- Greater - Matrices ----------------
+	//
+	template <typename T>
+	inline void inplace_less_equal_serial(utils::Matrix<T>& A, const T c) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_equal_serial: empty matrix");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) >= c){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_equal_rowwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_equal_rowwise_serial: empty matrix");
+		}
+
+		if (rows != v.size()){
+			throw std::runtime_error("inplace_less_equal_rowwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) >= v[i]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_equal_colwise_serial(utils::Matrix<T>& A, const utils::Vector<T>& v) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_equal_colwise_serial: empty matrix");
+		}
+
+		if (cols != v.size()){
+			throw std::runtime_error("inplace_less_equal_colwise_serial: Matrix-Vector size mismatch");
+		}
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) >= v[j]){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+	template <typename T>
+	inline void inplace_less_equal_serial(utils::Matrix<T>& A, const utils::Matrix<T>& B) {
+		const uint64_t rows = A.rows();
+		const uint64_t cols = A.cols();
+
+
+		if (rows == 0 || cols == 0){
+			throw std::runtime_error("inplace_less_equal_serial: empty matrix");
+		}
+
+		if (rows != B.rows() || cols != B.cols()){
+			throw std::runtime_error("inplace_less_equal_serial: Matrix size mismatch");
+		}
+
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				if (A(i,j) >= B(i,j)){
+					A(i,j) = T{0};
+				}
+				else{
+					A(i,j) = T{1};
+				}
+			}
+		}
+	}
+
+
+
+
+} // namespace numerics
+

include/numerics/detail/sign_serial.h

@@ -0,0 +1,48 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <cmath> // std::abs
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+
+namespace numerics::detail{
+	
+
+
+	template <typename T>
+	void inplace_sign_serial(T& c) {
+		if (c < T{0}){
+			c = T{-1};
+		}
+		else if(c > T{0}){
+			c = T{1};
+		}
+		else{
+			c = T{0};
+		}
+	}
+		
+	template <typename T>
+	void inplace_sign_serial(utils::Vector<T>& v) {
+		const uint64_t n = v.size();
+		for (uint64_t i = 0; i < n; ++i){
+			inplace_sign_serial(v[i]);
+		}
+	}
+
+
+	template <typename T>
+	void inplace_sign_serial(utils::Matrix<T>& A) {
+	    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){
+				inplace_sign_serial(A(i,j));
+			}
+		}
+	}
+
+
+} // namespace numerics
+

include/numerics/detail/standard_deviation_serial.h

@@ -0,0 +1,61 @@
+#pragma once
+
+#include <cstdint> //uint64_t
+#include <cmath> // std::abs
+
+#include "utils/vector.h"
+#include "utils/matrix.h"
+
+#include "numerics/mean.h"
+
+namespace numerics::detail{
+	
+		
+	template <typename T>
+	inline T standard_deviation_serial(const utils::Vector<T>& v) {
+		const uint64_t n = v.size();
+
+		if (n == 0){
+			throw std::runtime_error("std_serial: empty vector");
+		}
+
+		const T mean_value = numerics::mean(v);
+
+		T sum_squared_diff = T{0};
+
+
+		for (uint64_t i = 0; i < n; ++i){
+			const T diff = v[i] - mean_value;
+			sum_squared_diff += diff * diff;
+		}
+
+		const T variance = sum_squared_diff / static_cast<T>(n);
+
+		return numerics::sqrt(variance); 
+	}
+
+
+	template <typename T>
+	inline T standard_deviation_serial(const utils::Matrix<T>& A) {
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    const T mean_value = numerics::mean(A);
+
+	    T sum_squared_diff =  T{0};
+
+
+		for (uint64_t i = 0; i < rows; ++i){
+			for (uint64_t j = 0; j < cols; ++j){
+				const T diff = A(i,j) - mean_value;
+				sum_squared_diff += diff * diff;
+			}
+		}
+		const T variance = sum_squared_diff / static_cast<T>(rows*cols);
+
+		return numerics::sqrt(variance);
+	}
+
+
+} // namespace numerics
+

include/numerics/greater_or_less.h

@@ -0,0 +1,316 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/greater_or_less_serial.h"
+
+
+namespace numerics{
+
+
+	//
+	// ---------------- Greater - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_greater(utils::Vector<T>& v, const T c) {
+		detail::inplace_greater_serial(v, c);
+	}
+	template <typename T>
+	inline utils::Vector<T> greater(const utils::Vector<T>& v, const T c) {
+		utils::Vector<T> out = v;
+		inplace_greater(out, c);
+		return out;
+	}
+
+
+	template <typename T>
+	inline void inplace_greater(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_greater_serial(v, p);
+	}
+	template <typename T>
+	inline utils::Vector<T> greater(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+		utils::Vector<T> out = v;
+		inplace_greater(out, p);
+		return out;
+	}
+
+
+	//
+	// ---------------- Greater - Matrices  ----------------
+	//
+	template <typename T>
+	inline void inplace_greater(utils::Matrix<T>& A, const T c) {
+		detail::inplace_greater_serial(A, c);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater(const utils::Matrix<T>& A, const T c) {
+		utils::Matrix<T> out = A;
+		inplace_greater(out, c);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_greater_rowwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_greater_rowwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater_rowwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_greater_rowwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_greater_colwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_greater_colwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater_colwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_greater_colwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_greater(utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		detail::inplace_greater_serial(A, B);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater(const utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		utils::Matrix<T> out = A;
+		inplace_greater(out, B);
+		return out;
+	}
+
+
+	//
+	// ---------------- Greater or Equal - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_greater_equal(utils::Vector<T>& v, const T c) {
+		detail::inplace_greater_equal_serial(v, c);
+	}
+	template <typename T>
+	inline utils::Vector<T> greater_equal(const utils::Vector<T>& v, const T c) {
+		utils::Vector<T> out = v;
+		inplace_greater_equal(out, c);
+		return out;
+	}
+
+
+	template <typename T>
+	inline void inplace_greater_equal(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_greater_equal_serial(v, p);
+	}
+	template <typename T>
+	inline utils::Vector<T> greater_equal(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+		utils::Vector<T> out = v;
+		inplace_greater_equal(out, p);
+		return out;
+	}
+
+
+	//
+	// ---------------- Greater or Equal - Matrices  ----------------
+	//
+	template <typename T>
+	inline void inplace_greater_equal(utils::Matrix<T>& A, const T c) {
+		detail::inplace_greater_equal_serial(A, c);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater_equal(const utils::Matrix<T>& A, const T c) {
+		utils::Matrix<T> out = A;
+		inplace_greater_equal(out, c);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_greater_equal_rowwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_greater_equal_rowwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater_equal_rowwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_greater_equal_rowwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_greater_equal_colwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_greater_equal_colwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater_equal_colwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_greater_equal_colwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_greater_equal(utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		detail::inplace_greater_equal_serial(A, B);
+	}
+	template <typename T>
+	inline utils::Matrix<T> greater_euqal(const utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		utils::Matrix<T> out = A;
+		inplace_greater_equal(out, B);
+		return out;
+	}
+
+
+	//
+	// ---------------- Less - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_less(utils::Vector<T>& v, const T c) {
+		detail::inplace_less_serial(v, c);
+	}
+	template <typename T>
+	inline utils::Vector<T> less(const utils::Vector<T>& v, const T c) {
+		utils::Vector<T> out = v;
+		inplace_less(out, c);
+		return out;
+	}
+
+
+	template <typename T>
+	inline void inplace_less(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_less_serial(v, p);
+	}
+	template <typename T>
+	inline utils::Vector<T> less(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+		utils::Vector<T> out = v;
+		inplace_less(out, p);
+		return out;
+	}
+
+
+
+
+	//
+	// ---------------- Less - Matrices  ----------------
+	//
+	template <typename T>
+	inline void inplace_less(utils::Matrix<T>& A, const T c) {
+		detail::inplace_less_serial(A, c);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less(const utils::Matrix<T>& A, const T c) {
+		utils::Matrix<T> out = A;
+		inplace_less(out, c);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_less_rowwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_less_rowwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less_rowwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_less_rowwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_less_colwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_less_colwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less_colwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_less_colwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_less(utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		detail::inplace_less_serial(A, B);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less(const utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		utils::Matrix<T> out = A;
+		inplace_less(out, B);
+		return out;
+	}
+
+
+	//
+	// ---------------- Less or Equal - Vectors ----------------
+	//
+	template <typename T>
+	inline void inplace_less_equal(utils::Vector<T>& v, const T c) {
+		detail::inplace_less_equal_serial(v, c);
+	}
+	template <typename T>
+	inline utils::Vector<T> less_equal(const utils::Vector<T>& v, const T c) {
+		utils::Vector<T> out = v;
+		inplace_less_equal(out, c);
+		return out;
+	}
+
+
+	template <typename T>
+	inline void inplace_less_equal(utils::Vector<T>& v, const utils::Vector<T>& p) {
+		detail::inplace_less_equal_serial(v, p);
+	}
+	template <typename T>
+	inline utils::Vector<T> less_equal(const utils::Vector<T>& v, const utils::Vector<T>& p) {
+		utils::Vector<T> out = v;
+		inplace_less_equal(out, p);
+		return out;
+	}
+
+
+	//
+	// ---------------- Less or Equal - Matrices  ----------------
+	//
+	template <typename T>
+	inline void inplace_less_equal(utils::Matrix<T>& A, const T c) {
+		detail::inplace_less_equal_serial(A, c);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less_equal(const utils::Matrix<T>& A, const T c) {
+		utils::Matrix<T> out = A;
+		inplace_less_equal(out, c);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_less_equal_rowwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_less_equal_rowwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less_equal_rowwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_less_equal_rowwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_less_equal_colwise(utils::Matrix<T>& A, const utils::Vector<T> v) {
+		detail::inplace_less_equal_colwise_serial(A, v);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less_equal_colwise(const utils::Matrix<T>& A, const utils::Vector<T> v) {
+		utils::Matrix<T> out = A;
+		inplace_less_equal_colwise(out, v);
+		return out;
+	}
+
+	template <typename T>
+	inline void inplace_less_equal(utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		detail::inplace_less_equal_serial(A, B);
+	}
+	template <typename T>
+	inline utils::Matrix<T> less_euqal(const utils::Matrix<T>& A, const utils::Matrix<T> B) {
+		utils::Matrix<T> out = A;
+		inplace_less_equal(out, B);
+		return out;
+	}
+
+
+
+
+
+
+}

include/numerics/numerics.h

@@ -10,6 +10,7 @@
 #include "./numerics/dot.h"
 #include "./numerics/equal.h"
 #include "./numerics/exp.h"
+#include "./numerics/greater_or_less.h"
 #include "./numerics/interpolation1d.h"                       		// base
 #include "./numerics/inverse.h"
 #include "./numerics/isclose.h"
@@ -22,7 +23,9 @@
 #include "./numerics/neg.h"
 #include "./numerics/pow.h"
 #include "./numerics/random.h"
+#include "./numerics/sign.h"
 #include "./numerics/sqrt.h"
+#include "./numerics/standard_deviation.h"
 #include "./numerics/sub.h"
 #include "./numerics/sum.h"
 #include "./numerics/transpose.h"

include/numerics/sign.h

@@ -0,0 +1,51 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/sign_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline void inplace_sign(T& c) {
+		detail::inplace_sign_serial(c);
+	}
+	template <typename T>
+	inline T sign(const T c) {
+		T out = c;
+		inplace_sign(out);
+		return out;
+	}
+
+
+
+	template <typename T>
+	inline void inplace_sign(utils::Vector<T>& v) {
+		detail::inplace_sign_serial(v);
+	}
+
+	template <typename T>
+	inline utils::Vector<T> sign(const utils::Vector<T>& v) {
+	    utils::Vector<T> out = v;
+	    inplace_sign(out);
+	    return out;
+	}
+
+
+
+
+	template <typename T>
+	inline void inplace_sign(utils::Matrix<T>& A) {
+		detail::inplace_sign_serial(A);
+	}
+
+	template <typename T>
+	inline utils::Matrix<T> sign(const utils::Matrix<T>& A) {
+	    utils::Matrix<T> out = A;
+	    inplace_sign(out);
+	    return out;
+	}
+
+
+}

include/numerics/standard_deviation.h

@@ -0,0 +1,22 @@
+#pragma once
+
+#include "./core/omp_config.h"
+#include "detail/standard_deviation_serial.h"
+
+
+namespace numerics{
+
+	// ---------------- Elementwise ----------------
+	template <typename T>
+	inline T standard_deviation(const utils::Vector<T>& v) {
+		return detail::standard_deviation_serial(v);
+	}
+
+	template <typename T>
+	inline T standard_deviation(const utils::Matrix<T>& A){
+		return detail::standard_deviation_serial(A);
+	}
+
+
+
+}

include/random/detail/normal_serial.h

@@ -0,0 +1,57 @@
+#pragma once
+
+#include "random/engine.h"
+
+#include <random>
+//#include <type_traits>
+
+#include "utils/matrix.h"
+#include "utils/vector.h"
+
+namespace rng::detail {
+
+	//
+	// Base functions int
+	//
+    template <typename T>
+    typename std::enable_if<std::is_integral<T>::value, T>::type
+    normal_serial(const T mean, const T stddev) {
+        std::normal_distribution<T> dist(mean, stddev);
+        return dist(rng::engine());
+    }
+
+	//
+	// Base functions float
+	//
+    template <typename T>
+    typename std::enable_if<std::is_floating_point<T>::value, T>::type
+    normal_serial(const T mean, const T stddev) {
+        std::normal_distribution<T> dist(mean, stddev);
+        return dist(rng::engine());
+    }
+
+    template <typename T>
+    utils::Vector<T> normal_serial(const uint64_t size, const T mean, const T stddev){
+    	utils::Vector<T> v(size);
+    	for (uint64_t i = 0; i < size; ++i){
+    		v[i] = normal_serial(mean, stddev);
+    	}
+    	return v;
+    }
+
+    template <typename T>
+    utils::Matrix<T> normal_serial(const uint64_t rows, const uint64_t cols, const T mean, const T stddev){
+    	utils::Matrix<T> A(rows, cols);
+    	for (uint64_t i = 0; i < rows; ++i){
+    		for (uint64_t j = 0; j < cols; ++j){
+    			A(i,j) = normal_serial(mean, stddev);
+    		}
+    	}
+    	return A;
+    }
+
+
+
+
+
+}

include/random/normal.h

@@ -0,0 +1,28 @@
+#pragma once
+
+
+#include "detail/normal_serial.h"
+
+namespace rng {
+
+    template <typename T>
+    T normal(const T mean, const T stddev) {
+        return rng::detail::normal_serial(mean, stddev);
+    }
+
+    template <typename T>
+    utils::Vector<T> normal(const uint64_t size, const T mean, const T stddev) {
+        return rng::detail::normal_serial(size, mean, stddev);
+    }
+
+    template <typename T>
+    utils::Matrix<T> normal(const uint64_t rows, const uint64_t cols, const T mean, const T stddev) {
+        return rng::detail::normal_serial(rows, cols, mean, stddev);
+    }
+
+
+
+
+}
+
+

include/random/random.h

@@ -7,5 +7,6 @@
 
 #include "uniform.h"
 #include "binomial.h"
+#include "normal.h"