#include "matrix.hpp"
#include <stdint.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <cstdint>
#include <vector>
#include <algorithm>

std::vector<std::vector<uint64_t>> random_adjacency(const uint64_t vertex_count) {
    std::vector<std::vector<uint64_t>> adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    srand(time(NULL));

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (column_index != row_index) {
                adjacency_matrix[row_index][column_index] = rand() % 2;
            }
        }
    }

    return adjacency_matrix;
}

std::vector<std::vector<uint64_t>> calculate_distance_matrix(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> distance_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> power_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), power_matrix.begin());

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (row_index == column_index) {
                continue;
            } else if (adjacency_matrix[row_index][column_index] == 1) {
                distance_matrix[row_index][column_index] = 1;
            } else {
                distance_matrix[row_index][column_index] = UINT64_MAX;
            }
        }
    }

    for(uint64_t k = 2; k <= vertex_count; k++) {
        power_matrix = gemm_basic(adjacency_matrix, power_matrix);

        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
                if (power_matrix[row_index][column_index] != 0 && distance_matrix[row_index][column_index] == UINT64_MAX) {
                    distance_matrix[row_index][column_index] = k;
                }
            }
        }
    }

    return distance_matrix;
}

std::vector<uint64_t> get_eccentricities(const std::vector<std::vector<uint64_t>>& distance_matrix) {
    uint64_t vertex_count = distance_matrix.size();
    std::vector<uint64_t> eccentricities(vertex_count, UINT64_MAX);
    uint64_t eccentricity;

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        eccentricity = 0;

        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (distance_matrix[row_index][column_index] > eccentricity) {
                eccentricity = distance_matrix[row_index][column_index];
            }
        }

        if (eccentricity == 0) {
            break;
        }
        eccentricities[row_index] = eccentricity;
    }

    return eccentricities;
}

uint64_t get_radius(const std::vector<uint64_t>& eccentricities) {
    uint64_t radius = UINT64_MAX;

    for (uint64_t index = 0; index < eccentricities.size(); index++) {
        if (eccentricities[index] < radius) {
            radius = eccentricities[index];
        }
    }

    return radius;
}

uint64_t get_diameter(const std::vector<uint64_t>& eccentricities) {
    uint64_t diamter = 0;

    for (uint64_t index = 0; index < eccentricities.size(); index++) {
        if (eccentricities[index] > diamter) {
            diamter = eccentricities[index];
        }
    }

    return diamter;
}

std::vector<uint64_t> get_centre(const std::vector<uint64_t>& eccentricities) {
    std::vector<uint64_t> centre;
    uint64_t radius = get_radius(eccentricities);

    for (uint64_t index = 0; index < eccentricities.size(); index++) {
        if (eccentricities[index] == radius) {
            centre.push_back(index + 1);
        }
    }

    centre.shrink_to_fit();
    return centre;
}

std::vector<std::vector<uint64_t>> calculate_path_matrix(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> power_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), power_matrix.begin());

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (row_index == column_index || adjacency_matrix[row_index][column_index] == 1) {
                path_matrix[row_index][column_index] = 1;
            }
        }
    }

    for(uint64_t k = 2; k <= vertex_count; k++) {
        power_matrix = gemm_basic(adjacency_matrix, power_matrix);

        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
                if (power_matrix[row_index][column_index] != 0) {
                    path_matrix[row_index][column_index] = 1;
                }
            }
        }
    }

    return path_matrix;
}

std::vector<std::vector<uint64_t>> find_components_basic(const std::vector<std::vector<uint64_t>>& path_matrix) {
    uint64_t vertex_count = path_matrix.size();
    std::vector<std::vector<uint64_t>> components;
    std::vector<uint64_t> component;
    components.reserve(0);
    component.reserve(vertex_count);

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        component.clear();

        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (path_matrix[row_index][column_index] == 1) {
                component.push_back(column_index + 1);
            }
        }

        component.shrink_to_fit();
        std::sort(component.begin(), component.end());

        auto iterator = find(components.begin(), components.end(), component);

        if (iterator == components.end()) {
            components.push_back(component);
        }
    }

    components.shrink_to_fit();
    return components;
}

void dfs(const std::vector<std::vector<uint64_t>>& adjacency_matrix, const uint64_t vertex, std::vector<bool>& visited) {
    visited[vertex] = true;

    for (uint64_t neighbor_vertex = 0; neighbor_vertex < adjacency_matrix.size(); neighbor_vertex++) {
        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
            continue;
        }
        if (visited[neighbor_vertex]) {
            continue;
        }
        dfs(adjacency_matrix, neighbor_vertex, visited);
    }
}

std::vector<std::vector<uint64_t>> find_components_dfs(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> components;
    std::vector<uint64_t> component;
    std::vector<bool> visited(vertex_count, false);
    components.reserve(0);
    component.reserve(vertex_count);

    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
        component.clear();

        dfs(adjacency_matrix, vertex, visited);

        for (uint64_t index = 0; index < vertex_count; index++) {
            if (visited[index]) {
                component.push_back(index + 1);
            }
        }

        component.shrink_to_fit();
        std::sort(component.begin(), component.end());

        auto iterator = find(components.begin(), components.end(), component);

        if (iterator == components.end()) {
            components.push_back(component);
        }
    }

    components.shrink_to_fit();
    return components;
}

std::vector<std::vector<uint64_t>> find_bridges_basic(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> components = find_components_basic(path_matrix);
    std::vector<std::vector<uint64_t>> temp_components;
    std::vector<std::vector<uint64_t>> bridges;
    std::vector<uint64_t> bridge(2, 0);

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (row_index == column_index) {
                continue;
            }

            std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());

            temp_adjacency_matrix[row_index][column_index] = 0;
            temp_adjacency_matrix[column_index][row_index] = 0;

            path_matrix = calculate_path_matrix(temp_adjacency_matrix);
            temp_components = find_components_basic(path_matrix);

            if (temp_components.size() <= components.size()) {
                continue;
            }

            if (column_index < row_index) {
                bridge[0] = column_index + 1;
                bridge[1] = row_index + 1;
            } else {
                bridge[0] = row_index + 1;
                bridge[1] = column_index + 1;
            }

            auto iterator = find(bridges.begin(), bridges.end(), bridge);

            if (iterator == bridges.end()) {
                bridges.push_back(bridge);
            }
        }
    }

    bridges.shrink_to_fit();
    return bridges;
}

std::vector<std::vector<uint64_t>> find_bridges_dfs_v1(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> components = find_components_dfs(adjacency_matrix);
    std::vector<std::vector<uint64_t>> temp_components;
    std::vector<std::vector<uint64_t>> bridges;
    std::vector<uint64_t> bridge(2, 0);

    for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
        for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
            if (row_index == column_index) {
                continue;
            }

            std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());

            temp_adjacency_matrix[row_index][column_index] = 0;
            temp_adjacency_matrix[column_index][row_index] = 0;

            temp_components = find_components_dfs(temp_adjacency_matrix);

            if (temp_components.size() <= components.size()) {
                continue;
            }

            if (column_index < row_index) {
                bridge[0] = column_index + 1;
                bridge[1] = row_index + 1;
            } else {
                bridge[0] = row_index + 1;
                bridge[1] = column_index + 1;
            }

            auto iterator = find(bridges.begin(), bridges.end(), bridge);

            if (iterator == bridges.end()) {
                bridges.push_back(bridge);
            }
        }
    }

    bridges.shrink_to_fit();
    return bridges;
}

void dfs_bridges(const std::vector<std::vector<uint64_t>>& adjacency_matrix, const uint64_t vertex, const uint64_t parent_vertex, std::vector<bool>& visited,
                 uint64_t current_time, std::vector<uint64_t>& discovery_time, std::vector<uint64_t>& lowest_time, std::vector<std::vector<uint64_t>>& bridges) {
    current_time++;
    visited[vertex] = true;
    discovery_time[vertex] = current_time;
    lowest_time[vertex] = current_time;
    std::vector<uint64_t> bridge(2, 0);

    for (uint64_t neighbor_vertex = 0; neighbor_vertex < adjacency_matrix.size(); neighbor_vertex++) {
        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
            continue;
        }
        if (parent_vertex != neighbor_vertex && visited[neighbor_vertex]) {
            if (lowest_time[vertex] > discovery_time[neighbor_vertex]) {
                lowest_time[vertex] = discovery_time[neighbor_vertex];
            }
            continue;
        }
        if (visited[neighbor_vertex]) {
            continue;
        }

        dfs_bridges(adjacency_matrix, neighbor_vertex, vertex, visited, current_time, discovery_time, lowest_time, bridges);

        if (lowest_time[vertex] > lowest_time[neighbor_vertex]) {
            lowest_time[vertex] = lowest_time[neighbor_vertex];
        }
        if (discovery_time[vertex] >= lowest_time[neighbor_vertex]) {
            continue;
        }
        bridge[0] = vertex + 1;
        bridge[1] = neighbor_vertex + 1;

        auto iterator = find(bridges.begin(), bridges.end(), bridge);

        if (iterator == bridges.end()) {
            bridges.push_back(bridge);
        }
    }
}

std::vector<std::vector<uint64_t>> find_bridges_dfs_v2(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<bool> visited(vertex_count, false);
    uint64_t current_time = 0;
    std::vector<uint64_t> discovery_time(vertex_count, 0);
    std::vector<uint64_t> lowest_time(vertex_count, 0);
    std::vector<std::vector<uint64_t>> bridges;

    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
        if (!visited[vertex]) {
            dfs_bridges(adjacency_matrix, vertex, UINT64_MAX, visited, current_time, discovery_time, lowest_time, bridges);
        }
    }

    bridges.shrink_to_fit();
    return bridges;
}

std::vector<uint64_t> find_articulations_basic(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> components = find_components_basic(path_matrix);
    std::vector<std::vector<uint64_t>> temp_components;
    std::vector<uint64_t> articulations;

    for (uint64_t i = 0; i < vertex_count; i++) {
        std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());

        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
                temp_adjacency_matrix[row_index][i] = 0;
                temp_adjacency_matrix[i][column_index] = 0;
            }
        }

        path_matrix = calculate_path_matrix(temp_adjacency_matrix);
        temp_components = find_components_basic(path_matrix);

        // the + 1 is needed because I am not removing the vertex, I am just removing all of its edges
        // removing all of its edges, means it itself becomes a component, which needs to be accounted for
        if (temp_components.size() > components.size() + 1) {
            articulations.push_back(i + 1);
        }
    }

    articulations.shrink_to_fit();
    std::sort(articulations.begin(), articulations.end());
    return articulations;
}

std::vector<uint64_t> find_articulations_dfs_v1(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<std::vector<uint64_t>> path_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> temp_adjacency_matrix(vertex_count, std::vector<uint64_t>(vertex_count, 0));
    std::vector<std::vector<uint64_t>> components = find_components_dfs(adjacency_matrix);
    std::vector<std::vector<uint64_t>> temp_components;
    std::vector<uint64_t> articulations;

    for (uint64_t i = 0; i < vertex_count; i++) {
        std::copy(adjacency_matrix.begin(), adjacency_matrix.end(), temp_adjacency_matrix.begin());

        for (uint64_t row_index = 0; row_index < vertex_count; row_index++) {
            for (uint64_t column_index = 0; column_index < vertex_count; column_index++) {
                temp_adjacency_matrix[row_index][i] = 0;
                temp_adjacency_matrix[i][column_index] = 0;
            }
        }

        temp_components = find_components_dfs(temp_adjacency_matrix);

        // the + 1 is needed because I am not removing the vertex, I am just removing all of its edges
        // removing all of its edges, means it itself becomes a component, which needs to be accounted for
        if (temp_components.size() > components.size() + 1) {
            articulations.push_back(i + 1);
        }
    }

    articulations.shrink_to_fit();
    std::sort(articulations.begin(), articulations.end());
    return articulations;
}

void dfs_articulations(const std::vector<std::vector<uint64_t>>& adjacency_matrix, const uint64_t vertex, const uint64_t parent_vertex, std::vector<bool>& visited,
                 uint64_t current_time, std::vector<uint64_t>& discovery_time, std::vector<uint64_t>& lowest_time, std::vector<uint64_t>& articulations) {
    current_time++;
    visited[vertex] = true;
    discovery_time[vertex] = current_time;
    lowest_time[vertex] = current_time;
    uint64_t child_count = 0;
    bool is_articulation = false;

    for (uint64_t neighbor_vertex = 0; neighbor_vertex < adjacency_matrix.size(); neighbor_vertex++) {
        if (adjacency_matrix[vertex][neighbor_vertex] != 1) {
            continue;
        }
        if (visited[neighbor_vertex]) {
            if (lowest_time[vertex] > discovery_time[neighbor_vertex]) {
                lowest_time[vertex] = discovery_time[neighbor_vertex];
            }
            continue;
        }

        child_count++;

        dfs_articulations(adjacency_matrix, neighbor_vertex, vertex, visited, current_time, discovery_time, lowest_time, articulations);

        if (lowest_time[vertex] > lowest_time[neighbor_vertex]) {
            lowest_time[vertex] = lowest_time[neighbor_vertex];
        }

        if (parent_vertex != UINT64_MAX && discovery_time[vertex] <= lowest_time[neighbor_vertex]) {
            is_articulation = true;
        }
    }
    if (parent_vertex == UINT64_MAX && child_count > 1) {
        is_articulation = true;
    }
    if (is_articulation) {
        articulations.push_back(vertex + 1);
    }
}

std::vector<uint64_t> find_articulations_dfs_v2(const std::vector<std::vector<uint64_t>>& adjacency_matrix) {
    uint64_t vertex_count = adjacency_matrix.size();
    std::vector<bool> visited(vertex_count, false);
    uint64_t current_time = 0;
    std::vector<uint64_t> discovery_time(vertex_count, 0);
    std::vector<uint64_t> lowest_time(vertex_count, 0);
    std::vector<uint64_t> articulations;

    for (uint64_t vertex = 0; vertex < vertex_count; vertex++) {
        if (!visited[vertex]) {
            dfs_articulations(adjacency_matrix, vertex, UINT64_MAX, visited, current_time, discovery_time, lowest_time, articulations);
        }
    }

    articulations.shrink_to_fit();
    std::sort(articulations.begin(), articulations.end());
    return articulations;
}