CarSharingSimulation.cpp

#include <algorithm>
#include <chrono>
#include <iostream>

#include "CarSharingSimulation.hpp"

std::vector<double> sample_time_points(double end_time, double dt) {
    std::vector<double> time_points;
    
    for (double t = 0.0; t < end_time; t += dt) {
        time_points.push_back(t);
    }
    
    return time_points;
}

CarSharingSimulation::CarSharingSimulation(
    uint32_t num_stations,
    uint32_t cars_per_station,
    uint32_t spaces_per_station,
    double prob_class_1,
    double lam,
    double nu,
    double mu)
    : num_stations(num_stations),
      cars_per_station(cars_per_station),
      spaces_per_station(spaces_per_station),
      prob_class_1(prob_class_1),
      lam(lam),
      nu(nu),
      mu(mu),
      class1_reservations(num_stations, 0),
      class1_driving(num_stations, 0),
      class1_waiting(num_stations, 0),
      available_cars(num_stations, cars_per_station),
      class2_driving(num_stations, 0),
      unavailable_spaces(num_stations, cars_per_station),
      _class2_reservations(num_stations * num_stations, 0),
      stats(),
      transition_times(),
      num_reservations_1(0),
      num_reservations_2(0),
      num_driving_1(0),
      num_driving_2(0),
      num_waiting(0),
      num_stats(0),
      rng(static_cast<uint32_t>(
          std::chrono::high_resolution_clock::now().time_since_epoch().count()))
{
}


uint32_t CarSharingSimulation::uniform_random_station() {
    std::uniform_int_distribution<uint32_t> dist(0, num_stations - 1);
    return dist(rng);
}

/* 
 * Class 1 reservation: attempts to make a reservation at a station chosen uniformly at random.
 * The reservation is successful if and only if the station has an available car.
 * @return true if the reservation succeeded; false otherwise.
 */
bool CarSharingSimulation::reservation_update_1() {
    uint32_t station = uniform_random_station();
    if (available_cars.at(station) == 0) {
        return false;
    } else {
        available_cars.at(station) -= 1;
        class1_reservations.at(station) += 1;
        return true;
    }
}

/* 
 * Class 2 reservation: attempts to reserve a car at a starting station and a parking space at
 * a destination station, both chosen uniformly at random..
 * The reservation is successful if and only if the starting station has an available car and
 * the destination station has an available parking space.
 * @return true if the reservation succeeded; false otherwise.
 */
bool CarSharingSimulation::reservation_update_2() {
    uint32_t station = uniform_random_station();
    uint32_t destination = uniform_random_station();
    if (available_cars.at(station) == 0 || unavailable_spaces.at(destination) >= spaces_per_station) {
        return false;
    } else {
        available_cars.at(station) -= 1;
        class2_reservations(station, destination) += 1;
        unavailable_spaces.at(destination) += 1;
        return true;
    }
}

/*
 * Class 1 pickup: a class 1 user with a car reserved at a starting station
 * picks up the car and begins driving to a randomly selected destination station.
 */
void CarSharingSimulation::pickup_update_1() {
    uint32_t destination = uniform_random_station();
    auto start_dist = class1_reservation_distribution();
    uint32_t station = start_dist(rng);
    class1_reservations.at(station) -= 1;
    class1_driving.at(destination) += 1;
    if (class1_waiting.at(station) > 0) {
        class1_waiting.at(station) -= 1;
        num_waiting -= 1;
        available_cars.at(station) += 1;
    } else {
        unavailable_spaces.at(station) -= 1;
    }
}

/*
 * Class 1 pickup: a class 1 user with a car reserved at a starting station
 * picks up the car and begins driving to their previously reserved destination station.
 */
void CarSharingSimulation::pickup_update_2() {
    auto c2r_dist = class2_reservation_distribution();
    uint32_t index = c2r_dist(rng);
    uint32_t station = index / num_stations;
    uint32_t destination = index % num_stations;
    class2_reservations(station, destination) -= 1;
    class2_driving.at(destination) += 1;
    if (class1_waiting.at(station) > 0) {
        class1_waiting.at(station) -= 1;
        num_waiting -= 1;
        available_cars.at(station) += 1;
    } else {
        unavailable_spaces.at(station) -= 1;
    }
}

/*
 * Class 1 dropoff: a user of class 1 who is driving with a chosen destination
 * attempts to park at that destination. If there is a parking space available,
 * that user parks; otherwise, they are added to the queue of users waiting to
 * * park at that station.
 */
void CarSharingSimulation::dropoff_update_1() {
    auto dist = class1_driving_distribution();
    uint32_t destination = dist(rng);
    class1_driving.at(destination) -= 1;
    if (unavailable_spaces.at(destination) >= spaces_per_station) {
        class1_waiting.at(destination) += 1;
        num_waiting += 1;
    } else {
        available_cars.at(destination) += 1;
        unavailable_spaces.at(destination) += 1;
    }
}

/*
 * Class 2 dropoff: a user of class 2 who is driving toward their destination 
 * parks at that destination station (Note: the user already has a parking space
 * reserved, so this will never result in additional waiting).
 */
void CarSharingSimulation::dropoff_update_2() {
    auto dist = class2_driving_distribution();
    uint32_t destination = dist(rng);
    class2_driving.at(destination) -= 1;
    available_cars.at(destination) += 1;
}

/*
 * Determines the next transition using weights determined by the current
 * state of the system.
 * 
 * @return a value of type CarSharingSimulation::Transition specifying the next event.
 */
CarSharingSimulation::Transition CarSharingSimulation::next_transition() {
    std::discrete_distribution<uint32_t> dist({
        prob_class_1 * lam,
        (1 - prob_class_1) * lam,
        num_reservations_1 * nu,
        num_reservations_2 * nu,
        num_driving_1 * mu,
        num_driving_2 * mu
    });
    return static_cast<Transition>(dist(rng));
}

/*
 * Draws the time until the next event in the simulation.
 * The event holding time is sampled from an exponential distribution
 * with a rate equal to the sum of all possible event rates.
 *
 * @return Time until the next event.
 */
double CarSharingSimulation::next_holding_time() {
    std::exponential_distribution<double> dist(lam +
        (num_reservations_1 + num_reservations_2) * nu + 
        (num_driving_1 + num_driving_2) * mu);
    return dist(rng);
}

/**
 * Logs the current simulation time and a corresponding metric (reservations or drivers).
 * 
 * @param t The current simulation time.
 * @param column The metric to log:
 *        1 = class 1 reservations,
 *        2 = class 2 reservations,
 *        3 = class 1 drivers,
 *        4 = class 2 drivers.
 *
 * @throws std::invalid_argument if the column index is unsupported.
 */
void CarSharingSimulation::push_time_and_column(double t, uint32_t column) {
    transition_times.push_back(t);
    switch (column) {
    case 1:
        stats.push_back(num_reservations_1);
        break;
    case 2:
        stats.push_back(num_reservations_2);
        // std::cout << t << ' ' << num_reservations_2 << std::endl;
        break;
    case 3:
        stats.push_back(num_driving_1);
        break;
    case 4:
        stats.push_back(num_driving_2);
        break;
    case 5:
        stats.push_back(num_waiting);
        break;
    default:
        throw std::invalid_argument("Unsupported column index: " + std::to_string(column));
    }
}

/*
 * Main loop of the simulation.
 * 
 * @param end_time the time at which to stop running the simulation.
 * @param column The metric to log:
 *        1 = class 1 reservations,
 *        2 = class 2 reservations,
 *        3 = class 1 drivers,
 *        4 = class 2 drivers.
 */
void CarSharingSimulation::run(double end_time, uint32_t column) {
    double t = 0.0;
    while (t < end_time) {
        push_time_and_column(t, column);
        Transition transition = next_transition();
        double holding_time = next_holding_time();
        if (transition == Transition::CLASS1_RESERVATION) {
            if (!reservation_update_1()) {
                t += holding_time;
                continue;
            }
            num_reservations_1 += 1;
        } else if (transition == Transition::CLASS2_RESERVATION) {
            if (!reservation_update_2()) {
                t += holding_time;
                continue;
            }
            num_reservations_2 += 1;
        } else if (transition == Transition::CLASS1_PICKUP) {
            pickup_update_1();
            num_reservations_1 -= 1;
            num_driving_1 += 1;
        } else if (transition == Transition::CLASS2_PICKUP) {
            pickup_update_2();
            num_reservations_2 -= 1;
            num_driving_2 += 1;
        } else if (transition == Transition::CLASS1_DROPOFF) {
            dropoff_update_1();
            num_driving_1 -= 1;
        } else if (transition == Transition::CLASS2_DROPOFF) {
            dropoff_update_2();
            num_driving_2 -= 1;
        }
        t += holding_time;
    }
    ++num_stats;
    push_time_and_column(t, column);
}

void CarSharingSimulation::average(double end_time, double dt) {
    size_t transition_index = 0;
    for (size_t i = 0; i < cumulative_times.size(); i++) {
        while (transition_index + 1 < transition_times.size() &&
            cumulative_times.at(i) > transition_times.at(transition_index + 1))
        {
            ++transition_index;
        }
        double avg_stat = (cumulative_stats.at(i) * (num_stats - 1) +
            stats.at(transition_index)) / num_stats;
        cumulative_stats.at(i) = avg_stat;
    }
}

/*
 * Reset the simulation back to its initial state, preparing for another run.
 */
void CarSharingSimulation::reset_state() {
    std::fill(class1_reservations.begin(), class1_reservations.end(), 0);
    std::fill(_class2_reservations.begin(), _class2_reservations.end(), 0);
    std::fill(available_cars.begin(), available_cars.end(), cars_per_station);
    std::fill(class2_driving.begin(), class2_driving.end(), 0);
    std::fill(class1_waiting.begin(), class1_waiting.end(), 0);
    std::fill(unavailable_spaces.begin(), unavailable_spaces.end(), cars_per_station);

    num_reservations_1 = 0;
    num_reservations_2 = 0;
    num_driving_1 = 0;
    num_driving_2 = 0;

    stats.clear();
    transition_times.clear();

}

/**
 * Run the simulation repeatedly.
 * 
 * @param num_runs the number of times to run the simulation.
 * @param end_time the time at which to stop running each simulation.
 * @param dt the time step for averaging.
 * @param column the performance metric to output.
 * 
 * @return a vector of values from the specified metric at each time point.
 */
std::vector<double> CarSharingSimulation::repeated_runs(
    uint32_t num_runs,
    double end_time,
    double dt,
    uint32_t column
) {
    
    cumulative_times = sample_time_points(end_time, dt);
    cumulative_stats = std::vector<double>(cumulative_times.size(), 0.0);
    num_stats = 0;

    for (uint32_t i = 0; i < num_runs; i++) {
        reset_state();
        run(end_time, column);
        average(end_time, dt);
    }
    return cumulative_stats;
}

/**
 * Run the simulation repeatedly, averaging results.
 * This function is called by plotsim.py.
 * It instantiates a CarSharingSimulation object and calls
 * its repeated_runs method.
 */
std::vector<double> repeated_runs(
    uint32_t num_stations,
    uint32_t cars_per_station,
    uint32_t spaces_per_station,
    double prob_class_1,
    double lam,
    double nu,
    double mu,
    uint32_t num_runs,
    double end_time,
    double dt,
    uint32_t column
) {

    CarSharingSimulation sim(num_stations, cars_per_station, spaces_per_station,
        prob_class_1, lam, nu, mu);
    return sim.repeated_runs(num_runs, end_time, dt, column);
}