Find Depot-or-customer nodes shortest path.

Also cleanup input data in odiparpack/

Also rename sols/ to output/

1 files changed, 245 insertions, 1 deletions

diff --git a/example1.py b/example1.py
index 8b909bb..fdacd9d 100644
--- a/example1.py
+++ b/example1.py
@@ -1,6 +1,250 @@
 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'
@@ -17,4 +261,4 @@ def main():
 
 
 if __name__ == '__main__':
-    main()
+    lab()