Optimazation is done

include/math/sqrt.hpp

@@ -0,0 +1,163 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       math
+ * File Name:         sqrt.cpp
+ * Revision:          0.1.0
+ * Date:              06-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Functions to calculate the sqrt of numbers
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+#pragma once
+
+#include <tensor/vector.hpp> // for panic::vector
+#include <tensor/matrix.hpp> // for panic::matrix
+
+namespace panic{
+namespace math{
+
+
+
+/**
+ * @brief calculates the sqrt of a value.
+ *
+ * Computes:
+ * @code
+ * sqrt(x, y)
+ * @endcode
+ *
+ * @tparam T Numeric element type.
+ * @param x Value to take the sqrt of.
+ * @param y Result.
+ *
+ *
+ * @note This function is omp-friendly.
+ */
+template <typename T>
+bool sqrt(const T x, T& y);
+
+
+/**
+ * @brief calculates the sqrt of a value.
+ *
+ * Computes:
+ * @code
+ * result = sqrt(k)
+ * @endcode
+ *
+ * @tparam T Numeric element type.
+ * @param x Value to take the sqrt of.
+ *
+ * @return The calculated value
+ *
+ * @note This function is omp-friendly.
+ */
+template <typename T>
+T sqrt(const T x);
+
+
+/**
+ * @brief Calculates the sqrt elementwise in a vector
+ *
+ * Computes:
+ * @code
+ * c[i] = sqrt(a[i])
+ * @endcode
+ *
+ * @tparam T Numeric element type.
+ * @param a Input vector.
+ * @param c Output vector. Resized to match @p a.
+ *
+ * @return true if @p c was resized and filled successfully.
+ * @return false if resizing @p c failed.
+ *
+ * @note This overload writes the result into an existing vector to avoid
+ *       unnecessary temporary allocations.
+ */
+template <typename T>
+bool sqrt(const panic::tensor::vector<T>& a, panic::tensor::vector<T>& c);
+
+/**
+ * @brief Calculates the sqrt elementwise in a vector
+ *
+ * Computes:
+ * @code
+ * result[i] = sqrt(a[i])
+ * @endcode
+ *
+ * @tparam T Numeric element type.
+ * @param a Input vector.
+ *
+ * @return A new vector containing the result.
+ * @return An empty vector if the operation fails.
+ *
+ * @note This overload is convenient, but may allocate a new vector.
+ */
+template <typename T>
+panic::tensor::vector<T> sqrt(const panic::tensor::vector<T>& a);
+
+/**
+ * @brief Calculates the sqrt elementwise of a matrix
+ *
+ * Computes:
+ * @code
+ * C(i,j) = sqrt(A(i,j))
+ * @endcode
+ *
+ * @tparam T Numeric element type.
+ * @param A Input matrix.
+ * @param C Output Matrix. Resized to match @p A.
+ *
+ * @return true if @p C was resized and filled successfully.
+ * @return false if resizing @p C failed.
+ *
+ * @note This overload writes the result into an existing vector to avoid
+ *       unnecessary temporary allocations.
+ */
+template <typename T>
+bool sqrt(const panic::tensor::matrix<T>& A, panic::tensor::matrix<T>& C);
+
+/**
+ * @brief Returns the calculated sqrt elementwise of the matrix
+ *
+ * Computes:
+ * @code
+ * result(i,j) = sqrt(A(i,j))
+ * @endcode
+ *
+ * @tparam T Numeric element type.
+ * @param A Input matrix.
+ *
+ * @return A new matrix containing the result.
+ * @return An empty matrix if the operation fails.
+ *
+ * @note This overload is convenient, but may allocate a new vector.
+ */
+template <typename T>
+panic::tensor::matrix<T> sqrt(const panic::tensor::matrix<T>& A);
+
+} // namespace math
+} // namespace panic

include/neural_network/layer/trainable_layer.hpp

@@ -63,6 +63,24 @@ struct trainable_layer:layer{
     panic::tensor::real_matrix dweights;
     panic::tensor::real_vector dbiases;
 
+    /**
+     * @brief Previous parameter updates used by momentum SGD.
+     *
+     * These remain empty unless an optimizer using momentum
+     * initializes them.
+     */
+    panic::tensor::real_matrix weight_momentums;
+    panic::tensor::real_vector bias_momentums;
+
+    /**
+     * @brief Previous parameter updates used by AdaGrad.
+     *
+     * These remain empty unless an optimizer using cache
+     * initializes them.
+     */
+    panic::tensor::real_matrix weight_cache;
+    panic::tensor::real_vector bias_cache;
+
 
 	virtual ~trainable_layer() = default;
 

include/neural_network/model/model.hpp

@@ -323,9 +323,70 @@ struct model{
      *  
      *
      */
-    bool add_optimizer_sgd(const panic::types::real_t learning_rate = static_cast<panic::types::real_t>(1));
+    bool add_optimizer_sgd(const panic::types::real_t learning_rate = 1, 
+                            const panic::types::real_t decay = 0,
+                            const panic::types::real_t momentum = 0);
 
 
+    /**
+     * @brief Adds optimizer_adagrad to the model
+     *
+     * Computes:
+     * @code
+     * model.optimizer_adagrad(1e-4)
+     * @endcode
+     *
+     * @param learning_rate Learning rate for update_param (default 1e-3).
+     *
+     * @return true If looped and optimized every trainable layer.
+     *  
+     *
+     */
+    bool add_optimizer_adagrad(const panic::types::real_t learning_rate = 1, 
+                            const panic::types::real_t decay = 0,
+                            const panic::types::real_t epsilon = 1e-7);
+
+
+    /**
+     * @brief Adds optimizer_rmsprop to the model
+     *
+     * Computes:
+     * @code
+     * model.optimizer_adagrad(1e-4)
+     * @endcode
+     *
+     * @param learning_rate Learning rate for update_param (default 1e-3).
+     *
+     * @return true If looped and optimized every trainable layer.
+     *  
+     *
+     */
+    bool add_optimizer_rmsprop(const panic::types::real_t learning_rate = 0.001, 
+                                const panic::types::real_t decay = 0,
+                                const panic::types::real_t epsilon = 1e-7,
+                                const panic::types::real_t rho = 0.9);
+
+
+
+    /**
+     * @brief Adds optimizer_adam to the model
+     *
+     * Computes:
+     * @code
+     * model.optimizer_adam(1e-4)
+     * @endcode
+     *
+     * @param learning_rate Learning rate for update_param (default 1e-3).
+     *
+     * @return true If looped and optimized every trainable layer.
+     *  
+     *
+     */
+    bool add_optimizer_adam(const panic::types::real_t learning_rate = 0.001, 
+                                const panic::types::real_t decay = 0,
+                                const panic::types::real_t epsilon = 1e-7,
+                                const panic::types::real_t beta_1 = 0.9,
+                                const panic::types::real_t beta_2 = 0.999);
 
     /**
      * @brief Finalizes the model configuration.

include/neural_network/optimizers/optimizer.hpp

@@ -56,6 +56,8 @@ namespace panic{
  */
 struct optimizer{
 
+     panic::types::real_t current_learning_rate;
+
     /**
      * @brief Default de-constructor
      *
@@ -64,15 +66,33 @@ struct optimizer{
 
 
     /**
-     * @brief Virtual forward function for derivative layers
+     * @brief Virtual update parameters function for derivative optimizers
      *
-     * @param inputs Data matrix input for forward function.
+     * @param layer trainable layer to have their parameters updated
      * 
      * @Note It's equal to 0 because it make the derivative 
      *         object NEEDS to have these function to work.
      */
      virtual bool update_params(trainable_layer& layer) = 0;
 
+    /**
+     * @brief Virtual function to update internal parameters before update_params()
+     *
+     * 
+     */
+     virtual bool pre_update_params(){
+          return true;
+     }
+
+    /**
+     * @brief Virtual function to update internal parameters after update_params()
+     *
+     * 
+     */
+     virtual bool post_update_params(){
+          return true;
+     }
+
 };
 
 

include/neural_network/optimizers/optimizer_adagrad.hpp

@@ -0,0 +1,112 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       neural_network
+ * File Name:         optimizer_adagrad.hpp
+ * Revision:          0.1.0
+ * Date:              04-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Defines the optimizer_adagrad struct used in neural network
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+#pragma once
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//---------------------------------------------------------------------------------------------------------------------------
+#include <config/types.hpp>
+#include <neural_network/optimizers/optimizer.hpp>
+
+
+namespace panic{
+     namespace neural_network{
+
+
+
+/**
+ * @brief optimizer_adagrad for the rest of the neural network library to use
+ *
+ */
+struct optimizer_adagrad: optimizer{
+
+     panic::types::real_t learning_rate;
+
+     panic::types::real_t decay;
+
+     panic::types::real_t epsilon;
+
+     panic::types::uint_t iterations;
+
+
+
+    /**
+     * @brief Constructor
+     *
+     * @param learning_rate The learning rate for the optimization.
+     * @param decay The decay for the learning rate over the interations.
+     * 
+     */
+     optimizer_adagrad(const panic::types::real_t learning_rate = static_cast<panic::types::real_t>(1),
+                    const panic::types::real_t decay = static_cast<panic::types::real_t>(0),
+                    const panic::types::real_t epsilon = static_cast<panic::types::real_t>(1e-7));
+
+
+    /**
+     * @brief Default de-constructor
+     *
+     */
+     ~optimizer_adagrad() = default;
+
+
+    /**
+     * @brief Updates weights and biases in trainable  layers
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+     bool update_params(trainable_layer& layer) override;
+
+
+
+    /**
+     * @brief function to update internal parameters before update_params()
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+      bool pre_update_params() override;
+
+    /**
+     * @brief function to update internal parameters after update_params()
+     * 
+     */
+     bool post_update_params() override;
+
+};
+
+     } // namespace tensor
+} // namespace panic
+
+
+

include/neural_network/optimizers/optimizer_adam.hpp

@@ -0,0 +1,120 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       neural_network
+ * File Name:         optimizer_rmsprop.hpp
+ * Revision:          0.1.0
+ * Date:              04-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Defines the optimizer_adam struct used in neural network
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+#pragma once
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//---------------------------------------------------------------------------------------------------------------------------
+#include <config/types.hpp>
+#include <neural_network/optimizers/optimizer.hpp>
+
+
+namespace panic{
+     namespace neural_network{
+
+
+
+/**
+ * @brief optimizer_adam for the rest of the neural network library to use
+ *
+ */
+struct optimizer_adam: optimizer{
+
+     panic::types::real_t learning_rate;
+
+     panic::types::real_t decay;
+
+     panic::types::real_t epsilon;
+
+     panic::types::real_t beta_1;
+     panic::types::real_t beta_2;
+
+     panic::types::real_t beta_1_power;
+     panic::types::real_t beta_2_power;
+
+     panic::types::uint_t iterations;
+
+
+
+    /**
+     * @brief Constructor
+     *
+     * @param learning_rate The learning rate for the optimization.
+     * @param decay The decay for the learning rate over the interations.
+     * 
+     */
+     optimizer_adam(const panic::types::real_t learning_rate = static_cast<panic::types::real_t>(0.001),
+                    const panic::types::real_t decay = static_cast<panic::types::real_t>(0),
+                    const panic::types::real_t epsilon = static_cast<panic::types::real_t>(1e-7),
+                    const panic::types::real_t beta_1 = static_cast<panic::types::real_t>(0.9),
+                    const panic::types::real_t beta_2 = static_cast<panic::types::real_t>(0.999));
+
+
+    /**
+     * @brief Default de-constructor
+     *
+     */
+     ~optimizer_adam() = default;
+
+
+    /**
+     * @brief Updates weights and biases in trainable  layers
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+     bool update_params(trainable_layer& layer) override;
+
+
+
+    /**
+     * @brief function to update internal parameters before update_params()
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+      bool pre_update_params() override;
+
+    /**
+     * @brief function to update internal parameters after update_params()
+     * 
+     */
+     bool post_update_params() override;
+
+};
+
+     } // namespace tensor
+} // namespace panic
+
+
+

include/neural_network/optimizers/optimizer_rmsprop.hpp

@@ -0,0 +1,115 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       neural_network
+ * File Name:         optimizer_rmsprop.hpp
+ * Revision:          0.1.0
+ * Date:              04-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Defines the optimizer_rmsprop struct used in neural network
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+#pragma once
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//---------------------------------------------------------------------------------------------------------------------------
+#include <config/types.hpp>
+#include <neural_network/optimizers/optimizer.hpp>
+
+
+namespace panic{
+     namespace neural_network{
+
+
+
+/**
+ * @brief optimizer_rmsprop for the rest of the neural network library to use
+ *
+ */
+struct optimizer_rmsprop: optimizer{
+
+     panic::types::real_t learning_rate;
+
+     panic::types::real_t decay;
+
+     panic::types::real_t epsilon;
+
+     panic::types::real_t rho;
+
+     panic::types::uint_t iterations;
+
+
+
+    /**
+     * @brief Constructor
+     *
+     * @param learning_rate The learning rate for the optimization.
+     * @param decay The decay for the learning rate over the interations.
+     * 
+     */
+     optimizer_rmsprop(const panic::types::real_t learning_rate = static_cast<panic::types::real_t>(0.001),
+                    const panic::types::real_t decay = static_cast<panic::types::real_t>(0),
+                    const panic::types::real_t epsilon = static_cast<panic::types::real_t>(1e-7),
+                    const panic::types::real_t rho = static_cast<panic::types::real_t>(0.9));
+
+
+    /**
+     * @brief Default de-constructor
+     *
+     */
+     ~optimizer_rmsprop() = default;
+
+
+    /**
+     * @brief Updates weights and biases in trainable  layers
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+     bool update_params(trainable_layer& layer) override;
+
+
+
+    /**
+     * @brief function to update internal parameters before update_params()
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+      bool pre_update_params() override;
+
+    /**
+     * @brief function to update internal parameters after update_params()
+     * 
+     */
+     bool post_update_params() override;
+
+};
+
+     } // namespace tensor
+} // namespace panic
+
+
+

include/neural_network/optimizers/optimizer_sgd.hpp

@@ -52,13 +52,24 @@ struct optimizer_sgd: optimizer{
 
      panic::types::real_t learning_rate;
 
+     panic::types::real_t decay;
+
+     panic::types::real_t momentum;
+
+     panic::types::uint_t iterations;
+
+
+
     /**
      * @brief Constructor
      *
      * @param learning_rate The learning rate for the optimization.
+     * @param decay The decay for the learning rate over the interations.
      * 
      */
-     optimizer_sgd(const panic::types::real_t learning_rate = static_cast<panic::types::real_t>(1e-3));
+     optimizer_sgd(const panic::types::real_t learning_rate = static_cast<panic::types::real_t>(1),
+                    const panic::types::real_t decay = static_cast<panic::types::real_t>(0),
+                    const panic::types::real_t momentum = static_cast<panic::types::real_t>(0));
 
 
     /**
@@ -76,10 +87,24 @@ struct optimizer_sgd: optimizer{
      */
      bool update_params(trainable_layer& layer) override;
 
+
+
+    /**
+     * @brief function to update internal parameters before update_params()
+     *
+     * @param layer Trianable layer to update.
+     * 
+     */
+      bool pre_update_params() override;
+
+    /**
+     * @brief function to update internal parameters after update_params()
+     * 
+     */
+     bool post_update_params() override;
+
 };
 
-
-
      } // namespace tensor
 } // namespace panic
 

main.cpp

@@ -64,6 +64,7 @@
 #include <math/equal.hpp>
 #include <math/mean.hpp>
 #include <neural_network/activation_loss/activation_softmax_loss_categorical_crossentropy.hpp>
+#include <math/sqrt.hpp>
 
 #include <math.h>
 
@@ -709,8 +710,9 @@ int main(void) {
     panic::types::uint_t classes = 3;
 
 
-    // create spiral data
+    // create spiral data 
     panic::neural_network::spiral_data(samples, classes, X, y);
+    //panic::neural_network::vertical_data(samples, classes, X, y);
 
     // Initilise my model
     panic::neural_network::model mymodel;
@@ -741,10 +743,42 @@ int main(void) {
         return false;
     }
 
-    if (! mymodel.add_optimizer_sgd()){
+    /*
+    panic::types::real_t learning_rate = 1;
+    panic::types::real_t decay = 1e-3;
+    panic::types::real_t momentum = 0.9;
+    if (! mymodel.add_optimizer_sgd(learning_rate, decay, momentum)){
         return false;
     }
-   
+    */
+    /*
+    panic::types::real_t learning_rate = 1;
+    panic::types::real_t decay = 1e-4;
+    panic::types::real_t epsilon = 1e-7;
+    if (! mymodel.add_optimizer_adagrad(learning_rate, decay, epsilon)){
+        return false;
+    }
+    */
+    /*
+    panic::types::real_t learning_rate = 0.001;
+    panic::types::real_t decay = 1e-4;
+    panic::types::real_t epsilon = 1e-7;
+    panic::types::real_t rho = 0.999;
+    if (! mymodel.add_optimizer_rmsprop(learning_rate, decay, epsilon, rho)){
+        return false;
+    }
+    */
+
+    panic::types::real_t learning_rate = 0.02;
+    panic::types::real_t decay = 1e-5;
+    panic::types::real_t epsilon = 1e-7;
+    panic::types::real_t beta_1 = 0.9;
+    panic::types::real_t beta_2 = 0.999;
+    
+    if (! mymodel.add_optimizer_adam(learning_rate, decay, epsilon, beta_1, beta_2)){
+        return false;
+    }
+
     if (!mymodel.finalize()){
         return false;
     }
@@ -752,15 +786,13 @@ int main(void) {
     
     
     panic::types::uint_t epochs = 10000;
-    panic::types::uint_t print_every = 100;
+    panic::types::uint_t print_every = 250;
 
     if (!mymodel.train(X, y, epochs, print_every)){
         std::cout << "Training failed" << std::endl;
         return false;
     }
 
-    
-
 
 	return 0;
 }

src/math/sqrt.cpp

@@ -0,0 +1,322 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       math
+ * File Name:         sqrt.cpp
+ * Revision:          0.1.0
+ * Date:              06-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Functions to calculate the sqrt of numbers
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//-----------------------------------------------------------------------------------------------------
+
+#include <math/sqrt.hpp>
+#include <config/omp.hpp>
+#include <config/types.hpp>
+
+#include <tensor/vector.hpp> // for panic::vector
+#include <tensor/matrix.hpp> // for panic::matrix
+
+//---------------------------------------------------------------------------------------------------------------------------
+// PRIVATE CONSTANTS
+//---------------------------------------------------------------------------------------------------------------------------
+/**
+ * @brief Minimum number of element operations before using the OpenMP-enabled loop.
+ *
+ * Small vectors and matrices are kept serial because the overhead of starting
+ * worker threads can be larger than the work itself.
+ */
+static const panic::types::uint_t sqrt_omp_min_work = 250;
+//---------------------------------------------------------------------------------------------------------------------------
+// INPLEMENTATION
+//---------------------------------------------------------------------------------------------------------------------------
+
+
+namespace panic {
+	namespace math {
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::math::sqrt
+//
+// Description:
+//      Calculates the square root of a non-negative real number
+//      using Heron's method.
+//
+//      Returns 0 when:
+//          - input is negative
+//          - input is NaN
+//          - input is positive infinity
+//          - the method does not converge within the iteration limit
+//--------------------------------------------------------------------------------------------------------------------------
+template <typename T>
+bool sqrt(const T x, T& y){
+
+    const T zero = static_cast<T>(0);
+
+    const T one = static_cast<T>(1);
+
+    const T half = static_cast<T>(0.5);
+
+    const T tolerance = static_cast<T>(1e-6);
+
+    const panic::types::uint_t max_iterations = static_cast<panic::types::uint_t>(128);
+
+    // Reject negative numbers and NaN.
+    // NaN >= 0 evaluates to false.
+    if (!(x >= zero)){
+        return false;
+    }
+
+    // sqrt(0) is exactly 0.
+    // This special case also avoids division by zero.
+    if (x == zero){
+        y = T{0};
+        return true;
+    }
+
+    // Choose an initial estimate that works reasonably for
+    // values both above and below 1.
+    T current;
+
+    if (x >= one){
+        current = x;
+    }
+    else {
+        current = one;
+    }
+
+    for (panic::types::uint_t iteration = static_cast<panic::types::uint_t>(0); iteration < max_iterations; ++iteration){
+
+        const T next = half *(current + x / current);
+
+        T difference = next - current;
+
+        if (difference < zero){
+            difference = - difference;
+        }
+
+        /*
+         * next == current means floating-point rounding has
+         * prevented any further improvement.
+         */
+        if (next == current || difference <= tolerance * next){
+            y = next;
+            return true;
+        }
+        y = next;
+        current = next;
+    }
+
+
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// EXPLICIT TEMPLATE INSTANTIATION
+//
+// The implementation is in this .cpp file.
+// Build the overload for the official PANIC numeric types.
+//--------------------------------------------------------------------------------------------------------------------------
+template bool sqrt(const panic::types::real_t x,
+                    panic::types::real_t& y
+);
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::math::sqrt
+//
+// Description:
+//      Calculates the square root of a non-negative real number
+//      using Heron's method.
+//
+//      Returns false when:
+//          - input is negative
+//          - input is NaN
+//          - input is positive infinity
+//          - the method does not converge within the iteration limit
+//--------------------------------------------------------------------------------------------------------------------------
+template <typename T>
+T sqrt(const T x){
+    T y;
+    if (! sqrt(x,y)){
+        return T();
+    }
+    
+    return y;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// EXPLICIT TEMPLATE INSTANTIATION
+//
+// The implementation is in this .cpp file.
+// Build the overload for the official PANIC numeric types.
+//--------------------------------------------------------------------------------------------------------------------------
+template panic::types::real_t sqrt(const panic::types::real_t x
+);
+
+
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::math::sqrt
+//
+// Description:
+//		Calculates the sqrt elementwise of a vector
+//--------------------------------------------------------------------------------------------------------------------------
+template <typename T>
+bool sqrt(const panic::tensor::vector<T>& a, panic::tensor::vector<T>& c){
+
+	if (!c.resize(a.size())){
+		return false;
+	}
+
+    bool valid = true;
+
+    // valid remains true only if every sqrt is valid.
+    PANIC_OMP_PARALLEL_FOR_REDUCTION_IF( a.size() > sqrt_omp_min_work, &&, valid )
+    for (panic::types::uint_t i = 0; i < a.size(); ++i){
+        const bool nonzero = sqrt(a[i], c[i]);
+        valid = valid && nonzero;
+    }
+
+
+	return valid;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// EXPLICIT TEMPLATE INSTANTIATION
+//
+// The implementation is in this .cpp file.
+// Build the overload for the official PANIC numeric types.
+//--------------------------------------------------------------------------------------------------------------------------
+template bool sqrt<panic::types::real_t>(const panic::tensor::vector<panic::types::real_t>& a,
+										panic::tensor::vector<panic::types::real_t>& c
+);
+
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::math::sqrt
+//
+// Description:
+//		Calculates the sqrt elementwise for a vector
+//--------------------------------------------------------------------------------------------------------------------------
+template <typename T>
+panic::tensor::vector<T> sqrt(const panic::tensor::vector<T>& a){
+	panic::tensor::vector<T> c(a.size());
+
+	if (!sqrt(a, c)){
+		return panic::tensor::vector<T>();
+	}
+
+	return c;
+}
+//--------------------------------------------------------------------------------------------------------------------------
+// EXPLICIT TEMPLATE INSTANTIATION
+//
+// The implementation is in this .cpp file.
+// Build the overload for the official PANIC numeric types.
+//--------------------------------------------------------------------------------------------------------------------------
+template panic::tensor::vector<panic::types::real_t> 
+				sqrt(const panic::tensor::vector<panic::types::real_t>& a
+);
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::math::sqrt
+//
+// Description:
+//		calculates the natrual sqrt elementwise of a matrix
+//--------------------------------------------------------------------------------------------------------------------------
+template <typename T>
+bool sqrt(const panic::tensor::matrix<T>& A, panic::tensor::matrix<T>& C){
+
+	panic::types::uint_t rows = A.rows();
+	panic::types::uint_t cols = A.cols();
+	panic::types::uint_t work = rows*cols;
+
+	if ( !C.resize(rows, cols) ){
+		return false;
+	}
+
+    bool valid = true;
+
+    // valid remains true only if every sqrt is valid.
+	PANIC_OMP_PARALLEL_FOR_IF(work > sqrt_omp_min_work)
+	for (panic::types::uint_t i = 0; i < rows; ++i){
+		for (panic::types::uint_t j = 0; j < cols; ++j){
+            const bool nonzero = sqrt(A(i,j), C(i,j));
+            valid = valid && nonzero;
+		}
+	}
+
+	return valid;
+}
+//--------------------------------------------------------------------------------------------------------------------------
+// EXPLICIT TEMPLATE INSTANTIATION
+//
+// The implementation is in this .cpp file.
+// Build the overload for the official PANIC numeric types.
+//--------------------------------------------------------------------------------------------------------------------------
+template bool sqrt(const panic::tensor::matrix<panic::types::real_t>& A,
+					panic::tensor::matrix<panic::types::real_t>& C
+);
+
+
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::math::sqrt
+//
+// Description:
+//		Calculates the natrual sqrt element-wise of a matrix
+//--------------------------------------------------------------------------------------------------------------------------
+template <typename T>
+panic::tensor::matrix<T> sqrt(const panic::tensor::matrix<T>& A){
+	panic::tensor::matrix<T> C;
+
+	if (!sqrt(A, C)){
+		return panic::tensor::matrix<T>();
+	}
+
+	return C;
+}
+//--------------------------------------------------------------------------------------------------------------------------
+// EXPLICIT TEMPLATE INSTANTIATION
+//
+// The implementation is in this .cpp file.
+// Build the overload for the official PANIC numeric types.
+//--------------------------------------------------------------------------------------------------------------------------
+template panic::tensor::matrix<panic::types::real_t> 
+						sqrt(const panic::tensor::matrix<panic::types::real_t>& A
+);
+
+
+	} // namespace math
+} // namespace panic

src/math/sum.cpp

@@ -264,7 +264,7 @@ bool sum_colwise(const panic::tensor::matrix<T>& A, panic::tensor::vector<T>& b)
 	for (panic::types::uint_t i = 0; i < cols; ++i){
 		b[i] = T{0};
 		for (panic::types::uint_t j = 0; j < rows; ++j){
-			b[i] += A(i,j);
+			b[i] += A(j,i);
 		}
 		
 	}

src/neural_network/activation/activation_relu.cpp

@@ -51,9 +51,10 @@
  * Small vectors and matrices are kept serial because the overhead of starting
  * worker threads can be larger than the work itself.
  */
-static const panic::types::uint_t activation_ReLU_omp_min_size = 500;
+static const panic::types::uint_t activation_ReLU_omp_min_size = 250;
 
 static const panic::types::real_t almost_zero = static_cast<panic::types::real_t> (1e-7);
+static const panic::types::real_t zero = static_cast<panic::types::real_t> (0);
 //---------------------------------------------------------------------------------------------------------------------------
 // INPLEMENTATION
 //---------------------------------------------------------------------------------------------------------------------------
@@ -82,7 +83,7 @@ bool activation_relu::forward(const panic::tensor::real_matrix& input_data){
 
 	inputs = input_data;
 
-	panic::math::clip_lower(inputs, almost_zero, outputs);
+	panic::math::clip_lower(inputs, zero, outputs);
 
 
 	return true;
@@ -98,8 +99,28 @@ bool activation_relu::forward(const panic::tensor::real_matrix& input_data){
 //--------------------------------------------------------------------------------------------------------------------------
 bool activation_relu::backward(const panic::tensor::real_matrix& dvalues){
 
-	// Zero gradients where input values were negative
+    if (dvalues.rows() != inputs.rows() || dvalues.cols() != inputs.cols()){
-	dinputs = panic::math::clip_lower(dvalues, almost_zero);
+        return false;
+    }
+
+    dinputs = dvalues;
+
+    const panic::types::uint_t rows = dinputs.rows();
+    const panic::types::uint_t cols = dinputs.cols();
+    const panic::types::uint_t work = rows * cols;
+
+    PANIC_OMP_PARALLEL_FOR_IF(work > activation_ReLU_omp_min_size)
+    for (panic::types::uint_t i = 0; i < rows; ++i){
+        for (panic::types::uint_t j = 0; j < cols; ++j){
+
+            // ReLU output did not depend on this input,
+            // so no gradient passes backward.
+            if (inputs(i, j) <= static_cast<panic::types::real_t>(0)){
+                
+                dinputs(i, j) = static_cast<panic::types::real_t>(0);
+            }
+        }
+    }
 
 
 	return true;

src/neural_network/datasets/spiral_data.cpp

@@ -96,7 +96,7 @@ bool spiral_data(const panic::types::uint_t samples, const panic::types::uint_t
 	}
 
 	panic::tensor::matrix<T> random_matrix(samples*classes, 2);
-	panic::random::uniform(random_matrix, T{-0.15}, T{0.15});
+	panic::random::uniform(random_matrix, T{-0.0015}, T{0.0015});
 	
 	if (!panic::math::add(X, random_matrix, X)){
 		return false;

src/neural_network/layer/layer_dense.cpp

@@ -85,7 +85,7 @@ layer_dense::layer_dense() {
 //--------------------------------------------------------------------------------------------------------------------------
 layer_dense::layer_dense(panic::types::uint_t input_size, panic::types::uint_t neurons) {
 	weights.resize(input_size, neurons);
-	panic::random::uniform(weights);
+	panic::random::uniform(weights, static_cast<panic::types::real_t>(-1), static_cast<panic::types::real_t>(1));
 	panic::math::mul(weights, static_cast<panic::types::real_t>(0.01), weights);
 
 	biases.resize(neurons);

src/neural_network/model/model.cpp

@@ -49,7 +49,10 @@
 
 #include <neural_network/activation_loss/activation_softmax_loss_categorical_crossentropy.hpp>
 
-#include <neural_network/optimizers/optimizer_sgd.hpp>    
+#include <neural_network/optimizers/optimizer_sgd.hpp>
+#include <neural_network/optimizers/optimizer_adagrad.hpp>
+#include <neural_network/optimizers/optimizer_rmsprop.hpp>
+#include <neural_network/optimizers/optimizer_adam.hpp>
 
 #include <math/argmax.hpp>
 #include <math/equal.hpp>
@@ -545,9 +548,86 @@ bool model::add_loss_categorical_crossentropy(){
 // Example:
 //      model.add_optimizer_sgd(1e-4);
 //--------------------------------------------------------------------------------------------------------------------------
-bool model::add_optimizer_sgd(const panic::types::real_t learning_rate){
+bool model::add_optimizer_sgd(const panic::types::real_t learning_rate, const panic::types::real_t decay, const panic::types::real_t momentum){
 
-    optimizer_sgd* new_optimizer = new optimizer_sgd(learning_rate);
+    optimizer_sgd* new_optimizer = new optimizer_sgd(learning_rate, decay, momentum);
+
+    if (new_optimizer == 0){
+        return false;
+    }
+
+    delete optimizer_function;
+    optimizer_function = new_optimizer;
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::neural_network::model::add_optimizer_adagrad
+//
+// Description:
+//      Adds Stocastient Gradient Decent as an optimizer to the model.
+//
+// Example:
+//      model.add_optimizer_adagrad(1e-4);
+//--------------------------------------------------------------------------------------------------------------------------
+bool model::add_optimizer_adagrad(const panic::types::real_t learning_rate, const panic::types::real_t decay, const panic::types::real_t epsilon){
+
+    optimizer_adagrad* new_optimizer = new optimizer_adagrad(learning_rate, decay, epsilon);
+
+    if (new_optimizer == 0){
+        return false;
+    }
+
+    delete optimizer_function;
+    optimizer_function = new_optimizer;
+
+    return true;
+}
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::neural_network::model::add_optimizer_rmsprop
+//
+// Description:
+//      Adds Stocastient Gradient Decent as an optimizer to the model.
+//
+// Example:
+//      model.add_optimizer_rmsprop(1e-4);
+//--------------------------------------------------------------------------------------------------------------------------
+bool model::add_optimizer_rmsprop(const panic::types::real_t learning_rate, 
+                                    const panic::types::real_t decay, 
+                                    const panic::types::real_t epsilon,
+                                    const panic::types::real_t rho){
+
+    optimizer_rmsprop* new_optimizer = new optimizer_rmsprop(learning_rate, decay, epsilon, rho);
+
+    if (new_optimizer == 0){
+        return false;
+    }
+
+    delete optimizer_function;
+    optimizer_function = new_optimizer;
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Function Name : panic::neural_network::model::add_optimizer_adam
+//
+// Description:
+//      Adds Stocastient Gradient Decent as an optimizer to the model.
+//
+// Example:
+//      model.add_optimizer_adam(1e-4);
+//--------------------------------------------------------------------------------------------------------------------------
+bool model::add_optimizer_adam(const panic::types::real_t learning_rate, 
+                                    const panic::types::real_t decay, 
+                                    const panic::types::real_t epsilon,
+                                    const panic::types::real_t beta_1,
+                                    const panic::types::real_t beta_2){
+
+    optimizer_adam* new_optimizer = new optimizer_adam(learning_rate, decay, epsilon, beta_1, beta_2);
 
     if (new_optimizer == 0){
         return false;
@@ -618,7 +698,9 @@ bool model::optimize(){
         return false;
     }
 
-
+    if (! optimizer_function->pre_update_params()){
+        return false;
+    }
     
     for(panic::types::uint_t i = 0; i < trainable_layer_count; ++i){
         
@@ -627,10 +709,16 @@ bool model::optimize(){
             return false;
         }
 
+
+
         if (! optimizer_function->update_params(*trainable_layers[i]) ){
             return false;
         }
+    }
 
+
+    if (! optimizer_function->post_update_params()){
+        return false;
     }
 
 
@@ -654,7 +742,7 @@ bool model::train(const panic::tensor::real_matrix& X_train,
     panic::types::real_t accuracy;
     panic::tensor::uint_vector comparisons;
 
-    for (panic::types::uint_t epoch = 0; epoch < epochs; ++epoch){
+    for (panic::types::uint_t epoch = 0; epoch < epochs+1; ++epoch){
 
 
         forward(X_train);
@@ -688,8 +776,9 @@ bool model::train(const panic::tensor::real_matrix& X_train,
         // Its mean is therefore the classification accuracy.
         accuracy = panic::math::mean(comparisons);
 
-        if (epoch % print_every == static_cast<panic::types::uint_t>(0)){
+        if ((epoch) % print_every == static_cast<panic::types::uint_t>(0) ){
             std::cout << "Epoch: " << epoch;
+            std::cout << "      lr: " << optimizer_function->current_learning_rate;
             std::cout << "      data loss: " << loss_function->data_loss;
             std::cout << "      acc: " << accuracy << std::endl;
         }

src/neural_network/optimizers/optimizer_adagrad.cpp

@@ -0,0 +1,217 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       neural_network
+ * File Name:         optimizer_adagrad.cpp
+ * Revision:          0.1.0
+ * Date:              04-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Defines the optimizer_adagrad used in neural network
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//---------------------------------------------------------------------------------------------------------------------------
+#include <neural_network/optimizers/optimizer_adagrad.hpp>
+#include <config/omp.hpp>
+
+#include <math/mul.hpp>
+#include <math/add.hpp>
+#include <math/sub.hpp>
+#include <math/div.hpp>
+#include <math/sqrt.hpp>
+
+
+//---------------------------------------------------------------------------------------------------------------------------
+// PRIVATE CONSTANTS
+//---------------------------------------------------------------------------------------------------------------------------
+/**
+ * @brief Minimum number of element operations before using the OpenMP-enabled loop.
+ *
+ * Small vectors and matrices are kept serial because the overhead of starting
+ * worker threads can be larger than the work itself.
+ */
+static const panic::types::uint_t optimizer_adagrad_omp_min_size = 250;
+
+static const panic::types::real_t real_t_1 = static_cast<panic::types::real_t>(1);
+//---------------------------------------------------------------------------------------------------------------------------
+// INPLEMENTATION
+//---------------------------------------------------------------------------------------------------------------------------
+
+namespace panic{
+	namespace neural_network{
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::optimizer_adagrad
+//
+// Description:
+//		Constructor for optimizer_adagrad.
+//--------------------------------------------------------------------------------------------------------------------------
+optimizer_adagrad::optimizer_adagrad(const panic::types::real_t learning_rate, 
+                                        const panic::types::real_t decay,
+                                        const panic::types::real_t epsilon){
+
+	this->learning_rate = learning_rate;
+    this->decay = decay;
+    this->epsilon = epsilon;
+
+    iterations = static_cast<panic::types::uint_t>(0);
+    current_learning_rate = learning_rate;
+
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::update_params
+//
+// Description:
+//		Updates weights and biases in layer.
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_adagrad::update_params(trainable_layer& layer) {
+
+
+    // Gradients must match their corresponding parameters.
+    if (layer.weights.rows() != layer.dweights.rows() || layer.weights.cols() != layer.dweights.cols() || layer.biases.size() != layer.dbiases.size()){
+        return false;
+    }
+
+    if (layer.weights.rows() == static_cast<panic::types::uint_t>(0) || layer.weights.cols() == static_cast<panic::types::uint_t>(0) || layer.biases.size() == static_cast<panic::types::uint_t>(0) ){
+        return false;
+    }
+
+
+    // If layer doesn't contain cache arrays, create them
+    // and fill them with zeros
+    if(layer.weights.rows() != layer.weight_cache.rows() || 
+        layer.weights.cols() != layer.weight_cache.cols() ||
+        layer.biases.size() != layer.bias_cache.size()){
+        
+        if (! layer.weight_cache.resize(layer.weights.rows(), layer.weights.cols())){
+            return false;
+        }
+        layer.weight_cache.fill(0);
+
+        if (! layer.bias_cache.resize(layer.biases.size())){
+            return false;
+        }
+        layer.bias_cache.fill(0);
+    }
+
+    // Update cache with squared current gradients
+    if(!panic::math::add(
+            layer.bias_cache, 
+            panic::math::mul(layer.dbiases, layer.dbiases), 
+            layer.bias_cache)){
+        return false;
+    }
+
+
+    if(!panic::math::add(
+            layer.weight_cache, 
+            panic::math::mul(layer.dweights, layer.dweights), 
+            layer.weight_cache)){
+        return false;
+    }
+
+
+
+    panic::tensor::real_matrix weight_updates;
+    panic::tensor::real_vector bias_updates;
+
+    // Vanilla SGD parameter update + normalization
+    // with square rooted cache
+
+    // weight_updates = -current_learning_rate * dweights
+    if (!panic::math::mul(layer.dweights, -current_learning_rate, weight_updates)){
+    	return false;
+    }
+    // weight_updates = weight_updates / (sqrt(weight_cache) + epsilon)
+    if (!panic::math::div(weight_updates, panic::math::add(panic::math::sqrt(layer.weight_cache), epsilon), weight_updates)){
+        return false;
+    }
+    // Add weight updates to layer
+    if (! panic::math::add(layer.weights, weight_updates, layer.weights)){
+        return false;
+    }
+
+
+    // bias_updates = -current_learning_rate * dbiases
+    if (!panic::math::mul(layer.dbiases, -current_learning_rate, bias_updates)){
+        return false;
+    }
+    // bias_updates = bias_updates / (sqrt(bias_cache) + epsilon)
+    if (!panic::math::div(bias_updates, panic::math::add(panic::math::sqrt(layer.bias_cache), epsilon), bias_updates)){
+        return false;
+    }
+    // Add bias updates to layer
+    if (! panic::math::add(layer.biases, bias_updates, layer.biases)){
+        return false;
+    }
+
+
+
+    return true;
+	
+}
+
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::pre_update_params
+//
+// Description:
+//      function to update internal parameters before update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_adagrad::pre_update_params(){
+
+    if (decay){
+        current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations)));
+    }
+    
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::post_update_params
+//
+// Description:
+//      function to update internal parameters after update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_adagrad::post_update_params(){
+
+    iterations += static_cast<panic::types::uint_t>(1);
+
+    return true;
+}
+
+
+
+
+
+	} // namespace tensor
+} // namespace panic

src/neural_network/optimizers/optimizer_adam.cpp

@@ -0,0 +1,380 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       neural_network
+ * File Name:         optimizer_adam.cpp
+ * Revision:          0.1.0
+ * Date:              04-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Defines the optimizer_adam used in neural network
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//---------------------------------------------------------------------------------------------------------------------------
+#include <neural_network/optimizers/optimizer_adam.hpp>
+#include <config/omp.hpp>
+
+#include <math/mul.hpp>
+#include <math/add.hpp>
+#include <math/sub.hpp>
+#include <math/div.hpp>
+#include <math/sqrt.hpp>
+
+
+//---------------------------------------------------------------------------------------------------------------------------
+// PRIVATE CONSTANTS
+//---------------------------------------------------------------------------------------------------------------------------
+/**
+ * @brief Minimum number of element operations before using the OpenMP-enabled loop.
+ *
+ * Small vectors and matrices are kept serial because the overhead of starting
+ * worker threads can be larger than the work itself.
+ */
+static const panic::types::uint_t optimizer_adam_omp_min_size = 250;
+
+static const panic::types::real_t real_t_1 = static_cast<panic::types::real_t>(1);
+static const panic::types::real_t zero = static_cast<panic::types::real_t>(0);
+//---------------------------------------------------------------------------------------------------------------------------
+// INPLEMENTATION
+//---------------------------------------------------------------------------------------------------------------------------
+
+namespace panic{
+	namespace neural_network{
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::optimizer_adam
+//
+// Description:
+//		Constructor for optimizer_adam.
+//--------------------------------------------------------------------------------------------------------------------------
+optimizer_adam::optimizer_adam(const panic::types::real_t learning_rate, 
+                                        const panic::types::real_t decay,
+                                        const panic::types::real_t epsilon,
+                                        const panic::types::real_t beta_1,
+                                        const panic::types::real_t beta_2){
+
+	this->learning_rate = learning_rate;
+    this->decay = decay;
+    this->epsilon = epsilon;
+    this->beta_1 = beta_1;
+    this->beta_2 = beta_2;
+
+    this->beta_1_power = beta_1;
+    this->beta_2_power = beta_2;
+
+    iterations = static_cast<panic::types::uint_t>(0);
+    current_learning_rate = learning_rate;
+
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::update_params
+//
+// Description:
+//		Updates weights and biases in layer.
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_adam::update_params(trainable_layer& layer) {
+
+
+    // Gradients must match their corresponding parameters.
+    if (layer.weights.rows() != layer.dweights.rows() || layer.weights.cols() != layer.dweights.cols() || layer.biases.size() != layer.dbiases.size()){
+        return false;
+    }
+
+    if (layer.weights.rows() == static_cast<panic::types::uint_t>(0) || layer.weights.cols() == static_cast<panic::types::uint_t>(0) || layer.biases.size() == static_cast<panic::types::uint_t>(0) ){
+        return false;
+    }
+
+
+    // If layer doesn't contain cache arrays, create them
+    // and fill them with zeros
+    if(layer.weights.rows() != layer.weight_cache.rows() || 
+        layer.weights.cols() != layer.weight_cache.cols() ||
+        layer.biases.size() != layer.bias_cache.size()){
+        
+        if (! layer.weight_cache.resize(layer.weights.rows(), layer.weights.cols())){
+            return false;
+        }
+        layer.weight_cache.fill(0);
+
+        if (! layer.bias_cache.resize(layer.biases.size())){
+            return false;
+        }
+        layer.bias_cache.fill(0);
+    }
+
+    // If layer doesn't contain momentum arrays, create them
+    // and fill them with zeros
+    if(layer.weights.rows() != layer.weight_momentums.rows() || 
+        layer.weights.cols() != layer.weight_momentums.cols() ||
+        layer.biases.size() != layer.bias_momentums.size()){
+        
+        if (! layer.weight_momentums.resize(layer.weights.rows(), layer.weights.cols())){
+            return false;
+        }
+        layer.weight_momentums.fill(0);
+
+        if (! layer.bias_momentums.resize(layer.biases.size())){
+            return false;
+        }
+        layer.bias_momentums.fill(0);
+    }
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Update momentum calculations
+    //--------------------------------------------------------------------------------------------------------------------------
+    const panic::types::real_t one = static_cast<panic::types::real_t>(1);
+    const panic::types::real_t one_minus_beta_1 = one - beta_1;
+
+    panic::tensor::real_matrix retained_weight_momentum;
+    panic::tensor::real_matrix new_weight_gradient;
+
+
+    if (!panic::math::mul(layer.weight_momentums, beta_1, retained_weight_momentum)){
+        return false;
+    }
+
+    if (!panic::math::mul(layer.dweights, one_minus_beta_1, new_weight_gradient)){
+        return false;
+    }
+
+    if (!panic::math::add(retained_weight_momentum, new_weight_gradient, layer.weight_momentums)){
+        return false;
+    }
+
+
+    panic::tensor::real_vector retained_bias_momentum;
+    panic::tensor::real_vector new_bias_gradient;
+
+    if (!panic::math::mul(layer.bias_momentums, beta_1, retained_bias_momentum)){
+        return false;
+    }
+
+    if (!panic::math::mul(layer.dbiases, one_minus_beta_1, new_bias_gradient)){
+        return false;
+    }
+
+    if (!panic::math::add(retained_bias_momentum, new_bias_gradient, layer.bias_momentums)){
+        return false;
+    }
+
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Update momentum corrected
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    const panic::types::real_t momentum_correction = one - beta_1_power;
+
+    if (momentum_correction == zero){
+        return false;
+    }
+
+    panic::tensor::real_matrix corrected_weight_momentums;
+    panic::tensor::real_vector corrected_bias_momentums;
+
+    if (!panic::math::div(layer.weight_momentums, momentum_correction, corrected_weight_momentums)){
+        return false;
+    }
+
+    if (!panic::math::div(layer.bias_momentums, momentum_correction, corrected_bias_momentums)){
+        return false;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Update cache calculations
+    //--------------------------------------------------------------------------------------------------------------------------
+    const panic::types::real_t one_minus_beta_2 = one - beta_2;
+
+    panic::tensor::real_matrix retained_weight_cache;
+    panic::tensor::real_matrix squared_dweights;
+    panic::tensor::real_matrix new_weight_cache_values;
+
+    if (!panic::math::mul(layer.weight_cache, beta_2, retained_weight_cache)){
+        return false;
+    }
+
+    if (!panic::math::mul(layer.dweights, layer.dweights, squared_dweights)){
+        return false;
+    }
+
+    if (!panic::math::mul(squared_dweights, one_minus_beta_2, new_weight_cache_values)){
+        return false;
+    }
+
+    if (!panic::math::add(retained_weight_cache, new_weight_cache_values, layer.weight_cache)){
+        return false;
+    }
+
+
+    panic::tensor::real_vector retained_bias_cache;
+    panic::tensor::real_vector squared_dbiases;
+    panic::tensor::real_vector new_bias_cache_values;
+
+    if (!panic::math::mul(layer.bias_cache, beta_2, retained_bias_cache)){
+        return false;
+    }
+
+    if (!panic::math::mul(layer.dbiases, layer.dbiases, squared_dbiases)){
+        return false;
+    }
+
+    if (!panic::math::mul(squared_dbiases, one_minus_beta_2, new_bias_cache_values)){
+        return false;
+    }
+
+    if (!panic::math::add(retained_bias_cache, new_bias_cache_values, layer.bias_cache)){
+        return false;
+    }
+
+
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Update cache corrected
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    const panic::types::real_t cache_correction = one - beta_2_power;
+
+    if (cache_correction == zero){
+        return false;
+    }
+
+    panic::tensor::real_matrix corrected_weight_cache;
+    panic::tensor::real_vector corrected_bias_cache;
+
+    if (!panic::math::div(layer.weight_cache, cache_correction, corrected_weight_cache)){
+        return false;
+    }
+
+    if (!panic::math::div(layer.bias_cache, cache_correction, corrected_bias_cache)){
+        return false;
+    }
+
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Calculate Adam update
+    //--------------------------------------------------------------------------------------------------------------------------
+
+ 
+
+    panic::tensor::real_matrix sqrt_corrected_weight_cache;
+    panic::tensor::real_matrix weight_denominator;
+    panic::tensor::real_matrix scaled_weight_momentums;
+    panic::tensor::real_matrix weight_updates;
+
+    if (!panic::math::sqrt(corrected_weight_cache, sqrt_corrected_weight_cache)){
+        return false;
+    }
+
+    if (!panic::math::add(sqrt_corrected_weight_cache, epsilon, weight_denominator)){
+        return false;
+    }
+
+    if (!panic::math::mul(corrected_weight_momentums, -current_learning_rate, scaled_weight_momentums)){
+        return false;
+    }
+
+    if (!panic::math::div(scaled_weight_momentums, weight_denominator, weight_updates)){
+        return false;
+    }
+
+    if (!panic::math::add(layer.weights, weight_updates, layer.weights)){
+        return false;
+    }
+
+
+    panic::tensor::real_vector sqrt_corrected_bias_cache;
+    panic::tensor::real_vector bias_denominator;
+    panic::tensor::real_vector scaled_bias_momentums;
+    panic::tensor::real_vector bias_updates;
+
+    if (!panic::math::sqrt(corrected_bias_cache, sqrt_corrected_bias_cache)){
+        return false;
+    }
+
+    if (!panic::math::add(sqrt_corrected_bias_cache, epsilon, bias_denominator)){
+        return false;
+    }
+
+    if (!panic::math::mul(corrected_bias_momentums, -current_learning_rate, scaled_bias_momentums )){
+        return false;
+    }
+
+    if (!panic::math::div(scaled_bias_momentums, bias_denominator, bias_updates)){
+        return false;
+    }
+
+    if (!panic::math::add(layer.biases, bias_updates, layer.biases)){
+        return false;
+    }
+
+
+    return true;
+	
+}
+
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::pre_update_params
+//
+// Description:
+//      function to update internal parameters before update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_adam::pre_update_params(){
+
+    if (decay){
+        current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations)));
+    }
+    
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::post_update_params
+//
+// Description:
+//      function to update internal parameters after update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_adam::post_update_params(){
+
+    iterations += static_cast<panic::types::uint_t>(1);
+
+    beta_1_power *= beta_1;
+    beta_2_power *= beta_2;
+
+    return true;
+}
+
+
+
+
+
+	} // namespace tensor
+} // namespace panic

src/neural_network/optimizers/optimizer_rmsprop.cpp

@@ -0,0 +1,349 @@
+/**++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ *                                              PANIC
+ *                              Portable Algorithms and Numerics In C++
+ *
+ *                       Scientific computing from scratch, with feeling.
+ *
+ * Copyright (c) 2026 Michelle Bausager
+ *
+ * This file is part of PANIC.
+ *
+ * PANIC is free software licensed under the GNU General Public License v3.0 or later.
+ * You may redistribute and/or modify it under the terms of the GPL.
+ *
+ * PANIC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
+ * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the LICENSE file for the full license text.
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ *
+ * Project Name:      PANIC
+ * Module Name:       neural_network
+ * File Name:         optimizer_rmsprop.cpp
+ * Revision:          0.1.0
+ * Date:              04-08-2026
+ * Author:            Michelle Bausager
+ *
+ * Description:
+ * Defines the optimizer_rmsprop used in neural network
+ *
+ *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
+
+//---------------------------------------------------------------------------------------------------------------------------
+// INCLUDE DESCRIPTION
+//---------------------------------------------------------------------------------------------------------------------------
+#include <neural_network/optimizers/optimizer_rmsprop.hpp>
+#include <config/omp.hpp>
+
+#include <math/mul.hpp>
+#include <math/add.hpp>
+#include <math/sub.hpp>
+#include <math/div.hpp>
+#include <math/sqrt.hpp>
+
+
+//---------------------------------------------------------------------------------------------------------------------------
+// PRIVATE CONSTANTS
+//---------------------------------------------------------------------------------------------------------------------------
+/**
+ * @brief Minimum number of element operations before using the OpenMP-enabled loop.
+ *
+ * Small vectors and matrices are kept serial because the overhead of starting
+ * worker threads can be larger than the work itself.
+ */
+static const panic::types::uint_t optimizer_rmsprop_omp_min_size = 250;
+
+static const panic::types::real_t real_t_1 = static_cast<panic::types::real_t>(1);
+//---------------------------------------------------------------------------------------------------------------------------
+// INPLEMENTATION
+//---------------------------------------------------------------------------------------------------------------------------
+
+namespace panic{
+	namespace neural_network{
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::optimizer_rmsprop
+//
+// Description:
+//		Constructor for optimizer_rmsprop.
+//--------------------------------------------------------------------------------------------------------------------------
+optimizer_rmsprop::optimizer_rmsprop(const panic::types::real_t learning_rate, 
+                                        const panic::types::real_t decay,
+                                        const panic::types::real_t epsilon,
+                                        const panic::types::real_t rho){
+
+	this->learning_rate = learning_rate;
+    this->decay = decay;
+    this->epsilon = epsilon;
+    this->rho = rho;
+
+    iterations = static_cast<panic::types::uint_t>(0);
+    current_learning_rate = learning_rate;
+
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::update_params
+//
+// Description:
+//		Updates weights and biases in layer.
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_rmsprop::update_params(trainable_layer& layer) {
+
+
+    // Gradients must match their corresponding parameters.
+    if (layer.weights.rows() != layer.dweights.rows() || layer.weights.cols() != layer.dweights.cols() || layer.biases.size() != layer.dbiases.size()){
+        return false;
+    }
+
+    if (layer.weights.rows() == static_cast<panic::types::uint_t>(0) || layer.weights.cols() == static_cast<panic::types::uint_t>(0) || layer.biases.size() == static_cast<panic::types::uint_t>(0) ){
+        return false;
+    }
+
+
+    // If layer doesn't contain cache arrays, create them
+    // and fill them with zeros
+    if(layer.weights.rows() != layer.weight_cache.rows() || 
+        layer.weights.cols() != layer.weight_cache.cols() ||
+        layer.biases.size() != layer.bias_cache.size()){
+        
+        if (! layer.weight_cache.resize(layer.weights.rows(), layer.weights.cols())){
+            return false;
+        }
+        layer.weight_cache.fill(0);
+
+        if (! layer.bias_cache.resize(layer.biases.size())){
+            return false;
+        }
+        layer.bias_cache.fill(0);
+    }
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Update caches with squared current gradients
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    const panic::types::real_t one_minus_rho =
+        static_cast<panic::types::real_t>(1) - rho;
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Weight cache
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    panic::tensor::real_matrix retained_weight_cache;
+    panic::tensor::real_matrix squared_dweights;
+    panic::tensor::real_matrix weighted_dweights;
+
+    if (!panic::math::mul(
+            layer.weight_cache,
+            rho,
+            retained_weight_cache
+        )){
+        return false;
+    }
+
+    if (!panic::math::mul(
+            layer.dweights,
+            layer.dweights,
+            squared_dweights
+        )){
+        return false;
+    }
+
+    if (!panic::math::mul(
+            squared_dweights,
+            one_minus_rho,
+            weighted_dweights
+        )){
+        return false;
+    }
+
+    if (!panic::math::add(
+            retained_weight_cache,
+            weighted_dweights,
+            layer.weight_cache
+        )){
+        return false;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Bias cache
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    panic::tensor::real_vector retained_bias_cache;
+    panic::tensor::real_vector squared_dbiases;
+    panic::tensor::real_vector weighted_dbiases;
+
+    if (!panic::math::mul(
+            layer.bias_cache,
+            rho,
+            retained_bias_cache
+        )){
+        return false;
+    }
+
+    if (!panic::math::mul(
+            layer.dbiases,
+            layer.dbiases,
+            squared_dbiases
+        )){
+        return false;
+    }
+
+    if (!panic::math::mul(
+            squared_dbiases,
+            one_minus_rho,
+            weighted_dbiases
+        )){
+        return false;
+    }
+
+    if (!panic::math::add(
+            retained_bias_cache,
+            weighted_dbiases,
+            layer.bias_cache
+        )){
+        return false;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Calculate weight updates
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    panic::tensor::real_matrix weight_updates;
+
+    if (!panic::math::mul(
+            layer.dweights,
+            -current_learning_rate,
+            weight_updates
+        )){
+        return false;
+    }
+
+    panic::tensor::real_matrix sqrt_weight_cache;
+    panic::tensor::real_matrix weight_denominator;
+
+    if (!panic::math::sqrt(
+            layer.weight_cache,
+            sqrt_weight_cache
+        )){
+        return false;
+    }
+
+    if (!panic::math::add(
+            sqrt_weight_cache,
+            epsilon,
+            weight_denominator
+        )){
+        return false;
+    }
+
+    if (!panic::math::div(
+            weight_updates,
+            weight_denominator,
+            weight_updates
+        )){
+        return false;
+    }
+
+    if (!panic::math::add(
+            layer.weights,
+            weight_updates,
+            layer.weights
+        )){
+        return false;
+    }
+
+    //--------------------------------------------------------------------------------------------------------------------------
+    // Calculate bias updates
+    //--------------------------------------------------------------------------------------------------------------------------
+
+    panic::tensor::real_vector bias_updates;
+
+    if (!panic::math::mul(
+            layer.dbiases,
+            -current_learning_rate,
+            bias_updates
+        )){
+        return false;
+    }
+
+    panic::tensor::real_vector sqrt_bias_cache;
+    panic::tensor::real_vector bias_denominator;
+
+    if (!panic::math::sqrt(
+            layer.bias_cache,
+            sqrt_bias_cache
+        )){
+        return false;
+    }
+
+    if (!panic::math::add(
+            sqrt_bias_cache,
+            epsilon,
+            bias_denominator
+        )){
+        return false;
+    }
+
+    if (!panic::math::div(
+            bias_updates,
+            bias_denominator,
+            bias_updates
+        )){
+        return false;
+    }
+
+    if (!panic::math::add(
+            layer.biases,
+            bias_updates,
+            layer.biases
+        )){
+        return false;
+    }
+
+    return true;
+	
+}
+
+
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::pre_update_params
+//
+// Description:
+//      function to update internal parameters before update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_rmsprop::pre_update_params(){
+
+    if (decay){
+        current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations)));
+    }
+    
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::post_update_params
+//
+// Description:
+//      function to update internal parameters after update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_rmsprop::post_update_params(){
+
+    iterations += static_cast<panic::types::uint_t>(1);
+
+    return true;
+}
+
+
+
+
+
+	} // namespace tensor
+} // namespace panic

src/neural_network/optimizers/optimizer_sgd.cpp

@@ -40,6 +40,7 @@
 
 #include <math/mul.hpp>
 #include <math/add.hpp>
+#include <math/sub.hpp>
 
 //---------------------------------------------------------------------------------------------------------------------------
 // PRIVATE CONSTANTS
@@ -51,6 +52,8 @@
  * worker threads can be larger than the work itself.
  */
 static const panic::types::uint_t optimizer_sgd_omp_min_size = 250;
+
+static const panic::types::real_t real_t_1 = static_cast<panic::types::real_t>(1);
 //---------------------------------------------------------------------------------------------------------------------------
 // INPLEMENTATION
 //---------------------------------------------------------------------------------------------------------------------------
@@ -65,8 +68,18 @@ namespace panic{
 // Description:
 //		Constructor for optimizer_sgd.
 //--------------------------------------------------------------------------------------------------------------------------
-optimizer_sgd::optimizer_sgd(const panic::types::real_t learning_rate) {
+optimizer_sgd::optimizer_sgd(const panic::types::real_t learning_rate, 
+                                const panic::types::real_t decay,
+                                const panic::types::real_t momentum) {
+
 	this->learning_rate = learning_rate;
+    this->decay = decay;
+    this->momentum = momentum;
+
+    iterations = static_cast<panic::types::uint_t>(0);
+    current_learning_rate = learning_rate;
+
+
 }
 
 //--------------------------------------------------------------------------------------------------------------------------
@@ -87,28 +100,79 @@ bool optimizer_sgd::update_params(trainable_layer& layer) {
         return false;
     }
 
-
     panic::tensor::real_matrix weight_updates;
     panic::tensor::real_vector bias_updates;
 
-    // weight_updates = -learning_rate * dweights
+
-    if (!panic::math::mul(layer.dweights, -learning_rate, weight_updates)){
+    if(momentum){
-    	return false;
+
+        // If layer doesn't contain momentum arrays, create them
+        // and fill them with zeros
+        if(layer.weights.rows() != layer.weight_momentums.rows() || 
+            layer.weights.cols() != layer.weight_momentums.cols() ||
+            layer.biases.size() != layer.bias_momentums.size()){
+            
+            if (! layer.weight_momentums.resize(layer.weights.rows(), layer.weights.cols())){
+                return false;
+            }
+            layer.weight_momentums.fill(0);
+
+            if (! layer.bias_momentums.resize(layer.biases.size())){
+                return false;
+            }
+            layer.bias_momentums.fill(0);
+        }
+
+        // Build weight updates with momentum - take previous
+        // updates multiplied by retain factor and update with
+        // current gradient
+        if (!panic::math::sub(
+            panic::math::mul(layer.weight_momentums, momentum), 
+            panic::math::mul(layer.dweights, current_learning_rate), 
+            weight_updates)){
+
+            return false;
+        }
+
+        // Build bias updates
+        if (!panic::math::sub(
+            panic::math::mul(layer.bias_momentums, momentum), 
+            panic::math::mul(layer.dbiases, current_learning_rate), 
+            bias_updates)){
+
+            return false;
+        }
+
+        layer.weight_momentums = weight_updates;
+        layer.bias_momentums = bias_updates;
+
+
+    }
+    else{
+
+        // weight_updates = -current_learning_rate * dweights
+        if (!panic::math::mul(layer.dweights, -current_learning_rate, weight_updates)){
+        	return false;
+        }
+
+
+
+        // bias_updates = -current_learning_rate * dbiases
+        if (!panic::math::mul(layer.dbiases, -current_learning_rate, bias_updates)){
+        	return false;
+        }
+
+
     }
 
     // weights += weight_updates
     if (!panic::math::add(layer.weights, weight_updates, layer.weights)){
-    	return false;
+        return false;
-    }
-
-    // bias_updates = -learning_rate * dbiases
-    if (!panic::math::mul(layer.dbiases, -learning_rate, bias_updates)){
-    	return false;
     }
 
     // biases += bias_updates
     if (!panic::math::add(layer.biases, bias_updates, layer.biases)){
-    	return false;
+        return false;
     }
 
     return true;
@@ -116,5 +180,39 @@ bool optimizer_sgd::update_params(trainable_layer& layer) {
 }
 
 
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::pre_update_params
+//
+// Description:
+//      function to update internal parameters before update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_sgd::pre_update_params(){
+
+    if (decay){
+        current_learning_rate = learning_rate * (real_t_1 / (real_t_1 + (decay * iterations)));
+    }
+    
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------
+// Constructor Name : panic::neural_network::post_update_params
+//
+// Description:
+//      function to update internal parameters after update_params()
+//--------------------------------------------------------------------------------------------------------------------------
+bool optimizer_sgd::post_update_params(){
+
+    iterations += static_cast<panic::types::uint_t>(1);
+
+    return true;
+}
+
+
+
+
+
 	} // namespace tensor
 } // namespace panic