run.py

from multiprocessing import Pool
import math
import os
import random
import time
import numpy as np
import pandas as pd
from intervaltree import Interval, IntervalTree
import hashlib


def run_sim(comb):
    for i in range (30):
        # Create directory for logs
        condition = "sensorcount" + "-" + str(comb["sensor_count"]) + "_" + "channel" + "-" + str(comb["channel"]) + "_" + "rwdelay" + "-" +  str(comb["rw_delay"]) + "_" + "rwduration" + "-" +  str(comb["rw_duration"])
        directory = os.path.join(".", "logs", comb["param"].replace("_", ""), condition, hashlib.sha256(
            (str(int(time.time())) + str(random.randint(1000, 9999))).encode()).hexdigest())

        if not os.path.exists(directory):
            os.makedirs(directory)

        dataframes = []
        # Create gateway and sensors (dicts)  -> 0: idx, 1: x-coordinate, 2: y-coordinate, 3: spreading factor, 4: collision list
        gateway, sensors = create_sensors(comb)
        sensors = update_sf_collisions(sensors)

        # log_sensor_collisions(sensors, directory)

        ############################################ S I M U L A T I O N ############################################
        # General parameters
        sim_hours = 168  # 1 week in hours
        gateway_payload = 1  # bytes
        receive_window_delay = comb["rw_delay"]  # ms
        ack_delay = receive_window_delay  # ms
        receive_window_time = comb["rw_duration"]  # ms

        overlapping_packets = []  # Note down which packets overlap with the next hour to include them later

        for current_hour in range(sim_hours):
            hour_interval = (current_hour * 3600000, (current_hour + 1) * 3600000)

            # Append old overlapping packets and delete them
            packet_entries = overlapping_packets
            overlapping_packets = []

            # Generate 1 packet per sensor to gateway
            for idx in sensors.keys():
                sf = sensors[idx][3]

                # Generate time for packet_start, payload, and channel
                payload = np.random.randint(1,52)

                channel_low = 2
                channel_upper = channel_low + comb["channel"]
                channel = np.random.randint(channel_low, channel_upper)
                # channel = 2  # for testing
                packet_start = random.randint(hour_interval[0], hour_interval[1])
                packet_end = packet_start + payload_size_to_time(payload, sf)
                # Create a new packet entry for sensor without ACKs
                new_entry = {
                    'gateway_idx': 0,
                    'sensor_idx': idx,
                    'channel': channel,
                    'sf': sf,
                    'payload': payload,
                    'toa': payload_size_to_time(payload, sf),
                    'packet_start': packet_start,
                    'packet_end': packet_end,
                    'ack_delay': float('inf'),
                    'ack_start': float('inf'),
                    'ack_end': float('inf'),
                    'ack_received': float('inf')
                }
                packet_entries.append(new_entry)

                if packet_end >= (current_hour + 1) * 3600000:
                    overlapping_packets.append(new_entry)

            # Convert the list of entries to packet_df and append all at once
            df = pd.DataFrame(packet_entries)

            # Check for collisions per channel
            df = detect_collisions(df, sensors)

            # Create ACKs of packets which end in this hour
            df = create_acks(df, gateway_payload, ack_delay, receive_window_delay, receive_window_time, hour_interval )

            # Drop overlapping packets with inf
            df = df[df['ack_received'] != float('inf')]
            dataframes.append(df)

        # Save dfs in one big one
        big_dataframe = pd.concat(dataframes, ignore_index=True)
        filename = os.path.join(directory, f"transmission_df.parquet.gzip")
        big_dataframe.to_parquet(filename, compression='gzip')

        print(f"COMPLETED RUN {i+1}  |  Condition: {condition}")

def create_sensors(comb):
    # np.random.seed(42)
    taken_coordinates = []  # List of coordinates that are taken so sensors don't land on each other
    nr_sensors = int(comb["sensor_count"])

    # Create gateway first
    gateway_idx = 0
    gateway_x = 2171
    gateway_y = 2171
    gateway = [gateway_idx, gateway_x, gateway_y]
    taken_coordinates.append((gateway_x, gateway_y))

    # Sensor mapping with ID -> sensor_object
    sensors = dict()
    for idx in range(nr_sensors):
        sensor_idx = idx + 1
        while True:
            sensor_x = np.random.randint(gateway_x * 2)
            sensor_y = np.random.randint(gateway_y * 2)
            distance = np.sqrt((sensor_x - gateway_x) ** 2 + (sensor_y - gateway_y) ** 2)

            # Check if coordinates within circular range of gateway
            if distance <= gateway_x and (sensor_x, sensor_y) not in taken_coordinates:
                break

        sf = update_sf(gateway, sensor_x, sensor_y)
        sf_collisions = []
        sensors[sensor_idx] = [sensor_idx, sensor_x, sensor_y, sf, sf_collisions]
    return gateway, sensors

def update_sf(gateway, sensor_x, sensor_y):
    # dist is distance from sensor to gateway
    gateway_x = gateway[1]
    gateway_y = gateway[2]
    dist = math.sqrt((gateway_x - sensor_x) ** 2 + (gateway_y - sensor_y) ** 2)

    # Go for lowest sf possible
    currsf = 12
    for new_sf in range(12, 6, -1):
        hata_dist = hata(new_sf, 15, 1)
        if dist < hata_dist:
            currsf = new_sf
    return currsf

def hata(sf, gw_height, sens_height):
    # Returns distance in m
    # Represents Spreading Factor tolerances for Path Loss in dB - Credits to Noah's sim
    plrange = [131, 134, 137, 140, 141, 144]
    distance = 10 ** (-(69.55 + 76.872985 - 13.82 * math.log10(gw_height) - 3.2 * (math.log10(
        11.75 * sens_height) ** 2) + 4.97 - plrange[sf - 7]) / (44.9 - 6.55 * math.log10(gw_height)))
    return distance * 1000

def update_sf_collisions(sensors):
    updated_sensors = dict()
    # Precompute maximum sending distances for each sensor
    max_distances = {}
    for sensor_idx, sensor in sensors.items():
        sensor_sf = sensor[3]
        max_distances[sensor_idx] = max(hata(sensor_sf, 15, 1), 0)

    # Iterate over the sensors to detect collisions
    for sensor_idx, sensor in sensors.items():
        sensor_x = sensor[1]
        sensor_y = sensor[2]
        sensor_sf = sensor[3]
        sensor_sf_collisions = set()  # Avoid duplicate collisions

        # Check if this sensor collides with every other sensor
        for other_idx, other_sensor in sensors.items():
            if other_idx != sensor_idx:  # Avoid comparing sensor with itself
                other_x = other_sensor[1]
                other_y = other_sensor[2]

                # Get maximum sending distance of both sensors
                distance_sensor = max_distances[sensor_idx]
                distance_other_sensor = max_distances[other_idx]

                # Calculate Euclidean distance between the two sensors
                distance = math.sqrt((other_x - sensor_x) ** 2 + (other_y - sensor_y) ** 2)
                if distance <= max(distance_sensor, distance_other_sensor):
                    sensor_sf_collisions.add(other_idx)  # Add collision index to the set

        # Store the updated sensor with its collision list as a set
        updated_sensors[sensor_idx] = [sensor_idx, sensor_x, sensor_y, sensor_sf, list(sensor_sf_collisions)]
    return updated_sensors


def log_sensor_collisions(sensors, directory):
    # Log collisions list for each sensor to calculate collision probabilities
    sensor_entries = []

    for sensor in sensors.keys():
        # 0: idx, 1: x - coord, 2: y - coord, 3: sf, 4: sf_collisions
        new_entry = {
            'sensor_idx': sensor,
            'collisions': sensors[sensor][4]
        }
        sensor_entries.append(new_entry)

    sensor_df = pd.DataFrame(sensor_entries)
    filename = os.path.join(directory, f"sensor_df.parquet.gzip")
    sensor_df.to_parquet(filename, compression='gzip')


def payload_size_to_time(payload, sf):
    # Calculate time-on-air based on payload of packet and sf of sensor
    bw = 125
    pl = payload + 3
    cr = 4
    crc = 1
    h = 0
    de = 0
    sf = sf
    npreamble = 8

    rs = bw / (math.pow(2, sf))
    ts = 1 / rs
    symbol = 8 + max(math.ceil((8.0 * pl - 4.0 * sf + 28 + 16 * crc - 20.0 * h) /
                               (4.0 * (sf - 2.0 * de))) * (cr + 4), 0)
    tpreamble = (npreamble + 4.25) * ts
    tpayload = symbol * ts
    toa = tpreamble + tpayload
    return toa

def detect_collisions(packet_df, sensors):
    # Group by channel and sort by packet start
    for channel, group in packet_df.groupby('channel'):
        group = group.sort_values(by='packet_start')

        # Store the intervals of packets for the current channel
        tree = IntervalTree()

        # Iterate over packets in the current channel group
        for i, row_i in group.iterrows():
            packet_start = row_i['packet_start']
            packet_end = row_i['packet_end']
            sensor_idx = row_i['sensor_idx']

            # Create an interval from packet_start to packet_end & query for overlaps
            interval_i = Interval(packet_start, packet_end,
                                  data={'sensor_idx': sensor_idx, 'packet_end': packet_end, 'row_idx': i})
            overlapping_intervals = tree[packet_start:packet_end]

            # Check if overlapping intervals are colliding
            for interval in overlapping_intervals:
                row_j = interval.data
                if row_j['sensor_idx'] != sensor_idx:  # Skip comparing the same packet
                    # Check collision list
                    sensor_i_collisions = sensors[row_i['sensor_idx']][4]
                    sensor_j_collisions = sensors[row_j['sensor_idx']][4]

                    # Ensure the collision lists are of the correct type list/set
                    if isinstance(sensor_i_collisions, (list, set)) and isinstance(sensor_j_collisions,
                                                                                   (list, set)):
                        # Check if overlap in the collision list
                        if row_j['sensor_idx'] in sensor_i_collisions or row_i['sensor_idx'] in sensor_j_collisions:
                            # Mark both rows as colliding
                            packet_df.at[row_i.name, 'ack_received'] = 0
                            packet_df.at[row_j['row_idx'], 'ack_received'] = 0

            # Add the current packet's interval to the tree
            tree.add(interval_i)
    return packet_df

def create_acks(df, gateway_payload, ack_delay, receive_window_delay, receive_window_time, hour_interval):
    # Update ACKs
    df = df.sort_values(by='packet_end')
    max_ack_end = 0  # End of last ACK that was sent
    for i, row in df.iterrows():
        if row['ack_received'] != 0 and  hour_interval[0] <= row['packet_end'] < hour_interval[1]:
            ack_start = max(row['packet_end'] + ack_delay, max_ack_end)
            current_ack_delay = ack_start - row['packet_end']  # Standard ACK delay included
            # ack_end = ack_start + payload_size_to_time(gateway_payload, row['sf'])  # SF of ACK is senspr's SF
            ack_end = ack_start + payload_size_to_time(gateway_payload, 12) # SF of ACK is always 12
            df.at[i, 'ack_start'] = ack_start
            df.at[i, 'ack_delay'] = current_ack_delay
            df.at[i, 'ack_end'] = ack_end

            ack_received = -1
            if (row['packet_end'] + receive_window_delay) <= ack_start <= (
                    row['packet_end'] + receive_window_delay + receive_window_time):
                ack_received = 1

            df.at[i, 'ack_received'] = ack_received
            max_ack_end = max(max_ack_end, ack_end)

    return df

if __name__ == "__main__":
    # Tested parameters
    params = {
        "sensor_counts": [1500, 3000, 4500, 6000, 7500, 9000, 10500, 12000],
        "rw_delay": [0, 1000, 2000, 3000],
        "channels": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15],
        "rw_duration": [1000, 2000, 3000]
    }

    combinations = []
    for sensor_count in params["sensor_counts"]:
        combinations.append({"param": "sensor_count","sensor_count": sensor_count, "ack_delay": 1000, "rw_delay": 1000, "rw_duration": 1000, "channel": 15})

    for channel in params["channels"]:
        combinations.append({"param": "channel", "sensor_count": 7500, "ack_delay": 1000, "rw_delay": 1000,
                             "rw_duration": 1000, "channel": channel})

    for rw_delay in params["rw_delay"]:
        combinations.append({"param": "rw_delay", "sensor_count": 7500, "ack_delay": rw_delay, "rw_delay": rw_delay, "rw_duration": 1000, "channel": 15})

    for rw_duration in params["rw_duration"]:
        combinations.append({"param": "rw_duration", "sensor_count": 7500, "ack_delay": 1000, "rw_delay": 1000, "rw_duration": rw_duration, "channel": 15})


    with Pool(8) as p:
        results = p.map(run_sim, combinations)