NNFS p.135

Started on the iteration, not need to rewrite matadd.h and need to make matrandom.h.

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

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"
 

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_

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
+

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
+

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

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
+

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
 

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) {}

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 ---            #
 		//##########################################################

obj/main.d

@@ -1,23 +1,24 @@
 obj/main.o: src/main.cpp include/./core/omp_config.h \
  include/./utils/utils.h include/./utils/vector.h \
- include/./utils/matrix.h include/./utils/random.h \
+ include/./utils/random.h include/./utils/matrix.h \
  include/./utils/generators.h include/./utils/generators/linspace.h \
  include/utils/vector.h include/./utils/matcast.h \
  include/./numerics/numerics.h include/./numerics/max.h \
  include/./numerics/exp.h include/./numerics/log.h \
- include/./numerics/initializers/eye.h include/./numerics/matequal.h \
- include/./numerics/abs.h include/./numerics/transpose.h \
- include/./numerics/inverse.h \
- include/./numerics/inverse/inverse_gauss_jordan.h \
- include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
- include/./numerics/matmul.h include/./numerics/matmax.h \
- include/./numerics/matdiv.h include/./numerics/matvec.h \
- include/./numerics/matadd.h include/./numerics/matsubtract.h \
- include/./numerics/matsum.h include/./numerics/matclip.h \
- include/./numerics/matexp.h include/./numerics/matlog.h \
- include/./numerics/matdot.h include/./numerics/min.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/abs.h \
+ include/./numerics/transpose.h include/./numerics/inverse.h \
+ include/./numerics/inverse/inverse_gauss_jordan.h \
+ include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
+ include/./numerics/matmul.h include/./numerics/matscalar.h \
+ include/./numerics/matmax.h include/./numerics/matdiv.h \
+ include/./numerics/matvec.h include/./numerics/matadd.h \
+ include/./numerics/matsubtract.h include/./numerics/matsum.h \
+ include/./numerics/matclip.h include/./numerics/matexp.h \
+ include/./numerics/matlog.h include/./numerics/matdot.h \
+ include/./numerics/matargmax.h include/./numerics/min.h \
  include/./numerics/interpolation1d.h \
  include/./numerics/interpolation1d/interpolation1d_barycentric.h \
  include/./numerics/interpolation1d/interpolation1d_base.h \
@@ -31,6 +32,7 @@ obj/main.o: src/main.cpp include/./core/omp_config.h \
  include/utils/matrix.h \
  include/./modules/neural_networks/neural_networks.h \
  include/./modules/neural_networks/datasets/spiral.h \
+ include/./modules/neural_networks/datasets/vertical.h \
  include/./modules/neural_networks/layers/Dense_Layer.h \
  include/./modules/neural_networks/activation_functions/Activation_ReLU.h \
  include/./modules/neural_networks/activation_functions/Activation_Softmax.h \
@@ -40,8 +42,8 @@ obj/main.o: src/main.cpp include/./core/omp_config.h \
 include/./core/omp_config.h:
 include/./utils/utils.h:
 include/./utils/vector.h:
-include/./utils/matrix.h:
 include/./utils/random.h:
+include/./utils/matrix.h:
 include/./utils/generators.h:
 include/./utils/generators/linspace.h:
 include/utils/vector.h:
@@ -50,6 +52,11 @@ include/./numerics/numerics.h:
 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/abs.h:
@@ -59,6 +66,7 @@ include/./numerics/inverse/inverse_gauss_jordan.h:
 include/./numerics/inverse/inverse_lu.h:
 include/./decomp/lu.h:
 include/./numerics/matmul.h:
+include/./numerics/matscalar.h:
 include/./numerics/matmax.h:
 include/./numerics/matdiv.h:
 include/./numerics/matvec.h:
@@ -69,11 +77,8 @@ include/./numerics/matclip.h:
 include/./numerics/matexp.h:
 include/./numerics/matlog.h:
 include/./numerics/matdot.h:
+include/./numerics/matargmax.h:
 include/./numerics/min.h:
-include/./numerics/vecclip.h:
-include/./numerics/vecexp.h:
-include/./numerics/vecmax.h:
-include/./numerics/veclog.h:
 include/./numerics/interpolation1d.h:
 include/./numerics/interpolation1d/interpolation1d_barycentric.h:
 include/./numerics/interpolation1d/interpolation1d_base.h:
@@ -90,6 +95,7 @@ include/core/global_config.h:
 include/utils/matrix.h:
 include/./modules/neural_networks/neural_networks.h:
 include/./modules/neural_networks/datasets/spiral.h:
+include/./modules/neural_networks/datasets/vertical.h:
 include/./modules/neural_networks/layers/Dense_Layer.h:
 include/./modules/neural_networks/activation_functions/Activation_ReLU.h:
 include/./modules/neural_networks/activation_functions/Activation_Softmax.h:

src/main.cpp

@@ -25,29 +25,105 @@ int main(int argc, char const *argv[])
 
     utils::Mf X(10,2, 0);
     utils::Matrix<uint64_t> y(10,1, 0);
-    neural_networks::create_spital_data<float, uint64_t>(100, 3, X, y);
+    utils::Vector<uint64_t> class_targets;
+    float loss;
+    float accuracy;
+
+    //neural_networks::create_spital_data<float, uint64_t>(10000, 3, X, y);
+    neural_networks::create_vertical_data<float, uint64_t>(100, 3, X, y);
 
     neural_networks::Dense_Layer<float> dense1(2, 3);
-    
     neural_networks::Activation_ReLU<float> activation1;
     neural_networks::Dense_Layer<float> dense2(3, 3);
     neural_networks::Activation_Softmax<float> activation2;
     
     neural_networks::Loss_CategoricalCrossentrophy<float, uint64_t> loss_funtion;
 
+    float lowest_loss = 9999999;
+    utils::Mf best_dense_1_weights = dense1.weights;
+    utils::Vf best_dense_1_biases = dense1.biases;
+    utils::Mf best_dense_2_weights = dense2.weights;
+    utils::Vf best_dense_2_biases = dense2.biases;
+    utils::Mf matrRND;
+    utils::Vf vectRND;
+
+    utils::Vector<uint64_t> predections;
+
+
+
+
+
+        for (uint64_t i = 0; i < 10000; ++i){
+
+
+            // Generate a new set of weights for iteration
+            matrRND.random(dense1.weights.rows(), dense1.weights.cols(), -1, 1);
+            numerics::inplace_matscalar(matrRND,0.05f);
+            numerics::inplace_matadd(dense1.weights, matrRND);
+
+            dense1.biases.random(dense1.biases.size(), -1, 1);
+
+            matrRND.random(dense2.weights.rows(), dense2.weights.cols(), -1, 1);
+            numerics::inplace_matscalar(matrRND,0.05f);
+            dense2.weights.random(dense2.weights.rows(), dense2.weights.cols(), -1, 1);
+            dense2.biases.random(dense2.biases.size(), -1, 1);
+
+            // Perform a forward pass of the training data through this layer
             dense1.forward(X);
             activation1.forward(dense1.outputs);
             dense2.forward(activation1.outputs);
             activation2.forward(dense2.outputs);
 
+            // Perform a farward pass through activation function
+            // it takes the output of the second dense layer here and returns loss
+            loss = loss_funtion.calculate(activation2.outputs, y);
+
+            predections = numerics::matargmax_row<uint64_t, float>(activation2.outputs);
+
+            if (y.cols() < 1){
+                class_targets = numerics::matargmax_row<uint64_t, uint64_t>(y);
+            }else{
+                class_targets = y.get_col(0);
+            }
+
+            accuracy = numerics::vecmean_equal<float>(predections, class_targets);
+
+            if (loss < lowest_loss){
+                std::cout << "New set of weights found, iteration:" << i << ", loss:" << loss << ", acc:" << accuracy << std::endl;
+                best_dense_1_weights = dense1.weights;
+                best_dense_1_biases = dense1.biases;
+                best_dense_2_weights = dense2.weights;
+                best_dense_2_biases = dense2.biases;
+                lowest_loss = loss;
+            }
+
+
+        }
+
+
+
+
+
+
+
+
     for (int i = 0; i < 5; ++i){
         std::cout << activation2.outputs.get_row(i) << std::endl;
     }
 
-    float loss = loss_funtion.calculate(activation2.outputs, y);
     std::cout << loss << std::endl;
 
     
+    //predections.print();
+    
+    
+
+
+
+    //predections.print();
+    //class_targets.print();
+
+    std::cout << accuracy << std::endl;