#include "graph.h"
#include "matrix.h"
#include <string.h>
#include <time.h>
#include <stdlib.h>

void random_adjacency(const ulong vertex_count, ulong matrix[vertex_count][vertex_count]) {
    srand(time(NULL));

    for (ulong row_index = 0; row_index < vertex_count; row_index++) {
        for (ulong column_index = 0; column_index < vertex_count; column_index++) {
            if (column_index == row_index) {
                matrix[row_index][column_index] = 0;
            } else {
                matrix[row_index][column_index] = rand() % 2;
            }
        }
    }
}

void random_adjacency_struct(const Matrix *matrix) {
    if (matrix->row_length != matrix->column_length) {
        return;
    }

    srand(time(NULL));

    for (ulong row_index = 0; row_index < matrix->row_length; row_index++) {
        for (ulong column_index = 0; column_index < matrix->column_length; column_index++) {
            if (column_index == row_index) {
                matrix->values[row_index][column_index] = 0;
            } else {
                matrix->values[row_index][column_index] = rand() % 2;
            }
        }
    }
}

void calculate_distance_matrix(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count], ulong distance_matrix[vertex_count][vertex_count]) {
    ulong power_matrix[vertex_count][vertex_count];
    ulong temp_power_matrix[vertex_count][vertex_count];
    memcpy(power_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(ulong));

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

    for(ulong k = 2; k <= vertex_count; k++) {
        memcpy(temp_power_matrix, power_matrix, vertex_count * vertex_count * sizeof(ulong));
        gemm_basic(vertex_count, vertex_count, adjacency_matrix, vertex_count, vertex_count, temp_power_matrix, power_matrix);

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

int get_eccentricities(const ulong vertex_count, const ulong distance_matrix[vertex_count][vertex_count], ulong eccentricities[vertex_count]) {
    ulong eccentricity;

    // set all eccentricities to infinity in case this is a disconnected graph
    for (ulong index = 0; index < vertex_count; index++) {
        eccentricities[index] = -1;
    }

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

        for (ulong column_index = 0; column_index < vertex_count; column_index++) {
            if (distance_matrix[row_index][column_index] > eccentricity) {
                eccentricity = distance_matrix[row_index][column_index];
            }
        }
        // in case of a disconnected graph
        if (eccentricity == 0) {
            return 1;
        }
        eccentricities[row_index] = eccentricity;
    }

    return 0;
}

ulong get_radius(const ulong vertex_count, const ulong eccentricities[vertex_count]) {
    ulong radius = -1;

    for (ulong index = 0; index < vertex_count; index++) {
        if (eccentricities[index] < radius) {
            radius = eccentricities[index];
        }
    }

    return radius;
}

ulong get_diameter(const ulong vertex_count, const ulong eccentricities[vertex_count]) {
    ulong diamter = 0;

    for (ulong index = 0; index < vertex_count; index++) {
        if (eccentricities[index] > diamter) {
            diamter = eccentricities[index];
        }
    }

    return diamter;
}

void get_centre(const ulong vertex_count, const ulong eccentricities[vertex_count], const ulong radius, ulong centre[vertex_count]) {
    memset(centre, 0, vertex_count * sizeof(ulong));

    for (ulong index = 0; index < vertex_count; index++) {
        if (eccentricities[index] == radius) {
            centre[index] = 1;
        }
    }
}

void calculate_path_matrix(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count], ulong path_matrix[vertex_count][vertex_count]) {
    ulong power_matrix[vertex_count][vertex_count];
    ulong temp_power_matrix[vertex_count][vertex_count];
    memcpy(power_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(ulong));

    for (ulong row_index = 0; row_index < vertex_count; row_index++) {
        for (ulong 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;
            } else {
                path_matrix[row_index][column_index] = 0;
            }
        }
    }

    for(ulong k = 2; k <= vertex_count; k++) {
        memcpy(temp_power_matrix, power_matrix, vertex_count * vertex_count * sizeof(ulong));
        gemm_basic(vertex_count, vertex_count, adjacency_matrix, vertex_count, vertex_count, temp_power_matrix, power_matrix);

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

void find_components_basic(const ulong vertex_count, const ulong path_matrix[vertex_count][vertex_count], ulong components[vertex_count][vertex_count]) {
    ulong component[vertex_count];
    int contains_component;

    memset(components, 0, vertex_count * vertex_count * sizeof(ulong));

    for (ulong row_index = 0; row_index < vertex_count; row_index++) {
        memset(component, 0, vertex_count * sizeof(ulong));
        contains_component = 0;

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

        for (ulong index = 0; index < vertex_count; index++) {
            if (memcmp(components[index], component, vertex_count * sizeof(ulong)) == 0) {
                contains_component = 1;
            }
        }

        if (!contains_component) {
            for (ulong index = 0; index < vertex_count; index++) {
                components[row_index][index] = component[index];
            }
        }
    }
}

void dfs(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count], const ulong vertex, ulong visited[vertex_count]) {
    visited[vertex] = 1;

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

void find_components_dfs(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count], ulong components[vertex_count][vertex_count]) {
    ulong component[vertex_count];
    int contains_component = 0;

    memset(components, 0, vertex_count * vertex_count * sizeof(ulong));

    for (ulong vertex = 0; vertex < vertex_count; vertex++) {
        memset(component, 0, vertex_count * sizeof(ulong));
        contains_component = 0;

        dfs(vertex_count, adjacency_matrix, vertex, component);

        for (ulong index = 0; index < vertex_count; index++) {
            if (memcmp(components[index], component, vertex_count * sizeof(ulong)) == 0) {
                contains_component = 1;
            }
        }

        if (!contains_component) {
            for (ulong index = 0; index < vertex_count; index++) {
                components[vertex][index] = component[index];
            }
        }
    }
}

ulong amount_of_components(const ulong vertex_count, const ulong components[vertex_count][vertex_count]) {
    ulong amount_of_components = 0;

    for (ulong row_index = 0; row_index < vertex_count; row_index++) {
        for (ulong column_index = 0; column_index < vertex_count; column_index++) {
            if (components[row_index][column_index] != 0) {
                amount_of_components++;
                break;
            }
        }
    }

    return amount_of_components;
}

int contains_bridge(const ulong vertex_count, const ulong bridges[vertex_count][2], const ulong bridge[2]) {
    for (ulong index = 0; index < vertex_count; index++) {
        if (memcmp(bridges[index], bridge, 2 * sizeof(ulong)) == 0) {
            return 1;
        }
    }

    return 0;
}

void find_bridges_basic(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count],
                        const ulong components[vertex_count][vertex_count], ulong bridges[vertex_count][2]) {
    ulong path_matrix[vertex_count][vertex_count];
    ulong temp_adjacency_matrix[vertex_count][vertex_count];
    ulong temp_components[vertex_count][vertex_count];
    ulong bridge[2];

    memset(bridges, 0, vertex_count * 2 * sizeof(ulong));

    for (ulong row_index = 0; row_index < vertex_count; row_index++) {
        for (ulong column_index = 0; column_index < vertex_count; column_index++) {
            if (row_index == column_index) {
                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;
            }

            memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(ulong));

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

            calculate_path_matrix(vertex_count, temp_adjacency_matrix, path_matrix);
            find_components_basic(vertex_count, path_matrix, temp_components);

            if (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) && !contains_bridge(vertex_count, bridges, bridge)) {
                bridges[row_index][0] = bridge[0];
                bridges[row_index][1] = bridge[1];
            }
        }
    }
    print_matrix(vertex_count, 2, bridges);
}

void find_bridges_dfs_v1(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count],
                        const ulong components[vertex_count][vertex_count], ulong bridges[vertex_count][2]) {
    ulong temp_adjacency_matrix[vertex_count][vertex_count];
    ulong temp_components[vertex_count][vertex_count];
    ulong bridge[2];

    memset(bridges, 0, vertex_count * 2 * sizeof(ulong));

    for (ulong row_index = 0; row_index < vertex_count; row_index++) {
        for (ulong column_index = 0; column_index < vertex_count; column_index++) {
            if (row_index == column_index) {
                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;
            }

            memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(ulong));

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

            find_components_dfs(vertex_count, temp_adjacency_matrix, temp_components);

            if (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) && !contains_bridge(vertex_count, bridges, bridge)) {
                bridges[row_index][0] = bridge[0];
                bridges[row_index][1] = bridge[1];
            }
        }
    }
}

void dfs_bridges(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count], const ulong vertex, const ulong parent_vertex, ulong visited[vertex_count],
                 ulong current_time, ulong discovery_time[vertex_count], ulong lowest_time[vertex_count], ulong bridges[vertex_count][2]) {
    current_time++;
    visited[vertex] = 1;
    discovery_time[vertex] = current_time;
    lowest_time[vertex] = current_time;

    for (ulong neighbor_vertex = 0; neighbor_vertex < vertex_count; 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(vertex_count, 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]) {
            bridges[neighbor_vertex][0] = vertex + 1;
            bridges[neighbor_vertex][1] = neighbor_vertex + 1;
        }
    }
}

void find_bridges_dfs_v2(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count],
                         const ulong components[vertex_count][vertex_count], ulong bridges[vertex_count][2]) {
    ulong visited[vertex_count];
    ulong current_time = 0;
    ulong discovery_time[vertex_count];
    ulong lowest_time[vertex_count];

    memset(bridges, 0, vertex_count * 2 * sizeof(ulong));
    memset(visited, 0, vertex_count * sizeof(ulong));
    memset(discovery_time, 0, vertex_count * sizeof(ulong));
    memset(lowest_time, 0, vertex_count * sizeof(ulong));

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

void find_articulations_basic(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count],
                              const ulong components[vertex_count][vertex_count], ulong articulations[vertex_count]) {
    ulong path_matrix[vertex_count][vertex_count];
    ulong temp_adjacency_matrix[vertex_count][vertex_count];
    ulong temp_components[vertex_count][vertex_count];

    memset(articulations, 0, vertex_count * sizeof(ulong));

    for (ulong i = 0; i < vertex_count; i++) {
        memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(ulong));

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

        calculate_path_matrix(vertex_count, temp_adjacency_matrix, path_matrix);
        find_components_basic(vertex_count, path_matrix, temp_components);

        // 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 (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) + 1) {
            articulations[i] = i + 1;
        }
    }
}

void find_articulations_dfs_v1(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count],
                              const ulong components[vertex_count][vertex_count], ulong articulations[vertex_count]) {
    ulong temp_adjacency_matrix[vertex_count][vertex_count];
    ulong temp_components[vertex_count][vertex_count];

    memset(articulations, 0, vertex_count * sizeof(ulong));

    for (ulong i = 0; i < vertex_count; i++) {
        memcpy(temp_adjacency_matrix, adjacency_matrix, vertex_count * vertex_count * sizeof(ulong));

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

        find_components_dfs(vertex_count, temp_adjacency_matrix, temp_components);

        // 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 (amount_of_components(vertex_count, temp_components) > amount_of_components(vertex_count, components) + 1) {
            articulations[i] = i + 1;
        }
    }
}

void dfs_articulations(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count], const ulong vertex, const ulong parent_vertex,
                       ulong visited[vertex_count], ulong current_time, ulong discovery_time[vertex_count], ulong lowest_time[vertex_count],
                       ulong articulations[vertex_count]) {
    current_time++;
    visited[vertex] = 1;
    discovery_time[vertex] = current_time;
    lowest_time[vertex] = current_time;
    ulong child_count = 0;
    ulong is_articulation = 0;

    for (ulong neighbor_vertex = 0; neighbor_vertex < vertex_count; 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(vertex_count, 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 != -1 && discovery_time[vertex] <= lowest_time[neighbor_vertex]) {
            is_articulation = 1;
        }
    }
    if (parent_vertex == -1 && child_count > 1) {
        is_articulation = 1;
    }
    if (is_articulation) {
        articulations[vertex] = vertex + 1;
    }
}

void find_articulations_dfs_v2(const ulong vertex_count, const ulong adjacency_matrix[vertex_count][vertex_count],
                         const ulong components[vertex_count][vertex_count], ulong articulations[vertex_count]) {
    ulong visited[vertex_count];
    ulong current_time = 0;
    ulong discovery_time[vertex_count];
    ulong lowest_time[vertex_count];

    memset(articulations, 0, vertex_count * sizeof(ulong));
    memset(visited, 0, vertex_count * sizeof(ulong));
    memset(discovery_time, 0, vertex_count * sizeof(ulong));
    memset(lowest_time, 0, vertex_count * sizeof(ulong));

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