/**
 * @file
 * @brief [Gram Schmidt Orthogonalisation
 * Process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process)
 *
 * @details
 * Takes the input of Linearly Independent Vectors,
 * returns vectors orthogonal to each other.
 *
 * ### Algorithm
 * Take the first vector of given LI vectors as first vector of Orthogonal
 * vectors. Take projection of second input vector on the first vector of
 * Orthogonal vector and subtract it from the 2nd LI vector. Take projection of
 * third vector on the second vector of Othogonal vectors and subtract it from
 * the 3rd LI vector. Keep repeating the above process until all the vectors in
 * the given input array are exhausted.
 *
 * For Example:
 * In R2,
 * Input LI Vectors={(3,1),(2,2)}
 * then Orthogonal Vectors= {(3, 1),(-0.4, 1.2)}
 *
 *  Have defined maximum dimension of vectors to be 10 and number of vectors
 *  taken is 20.
 *  Please do not give linearly dependent vectors
 *
 *
 * @author [Akanksha Gupta](https://github.com/Akanksha-Gupta920)
 */

#include <array>     /// for std::array
#include <cassert>   /// for assert
#include <cmath>     /// for fabs
#include <iostream>  /// for io operations

#include "math.h"

/**
 * @namespace numerical_methods
 * @brief Numerical Methods algorithms
 */
namespace numerical_methods {
/**
 * @namespace gram_schmidt
 * @brief Functions for [Gram Schmidt Orthogonalisation
 * Process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process)
 */
namespace gram_schmidt {
/**
 * Dot product function.
 * Takes 2 vectors along with their dimension as input and returns the dot
 * product.
 * @param x vector 1
 * @param y vector 2
 * @param c dimension of the vectors
 *
 * @returns sum
 */
double dot_product(const std::array<double, 10>& x,
                   const std::array<double, 10>& y, const int& c) {
    double sum = 0;
    for (int i = 0; i < c; ++i) {
        sum += x[i] * y[i];
    }
    return sum;
}

/**
 * Projection Function
 * Takes input of 2 vectors along with their dimension and evaluates their
 * projection in temp
 *
 * @param x Vector 1
 * @param y Vector 2
 * @param c dimension of each vector
 *
 * @returns factor
 */
double projection(const std::array<double, 10>& x,
                  const std::array<double, 10>& y, const int& c) {
    double dot =
        dot_product(x, y, c);  /// The dot product of two vectors is taken
    double anorm =
        dot_product(y, y, c);  /// The norm of the second vector is taken.
    double factor =
        dot /
        anorm;  /// multiply that factor with every element in a 3rd vector,
                /// whose initial values are same as the 2nd vector.
    return factor;
}

/**
 * Function to print the orthogonalised vector
 *
 * @param r number of vectors
 * @param c dimenaion of vectors
 * @param B stores orthogonalised vectors
 *
 * @returns void
 */
void display(const int& r, const int& c,
             const std::array<std::array<double, 10>, 20>& B) {
    for (int i = 0; i < r; ++i) {
        std::cout << "Vector " << i + 1 << ": ";
        for (int j = 0; j < c; ++j) {
            std::cout << B[i][j] << " ";
        }
        std::cout << '\n';
    }
}

/**
 * Function for the process of Gram Schimdt Process
 * @param r number of vectors
 * @param c dimension of vectors
 * @param A stores input of given LI vectors
 * @param B stores orthogonalised vectors
 *
 * @returns void
 */
void gram_schmidt(int r, const int& c,
                  const std::array<std::array<double, 10>, 20>& A,
                  std::array<std::array<double, 10>, 20> B) {
    if (c < r) {  /// we check whether appropriate dimensions are given or not.
        std::cout << "Dimension of vector is less than number of vector, hence "
                     "\n first "
                  << c << " vectors are orthogonalised\n";
        r = c;
    }

    int k = 1;

    while (k <= r) {
        if (k == 1) {
            for (int j = 0; j < c; j++)
                B[0][j] = A[0][j];  /// First vector is copied as it is.
        }

        else {
            std::array<double, 10>
                all_projection{};  /// array to store projections
            for (int i = 0; i < c; ++i) {
                all_projection[i] = 0;  /// First initialised to zero
            }

            int l = 1;
            while (l < k) {
                std::array<double, 10>
                    temp{};           /// to store previous projected array
                double factor = NAN;  /// to store the factor by which the
                                      /// previous array will change
                factor = projection(A[k - 1], B[l - 1], c);
                for (int i = 0; i < c; ++i) {
                    temp[i] = B[l - 1][i] * factor;  /// projected array created
                }
                for (int j = 0; j < c; ++j) {
                    all_projection[j] =
                        all_projection[j] +
                        temp[j];  /// we take the projection with all the
                                  /// previous vector and add them.
                }
                l++;
            }
            for (int i = 0; i < c; ++i) {
                B[k - 1][i] =
                    A[k - 1][i] -
                    all_projection[i];  /// subtract total projection vector
                                        /// from the input vector
            }
        }
        k++;
    }
    display(r, c, B);  // for displaying orthogoanlised vectors
}
}  // namespace gram_schmidt
}  // namespace numerical_methods
/**
 * Test Function. Process has been tested for 3 Sample Inputs
 * @returns void
 */
static void test() {
    std::array<std::array<double, 10>, 20> a1 = {
        {{1, 0, 1, 0}, {1, 1, 1, 1}, {0, 1, 2, 1}}};
    std::array<std::array<double, 10>, 20> b1 = {{0}};
    double dot1 = 0;
    numerical_methods::gram_schmidt::gram_schmidt(3, 4, a1, b1);
    int flag = 1;
    for (int i = 0; i < 2; ++i) {
        for (int j = i + 1; j < 3; ++j) {
            dot1 = fabs(
                numerical_methods::gram_schmidt::dot_product(b1[i], b1[j], 4));
            if (dot1 > 0.1) {
                flag = 0;
                break;
            }
        }
    }
    if (flag == 0)
        std::cout << "Vectors are linearly dependent\n";
    assert(flag == 1);
    std::cout << "Passed Test Case 1\n ";

    std::array<std::array<double, 10>, 20> a2 = {{{3, 1}, {2, 2}}};
    std::array<std::array<double, 10>, 20> b2 = {{0}};
    double dot2 = 0;
    numerical_methods::gram_schmidt::gram_schmidt(2, 2, a2, b2);
    flag = 1;
    for (int i = 0; i < 1; ++i) {
        for (int j = i + 1; j < 2; ++j) {
            dot2 = fabs(
                numerical_methods::gram_schmidt::dot_product(b2[i], b2[j], 2));
            if (dot2 > 0.1) {
                flag = 0;
                break;
            }
        }
    }
    if (flag == 0)
        std::cout << "Vectors are linearly dependent\n";
    assert(flag == 1);
    std::cout << "Passed Test Case 2\n";

    std::array<std::array<double, 10>, 20> a3 = {{{1, 2, 2}, {-4, 3, 2}}};
    std::array<std::array<double, 10>, 20> b3 = {{0}};
    double dot3 = 0;
    numerical_methods::gram_schmidt::gram_schmidt(2, 3, a3, b3);
    flag = 1;
    for (int i = 0; i < 1; ++i) {
        for (int j = i + 1; j < 2; ++j) {
            dot3 = fabs(
                numerical_methods::gram_schmidt::dot_product(b3[i], b3[j], 3));
            if (dot3 > 0.1) {
                flag = 0;
                break;
            }
        }
    }
    if (flag == 0)
        std::cout << "Vectors are linearly dependent\n";
    assert(flag == 1);
    std::cout << "Passed Test Case 3\n";
}

/**
 * @brief Main Function
 * @return 0 on exit
 */
int main() {
    int r = 0, c = 0;
    test();  // perform self tests
    std::cout << "Enter the dimension of your vectors\n";
    std::cin >> c;
    std::cout << "Enter the number of vectors you will enter\n";
    std::cin >> r;

    std::array<std::array<double, 10>, 20>
        A{};  /// a 2-D array for storing all vectors
    std::array<std::array<double, 10>, 20> B = {
        {0}};  /// a 2-D array for storing orthogonalised vectors
    /// storing vectors in array A
    for (int i = 0; i < r; ++i) {
        std::cout << "Enter vector " << i + 1
                  << '\n';  /// Input of vectors is taken
        for (int j = 0; j < c; ++j) {
            std::cout << "Value " << j + 1 << "th of vector: ";
            std::cin >> A[i][j];
        }
        std::cout << '\n';
    }

    numerical_methods::gram_schmidt::gram_schmidt(r, c, A, B);

    double dot = 0;
    int flag = 1;  /// To check whether vectors are orthogonal or  not
    for (int i = 0; i < r - 1; ++i) {
        for (int j = i + 1; j < r; ++j) {
            dot = fabs(
                numerical_methods::gram_schmidt::dot_product(B[i], B[j], c));
            if (dot > 0.1)  /// take make the process numerically stable, upper
                            /// bound for the dot product take 0.1
            {
                flag = 0;
                break;
            }
        }
    }
    if (flag == 0)
        std::cout << "Vectors are linearly dependent\n";
    return 0;
}