osm3000/kebab_analysis.py

## kebab_analysis.py
"""
The core code for the kebab restaurant analysis project.

Notes:
- Given that I couldn't get the Google Maps python SDK to work, but the curl requests work perfectly (!!), this is a wrapper with "requests" to make the requests to the Google Maps API.
"""

from dotenv import load_dotenv

load_dotenv()
import requests
import json
import os
from time import sleep
import osmnx as ox
from geopandas import GeoDataFrame
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import json
import requests
import geopandas
from sklearn.model_selection import cross_validate
from sklearn.ensemble import RandomForestClassifier
from scipy.stats import pearsonr

# Constants
PROJECTION = 32631
DISTANCE_THRESHOLD = 1000
# DISTANCE_THRESHOLD = 750
# DISTANCE_THRESHOLD = 500

def get_stations_from_file(city: str):
    """
    Load the "liste-des-gares.geojson" file as a GeoDataFrame.
	You can download it from here:
	https://ressources.data.sncf.com/explore/dataset/liste-des-gares/export/?location=6,46.90879,1.85167&basemap=jawg.transports
    """
    stations = GeoDataFrame.from_file("liste-des-gares.geojson")
    city = city.upper()
    # Print the shape of the GeoDataFrame
    # Print all the headers, one by one
    # Make sure the city is in the "COMMUNE" column
    assert city in stations["commune"].str.upper().values
    # Filter the stations on the city
    stations = stations[stations["commune"].str.upper().str.contains(city)]
    # Filter on VOYAGEURS == "O"
    stations = stations[stations["voyageurs"] == "O"]
    # Filter on FRET == "N"
    stations = stations[stations["fret"] == "N"]
    print(f"Number of stations in {city}: {stations.shape[0]}")

    return stations


def get_data(text_query: str, lat: float, lon: float, radius: float = 500.0):
    """
    If the response has "nextPageToken" in it, then it will make another request to get the next page of results.
    """
    request_body = {
        "textQuery": text_query,
        "pageSize": 20,
        "locationBias": {
            "circle": {
                "center": {"latitude": lat, "longitude": lon},
                "radius": radius,
            }
        },
    }

    """
    Example of the request header:
    -H 'Content-Type: application/json' -H "X-Goog-Api-Key: ADD-KEY-HERE" \
    -H "X-Goog-FieldMask: places.id,places.displayName,places.formattedAddress,places.types,places.location"
    """
    request_headers = {
        "Content-Type": "application/json",
        "X-Goog-Api-Key": os.getenv("GMAPS_API_KEY"),
        "X-Goog-FieldMask": "places.id,places.displayName,places.formattedAddress,places.types,places.location,places.rating,places.reviews",
    }

    url = "https://places.googleapis.com/v1/places:searchText"

    # Make the request
    all_responses = []

    response = requests.post(url, headers=request_headers, json=request_body)
    response_data = json.loads(response.text)
    all_responses += response_data["places"]

    # While there is a nextPageToken, make another request
    while "nextPageToken" in response_data:
        request_body["pageToken"] = response_data["nextPageToken"]
        response = requests.post(url, headers=request_headers, json=request_body)
        response_data = json.loads(response.text)
        # all_responses.append(response_data['places'])
        all_responses += response_data["places"]
        sleep(2)
    sleep(2)

    return all_responses


def clean_gmap_results(data: list):
    """
    Given the data from the API, which is a list of dictionaries, we want to return a one dictionary, with "id" of each place as the key, and the rest of the data as the value.
    """
    clean_data = {}
    for item in data:
        id = item["id"]
        clean_data[id] = item
        # Remove the "id" key from the dictionary
        del clean_data[id]["id"]
        # Add the latitude and longitude keys
        clean_data[id]["latitude"] = clean_data[id]["location"]["latitude"]
        clean_data[id]["longitude"] = clean_data[id]["location"]["longitude"]
        del clean_data[id]["location"]
        # Add display name
        clean_data[id]["display_name"] = clean_data[id]["displayName"]["text"]
        del clean_data[id]["displayName"]
        # Convert the types (a list) into a string (to be able to store it in the database)
        clean_data[id]["types"] = ", ".join(clean_data[id]["types"])
        # Get the number of reviews, and the average rating. Then remove the reviews key.
        if "reviews" in clean_data[id]:
            reviews = clean_data[id]["reviews"]
            clean_data[id]["all_ratings"] = [
                float(review["rating"]) for review in reviews
            ]
            del clean_data[id]["reviews"]
        else:
            clean_data[id]["all_ratings"] = 0
    return clean_data

def get_all_restaurant_data(city: str = "Paris", fresh_data: bool = False):
    """
    Given the stations, we want to get the restaurant data for each station.
    """
    # Check if the file exists
    force_fresh_data = False
    if not os.path.exists(f"{city.lower()}_all_restaurant_data.json"):
        # raise FileNotFoundError(f"{city.lower()}_all_restaurant_data.json not found.")
        force_fresh_data = True
    all_restaurant_data = {}
    if fresh_data or force_fresh_data:
        stations = get_stations_from_file(city)
        all_restaurant_data = {}
        for i, row in stations.iterrows():
            lat, lon = row.geometry.y, row.geometry.x
            print(row["libelle"] + ":", lat, lon)
            data = get_data("kebab", lat, lon, radius=1000.0)
            clean_data = clean_gmap_results(data)
            all_restaurant_data.update(clean_data)

        # Save the data to a file
        with open(f"{city.lower()}_all_restaurant_data.json", "w") as f:
            json.dump(all_restaurant_data, f, indent=4)
    else:
        # Load the data from the file
        with open(f"{city.lower()}_all_restaurant_data.json", "r") as f:
            all_restaurant_data = json.load(f)

    return all_restaurant_data


def calculate_distances(city: str):
    """
    Load the stations from the file as GeoDataFrame, and load the restaurant data from the file as a dictionary.
    Hash the name of the station (libelle) to the station data.
    Convert the restaurant data to a GeoDataFrame.
    Calculate the distance between each station and each restaurant.
    """
    stations = get_stations_from_file(city)
    all_restaurant_data = get_all_restaurant_data(city=city, fresh_data=False)
    # Hash the name of the station to the station data
    station_hash = {}
    for i, row in stations.iterrows():
        station_hash[row["libelle"]] = row

    # Convert the restaurant data to a GeoDataFrame
    restaurant_data = []
    for item in all_restaurant_data:
        restaurant = all_restaurant_data[item]
        restaurant_data.append(
            {
                "name": restaurant.get("display_name", ""),
                "latitude": restaurant.get("latitude", 0),
                "longitude": restaurant.get("longitude", 0),
                "types": restaurant.get("types", ""),
                "rating": restaurant.get("rating", 0),
            }
        )
    restaurant_data = pd.DataFrame(restaurant_data)
    # Remove the restaurants that have a rating of 0
    restaurant_data = restaurant_data[restaurant_data["rating"] > 0]
    # Convert the restaurant data to a GeoDataFrame
    restaurant_data = GeoDataFrame(
        restaurant_data,
        geometry=geopandas.points_from_xy(
            restaurant_data.longitude, restaurant_data.latitude
        ),
        crs="EPSG:4326",
    )
    # Project the GeoDataFrame to the same CRS as the stations
    restaurant_data = restaurant_data.to_crs(PROJECTION)
    # Project the GeoDataFrame stations to the same CRS as the restaurant data
    stations = stations.to_crs(PROJECTION)
    # Calculate the euclidean distance between each station and each restaurant
    all_distances = []
    for i, single_station in stations.iterrows():
        # print(single_station["libelle"], single_station.geometry)
        for j, single_restaurant in restaurant_data.iterrows():
            # print(single_station["libelle"], single_restaurant["name"], single_restaurant.geometry)
            distance = single_station.geometry.distance(single_restaurant.geometry)
            # Calculate the angle between the station and the restaurant: do the math from the station to the restaurant
            angle = np.arctan2(
                single_restaurant.geometry.y - single_station.geometry.y,
                single_restaurant.geometry.x - single_station.geometry.x,
            )
            all_distances.append(
                {
                    "station": single_station["libelle"],
                    "restaurant": single_restaurant["name"],
                    "distance": distance,
                    "angle": angle,
                    "restaurant_rating": single_restaurant["rating"],
                    "train_station_latitude": single_station.geometry.y,
                    "train_station_longitude": single_station.geometry.x,
                    "restaurant_latitude": single_restaurant.geometry.y,
                    "restaurant_longitude": single_restaurant.geometry.x,
                }
            )
    all_distances = pd.DataFrame(all_distances)
    print(
        f"Number of restaurants X stations, before filtering: {all_distances.shape[0]}"
    )
    # Remove the restaurants that are more than 1000m away from the station
    all_distances = all_distances[all_distances["distance"] <= DISTANCE_THRESHOLD]
    print(
        f"Number of restaurants X stations, after filtering: {all_distances.shape[0]}"
    )
    # Get median of ratings, and create a new column for the ratings, thresholded by the median
    median_rating = all_distances["restaurant_rating"].median()
    print(f"Median rating: {median_rating}")
    all_distances["rating_threshold"] = all_distances["restaurant_rating"].apply(
        lambda x: 1 if x >= median_rating else 0
    )

    # Plot the ratings histogram
    # plt.hist(all_distances["restaurant_rating"], bins=20)
    # plt.show()
    # # Plot the distance histogram
    # plt.hist(all_distances["rating_threshold"], bins=20)
    # plt.show()
    # # Plot rating_threshold vs distance
    # plt.scatter(all_distances["distance"], all_distances["rating_threshold"])
    # plt.xlabel("Distance (m)")
    # plt.ylabel("Rating threshold")
    # plt.show()

    # Calculate the correlation between distance and restaurant_rating
    correlation = np.corrcoef(
        1 / all_distances["distance"], all_distances["restaurant_rating"]
    )
    # correlation = np.corrcoef(all_distances["distance"], all_distances["restaurant_rating"])
    correlation = np.round(correlation, 3)
    print(f"Correlation between distance and restaurant_rating: {correlation[0, 1]}")

    # Perform the Wilcoxon test
    # Calculate the Wilcoxon test
    try:
        stat, p = pearsonr(
            # 1 / all_distances["distance"], all_distances["restaurant_rating"]
            1 / all_distances["distance"],
            all_distances["rating_threshold"],
        )
        print(f"Statistics={stat}, p = {np.round(p, 4)}")
        # Interpret the results
        alpha = 0.05
        if p > alpha:
            print("Same distribution (fail to reject H0)")
        else:
            print("Different distribution (reject H0)")
    except Exception as e:
        print(f"Data is too small to apply the statistical test: {e}")

    return all_distances


def predictive_power_analysis(
    all_distances: pd.DataFrame, city: str, use_angle: bool = False
):
    """
    We went to find out if there is a non-linear relationship between the distance and the rating of the restaurant.
    We will use a non-linear model (random forest) to predict the rating of the restaurant based on the distance.
    Apply cross_validate to get the accuracy and f1 score.
    """
    # Split the data into training and testing data
    if use_angle:
        X = all_distances[["distance", "angle"]]
    else:
        X = all_distances[["distance"]]
    y = all_distances["rating_threshold"]

    # Count and print the number of 1s and 0s in the rating_threshold column
    print(
        f"Percentage of labels (ratings split based on the median rating): {y.value_counts(normalize=True).round(3)}"
    )
    # Train the model
    model = RandomForestClassifier()
    ignore_ml_model = False
    try:
        cross_val = cross_validate(model, X, y, cv=5, scoring=["accuracy", "f1"])
    except Exception as e:
        print(f"Data is too small to apply cross_validate: {e}")
        ignore_ml_model = True

    if not ignore_ml_model:
        print("Cross validation results:")
        for key in cross_val:
            if "time" in key:
                continue
            mean, std = np.mean(cross_val[key]), np.std(cross_val[key])
            # Round the mean and std to 3 decimal places
            mean, std = np.round(mean, 3), np.round(std, 3)
            print(key, mean, std)

        # Fit the model
        model.fit(X, y)
        plt.figure()

        if use_angle:
            ########################################################
            # Plot the decision boundary - In case of 2 features (distance and angle)
            ########################################################
            # Create a meshgrid
            x_min, x_max = X["distance"].min() - 1, X["distance"].max() + 1
            y_min, y_max = X["angle"].min() - 1, X["angle"].max() + 1
            xx, yy = np.meshgrid(
                np.arange(x_min, x_max, 0.1), np.arange(y_min, y_max, 0.1)
            )
            # Predict the values
            Z = model.predict(np.c_[xx.ravel(), yy.ravel()])
            Z = Z.reshape(xx.shape)
            # Plot the decision boundary
            plt.contourf(xx, yy, Z, alpha=0.4)
            plt.scatter(X["distance"], X["angle"], c=y, s=20, edgecolor="k")
            plt.xlabel("Distance (m)")
            plt.ylabel("Angle")
            # Colorbar
            plt.colorbar()
        else:
            ########################################################
            # Plot the decision boundary - In case of 1 feature (distance)
            ########################################################
            # Create a meshgrid
            x_min, x_max = X["distance"].min() - 1, X["distance"].max() + 1
            xx = np.linspace(x_min, x_max, 1000)
            # Predict the values
            Z = model.predict(xx.reshape(-1, 1))
            # Plot the decision boundary
            plt.plot(xx, Z, color="r")
            plt.scatter(X["distance"], y, s=20, edgecolor="k")
            plt.xlabel("Distance (m)")
            plt.ylabel("Rating threshold")

        plt.title(
            f"Decision boundary - {city} - Model F1 score: {np.round(cross_val['test_f1'].mean(), 3)} - Distance threshold: {DISTANCE_THRESHOLD}m"
        )
        # plt.show()
        plt.savefig(
            f"{city}_decision_boundary_angle:{use_angle}_{DISTANCE_THRESHOLD}m.png"
        )


if __name__ == "__main__":
    data = get_all_restaurant_data(city="Paris", fresh_data=True)
    # print(data)
    for city in [
        "Paris",
        "Toulouse",
        "Lyon",
        "Marseille",
        "Lille",
        "Bordeaux",
        "Nantes",
        "Strasbourg",
        "Rennes",
        "Montpellier",
    ]:
        print(f"Calculating distances for {city}...")
        all_distances = calculate_distances(city)
        print("*" * 50)
        print("Let's discover non-linear relationships...")
        print("*" * 50)
        for use_angle in [False, True]:
            print(f"Using angle: {use_angle}")
            predictive_power_analysis(all_distances, city, use_angle=use_angle)
        print("-" * 100)
        print("\n")
        # break
	"""
	The core code for the kebab restaurant analysis project.

	Notes:
	- Given that I couldn't get the Google Maps python SDK to work, but the curl requests work perfectly (!!), this is a wrapper with "requests" to make the requests to the Google Maps API.
	"""

	from dotenv import load_dotenv

	load_dotenv()
	import requests
	import json
	import os
	from time import sleep
	import osmnx as ox
	from geopandas import GeoDataFrame
	import pandas as pd
	import numpy as np
	import matplotlib.pyplot as plt
	import json
	import requests
	import geopandas
	from sklearn.model_selection import cross_validate
	from sklearn.ensemble import RandomForestClassifier
	from scipy.stats import pearsonr

	# Constants
	PROJECTION = 32631
	DISTANCE_THRESHOLD = 1000
	# DISTANCE_THRESHOLD = 750
	# DISTANCE_THRESHOLD = 500

	def get_stations_from_file(city: str):
	"""
	Load the "liste-des-gares.geojson" file as a GeoDataFrame.
	You can download it from here:
	https://ressources.data.sncf.com/explore/dataset/liste-des-gares/export/?location=6,46.90879,1.85167&basemap=jawg.transports
	"""
	stations = GeoDataFrame.from_file("liste-des-gares.geojson")
	city = city.upper()
	# Print the shape of the GeoDataFrame
	# Print all the headers, one by one
	# Make sure the city is in the "COMMUNE" column
	assert city in stations["commune"].str.upper().values
	# Filter the stations on the city
	stations = stations[stations["commune"].str.upper().str.contains(city)]
	# Filter on VOYAGEURS == "O"
	stations = stations[stations["voyageurs"] == "O"]
	# Filter on FRET == "N"
	stations = stations[stations["fret"] == "N"]
	print(f"Number of stations in {city}: {stations.shape[0]}")

	return stations


	def get_data(text_query: str, lat: float, lon: float, radius: float = 500.0):
	"""
	If the response has "nextPageToken" in it, then it will make another request to get the next page of results.
	"""
	request_body = {
	"textQuery": text_query,
	"pageSize": 20,
	"locationBias": {
	"circle": {
	"center": {"latitude": lat, "longitude": lon},
	"radius": radius,
	}
	},
	}

	"""
	Example of the request header:
	-H 'Content-Type: application/json' -H "X-Goog-Api-Key: ADD-KEY-HERE" \
	-H "X-Goog-FieldMask: places.id,places.displayName,places.formattedAddress,places.types,places.location"
	"""
	request_headers = {
	"Content-Type": "application/json",
	"X-Goog-Api-Key": os.getenv("GMAPS_API_KEY"),
	"X-Goog-FieldMask": "places.id,places.displayName,places.formattedAddress,places.types,places.location,places.rating,places.reviews",
	}

	url = "https://places.googleapis.com/v1/places:searchText"

	# Make the request
	all_responses = []

	response = requests.post(url, headers=request_headers, json=request_body)
	response_data = json.loads(response.text)
	all_responses += response_data["places"]

	# While there is a nextPageToken, make another request
	while "nextPageToken" in response_data:
	request_body["pageToken"] = response_data["nextPageToken"]
	response = requests.post(url, headers=request_headers, json=request_body)
	response_data = json.loads(response.text)
	# all_responses.append(response_data['places'])
	all_responses += response_data["places"]
	sleep(2)
	sleep(2)

	return all_responses


	def clean_gmap_results(data: list):
	"""
	Given the data from the API, which is a list of dictionaries, we want to return a one dictionary, with "id" of each place as the key, and the rest of the data as the value.
	"""
	clean_data = {}
	for item in data:
	id = item["id"]
	clean_data[id] = item
	# Remove the "id" key from the dictionary
	del clean_data[id]["id"]
	# Add the latitude and longitude keys
	clean_data[id]["latitude"] = clean_data[id]["location"]["latitude"]
	clean_data[id]["longitude"] = clean_data[id]["location"]["longitude"]
	del clean_data[id]["location"]
	# Add display name
	clean_data[id]["display_name"] = clean_data[id]["displayName"]["text"]
	del clean_data[id]["displayName"]
	# Convert the types (a list) into a string (to be able to store it in the database)
	clean_data[id]["types"] = ", ".join(clean_data[id]["types"])
	# Get the number of reviews, and the average rating. Then remove the reviews key.
	if "reviews" in clean_data[id]:
	reviews = clean_data[id]["reviews"]
	clean_data[id]["all_ratings"] = [
	float(review["rating"]) for review in reviews
	]
	del clean_data[id]["reviews"]
	else:
	clean_data[id]["all_ratings"] = 0
	return clean_data

	def get_all_restaurant_data(city: str = "Paris", fresh_data: bool = False):
	"""
	Given the stations, we want to get the restaurant data for each station.
	"""
	# Check if the file exists
	force_fresh_data = False
	if not os.path.exists(f"{city.lower()}_all_restaurant_data.json"):
	# raise FileNotFoundError(f"{city.lower()}_all_restaurant_data.json not found.")
	force_fresh_data = True
	all_restaurant_data = {}
	if fresh_data or force_fresh_data:
	stations = get_stations_from_file(city)
	all_restaurant_data = {}
	for i, row in stations.iterrows():
	lat, lon = row.geometry.y, row.geometry.x
	print(row["libelle"] + ":", lat, lon)
	data = get_data("kebab", lat, lon, radius=1000.0)
	clean_data = clean_gmap_results(data)
	all_restaurant_data.update(clean_data)

	# Save the data to a file
	with open(f"{city.lower()}_all_restaurant_data.json", "w") as f:
	json.dump(all_restaurant_data, f, indent=4)
	else:
	# Load the data from the file
	with open(f"{city.lower()}_all_restaurant_data.json", "r") as f:
	all_restaurant_data = json.load(f)

	return all_restaurant_data


	def calculate_distances(city: str):
	"""
	Load the stations from the file as GeoDataFrame, and load the restaurant data from the file as a dictionary.
	Hash the name of the station (libelle) to the station data.
	Convert the restaurant data to a GeoDataFrame.
	Calculate the distance between each station and each restaurant.
	"""
	stations = get_stations_from_file(city)
	all_restaurant_data = get_all_restaurant_data(city=city, fresh_data=False)
	# Hash the name of the station to the station data
	station_hash = {}
	for i, row in stations.iterrows():
	station_hash[row["libelle"]] = row

	# Convert the restaurant data to a GeoDataFrame
	restaurant_data = []
	for item in all_restaurant_data:
	restaurant = all_restaurant_data[item]
	restaurant_data.append(
	{
	"name": restaurant.get("display_name", ""),
	"latitude": restaurant.get("latitude", 0),
	"longitude": restaurant.get("longitude", 0),
	"types": restaurant.get("types", ""),
	"rating": restaurant.get("rating", 0),
	}
	)
	restaurant_data = pd.DataFrame(restaurant_data)
	# Remove the restaurants that have a rating of 0
	restaurant_data = restaurant_data[restaurant_data["rating"] > 0]
	# Convert the restaurant data to a GeoDataFrame
	restaurant_data = GeoDataFrame(
	restaurant_data,
	geometry=geopandas.points_from_xy(
	restaurant_data.longitude, restaurant_data.latitude
	),
	crs="EPSG:4326",
	)
	# Project the GeoDataFrame to the same CRS as the stations
	restaurant_data = restaurant_data.to_crs(PROJECTION)
	# Project the GeoDataFrame stations to the same CRS as the restaurant data
	stations = stations.to_crs(PROJECTION)
	# Calculate the euclidean distance between each station and each restaurant
	all_distances = []
	for i, single_station in stations.iterrows():
	# print(single_station["libelle"], single_station.geometry)
	for j, single_restaurant in restaurant_data.iterrows():
	# print(single_station["libelle"], single_restaurant["name"], single_restaurant.geometry)
	distance = single_station.geometry.distance(single_restaurant.geometry)
	# Calculate the angle between the station and the restaurant: do the math from the station to the restaurant
	angle = np.arctan2(
	single_restaurant.geometry.y - single_station.geometry.y,
	single_restaurant.geometry.x - single_station.geometry.x,
	)
	all_distances.append(
	{
	"station": single_station["libelle"],
	"restaurant": single_restaurant["name"],
	"distance": distance,
	"angle": angle,
	"restaurant_rating": single_restaurant["rating"],
	"train_station_latitude": single_station.geometry.y,
	"train_station_longitude": single_station.geometry.x,
	"restaurant_latitude": single_restaurant.geometry.y,
	"restaurant_longitude": single_restaurant.geometry.x,
	}
	)
	all_distances = pd.DataFrame(all_distances)
	print(
	f"Number of restaurants X stations, before filtering: {all_distances.shape[0]}"
	)
	# Remove the restaurants that are more than 1000m away from the station
	all_distances = all_distances[all_distances["distance"] <= DISTANCE_THRESHOLD]
	print(
	f"Number of restaurants X stations, after filtering: {all_distances.shape[0]}"
	)
	# Get median of ratings, and create a new column for the ratings, thresholded by the median
	median_rating = all_distances["restaurant_rating"].median()
	print(f"Median rating: {median_rating}")
	all_distances["rating_threshold"] = all_distances["restaurant_rating"].apply(
	lambda x: 1 if x >= median_rating else 0
	)

	# Plot the ratings histogram
	# plt.hist(all_distances["restaurant_rating"], bins=20)
	# plt.show()
	# # Plot the distance histogram
	# plt.hist(all_distances["rating_threshold"], bins=20)
	# plt.show()
	# # Plot rating_threshold vs distance
	# plt.scatter(all_distances["distance"], all_distances["rating_threshold"])
	# plt.xlabel("Distance (m)")
	# plt.ylabel("Rating threshold")
	# plt.show()

	# Calculate the correlation between distance and restaurant_rating
	correlation = np.corrcoef(
	1 / all_distances["distance"], all_distances["restaurant_rating"]
	)
	# correlation = np.corrcoef(all_distances["distance"], all_distances["restaurant_rating"])
	correlation = np.round(correlation, 3)
	print(f"Correlation between distance and restaurant_rating: {correlation[0, 1]}")

	# Perform the Wilcoxon test
	# Calculate the Wilcoxon test
	try:
	stat, p = pearsonr(
	# 1 / all_distances["distance"], all_distances["restaurant_rating"]
	1 / all_distances["distance"],
	all_distances["rating_threshold"],
	)
	print(f"Statistics={stat}, p = {np.round(p, 4)}")
	# Interpret the results
	alpha = 0.05
	if p > alpha:
	print("Same distribution (fail to reject H0)")
	else:
	print("Different distribution (reject H0)")
	except Exception as e:
	print(f"Data is too small to apply the statistical test: {e}")

	return all_distances


	def predictive_power_analysis(
	all_distances: pd.DataFrame, city: str, use_angle: bool = False
	):
	"""
	We went to find out if there is a non-linear relationship between the distance and the rating of the restaurant.
	We will use a non-linear model (random forest) to predict the rating of the restaurant based on the distance.
	Apply cross_validate to get the accuracy and f1 score.
	"""
	# Split the data into training and testing data
	if use_angle:
	X = all_distances[["distance", "angle"]]
	else:
	X = all_distances[["distance"]]
	y = all_distances["rating_threshold"]

	# Count and print the number of 1s and 0s in the rating_threshold column
	print(
	f"Percentage of labels (ratings split based on the median rating): {y.value_counts(normalize=True).round(3)}"
	)
	# Train the model
	model = RandomForestClassifier()
	ignore_ml_model = False
	try:
	cross_val = cross_validate(model, X, y, cv=5, scoring=["accuracy", "f1"])
	except Exception as e:
	print(f"Data is too small to apply cross_validate: {e}")
	ignore_ml_model = True

	if not ignore_ml_model:
	print("Cross validation results:")
	for key in cross_val:
	if "time" in key:
	continue
	mean, std = np.mean(cross_val[key]), np.std(cross_val[key])
	# Round the mean and std to 3 decimal places
	mean, std = np.round(mean, 3), np.round(std, 3)
	print(key, mean, std)

	# Fit the model
	model.fit(X, y)
	plt.figure()

	if use_angle:
	########################################################
	# Plot the decision boundary - In case of 2 features (distance and angle)
	########################################################
	# Create a meshgrid
	x_min, x_max = X["distance"].min() - 1, X["distance"].max() + 1
	y_min, y_max = X["angle"].min() - 1, X["angle"].max() + 1
	xx, yy = np.meshgrid(
	np.arange(x_min, x_max, 0.1), np.arange(y_min, y_max, 0.1)
	)
	# Predict the values
	Z = model.predict(np.c_[xx.ravel(), yy.ravel()])
	Z = Z.reshape(xx.shape)
	# Plot the decision boundary
	plt.contourf(xx, yy, Z, alpha=0.4)
	plt.scatter(X["distance"], X["angle"], c=y, s=20, edgecolor="k")
	plt.xlabel("Distance (m)")
	plt.ylabel("Angle")
	# Colorbar
	plt.colorbar()
	else:
	########################################################
	# Plot the decision boundary - In case of 1 feature (distance)
	########################################################
	# Create a meshgrid
	x_min, x_max = X["distance"].min() - 1, X["distance"].max() + 1
	xx = np.linspace(x_min, x_max, 1000)
	# Predict the values
	Z = model.predict(xx.reshape(-1, 1))
	# Plot the decision boundary
	plt.plot(xx, Z, color="r")
	plt.scatter(X["distance"], y, s=20, edgecolor="k")
	plt.xlabel("Distance (m)")
	plt.ylabel("Rating threshold")

	plt.title(
	f"Decision boundary - {city} - Model F1 score: {np.round(cross_val['test_f1'].mean(), 3)} - Distance threshold: {DISTANCE_THRESHOLD}m"
	)
	# plt.show()
	plt.savefig(
	f"{city}_decision_boundary_angle:{use_angle}_{DISTANCE_THRESHOLD}m.png"
	)


	if __name__ == "__main__":
	data = get_all_restaurant_data(city="Paris", fresh_data=True)
	# print(data)
	for city in [
	"Paris",
	"Toulouse",
	"Lyon",
	"Marseille",
	"Lille",
	"Bordeaux",
	"Nantes",
	"Strasbourg",
	"Rennes",
	"Montpellier",
	]:
	print(f"Calculating distances for {city}...")
	all_distances = calculate_distances(city)
	print("" 50)
	print("Let's discover non-linear relationships...")
	print("" 50)
	for use_angle in [False, True]:
	print(f"Using angle: {use_angle}")
	predictive_power_analysis(all_distances, city, use_angle=use_angle)
	print("-" * 100)
	print("\n")
	# break
No results found