<jupyter_start><jupyter_text>Discrete Simulation Final ProjectAuthors: [<NAME>](https://github.com/FelipOliveira), [<NAME>](https://github.com/jvcanavarro) and [<NAME>](https://github.com/VitorCantao). *** Metro Bike [Dataset](https://bikeshare.metro.net/about/data/)*The Metro Bike Share system makes bikes available 24/7, 365 days a year in Downtown LA, Central LA, and North Hollywood (...). Metro Bike Share offers convenient round-the-clock access to a fleet of bicycles for short trips. Metro Bike Share is one of LA Metro's multiple public transportation options for Angelenos and visitors to get around.*<jupyter_code>import matplotlib.pyplot as plt
import geopy.distance as gd
import pandas as pd
import numpy as np<jupyter_output><empty_output><jupyter_text>Dataset Overview<jupyter_code>df = pd.read_csv('../datasets/metrobike-2020.csv')
df.head()
# Number of rides in 2020
df.shape[0]
df.columns
# Remove unecessary columns
columns = ['duration', 'start_time', 'end_time', 'start_station', 'start_lat', 'start_lon', 'end_station', 'end_lat', 'end_lon', 'bike_id', 'trip_route_category', 'bike_type']
df = df[columns]
df.head()<jupyter_output><empty_output><jupyter_text>According to Stations [Dataset](https://bikeshare.metro.net/wp-content/uploads/2021/01/metro-bike-share-stations-2021-01-01.csv):* `4285,Metro Bike Share Free Bikes,2/27/2019,Free Bikes,Active`* `4286,Metro Bike Share Out of Service Area Smart Bike,2/27/2019,Free Bikes,Active`* `3000,Virtual Station,7/7/2016,N/A,Active`Those stations doesn't possess geolocation information, which leads us to ignore then when calculating rides distances<jupyter_code>df_null = df[df.isnull().any(axis=1)]
nan_stations = df_null.start_station[df_null['start_lat'].isnull()].unique()
nan_stations
# Remove stations without geolocation
df = df.dropna()
# Remove trips starting and ending at the same station
df = df.loc[df.trip_route_category == 'One Way']
# Total bikes within all stations and Nº of stations
print('Number of bicycles: ', len(df['bike_id'].unique()))
print('Number of stations: ', df['start_station'].unique().shape[0])
# Percentage of rides by each model of bicycle
df['bike_type'].value_counts(normalize=True)
# Mean duration of a bike travel
df['duration'].mean()<jupyter_output><empty_output><jupyter_text>Calculate distance using geolocation coordenates<jupyter_code>def get_distance(coords):
    return gd.geodesic((coords[0], coords[1]), (coords[2], coords[3])).km
# Create new distance column
columns = ['start_lat', 'start_lon', 'end_lat', 'end_lon']
df['distance'] = df.apply(lambda x: get_distance(x[columns]), axis=1)
df.head()<jupyter_output><empty_output><jupyter_text>A better option was to filter unique pairs of start and end stations, then calculate the distances. It may be implemented in the future.<jupyter_code>df[['duration', 'distance']].describe()<jupyter_output><empty_output><jupyter_text>The duration and distance values of rides are very discrepant (see the standard deviation). Eventually, there are cases where people stay with the bike an entire day and then return then to a station, which doesn't reflect the real time they spending at the ride.In this sense, we decided to ignore those cases when calculating the mean speed of each bycicle model.  Calculate the average speed of each bicycle model with different intervals of time and space<jupyter_code># Select only necessary columns
time_space = df.iloc[:, [0, 11, 12]]
time_space.head()
def get_speed(space, time, unit):

    if unit == 'km/min': return space / time
    if unit == 'km/h'  : return space / (time / 60)
    if unit == 'm/s'   : return (space * 1000) / (time * 60)
    
unit = 'km/min'
intervals = [(0, 0), (15, 3), (30, 6), (45, 9), (120, 12)]
distance = 0
avg_speed = pd.DataFrame()

for x, y in zip(intervals, intervals[1:]):
    title = f'{x[0]} - {y[0]} Min'
    # Loc the intervals of duration and distance in order to calculate the speed of each model at these situations
    time_space_filtered = time_space.loc[(df.duration > x[0]) & (df.duration <= y[0]) 
                                           & (df.distance > x[1]) & (df.distance <= y[1])]
    time_space_filtered['speed'] = time_space_filtered.apply(lambda x: get_speed(x['distance'], x['duration'], unit=unit), axis=1)
    
    bike_speed = time_space_filtered.groupby('bike_type').mean()
    avg_speed[title] = bike_speed['speed']
    distance += 3

avg_speed
# Calculate average speed of each bike model

time_space_filtered = time_space.loc[(df.duration > 0) & (df.duration <= 180) 
                                           & (df.distance > 0) & (df.distance <= 12)]

time_space['speed'] = time_space_filtered.apply(lambda x: get_speed(x.distance, x.duration, unit), axis=1)
bike_speed = time_space.groupby('bike_type').mean()
avg_speed['Overall'] = bike_speed['speed']<jupyter_output><empty_output><jupyter_text>We are only considering trips with duration and distance up to 120 minutes and 12 kilometers, respectively.Idea: Create a 1st histogram that represents the normal distribuition of trips durationThen, create a 2nd histogram that represents each bike model speed based on trips duration<jupyter_code>title = 'Average speed of bicycle models for different intervals of time and travel distances'
xlabel = 'Duration(min) | Distance(km) Intervals'

ax = avg_speed.T.plot(kind='bar', figsize=(17, 8), xlabel=xlabel, ylabel='Speed in ' + unit, colormap='tab20b', rot=0, title=title)
for p in ax.patches:                 
    ax.annotate(round(p.get_height(),3), (p.get_x()+p.get_width()/2., p.get_height()), ha='center', va='center', xytext=(0, 10), textcoords='offset points')
ax.grid(linestyle=':', axis='y')<jupyter_output><empty_output><jupyter_text>it is noticeable that electric bicycles have a certain "advantage" over ordinary ones as the distance and duration of the rides increases  Normalize the frequency of rides for each unique pair of (start|end) station<jupyter_code># Ignore rides that end at same station and calculate total ocurrences by each ride
group = ['start_station', 'end_station', 'distance']
stations = df.groupby(group).size().reset_index(name='occurs')

# Normalize weights by each group (i.e. each start_station)
stations['norm_rides'] = stations['occurs'] / stations.groupby('start_station')['occurs'].transform('sum')
stations = stations.drop_duplicates(['start_station', 'end_station'])
stations<jupyter_output><empty_output><jupyter_text>Create frequency and distance matrices<jupyter_code># Pivot to create rides frequency matrix
freq_matrix = stations.pivot(index='start_station', columns='end_station', values='norm_rides').fillna(0)
freq_matrix
# Stations distance matrix
distance_matrix = stations.pivot(index='start_station', columns='end_station', values='distance').fillna(0)
distance_matrix.head()
# Normalize start stations
start_station_matrix = df.start_station.value_counts(normalize=True)
start_station_matrix.head()
avg_speed.at['standard', (15, 3)]<jupyter_output><empty_output><jupyter_text>*** Simulation<jupyter_code>def reset():
    stations = []
    total_clients = 0
    total_wait_time = 0
    total_overflow = 0
from enum import Enum
import simpy
import random

STATIONS = freq_matrix.index
destinations_matrix = freq_matrix.loc[STATIONS, STATIONS]
BIKE_TYPE = df['bike_type'].unique()

NUMBER_OF_RIDES = 209974
RIDES_PER_DAY = round(NUMBER_OF_RIDES / 365)
NUMBER_OF_BIKES = 3908
NUMBER_OF_STATIONS = len(STATIONS)
NUMBER_OF_BIKES_PER_STATION = NUMBER_OF_BIKES // NUMBER_OF_STATIONS
INTERVALS = [(0, 0), (15, 3), (30, 6), (45, 9), (120, 12)]

SHOULD_PRINT_EVENT = False

def random_weighted(vector):
    PRECISION = 10**6
    integer_vector = round(vector * PRECISION)
    pool_size = integer_vector.sum()

    random_number = random.randrange(0, pool_size)

    accumulated_probability = 0
    for index, item in enumerate(integer_vector):
        accumulated_probability += item
        if (random_number < accumulated_probability):
            return index

class BikeType(Enum):
    STANDARD = 1
    ELECTRIC = 2
    SMART = 3

class BikeDock:
    def __init__(self, env, standard_bikes, electric_bikes, smart_bikes, capacity):
        self.env = env
        self._standard_dock = simpy.Container(self.env, init=standard_bikes, capacity=capacity)
        self._electric_dock = simpy.Container(self.env, init=electric_bikes, capacity=capacity)
        self._smart_dock = simpy.Container(self.env, init=smart_bikes, capacity=capacity)
        self._type_to_dock = {BikeType.STANDARD: self._standard_dock, BikeType.ELECTRIC: self._electric_dock, BikeType.SMART: self._smart_dock}

    def capacity(self, bike_type):
        return self._type_to_dock[bike_type].capacity

    def available_bikes(self, bike_type):
        return self._type_to_dock[bike_type].level

    def has_available(self, bike_type):
        return self._type_to_dock[bike_type].level > 0

    def is_full(self, bike_type):
        return self._type_to_dock[bike_type].level >= self._type_to_dock[bike_type].capacity

    def rent_bike(self, bike_type):
        return self._type_to_dock[bike_type].get(1)

    def return_bike(self, bike_type):
        return self._type_to_dock[bike_type].put(1)


class Ride:
    def __init__(self, env, start_station, end_station, bike_type):
        self.start_station = start_station
        self.end_station = end_station
        self.bike_type = bike_type
        self._type_to_speed = {BikeType.STANDARD: 0.127, BikeType.ELECTRIC: 0.12, BikeType.SMART: 0.11}
        self.env = env

    @property
    def distance(self):
        return distance_matrix.at[self.start_station.id, self.end_station.id]

    @property
    def duration(self):
        return round(self.distance / self._bike_speed)

    def start(self):
        env.process(self.start_station.get_bike(self.bike_type))
        yield env.timeout(self.duration)
        env.process(self.end_station.put_bike(self.bike_type))

    @property
    def _bike_speed(self):
        bike_model = self.bike_type.name.lower()
        for x, y in zip(INTERVALS, INTERVALS[1:]):
            if x[0] < self.distance <= y[0]:
                return avg_speed.at[bike_model, y]
        return avg_speed.at[bike_model, 'overall']


class Station(object):
    def __init__(self, env, station_id):
        self.env = env
        self.id = station_id

        standard_bikes = round(NUMBER_OF_BIKES_PER_STATION * 0.4)
        electric_bikes = round(NUMBER_OF_BIKES_PER_STATION * 0.4)
        smart_bikes = round(NUMBER_OF_BIKES_PER_STATION * 0.2)
        self.dock = BikeDock(env, standard_bikes, electric_bikes, smart_bikes, 
                             capacity=NUMBER_OF_BIKES_PER_STATION * 1.05)

        self.overflow_count = 0
        self.client_wait_time = 0
        self.client_giveup = 0
        self.rides = 0

    def get_bike(self, bike_type):
        wait_time = 0
        while (wait_time < 10):
            if (self.dock.has_available(bike_type)):
                if (SHOULD_PRINT_EVENT): print(f'''Bike requested at station {self.id} - 
                {self.dock.available_bikes(bike_type) - 1} / {self.dock.capacity(bike_type)}''')

                self.rides += 1
                return self.dock.rent_bike(bike_type)
            else: 
                wait_time += 1
                self.client_wait_time += 1
                if (SHOULD_PRINT_EVENT): print(f'No more bikes at station {self.id}. Waiting {wait_time}...')
                yield env.timeout(1)
        self.client_giveup += 1
        if (SHOULD_PRINT_EVENT): print('#-#-#  I\'m out of here! Trash company!  #-#-#')
        

    def put_bike(self, bike_type):
        if (self.dock.is_full(bike_type)):
            self.overflow_count += 1

        if (SHOULD_PRINT_EVENT): print(f'''Bike arrived at station {self.id} - 
        {self.dock.available_bikes(bike_type) + 1} / {self.dock.capacity(bike_type)}''')

        yield self.dock.return_bike(bike_type)

def cyclist_arrivals(env, stations):
    random_station_index = random.randint(0, len(stations) - 1)
    start_station = stations[random_station_index]

    random_destination_index = random_weighted(destinations_matrix.loc[start_station.id])
    final_station = stations[random_destination_index]
    if (SHOULD_PRINT_EVENT): print(f'New Ciclist arrives! {start_station.id} -> {final_station.id}')

    bike_probability = df['bike_type'].value_counts(normalize=True)
    preferred_bike_type = BikeType(random_weighted(bike_probability) + 1)
    if (SHOULD_PRINT_EVENT): print(f'Preferred Bike Type: {preferred_bike_type}')

    ride = Ride(env, start_station, final_station, preferred_bike_type)
    yield env.process(ride.start())


def setup(env):
    global total_clients
    for station_id in STATIONS:
        stations.append(Station(env, station_id))
    while True:
        yield env.timeout(random.randint(1, 2))
        total_clients += 1
        env.process(cyclist_arrivals(env, stations))



total_clients = 0
total_wait_time = 0
total_overflow = 0
for i in range(30):
    print('Iteration: ',i+1)
    stations = []
    
    env = simpy.Environment()
    env.process(setup(env))
    env.run(until=1440)

    for station in stations:
        total_wait_time += station.client_wait_time
        total_overflow += station.overflow_count


    print(f'Total wait time: {total_wait_time}')
    print(f'Overflows: {total_overflow}')
    print(f'Trips: {total_clients}')
# Initial Simulation
ITERATIONS = 30

hipos = []
h1 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.5, .3, .2)
arrival_interval = 2.5

h1.append(mean_trips)
h1.append(n_bikes)
h1.append(mean_overflows)
h1.append(iterations)
h1.append(mean_w_time)
h1.append(arrival_interval)
h1.append(distr_bikes)

h1
hipos = []
hipos.append(h1)
# Reducing nº of bikes by half
h2 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.5, .3, .2)
arrival_interval = 2.5

h2.append(mean_trips)
h2.append(n_bikes)
h2.append(mean_overflows)
h2.append(iterations)
h2.append(mean_w_time)
h2.append(arrival_interval)
h2.append(distr_bikes)

hipos.append(h2)
print(h1)
print(h2)
print(h3)
print(h4)
print(h5)
print(h6)
print(h7)
print(h8)
print(h9)
print(h10)
print(h11)
print(h12)
print(h13)
# Reducing wait time with reduced number of bikes by changing bike model distr
h3 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.4, .4, .2)
arrival_interval = 2.5

h3.append(mean_trips)
h3.append(n_bikes)
h3.append(mean_overflows)
h3.append(iterations)
h3.append(mean_w_time)
h3.append(arrival_interval)
h3.append(distr_bikes)

hipos.append(h3)
# Reducing wait time with reduced number of bikes by changing bike model distr (2)
h4 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.3, .5, .2)
arrival_interval = 2.5

h4.append(mean_trips)
h4.append(n_bikes)
h4.append(mean_overflows)
h4.append(iterations)
h4.append(mean_w_time)
h4.append(arrival_interval)
h4.append(distr_bikes)

hipos.append(h4)
# Reducing wait time with reduced number of bikes by changing bike model distr and capacity by 20 %(3)
h5 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.3, .5, .2)
arrival_interval = 2.5

h5.append(mean_trips)
h5.append(n_bikes)
h5.append(mean_overflows)
h5.append(iterations)
h5.append(mean_w_time)
h5.append(arrival_interval)
h5.append(distr_bikes)

hipos.append(h5)
# Reducing wait time with reduced number of bikes by changing bike model distr and reducing arrival interval
h6 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.4, .4, .2)
arrival_interval = 2.0

h6.append(mean_trips)
h6.append(n_bikes)
h6.append(mean_overflows)
h6.append(iterations)
h6.append(mean_w_time)
h6.append(arrival_interval)
h6.append(distr_bikes)

hipos.append(h6)
# Reducing wait time with reduced number of bikes by changing bike model distr and reducing arrival interval (2)
h7 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.4, .4, .2)
arrival_interval = 1.5

h7.append(mean_trips)
h7.append(n_bikes)
h7.append(mean_overflows)
h7.append(iterations)
h7.append(mean_w_time)
h7.append(arrival_interval)
h7.append(distr_bikes)

hipos.append(h7)
# Reducing wait time with reduced number of bikes by changing bike model distr by incrementing smart type
h8 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.3, .4, .3)
arrival_interval = 2.5

h8.append(mean_trips)
h8.append(n_bikes)
h8.append(mean_overflows)
h8.append(iterations)
h8.append(mean_w_time)
h8.append(arrival_interval)
h8.append(distr_bikes)

hipos.append(h8)
# System Capacity: Reducing arrival interval (1.5)
h9 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.5, .3, .2)
arrival_interval = 1.5

h9.append(mean_trips)
h9.append(n_bikes)
h9.append(mean_overflows)
h9.append(iterations)
h9.append(mean_w_time)
h9.append(arrival_interval)
h9.append(distr_bikes)

hipos.append(h9)
# System Capacity: Reducing arrival interval (1)
h10 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.5, .3, .2)
arrival_interval = 1.0

h10.append(mean_trips)
h10.append(n_bikes)
h10.append(mean_overflows)
h10.append(iterations)
h10.append(mean_w_time)
h10.append(arrival_interval)
h10.append(distr_bikes)

hipos.append(h10)
# System Capacity: Reducing arrival interval (1) with the best bike distr
h11 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.4, .4, .2)
arrival_interval = 1.0

h11.append(mean_trips)
h11.append(n_bikes)
h11.append(mean_overflows)
h11.append(iterations)
h11.append(mean_w_time)
h11.append(arrival_interval)
h11.append(distr_bikes)

hipos.append(h11)
# System Capacity: Reducing arrival interval (1.5) with the best bike distr
h13 = []

mean_trips = total_clients / ITERATIONS
n_bikes = NUMBER_OF_BIKES
mean_overflows = total_overflow / ITERATIONS
iterations = ITERATIONS
mean_w_time = total_wait_time / ITERATIONS
distr_bikes = (.4, .4, .2)
arrival_interval = 1.5

h13.append(mean_trips)
h13.append(n_bikes)
h13.append(mean_overflows)
h13.append(iterations)
h13.append(mean_w_time)
h13.append(arrival_interval)
h13.append(distr_bikes)

hipos.append(h13)
# Plots: Redução pela metade -> Distribuições diferentes
# Plots: Redução pela metade -> Melhor distribuição reduzindo o intervalo de chegada
# Plots: Sem redução -> Reduzindo os intervalos de chegada
# Plots: Sem redução -> Reduzindo os intervalos de chegada com a melhor distribuição encontrada
# Comentario: aumento de capacidade em 10% com e sem redução
columns = ['mean_trips', 'n_bikes', 'mean_overflows', 'days', 'mean_wtime', 'arrival_interval', 'bike_distr']
hypoteses = [h1, h2, h3, h4, h8, h6, h7, h5, h9, h10, h11, h12, h13]
hps = pd.DataFrame(hypoteses, columns=columns)
hps
labels = ['Initial\n Half Bikes\n(50-30-20(%))', '40-40-20(%)', '30-50-20(%)', '30-40-30(%)', '2.0 min Interval\n + Best Distr', '1.5 min Interval + \nBest Distr', 'Initial + \n5% Capacity']
title = 'Média de viagens por dia para cada hipótese com metade da frota'

ax = hps[['mean_wtime', 'mean_overflows']].loc[hps.n_bikes < 3908].plot(kind='bar', figsize=(15, 10), rot=0, title=title)
for p in ax.patches:                 
    ax.annotate(round(p.get_height(),3), (p.get_x()+p.get_width()/2., p.get_height()), ha='center', va='center', xytext=(0, 10), textcoords='offset points')
ax.grid(linestyle=':', axis='y')
ax.set_xticklabels(labels)
labels = ['Initial', '1.5 min Interval', '1.0 min Interval', '1.0 min Interval + \nBest Distr', '1.0 min Interval + \nBest Distr + \n5% Capacity', '1.5 min Interval + \nBest Distr + \n5% Capacity']
title = 'Média de tempo de espera e overflow para cada hipótese sem frota reduzida'

ax = hps[['mean_trips', 'mean_overflows']].loc[hps.n_bikes == 3908].plot(kind='bar', figsize=(15, 10), rot=0, title=title)
ax.set_xticklabels(labels)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11649936377774636, 'fraction_numerical': 0.02559027286865545, 'mean_word_length': 4.672280139000267, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 55, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '618052', 'license': 'MIT', 'n_tokens_mistral': 8046, 'n_tokens_neox': 7492, 'n_words': 1718, 'path': 'notebooks/main.ipynb', 'repo_name': 'jvcanavarro/MetroBike', 'size': '21238'}