import numpy as np
from geopy.distance import geodesic
    
# 車両数
nvehicle = 2
    
# 都市の数
ncity = 10
    
dimension = ncity + 1
avg_cities_per_vehicle = dimension // nvehicle
    
# 各都市における配送需要（重量）を乱数で決定
demand = np.random.randint(1, 100, size=dimension)
demand[0] = 0
demand_max = max(demand)
demand_mean = sum(demand) / len(demand)
    
# 全体的な需要に合わせ、車両の積載可能量 Q を設定する。
capacity = demand_max + int(demand_mean) * avg_cities_per_vehicle

# 座標の取り得る範囲を設定
lat_range = [35.7014, 35.968]
lon_range = [139.34, 140.04]
    
# 各都市の座標
ind2coord = {
    i: (
        np.random.uniform(lon_range[0], lon_range[1]),
        np.random.uniform(lat_range[0], lat_range[1]),
    )
    for i in range(dimension)
}
    
# 2都市間の座標距離行列 D
distance_matrix = np.array(
    [
        [geodesic(coord_i[::-1], coord_j[::-1]).m for coord_j in ind2coord.values()]
        for coord_i in ind2coord.values()
    ]
)

_colors = [
    "green",
    "orange",
    "blue",
    "red",
    "purple",
    "pink",
    "darkblue",
    "cadetblue",
    "darkred",
    "lightred",
    "darkgreen",
    "lightgreen",
    "lightblue",
    "darkpurple",
]

import folium
    
    
def plot_solution(coord: dict, title: str, best_tour: dict = dict()):
    l = len(coord)
    center = [
        sum(lat for _, lat in coord.values()) / l,
        sum(lon for lon, _ in coord.values()) / l,
    ]
    m = folium.Map(center, tiles="Stamen Toner", zoom_start=10.5)
    folium.Marker(
        location=coord[0][::-1],
        popup=f"depot",
        icon=folium.Icon(icon="car", prefix="fa"),
    ).add_to(m)
    
    _color = _colors[1]
    if best_tour:
        for k, tour in best_tour.items():
            _color = _colors[k % len(_colors)]
            for city in tour:
                if city == 0:
                    continue
    
                folium.Marker(
                    location=coord[city][::-1],
                    popup=f"person{k}",
                    icon=folium.Icon(
                        icon="user", prefix="fa", color="white", icon_color=_color
                    ),
                ).add_to(m)
            folium.vector_layers.PolyLine(
                locations=[coord[city][::-1] for city in tour], color=_color, weight=3
            ).add_to(m)
    else:
        for k, node in coord.items():
            if k == 0:
                continue
            folium.Marker(
                location=node[::-1],
                popup=f"customer{k}",
                icon=folium.Icon(
                    icon="user", prefix="fa", color="white", icon_color=_color
                ),
            ).add_to(m)
    
    title = f"<h4>{title}</h4>"
    m.get_root().html.add_child(folium.Element(title))
    
    return m

# 埼玉・千葉・東京近辺
title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}"
plot_solution(ind2coord, title)

def upperbound_of_tour(capacity: int, demand: np.ndarray) -> int:
    max_tourable_cities = 1
    for w in sorted(demand):
        capacity -= w
        if capacity >= 0:
            max_tourable_cities += 1
        else:
            return max_tourable_cities
    return max_tourable_cities

from amplify import BinaryPoly, BinarySymbolGenerator
    
gen = BinarySymbolGenerator()
    
# 積載可能量から1台の車両が訪問できる都市の最大数
max_tourable_cities = upperbound_of_tour(capacity, demand)
    
x = gen.array(max_tourable_cities, dimension, nvehicle)

x[0][1:][:] = 0
x[-1][1:][:] = 0
    
x[0][0][:] = 1
x[-1][0][:] = 1

from amplify.constraint import one_hot
from amplify import sum_poly
    
    
def make_onetrip_constraints(x) -> list:
    max_tourable_cities = x.shape[0]
    dimension = x.shape[1]
    nvehicle = x.shape[2]
    
    constraints = [
        one_hot(sum_poly(dimension, lambda j: x[i][j][k]))
        for i in range(max_tourable_cities)
        for k in range(nvehicle)
    ]
    return constraints

def make_onevisit_constraints(x) -> list:
    max_tourable_cities = x.shape[0]
    dimension = x.shape[1]
    nvehicle = x.shape[2]
    
    constraints = [
        one_hot(
            sum_poly(
                max_tourable_cities, lambda i: sum_poly(nvehicle, lambda k: x[i][j][k])
            )
        )
        for j in range(1, dimension)
    ]
    return constraints

from amplify.constraint import less_equal
from amplify import InequalityFormulation
    
    
def make_capacity_constraints(x, demand: np.ndarray, capacity: int) -> list:
    max_tourable_cities = x.shape[0]
    dimension = x.shape[1]
    nvehicle = x.shape[2]
    
    constraints = [
        less_equal(
            demand * x[:, :, k],
            capacity,
            method=InequalityFormulation.RelaxationQuadra,
        )
        / capacity
        / capacity
        for k in range(nvehicle)
    ]
    
    return constraints

from amplify import einsum
    
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
    
xr = x.roll(-1, axis=0)
# 経路の総距離
objective = einsum("pq,ipk,iqk", distance_matrix, x, xr)

constraints1 = make_onetrip_constraints(x)
constraints2 = make_onevisit_constraints(x)
constraints3 = make_capacity_constraints(x, demand, capacity)
constraints = sum(constraints1) + sum(constraints2) + sum(constraints3)
    
model = objective + constraints * np.max(distance_matrix)

from amplify.client import FixstarsClient
    
client = FixstarsClient()
client.parameters.timeout = 2000  # タイムアウト2秒
# client.token = "API トークンを入力してください"

from amplify import Solver
    
solver = Solver(client)
    
result = solver.solve(model)
if len(result.solutions) == 0:
    raise RuntimeError("Some of the constraints are not satisfied.")
    
x_values = result.solutions[0].values  # 変数値の取り出し

def onehot2_sequence(x_values) -> dict[int, list]:
    nvehicle = x_values.shape[2]
    sequence = dict()
    for k in range(nvehicle):
        sequence[k] = np.where(x_values[:, :, k])[1]
    return sequence

def process_sequence(sequence: dict[int, list]) -> dict[int, list]:
    new_seq = dict()
    for k, v in sequence.items():
        v = np.append(v, v[0])
        mask = np.concatenate(([True], np.diff(v) != 0))
        new_seq[k] = v[mask]
    return new_seq

solution = x.decode(x_values)
sequence = onehot2_sequence(
    solution
)  # one-hot な変数テーブルを辞書に変換。key：車両インデックス, value：各車両が訪問した都市の順番が入ったリスト
best_tour = process_sequence(sequence)  # 上の辞書からデポへの余計な訪問を取り除く
print(result.solutions[0].energy)  # 目的関数値を表示
print(*best_tour.items(), sep="\n")  # 求めた解を表示

title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}, cost={result.solutions[0].energy:.2f}"
plot_solution(ind2coord, title, best_tour)

$x_{i,j,1}$	depot	A	C	E	F
0(出発)	1	0	0	0	0
1	0	1	0	0	0
2	0	0	1	0	0
3	0	0	0	1	0
4	0	0	0	0	1
5	1	0	0	0	0
6	1	0	0	0	0
7	1	0	0	0	0
8	1	0	0	0	0
9(到着)	1	0	0	0	0

$x_{i,j,2}$	depot	B	G	I
0(出発)	1	0	0	0
1	1	0	0	0
2	1	0	0	0
3	1	0	0	0
4	0	1	0	0
5	0	0	1	0
6	0	0	0	1
7	1	0	0	0
8	1	0	0	0
9(到着)	1	0	0	0

$x_{i,j,3}$	depot	D	H
0(出発)	1	0	0
1	0	1	0
2	0	0	1
3	1	0	0
4	1	0	0
5	1	0	0
6	1	0	0
7	1	0	0
8	1	0	0
9(到着)	1	0	0

容量制約つき運搬経路問題(Capacitated Vehicle Routing Problem, CVRP)¶

定式化¶

定数¶

変数¶

制約条件¶

目的関数¶

具体例を使った制約式の説明¶

実装¶

問題の作成・可視化¶

訪問都市の上限の計算¶

二次多項式モデルの構築¶

制約条件①：全ての車両はデポ（都市 0）から出発し、デポに到着する¶

制約条件②：車両$k$が$i$番目に訪問する場所は一箇所だけ¶

制約条件③：都市$j$はいずれかの車両によってちょうど一回だけ訪問される¶

制約条件④：車両 $k$ が運ぶ荷物の総重量は $Q$ 以下¶

目的関数¶

Amplify クライアントの設定¶

結果取得・可視化¶

容量制約つき運搬経路問題(Capacitated Vehicle Routing Problem, CVRP)¶

定式化¶

定数¶

変数¶

制約条件¶

目的関数¶

具体例を使った制約式の説明¶

実装¶

問題の作成・可視化¶

訪問都市の上限の計算¶

二次多項式モデルの構築¶

制約条件①： 全ての車両はデポ（都市 0）から出発し、デポに到着する¶

制約条件②： 車両$k$が$i$番目に訪問する場所は一箇所だけ¶

制約条件③： 都市$j$はいずれかの車両によってちょうど一回だけ訪問される¶

制約条件④： 車両 $k$ が運ぶ荷物の総重量は $Q$ 以下¶

目的関数¶

Amplify クライアントの設定¶

結果取得・可視化¶

制約条件①：全ての車両はデポ（都市 0）から出発し、デポに到着する¶

制約条件②：車両$k$が$i$番目に訪問する場所は一箇所だけ¶

制約条件③：都市$j$はいずれかの車両によってちょうど一回だけ訪問される¶

制約条件④：車両 $k$ が運ぶ荷物の総重量は $Q$ 以下¶