Merge branch 'master' into is_graph_bipartite

CMakeLists.txt

@@ -40,7 +40,8 @@ add_subdirectory(backtracking)
 add_subdirectory(data_structures)
 add_subdirectory(machine_learning)
 add_subdirectory(numerical_methods)
-    
+add_subdirectory(graph)
+
 cmake_policy(SET CMP0054 NEW)
 cmake_policy(SET CMP0057 NEW)
 find_package(Doxygen OPTIONAL_COMPONENTS dot dia)

CONTRIBUTING.md

@@ -23,7 +23,8 @@ We are very happy that you consider implementing algorithms and data structures
 - Please use the directory structure of the repository.
 - File extension for code should be *.h *.cpp.
 - Don't use **bits/stdc++.h** because this is quite Linux specific and slows down the compilation process.
-- Avoid using **struct** and instead use the **class** keyword.
+- Organize your code using **`struct`**, **`class`** and/or **`namespace`** keywords
+- If an implementation of the algorithm already exists, please refer to the [file-name section below](#new-file-name-guidelines).
 - You can suggest reasonable changes to existing algorithms.
 - Strictly use snake_case (underscore_separated) in filenames.
 - If you have added or modified code, please make sure the code compiles before submitting.
@@ -116,6 +117,7 @@ my_new_cpp_class.cpp    is correct format
 ```
 - It will be used to dynamically create a directory of files and implementation.
 - File name validation will run on docker to ensure the validity.
+- If an implementation of the algorithm already exists and your version is different from that implemented, please use incremental numeric digit as a suffix. For example, if `median_search.cpp` already exists in the `search` folder and you are contributing a new implementation, the filename should be `median_search2.cpp` and for a third implementation, `median_search3.cpp`.
 
 #### New Directory guidelines 
 - We recommend adding files to existing directories as much as possible.

DIRECTORY.md

@@ -67,9 +67,11 @@
   * [Longest Increasing Subsequence (Nlogn)](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/dynamic_programming/longest_increasing_subsequence_(nlogn).cpp)
   * [Matrix Chain Multiplication](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/dynamic_programming/matrix_chain_multiplication.cpp)
   * [Searching Of Element In Dynamic Array](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/dynamic_programming/searching_of_element_in_dynamic_array.cpp)
+  * [Shortest Common Supersequence](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/dynamic_programming/shortest_common_supersequence.cpp)
   * [Tree Height](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/dynamic_programming/tree_height.cpp)
 
 ## Geometry
+  * [Jarvis Algorithm](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/geometry/jarvis_algorithm.cpp)
   * [Line Segment Intersection](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/geometry/line_segment_intersection.cpp)
 
 ## Graph
@@ -78,8 +80,8 @@
   * [Connected Components](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/connected_components.cpp)
   * [Connected Components With Dsu](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/connected_components_with_dsu.cpp)
   * [Cycle Check Directed Graph](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/cycle_check_directed_graph.cpp)
-  * [Dfs](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/dfs.cpp)
-  * [Dfs With Stack](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/dfs_with_stack.cpp)
+  * [Depth First Search](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/depth_first_search.cpp)
+  * [Depth First Search With Stack](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/depth_first_search_with_stack.cpp)
   * [Dijkstra](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/dijkstra.cpp)
   * [Hamiltons Cycle](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/hamiltons_cycle.cpp)
   * [Is Graph Bipartite](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/graph/is_graph_bipartite.cpp)
@@ -111,7 +113,9 @@
   * [Adaline Learning](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/adaline_learning.cpp)
   * [Kohonen Som Topology](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/kohonen_som_topology.cpp)
   * [Kohonen Som Trace](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/kohonen_som_trace.cpp)
+  * [Neural Network](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/neural_network.cpp)
   * [Ordinary Least Squares Regressor](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/ordinary_least_squares_regressor.cpp)
+  * [Vector Ops](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/vector_ops.hpp)
 
 ## Math
   * [Armstrong Number](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/math/armstrong_number.cpp)
@@ -180,7 +184,7 @@
   * [Decimal To Binary](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/decimal_to_binary.cpp)
   * [Decimal To Hexadecimal](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/decimal_to_hexadecimal.cpp)
   * [Decimal To Roman Numeral](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/decimal_to_roman_numeral.cpp)
-  * [Fast Interger Input](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/fast_interger_input.cpp)
+  * [Fast Integer Input](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/fast_integer_input.cpp)
   * [Happy Number](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/happy_number.cpp)
   * [Matrix Exponentiation](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/matrix_exponentiation.cpp)
   * [Palindrome Of Number](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/others/palindrome_of_number.cpp)
@@ -202,6 +206,7 @@
 
 ## Range Queries
   * [Fenwick Tree](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/range_queries/fenwick_tree.cpp)
+  * [Heavy Light Decomposition](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/range_queries/heavy_light_decomposition.cpp)
   * [Mo](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/range_queries/mo.cpp)
   * [Segtree](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/range_queries/segtree.cpp)
 
@@ -236,6 +241,7 @@
   * [Non Recursive Merge Sort](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/non_recursive_merge_sort.cpp)
   * [Numeric String Sort](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/numeric_string_sort.cpp)
   * [Odd Even Sort](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/odd_even_sort.cpp)
+  * [Pigeonhole Sort](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/pigeonhole_sort.cpp)
   * [Quick Sort](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/quick_sort.cpp)
   * [Quick Sort 3](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/quick_sort_3.cpp)
   * [Radix Sort](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/sorting/radix_sort.cpp)

dynamic_programming/shortest_common_supersequence.cpp

@@ -0,0 +1,179 @@
+/**
+ * @file
+ * @brief SCS is a string Z which is the shortest supersequence of strings X and Y (may not be continuous in Z, but order is maintained).
+ *
+ * @details
+ * The idea is to use lookup table method as used in LCS.
+ * For example: example 1:-
+ * X: 'ABCXYZ', Y: 'ABZ' then Z will be 'ABCXYZ' (y is not continuous but in order)
+ * 
+ * For example: example 2:-
+ * X: 'AGGTAB', Y: 'GXTXAYB' then Z will be 'AGGXTXAYB'
+ * @author [Ridhish Jain](https://github.com/ridhishjain)
+ * @see more on [SCS](https://en.wikipedia.org/wiki/Shortest_common_supersequence_problem)
+ * @see related problem [Leetcode](https://leetcode.com/problems/shortest-common-supersequence/)
+*/
+
+// header files
+#include <iostream>
+#include <string>
+#include <vector>
+#include <algorithm>
+#include <cassert>
+
+/**
+ * @namespace dynamic_programming
+ * @brief Dynamic Programming algorithms
+*/
+namespace dynamic_programming {
+
+    /**
+    * @namespace shortest_common_supersequence
+    * @brief Shortest Common Super Sequence algorithm
+    */
+    namespace shortest_common_supersequence {
+        
+        /**
+         * Function implementing Shortest Common Super-Sequence algorithm using look-up table method.
+         * @param str1 first string 'X'
+         * @param str2 second string 'Y'
+         * @returns string 'Z', superSequence of X and Y 
+        */
+        std::string scs(const std::string &str1, const std::string &str2) {
+
+            // Edge cases
+            // If either str1 or str2 or both are empty
+            if(str1.empty() && str2.empty()) {
+                return "";
+            }
+            else if(str1.empty()) {
+                return str2;
+            }
+            else if(str2.empty()) {
+                return str1;
+            }
+
+            // creating lookup table
+            std::vector <std::vector <int>> lookup(str1.length() + 1, std::vector <int> (str2.length() + 1, 0));
+      
+            for(int i=1; i <= str1.length(); i++) {
+                for(int j=1; j <= str2.length(); j++) {
+                    if(str1[i-1] == str2[j-1]) {
+                        lookup[i][j] = lookup[i-1][j-1] + 1;
+                    }
+                    else {
+                        lookup[i][j] = std::max(lookup[i-1][j], lookup[i][j-1]);
+                    }
+                }
+            }
+
+            // making supersequence
+            // i and j are initially pointed towards end of strings
+            // Super-sequence will be constructed backwards
+            int i=str1.length();
+            int j=str2.length();
+            std::string s;
+      
+            while(i>0 && j>0) {
+
+                // If the characters at i and j of both strings are same
+                // We only need to add them once in s
+                if(str1[i-1] == str2[j-1]) {
+                    s.push_back(str1[i-1]);
+                    i--;
+                    j--;
+                }
+                // otherwise we check lookup table for recurrences of characters
+                else {
+                    if(lookup[i-1][j] > lookup[i][j-1]) {
+                        s.push_back(str1[i-1]);
+                        i--;
+                    }
+                    else {
+                        s.push_back(str2[j-1]);
+                        j--;
+                    }
+                }
+            }
+
+            // copying remaining elements
+            // if j becomes 0 before i
+            while(i > 0) {
+                s.push_back(str1[i-1]);
+                i--;
+            }
+
+            // if i becomes 0 before j
+            while(j > 0) {
+                s.push_back(str2[j-1]);
+                j--;
+            }
+
+            // As the super sequence is constructd backwards
+            // reversing the string before returning gives us the correct output  
+            reverse(s.begin(), s.end());
+            return s;
+        }
+    } // namespace shortest_common_supersequence
+} // namespace dynamic_programming
+
+/** 
+ * Test Function
+ * @return void 
+*/
+static void test() {
+    // custom input vector
+    std::vector <std::vector <std::string>> scsStrings {
+        {"ABCXYZ", "ABZ"},
+        {"ABZ", "ABCXYZ"},
+        {"AGGTAB", "GXTXAYB"},
+        {"X", "Y"},
+    };
+
+    // calculated output vector by scs function
+    std::vector <std::string> calculatedOutput(4, "");
+    int i=0;
+    for(auto & scsString : scsStrings) {
+        
+        calculatedOutput[i] = dynamic_programming::shortest_common_supersequence::scs(
+            scsString[0], scsString[1]
+        );
+        i++;
+    }
+
+    // expected output vector acc to problem statement
+    std::vector <std::string> expectedOutput {
+        "ABCXYZ",
+        "ABCXYZ",
+        "AGGXTXAYB",
+        "XY"
+    };
+
+    // Testing implementation via assert function
+    // It will throw error if any of the expected test fails
+    // Else it will give nothing
+    for(int i=0; i < scsStrings.size(); i++) {
+        assert(expectedOutput[i] == calculatedOutput[i]);
+    }
+
+    std::cout << "All tests passed successfully!\n";
+    return;
+}
+
+/** Main function (driver code)*/
+int main() {
+    // test for implementation
+    test();
+
+    // user input
+    std::string s1, s2;
+    std::cin >> s1;
+    std::cin >> s2;
+
+    std::string ans;
+
+    // user output
+    ans = dynamic_programming::shortest_common_supersequence::scs(s1, s2);
+    std::cout << ans;
+    return 0;
+}

geometry/jarvis_algorithm.cpp

@@ -0,0 +1,179 @@
+/**
+ * @file
+ * @brief Implementation of [Jarvis’s](https://en.wikipedia.org/wiki/Gift_wrapping_algorithm) algorithm.
+ *
+ * @details
+ * Given a set of points in the plane. the convex hull of the set
+ * is the smallest convex polygon that contains all the points of it.
+ *
+ * ### Algorithm
+ * The idea of Jarvis’s Algorithm is simple, we start from the leftmost point
+ * (or point with minimum x coordinate value) and we
+ * keep wrapping points in counterclockwise direction.
+ *
+ * The idea is to use orientation() here. Next point is selected as the
+ * point that beats all other points at counterclockwise orientation, i.e.,
+ * next point is q if for any other point r,
+ * we have “orientation(p, q, r) = counterclockwise”.
+ *
+ * For Example,
+ * If points = {{0, 3}, {2, 2}, {1, 1}, {2, 1},
+                      {3, 0}, {0, 0}, {3, 3}};
+ *
+ * then the convex hull is
+ * (0, 3), (0, 0), (3, 0), (3, 3)
+ *
+ * @author [Rishabh Agarwal](https://github.com/rishabh-997)
+ */
+
+#include <vector>
+#include <cassert>
+#include <iostream>
+
+/**
+ *  @namespace geometry
+ *  @brief Geometry algorithms
+ */
+namespace geometry {
+    /**
+     * @namespace jarvis
+     * @brief Functions for [Jarvis’s](https://en.wikipedia.org/wiki/Gift_wrapping_algorithm) algorithm
+     */
+    namespace jarvis {
+        /**
+         * Structure defining the x and y co-ordinates of the given
+         * point in space
+         */
+        struct Point {
+            int x, y;
+        };
+
+        /**
+         * Class which can be called from main and is globally available
+         * throughout the code
+         */
+        class Convexhull {
+            std::vector<Point> points;
+            int size;
+
+        public:
+            /**
+             * Constructor of given class
+             *
+             * @param pointList list of all points in the space
+             * @param n number of points in space
+             */
+            explicit Convexhull(const std::vector<Point> &pointList) {
+                points = pointList;
+                size = points.size();
+            }
+
+            /**
+             * Creates convex hull of a set of n points.
+             * There must be 3 points at least for the convex hull to exist
+             *
+             * @returns an vector array containing points in space
+             * which enclose all given points thus forming a hull
+             */
+            std::vector<Point> getConvexHull() const {
+                // Initialize Result
+                std::vector<Point> hull;
+
+                // Find the leftmost point
+                int leftmost_point = 0;
+                for (int i = 1; i < size; i++) {
+                    if (points[i].x < points[leftmost_point].x) {
+                        leftmost_point = i;
+                    }
+                }
+                // Start from leftmost point, keep moving counterclockwise
+                // until reach the start point again.  This loop runs O(h)
+                // times where h is number of points in result or output.
+                int p = leftmost_point, q = 0;
+                do {
+                    // Add current point to result
+                    hull.push_back(points[p]);
+
+                    // Search for a point 'q' such that orientation(p, x, q)
+                    // is counterclockwise for all points 'x'. The idea
+                    // is to keep track of last visited most counter clock-
+                    // wise point in q. If any point 'i' is more counter clock-
+                    // wise than q, then update q.
+                    q = (p + 1) % size;
+                    for (int i = 0; i < size; i++) {
+                        // If i is more counterclockwise than current q, then
+                        // update q
+                        if (orientation(points[p], points[i], points[q]) == 2) {
+                            q = i;
+                        }
+                    }
+
+                    // Now q is the most counterclockwise with respect to p
+                    // Set p as q for next iteration, so that q is added to
+                    // result 'hull'
+                    p = q;
+
+                } while (p != leftmost_point);        // While we don't come to first point
+
+                return hull;
+            }
+
+            /**
+             * This function returns the geometric orientation for the three points
+             * in a space, ie, whether they are linear ir clockwise or
+             * anti-clockwise
+             * @param p first point selected
+             * @param q adjacent point for q
+             * @param r adjacent point for q
+             *
+             * @returns 0 -> Linear
+             * @returns 1 -> Clock Wise
+             * @returns 2 -> Anti Clock Wise
+             */
+            static int orientation(const Point &p, const Point &q, const Point &r) {
+                int val = (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y);
+
+                if (val == 0) {
+                    return 0;
+                }
+                return (val > 0) ? 1 : 2;
+            }
+
+        };
+
+    } // namespace jarvis
+} // namespace geometry
+
+/**
+ * Test function
+ * @returns void
+ */
+static void test() {
+    std::vector<geometry::jarvis::Point> points = {{0, 3},
+                                                   {2, 2},
+                                                   {1, 1},
+                                                   {2, 1},
+                                                   {3, 0},
+                                                   {0, 0},
+                                                   {3, 3}
+    };
+    geometry::jarvis::Convexhull hull(points);
+    std::vector<geometry::jarvis::Point> actualPoint;
+    actualPoint = hull.getConvexHull();
+
+    std::vector<geometry::jarvis::Point> expectedPoint = {{0, 3},
+                                                          {0, 0},
+                                                          {3, 0},
+                                                          {3, 3}};
+    for (int i = 0; i < expectedPoint.size(); i++) {
+        assert(actualPoint[i].x == expectedPoint[i].x);
+        assert(actualPoint[i].y == expectedPoint[i].y);
+    }
+    std::cout << "Test implementations passed!\n";
+}
+
+/** Driver Code */
+int main() {
+    test();
+    return 0;
+}

graph/CMakeLists.txt

@@ -0,0 +1,20 @@
+# If necessary, use the RELATIVE flag, otherwise each source file may be listed
+# with full pathname. RELATIVE may makes it easier to extract an executable name
+# automatically.
+file( GLOB APP_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.cpp )
+# file( GLOB APP_SOURCES ${CMAKE_SOURCE_DIR}/*.c )
+# AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR} APP_SOURCES)
+foreach( testsourcefile ${APP_SOURCES} )
+    # I used a simple string replace, to cut off .cpp.
+    string( REPLACE ".cpp" "" testname ${testsourcefile} )
+    add_executable( ${testname} ${testsourcefile} )
+
+    set_target_properties(${testname} PROPERTIES
+        LINKER_LANGUAGE CXX
+    )
+    if(OpenMP_CXX_FOUND)
+        target_link_libraries(${testname} OpenMP::OpenMP_CXX)
+    endif()
+    install(TARGETS ${testname} DESTINATION "bin/graph")
+
+endforeach( testsourcefile ${APP_SOURCES} )

graph/cycle_check_directed_graph.cpp

@@ -53,7 +53,7 @@ using AdjList = std::map<unsigned int, std::vector<unsigned int>>;
  */
 class Graph {
  public:
-    Graph() : m_vertices(0), m_adjList({}) {}
+    Graph() : m_adjList({}) {}
     ~Graph() = default;
     Graph(Graph&&) = default;
     Graph& operator=(Graph&&) = default;
@@ -65,8 +65,8 @@ class Graph {
      * @param vertices specify the number of vertices the graph would contain.
      * @param adjList is the adjacency list representation of graph.
      */
-    Graph(unsigned int vertices, AdjList const& adjList)
-        : m_vertices(vertices), m_adjList(adjList) {}
+    Graph(unsigned int vertices, AdjList adjList)
+        : m_vertices(vertices), m_adjList(std::move(adjList)) {}
 
     /** Create a graph from vertices and adjacency list.
      *
@@ -142,7 +142,7 @@ class Graph {
     }
 
  private:
-    unsigned int m_vertices;
+    unsigned int m_vertices = 0;
     AdjList m_adjList;
 };
 

graph/depth_first_search.cpp

@@ -0,0 +1,133 @@
+/**
+ *
+ * \file
+ * \brief [Depth First Search Algorithm
+ * (Depth First Search)](https://en.wikipedia.org/wiki/Depth-first_search)
+ *
+ * \author [Ayaan Khan](http://github.com/ayaankhan98)
+ *
+ * \details
+ * Depth First Search also quoted as DFS is a Graph Traversal Algorithm.
+ * Time Complexity O(|V| + |E|) where V is number of vertices and E
+ * is number of edges in graph.
+ *
+ * Application of Depth First Search are
+ *
+ * 1. Finding connected components
+ * 2. Finding 2-(edge or vertex)-connected components.
+ * 3. Finding 3-(edge or vertex)-connected components.
+ * 4. Finding the bridges of a graph.
+ * 5. Generating words in order to plot the limit set of a group.
+ * 6. Finding strongly connected components.
+ *
+ * And there are many more...
+ *
+ * <h4>Working</h4>
+ * 1. Mark all vertices as unvisited first
+ * 2. start exploring from some starting vertex.
+ *
+ *      While exploring vertex we mark the vertex as visited
+ *      and start exploring the vertices connected to this
+ *      vertex in recursive way.
+ *
+ */
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+
+/**
+ *
+ * \namespace graph
+ * \brief Graph Algorithms
+ *
+ */
+namespace graph {
+/**
+ * \brief
+ * Adds and edge between two vertices of graph say u and v in this
+ * case.
+ *
+ * @param adj Adjacency list representation of graph
+ * @param u first vertex
+ * @param v second vertex
+ *
+ */
+void addEdge(std::vector<std::vector<size_t>> *adj, size_t u, size_t v) {
+    /*
+     *
+     * Here we are considering undirected graph that's the
+     * reason we are adding v to the adjacency list representation of u
+     * and also adding u to the adjacency list representation of v
+     *
+     */
+    (*adj)[u - 1].push_back(v - 1);
+    (*adj)[v - 1].push_back(u - 1);
+}
+
+/**
+ *
+ * \brief
+ * Explores the given vertex, exploring a vertex means traversing
+ * over all the vertices which are connected to the vertex that is
+ * currently being explored.
+ *
+ * @param adj garph
+ * @param v vertex to be explored
+ * @param visited already visited vertices
+ *
+ */
+void explore(const std::vector<std::vector<size_t>> &adj, size_t v,
+             std::vector<bool> *visited) {
+    std::cout << v + 1 << " ";
+    (*visited)[v] = true;
+    for (auto x : adj[v]) {
+        if (!(*visited)[x]) {
+            explore(adj, x, visited);
+        }
+    }
+}
+
+/**
+ * \brief
+ * initiates depth first search algorithm.
+ *
+ * @param adj adjacency list of graph
+ * @param start vertex from where DFS starts traversing.
+ *
+ */
+void depth_first_search(const std::vector<std::vector<size_t>> &adj,
+                        size_t start) {
+    size_t vertices = adj.size();
+
+    std::vector<bool> visited(vertices, false);
+    explore(adj, start, &visited);
+}
+}  // namespace graph
+
+/** Main function */
+int main() {
+    size_t vertices = 0, edges = 0;
+    std::cout << "Enter the Vertices : ";
+    std::cin >> vertices;
+    std::cout << "Enter the Edges : ";
+    std::cin >> edges;
+
+    /// creating graph
+    std::vector<std::vector<size_t>> adj(vertices, std::vector<size_t>());
+
+    /// taking input for edges
+    std::cout << "Enter the vertices which have edges between them : "
+              << std::endl;
+    while (edges--) {
+        size_t u = 0, v = 0;
+        std::cin >> u >> v;
+        graph::addEdge(&adj, u, v);
+    }
+
+    /// running depth first search over graph
+    graph::depth_first_search(adj, 2);
+
+    std::cout << std::endl;
+    return 0;
+}

graph/depth_first_search_with_stack.cpp

@@ -0,0 +1,46 @@
+#include <iostream>
+#include <list>
+#include <vector>
+#include <stack>
+
+constexpr int WHITE = 0;
+constexpr int GREY = 1;
+constexpr int BLACK = 2;
+constexpr int INF = 99999;
+
+void dfs(const std::vector< std::list<int> > &graph, int start) {
+    std::vector<int> checked(graph.size(), WHITE);
+    checked[start] = GREY;
+    std::stack<int> stack;
+    stack.push(start);
+    while (!stack.empty()) {
+        int act = stack.top();
+        stack.pop();
+
+        if (checked[act] == GREY) {
+            std::cout << act << ' ';
+            for (auto it : graph[act]) {
+                stack.push(it);
+                if (checked[it] != BLACK) {
+                    checked[it] = GREY;
+                }
+            }
+            checked[act] = BLACK;  // nodo controllato
+        }
+    }
+}
+
+int main() {
+    int n = 0;
+    std::cin >> n;
+    std::vector< std::list<int> > graph(INF);
+    for (int i = 0; i < n; ++i) {
+        int u = 0, w = 0;
+        std::cin >> u >> w;
+        graph[u].push_back(w);
+    }
+
+    dfs(graph, 0);
+
+    return 0;
+}

graph/dfs.cpp

@@ -1,26 +0,0 @@
-#include <iostream>
-using namespace std;
-int v = 4;
-void DFSUtil_(int graph[4][4], bool visited[], int s) {
-    visited[s] = true;
-    cout << s << " ";
-    for (int i = 0; i < v; i++) {
-        if (graph[s][i] == 1 && visited[i] == false) {
-            DFSUtil_(graph, visited, i);
-        }
-    }
-}
-
-void DFS_(int graph[4][4], int s) {
-    bool visited[v];
-    memset(visited, 0, sizeof(visited));
-    DFSUtil_(graph, visited, s);
-}
-
-int main() {
-    int graph[4][4] = {{0, 1, 1, 0}, {0, 0, 1, 0}, {1, 0, 0, 1}, {0, 0, 0, 1}};
-    cout << "DFS: ";
-    DFS_(graph, 2);
-    cout << endl;
-    return 0;
-}

graph/dfs_with_stack.cpp

@@ -1,51 +0,0 @@
-#include <iostream>
-#include <list>
-#include <stack>
-
-#define WHITE 0
-#define GREY 1
-#define BLACK 2
-#define INF 99999
-
-using namespace std;
-
-int checked[999] = {WHITE};
-
-void dfs(const list<int> lista[], int start) {
-    stack<int> stack;
-
-    int checked[999] = {WHITE};
-
-    stack.push(start);
-
-    checked[start] = GREY;
-    while (!stack.empty()) {
-        int act = stack.top();
-        stack.pop();
-
-        if (checked[act] == GREY) {
-            cout << act << ' ';
-            for (auto it = lista[act].begin(); it != lista[act].end(); ++it) {
-                stack.push(*it);
-                if (checked[*it] != BLACK)
-                    checked[*it] = GREY;
-            }
-            checked[act] = BLACK;  // nodo controllato
-        }
-    }
-}
-
-int main() {
-    int u, w;
-    int n;
-    cin >> n;
-    list<int> lista[INF];
-    for (int i = 0; i < n; ++i) {
-        cin >> u >> w;
-        lista[u].push_back(w);
-    }
-
-    dfs(lista, 0);
-
-    return 0;
-}

graph/kosaraju.cpp

@@ -4,77 +4,84 @@
 
 #include <iostream>
 #include <vector>
+#include <stack>
 
-using namespace std;
 
 /**
  * Iterative function/method to print graph:
- * @param a[] : array of vectors (2D)
- * @param V : vertices
+ * @param a adjacency list representation of the graph
+ * @param V number of vertices
  * @return void
  **/
-void print(vector<int> a[], int V) {
+void print(const std::vector< std::vector<int> > &a, int V) {
     for (int i = 0; i < V; i++) {
-        if (!a[i].empty())
-            cout << "i=" << i << "-->";
-        for (int j = 0; j < a[i].size(); j++) cout << a[i][j] << " ";
-        if (!a[i].empty())
-            cout << endl;
+        if (!a[i].empty()) {
+            std::cout << "i=" << i << "-->";
+        }
+        for (int j : a[i]) {
+            std::cout << j << " ";
+        }
+        if (!a[i].empty()) {
+            std::cout << std::endl;
+        }
     }
 }
 
 /**
  * //Recursive function/method to push vertices into stack passed as parameter:
- * @param v : vertices
- * @param &st : stack passed by reference
- * @param vis[] : array to keep track of visited nodes (boolean type)
- * @param adj[] : array of vectors to represent graph
+ * @param v vertices
+ * @param st stack passed by reference
+ * @param vis array to keep track of visited nodes (boolean type)
+ * @param adj adjacency list representation of the graph
  * @return void
  **/
-void push_vertex(int v, stack<int> &st, bool vis[], vector<int> adj[]) {
-    vis[v] = true;
+void push_vertex(int v, std::stack<int> *st, std::vector<bool> *vis, const std::vector< std::vector<int> > &adj) {
+    (*vis)[v] = true;
     for (auto i = adj[v].begin(); i != adj[v].end(); i++) {
-        if (vis[*i] == false)
+        if ((*vis)[*i] == false) {
             push_vertex(*i, st, vis, adj);
+        }
     }
-    st.push(v);
+    st->push(v);
 }
 
 /**
  * //Recursive function/method to implement depth first traversal(dfs):
- * @param v : vertices
- * @param vis[] : array to keep track of visited nodes (boolean type)
- * @param grev[] : graph with reversed edges
+ * @param v vertices
+ * @param vis array to keep track of visited nodes (boolean type)
+ * @param grev graph with reversed edges
  * @return void
  **/
-void dfs(int v, bool vis[], vector<int> grev[]) {
-    vis[v] = true;
+void dfs(int v, std::vector<bool> *vis, const std::vector< std::vector<int> > &grev) {
+    (*vis)[v] = true;
     // cout<<v<<" ";
     for (auto i = grev[v].begin(); i != grev[v].end(); i++) {
-        if (vis[*i] == false)
+        if ((*vis)[*i] == false) {
             dfs(*i, vis, grev);
+        }
     }
 }
 
 // function/method to implement Kosaraju's Algorithm:
 /**
 * Info about the method
-* @param V : vertices in graph
-* @param adj[] : array of vectors that represent a graph (adjacency list/array)
+* @param V vertices in graph
+* @param adj array of vectors that represent a graph (adjacency list/array)
 * @return int ( 0, 1, 2..and so on, only unsigned values as either there can be
 no SCCs i.e. none(0) or there will be x no. of SCCs (x>0)) i.e. it returns the
 count of (number of) strongly connected components (SCCs) in the graph.
 (variable 'count_scc' within function)
 **/
-int kosaraju(int V, vector<int> adj[]) {
-    bool vis[V] = {};
-    stack<int> st;
+int kosaraju(int V, const std::vector< std::vector<int> > &adj) {
+    std::vector<bool> vis(V, false);
+    std::stack<int> st;
     for (int v = 0; v < V; v++) {
-        if (vis[v] == false)
-            push_vertex(v, st, vis, adj);
+        if (vis[v] == false) {
+            push_vertex(v, &st, &vis, adj);
+        }
     }
     // making new graph (grev) with reverse edges as in adj[]:
-    vector<int> grev[V];
+    std::vector< std::vector<int> > grev(V);
     for (int i = 0; i < V + 1; i++) {
         for (auto j = adj[i].begin(); j != adj[i].end(); j++) {
             grev[*j].push_back(i);
@@ -89,7 +96,7 @@ int kosaraju(int V, vector<int> adj[]) {
         int t = st.top();
         st.pop();
         if (vis[t] == false) {
-            dfs(t, vis, grev);
+            dfs(t, &vis, grev);
             count_scc++;
         }
     }
@@ -101,21 +108,21 @@ int kosaraju(int V, vector<int> adj[]) {
 // All critical/corner cases have been taken care of.
 // Input your required values: (not hardcoded)
 int main() {
-    int t;
-    cin >> t;
+    int t = 0;
+    std::cin >> t;
     while (t--) {
-        int a, b;  // a->number of nodes, b->directed edges.
-        cin >> a >> b;
-        int m, n;
-        vector<int> adj[a + 1];
+        int a = 0, b = 0;  // a->number of nodes, b->directed edges.
+        std::cin >> a >> b;
+        int m = 0, n = 0;
+        std::vector< std::vector<int> > adj(a + 1);
         for (int i = 0; i < b; i++)  // take total b inputs of 2 vertices each
                                      // required to form an edge.
         {
-            cin >> m >> n;  // take input m,n denoting edge from m->n.
+            std::cin >> m >> n;  // take input m,n denoting edge from m->n.
             adj[m].push_back(n);
         }
         // pass number of nodes and adjacency array as parameters to function:
-        cout << kosaraju(a, adj) << endl;
+        std::cout << kosaraju(a, adj) << std::endl;
     }
     return 0;
 }

graph/kruskal.cpp

@@ -1,73 +1,21 @@
 #include <iostream>
+#include <vector>
+#include <algorithm>
+#include <array>
 //#include <boost/multiprecision/cpp_int.hpp>
 // using namespace boost::multiprecision;
 const int mx = 1e6 + 5;
-const long int inf = 2e9;
-typedef long long ll;
-#define rep(i, n) for (i = 0; i < n; i++)
-#define repp(i, a, b) for (i = a; i <= b; i++)
-#define pii pair<int, int>
-#define vpii vector<pii>
-#define vi vector<int>
-#define vll vector<ll>
-#define r(x) scanf("%d", &x)
-#define rs(s) scanf("%s", s)
-#define gc getchar_unlocked
-#define pc putchar_unlocked
-#define mp make_pair
-#define pb push_back
-#define lb lower_bound
-#define ub upper_bound
-#define endl "\n"
-#define fast                          \
-    ios_base::sync_with_stdio(false); \
-    cin.tie(NULL);                    \
-    cout.tie(NULL);
-using namespace std;
-void in(int &x) {
-    register int c = gc();
-    x = 0;
-    int neg = 0;
-    for (; ((c < 48 || c > 57) && c != '-'); c = gc())
-        ;
-    if (c == '-') {
-        neg = 1;
-        c = gc();
-    }
-    for (; c > 47 && c < 58; c = gc()) {
-        x = (x << 1) + (x << 3) + c - 48;
-    }
-    if (neg)
-        x = -x;
-}
-void out(int n) {
-    int N = n, rev, count = 0;
-    rev = N;
-    if (N == 0) {
-        pc('0');
-        return;
-    }
-    while ((rev % 10) == 0) {
-        count++;
-        rev /= 10;
-    }
-    rev = 0;
-    while (N != 0) {
-        rev = (rev << 3) + (rev << 1) + N % 10;
-        N /= 10;
-    }
-    while (rev != 0) {
-        pc(rev % 10 + '0');
-        rev /= 10;
-    }
-    while (count--) pc('0');
-}
-ll parent[mx], arr[mx], node, edge;
-vector<pair<ll, pair<ll, ll>>> v;
+using ll = int64_t;
+
+std::array<ll, mx> parent;
+ll node, edge;
+std::vector<std::pair<ll, std::pair<ll, ll>>> edges;
 void initial() {
-    int i;
-    rep(i, node + edge) parent[i] = i;
+    for (int i = 0; i < node + edge; ++i) {
+      parent[i] = i;
+    }
 }
+
 int root(int i) {
     while (parent[i] != i) {
         parent[i] = parent[parent[i]];
@@ -75,41 +23,42 @@ int root(int i) {
     }
     return i;
 }
+
 void join(int x, int y) {
     int root_x = root(x);  // Disjoint set union by rank
     int root_y = root(y);
     parent[root_x] = root_y;
 }
+
 ll kruskal() {
-    ll mincost = 0, i, x, y;
-    rep(i, edge) {
-        x = v[i].second.first;
-        y = v[i].second.second;
+    ll mincost = 0;
+    for (int i = 0; i < edge; ++i) {
+        ll x = edges[i].second.first;
+        ll y = edges[i].second.second;
         if (root(x) != root(y)) {
-            mincost += v[i].first;
+            mincost += edges[i].first;
             join(x, y);
         }
     }
     return mincost;
 }
+
 int main() {
-    fast;
-    while (1) {
-        int i, j, from, to, cost, totalcost = 0;
-        cin >> node >> edge;  // Enter the nodes and edges
-        if (node == 0 && edge == 0)
+    while (true) {
+        int from = 0, to = 0, cost = 0, totalcost = 0;
+        std::cin >> node >> edge;  // Enter the nodes and edges
+        if (node == 0 && edge == 0) {
             break;  // Enter 0 0 to break out
+        }
         initial();  // Initialise the parent array
-        rep(i, edge) {
-            cin >> from >> to >> cost;
-            v.pb(mp(cost, mp(from, to)));
+        for (int i = 0; i < edge; ++i) {
+            std::cin >> from >> to >> cost;
+            edges.emplace_back(make_pair(cost, std::make_pair(from, to)));
             totalcost += cost;
         }
-        sort(v.begin(), v.end());
-        // rep(i,v.size())
-        // 	cout<<v[i].first<<"  ";
-        cout << kruskal() << endl;
-        v.clear();
+        sort(edges.begin(), edges.end());
+        std::cout << kruskal() << std::endl;
+        edges.clear();
     }
     return 0;
 }

machine_learning/iris.csv

@@ -0,0 +1,152 @@
+https://archive.ics.uci.edu/ml/datasets/iris
+sepal length in cm,sepal width in cm,petal length in cm,petal width in cm
+5.1,3.5,1.4,.2,0
+4.9,3,1.4,.2,0
+4.7,3.2,1.3,.2,0
+4.6,3.1,1.5,.2,0
+5,3.6,1.4,.2,0
+5.4,3.9,1.7,.4,0
+4.6,3.4,1.4,.3,0
+5,3.4,1.5,.2,0
+4.4,2.9,1.4,.2,0
+4.9,3.1,1.5,.1,0
+5.4,3.7,1.5,.2,0
+4.8,3.4,1.6,.2,0
+4.8,3,1.4,.1,0
+4.3,3,1.1,.1,0
+5.8,4,1.2,.2,0
+5.7,4.4,1.5,.4,0
+5.4,3.9,1.3,.4,0
+5.1,3.5,1.4,.3,0
+5.7,3.8,1.7,.3,0
+5.1,3.8,1.5,.3,0
+5.4,3.4,1.7,.2,0
+5.1,3.7,1.5,.4,0
+4.6,3.6,1,.2,0
+5.1,3.3,1.7,.5,0
+4.8,3.4,1.9,.2,0
+5,3,1.6,.2,0
+5,3.4,1.6,.4,0
+5.2,3.5,1.5,.2,0
+5.2,3.4,1.4,.2,0
+4.7,3.2,1.6,.2,0
+4.8,3.1,1.6,.2,0
+5.4,3.4,1.5,.4,0
+5.2,4.1,1.5,.1,0
+5.5,4.2,1.4,.2,0
+4.9,3.1,1.5,.2,0
+5,3.2,1.2,.2,0
+5.5,3.5,1.3,.2,0
+4.9,3.6,1.4,.1,0
+4.4,3,1.3,.2,0
+5.1,3.4,1.5,.2,0
+5,3.5,1.3,.3,0
+4.5,2.3,1.3,.3,0
+4.4,3.2,1.3,.2,0
+5,3.5,1.6,.6,0
+5.1,3.8,1.9,.4,0
+4.8,3,1.4,.3,0
+5.1,3.8,1.6,.2,0
+4.6,3.2,1.4,.2,0
+5.3,3.7,1.5,.2,0
+5,3.3,1.4,.2,0
+7,3.2,4.7,1.4,1
+6.4,3.2,4.5,1.5,1
+6.9,3.1,4.9,1.5,1
+5.5,2.3,4,1.3,1
+6.5,2.8,4.6,1.5,1
+5.7,2.8,4.5,1.3,1
+6.3,3.3,4.7,1.6,1
+4.9,2.4,3.3,1,1
+6.6,2.9,4.6,1.3,1
+5.2,2.7,3.9,1.4,1
+5,2,3.5,1,1
+5.9,3,4.2,1.5,1
+6,2.2,4,1,1
+6.1,2.9,4.7,1.4,1
+5.6,2.9,3.6,1.3,1
+6.7,3.1,4.4,1.4,1
+5.6,3,4.5,1.5,1
+5.8,2.7,4.1,1,1
+6.2,2.2,4.5,1.5,1
+5.6,2.5,3.9,1.1,1
+5.9,3.2,4.8,1.8,1
+6.1,2.8,4,1.3,1
+6.3,2.5,4.9,1.5,1
+6.1,2.8,4.7,1.2,1
+6.4,2.9,4.3,1.3,1
+6.6,3,4.4,1.4,1
+6.8,2.8,4.8,1.4,1
+6.7,3,5,1.7,1
+6,2.9,4.5,1.5,1
+5.7,2.6,3.5,1,1
+5.5,2.4,3.8,1.1,1
+5.5,2.4,3.7,1,1
+5.8,2.7,3.9,1.2,1
+6,2.7,5.1,1.6,1
+5.4,3,4.5,1.5,1
+6,3.4,4.5,1.6,1
+6.7,3.1,4.7,1.5,1
+6.3,2.3,4.4,1.3,1
+5.6,3,4.1,1.3,1
+5.5,2.5,4,1.3,1
+5.5,2.6,4.4,1.2,1
+6.1,3,4.6,1.4,1
+5.8,2.6,4,1.2,1
+5,2.3,3.3,1,1
+5.6,2.7,4.2,1.3,1
+5.7,3,4.2,1.2,1
+5.7,2.9,4.2,1.3,1
+6.2,2.9,4.3,1.3,1
+5.1,2.5,3,1.1,1
+5.7,2.8,4.1,1.3,1
+6.3,3.3,6,2.5,2
+5.8,2.7,5.1,1.9,2
+7.1,3,5.9,2.1,2
+6.3,2.9,5.6,1.8,2
+6.5,3,5.8,2.2,2
+7.6,3,6.6,2.1,2
+4.9,2.5,4.5,1.7,2
+7.3,2.9,6.3,1.8,2
+6.7,2.5,5.8,1.8,2
+7.2,3.6,6.1,2.5,2
+6.5,3.2,5.1,2,2
+6.4,2.7,5.3,1.9,2
+6.8,3,5.5,2.1,2
+5.7,2.5,5,2,2
+5.8,2.8,5.1,2.4,2
+6.4,3.2,5.3,2.3,2
+6.5,3,5.5,1.8,2
+7.7,3.8,6.7,2.2,2
+7.7,2.6,6.9,2.3,2
+6,2.2,5,1.5,2
+6.9,3.2,5.7,2.3,2
+5.6,2.8,4.9,2,2
+7.7,2.8,6.7,2,2
+6.3,2.7,4.9,1.8,2
+6.7,3.3,5.7,2.1,2
+7.2,3.2,6,1.8,2
+6.2,2.8,4.8,1.8,2
+6.1,3,4.9,1.8,2
+6.4,2.8,5.6,2.1,2
+7.2,3,5.8,1.6,2
+7.4,2.8,6.1,1.9,2
+7.9,3.8,6.4,2,2
+6.4,2.8,5.6,2.2,2
+6.3,2.8,5.1,1.5,2
+6.1,2.6,5.6,1.4,2
+7.7,3,6.1,2.3,2
+6.3,3.4,5.6,2.4,2
+6.4,3.1,5.5,1.8,2
+6,3,4.8,1.8,2
+6.9,3.1,5.4,2.1,2
+6.7,3.1,5.6,2.4,2
+6.9,3.1,5.1,2.3,2
+5.8,2.7,5.1,1.9,2
+6.8,3.2,5.9,2.3,2
+6.7,3.3,5.7,2.5,2
+6.7,3,5.2,2.3,2
+6.3,2.5,5,1.9,2
+6.5,3,5.2,2,2
+6.2,3.4,5.4,2.3,2
+5.9,3,5.1,1.8,2

machine_learning/neural_network.cpp

@@ -0,0 +1,790 @@
+/**
+ * @file 
+ * @author [Deep Raval](https://github.com/imdeep2905)
+ * 
+ * @brief Implementation of [Multilayer Perceptron] (https://en.wikipedia.org/wiki/Multilayer_perceptron). 
+ * 
+ * @details
+ * A multilayer perceptron (MLP) is a class of feedforward artificial neural network (ANN). The term MLP is used ambiguously, 
+ * sometimes loosely to any feedforward ANN, sometimes strictly to refer to networks composed of multiple layers of perceptrons 
+ * (with threshold activation). Multilayer perceptrons are sometimes colloquially referred to as "vanilla" neural networks, 
+ * especially when they have a single hidden layer.
+ * 
+ * An MLP consists of at least three layers of nodes: an input layer, a hidden layer and an output layer. Except for the 
+ * input nodes, each node is a neuron that uses a nonlinear activation function. MLP utilizes a supervised learning technique 
+ * called backpropagation for training. Its multiple layers and non-linear activation distinguish MLP from a linear 
+ * perceptron. It can distinguish data that is not linearly separable.
+ * 
+ * See [Backpropagation](https://en.wikipedia.org/wiki/Backpropagation) for training algorithm.
+ * 
+ * \note This implementation uses mini-batch gradient descent as optimizer and MSE as loss function. Bias is also not included.
+ */
+
+#include "vector_ops.hpp" // Custom header file for vector operations
+
+#include <iostream>
+#include <valarray>
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <chrono>
+#include <string>
+#include <fstream>
+#include <sstream>
+#include <cassert>
+
+/** \namespace machine_learning
+ * \brief Machine learning algorithms
+ */
+namespace machine_learning {
+    /** \namespace neural_network
+     * \brief Neural Network or Multilayer Perceptron
+     */   
+    namespace neural_network {
+        /** \namespace activations
+         * \brief Various activation functions used in Neural network
+         */         
+        namespace activations {
+            /**
+             * Sigmoid function
+             * @param X Value 
+             * @return Returns sigmoid(x)
+             */  
+            double sigmoid (const double &x) {
+                return 1.0 / (1.0 + std::exp(-x));
+            }
+
+            /**
+             * Derivative of sigmoid function
+             * @param X Value 
+             * @return Returns derivative of sigmoid(x)
+             */  
+            double dsigmoid (const double &x) {
+                return x * (1 - x);
+            }
+
+            /**
+             * Relu function
+             * @param X Value 
+             * @returns relu(x)
+             */  
+            double relu (const double &x) {
+                return std::max(0.0, x);
+            }
+
+            /**
+             * Derivative of relu function
+             * @param X Value 
+             * @returns derivative of relu(x)
+             */  
+            double drelu (const double &x) {
+                return x >= 0.0 ? 1.0 : 0.0;
+            }
+
+            /**
+             * Tanh function
+             * @param X Value 
+             * @return Returns tanh(x)
+             */  
+            double tanh (const double &x) {
+                return 2 / (1 + std::exp(-2 * x)) - 1;
+            }
+
+            /**
+             * Derivative of Sigmoid function
+             * @param X Value 
+             * @return Returns derivative of tanh(x)
+             */  
+            double dtanh (const double &x) {
+                return 1 - x * x;
+            }
+        } // namespace activations
+        /** \namespace util_functions
+         * \brief Various utility functions used in Neural network
+         */      
+        namespace util_functions {
+            /**
+             * Square function
+             * @param X Value 
+             * @return Returns x * x
+             */  
+            double square(const double &x) {
+                return x * x;
+            }
+            /**
+             * Identity function
+             * @param X Value 
+             * @return Returns x 
+             */  
+            double identity_function(const double &x) {
+                return x;
+            }
+        } // namespace util_functions
+        /** \namespace layers
+         * \brief This namespace contains layers used 
+         * in MLP.
+         */   
+        namespace layers {
+            /**
+             * neural_network::layers::DenseLayer class is used to store all necessary information about 
+             * the layers (i.e. neurons, activation and kernal). This class 
+             * is used by NeuralNetwork class to store layers.
+             * 
+             */
+            class DenseLayer {
+                public:
+                    // To store activation function and it's derivative
+                    double (*activation_function)(const double &); 
+                    double (*dactivation_function)(const double &);
+                    int neurons; // To store number of neurons (used in summary)
+                    std::string activation; // To store activation name (used in summary)
+                    std::vector <std::valarray <double>> kernal; // To store kernal (aka weights)
+                    
+                    /**
+                     * Constructor for neural_network::layers::DenseLayer class
+                     * @param neurons number of neurons
+                     * @param activation activation function for layer
+                     * @param kernal_shape shape of kernal 
+                     * @param random_kernal flag for whether to intialize kernal randomly
+                     */
+                    DenseLayer(const int &neurons, 
+                          const std::string &activation,
+                          const std::pair<size_t, size_t> &kernal_shape,
+                          const bool &random_kernal) {
+                        // Choosing activation (and it's derivative)
+                        if (activation == "sigmoid") {
+                            activation_function = neural_network::activations::sigmoid;
+                            dactivation_function = neural_network::activations::sigmoid;
+                        }
+                        else if (activation == "relu") {
+                            activation_function = neural_network::activations::relu;
+                            dactivation_function = neural_network::activations::drelu;
+                        }
+                        else if (activation == "tanh") {
+                            activation_function = neural_network::activations::tanh;
+                            dactivation_function = neural_network::activations::dtanh;
+                        }
+                        else if (activation == "none") {
+                            // Set identity function in casse of none is supplied
+                            activation_function = neural_network::util_functions::identity_function;
+                            dactivation_function = neural_network::util_functions::identity_function;
+                        }
+                        else {
+                            // If supplied activation is invalid
+                            std::cerr << "ERROR: Invalid argument for layer -> constructor -> activation, ";
+                            std::cerr << "Expected from {none, sigmoid, relu, tanh} got ";
+                            std::cerr << activation << std::endl;
+                            std::exit(EXIT_FAILURE);
+                        }
+                        this -> activation = activation; // Setting activation name
+                        this -> neurons = neurons; // Setting number of neurons
+                        // Initialize kernal according to flag
+                        if(random_kernal) {
+                            uniform_random_initialization(kernal, kernal_shape, -1.0, 1.0);
+                        }
+                        else {
+                            unit_matrix_initialization(kernal, kernal_shape);
+                        }
+                    }
+                    /**
+                     * Constructor for neural_network::layers::DenseLayer class
+                     * @param neurons number of neurons
+                     * @param activation activation function for layer
+                     * @param kernal values of kernal (useful in loading model)
+                     */
+                    DenseLayer (const int &neurons, 
+                           const std::string &activation, 
+                           const std::vector <std::valarray<double>> &kernal) {
+                        // Choosing activation (and it's derivative)
+                        if (activation == "sigmoid") {
+                            activation_function = neural_network::activations::sigmoid;
+                            dactivation_function = neural_network::activations::sigmoid;
+                        }
+                        else if (activation == "relu") {
+                            activation_function = neural_network::activations::relu;
+                            dactivation_function = neural_network::activations::drelu;
+                        }
+                        else if (activation == "tanh") {
+                            activation_function = neural_network::activations::tanh;
+                            dactivation_function = neural_network::activations::dtanh;
+                        }
+                        else if (activation == "none") {
+                            // Set identity function in casse of none is supplied
+                            activation_function = neural_network::util_functions::identity_function;
+                            dactivation_function = neural_network::util_functions::identity_function;
+                        }
+                        else {
+                            // If supplied activation is invalid
+                            std::cerr << "ERROR: Invalid argument for layer -> constructor -> activation, ";
+                            std::cerr << "Expected from {none, sigmoid, relu, tanh} got ";
+                            std::cerr << activation << std::endl;
+                            std::exit(EXIT_FAILURE);
+                        }
+                        this -> activation = activation; // Setting activation name
+                        this -> neurons = neurons; // Setting number of neurons
+                        this -> kernal = kernal;  // Setting supplied kernal values                 
+                    }
+                                       
+                /**
+                 * Copy Constructor for class DenseLayer. 
+                 * 
+                 * @param model instance of class to be copied.
+                 */
+                DenseLayer(const DenseLayer &layer) = default;   
+
+                /**
+                 * Destructor for class DenseLayer. 
+                 */                   
+                ~DenseLayer() = default;
+
+                /**
+                 * Copy assignment operator for class DenseLayer
+                 */
+                DenseLayer& operator = (const DenseLayer &layer) = default;
+
+                /**
+                 * Move constructor for class DenseLayer
+                 */
+                DenseLayer(DenseLayer &&) = default;
+
+                /**
+                 * Move assignment operator for class DenseLayer
+                 */
+                DenseLayer& operator = (DenseLayer &&) = default;
+            };
+        } // namespace layers
+        /**
+         * NeuralNetwork class is implements MLP. This class is
+         * used by actual user to create and train networks.
+         * 
+         */
+        class NeuralNetwork {
+            private:
+                std::vector <neural_network::layers::DenseLayer> layers; // To store layers
+                /**
+                 * Private Constructor for class NeuralNetwork. This constructor
+                 * is used internally to load model.
+                 * @param config vector containing pair (neurons, activation)
+                 * @param kernals vector containing all pretrained kernals
+                 */                
+                NeuralNetwork(const std::vector <std::pair<int, std::string>> &config,
+                              const std::vector <std::vector<std::valarray<double>>> &kernals) {
+                    // First layer should not have activation
+                    if(config.begin() -> second != "none") {
+                        std::cerr << "ERROR: First layer can't have activation other than none";
+                        std::cerr << std::endl;
+                        std::exit(EXIT_FAILURE);
+                    }
+                    // Network should have atleast two layers
+                    if(config.size() <= 1) {
+                        std::cerr << "ERROR: Invalid size of network, ";
+                        std::cerr << "Atleast two layers are required";
+                        std::exit(EXIT_FAILURE);
+                    }
+                    // Reconstructing all pretrained layers
+                    for(size_t i = 0; i < config.size(); i++) {
+                        layers.emplace_back(neural_network::layers::DenseLayer(config[i].first, 
+                                               config[i].second,
+                                               kernals[i])); 
+                    } 
+                    std::cout << "INFO: Network constructed successfully" << std::endl;                  
+                }
+                /**
+                 * Private function to get detailed predictions (i.e.
+                 * activated neuron values). This function is used in
+                 * backpropagation, single predict and batch predict.
+                 * @param X input vector
+                 */
+                std::vector<std::vector<std::valarray <double>>> 
+                __detailed_single_prediction (const std::vector<std::valarray <double>> &X) {
+                    std::vector<std::vector < std::valarray <double> >> details;
+                    std::vector < std::valarray <double> > current_pass = X;
+                    details.emplace_back(X);
+                    for(const auto &l : layers) {
+                        current_pass = multiply(current_pass, l.kernal);
+                        current_pass = apply_function(current_pass, l.activation_function);
+                        details.emplace_back(current_pass);
+                    }
+                    return details;
+                }
+            public:
+                /**
+                 * Default Constructor for class NeuralNetwork. This constructor
+                 * is used to create empty variable of type NeuralNetwork class.
+                 */   
+                NeuralNetwork() = default;
+
+                /**
+                 * Constructor for class NeuralNetwork. This constructor
+                 * is used by user.
+                 * @param config vector containing pair (neurons, activation)
+                 */   
+                explicit NeuralNetwork(const std::vector <std::pair<int, std::string>> &config) {
+                    // First layer should not have activation
+                    if(config.begin() -> second != "none") {
+                        std::cerr << "ERROR: First layer can't have activation other than none";
+                        std::cerr << std::endl;
+                        std::exit(EXIT_FAILURE);
+                    }
+                    // Network should have atleast two layers
+                    if(config.size() <= 1) {
+                        std::cerr << "ERROR: Invalid size of network, ";
+                        std::cerr << "Atleast two layers are required";
+                        std::exit(EXIT_FAILURE);
+                    }
+                    // Separately creating first layer so it can have unit matrix 
+                    // as kernal.
+                    layers.push_back(neural_network::layers::DenseLayer(config[0].first, 
+                                           config[0].second, 
+                                           {config[0].first, config[0].first},
+                                           false));
+                    // Creating remaining layers
+                    for(size_t i = 1; i < config.size(); i++) {
+                        layers.push_back(neural_network::layers::DenseLayer(config[i].first, 
+                                               config[i].second,
+                                               {config[i - 1].first, config[i].first},
+                                               true));
+                    }
+                    std::cout << "INFO: Network constructed successfully" << std::endl;
+                }
+
+                /**
+                 * Copy Constructor for class NeuralNetwork. 
+                 * 
+                 * @param model instance of class to be copied.
+                 */
+                NeuralNetwork(const NeuralNetwork &model) = default;   
+
+                /**
+                 * Destructor for class NeuralNetwork. 
+                 */                   
+                ~NeuralNetwork() = default;
+
+                /**
+                 * Copy assignment operator for class NeuralNetwork
+                 */
+                NeuralNetwork& operator = (const NeuralNetwork &model) = default;
+
+                /**
+                 * Move constructor for class NeuralNetwork
+                 */
+                NeuralNetwork(NeuralNetwork &&) = default;
+
+                /**
+                 * Move assignment operator for class NeuralNetwork
+                 */
+                NeuralNetwork& operator = (NeuralNetwork &&) = default;
+
+                /**
+                 * Function to get X and Y from csv file (where X = data, Y = label)
+                 * @param file_name csv file name
+                 * @param last_label flag for whether label is in first or last column
+                 * @param normalize flag for whether to normalize data 
+                 * @param slip_lines number of lines to skip
+                 * @return returns pair of X and Y
+                 */  
+                std::pair<std::vector<std::vector<std::valarray<double>>>, std::vector<std::vector<std::valarray<double>>>>
+                get_XY_from_csv(const std::string &file_name, 
+                                     const bool &last_label, 
+                                     const bool &normalize,
+                                     const int &slip_lines = 1) {
+                    std::ifstream in_file; // Ifstream to read file
+                    in_file.open(file_name.c_str(), std::ios::in); // Open file
+                    std::vector <std::vector<std::valarray<double>>> X, Y; // To store X and Y
+                    std::string line; // To store each line
+                    // Skip lines
+                    for(int i = 0; i < slip_lines; i ++) {
+                        std::getline(in_file, line, '\n'); // Ignore line
+                    }
+                    // While file has information
+                    while(!in_file.eof() && std::getline(in_file, line, '\n'))
+                    {
+                        std::valarray <double> x_data, y_data; // To store single sample and label
+                        std::stringstream ss(line); // Constructing stringstream from line
+                        std::string token; // To store each token in line (seprated by ',')
+                        while(std::getline(ss, token, ',')) { // For each token
+                            // Insert numerical value of token in x_data
+                            x_data = insert_element(x_data, std::stod(token));
+                        } 
+                        // If label is in last column
+                        if(last_label) {
+                            y_data.resize(this -> layers.back().neurons);
+                            // If task is classification
+                            if(y_data.size() > 1) {
+                                y_data[x_data[x_data.size() - 1]] = 1;
+                            }
+                            // If task is regrssion (of single value)
+                            else {
+                                y_data[0] = x_data[x_data.size() - 1];
+                            } 
+                            x_data = pop_back(x_data); // Remove label from x_data
+                        }
+                        else {
+                            y_data.resize(this -> layers.back().neurons);
+                            // If task is classification
+                            if(y_data.size() > 1) {
+                                y_data[x_data[x_data.size() - 1]] = 1;
+                            }
+                            // If task is regrssion (of single value)
+                            else {
+                                y_data[0] = x_data[x_data.size() - 1];
+                            } 
+                            x_data = pop_front(x_data); // Remove label from x_data
+                        }
+                        // Push collected X_data and y_data in X and Y
+                        X.push_back({x_data});
+                        Y.push_back({y_data});
+                    }
+                    in_file.close();
+                    // Normalize training data if flag is set
+                    if(normalize) {
+                        // Scale data between 0 and 1 using min-max scaler
+                        X = minmax_scaler(X, 0.01, 1.0);
+                    }
+                    return make_pair(X, Y); // Return pair of X and Y
+                }
+
+                /**
+                 * Function to get prediction of model on single sample.
+                 * @param X array of feature vectors
+                 * @return returns predictions as vector
+                 */  
+                std::vector<std::valarray <double>> 
+                single_predict (const std::vector<std::valarray <double>> &X) {
+                    // Get activations of all layers
+                    auto activations = this -> __detailed_single_prediction(X);
+                    // Return activations of last layer (actual predicted values)
+                    return activations.back();
+                }
+
+                /**
+                 * Function to get prediction of model on batch
+                 * @param X array of feature vectors
+                 * @return returns predicted values as vector
+                 */  
+                std::vector < std::vector <std::valarray<double>>>
+                batch_predict (const std::vector <std::vector <std::valarray <double>>> &X) {
+                    // Store predicted values
+                    std::vector < std::vector <std::valarray<double>>> predicted_batch(X.size());
+                    for(size_t i = 0; i < X.size(); i++) { // For every sample
+                        // Push predicted values
+                        predicted_batch[i] = this -> single_predict(X[i]);
+                    }
+                    return predicted_batch; // Return predicted values
+                }
+
+                /**
+                 * Function to fit model on supplied data
+                 * @param X array of feature vectors
+                 * @param Y array of target values
+                 * @param epochs number of epochs (default = 100)
+                 * @param learning_rate learning rate (default = 0.01)
+                 * @param batch_size batch size for gradient descent (default = 32)
+                 * @param shuffle flag for whether to shuffle data (default = true)
+                 */  
+                void fit(const std::vector < std::vector <std::valarray<double>>>  &X_, 
+                         const std::vector < std::vector <std::valarray<double>>>  &Y_, 
+                         const int &epochs = 100, 
+                         const double &learning_rate = 0.01,
+                         const size_t &batch_size = 32,
+                         const bool &shuffle = true) {
+                    std::vector < std::vector <std::valarray<double>>> X = X_, Y = Y_;
+                    // Both label and input data should have same size
+                    if (X.size() != Y.size()) {
+                        std::cerr << "ERROR : X and Y in fit have different sizes" << std::endl;
+                        std::exit(EXIT_FAILURE);
+                    }
+                    std::cout << "INFO: Training Started" << std::endl;
+                    for (int epoch = 1; epoch <= epochs; epoch++) { // For every epoch
+                        // Shuffle X and Y if flag is set
+                        if(shuffle) {
+                            equal_shuffle(X, Y);
+                        }
+                        auto start = std::chrono::high_resolution_clock::now(); // Start clock
+                        double loss = 0, acc = 0; // Intialize performance metrics with zero
+                        // For each starting index of batch
+                        for(size_t batch_start = 0; batch_start < X.size(); batch_start += batch_size) {
+                            for(size_t i = batch_start; i < std::min(X.size(), batch_start + batch_size); i++) {
+                                    std::vector <std::valarray<double>> grad, cur_error, predicted;
+                                    auto activations = this -> __detailed_single_prediction(X[i]);
+                                    // Gradients vector to store gradients for all layers
+                                    // They will be averaged and applied to kernal
+                                    std::vector<std::vector<std::valarray<double>>> gradients;
+                                    gradients.resize(this -> layers.size());
+                                    // First intialize gradients to zero
+                                    for(size_t i = 0; i < gradients.size(); i++) {
+                                        zeroes_initialization(gradients[i], get_shape(this -> layers[i].kernal));
+                                    }
+                                    predicted = activations.back(); // Predicted vector
+                                    cur_error = predicted - Y[i]; // Absoulute error
+                                    // Calculating loss with MSE
+                                    loss += sum(apply_function(cur_error, neural_network::util_functions::square));
+                                    // If prediction is correct
+                                    if(argmax(predicted) == argmax(Y[i])) {
+                                        acc += 1;
+                                    }
+                                    // For every layer (except first) starting from last one
+                                    for(size_t j = this -> layers.size() - 1; j >= 1; j--) {
+                                        // Backpropogating errors
+                                        cur_error = hadamard_product(cur_error, 
+                                                                    apply_function(activations[j + 1], 
+                                                                    this -> layers[j].dactivation_function));
+                                        // Calculating gradient for current layer
+                                        grad = multiply(transpose(activations[j]), cur_error);
+                                        // Change error according to current kernal values
+                                        cur_error = multiply(cur_error, transpose(this -> layers[j].kernal));
+                                        // Adding gradient values to collection of gradients
+                                        gradients[j] = gradients[j] + grad / double(batch_size);
+                                    }
+                                    // Applying gradients
+                                    for(size_t j = this -> layers.size() - 1; j >= 1; j--) {
+                                        // Updating kernal (aka weights)
+                                        this -> layers[j].kernal = this -> layers[j].kernal -
+                                                            gradients[j] * learning_rate;
+                                    }
+                                }
+                        }
+                        auto stop = std::chrono::high_resolution_clock::now(); // Stoping the clock
+                        // Calculate time taken by epoch
+                        auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start);
+                        loss /= X.size(); // Averaging loss
+                        acc /= X.size(); // Averaging accuracy
+                        std::cout.precision(4); // set output precision to 4 
+                        // Printing training stats
+                        std::cout << "Training: Epoch " << epoch << '/' << epochs; 
+                        std::cout << ", Loss: " << loss;
+                        std::cout << ", Accuracy: " << acc;
+                        std::cout << ", Taken time: " << duration.count() / 1e6 << " seconds";
+                        std::cout << std::endl;
+                    }
+                    return;
+                }
+
+                /**
+                 * Function to fit model on data stored in csv file
+                 * @param file_name csv file name
+                 * @param last_label flag for whether label is in first or last column
+                 * @param epochs number of epochs
+                 * @param learning_rate learning rate 
+                 * @param normalize flag for whether to normalize data 
+                 * @param slip_lines number of lines to skip
+                 * @param batch_size batch size for gradient descent (default = 32)
+                 * @param shuffle flag for whether to shuffle data (default = true)
+                 */                
+                void fit_from_csv (const std::string &file_name, 
+                                     const bool &last_label,
+                                     const int &epochs, 
+                                     const double &learning_rate,
+                                     const bool &normalize,
+                                     const int &slip_lines = 1,
+                                     const size_t &batch_size = 32,
+                                     const bool &shuffle = true) {
+                    // Getting training data from csv file                    
+                    auto data = this -> get_XY_from_csv(file_name, last_label, normalize, slip_lines);
+                    // Fit the model on training data
+                    this -> fit(data.first, data.second, epochs, learning_rate, batch_size, shuffle);
+                    return;
+                }
+
+                /**
+                 * Function to evaluate model on supplied data
+                 * @param X array of feature vectors (input data)
+                 * @param Y array of target values (label)
+                 */  
+                void evaluate(const std::vector< std::vector <std::valarray <double>>> &X, 
+                                const std::vector< std::vector <std::valarray <double>>> &Y) {
+                    std::cout << "INFO: Evaluation Started" << std::endl;
+                    double acc = 0, loss = 0; // intialize performance metrics with zero
+                    for(size_t i = 0; i < X.size(); i++) { // For every sample in input
+                        // Get predictions
+                        std::vector<std::valarray<double>> pred = this -> single_predict(X[i]);
+                        // If predicted class is correct 
+                        if(argmax(pred) == argmax(Y[i])) {
+                            acc += 1; // Increment accuracy
+                        }
+                        // Calculating loss - Mean Squared Error
+                        loss += sum(apply_function((Y[i] - pred), 
+                                    neural_network::util_functions::square) * 0.5);
+                    }
+                    acc /= X.size(); // Averaging accuracy
+                    loss /= X.size(); // Averaging loss
+                    // Prinitng performance of the model
+                    std::cout << "Evaluation: Loss: " << loss;
+                    std::cout << ", Accuracy: " << acc << std::endl;
+                    return;
+                }
+
+                /**
+                 * Function to evaluate model on data stored in csv file
+                 * @param file_name csv file name
+                 * @param last_label flag for whether label is in first or last column
+                 * @param normalize flag for whether to normalize data 
+                 * @param slip_lines number of lines to skip
+                 */                
+                void evaluate_from_csv (const std::string &file_name, 
+                                     const bool &last_label,
+                                     const bool &normalize,
+                                     const int &slip_lines = 1) {
+                    // Getting training data from csv file
+                    auto data = this -> get_XY_from_csv(file_name, last_label, normalize, slip_lines);
+                    // Evaluating model
+                    this -> evaluate(data.first, data.second);
+                    return;
+                }                
+
+                /**
+                 * Function to save current model.
+                 * @param file_name file name to save model (*.model)
+                 */
+                void save_model (const std::string &_file_name) {
+                    std::string file_name = _file_name;
+                    // Adding ".model" extension if it is not already there in name
+                    if(file_name.find(".model") == file_name.npos) {
+                        file_name += ".model";
+                    }
+                    std::ofstream out_file; // Ofstream to write in file
+                    // Open file in out|trunc mode
+                    out_file.open(file_name.c_str(), std::ofstream::out | std::ofstream::trunc);
+                    /**
+                        Format in which model is saved:
+
+                        total_layers
+                        neurons(1st neural_network::layers::DenseLayer) activation_name(1st neural_network::layers::DenseLayer)
+                        kernal_shape(1st neural_network::layers::DenseLayer)
+                        kernal_values
+                        .
+                        .
+                        .
+                        neurons(Nth neural_network::layers::DenseLayer) activation_name(Nth neural_network::layers::DenseLayer)
+                        kernal_shape(Nth neural_network::layers::DenseLayer)
+                        kernal_value
+
+                        For Example, pretrained model with 3 layers:
+                        <pre>
+                        3
+                        4 none
+                        4 4
+                        1 0 0 0 
+                        0 1 0 0 
+                        0 0 1 0 
+                        0 0 0 1 
+                        6 relu
+                        4 6
+                        -1.88963 -3.61165 1.30757 -0.443906 -2.41039 -2.69653 
+                        -0.684753 0.0891452 0.795294 -2.39619 2.73377 0.318202 
+                        -2.91451 -4.43249 -0.804187 2.51995 -6.97524 -1.07049 
+                        -0.571531 -1.81689 -1.24485 1.92264 -2.81322 1.01741 
+                        3 sigmoid
+                        6 3
+                        0.390267 -0.391703 -0.0989607 
+                        0.499234 -0.564539 -0.28097 
+                        0.553386 -0.153974 -1.92493 
+                        -2.01336 -0.0219682 1.44145 
+                        1.72853 -0.465264 -0.705373 
+                        -0.908409 -0.740547 0.376416 
+                        </pre>
+                    */
+                    // Saving model in the same format
+                    out_file << layers.size();
+                    out_file << std::endl;
+                    for(const auto &layer : this -> layers) {
+                        out_file << layer.neurons << ' ' << layer.activation << std::endl;
+                        const auto shape = get_shape(layer.kernal);
+                        out_file << shape.first << ' ' << shape.second << std::endl;
+                        for(const auto &row : layer.kernal) {
+                            for(const auto &val : row) {
+                                out_file << val << ' ';
+                            }
+                            out_file << std::endl;
+                        }
+                    }
+                    std::cout << "INFO: Model saved successfully with name : ";
+                    std::cout << file_name << std::endl;
+                    return;
+                }
+
+                /**
+                 * Function to load earlier saved model.
+                 * @param file_name file from which model will be loaded (*.model)
+                 * @return instance of NeuralNetwork class with pretrained weights
+                 */
+                NeuralNetwork load_model (const std::string &file_name) {
+                    std::ifstream in_file; // Ifstream to read file
+                    in_file.open(file_name.c_str()); // Openinig file
+                    std::vector <std::pair<int, std::string>> config; // To store config
+                    std::vector <std::vector<std::valarray<double>>> kernals; // To store pretrained kernals
+                    // Loading model from saved file format
+                    size_t total_layers = 0; 
+                    in_file >> total_layers;
+                    for(size_t i = 0; i < total_layers; i++) {
+                        int neurons = 0;
+                        std::string activation;
+                        size_t shape_a = 0, shape_b = 0;
+                        std::vector<std::valarray<double>> kernal;
+                        in_file >> neurons >> activation >> shape_a >> shape_b;
+                        for(size_t r = 0; r < shape_a; r++) {
+                            std::valarray<double> row(shape_b);
+                            for(size_t c = 0; c < shape_b; c++) {
+                                in_file >> row[c];
+                            }
+                            kernal.push_back(row);
+                        }
+                        config.emplace_back(make_pair(neurons, activation));;
+                        kernals.emplace_back(kernal);
+                    }
+                    std::cout << "INFO: Model loaded successfully" << std::endl;
+                    return NeuralNetwork(config, kernals); // Return instance of NeuralNetwork class
+                }
+
+                /**
+                 * Function to print summary of the network.
+                 */
+                void summary () {
+                    // Printing Summary 
+                    std::cout << "===============================================================" << std::endl;
+                    std::cout << "\t\t+ MODEL SUMMARY +\t\t\n";
+                    std::cout << "===============================================================" << std::endl;
+                    for(size_t i = 1; i <= layers.size(); i++) { // For every layer
+                        std::cout << i << ")";
+                        std::cout << " Neurons : " << layers[i - 1].neurons; // number of neurons
+                        std::cout << ", Activation : " << layers[i - 1].activation; // activation
+                        std::cout << ", Kernal Shape : " << get_shape(layers[i - 1].kernal); // kernal shape
+                        std::cout << std::endl;
+                    }
+                    std::cout << "===============================================================" << std::endl;
+                    return;
+                }
+        
+        };
+    } // namespace neural_network
+} // namespace machine_learning
+
+/**
+ * Function to test neural network
+ * @returns none
+ */
+static void test() {
+    // Creating network with 3 layers for "iris.csv"
+    machine_learning::neural_network::NeuralNetwork myNN =
+    machine_learning::neural_network::NeuralNetwork({
+        {4, "none"}, // First layer with 3 neurons and "none" as activation
+        {6, "relu"}, // Second layer with 6 neurons and "relu" as activation
+        {3, "sigmoid"} // Third layer with 3 neurons and "sigmoid" as activation
+    });
+    // Printing summary of model
+    myNN.summary();
+    // Training Model
+    myNN.fit_from_csv("iris.csv", true, 100, 0.3, false, 2, 32, true);
+    // Testing predictions of model
+    assert(machine_learning::argmax(myNN.single_predict({{5,3.4,1.6,0.4}})) == 0);
+    assert(machine_learning::argmax(myNN.single_predict({{6.4,2.9,4.3,1.3}})) == 1);
+    assert(machine_learning::argmax(myNN.single_predict({{6.2,3.4,5.4,2.3}})) == 2);
+    return;
+}
+
+/** Driver Code */
+int main() {
+    // Testing
+    test();
+    return 0;
+}

machine_learning/vector_ops.hpp

@@ -0,0 +1,484 @@
+/**
+ * @file vector_ops.hpp
+ * @author [Deep Raval](https://github.com/imdeep2905)
+ * 
+ * @brief Various functions for vectors associated with [NeuralNetwork (aka Multilayer Perceptron)] 
+ * (https://en.wikipedia.org/wiki/Multilayer_perceptron).
+ * 
+ */
+#ifndef VECTOR_OPS_FOR_NN
+#define VECTOR_OPS_FOR_NN
+
+#include <iostream>
+#include <algorithm>
+#include <vector>
+#include <valarray>
+#include <chrono>
+#include <random>
+
+/**
+ * @namespace machine_learning
+ * @brief Machine Learning algorithms
+ */
+namespace machine_learning {
+/**
+ * Overloaded operator "<<" to print 2D vector
+ * @tparam T typename of the vector
+ * @param out std::ostream to output
+ * @param A 2D vector to be printed
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out,
+                         std::vector<std::valarray<T>> const &A) {
+    // Setting output precision to 4 in case of floating point numbers
+    out.precision(4); 
+    for(const auto &a : A) { // For each row in A
+        for(const auto &x : a) { // For each element in row
+            std::cerr << x << ' '; // print element 
+        }
+        std::cerr << std::endl;
+    }
+    return out;
+}
+
+/**
+ * Overloaded operator "<<" to print a pair 
+ * @tparam T typename of the pair
+ * @param out std::ostream to output
+ * @param A Pair to be printed
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out, const std::pair<T, T> &A) {
+    // Setting output precision to 4 in case of floating point numbers
+    out.precision(4);
+    // printing pair in the form (p, q)
+    std::cerr << "(" << A.first << ", " << A.second << ")";
+    return out;
+}
+
+/**
+ * Overloaded operator "<<" to print a 1D vector
+ * @tparam T typename of the vector
+ * @param out std::ostream to output
+ * @param A 1D vector to be printed
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out, const std::valarray<T> &A) {
+    // Setting output precision to 4 in case of floating point numbers
+    out.precision(4);
+    for(const auto &a : A) { // For every element in the vector.
+        std::cerr << a << ' '; // Print element
+    }
+    std::cerr << std::endl;
+    return out;
+}
+
+/**
+ * Function to insert element into 1D vector
+ * @tparam T typename of the 1D vector and the element
+ * @param A 1D vector in which element will to be inserted
+ * @param ele element to be inserted
+ * @return new resultant vector
+ */
+template <typename T>
+std::valarray<T> insert_element(const std::valarray <T> &A, const T &ele) {
+    std::valarray <T> B; // New 1D vector to store resultant vector
+    B.resize(A.size() + 1); // Resizing it accordingly
+    for(size_t i = 0; i < A.size(); i++) { // For every element in A
+        B[i] = A[i]; // Copy element in B
+    }
+    B[B.size() - 1] = ele; // Inserting new element in last position
+    return B; // Return resultant vector
+}
+
+/**
+ * Function to remove first element from 1D vector
+ * @tparam T typename of the vector
+ * @param A 1D vector from which first element will be removed
+ * @return new resultant vector
+ */
+template <typename T>
+std::valarray <T> pop_front(const std::valarray<T> &A) {
+    std::valarray <T> B; // New 1D vector to store resultant vector
+    B.resize(A.size() - 1); // Resizing it accordingly
+    for(size_t i = 1; i < A.size(); i ++) { // // For every (except first) element in A 
+        B[i - 1] = A[i]; // Copy element in B with left shifted position
+    }
+    return B; // Return resultant vector
+}
+
+/**
+ * Function to remove last element from 1D vector
+ * @tparam T typename of the vector
+ * @param A 1D vector from which last element will be removed
+ * @return new resultant vector
+ */
+template <typename T>
+std::valarray <T> pop_back(const std::valarray<T> &A) {
+    std::valarray <T> B; // New 1D vector to store resultant vector
+    B.resize(A.size() - 1); // Resizing it accordingly
+    for(size_t i = 0; i < A.size() - 1; i ++) { // For every (except last) element in A 
+        B[i] = A[i]; // Copy element in B
+    }
+    return B; // Return resultant vector
+}
+
+/**
+ * Function to equally shuffle two 3D vectors (used for shuffling training data)
+ * @tparam T typename of the vector
+ * @param A First 3D vector
+ * @param B Second 3D vector
+ */
+template <typename T>
+void equal_shuffle(std::vector < std::vector <std::valarray<T>> >  &A, 
+                   std::vector < std::vector <std::valarray<T>> >  &B) {
+    // If two vectors have different sizes
+    if(A.size() != B.size())
+    {
+        std::cerr << "ERROR : Can not equally shuffle two vectors with different sizes: ";
+        std::cerr << A.size() << " and " << B.size() << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    for(size_t i = 0; i < A.size(); i++) { // For every element in A and B
+        // Genrating random index < size of A and B
+        std::srand(std::chrono::system_clock::now().time_since_epoch().count());
+        size_t random_index = std::rand() % A.size();
+        // Swap elements in both A and B with same random index
+        std::swap(A[i], A[random_index]);
+        std::swap(B[i], B[random_index]);
+    }
+    return;
+}
+
+/**
+ * Function to initialize given 2D vector using uniform random initialization
+ * @tparam T typename of the vector
+ * @param A 2D vector to be initialized
+ * @param shape required shape
+ * @param low lower limit on value
+ * @param high upper limit on value
+ */
+template <typename T>
+void uniform_random_initialization(std::vector<std::valarray<T>> &A, 
+                   const std::pair<size_t, size_t> &shape, 
+                   const T &low, 
+                   const T &high) {
+    A.clear(); // Making A empty 
+    // Uniform distribution in range [low, high]
+    std::default_random_engine generator(std::chrono::system_clock::now().time_since_epoch().count());
+    std::uniform_real_distribution <T> distribution(low, high);
+    for(size_t i = 0; i < shape.first; i++) { // For every row 
+        std::valarray <T> row; // Making empty row which will be inserted in vector
+        row.resize(shape.second);
+        for(auto &r : row) { // For every element in row
+            r = distribution(generator); // copy random number 
+        }  
+        A.push_back(row); // Insert new row in vector
+    }
+    return;
+}
+
+
+/**
+ * Function to Intialize 2D vector as unit matrix
+ * @tparam T typename of the vector
+ * @param A 2D vector to be initialized
+ * @param shape required shape
+ */
+template <typename T>
+void unit_matrix_initialization(std::vector<std::valarray<T>> &A, 
+                   const std::pair<size_t, size_t> &shape
+                   ) {
+    A.clear(); // Making A empty 
+    for(size_t i = 0; i < shape.first; i++) {
+        std::valarray <T> row; // Making empty row which will be inserted in vector
+        row.resize(shape.second);
+        row[i] = T(1); // Insert 1 at ith position 
+        A.push_back(row); // Insert new row in vector
+    }
+    return;
+}
+
+/**
+ * Function to Intialize 2D vector as zeroes
+ * @tparam T typename of the vector
+ * @param A 2D vector to be initialized
+ * @param shape required shape
+ */
+template <typename T>
+void zeroes_initialization(std::vector<std::valarray<T>> &A, 
+                   const std::pair<size_t, size_t> &shape
+                   ) {
+    A.clear(); // Making A empty 
+    for(size_t i = 0; i < shape.first; i++) {
+        std::valarray <T> row; // Making empty row which will be inserted in vector
+        row.resize(shape.second); // By default all elements are zero
+        A.push_back(row); // Insert new row in vector
+    }
+    return;
+}
+
+/**
+ * Function to get sum of all elements in 2D vector
+ * @tparam T typename of the vector
+ * @param A 2D vector for which sum is required
+ * @return returns sum of all elements of 2D vector
+ */
+template <typename T>
+T sum(const std::vector<std::valarray<T>> &A) {
+    T cur_sum = 0; // Initially sum is zero
+    for(const auto &a : A) { // For every row in A
+        cur_sum += a.sum(); // Add sum of that row to current sum
+    }
+    return cur_sum; // Return sum
+}
+
+/**
+ * Function to get shape of given 2D vector
+ * @tparam T typename of the vector
+ * @param A 2D vector for which shape is required
+ * @return shape as pair
+ */
+template <typename T>
+std::pair<size_t, size_t> get_shape(const std::vector<std::valarray<T>> &A) {
+    const size_t sub_size = (*A.begin()).size();
+    for(const auto &a : A) {
+        // If supplied vector don't have same shape in all rows
+        if(a.size() != sub_size) {
+            std::cerr << "ERROR: (get_shape) Supplied vector is not 2D Matrix" << std::endl;
+            std::exit(EXIT_FAILURE);
+        }
+    }
+    return std::make_pair(A.size(), sub_size); // Return shape as pair
+}
+
+/**
+ * Function to scale given 3D vector using min-max scaler
+ * @tparam T typename of the vector
+ * @param A 3D vector which will be scaled
+ * @param low new minimum value 
+ * @param high new maximum value
+ * @return new scaled 3D vector
+ */
+template <typename T>
+std::vector<std::vector<std::valarray<T>>>
+minmax_scaler(const std::vector<std::vector<std::valarray<T>>> &A, const T &low, const T &high) {
+    std::vector<std::vector<std::valarray<T>>> B = A; // Copying into new vector B
+    const auto shape = get_shape(B[0]); // Storing shape of B's every element
+    // As this function is used for scaling training data vector should be of shape (1, X)
+    if(shape.first != 1) {
+        std::cerr << "ERROR: (MinMax Scaling) Supplied vector is not supported for minmax scaling, shape: ";
+        std::cerr << shape << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    for(size_t i = 0; i < shape.second; i++) {
+        T min = B[0][0][i], max = B[0][0][i]; 
+        for(size_t j = 0; j < B.size(); j++) {
+            // Updating minimum and maximum values
+            min = std::min(min, B[j][0][i]);
+            max = std::max(max, B[j][0][i]);
+        }
+        for(size_t j = 0; j < B.size(); j++) {
+            // Applying min-max scaler formula
+            B[j][0][i] = ((B[j][0][i] - min) / (max - min)) * (high - low) + low;
+        }
+    }
+    return B; // Return new resultant 3D vector
+}
+
+/**
+ * Function to get index of maximum element in 2D vector 
+ * @tparam T typename of the vector
+ * @param A 2D vector for which maximum index is required
+ * @return index of maximum element
+ */
+template <typename T>
+size_t argmax(const std::vector<std::valarray<T>> &A) {
+    const auto shape = get_shape(A);
+    // As this function is used on predicted (or target) vector, shape should be (1, X)    
+    if(shape.first != 1) {
+        std::cerr << "ERROR: (argmax) Supplied vector is ineligible for argmax" << std::endl;
+        std::exit(EXIT_FAILURE);        
+    }
+    // Return distance of max element from first element (i.e. index)
+    return std::distance(std::begin(A[0]), std::max_element(std::begin(A[0]), std::end(A[0])));
+}
+
+/**
+ * Function which applys supplied function to every element of 2D vector
+ * @tparam T typename of the vector
+ * @param A 2D vector on which function will be applied
+ * @param func Function to be applied
+ * @return new resultant vector
+ */
+template <typename T> 
+std::vector <std::valarray <T>> apply_function(const std::vector <std::valarray <T>> &A, 
+                                               T (*func) (const T &)) {
+    std::vector<std::valarray<double>> B = A; // New vector to store resultant vector
+    for(auto &b : B) { // For every row in vector
+        b = b.apply(func); // Apply function to that row
+    }
+    return B; // Return new resultant 2D vector
+}
+
+/**
+ * Overloaded operator "*" to multiply given 2D vector with scaler 
+ * @tparam T typename of both vector and the scaler
+ * @param A 2D vector to which scaler will be multiplied
+ * @param val Scaler value which will be multiplied
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T> > operator * (const std::vector<std::valarray<T>> &A, const T& val) {
+    std::vector<std::valarray<double>> B = A; // New vector to store resultant vector
+    for(auto &b : B) { // For every row in vector
+        b = b * val; // Multiply row with scaler
+    }
+    return B; // Return new resultant 2D vector
+}
+
+/**
+ * Overloaded operator "/" to divide given 2D vector with scaler 
+ * @tparam T typename of the vector and the scaler
+ * @param A 2D vector to which scaler will be divided
+ * @param val Scaler value which will be divided
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T> > operator / (const std::vector<std::valarray<T>> &A, const T& val) {
+    std::vector<std::valarray<double>> B = A; // New vector to store resultant vector
+    for(auto &b : B) { // For every row in vector
+        b = b / val; // Divide row with scaler
+    }
+    return B; // Return new resultant 2D vector
+}
+
+/**
+ * Function to get transpose of 2D vector
+ * @tparam T typename of the vector
+ * @param A 2D vector which will be transposed
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T>> transpose(const std::vector<std::valarray<T>> &A) {
+    const auto shape = get_shape(A); // Current shape of vector
+    std::vector <std::valarray <T> > B; // New vector to store result
+    // Storing transpose values of A in B
+    for(size_t j = 0; j < shape.second; j++) { 
+        std::valarray <T> row; 
+        row.resize(shape.first);
+        for(size_t i = 0; i < shape.first; i++) {
+            row[i] = A[i][j];
+        }
+        B.push_back(row);
+    }
+    return B; // Return new resultant 2D vector
+}
+
+/**
+ * Overloaded operator "+" to add two 2D vectors
+ * @tparam T typename of the vector
+ * @param A First 2D vector 
+ * @param B Second 2D vector
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T> > operator + (const std::vector<std::valarray<T>> &A, const std::vector<std::valarray<T>> &B) {
+    const auto shape_a = get_shape(A);
+    const auto shape_b = get_shape(B);
+    // If vectors don't have equal shape
+    if(shape_a.first != shape_b.first || shape_a.second != shape_b.second) {
+        std::cerr << "ERROR: (vector addition) Supplied vectors have different shapes ";
+        std::cerr << shape_a << " and " << shape_b << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    std::vector<std::valarray <T>> C;
+    for(size_t i = 0; i < A.size(); i++) { // For every row
+        C.push_back(A[i] + B[i]); // Elementwise addition
+    }
+    return C; // Return new resultant 2D vector
+}
+
+/**
+ * Overloaded operator "-" to add subtract 2D vectors
+ * @tparam T typename of the vector
+ * @param A First 2D vector 
+ * @param B Second 2D vector
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T>> operator - (const std::vector<std::valarray<T>> &A, const std::vector<std::valarray<T>> &B) {
+    const auto shape_a = get_shape(A);
+    const auto shape_b = get_shape(B);
+    // If vectors don't have equal shape
+    if(shape_a.first != shape_b.first || shape_a.second != shape_b.second) {
+        std::cerr << "ERROR: (vector subtraction) Supplied vectors have different shapes ";
+        std::cerr << shape_a << " and " << shape_b << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    std::vector<std::valarray<T>> C; // Vector to store result
+    for(size_t i = 0; i < A.size(); i++) { // For every row
+        C.push_back(A[i] - B[i]); // Elementwise substraction
+    }
+    return C; // Return new resultant 2D vector
+}
+
+/**
+ * Function to multiply two 2D vectors
+ * @tparam T typename of the vector
+ * @param A First 2D vector 
+ * @param B Second 2D vector
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T>> multiply(const std::vector<std::valarray<T>> &A, const std::vector<std::valarray<T>> &B) {
+    const auto shape_a = get_shape(A);
+    const auto shape_b = get_shape(B);
+    // If vectors are not eligible for multiplication
+    if(shape_a.second != shape_b.first ) {
+        std::cerr << "ERROR: (multiply) Supplied vectors are not eligible for multiplication ";
+        std::cerr << shape_a << " and " << shape_b << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    std::vector<std::valarray<T>> C; // Vector to store result
+    // Normal matrix multiplication 
+    for (size_t i = 0; i < shape_a.first; i++) {
+        std::valarray<T> row;
+        row.resize(shape_b.second);
+        for(size_t j = 0; j < shape_b.second; j++) {
+            for(size_t k = 0; k < shape_a.second; k++) {
+                row[j] += A[i][k] * B[k][j];
+            }
+        }
+        C.push_back(row);
+    }
+    return C; // Return new resultant 2D vector
+}
+
+/**
+ * Function to get hadamard product of two 2D vectors
+ * @tparam T typename of the vector
+ * @param A First 2D vector 
+ * @param B Second 2D vector
+ * @return new resultant vector
+ */
+template <typename T>
+std::vector <std::valarray <T>> hadamard_product(const std::vector<std::valarray<T>> &A, const std::vector<std::valarray<T>> &B) {
+    const auto shape_a = get_shape(A);
+    const auto shape_b = get_shape(B);
+    // If vectors are not eligible for hadamard product
+    if(shape_a.first != shape_b.first || shape_a.second != shape_b.second) {
+        std::cerr << "ERROR: (hadamard_product) Supplied vectors have different shapes ";
+        std::cerr << shape_a << " and " << shape_b << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    std::vector<std::valarray<T>> C; // Vector to store result
+    for(size_t i = 0; i < A.size(); i++) {
+        C.push_back(A[i] * B[i]); // Elementwise multiplication
+    }
+    return C; // Return new resultant 2D vector
+}
+}  // namespace machine_learning
+
+
+#endif

math/fibonacci.cpp

@@ -13,8 +13,10 @@
 
 /**
  * Recursively compute sequences
+ * @param n input
+ * @returns n-th element of the Fbinacci's sequence
  */
-unsigned int fibonacci(unsigned int n) {
+uint64_t fibonacci(uint64_t n) {
     /* If the input is 0 or 1 just return the same
        This will set the first 2 values of the sequence */
     if (n <= 1) {
@@ -31,27 +33,27 @@ unsigned int fibonacci(unsigned int n) {
  * @returns `void`
 */
 static void test() {
-    unsigned int test_case_1 = fibonacci(0);
+    uint64_t test_case_1 = fibonacci(0);
     assert(test_case_1 == 0);
     std::cout << "Passed Test 1!" << std::endl;
 
-    unsigned int test_case_2 = fibonacci(1);
+    uint64_t test_case_2 = fibonacci(1);
     assert(test_case_2 == 1);
     std::cout << "Passed Test 2!" << std::endl;
 
-    unsigned int test_case_3 = fibonacci(2);
+    uint64_t test_case_3 = fibonacci(2);
     assert(test_case_3 == 1);
     std::cout << "Passed Test 3!" << std::endl;
 
-    unsigned int test_case_4 = fibonacci(3);
+    uint64_t test_case_4 = fibonacci(3);
     assert(test_case_4 == 2);
     std::cout << "Passed Test 4!" << std::endl;
 
-    unsigned int test_case_5 = fibonacci(4);
+    uint64_t test_case_5 = fibonacci(4);
     assert(test_case_5 == 3);
     std::cout << "Passed Test 5!" << std::endl;
 
-    unsigned int test_case_6 = fibonacci(15);
+    uint64_t test_case_6 = fibonacci(15);
     assert(test_case_6 == 610);
     std::cout << "Passed Test 6!" << std::endl << std::endl;
 }

range_queries/heavy_light_decomposition.cpp

@@ -0,0 +1,639 @@
+/**
+ * @file
+ * @brief [Heavy Light
+ * Decomposition](https://en.wikipedia.org/wiki/Heavy_path_decomposition)
+ * implementation
+ * @author [Aniruthan R](https://github.com/aneee004)
+ *
+ * @details
+ * Heavy-Light Decomposition is a technique on trees, that supports the
+ * following:
+ * 1. Update node s, with a value v
+ * 2. Return the (sum) of all node values on the simple path from a to b
+ * (sum) can also be replced with XOR, OR, AND, min, or max
+ *
+ * The update is done in O(log n) time, and
+ * the query is done in O(log^2 n) time with HLD
+ * where, n is the number of nodes
+ *
+ * The template type is the data type of the value stored in the nodes.
+ * If a non-primitive data-type is used as a template,
+ * the coressponding operators must be overloaded.
+ *
+ * An HLD object can only be created with a constant number of nodes, and
+ * it cannot be changed later. Creaty an empty instance is not supported.
+ *
+ * To start answering updates and queries,
+ * 1. Create an instance of HLD<X> object (obj), with the required data type.
+ * 2. Read in the edge/parent information and update it with obj.add_edge().
+ * Note: The edges addes must be 0 indexed.
+ * 3. Create a vector with initial node values, and call set_node_val() with it.
+ * 4. Call obj.init() to populate the required information for supporting
+ * operations.
+ * 5. Call obj.update(node, new_val), to update the value at index 'node' to the
+ * new value. Note: node must be 0 indexed
+ * 6. Call obj.query(a, b) to get the (sum) of node values in the simple path
+ * from a to b. Note: a and b, must be 0 indexed.
+ *
+ * Sample I/O at the bottom.
+ * @todo Support edge weight queries, by storing the edge weight value in it's
+ * child algorithm verified by testing in CSES path queries:
+ * https://cses.fi/problemset/task/1138
+ */
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstring>
+#include <iostream>
+#include <list>
+#include <numeric>
+#include <string>
+#include <vector>
+
+/**
+ * @namespace range_queries
+ * @brief Algorithms and Data Structures that support range queries and updates.
+ */
+namespace range_queries {
+/**
+ * @namespace heavy_light_decomposition
+ * @brief Heavy light decomposition algorithm
+ */
+namespace heavy_light_decomposition {
+/**
+ * @brief A Basic Tree, which supports binary lifting
+ * @tparam the data type of the values stored in the tree nodes
+ * @details Deleting the default constructor
+ * An instance can only be created with the number of nodes
+ * Defaults:
+ * t_node indexing are zero based
+ * t_root is 0
+ * depth of root_node is 0
+ * Supports:
+ * lift :- lift a node k units up the tree
+ * kth_ancestor :- returns the kth ancestor
+ * lca :- returns the least common ancestor
+ */
+template <typename X> class Tree {
+  //
+
+private:
+  std::vector<std::list<int>>
+      t_adj;         ///< an adjacency list to stores the tree edges
+  const int t_nodes, ///< number of nodes
+      t_maxlift;     ///< maximum possible height of the tree
+  std::vector<std::vector<int>>
+      t_par;                ///< a matrix to store every node's 2^kth parent
+  std::vector<int> t_depth, ///< a vector to store the depth of a node,
+      t_size; ///< a vector to store the subtree size rooted at node
+
+  int t_root;           ///< the root of the tree
+  std::vector<X> t_val; ///< values of nodes
+  template <typename T> friend class HLD;
+
+  /**
+   * @brief Utility function to compute sub-tree sizes
+   * @param u current dfs node
+   * @param p the parent of node @param u
+   * @returns void
+   */
+  void dfs_size(int u, int p = -1) {
+    for (const int &v : t_adj[u]) {
+      if (v ^ p) {
+        dfs_size(v, u);
+        t_size[u] += t_size[v];
+      }
+    }
+  }
+
+  /**
+   * @brief Utility function to populate the t_par vector
+   * @param u current dfs node
+   * @param p the parent of node u
+   * @returns void
+   */
+  void dfs_lca(int u, int p = -1) {
+    t_par[u][0] = p;
+    if (p != -1) {
+      t_depth[u] = 1 + t_depth[p];
+    }
+    for (int k = 1; k < t_maxlift; k++) {
+      if (t_par[u][k - 1] != -1) {
+        t_par[u][k] = t_par[t_par[u][k - 1]][k - 1];
+      }
+    }
+
+    for (const int &v : t_adj[u]) {
+      if (v ^ p) {
+        dfs_lca(v, u);
+      }
+    }
+  }
+
+public:
+  /**
+   * @brief Class parameterized constructor, resizes the and initializes the
+   * data members
+   * @param nodes the total number of nodes in the tree
+   */
+  explicit Tree(int nodes)
+      : t_nodes(nodes), t_maxlift(static_cast<int>(floor(log2(nodes))) + 1) {
+    /* Initialize and resize all the vectors */
+    t_root = 0; /* Default */
+    t_adj.resize(t_nodes);
+    t_par.assign(t_nodes, std::vector<int>(t_maxlift, -1));
+    t_depth.assign(t_nodes, 0);
+    t_size.assign(t_nodes, 1);
+    t_val.resize(t_nodes);
+  }
+
+  /**
+   * @brief Adds an undirected edge from node u to node v in the tree
+   * @param u the node where the edge is from
+   * @param v the node where the edge is to
+   * @returns void
+   */
+  void add_edge(const int u, const int v) {
+    t_adj[u].push_back(v);
+    t_adj[v].push_back(u);
+  }
+
+  /**
+   * @brief Set the root for the tree
+   * @param new_root the new root
+   * @returns void
+   */
+  void change_root(int new_root) { t_root = new_root; }
+
+  /**
+   * @brief Set the values for all the nodes
+   * @param node_val a vector of size n, with all the node values where, n is
+   * the number of nodes
+   * @returns void
+   */
+  void set_node_val(const std::vector<X> &node_val) {
+    assert(static_cast<int>(node_val.size()) == t_nodes);
+    t_val = node_val;
+  }
+
+  /**
+   * @brief This function must be called after the tree adjacency list and node
+   * values are populated The function initializes the required parameters, and
+   * populates the segment tree
+   * @returns void
+   */
+  void init() {
+    assert(t_nodes > 0);
+    dfs_size(t_root);
+    dfs_lca(t_root);
+  }
+
+  /**
+   * @brief The function lifts a node, k units up the tree.
+   * The lifting is done in place, and the result is stored in the address
+   * pointed by p.
+   * @param p a pointer to the variable that stores the node id
+   * @param dist the distance to move up the tree
+   * @returns void
+   */
+  void lift(int *const p, int dist) {
+    for (int k = 0; k < t_maxlift; k++) {
+      if (*p == -1) {
+        return;
+      }
+      if (dist & 1) {
+        *p = t_par[*p][k];
+      }
+      dist >>= 1;
+    }
+  }
+
+  /**
+   * @brief The function returns the kth ancestor of a node
+   * @param p the node id whose kth ancestor is to be found
+   * @param dist the distance to move up the tree
+   * @returns the kth ancestor of node
+   */
+  int kth_ancestor(int p, const int &dist) {
+    lift(&p, dist);
+    return p;
+  }
+
+  /**
+   * @brief The function returns the least common ancestor of two nodes
+   * @param a node id_1
+   * @param b node id_2
+   * @returns the least common ancestor of node a, and node b
+   */
+  int lca(int a, int b) {
+    assert(a >= 0 and b >= 0 and a < t_nodes and b < t_nodes);
+    if (t_depth[a] > t_depth[b]) {
+      lift(&a, t_depth[a] - t_depth[b]);
+    }
+    if (t_depth[b] > t_depth[a]) {
+      lift(&b, t_depth[b] - t_depth[a]);
+    }
+    if (a == b) {
+      return a;
+    }
+    for (int k = t_maxlift - 1; k >= 0; k--) {
+      if (t_par[a][k] != t_par[b][k]) {
+        a = t_par[a][k];
+        b = t_par[b][k];
+      }
+    }
+    return t_par[a][0];
+  }
+};
+
+/**
+ * @brief Segment Tree, to store heavy chains
+ * @tparam the data type of the values stored in the tree nodes
+ */
+template <typename X> class SG {
+private:
+  /**
+   * @brief Everything here is private,
+   * and can only be accessed through the methods,
+   * in the derived class (HLD)
+   */
+
+  std::vector<X> s_tree; ///< the segment tree, stored as a vector
+  int s_size;            ///< number of leaves in the segment tree
+  X sret_init = 0;       ///< inital query return value
+  template <typename T> friend class HLD;
+
+  /**
+   * @brief Function that specifies the type of operation involved when segments
+   * are combined
+   * @param lhs the left segment
+   * @param rhs the right segment
+   * @returns the combined result
+   */
+  X combine(X lhs, X rhs) { return lhs + rhs; }
+
+  /**
+   * @brief Class parameterized constructor. Resizes the and initilizes the data
+   * members.
+   * @param nodes the total number of nodes in the tree
+   * @returns void
+   */
+  explicit SG(int size) {
+    s_size = size;
+    s_tree.assign(2 * s_size, 0ll);
+  }
+
+  /**
+   * @brief Update the value at a node
+   * @param p the node to be udpated
+   * @param v the update value
+   * @returns void
+   */
+  void update(int p, X v) {
+    for (p += s_size; p > 0; p >>= 1) {
+      s_tree[p] += v;
+    }
+  }
+
+  /**
+   * @brief Make a range query from node label l to node label r
+   * @param l node label where the path starts
+   * @param r node label where the path ends
+   * @returns void
+   */
+  X query(int l, int r) {
+    X lhs = sret_init, rhs = sret_init;
+    for (l += s_size, r += s_size + 1; l < r; l >>= 1, r >>= 1) {
+      if (l & 1) {
+        lhs = combine(lhs, s_tree[l++]);
+      }
+      if (r & 1) {
+        rhs = combine(s_tree[--r], rhs);
+      }
+    }
+    return combine(lhs, rhs);
+  }
+
+  /**
+   * @brief Set the initialization for the query data type, based on requirement
+   *
+   * @details
+   * Change the sret_init, based on requirement:
+   * * Sum Query: 0 (Default)
+   * * XOR Query: 0 (Default)
+   * * Min Query: Infinity
+   * * Max Query: -Infinity
+   * @param new_sret_init the new init
+   */
+  void set_sret_init(X new_sret_init) { sret_init = new_sret_init; }
+};
+
+/**
+ * @brief The Heavy-Light Decomposition class
+ * @tparam the data type of the values stored in the tree nodes
+ */
+template <typename X> class HLD : public Tree<X>, public SG<X> {
+private:
+  int label;                ///< utility member to assign labels in dfs_labels()
+  std::vector<int> h_label, ///< stores the label of a node
+      h_heavychlid,         ///< stores the heavy child of a node
+      h_parent;             ///< stores the top of the heavy chain from a node
+
+  /**
+   * @brief Utility function to assign heavy child to each node (-1 for a leaf
+   * node)
+   * @param u current dfs node
+   * @param p the parent of node u
+   * @returns void
+   */
+  void dfs_hc(int u, int p = -1) {
+    int hc_size = -1, hc_id = -1;
+    for (const int &v : Tree<X>::t_adj[u]) {
+      if (v ^ p) {
+        dfs_hc(v, u);
+        if (Tree<X>::t_size[v] > hc_size) {
+          hc_size = Tree<X>::t_size[v];
+          hc_id = v;
+        }
+      }
+    }
+    h_heavychlid[u] = hc_id;
+  }
+
+  /**
+   * @brief Utility function to assign highest parent that can be reached though
+   * heavy chains
+   * @param u current dfs node
+   * @param p the parent of node u
+   * @returns void
+   */
+  void dfs_par(int u, int p = -1) {
+    if (h_heavychlid[u] != -1) {
+      h_parent[h_heavychlid[u]] = h_parent[u];
+      dfs_par(h_heavychlid[u], u);
+    }
+    for (const int &v : Tree<X>::t_adj[u]) {
+      if (v ^ p and v ^ h_heavychlid[u]) {
+        dfs_par(v, u);
+      }
+    }
+  }
+
+  /**
+   * @brief Utility function to lable the nodes so that heavy chains have a
+   * contigous lable
+   * @param u current dfs node
+   * @param p the parent of node u
+   * @returns void
+   */
+  void dfs_labels(int u, int p = -1) {
+    h_label[u] = label++;
+    if (h_heavychlid[u] != -1) {
+      dfs_labels(h_heavychlid[u], u);
+    }
+    for (const int &v : Tree<X>::t_adj[u]) {
+      if (v ^ p and v ^ h_heavychlid[u]) {
+        dfs_labels(v, u);
+      }
+    }
+  }
+
+  /**
+   * @brief Utility function to break down a path query into two chain queries
+   * @param a node where the path starts
+   * @param b node where the path ends
+   * a and b must belong to a single root to leaf chain
+   * @returns the sum of ndoe values in the simple path from a to b
+   */
+  X chain_query(int a, int b) {
+    X ret = SG<X>::sret_init;
+    if (Tree<X>::t_depth[a] < Tree<X>::t_depth[b]) {
+      std::swap(a, b);
+    }
+    while (Tree<X>::t_depth[a] >= Tree<X>::t_depth[b]) {
+      int l = h_label[h_parent[a]];
+      int r = h_label[a];
+      if (Tree<X>::t_depth[h_parent[a]] < Tree<X>::t_depth[b]) {
+        l += Tree<X>::t_depth[b] - Tree<X>::t_depth[h_parent[a]];
+      }
+      ret = SG<X>::combine(ret, SG<X>::query(l, r));
+      a = Tree<X>::t_par[h_parent[a]][0];
+      if (a == -1) {
+        break;
+      }
+    }
+    return ret;
+  }
+
+public:
+  /**
+   * @brief Class parameterized constructor. Resizes the and initilizes the data
+   * members.
+   * @param nodes the total number of nodes in the tree
+   */
+  explicit HLD<X>(int nodes) : Tree<X>(nodes), SG<X>(nodes) {
+    /* Initialization and resize vectors */
+    label = 0;
+    h_label.assign(Tree<X>::t_nodes, -1);
+    h_heavychlid.assign(Tree<X>::t_nodes, -1);
+    h_parent.resize(Tree<X>::t_nodes);
+    iota(h_parent.begin(), h_parent.end(), 0);
+  }
+
+  /**
+   * @brief This function must be called after the tree adjacency list and node
+   * values are populated The function initializes the required parametes, and
+   * populates the segment tree
+   * @returns void
+   */
+  void init() {
+    Tree<X>::init();
+
+    // Fill the heavy child, greatest parent, and labels
+    label = 0;
+    dfs_hc(Tree<X>::t_root);
+    dfs_par(Tree<X>::t_root);
+    dfs_labels(Tree<X>::t_root);
+
+    // Segment Tree Initialization
+    for (int i = 0; i < Tree<X>::t_nodes; i++) {
+      SG<X>::s_tree[h_label[i] + Tree<X>::t_nodes] = Tree<X>::t_val[i];
+    }
+    for (int i = Tree<X>::t_nodes - 1; i > 0; i--) {
+      SG<X>::s_tree[i] =
+          SG<X>::combine(SG<X>::s_tree[i << 1], SG<X>::s_tree[i << 1 | 1]);
+    }
+  }
+
+  /**
+   * @brief This function updates the value at node with val
+   * @param node the node where the update is done
+   * @param val the value that is being updated
+   * @returns void
+   */
+  void update(int node, X val) {
+    X diff = val - Tree<X>::t_val[node];
+    SG<X>::update(h_label[node], diff);
+    Tree<X>::t_val[node] = val;
+  }
+
+  /**
+   * @brief This function returns the sum of node values in the simple path from
+   * from node_1 to node_2
+   * @param a the node where the simple path starts
+   * @param b the node where the simple path ends
+   * (parameters are interchangeable, i.e., the function is commutative)
+   * @returns the sum of node values in the simple path from a to b
+   */
+  X query(int a, int b) {
+    int lc = Tree<X>::lca(a, b);
+    X ret = SG<X>::sret_init;
+    assert(lc != -1);
+    ret += chain_query(a, lc);
+    ret += chain_query(b, lc);
+    return ret - Tree<X>::t_val[lc];
+  }
+};
+} // namespace heavy_light_decomposition
+} // namespace range_queries
+
+/**
+ * Test implementations
+ * @returns none
+ */
+static void test_1() {
+  std::cout << "Test 1:\n";
+
+  // Test details
+  int n = 5;
+  std::vector<int64_t> node_values = {4, 2, 5, 2, 1};
+  std::vector<std::vector<int>> edges = {{1, 2}, {1, 3}, {3, 4}, {3, 5}};
+  std::vector<std::vector<int>> queries = {
+      {2, 1, 4},
+      {1, 3, 2},
+      {2, 1, 4},
+  };
+  std::vector<int> expected_result = {11, 8};
+  std::vector<int> code_result;
+
+  range_queries::heavy_light_decomposition::HLD<int64_t> hld(n);
+  hld.set_node_val(node_values);
+  for (int i = 0; i < n - 1; i++) {
+    int u = edges[i][0], v = edges[i][1];
+    hld.add_edge(u - 1, v - 1);
+  }
+  hld.init();
+  for (const auto &q : queries) {
+    int type = q[0];
+    if (type == 1) {
+      int p = q[1], x = q[2];
+      hld.update(p - 1, x);
+    } else if (type == 2) {
+      int a = q[1], b = q[2];
+      code_result.push_back(hld.query(a - 1, b - 1));
+    } else {
+      continue;
+    }
+  }
+  for (int i = 0; i < static_cast<int>(expected_result.size()); i++) {
+    assert(expected_result[i] == code_result[i]);
+  }
+  std::cout << "\nTest 1 passed!\n";
+}
+
+/**
+ * Second test implementations
+ * @returns void
+ */
+static void test_2() {
+  std::cout << "Test 2:\n";
+
+  // Test details (Bamboo)
+  int n = 10;
+  std::vector<int64_t> node_values = {1, 8, 6, 8, 6, 2, 9, 2, 3, 2};
+  std::vector<std::vector<int>> edges = {
+      {10, 5}, {6, 2}, {10, 7}, {5, 2}, {3, 9}, {8, 3}, {1, 4}, {6, 4}, {8, 7}};
+  std::vector<std::vector<int>> queries = {
+      {2, 1, 10}, {2, 1, 6}, {1, 3, 4}, {2, 1, 9}, {1, 5, 3},
+      {1, 7, 8},  {2, 1, 4}, {2, 1, 8}, {1, 1, 4}, {1, 2, 7}};
+  std::vector<int> expected_result = {27, 11, 45, 9, 34};
+  std::vector<int> code_result;
+
+  range_queries::heavy_light_decomposition::HLD<int64_t> hld(n);
+  hld.set_node_val(node_values);
+  for (int i = 0; i < n - 1; i++) {
+    int u = edges[i][0], v = edges[i][1];
+    hld.add_edge(u - 1, v - 1);
+  }
+  hld.init();
+  for (const auto &q : queries) {
+    int type = q[0];
+    if (type == 1) {
+      int p = q[1], x = q[2];
+      hld.update(p - 1, x);
+    } else if (type == 2) {
+      int a = q[1], b = q[2];
+      code_result.push_back(hld.query(a - 1, b - 1));
+    } else {
+      continue;
+    }
+  }
+  for (int i = 0; i < static_cast<int>(expected_result.size()); i++) {
+    assert(expected_result[i] == code_result[i]);
+  }
+  std::cout << "\nTest2 passed!\n";
+}
+
+/**
+ * Third test implementations
+ * @returns void
+ */
+static void test_3() {
+  std::cout << "Test 3:\n";
+
+  // Test details
+  int n = 8;
+  std::vector<int64_t> node_values = {1, 8, 6, 8, 6, 2, 9, 2};
+  std::vector<std::vector<int>> edges = {{1, 2}, {2, 3}, {3, 4}, {1, 5},
+                                         {6, 3}, {7, 5}, {8, 7}};
+  std::vector<std::vector<int>> queries = {
+      {2, 6, 8}, {2, 3, 6}, {1, 3, 4}, {2, 7, 1}, {1, 5, 3},
+      {1, 7, 8}, {2, 6, 4}, {2, 7, 8}, {1, 1, 4}, {1, 2, 7}};
+  std::vector<int> expected_result = {34, 8, 16, 14, 10};
+  std::vector<int> code_result;
+
+  range_queries::heavy_light_decomposition::HLD<int64_t> hld(n);
+  hld.set_node_val(node_values);
+  for (int i = 0; i < n - 1; i++) {
+    int u = edges[i][0], v = edges[i][1];
+    hld.add_edge(u - 1, v - 1);
+  }
+  hld.init();
+  for (const auto &q : queries) {
+    int type = q[0];
+    if (type == 1) {
+      int p = q[1], x = q[2];
+      hld.update(p - 1, x);
+    } else if (type == 2) {
+      int a = q[1], b = q[2];
+      code_result.push_back(hld.query(a - 1, b - 1));
+    } else {
+      continue;
+    }
+  }
+  for (int i = 0; i < static_cast<int>(expected_result.size()); i++) {
+    assert(expected_result[i] == code_result[i]);
+  }
+  std::cout << "\nTest3 passed!\n";
+}
+
+/**
+ * Main function
+ */
+int main() {
+  test_1();
+  test_2();
+  test_3();
+  return 0;
+}

search/median_search.cpp

@@ -1,78 +1,156 @@
 /**
- * \file
- * \brief [Median search](https://en.wikipedia.org/wiki/Median_search) algorithm
- * \warning This core is erroneous and gives invorrect answers. Tested using
- * cases from [here](https://brilliant.org/wiki/median-finding-algorithm/)
- * \ingroup median search
- * \{
+ * @file median_search.cpp
+ * @brief Implementation of [Median search](https://en.wikipedia.org/wiki/Median_of_medians) algorithm.
+ * @cases from [here](https://brilliant.org/wiki/median-finding-algorithm/)
+ *
+ * @details
+ * Given an array A[1,...,n] of n numbers and an index i, where 1 ≤ i ≤ n, find the i-th smallest element of A.
+ * median_of_medians(A, i):
+ *  #divide A into sublists of len 5
+ *  sublists = [A[j:j+5] for j in range(0, len(A), 5)]
+ *  medians = [sorted(sublist)[len(sublist)/2] for sublist in sublists]
+ *  if len(medians) <= 5:
+ *	pivot = sorted(medians)[len(medians)/2]
+ *  else:
+ *      #the pivot is the median of the medians
+ *      pivot = median_of_medians(medians, len(medians)/2)
+ *  #partitioning step
+ *  low = [j for j in A if j < pivot]
+ *  high = [j for j in A if j > pivot]
+ *  k = len(low)
+ *   if i < k:
+ *      return median_of_medians(low,i)
+ *   elif i > k:
+ *      return median_of_medians(high,i-k-1)
+ *  else: #pivot = k
+ *       return pivot
+ *
+ * \note this algorithm implements median search for only arrays which have distinct elements
+ *
+ * Here are some example lists you can use to see how the algorithm works
+ * A = [1,2,3,4,5,1000,8,9,99] (Contain Unique Elements)
+ * B = [1,2,3,4,5,6] (Contains Unique Elements)
+ * print median_of_medians(A, 0) #should be 1
+ * print median_of_medians(A,7) #should be 99
+ * print median_of_medians(B,4) #should be 5
+ *
+ * @author Unknown author
+ * @author [Sushil Kumar](https://github.com/Rp-sushil)
  */
-#include <algorithm>
+
 #include <iostream>
+#include <algorithm>
 #include <vector>
+#include <cassert>
 
 /**
- * @todo add documentation
+ * @namespace search
+ * @brief Search algorithms
  */
-template <class X>
-void comp(X x, std::vector<int> *s1, std::vector<int> *s2,
-          std::vector<int> *s3) {
-    if (s1->size() >= x && s1->size() + s2->size() < x) {
-        std::cout << (*s2)[0] << " is the " << x + 1 << "th element from front";
-    } else if (s1->size() > x) {
-        std::sort(s1->begin(), s1->end());
-        std::cout << (*s1)[x] << " is the " << x + 1 << "th element from front";
-    } else if (s1->size() + s2->size() <= x && s3->size() > x) {
-        std::sort(s3->begin(), s3->end());
-        std::cout << (*s3)[x - s1->size() - s2->size()] << " is the " << x + 1
-                  << "th element from front";
-    } else {
-        std::cout << x + 1 << " is invalid location";
-    }
+namespace search {
+/**
+ * @namespace median_search
+ * @brief Functions for [Median search](https://en.wikipedia.org/wiki/Median_search) algorithm
+ */
+namespace median_search {
+/**
+* This function search the element in an array for the given index.
+* @param A array where numbers are saved
+* @param idx current index in array
+* @returns corresponding element which we want to search.
+*/  
+int median_of_medians(const std::vector<int>& A,  const int& idx) {
+	int pivot = 0;					// initialized with zero
+	std::vector<int> a(A.begin(), A.end());
+	std::vector<int> m;
+	int r = a.size();
+	for(int i = 0; i < r; i += 5){
+		std::sort(a.begin() + i, a.begin() + std::min(r, i + 5));
+		int mid = (i + std::min(r, i + 5)) / 2;
+		m.push_back(a[mid]);
+	}
+	int sz = int(m.size());
+	if(sz <= 5){
+		std::sort(m.begin(), m.end());
+		pivot = m[(sz- 1) / 2];
+	}
+	else{
+		pivot = median_of_medians(m, idx);
+	}
+	std::vector<int> low;
+	std::vector<int> high;
+	for(int i = 0; i < r; i++){
+		if(a[i] < pivot){
+			low.push_back(a[i]);
+		}
+		else if(a[i] > pivot){
+			high.push_back(a[i]);
+		}
+	}
+	int k = int(low.size());
+	if(idx < k){
+		return median_of_medians(low, idx);
+	}
+	else if(idx > k){
+		return median_of_medians(high, idx-k-1);
+	}
+	else{
+		return pivot;
+	}
 }
+}  // namespace median_search
+}  // namespace search
 
-#define MAX_NUM 20  ///< maximum number of values to sort from
+/**
+ * Function to test above algorithm
+ */
+void test(){
+	std::vector<int> A{25,21,98,100,76,22,43,60,89,87};
+	int i = 3;
+	assert(A[6] == search::median_search::median_of_medians(A, i));		// A[6]  = 43, is the fourth smallest element.
+	std::cout << "test case:1 passed\n";
+	
+	std::vector<int> B{1,2,3,4,5,6};
+	int j = 4;
+	assert(B[4] == search::median_search::median_of_medians(B, j));		// B[4] = 5, is the fifth smallest element.
+	std::cout << "test case:2 passed\n";
+	
+	std::vector<int> C{1,2,3,4,5,1000,8,9,99};
+	int k = 3;
+	assert(C[3] == search::median_search::median_of_medians(C, k)); 	// C[3] = 4, is the fourth smallest element.
+	std::cout << "test case:3 passed\n";
+	std::cout << "--All tests passed--\n";
+}
 
 /**
  * Main function
  */
-int main() {
-    std::vector<int> v{25, 21, 98, 100, 76, 22, 43, 60, 89, 87};
-    std::vector<int> s1;
-    std::vector<int> s2;
-    std::vector<int> s3;
-
-    // creates an array of random numbers
-    // for (int i = 0; i < MAX_NUM; i++) {
-    //     int r = std::rand() % 1000;
-    //     v.push_back(r);
-    //     std::cout << r << " ";
-    // }
-    for (int r : v) std::cout << r << " ";
-
-    int median = std::rand() % 1000;  // initialize to a random numnber
-
-    std::cout << "\nmedian=" << median << std::endl;
-    int avg1, avg2, avg3, sum1 = 0, sum2 = 0, sum3 = 0;
-
-    for (int i = 0; i < v.size(); i++) {  // iterate through all numbers
-        if (v.back() == v[median]) {
-            avg1 = sum1 + v.back();
-            s2.push_back(v.back());
-        } else if (v.back() < v[median]) {
-            avg2 = sum2 + v.back();
-            s1.push_back(v.back());
-        } else {
-            avg3 = sum3 + v.back();
-            s3.push_back(v.back());
-        }
-        v.pop_back();
-    }
-
-    int x;
-    std::cout << "enter the no. to be searched form begining:- ";
-    std::cin >> x;
-    comp(x - 1, &s1, &s2, &s3);
-
-    return 0;
+int main()
+{
+	test();
+	int n = 0;
+	std::cout << "Enter Size of Array: ";
+	std::cin >> n;
+	std::vector<int> a(n);
+	std::cout << "Enter Array: ";
+	for(int i = 0; i < n; i++){
+		std::cin >> a[i];
+	}
+	std::cout << "Median: ";			// Median defination: https://en.wikipedia.org/wiki/Median
+	int x = search::median_search::median_of_medians(a,  (n - 1) / 2);
+	if(n % 2 == 0){
+		int y = search::median_search::median_of_medians(a, n / 2);
+		std::cout << (float(x) + float(y))/2.0;
+	}
+	else{
+		std::cout << x;
+	}
+	std::cout << "\nTo find i-th smallest element ";
+       	std::cout << "\nEnter i: ";
+	int idx = 0;
+	std::cin >> idx;
+	idx--;
+	std::cout << idx + 1<< "-th smallest element: " << search::median_search::median_of_medians(a, idx) << '\n';
+	return 0;
 }
-/// }
+

sorting/pigeonhole_sort.cpp

@@ -0,0 +1,133 @@
+/**
+ * @file
+ * @brief Implementation of [Pigeonhole Sort algorithm]
+ * (https://en.wikipedia.org/wiki/Pigeonhole_sort)
+ * @author [Lownish](https://github.com/Lownish)
+ * @details
+ * Pigeonhole sorting is a sorting algorithm that is suitable for sorting lists
+ * of elements where the number of elements and the number of possible key
+ * values are approximately the same. It requires O(n + Range) time where n is
+ * number of elements in input array and ‘Range’ is number of possible values in
+ * array.
+ *
+ * The time Complexity of the algorithm is \f$O(n+N)\f$.
+ */
+
+#include <algorithm>  //for std::is_sorted
+#include <array>      //for std::array
+#include <cassert>    //for assert
+#include <iostream>   //for io operations
+
+/**
+ * @namespace sorting
+ * @brief Sorting algorithms
+ */
+namespace sorting {
+
+/**
+ * Pigeonhole sorting of array of size n
+ * The function will sort the array through Pigeonhole algorithm and print
+ * @param arr unsorted array of elements
+ * @returns sorted array of elements
+ */
+template <std::size_t N>
+std::array<int, N> pigeonSort(std::array<int, N> arr) {
+    // Finding min and max*
+    auto min = std::min_element(std::begin(arr), std::end(arr));
+    auto max = std::max_element(std::begin(arr), std::end(arr));
+
+    // Range refers to the number of holes required
+    int range = *max - *min + 1;
+    int *hole = new int[range]();
+
+    // Copying all array values to pigeonhole
+    for (int i = 0; i < N; i++) {
+        hole[arr[i] - *min] = arr[i];
+    }
+
+    // Deleting elements from list and storing to original array
+    int count = 0;
+    for (int i = 0; i < range; i++) {
+        while (hole[i] != '\0') {
+            arr[count] = hole[i];
+            hole[i] = {};
+            count++;
+        }
+    }
+    delete[] hole;
+
+    return arr;
+}
+}  // namespace sorting
+
+/**
+ * Test function 1 with unsorted array
+ * {8, 3, 2, 7, 4, 6, 8}
+ * @returns none
+ */
+static void test_1() {
+    const int n = 7;
+    std::array<int, n> test_array = {8, 3, 2, 7, 4, 6, 8};
+
+    test_array = sorting::pigeonSort<n>(test_array);
+
+    assert(std::is_sorted(std::begin(test_array), std::end(test_array)));
+
+    // Printing sorted array
+    for (int i = 0; i < n; i++) {
+        std::cout << test_array.at(i) << " ";
+    }
+    std::cout << "\nPassed\n";
+}
+
+/**
+ * Test function 2 with unsorted array
+ * {802, 630, 20, 745, 52, 300, 612, 932, 78, 187}
+ * @returns none
+ */
+static void test_2() {
+    const int n = 10;
+    std::array<int, n> test_array = {802, 630, 20,  745, 52,
+                                     300, 612, 932, 78,  187};
+
+    test_array = sorting::pigeonSort<n>(test_array);
+
+    assert(std::is_sorted(std::begin(test_array), std::end(test_array)));
+
+    // Printing sorted array
+    for (int i = 0; i < n; i++) {
+        std::cout << test_array.at(i) << " ";
+    }
+    std::cout << "\nPassed\n";
+}
+
+/**
+ * Test function 1 with unsorted array
+ * {11,13,12,14}
+ * @returns none
+ */
+static void test_3() {
+    const int n = 4;
+    std::array<int, n> test_array = {11, 13, 12, 14};
+
+    test_array = sorting::pigeonSort<n>(test_array);
+
+    assert(std::is_sorted(std::begin(test_array), std::end(test_array)));
+
+    // Printing sorted array
+    for (int i = 0; i < n; i++) {
+        std::cout << test_array.at(i) << " ";
+    }
+    std::cout << "\nPassed\n";
+}
+
+/**
+ * Main function
+ */
+int main() {
+    test_1();
+    test_2();
+    test_3();
+
+    return 0;
+}