Matrix.hpp

/******************************************************************************
* Copyright (c) ATTX LLC 2024. All Rights Reserved.
*
* This software and associated documentation (the "Software") are the 
* proprietary and confidential information of ATTX, LLC. The Software is 
* furnished under a license agreement between ATTX and the user organization 
* and may be used or copied only in accordance with the terms of the agreement.
* Refer to 'license/attx_license.adoc' for standard license terms.
*
* EXPORT CONTROL NOTICE: THIS SOFTWARE MAY INCLUDE CONTENT CONTROLLED UNDER THE
* INTERNATIONAL TRAFFIC IN ARMS REGULATIONS (ITAR) OR THE EXPORT ADMINISTRATION 
* REGULATIONS (EAR99). No part of the Software may be used, reproduced, or 
* transmitted in any form or by any means, for any purpose, without the express 
* written permission of ATTX, LLC.
******************************************************************************/
/*
Matrix math header file
 
Author: Alex Reynolds
*/
#ifndef MATRIX_HPP
#define MATRIX_HPP
 
#include <array>
#include <iostream>
 
#include "clockwerkerrors.h"
#include "safemath.hpp"
 
namespace clockwerk {
 
    /// @brief   Matrix math implementation
    /// 
    /// This file defines a base class for computational matrix math.
    /// It is designed for use in real-time, embedded mode or with
    /// simulation. As such, it contains certain limitations on dynamic
    /// allocation to make it real-time compliant. 
    /// 
    /// Note: The naming convention observed for matrix math is as follows:
    /// - Capital letters represent matrices
    /// - Lowercase letters represent scalars
    /// - All vectors are 1-dimensional matrices for the purpose of math
    /// - Pass by reference return values are all named "result"
    /// - The letter used in math indicates the order of operations - A/a
    ///   comes first, then B/b. A/a or B/b always represent the "other"
    ///   value in an operation
    /// 
    /// Important: The matrix class is templated such that it can be used
    /// with a generic type T, but should only be used on native types,
    /// as other implementations could result in undefined behavior.
    /// 
    /// TODO: Flattening and flattened constructor
    template <typename T, unsigned int R, unsigned int C>
    class Matrix {
    public:
        /// Default constructor to initialize a matrix to all zeroes for ease of use
        Matrix<T, R, C>();
 
        /// Constructor to initialize a matrix with all the same element
        Matrix<T, R, C>(T elements);
 
        /// Constructor for Matrix class with initialization. 
        /// Initializes matrix to values passed in via array
        Matrix<T, R, C>(const T(&initial)[R][C]);
 
        /// Copy constructor for Matrix class. Copies data from matrix
        /// object to the current instance
        Matrix<T, R, C>(const Matrix<T, R, C> &initial);
 
        /// Array constructor for the matrix class. Initializes from std::array
        Matrix<T, R, C>(const std::array<std::array<T, C>, R> &initial);
 
        /// Destructor -- doesn't do anything because we don't dynamically
        /// allocate
        ~Matrix<T, R, C>() {}
 
        /// @brief   Function to dump information on matrix
        void dump() const;
 
        /// @brief   Function to dump information on matrix to string
        std::string str() const;
 
        /// ----------------------------------------------------------
        /// Getters and setters
        /// ----------------------------------------------------------
        /// @brief   Function to set a single value in the matrix
        /// @param   row     The row index
        /// @param   col     The column index
        /// @param   value   The value to set the matrix to
        /// @return  Integer value corresponding to error codes in clockwerkerrrors.h
        int set(const unsigned int &row, const unsigned int &col, const T &value);
 
        /// @brief   Function to get a single value in the matrix
        /// @param   row     The row index
        /// @param   col     The column index
        /// @param   value   PBR return of the value in the matrix
        /// @return  Integer value corresponding to error codes in clockwerkerrrors.h
        int get(const unsigned int &row, const unsigned int &col, T &result) const;
 
        /// @brief Function to set the values of the matrix row-wise
        /// @param start_ptr The data address at which read should begin. 
        /// @note start_ptr must represent an array or similar type of at least R*C values
        /// @note Not protected. May segfault if not used properly
        /// @note The array [1 2 3 4] would read into a matrix as [1 2; 3 4]
        void setFromArray(const T* start_ptr);
 
        /// @brief Function to get the values of the matrix row-wise
        /// @param start_ptr The data address at which write should begin. 
        /// @note start_ptr must represent an array or similar type of at least R*C values
        /// @note Not protected. May segfault if not used properly
        /// @note The matrix [1 2; 3 4] would output as an array [1 2 3 4]
        void getAsArray(T* start_ptr) const;
 
        /// @brief   Function to get a single value in the matrix
        /// @param   row     The row index
        /// @param   col     The column index
        /// @return  Matrix value at index
        /// @note Does not bounds check and is unsafe for embedded operations
        T get(const unsigned int &row, const unsigned int &col) const;
 
        /// @brief   Function to get a copy of the matrix
        /// @param   result  PBR return of a copy of this matrix
        void getCopy(Matrix<T, R, C> &result) const;
 
    /// @brief    Equals operator overload for matrix 
        Matrix<T, R, C>& operator=(const Matrix<T, R, C>& other);
 
        /// @brief Function to return a matrix row or vector value
        /// @param idx The row index to return
        /// @return A pointer to the row element of the matrix. Nullptr if out of range
        T* operator [](unsigned int idx);
 
        /// ----------------------------------------------------------
        /// Operations to get matrix info
        /// ----------------------------------------------------------
 
        /// Function to get the size of the matrix
        std::pair<unsigned int, unsigned int> size() {return std::pair<unsigned int, unsigned int>(R, C);}
 
        /// @brief   Function to return the maximum value in the matrix
        /// @param   result  PBR return of maximum value
        void max(T &result, std::pair<unsigned int, unsigned int> &index) const;
 
        /// @brief   Function to return the minimum value in the matrix
        /// @param   result  PBR return of minimum value
        void min(T &result, std::pair<unsigned int, unsigned int> &index) const;
 
        /// ----------------------------------------------------------
        /// Operations to return important matrix information
        /// ----------------------------------------------------------
 
        /// @brief   Function to return the determinant of the matrix
        /// @param   result  PBR return of determinant
        /// @return  Integer value corresponding to error codes in clockwerkerrrors.h
        int det(T &result) const;
 
        /// @brief   Function to return the inverse of the matrix
        /// @param   result  PBR return of matrix inverse
        /// @return  Integer value corresponding to error codes in clockwerkerrrors.h
        int inverse(Matrix<T, R, C> &result) const;
        Matrix<T, R, C> inverse() const {Matrix<T, R, C> tmp; inverse(tmp); return tmp;}
 
        /// @brief   Function to return the transpose of the matrix
        /// @param   result  PBR return of matrix transpose
        /// @note    Overloaded for return as well
        void transpose(Matrix<T, C, R> &result) const;
        Matrix<T, C, R> transpose() const;
 
        /// @brief   Function to return the trace of the matrix
        /// @param   result  PBR return of matrix trace
        /// @return  Integer value corresponding to error codes in clockwerkerrrors.h
        int trace(T &result) const;
 
        /// @brief   Function to set all elements of the matrix to zero
        void setToZeros();
 
        /// @brief Function to set matrix to identity, if it is a square matrix
        /// @return Error code corresponding to success/failure
        int identity();
        int eye() {return identity();}
 
        /// The actual values held by the matrix -- a two dimensional 
        /// array of values indexed as (row, column). 
        /// NOTE: Public for ease of access and speed (no need to use setter/getter),
        /// and functions in Matrix and Safemath libraries are tested safe. Behavior
        /// using matrices outside of clockwerk libraries should use setter/getter, rather
        /// than direct access, for safety.
        std::array<std::array<T, C>, R> values;
 
    protected:
        /// Function to check, given a set of submatrix boundaries, that those boundaries are
        /// valid for the current matrix
        int _checkLookupBoundaries(const unsigned int &start_r, const unsigned int &end_r,
                                const unsigned int &start_c, const unsigned int &end_c) const;
 
        /// @brief   Function to take a 2-d matrix represented by A and decompose it into LU form
        /// @param   A   A 2-d matrix represented as a double pointer
        ///              After ops A contains a copy of both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
        /// @param   P   The permutation matrix to hold info on matrix changes
        /// @return  Error code corresponding to codes in clockwerkerrors.h
        /// @note    Code borrowed/updated from sample C code on wikipedia here: 
        ///          https://en.wikipedia.org/wiki/LU_decomposition
        int _LUPDecompose(T *A[R], unsigned int P[R+1]) const;
 
    };
 
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C>::Matrix() {
        unsigned int i, j;
        /// Explicitly set initialized values to element
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[i][j] = T();
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C>::Matrix(T elements) {
        unsigned int i, j;
        /// Explicitly set initialized values to element
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[i][j] = elements;
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C>::Matrix(const T(&initial)[R][C]) {
        /// Note -- one absolutely essential check that's occurring here is that the size
        /// of the initializer list matches the size of the matrix. This is occurring implicitly
        /// at compile time via templating and is the reason why an array is used rather than
        /// a std::initializer_list
 
        /// Initialize to values from initializer list
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[i][j] = initial[i][j];
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C>::Matrix(const Matrix<T, R, C> &initial) {
        unsigned int i, j;
        /// Explicitly set initialized values to zero
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[i][j] = initial.values[i][j];
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C>::Matrix(const std::array<std::array<T, C>, R> &initial) {
        unsigned int i, j;
        /// Loop through matrix and copy our values
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[i][j] = initial[i][j];
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    T* Matrix<T, R, C>::operator [](unsigned int idx) {
        // Return pointer to our array element
        if(idx < R) {
            return &values[idx][0];
        } else {
            return nullptr;
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::dump() const {
        unsigned int i, j;
        /// Loop through matrix and write out values
        for(i = 0; i < R; i++) {
            std::cout<<"[";
            for(j = 0; j < C; j++) {
                if(j > 0) {
                    std::cout<<", ";
                }
                std::cout<<values[i][j];
            }
            std::cout<<"]"<<std::endl;
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    std::string Matrix<T, R, C>::str() const {
        unsigned int i, j;
        std::string out = "[";
        /// Loop through matrix and write out values
        for(i = 0; i < R; i++) {
            out += "[";
            for(j = 0; j < C; j++) {
                if(j > 0) {
                    out += ",";
                }
                out += std::to_string(values[i][j]);
            }
            out += "]";
        }
 
        out += "]";
        return out;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::set(const unsigned int &row, const unsigned int &col, const T &value){
        /// Check boundaries on lookup values to ensure validity 
        if(_checkLookupBoundaries(row, row, col, col)) {
            /// If we returned true, we had a dimension mismatch. Return error
            return ERROR_DIMENSIONS;
        }
 
        /// Set value
        values[row][col] = value;
        return NO_ERROR;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::get(const unsigned int &row, const unsigned int &col, T &result) const{
        /// Check boundaries on lookup values to ensure validity 
        if(_checkLookupBoundaries(row, row, col, col)) {
            /// If we returned true, we had a dimension mismatch. Return error
            return ERROR_DIMENSIONS;
        }
 
        /// Set value
        result = values[row][col];
        return NO_ERROR;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    T Matrix<T, R, C>::get(const unsigned int &row, const unsigned int &col) const{
        return values[row][col];
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::setFromArray(const T* start_ptr) {
        const T* ptr = start_ptr;
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[i][j] = *ptr;
                ptr++;
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::getAsArray(T* start_ptr) const {
        T* ptr = start_ptr;
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                *ptr = values[i][j];
                ptr++;
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::getCopy(Matrix<T, R, C> &result) const {
        /// NOTE: This function does not need to check lookup boundaries because
        /// templating will enforce. Code will not compile if dimensions are incorrect
 
        unsigned int i, j;
        /// Loop through matrix and search for maximum value
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[i][j] = values[i][j];
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C>& Matrix<T, R, C>::operator=(const Matrix<T, R, C>& other)
    {
        /// Guard self assignment
        if (this == &other)
            return *this;
 
        /// Copy values from current matrix to return
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                this->values[i][j] = other.values[i][j];
            }
        }
        return *this;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::max(T &result, std::pair<unsigned int, unsigned int> &index) const {
        unsigned int i, j;
        index.first = 0;
        index.second = 0;
        result = values[0][0];      /// Set our initial index to first value
        /// Loop through matrix and search for maximum value
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                if(values[i][j] > result) {
                    result = values[i][j];
                    index.first = i;
                    index.second = j;
                }
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::min(T &result, std::pair<unsigned int, unsigned int> &index) const {
        unsigned int i, j;
        index.first = 0;
        index.second = 0;
        result = values[0][0];      /// Set our initial index to first value
        /// Loop through matrix and search for minimum value
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                if(values[i][j] < result) {
                    result = values[i][j];
                    index.first = i;
                    index.second = j;
                }
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::det(T &result) const {
        /// Note: This function needs to check that we have a square matrix to calculate
        /// determinant. That check is perfomed in the function _LUPDecompose, and therefore
        /// not performed here.
 
        /// Initialize some tracking variables for this function
        unsigned int i, j, err;
 
        /// Now that we know we have a square matrix, initialize a copy of our matrix as a
        /// 2-d array. We need to do this because we're going to decompose the matrix into
        /// upper triangular form in order to take its determinant, but we don't want to 
        /// actually modify the base matrix. Thus, the copy.
        T A[R][C];
        T *decomposed[R];
        for(i = 0; i < R; i++) {
            decomposed[i] = A[i];
            for(j = 0; j < C; j++) {
                A[i][j] = values[i][j];
            }
        }
 
        /// Define our permutation matrix, which exists purely to track changes to the base
        /// matrix. This is tracked as a single array.
        /// Note: from this point forward using R to indicate size. R and C are equivalent
        /// or there will be an error thrown, so which used doesn't matter.
        unsigned int P[R+1];
 
        /// Now decompose our matrix into its L*U form for easy determinant calculation.
        /// If our function returns that there was a poorly conditioned matrix (rank deficient),
        /// we know that the determinant is zero and didn't actually have an error. Otherwise,
        /// if our function returns that our matrix is not square, return error
        err = _LUPDecompose(decomposed, P);
        if(err == ERROR_POORLY_CONDITIONED) {
            result = 0;
            return NO_ERROR;
        } else if(err) {
            return err;
        }
 
        /// Now, calculate our determinant simply using known relationship for upper triangular matrices
        result = decomposed[0][0];
        for(i = 1; i < R; i++) {
            result *= decomposed[i][i];
        }
        if((P[R] - R)%2) {
            result = -result;
        }
 
        return NO_ERROR;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::inverse(Matrix<T, R, C> &result) const {
        /// Note: This function needs to check that we have a square matrix to calculate
        /// determinant. That check is perfomed in the function _LUPDecompose, and therefore
        /// not performed here.
 
        /// Note: Need to check that we have matching sizes between matrix and output -- this
        /// is handled automatically by compiler and doesn't need to be checked here
 
        /// Initialize some tracking variables for this function
        unsigned int i, j, k, idx, err;
 
        /// Now that we know we have a square matrix, initialize a copy of our matrix as a
        /// 2-d array. We need to do this because we're going to decompose the matrix into
        /// upper triangular form in order to take its determinant, but we don't want to 
        /// actually modify the base matrix. Thus, the copy.
        T A[R][C];
        T *decomposed[R];
        for(i = 0; i < R; i++) {
            decomposed[i] = A[i];
            for(j = 0; j < C; j++) {
                A[i][j] = values[i][j];
            }
        }
 
        /// Define our permutation matrix, which exists purely to track changes to the base
        /// matrix. This is tracked as a single array.
        /// Note: from this point forward using R to indicate size. R and C are equivalent
        /// or there will be an error thrown, so which used doesn't matter.
        unsigned int P[R+1];
 
        /// Now decompose our matrix into its L*U form for easy determinant calculation. Check
        /// that there was not an error returned due to poorly conditioned matrix
        err = _LUPDecompose(decomposed, P);
        if(err) {
            return err;
        }
 
        /// Now, calculate our inverse and assign to our return value
        for (j = 0; j < R; j++) {
            for (i = 0; i < R; i++) {
                result.values[i][j] = P[i] == j ? 1.0 : 0.0;
                for (k = 0; k < i; k++) {
                    result.values[i][j] -= decomposed[i][k]*result.values[k][j];
                }
            }
 
 
            for(idx = 1; idx < R+1; idx++) {
                i = R - idx;
                for (k = i+1; k < R; k++) {
                    result.values[i][j] -= decomposed[i][k]*result.values[k][j];
                }
                err = safeDivide(result.values[i][j], decomposed[i][i], result.values[i][j]);
                if(err) {;
                    return err;
                }
            }
        }
 
        return NO_ERROR;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::transpose(Matrix<T, C, R> &result) const {
        /// NOTE: No need to enforce size constraints on matrix here because
        /// they will be enforced by the compiler. Code will not compile
        /// if matrix of the wrong size is fed to this function
 
        unsigned int i, j;
        /// Loop through matrix and write out values in transposed form
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[j][i] = values[i][j];
            }
        }
    }
    template <typename T, unsigned int R, unsigned int C>
    Matrix<T, C, R> Matrix<T, R, C>::transpose() const {
        /// NOTE: No need to enforce size constraints on matrix here because
        /// they will be enforced by the compiler. Code will not compile
        /// if matrix of the wrong size is fed to this function
        Matrix<T, C, R> tmp;
        transpose(tmp);
        return tmp;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::trace(T &result) const {
        /// Ensure we have a square matrix -- error if not
        /// TODO: See if this should be checked at compile time for templating speed
        if(R != C) {
            return ERROR_DIMENSIONS;
        }
 
        result = T();
        unsigned int i;
        /// Now calculate and return trace
        for(i = 0; i < R; i++) {
            result += values[i][i];
        }
        return NO_ERROR;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    void Matrix<T, R, C>::setToZeros() {
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                values[j][i] = 0;
            }
        }
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::identity() {
        if(C != R) {return ERROR_DIMENSIONS;}
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                if(i == j) {
                    values[i][i] = 1.0;
                } else {
                    values[i][j] = 0.0;
                }
            }
        }
 
        return NO_ERROR;
    }
 
    template <typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::_checkLookupBoundaries(const unsigned int &start_r, const unsigned int &end_r,
                                                const unsigned int &start_c, const unsigned int &end_c) const{
        /// Check that start is not greater than end
        if(start_r > end_r || start_c > end_c) {
            return ERROR_DIMENSIONS;
        }
        /// Check that end does not exceed matrix size. Note we do not need to check start
        /// because it is an unsigned int where we're checking start !> end
        if(end_r >= R || end_c >= C) {
            return ERROR_DIMENSIONS;
        }
        /// We passed our checks - return OK
        return NO_ERROR;
    }
 
    template<typename T, unsigned int R, unsigned int C>
    int Matrix<T, R, C>::_LUPDecompose(T *A[R], unsigned int P[R+1]) const {
        /// Check to confirm that our matrix is square. If it is not, return with error.
        /// Note: afterwards we know R and C are equivalent and interchangeable. We'll use
        /// just R moving forward
        /// TODO: Figure out if this should be done via macro at compile time due to templating
        if(R != C) {
            return ERROR_DIMENSIONS;
        }
 
        /// Initialize our iterators, etc. for the function
        unsigned int i, j, k, imax;
        T max_A, *ptr, abs_A;
 
        /// Initialize our unit permutation matrix. P[N] is initialized with N
        for (k = 0; k <= R; k++) {
            P[k] = k;
        }
 
        /// Loop through matrix columnwise 
        for (i = 0; i < R; i++) {
            /// Reset max_A and imax to default to start loop. Explicitly invoke
            /// default constructor on max_A for consistency
            max_A = T();
            imax = i;
 
            /// Loop through elements of column N to identify row with maximum value. This
            /// will be identified by imax and used in a later pivot operation. The pivot is
            /// not necessary to take the determinant, but provides better numerical stability
            for (k = i; k < R; k++) {
                abs_A = std::abs(A[k][i]);      /// Get abs value of individual element
                /// If absolute value is greater than previous max, set max and index
                if (abs_A > max_A) {
                    max_A = abs_A;
                    imax = k;
                }
            }
 
            /// If it was not possible to find a value greater than zero (within tolerance)
            /// in the current column, we know our matrix is rank deficient and should return
            /// as such. This ensures a zero determinant as well.
            if (max_A < TOLERANCE_CONDITIONING)
                return ERROR_POORLY_CONDITIONED; //failure, matrix is degenerate
 
            /// If the first value of the current row is not the highest value, we will pivot
            /// for numerical stability in the operation. Swap pointers to current and pivot
            /// row and record change in Permutation matrix
            if (imax != i) {
                /// Count pivots that occurred. This tracking is for determinant purposes --
                /// Each pivot flips the sign of the determinant and must be tracked
                P[R]++;
 
                /// Track pivot in P "matrix"
                j = P[i];
                P[i] = P[imax];
                P[imax] = j;
 
                /// Pivot the rows of A by rearranging pointers for efficiency
                ptr = A[i];
                A[i] = A[imax];
                A[imax] = ptr;
            }
 
            /// Now loop through our matrix and perform the actual gaussian elimination
            for (j = i+1; j < R; j++) {
                safeDivide(A[j][i], A[i][i], A[j][i]);
                for (k = i+1; k < R; k++) {
                    A[j][k] -= A[j][i] * A[i][k];
                }
            }
        }
 
        /// We reached the end and did not hit an error -- return good job
        return NO_ERROR;
    }
 
}
 
#endif