Получите наиболее оптимальную комбинацию мест из списка мест с помощью алгоритма графа

Я написал короткий пример, используя networkx (основы которого я изучил для целей этого ответа); код, вероятно, неуклюж и его можно сократить, но он должен предоставить хотя бы основу для написания вашего кода.

Я построил график на основе этого рисунка:

import networkx as nx
import math

nodes = {
    "A": [("A_1", (0, 2)), ("A_2", (0, 5)), ("A_3", (0, 8)), ],
    "B": [("B_1", (1, 1)), ("B_2", (1, 3)), ("B_3", (1, 5)), ("B_4", (1, 7)), ],
    "C": [("C_1", (2, 0)), ("C_2", (2, 5)), ("C_3", (2, 9)), ],
    "D": [("D_1", (3, 1)), ("D_2", (3, 5)), ("D_3", (3, 7)), ],
}

G = nx.Graph()

G.add_nodes_from(["in", "out"])


def add_edge(n1, n2):
    origin = G._node[n1]["coords"]
    extremity = G._node[n2]["coords"]
    G.add_edge(n1, n2, weight=math.sqrt((origin[0] - extremity[0])**2 + (origin[1] - extremity[1])**2))


for node in nodes["A"]:
    G.add_node(node[0], coords=node[1])
    G.add_edge("in", node[0], weight=0)

for node in nodes["B"]:
    G.add_node(node[0], coords=node[1])
    for a_node in nodes["A"]:
        add_edge(a_node[0], node[0])

for node in nodes["C"]:
    G.add_node(node[0], coords=node[1])
    for b_node in nodes["B"]:
        add_edge(b_node[0], node[0])

for node in nodes["D"]:
    G.add_node(node[0], coords=node[1])
    G.add_edge(node[0], "out", weight=0)
    for c_node in nodes["C"]:
        add_edge(c_node[0], node[0])

print(nx.shortest_path(G, "in", "out", weight = "weight"))

Результат: ['in', 'A_2', 'B_3', 'C_2', 'D_2', 'out']

Одним из вариантов может быть создание KDTree для каждой категории и запроса чередующегося ближайшего места:

from scipy.spatial import KDTree

kdtrees = {cat: KDTree(list(coords.values())) for cat, coords in places.items()}


def shortest_comb(places, categories, kdtrees):
    min_dis = float("inf")

    for init_name, init_coo in places[categories[0]].items():
        route, total_dis = [init_name], 0
        curr = init_coo

        for next_cat in categories[1:]: # from B and further..
            distance, index = kdtrees[next_cat].query(curr)
            following = list(places[next_cat].keys())[index]
            curr = places[next_cat][following]
            route.append(following)
            total_dis += distance

        if total_dis < min_dis:
            min_dis = total_dis
            mpr = route

    return mpr, min_dis

mpr, min_dis = shortest_comb(places, categories, kdtrees)

print(f"most plausible route: {'->'.join(mpr)}, distance = {min_dis/1000:.2f}km")
# most plausible route: A6->B3->C3->D4, distance=43.50km

Используемые входы:

import random

random.seed(1)


def rnd_points(n, cat, bbox=(21, -32, 22, -31), as_latlon=False):
    min_lon, min_lat, max_lon, max_lat = bbox
    points = [
        (
            random.uniform(min_lat, max_lat),
            random.uniform(min_lon, max_lon)
        ) for _ in range(n)
    ]
    if not as_latlon:
        from pyproj import Transformer
        transformer = Transformer.from_crs(4326, 32734)
        points = [pt for pt in transformer.itransform(points)]
    return {f"{cat}{i}": coo for i, coo in enumerate(points, 1)}


categories = "ABCD"
places = {cat: rnd_points(n, cat) for cat, n in zip(categories, (9, 7, 5, 6))}

Допустим, у вас есть только категории A и B. Тогда поиск оптимального A для каждого B будет стоить вам N(A) * N(B) (в худшем случае — более сложная структура данных вполне может обойтись). Теперь добавьте категорию C. За время N(B) * N( C), для каждого C вы можете попробовать каждый B и объединить его с оптимальным путем к этому B (уже решенным). Теперь для каждого C у вас есть оптимальная тройка A->B->C.

И так далее. Существует количество (в худшем случае) шагов за квадратичное время, на единицу меньше числа категорий. Работая на этапе (скажем) K, вам нужно запомнить N(K) оптимальных путей, по одному для каждого K: для этого K попробуйте каждый J, а затем соедините его с уже вычисленным оптимальным путем, заканчивающимся этим J.

Это форма «динамического программирования».

КОД

Вот полный код, показывающий грубую силу и вышеописанное; предложение намного быстрее:

Обновлено: заменено quicker на более экономичную перезапись:

from random import random, seed
from math import dist, inf

seed(12345678987654321)
def getnum():
    return round(random(), 3)

def getpoints(n):
    return [(getnum(), getnum())
            for i in range(n)]

NCAT = 5
NP = 30
categories = [getpoints(NP) for i in range(NCAT)]

def brute(categories):
    from itertools import product

    smallest = inf
    for t in product(*categories):
        tot = 0.0
        it = iter(t)
        last = next(it)
        for this in it:
            tot += dist(last, this)
            last = this
        if tot < smallest:
            smallest = tot
            best = t
    return smallest, best

def quicker(categories):
    old_dists = [0.0] * len(categories[0])
    old_paths = [[p] for p in categories[0]]
    for i in range(1, len(categories)):
        new_dists = []
        new_paths = []
        for new in categories[i]:
            smallest = inf
            for oldi, (old, old_dist) in enumerate(zip(categories[i-1],
                                                       old_dists)):
                test = dist(old, new) + old_dist
                if test < smallest:
                    smallest = test
                    best_old = oldi
            new_dists.append(smallest)
            new_paths.append(old_paths[best_old] + [new])
        old_dists = new_dists
        old_paths = new_paths
    winner = min(range(len(old_dists)), key=old_dists.__getitem__)
    return old_dists[winner], old_paths[winner]

from time import perf_counter as now
for f in brute, quicker:
    start = now()
    x = f(categories)
    finish = now()
    print(f.__name__, "ran in", finish - start)
    print("   ", x)

Обратите внимание, что я использовал случайное начальное число, поэтому вы должны получить тот же результат, что и я (хотя, конечно, время может отличаться):

brute ran in 18.69626900000003
    (0.17976537535418252, ((0.339, 0.496), (0.333, 0.436),
        (0.333, 0.401), (0.318, 0.477), (0.318, 0.47)))
quicker ran in 0.0008059000028879382
    (0.17976537535418252, [(0.339, 0.496), (0.333, 0.436),
        (0.333, 0.401), (0.318, 0.477), (0.318, 0.47)])

Результат тот же, но quicker на много порядков быстрее.

ОБРАТНЫЙ ТРЕК

Я советую не использовать общий алгоритм графа «кратчайшего пути». В вашем графике огромное количество структур, которые, как указано выше, легко использовать. Другой способ использовать это — выполнить поиск с возвратом по всему пространству продукта, но прерывая его раньше, когда расстояние, для которого это значение, >= лучшего результата, наблюдаемого на данный момент. В моих небольших тестовых случаях он работает на удивление хорошо. YMMV — очень зависит от данных.

def backtrack(categories):
    path = []
    best_path = []
    smallest = inf
    nc = len(categories)

    def inner(i, sofar):
        nonlocal smallest, best_path
        if i >= nc:
            smallest = sofar
            best_path = path[:]
            return
        if i:
            last = path[-1]
        path.append(None)
        for new in categories[i]:
            if i:
                test = sofar + dist(last, new)
            else:
                test = 0;0
            if test < smallest:
                path[-1] = new
                inner(i+1, test)
        path.pop()

    inner(0, 0.0)
    return smallest, best_path

Другой поворот заключается в том, чтобы сначала использовать другой подход, чтобы найти «довольно хороший» путь, чтобы backtrack() мог раньше прервать бесполезные поиски.

Например, храните категории в виде k-d деревьев. Выберите наугад какую-нибудь точку P. Затем постройте путь от ближайшей точки к P в каждой из категорий (дерево k-d может возвращать самое близкое по логарифму времени, и для этой цели должно быть достаточно использовать опцию более быстрого поиска дерева «согласовать аппроксимацию»).

Если повезет, поскольку все точки пути «близки» к P, они также будут близки друг к другу. Вспеньте, промойте и повторяйте, пока не перестанете добиваться «разумного» прогресса. Тогда начните backtrack() с лучшего, что вы видели на данный момент.

РЕАЛИЗМ

Если вы хотите выполнить поиск в литературе, ваша проблема — это пример кратчайшего пути на «направленных, взвешенных, многоступенчатых» графах. «Многоступенчатость» там самое главное качество.

Динамическое программирование — это предпочтительное решение, и я уже привел для него полный код. Это чрезвычайно эффективно, сокращая время выполнения с экспоненциального до (в основном) квадратичного.

Но в комментариях вы говорили о «сотнях миллионов» очков. В этом случае проблема меняется с «сложной» на «ничего страшного — вы никогда не найдете ничего, что могло бы решить столь большие проблемы точно и быстро». quicker() просто непрактично, если у вас есть хотя бы миллион баллов в какой-то категории (квадрат этого числа равен триллиону).

Если вы хотите получить результаты в разумные сроки, вам придется согласиться на приблизительные значения. Насколько я могу судить, нет полезной литературы по аппроксимации кратчайшего пути для многоэтапных графов.

Если вы хотите продолжить это, я уже дал несколько советов. Однако вам придется сделать гораздо больше. На вашем месте я бы запустил два процесса. Один бежит backtrack(). Другой (или более одного) систематически выбирает точки P из сетки, покрывающей пространство XY, и строит пути-кандидаты из ближайших (к P) точек в каждой категории. Если они находят новое лучшее, они убивают процесс backtrack() и перезапускают его с новым лучшим. Процесс, который привел к выигрышу P, возможно, пожелает начать на новой, меньшей сетке с центром в P (приняв своего рода эвристику восхождения на холм).

backtrack() в таких случаях, вероятно, никогда не закончишь, пока у тебя не кончится терпение. В этот момент лучшее, что вы можете сделать, — это согласиться на лучшее, что удалось найти на данный момент. Вы не будете знать, лучшее ли это или даже «близкое». Такие вот дела. Неразрешимое есть неразрешимое :-)

Такие трюки, как генетическое программирование или имитация отжига, мне не кажутся многообещающими. Локальных минимумов предостаточно, а «глубокие» локальные минимумы разделены высокими холмами. «Невозможно добраться отсюда» через правдоподобную последовательность небольших локальных изменений.

ЖАДНЫЙ

Вероятный способ получить «довольно хороший» путь — использовать жадный алгоритм. Выберите один из категории A. Затем выберите тот B, который ближе всего к этому. Затем выберите тот C, который ближе всего к выбранному B. И так далее.

На самом деле это работает не очень хорошо. Однако это можно сделать намного лучше, переписав quicker() так, чтобы после обработки каждой категории запоминались только лучшие LIM пути. Так, например, при работе на этапе K вместо рассмотрения возможностей len(J) * len(K) мы рассматриваем не более LIM * len(K) и запоминаем максимум только LIM лучшие из них. Это делает quicker() линейное время (а не квадратичным), но тогда нет никакой гарантии оптимальности.

Кажется, это хороший способ найти хороший путь, с которого можно начать backtrack(). Поскольку он рандомизирует точки LIM категории A, с которых он начинается, его многократное выполнение может привести к разным путям.

def quicker_lim(categories, LIM=100):
    from heapq import nsmallest
    from random import sample

    old_info = [(0.0, [p])
                for p in sample(categories[0],
                                min(len(categories[0]), LIM))]
    for i in range(1, len(categories)):
        new_info = []
        for new in categories[i]:
            smallest = inf
            for old_dist, p in old_info:
                test = dist(p[-1], new) + old_dist
                if test < smallest:
                    smallest = test
                    bestp = p
            new_info.append((smallest, bestp + [new]))
        old_info = nsmallest(LIM, new_info)
    return old_info[0]

ЖАДНЫЙ к-д

Это наиболее многообещающая схема аппроксимации, которую я придумал. Как и прежде, действуйте поэтапно, по одной категории за раз, сохраняя не более LIM лучшего с каждого этапа.

На следующем этапе для каждого наилучшего пути попробуйте добавить FANOUT ближайшие (к концу пути) точки в текущей категории. Деревья kd делают это быстро и легко, при условии, что дерево kd уже построено. Приведенный ниже код создает их «на лету» и тратит большую часть времени на построение деревьев (для категории с n баллами затраты составляют примерно O(n log n)).

Я использую очень короткий класс на чистом Python KDTree (код здесь не приведен) из

https://github.com/Vectorized/Python-KD-Tree

Это прототип просто для проверки идеи. В производстве, если это вообще возможно, было бы лучше построить деревья заранее и использовать их повторно.

def kdgreedy(categories, LIM=100, FANOUT=3):
    from heapq import nsmallest
    from random import sample
    old_info = [(0.0, [p])
                for p in sample(categories[0],
                                min(len(categories[0]), LIM))]
    for i in range(1, len(categories)):
        new_info = []
        print("cat", i, end = " ")
        kdcat = KDTree(categories[i], 2)
        print("built")
        for old_dist, path in old_info:
            p = path[-1]
            for q in kdcat.get_knn(p, FANOUT, False):
                new_info.append((old_dist + dist(p, q),
                                 path + [q]))
        old_info = nsmallest(LIM, new_info)
    return old_info[0]

Самая большая «дыра» в том, что он пытается набрать максимум LIM очков только из первой категории. Если оптимальный путь не начинается с одного из них, он не будет найден. При повторном запуске будет выбран другой (хотя, возможно, и перекрывающийся) набор отправных точек. Опять же, это предназначено для поиска достойного пути для backtrack() для начала.

И ЖАДНЫЙ К-Д ОТСЛЕЖИВАНИЕ

Идея: сделать возврат, но на каждом уровне стараться не больше, чем LIM ближайших соседей (через k-d деревья). Это не гарантирует оптимальности, но во всех случаях, которые я пробовал до сих пор, оптимальный путь находится гораздо быстрее, чем полный возврат.

kds хранит предварительно вычисленные деревья kd при входе; kds[i] это KDTree(categories[i[, 2).

def backtrack_lim(categories, kds, LIM=5):
    best_path = None
    smallest = inf
    nc = len(kds)
    path = [None] * nc

    def inner(i, sofar):
        nonlocal smallest, best_path
        if i >= nc:
            smallest = sofar
            best_path = path[:]
            print("new smallest", smallest)
            return
        last = path[i-1]
        for new in kds[i].get_knn(last, LIM, False):
            test = sofar + dist(last, new)
            if test < smallest:
                path[i] = new
                inner(i+1, test)

    n = len(categories[0])
    for j, path[0] in enumerate(categories[0]):
        if not j & 0x1ff:
            print(f"{j:_} {j/n:.2%}", end = "\r")
        inner(1, 0.0)
    print()
    return smallest, best_path

Его можно использовать как есть или, например, переработать в генератор, чтобы с самого начала находить всё более лучшие пути для полноценного backtrack().

АСИМПТОТИКА

Это может помочь запомнить следующее:

Во всех случаях brute() занимает время, пропорциональное

prod(map(len, categories))

Катастрофа.

Во всех случаях quicker() занимает время, пропорциональное

sum(len(a) * len(b) for a, b in zip(categories, categories[1:]))

Не только намного быстрее, но и оптимально в худшем случае. Он смотрит на каждое ребро ровно один раз. Чтобы добиться большего, алгоритм должен быть в состоянии сделать вывод, что на некоторые ребра никогда не нужно смотреть. Но в худшем случае все ребра все-таки придется осмотреть.

backtrack() сильно зависит от данных. В худшем случае это может занять до brute() времени, хотя я не уверен, что вообще можно построить такие случаи. По мере увеличения количества баллов в категориях он, скорее всего, будет работать быстрее, чем quicker(). В лучшем случае оно может завершиться за время, пропорциональное len(categories) + len(categories[0]) * len(categories[1]). Например, рассмотрим 100 категорий, которые попарно не пересекаются, за исключением того, что первая точка в каждом списке категорий — p. Первый путь, который он пробует, — это [p] * len(categories), общая длина которого равна 0. Затем все остальные пути отсекаются, как только будут выбраны первые две точки. Алгоритм можно было бы изменить так, чтобы он немедленно останавливался, если бы он обнаружил, что общая длина равна 0, но это почти наверняка будет контрпродуктивно, поскольку циклы записи будут проверяться на то, что почти наверняка никогда не произойдет.

Все остальные методы, приведенные здесь, не гарантируют оптимальности.

ПОСЛЕДУЮЩИЕ ДЕЙСТВИЯ

Во-первых, после дальнейшего тестирования backtrack_lim() является явным победителем, поскольку относительно быстро получает отличные результаты. Он остается эффективным, даже если LIM сократить с 5 до 2.

Во-вторых, большой сюрприз в PyPy 3.10.12: он math.dist намного медленнее, чем CPython, примерно в 8 раз медленнее. PyPy почти всегда значительно быстрее CPython, но не в этом случае. В этом случае это можно исправить, используя вместо этого следующее:

def dist(x, y):
    return sqrt((t := x[0] - y[0]) * t + (t := x[1] - y[1]) * t)

хотя это ограничено 2D-точками и не защищено от ложного переполнения/недополнения для очень больших/маленьких чисел с плавающей запятой. Обратите внимание, что ** 2 вместо умножения t*t работает медленнее - универсальный int.__pow() требует значительных накладных расходов, чтобы добраться до специального (для него) случая «о, нужно только одно умножение».