path: root/example1.py

from vrptw_base import VrptwGraph
from multiple_ant_colony_system import MultipleAntColonySystem

import numpy as np
import csv
import math
import random
import networkx as nx
import matplotlib.pyplot as plt


class Node:
    """
    Attributes
    ----------
    index : int
        index
    ubigeo : str
        6 digits
    lat : float
        Latitud (angulo de los Paralelos) en grados sexagesimales
    lon : float
        Longitud (angulo de los Meridianos) en grados sexagesimales
    is_depot : bool
    demand : int
    ready_time : float
    due_time : float
    service_time : float
    x : float
        Aproximacion local, en Km respecto
    y : float
        Aproximacion local, en Km respecto ("Lima")

    Notes
    -----
    Web Mercator projection (BAD):

        x = floor(256 / (2 * math.pi) * 2**(zoom_level) * (lon + math.pi))
        y = floor(265 / (2 * math.pi) * 2**(zoom_level) * (math.pi - math.ln(math.tan( math.pi / 4 + lat / 2 ))))

        x = R * lon
        y = R * ln(tan(pi/4 + lat/2)

    Both `lon` and `lat` in radians.

    "Lima": -12.04591952,-77.03049615  (lat, long)
    """

    def __init__(self, index: int, ubigeo, lat, lon, is_depot,
                 demand, ready_time, due_time, service_time):
        super()
        self.index = index
        self.ubigeo = ubigeo

        if is_depot:
            self.is_depot = True
        else:
            self.is_depot = False

        earth_radius_km = 6371  # Avg. radius
        lima_lat = -12.04591952
        lima_lon = -77.03049615

        self.lat = lat
        self.lon = lon
        self.x = (lon - lima_lon) * (math.pi / 180) * earth_radius_km
        self.y = (lat - lima_lat) * (math.pi / 180) * earth_radius_km
        self.x = round(self.x, 3)
        self.y = round(self.y, 3)
        self.demand = demand
        self.ready_time = 0  # ready_time
        self.due_time = due_time
        self.service_time = service_time


def _deg2rad(degrees):
    return degrees * math.pi / 180


def distance_between_coordinates(lat1, lon1, lat2, lon2):
    """Returns the distance in km"""
    earth_radius_km = 6371        # Avg. radius

    lat1 = _deg2rad(lat1)
    lat2 = _deg2rad(lat2)

    d_lat = lat2 - lat1
    d_lon = lon2 - lon1

    a = math.sin(d_lat/2) * math.sin(d_lat/2) \
        + math.sin(d_lon/2) * math.sin(d_lon/2) * math.cos(lat1) * math.cos(lat2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    return earth_radius_km * c


def _read_nodes(depot_num):
    """
    Read nodes from input file.

    Attributes
    ----------
    depot_num : int
        Number of depots

    Returns
    -------
    list
        List of Nodes.  The first `depot_num` nodes are the depots.
    """
    depots = []
    no_depots = []
    with open('odiparpack/inf226.oficinas_mod.csv', newline='') as csvfile:
        orders = csv.reader(csvfile)
        count = 1
        index_customer = depot_num
        index_depot = 0
        for row in orders:
            if count == 1:
                count += 1
                continue

            ubigeo, dept, prov, lat, lon, region_natural, is_depot = row[:7]
            demand, ready_time, due_time, service_time = row[7:]

            n = Node(
                -1, ubigeo, float(lat), float(lon),
                int(is_depot),
                demand=float(demand), ready_time=0,
                due_time=float(due_time), service_time=60
            )

            if n.is_depot:
                n.index = index_depot
                depots.append(n)
                index_depot += 1
            else:
                n.index = index_customer
                no_depots.append(n)
                index_customer += 1

            count += 1

    return depots + no_depots


def _read_tramos():
    """
    Lee archivo de tramos

    Returns
    -------
    list
        Lista de tuplas con los tramos  (ubigeo1, ubigeo2)
    """
    tramos = []
    with open('odiparpack/inf226.tramos.v.2.0.csv', newline='') as f:
        rows = csv.reader(f)
        count = 1
        for row in rows:
            if count >= 2:
                tramos.append([row[0], row[1]])
            count += 1
    return tramos


def make_complete_customer_node_graph(nodes, tramos):
    """
    Nodes + Tramos make a network.  An undirected graph.  But we only care
    about the nodes that are either have demand (aka. customers), or are depots.

    This function creates a sub-graph of only depot-or-customer nodes, and
    finds the nearest path + path_lenght between them (on the original network).

    Parameters
    ----------
    nodes : list
        List of Node objects
    tramos : list
        List of tuples of "ubigeo"

    Returns
    -------
    list
        List of depot-or-customer nodes,
        Shortest distance matrix (np.narray),
        Dict with shortest distance and path object for all pairs of nodes
    """
    elist = []

    # Load sparse graph
    nodes_d = {n.ubigeo: n for n in nodes}
    for ubigeo1, ubigeo2 in tramos:
        n1 = nodes_d[ubigeo1]
        n2 = nodes_d[ubigeo2]
        dist = round(distance_between_coordinates(n1.lat, n1.lon,
                                                  n2.lat, n2.lon), 0)
        elist.append((ubigeo1, ubigeo2, dist))

    G = nx.Graph()
    G.add_weighted_edges_from(elist)

    # Calculate missing edges using the shortest path
    # (Thus making the graph complete)
    len_path = dict(nx.all_pairs_dijkstra(G, weight='weight'))
    # for node, (dist, path) in len_path:
    #    print(f"{node}: {dist} [{path}]")

    # Prune all nodes that have no demand and "reset index" of nodes
    depot_customer_nodes = []
    i = 0
    for n in nodes:
        if n.demand > 0 or n.is_depot:
            n.index = i
            depot_customer_nodes.append(n)
            i += 1
    node_num = i

    # Create distance matrix
    node_dist_mat = np.zeros((node_num, node_num))
    for n1 in depot_customer_nodes:
        for n2 in depot_customer_nodes:
            node_dist_mat[n1.index, n2.index] = \
                len_path[n1.ubigeo][0][n2.ubigeo]

    return depot_customer_nodes, node_dist_mat, len_path


def lab():
    depot_num = 3
    nodes = _read_nodes(depot_num)
    tramos = _read_tramos()
    nodes, node_dist_mat, len_path = \
        make_complete_customer_node_graph(nodes, tramos)

    # nodes with demand
    for n in nodes[:8]:
        print(n.__dict__)

    # shortest distance between nodes in network
    print(node_dist_mat[:8, :8])

    # shortest path
    for n1 in nodes[:1]:
        for n2 in nodes[:8]:
            shortest_path = len_path[n1.ubigeo][1][n2.ubigeo]
            print(f"{n1.ubigeo} => {n2.ubigeo}: [{shortest_path}]")


def main():
    #file_path = './solomon-100/c101.txt'
    file_path = './odiparpack/pedidosperu195.txt'
    ants_num = 10
    beta = 2  # 5
    q0 = 0.1  # 0.5 ?
    show_figure = True

    graph = VrptwGraph(file_path)
    macs = MultipleAntColonySystem(graph, ants_num=ants_num, beta=beta, q0=q0,
                                   whether_or_not_to_show_figure=show_figure)
    macs.run_multiple_ant_colony_system()


if __name__ == '__main__':
    lab()