main.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <list>
#include <chrono>
#include <cmath>
#include <queue>
#include <random>
#include <unordered_set>
#include <unordered_map>
#include <algorithm>

using namespace std;
using pii = pair<int,int>;
using Clock = chrono::high_resolution_clock;
using ll = long long;

// Inputs: students, classes, clubs, dorms, school, relationships
const int students = 3000;
const int classCap = 20;
const int avgClasses = 4;
const int avgClubs = 2;
const int dormAmount = 31;
const int NumClubs = 130;
const int NumSchools = 3; // Arts&Sciences, Business, Law
const int SampleSeedsForPaths = 100; // For avg path length estimation
const int repeatRuns = 3; // runs once by default (increase for averaging metrics)


// Outputs: unweighted and weighted connections
 using UnweightedGraph = vector<vector<int>>;
 using WeightedGraph = vector<vector<pair<int, int>>>;

// Relationship Weights
unordered_map<string, int> W = {
    {"Friends", 1},
    {"Classess", 2},
    {"Clubs", 3},
    {"Dorms", 4},
    {"Schools", 5}
};



// initialize random number generator
 std::mt19937_64 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

 int randint(int a, int b) {
    std::uniform_int_distribution<int> dist(a, b);
    return dist(rng);

}

double rand01() {
    std::uniform_real_distribution<double> dist(0.0, 1.0);
    return dist(rng);
}


// Intialize parameters
vector<vector<int>> friends;
vector<vector<int>> classes;
vector<vector<int>> clubs(NumClubs);
vector<vector<int>> dorms(dormAmount );
vector<vector<int>> schools(NumSchools);


// Assign ~4 classes to each student
vector<vector<int>> assignClasses(int students, int classSize, int avgClasses) {
    int totalEnrollments = students * avgClasses;
    int classCountSize = totalEnrollments / classSize;
    int classCrowdFactor = 1.0; // Adjust to increase/decrease number of classes
    int NumClasses = ceil((totalEnrollments / 20) * classCrowdFactor); // Ensure enough classes to accommodate all students

    vector<vector<int>> classes(NumClasses);

    std::random_device rd;
    std::mt19937 rng(rd());
    std::uniform_int_distribution<int> dis(0, NumClasses - 1);

    for (int i = 0; i < students; i++) {
        int enrollments = 0;

        while (enrollments < avgClasses) {
            int c = dis(rng);
            if (classes[c].size() < classSize) {
                classes[c].push_back(i);
                
            }
            enrollments++;
        }

        
    }
    return classes;
};

// Assign dorms evenly
vector<vector<int>> assignDorms( int students, int dormAmount) {
    vector<vector<int>> dorms(dormAmount);

    int baseSize = students / dormAmount;
    int remainder = students % dormAmount;

    int studentPos = 0;

    for (int d = 0; d < dormAmount; d++) {
        int size = (d < remainder) ? baseSize + 1 : baseSize;
    
        for (int r = 0; r < size; r++) {
            dorms[d].push_back(studentPos);
            studentPos++;
        }
    }

    return dorms;
};

// Assign schools
vector<vector<int>> assignSchools(int students, int NumSchools) {
    vector<vector<int>> schools(NumSchools);

    int baseSize = students / NumSchools;
    int remainder = students % NumSchools;

    int studentPos = 0;

    for (int sc = 0; sc < NumSchools; sc++) {
        int size = (sc < remainder) ? baseSize + 1 : baseSize;

        
    
        while (studentPos < students - 1) {
            int c = randint(0, 2);
            schools[c].push_back(studentPos);
            studentPos++;
        };
    }
    
    return schools;
};

// Watts-Strogatz small-world model
// N = number of nodes
// k = each node connects to k/2 nearest neighbors (k must be even)
// beta = rewiring probability (0.0 = lattice, 1.0 = random)

vector<pair<int,int>> WS_model(int N, int k, double beta) {

    vector<unordered_set<int>> adj(N);

    // Create ring latticeWS
    int half = k / 2;
    for (int i = 0; i < N; i++) {
        for (int j = 1; j <= half; j++) {
            int v = (i + j) % N;
            adj[i].insert(v);
            adj[v].insert(i);
        }
    }

    // Rewire edges w/ probability beta
    std::mt19937_64 rng(time(NULL));
    uniform_real_distribution<double> dist(0.0, 1.0);

    for (int i = 0; i < N; i++) {
        vector<int> neighbors(adj[i].begin(), adj[i].end());
        for (int current : neighbors) {
            if (current <= i) continue; // avoid double processing
            if (dist(rng) < beta) {
                // Remove existing edge
                adj[i].erase(current);
                adj[current].erase(i);

                // Pick new neighbor
                int newNeighbor = i;
                int attempts = 0;
                do { 
                    newNeighbor = randint(0, N - 1);
                    attempts++;
                }
                while ((newNeighbor == i || adj[i].count(newNeighbor)) && attempts < 50);
                    
                // Add new edge if valid
                if (newNeighbor != i && adj[i].count(newNeighbor)) {
                    adj[i].insert(newNeighbor);
                    adj[newNeighbor].insert(i);
                } else {
                    // Re-add original edge if no valid new neighbor found
                    adj[i].insert(current);
                    adj[current].insert(i);
                }
            }
        }
    }

    // Convert to edge list
    vector<pair<int,int>> edges;
    for (int v = 0; v < N; v++) {
        for (int e : adj[v]) {
            if (e > v) edges.emplace_back(v, e);
        }
    }
    return edges;

};

void addSchoolEdges_WS(UnweightedGraph &unwg,
    WeightedGraph &wg,
    const vector<vector<int>> &schoolGroups,
    int k,
    double beta,
    int weight ) {
        

    for (const auto &group : schoolGroups) {


        int size = group.size();
        if (size < 2) continue;


        // Create local WS Adjacency for [0..size-1]
        vector<pair<int,int>>edges = WS_model(size, k, beta);

        // Map edges
        for (auto &u : edges) {
            int e = group[u.first];
            int v = group[u.second];


            // Unweighted Graph
            unwg[e].push_back(v);
            unwg[v].push_back(e);


            //Weighted with minimization
            bool foundVE = false;
            for (auto &p : wg[e]) {
                if (p.first == v) {
                    p.second = min(p.second, weight);
                    foundVE = true;
                    break;
                }
            }
            if (!foundVE) {
                wg[e].push_back({v, weight});
                wg[v].push_back({e, weight});
            }
        }

    }
    };

// Assign clubs 0-4 per student
vector<vector<int>> assignClubs(int students, int NumClubs, int avgClubs) {
    vector<vector<int>> clubs(NumClubs);
    int joinedClubs;
    
    std::random_device rd;
    std::mt19937 rng(rd());

   

    for (int s = 0; s < students; s++) {
        int joinedClubs = randint(0, 4); // Club joining variability
    
   int assigned = 0;
   while (assigned < joinedClubs) {
        int c = randint(0, NumClubs - 1);

        //Search for duplicate club assignment
        if (std::find(clubs[c].begin(), clubs[c].end(), c) == clubs[c].end()) clubs[c].push_back(s);
           
        
        
            
        assigned++;
   }
}
    return clubs;
};


vector<pair<int,int>> generateFriendEdges(int students, 
    vector<vector<int>> dorms, vector<vector<int>> clubs,
     vector<vector<int>> classes, vector<vector<int>> schools, double scaleFactor = 0.03,
    double wsRewireProb = 0.1, int k_ws = 4, int ba_extra = 2) {

    // Initialize Watts-Strogatz friends graph: each node connects to k_ws neighors on each side
     vector<vector<int>> friends(students);
     
    // Build ring lattice
    for (int i = 0; i < students; i++) {
        for(int j = 1; j <= k_ws; j++) {
            int neighbor = (i + j) % students;
            friends[i].push_back(neighbor);
            friends[neighbor].push_back(i);
            int neighbor2 = (i - j + students) % students;
            friends[i].push_back(neighbor2);
            friends[neighbor2].push_back(i);
        }
    }
    // Rewire edges with probability wsRewireProb
    for (int i = 0; i < students; i++) {
        vector<int> neighbors(friends[i].begin(), friends[i].end());
        for (int neighbor : neighbors) {
            if (neighbor <= i) continue; // Avoid double processing
            if (rand01() < wsRewireProb) {
               // remove existing edge (erase by value using remove-erase idiom)
               friends[i].erase(std::remove(friends[i].begin(), friends[i].end(), neighbor), friends[i].end());
               friends[neighbor].erase(std::remove(friends[neighbor].begin(), friends[neighbor].end(), i), friends[neighbor].end());

               // Pick new neighbor (not including self and current neighbors)
               int newEdge = randint(0, students - 1);
               int attempts = 0;
               while ((newEdge == i || std::find(friends[i].begin(), friends[i].end(), newEdge) != friends[i].end()) && attempts < 50) {
                    newEdge = randint(0, students - 1);
                    attempts++;
               }
               if (newEdge != i && std::find(friends[i].begin(), friends[i].end(), newEdge) != friends[i].end()) {
                    friends[i].push_back(newEdge);
                    friends[newEdge].push_back(i);
            } else {
                // Re-add original edge if no valid new edge found
                friends[i].push_back(neighbor);
                friends[neighbor].push_back(i);
               }
            }
        }
    }

    // Preferential attachment based on relationships (Barabasi-Albert)
    vector<double> degrees(students);
for (int i = 0; i < students; i++) degrees[i] = (double)friends[i].size();


// Compute cumulative degree distribution
for (int t = 0; t < students; t++) {
    // Add a few random edges based on preferential attachment
    for (int e = 0; e < ba_extra; e++) {
        double sumdeg = 0.0;
        for (double d : degrees) sumdeg += d + 1.0; // +1 to avoid zero degree
        double r = rand01() * sumdeg;
        double acc = 0.0;
        int target = 0;
        for (int i = 0; i < students; i++) {
            acc += degrees[i] + 1.0;
            if (acc >= r) {
                target = i;
                break;
            }
        }
        int a = randint(0, students - 1);
        if (a != target && std::find(friends[a].begin(), friends[a].end(), target) != friends[a].end()) {
            friends[a].push_back(target);
            friends[target].push_back(a);
            degrees[a] += 1.0;
            degrees[target] += 1.0;
        }
    }

}

    // Convert to edges
    vector<pair<int,int>> edges;
    edges.reserve(students * 4);
    for (int v = 0; v < students; v++) {
        for (int e : friends[v]) {
            if (e > v) edges.emplace_back(v, e);
        }
    }
    return edges;
};

void addRelationshipGroupEdges(UnweightedGraph &unweighted,
    WeightedGraph &weighted,
    const vector<vector<int>> &groups,
    const string &relationType) {
    // Relationship weights
    int weight = W[relationType];
    int N = (int)unweighted.size();

    // For each group (dorm, club, class, school)
    for (const auto& group : groups) {
        int groupSize = group.size();

        // Add edge for all pairs (v,e) in the group
        for (int i = 0; i < groupSize; i++) {
            for (int j = i + 1; j < groupSize; j++) {
                int v = group[i];
                int e = group[j];

                // Unweighted graph
                unweighted[v].push_back(e);
                unweighted[e].push_back(v);

                // Weighted with minimization
                bool foundVE = false;
                for (auto &p : weighted[v]) {
                    if (p.first == e) {
                        p.second = min(p.second, weight);
                        foundVE = true;
                        break;
                    }
                }
                if (!foundVE) {
                    weighted[v].push_back({e, weight});
                    weighted[e].push_back({v, weight});
                }
            }
        }
    }
}

// Post Processing: remove duplicates and sort
void unique_and_sort(UnweightedGraph &G) {
    for (auto &neighbors : G) {
        sort(neighbors.begin(), neighbors.end());
        neighbors.erase(std::unique(neighbors.begin(), neighbors.end()), neighbors.end());
    }
}

void unique_and_minimize(WeightedGraph &G) {
    for (auto &neighbors : G) {
        if (neighbors.empty()) continue;
        // Sort by neighbor index
        sort(neighbors.begin(), neighbors.end(),
            [](const pair<int,int> &a, const pair<int,int> &b) {
                return a.first < b.first;
            });

        // Remove duplicates, keep minimum weight
        vector<pair<int,int>> uniqueNeighbors;
        uniqueNeighbors.reserve(neighbors.size());
        
        int curN = neighbors[0].first;
        int curW = neighbors[0].second;
        for (size_t i = 1; i < neighbors.size(); i++) {
            if (neighbors[i].first == curN) {
                curW = min(curW, neighbors[i].second);
            } else {
                uniqueNeighbors.push_back({curN, curW});
                curN = neighbors[i].first;
                curW = neighbors[i].second;
            }
        }
        uniqueNeighbors.push_back({curN, curW});
        neighbors.swap(uniqueNeighbors);
    }
}

void addFriendEdges(UnweightedGraph &unweighted,
    WeightedGraph &weighted,
    const vector<pair<int,int>> &friendEdges) {
    int weight = W["Friends"];
    
    for (const auto& edge : friendEdges) {
        int u = edge.first;
        int v = edge.second;
        
        // Unweighted graph
        unweighted[u].push_back(v);
        unweighted[v].push_back(u);
        
        // Weighted with minimization
        bool foundUV = false;
        for (auto &p : weighted[u]) {
            if (p.first == v) {
                p.second = min(p.second, weight);
                foundUV = true;
                break;
            }
        }
        if (!foundUV) {
            weighted[u].push_back({v, weight});
            weighted[v].push_back({u, weight});
        }
    }
}

void write_edges(const string &filename, const vector<pair<int,int>> &edges) {
    ofstream out(filename);
    for (const auto& edge : edges) {
        out << edge.first << "," << edge.second << "\n";
    }
    out.close();
}

// BFS (single source) - returns vector<int> dist (int = -1) and parent for path building
pair<vector<int>, vector<int>> bfs_single_source(const UnweightedGraph &Graph, int source) {
    int n = Graph.size();
    vector<int> dist(n, -1), parent(n, -1);
    queue<int> q;
    dist[source] = 0;
    q.push(source);
    while (!q.empty()) {
        int e = q.front(); 
        q.pop();

        for (int v : Graph[e]) {
            if (dist[v] == -1) {
                dist[v] = dist[e] + 1;
                parent[v] = e;
                q.push(v);
            }
        }
    }
    return {dist, parent};

};

// Reconstruct path from source to target using parent array
vector<int> reconstruct_path(int target, const vector<int> &parent) {
    vector<int> path;

if (parent.empty()) return path;

    int cur = target;
    if (cur < 0) return path;

    // If parent[cur] == -1 and cur != source, no path
    while (cur != -1) {
        path.push_back(cur);
        cur = parent[cur];
    }
    reverse(path.begin(), path.end());
    return path;
};


// Dijkstra

pair<vector<double>, vector<int>> dijkstra_single_source(const WeightedGraph &Graph, int source) {
    int n = (int)Graph.size();
    const double INF = 1e18;
    vector<double> dist(n, INF);
    vector<int> parent(n, -1);
    using PDD = pair<double, int>; // (distance, node)
    priority_queue<PDD, vector<PDD>, greater<PDD>> pq;

    dist[source] = 0.0;
    pq.push({0.0, source});

    while (!pq.empty()) {
        auto [curDist, e] = pq.top();
        pq.pop();

        if (curDist > dist[e]) continue; // stale entry

        for (const auto& [v, weight] : Graph[e]) {
            double newDist = curDist + (double)weight;
            if (newDist < dist[v]) {
                dist[v] = newDist;
                parent[v] = e;
                pq.push({newDist, v});
            }
        }
    }
    return {dist, parent};
};

// DKA on induced subgraph when given a set of nodes 
pair<vector<double>, vector<int>> dijkstra_subgraph(const WeightedGraph &Graph, int source, int target, const vector<int> &nodes) {
    int n = (int)nodes.size();
    unordered_map<int,int> nodeToIndex;
    nodeToIndex.reserve(n*2);
    for (int i = 0; i < n; i++) nodeToIndex[nodes[i]] = i;
    
    // Initialize subgraph
    WeightedGraph subG(n);

    // Build subgraph
    for (int i = 0; i < n; i++) {
        int u = nodes[i];
        for (auto &[v, weight] : Graph[u]) {
            if (nodeToIndex.find(v) != nodeToIndex.end()) {
                subG[i].push_back({nodeToIndex[v], weight});
            }
        }
    }
    int s = nodeToIndex.count(source) ? nodeToIndex[source] : -1;
    int t = nodeToIndex.count(target) ? nodeToIndex[target] : -1;
    double INF = 1e18;
    vector<double> dist(n, INF);
    vector<int> parent(n, -1);
    using PDD = pair<double, int>; // (distance, node)
    priority_queue<PDD, vector<PDD>, greater<PDD>> pq;
    if (s == -1) return {vector<double>(), vector<int>()}; // Start outside of induced set
    dist[s] = 0.0; pq.push({0.0, s});
    while (!pq.empty()) {
         auto [curDist, e] = pq.top(); pq.pop();
        if (curDist > dist[e]) continue; // stale entry
        for (auto& [v, weight] : subG[e]) {
            double newDist = curDist + (double)weight;
            if (newDist < dist[v]) {
                dist[v] = newDist;
                parent[v] = e;
                pq.push({newDist, v});
            }
        }
    }

    return {dist, parent};

}
// Hybrid shortest path

// Metrics
 // Approx. average path length via BFS sampling from random seeds
double approximateAvgPathLength(const UnweightedGraph &Graph, int sampleSeeds) {
    int n = (int)Graph.size();
    if (n== 0) return 0.0;
    long long total = 0;
    long long count = 0;
    

    for (int i = 0; i < sampleSeeds; i++) {
        int src = randint(0, n - 1);
        auto [dist, parent] = bfs_single_source(Graph, src);
        for (int v = 0; v < n; v++) {
            if (v == src) continue;
            if (dist[v] >= 0) {
                 total =+ dist[v];
                  count++; 
            }
        }
    }
    if (count == 0) return 0.0;
    return (double)total / (double)count;
};

// Approx. diameter via double sweep
int approximateDiameter(const UnweightedGraph &Graph) {
    int n = (int)Graph.size();
    if (n == 0) return 0;
    int source = randint(0, n - 1);
    auto [dist1, parent1] = bfs_single_source(Graph, source);
    // Find farthest from source
    int farthest = source;
    for (int i = 0; i < n; i ++) {
        if (dist1[i] > dist1[farthest]) {
            farthest = i;
        }
    }
    auto [dist2, parent2] = bfs_single_source(Graph, farthest);
    int farthest2 = farthest;
    for (int i = 0; i < n; i++) {
        if (dist2[i] > dist2[farthest2]) {
            farthest2 = i;
        }
    }
    return dist2[farthest2];
    };

// Average degree and degree distribution
pair<double, vector<int>> degree_stats(const UnweightedGraph &Graph) {
    int n = (int)Graph.size();
    vector<int> deg(n);
    long long sum = 0;
    for (int i = 0; i < n; i++) {
        deg[i] = (int)Graph[i].size();
        sum += deg[i];
    }
    double avg = (double)sum / (double)n;
    cout << "n = " << n << ". \n" << "sum = " << sum << " . \n";
    return {avg, deg};
};

// Clustering coefficient (average local)
double clustering_coefficient(const UnweightedGraph &Graph) {
    int n = (int)Graph.size();
    double totalC = 0.0;
    for (int e = 0; e < n; e++) {
        int k = (int)Graph[e].size();
        if (k < 2) continue; // C(e) = 0 if degree < 2
        // Count edges between neighbors
        unordered_set<int> neighborSet(Graph[e].begin(), Graph[e].end());
        int links = 0;
        for (int v : Graph[e]) {
            for (int w : Graph[v]) {
                if (v < w && neighborSet.count(w)) {
                    links++; // each link is counted once if v<w
                }
            }
        }
        // C(e) = 2*links / (k*(k-1))
        double possible = (double)k * (k - 1) / 2.0;

        int actual = 0;
        for (int i = 0; i < k; i++) {
            for (int j = i + 1; j < k; j++) {
                int a = Graph[e][i];
                int b = Graph[e][j];
                // If a and b are connected
                // check smaller degree first
                // brute force
                if(find(Graph[a].begin(), Graph[a].end(), b) != Graph[a].end()) {
                    actual++;
                }
            }         
        }
            totalC += (possible > 0.0) ? (double)actual / possible : 0.0;
    }

    return totalC / (double)n;
};

vector<double> betweenness_centrality(const UnweightedGraph &Graph) {
    int n = (int)Graph.size();
    vector<double> CB(n, 0.0);
    for (int s = 0; s < n; s++) {
        // Single-source shortest paths
        vector<vector<int>> pred(n);
        vector<int> dist(n, -1);
        vector<double> sigma(n, 0.0);
        queue<int> q;
        vector<int> stack;
        dist[s] = 0;
        sigma[s] = 1.0;
        q.push(s);

        while (!q.empty()) {
            int v = q.front(); q.pop();
            stack.push_back(v);
            for (int w : Graph[v]) {
                // Found for the first time
                if (dist[w] < 0) {
                    dist[w] = dist[v] + 1;
                    q.push(w);
                }
                // Found shortest path to w via v
                if (dist[w] == dist[v] + 1) {
                    sigma[w] += sigma[v];
                    pred[w].push_back(v);
                }
            }
        }
        vector <double> delta(n, 0.0);
        for (int i = (int)stack.size() - 1; i >= 0; i--) {
            int w = stack[i];
            for (int v : pred[w]) {
                if (sigma[w] != 0) {
                    delta[v] += (sigma[v] / sigma[w]) * (1.0 + delta[w]);
                }
            }
            if (w != s) CB[w] += delta[w];
        }
    }
    return CB;
};


int main() {    

    cout << "Campus Network Simulation\n";
    cout << "N= " << students << ", Dorms= " << dormAmount << ", Class cap= " << classCap <<", Avg classes= " << avgClasses << "\n";
    cout << "Generating relationship groups...\n";

    // Assign relationships
    auto classes = assignClasses(students, classCap, avgClasses);
    auto dorms = assignDorms(students, dormAmount);
    auto schools = assignSchools(students, NumSchools);
    auto clubs = assignClubs(students, NumClubs, avgClubs);
    auto friends = WS_model(students, 6, 0.2);
    
    

    // CSV outputs for verification
   // write_edges("friend_edges.csv", friends);

    // Initialize graphs
    UnweightedGraph unwg(students);
    WeightedGraph wg(students);
    // Add groups to graphs

    // Print a summary of the schools vector (cannot stream vector-of-vectors directly)
    cout << "schools: " << schools.size() << " groups\n";
    for (size_t si = 0; si < schools.size(); ++si) {
        cout << "  school " << si << " size: " << schools[si].size() << "\n";
    }
    
    addRelationshipGroupEdges(unwg, wg, dorms, "Dorms");
    addRelationshipGroupEdges(unwg, wg, classes, "Classes");
    addRelationshipGroupEdges(unwg, wg, clubs, "Clubs");
    addSchoolEdges_WS(unwg, wg, schools,
    6,        // k = each student linked to 6 neighbors (3 on each side)
    0.1,      // beta = 0.1 rewiring probability
    5         // school relationship weight
);

    //addRelationshipGroupEdges(unwg, wg, schools, "Schools");
    addFriendEdges(unwg, wg, friends);
    

    // Deduplicate adjacency lists
    unique_and_sort(unwg);
    unique_and_minimize(wg);

    cout << "Graph construction complete. Nodes: " << students << "\n";

    // Metrics
    auto [avgDeg, deg] = degree_stats(unwg);
    cout << "Average Degree (unweighted): " << avgDeg << "\n";

    // Write degree distribution (CSV)


    // Clustering Coefficient
    double clustering = clustering_coefficient(unwg);
    cout << "Avg clustering coefficient: " << clustering << "\n";


    // Approx. Average Path Length
    double avgPath = approximateAvgPathLength(unwg, SampleSeedsForPaths);
    cout << "Approx diameter: " << avgPath << "\n";


    // Approx Diameter
    int diameter = approximateDiameter(unwg);
    cout << "Approximate diameter: " << diameter << "\n";

    // Betweenness centrality



    // Community detection (label propagation & modularity)


    // Runtime comparisons


    // Write summary CSV
    write_edges("../StetsonSmallWorld/data/watts_strogatz_edges.csv", friends);
    
};