graph.py

# coding: utf-8
# pylint: disable=missing-docstring, invalid-name, attribute-defined-outside-init, no-member

import math
import os
import pickle
from decimal import Decimal
import numpy as np
from matplotlib import pyplot as plt

from log import log
from worker import SharedMatMain

# get mantissa & exponent of a number
def getManExp(number):
    d = Decimal(number)
    _, digits, exponent = d.as_tuple()
    exp = len(digits) + exponent - 1
    man = d.scaleb(-exp).normalize()
    return man, exp

# Tour object for storing tours
class Tour:

    __slots__ = ('path', 'length') # faster acces
    graph = None # shared by all tours instances

    def __init__(self, tour = None, length = np.Inf, path = None):
        if tour is None:
            self.length = length
            self.path = path
        else:
            self.length = tour.length
            self.path = tour.path.copy() # need a copy since tour.path might be changed

    def load(self, fname):
        if os.path.exists(fname):
            with open(fname, 'rb') as f:
                self.length, self.path = pickle.load(f)
                return True
        return False

    def save(self, fname):
        if self.isvalid():
            with open(fname, 'wb') as f:
                pickle.dump((self.length, self.path), f)
                return True
        return False

    def isvalid(self):
        return self.length != np.Inf

    @staticmethod
    def init_graph(graph):
        Tour.graph = graph

    def add_phero(self, delta):
        Tour.graph.add_phero_path(self.path, delta / self.length)


# symmetric graph to solve TSP
class Graph:
    isvalid = True

    bestdir = 'best'
    bestFoundTour = None
    optTour = None

    def __init__(self, nNodes):
        self.nNodes = nNodes

        # names
        self.get_names()
        log (f'{self.pname} created with {nNodes} nodes')

        # how many permutations for this TSP ?
        try:
            nPermutations = .5 * math.factorial(nNodes-1)
            self.permutations = getManExp(nPermutations)
        except OverflowError:
            self.permutations = None

        # all matrices are self.nNodes * self.nNodes and symmetrical
        # (i,j) is the edge from node i to node j
        self.matShape = (nNodes, nNodes)

        # create shared dist & weight matrices for the workers use
        self.sharedMatDist = SharedMatMain(self.matShape, 'd')
        self.matDist = self.sharedMatDist.get()

        self.sharedMatWeight = SharedMatMain(self.matShape, 'd')
        self.matWeight = self.sharedMatWeight.get()

        # get matDist
        matDist = self.get_mat_dist()

        # fill diagonal with +infinity to avoid /0 in visibility computation
        np.fill_diagonal(matDist, np.Inf)

        # update the SharedMatWeight for the workers use
        self.sharedMatDist.update(matDist)

        # best found tour
        self.load_best_found()

    # init pheros on the whole graph
    def init_phero(self, phero):
        self.matPhero = np.full(self.matShape, phero) # diagonal is not used

    # init visibility on the whole graph
    def init_visibility(self, beta):
        self.matVisib = (1. / self.matDist) # avoid \0 by having matDist diagonal filled with np.Inf
        self.matVisib **= beta

    # compute weight for each edge
    def compute_weights(self, alpha):
        matWeight = self.matPhero ** alpha
        matWeight *= self.matVisib

        # update the SharedMatWeight for the workers use
        self.sharedMatWeight.update(matWeight)

    # evaporate phero on the whole graph
    def evaporate(self, evaporation):
        self.matPhero *= evaporation

    # clamp phero on the whole graph
    def clamp(self, minPhero, maxPhero):
        self.matPhero.clip(minPhero, maxPhero, out = self.matPhero)

    # add phero on the edges of a path
    def add_phero_path(self, path, delta):
        rows, cols = path[:-1], path[1:]
        self.matPhero[rows, cols] += delta

        # don't forget it's symmetrical !
        self.matPhero[cols, rows] += delta

    # don't need to save if already an optimal tour
    def load_best_found(self):
        if not os.path.exists(self.bestdir):
            os.makedirs(self.bestdir)

        if self.optTour is None:
            self.bestFoundTour = Tour()

            self.bestfile = os.path.join(self.bestdir, self.fname + '.best')
            if self.bestFoundTour.load(self.bestfile):
                log (f'best found tour so far is {self.bestFoundTour.length}km')

    # don't need to save if already an optimal tour
    def save_best_found(self, found):
        if self.optTour is None:
            if found.length < self.bestFoundTour.length:
                self.bestFoundTour = Tour(found)

                if self.bestFoundTour.save(self.bestfile):
                    log (f'best found tour is now {found.length}km')


# graph with matplotlib regular plot using Xs, Ys and nodeLabels
class GraphXY(Graph):

    def __init__(self, nNodes):
        self.nNodes = nNodes

        # nodes coords in a xmax * ymax plane
        self.Xs, self.Ys = self.get_nodes()

        # bbox and labels for plotting
        self.bbox = (self.Xs.min(), self.Ys.min(), self.Xs.max(), self.Ys.max())
        self.nodeLabels = self.get_node_labels()

        super(GraphXY, self).__init__(nNodes)

    # plot a tour
    def plot_tour(self, tour):
        bbox = self.bbox
        ax, _ = plt.gca(), plt.gcf()

        # keep the original ratio
        ax.set_aspect('equal')

        # border around graph in %
        border = .1
        ratio = 1 / (1 + 2 * border) # to map % on [0,1] axes

        # show scale
        scale = .25 # in %
        scalecolor = '.5'
        ax.annotate('', xy=(border * ratio, 0), xycoords='axes fraction',
                    xytext=((border + scale) * ratio, 0),
                    arrowprops=dict(arrowstyle="<->", color=scalecolor, shrinkA=.0, shrinkB=.0))

        xmin, ymin, xmax, ymax = bbox
        mult = {'m': 1000, 'km': 1}
        w, h = xmax - xmin, ymax - ymin

        xscale, yscale = w * scale, h * scale

        unit = 'm' if xscale < 1. else 'km'
        ax.annotate(f'{xscale * mult[unit]} {unit}',
                    xy=((border + scale * .5) * ratio, 0), xycoords='axes fraction',
                    xytext=(0, -10), textcoords='offset pixels',
                    va='center', ha='center', color=scalecolor)

        # background color
        ax.set_facecolor('.98')

        # limits of the plot
        plt.xlim((xmin - border * w, xmax + border * w))
        plt.ylim((ymin - border * h, ymax + border * h))

        # ticks step
        plt.xticks(np.arange(xmin, xmax + xscale * .5, xscale))
        plt.yticks(np.arange(ymin, ymax + yscale * .5, xscale))

        # grid & axis color
        gridcolor = '.93'
        plt.grid(color = gridcolor)
        for spine in ax.spines:
            ax.spines[spine].set_color(gridcolor)

        # no labels & ticks on axis
        ax.tick_params(axis='both', labelbottom=False, labelleft=False, bottom=False, left=False)

        #  optimal or bestFound tour available  ?
        bestTour = None
        if self.optTour:
            bestTour, bestTxt = self.optTour, 'optimal'
        elif self.bestFoundTour.isvalid():
            bestTour, bestTxt = self.bestFoundTour, r'best_{found}'

        # plot title
        title = r'$\bf{%s} TSP$, $n_{nodes}=%d$' %(self.name, self.nNodes) + '\n'

        # plot nodes edges & legends
        txt = r'$found=%gkm$'%tour.length

        # enrich if bestTour
        if bestTour:
            bestLength =  round((1. - bestTour.length / tour.length) * 100, 2)
            txt += r' i.e. $\it%s$ $\bf%+g$%%' %(bestTxt, bestLength)
            title += '\n'

        self.plot_path(tour.path, '-g.', 'blue', txt, annotateNodes = True)
        plt.title(title)

        # plot bestTour
        if bestTour:
            label = r'$%s=%gkm$'%(bestTxt, bestTour.length)
            self.plot_path(bestTour.path, '--', '.6', label, w*.01, h*.01)

        # show legend
        plt.legend(bbox_to_anchor=(0.5, 1), loc='lower center', borderaxespad=0.,
                   ncol = 1 if bestTour else 1, frameon = False, labelspacing= 0)

    # plot all edges of a path
    def plot_path(self, path, linestyle, color, label, dx = 0, dy = 0, annotateNodes = False):
        # plot nodes edges
        Xpath, Ypath = self.Xs[path], self.Ys[path]
        if dx > .0:
            Xpath += dx
        if dy > .0:
            Ypath += dy

        plt.plot(Xpath, Ypath, linestyle, linewidth = .5, color = color, label = label)

        if annotateNodes:
            for p, x, y in zip(path, Xpath, Ypath):
                plt.annotate(f'{self.nodeLabels[p]}', (x, y),
                             xytext = (2, 0), textcoords = 'offset pixels',
                             size = 8, color = 'cornflowerblue')


# graph random generated
class GraphRND(GraphXY):

    def __init__(self, nNodes = 50, seed = 1000, xmax = 1000, ymax = 1000):
        self.seed = seed
        self.xmax, self.ymax = xmax, ymax
        super(GraphRND, self).__init__(nNodes)

    def get_names(self):
        self.pname = f'RndSeed_{self.seed}_{self.nNodes}'           # for print
        self.name = r'RndSeed_{%d, %d}'%(self.seed, self.nNodes)    # for matplotlib
        self.fname = self.pname + f'_({self.xmax}, {self.ymax})'    # for file

    def get_node_labels(self):
        return np.array([i+1 for i in range(self.nNodes)])

    # get random nodes Xs & Ys in a xmax * ymax plane
    def get_nodes(self):
        # same seed gives the same graph
        np.random.seed(self.seed)

        if self.xmax == self.ymax:
            XY = np.random.uniform(0, self.xmax, 2 * self.nNodes)
            return XY[::2], XY[1::2]

        X = np.random.uniform(0, self.xmax, self.nNodes)
        Y = np.random.uniform(0, self.ymax, self.nNodes)
        return  X, Y

    # get mat distance on the whole graph
    def get_mat_dist(self):
        # vector of complex nodes coords reshaped as an 1 * nNodes matrix so we can transpose it
        coords = (self.Xs + self.Ys * 1j)[np.newaxis]

        # np.abs is the euclidian norm on complex
        return np.abs(coords.T - coords)