Ride sharing¶

import numpy as np
from geopy.distance import geodesic


def generate_problem(
    lon_range,
    lat_range,
    parking,
    ncars=None,
    npeople=None,
    C=None,
    lb=1,
    ub=160,
    seed=1,
):
    """
    車の数、人の数、車の定員をランダムに決定した後、
    車の数＋人の数の点の座標を生成し、座標を元に距離行列を生成する関数
    Params
    ------
    lon_range : list
        車と人の座標をサンプルする際の経度の範囲
    lat_range : list
        車と人の座標をサンプルする際の緯度の範囲
    parking : list
        使用可能な駐車場の座標が入ったリスト
    ncars : int (default=None)
        実際に使用する駐車場の数を指定する整数
        ncars <= len(parking)を満たすように設定しなければならない
    npeople : int (default=None)
        ライドシェアの利用者数
    C : int (default=None)
        車の乗車可能人数の上限
    lb : int (default=1)
        ランダムな整数値を生成する時の下限
    seed : int (default=None)
        乱数のシード

    Retuens
    -------
    ncars : int
        実際に使用する駐車場の数を指定する整数
        ncars <= len(parking)を満たすように設定しなければならない
    npeople : int
        ライドシェアの利用者数
    D : np.ndarray
        利用者、車の間の距離を成分に持つ行列
    C : int
        車の乗車可能人数の上限
    ind2coord : dict
        車(人)のインデックスと座標を紐づける辞書
        key : int
            車(人)のインデックス
        value : list
            座標を表すリスト
    """
    np.random.seed(seed)
    if ncars is None or (isinstance(ncars, int) and ncars > len(parking)):
        if isinstance(ncars, int) and ncars > len(parking):
            print(
                f"Maximum value of ncars is {len(parking)}.\n ncars : {ncars} -> {len(parking)}."
            )
        ncars = len(parking)
    if npeople is None:
        npeople = np.random.randint(lb, ub)
    if C is None:
        C = np.random.randint(npeople // ncars + 1, npeople + 2)
    if ncars * C < npeople:
        print("Fail to create valid problem.\nPlease retry after changing random seed.")
        return None, None, None, None, None
    n = ncars + npeople
    ind2coord = dict()
    tmp = [
        parking[i][::-1] for i in np.random.choice(len(parking), ncars, replace=False)
    ]
    for i in range(ncars):
        ind2coord[i] = (tmp[i][0], tmp[i][1])
    for i in range(ncars, n):
        lon = np.random.uniform(lon_range[0], lon_range[1])
        lat = np.random.uniform(lat_range[0], lat_range[1])
        tmp.append((lon, lat))
        ind2coord[i] = (lon, lat)

    D = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            D[i, j] = geodesic(tmp[i][::-1], tmp[j][::-1]).m
    return ncars, npeople, D, C, ind2coord

import folium

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


def simple_plot(coord, ncars):
    m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
    tmp = list(coord.items())
    for j, x in enumerate(tmp):
        if j < ncars:
            folium.Marker(
                location=x[1][::-1],
                icon=folium.Icon(icon="car", prefix="fa", color=_colors[0]),
            ).add_to(m)
        else:
            folium.Marker(
                location=x[1][::-1],
                popup="person",
                icon=folium.Icon(icon="user", prefix="fa", color=_colors[1]),
            ).add_to(m)
    return m

# 船橋駅周辺
lon = (139.9, 140.08)
lat = (35.675500, 35.76)
# 9箇所
parking = [
    (35.67699938102926, 140.0434199237448),
    (35.68494726920934, 139.99303731029542),
    (35.68604762650153, 140.01831984588475),
    (35.69720660219214, 139.98034538800417),
    (35.6981824540223, 140.00360550271415),
    (35.698774929464875, 139.9982410856558),
    (35.700029569368, 139.98558105961536),
    (35.70599837320516, 139.93269833544272),
    (35.71199204224218, 140.0415316476293),
]

ncars, npeople, D, C, index2coordinate = generate_problem(lon, lat, parking, seed=0)
simple_plot(index2coordinate, ncars)

print(ncars, npeople, C)

from amplify import BinarySymbolGenerator

gen = BinarySymbolGenerator()
q = gen.array(npeople, ncars)

from math import sqrt


def regularizeDistance(D):
    average = D.mean(axis=0, keepdims=True)
    std = D.std(axis=0, keepdims=True, ddof=0)
    return (D - average) / std

D = regularizeDistance(D)

from amplify import sum_poly, einsum


def setObjective(q, ncars, npeople, D, C, alpha=1):
    """目的関数"""
    # 各利用者の移動距離に関係する項
    distance_cost = sum_poly(D[ncars:, :ncars] * q)

    # 各車の乗車率に関係する項
    ride_rate_cost = ((q.sum(axis=0) / C) ** 2).sum()

    cost = distance_cost - alpha * ride_rate_cost
    return cost

from amplify.constraint import one_hot, less_equal


def setConstraints(q, ncars, npeople, C, k1=None, k2=None, alpha=1):
    """小規模問題の制約式を設定する関数"""
    if k2 is None:
        k2 = 2 * alpha / C + 1
    if k1 is None:
        k1 = (2 + 2 * alpha / C) + 1

    # 一人一台の車に乗車するという制約(1)
    allocate_constraints = [one_hot(q[i]) for i in range(npeople)]

    # 一台の車に乗れるのはC人以下だという制約(2)
    capacity_constraints = [less_equal(sum_poly(q[:, j]), C) for j in range(ncars)]

    constraints = k1 * sum(allocate_constraints) + k2 * sum(capacity_constraints)
    return constraints

cost = setObjective(q, ncars, npeople, D, C)
constraints = setConstraints(q, ncars, npeople, C)

model1 = cost + constraints

from amplify import Solver
from amplify.client import FixstarsClient

client = FixstarsClient()
client.parameters.timeout = 2000  # 制限時間
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境で使用する場合は、Fixstars Amplify AEのアクセストークンを入力してください

solver = Solver(client)

result = solver.solve(model1)

if len(result) == 0:
    # 実行可能解が見つかっていなければ例外を投げる
    raise RuntimeError("No feasible solution was found.")

q_values = q.decode(result[0].values)

def plot_result(coord, q_values):
    m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
    npeople = len(q_values)
    ncars = len(q_values[0])
    columns = ["latitude", "longitude", "size", "name"]
    data = {label: list() for label in columns}
    answer = dict()
    for i in range(npeople):
        car = np.where(np.array(q_values[i]) == 1)[0][-1]
        if car not in answer:
            answer[car] = []
        answer[car].append(i + ncars)

    for k in range(ncars):
        status = "active"
        car_loc = coord[k]
        if k in answer:
            tmp = answer[k]
            x = [coord[p][0] for p in tmp] + [car_loc[0]]
            y = [coord[p][1] for p in tmp] + [car_loc[1]]
        else:
            x = car_loc[:1]
            y = car_loc[1:]
            status = "empty"
        folium.Marker(
            location=[y[-1], x[-1]],
            popup=f"cluster{k}",
            icon=folium.Icon(icon="car", prefix="fa", color=_colors[k % len(_colors)]),
        ).add_to(m)
        for a, b in zip(y[:-1], x[:-1]):
            folium.Marker(
                location=[a, b],
                popup=f"person{k}",
                icon=folium.Icon(
                    icon="user",
                    prefix="fa",
                    color="white",
                    icon_color=_colors[k % len(_colors)],
                ),
            ).add_to(m)
    return m

plot_result(index2coordinate, q_values)

n = ncars + npeople

q = BinarySymbolGenerator().array(n)

print(D.shape)
print(q.shape)

from amplify import einsum

cost = einsum("ij,i,j->", D, (2 * q - 1), (2 * q - 1))

from amplify.constraint import equal_to


# 分割後の台数が元の台数の半分になるための制約
car_constraints = equal_to(sum_poly(q[:ncars]), ncars // 2)

# 分割後の人数が元の台数の半分になるための制約
people_constraints = equal_to(sum_poly(q[ncars:n]), npeople // 2)

# 制約の強さを設定
k1 = 2 * int(D.sum(axis=1).max()) + 3

constraints = car_constraints + people_constraints
constraints *= k1

model_split = cost + constraints

result = solver.solve(model_split)

if len(result) == 0:
    # 実行可能解が見つからなければ例外を投げる
    raise RuntimeError("No feasible solution was found.")
else:
    # 実行可能解が見つかればそれらの目的関数値を順番に表示する
    for solution in result:
        energy = solution.energy
        values = solution.values
        print(f"energy = {energy}")

def devide(q, D, coord, result):
    """クラスタリングの結果を分割する関数"""
    energy, values = result[0].energy, result[0].values
    q_values = q.decode(values)
    cluster1 = np.where(np.array(q_values) == 1)[0]
    cluster2 = np.where(np.array(q_values) == 0)[0]
    nc1 = len(cluster1)
    nc2 = len(cluster2)
    D1 = np.zeros((nc1, nc1))
    D2 = np.zeros((nc2, nc2))

    C1 = dict()
    C2 = dict()
    for i in range(nc1):
        C1[i] = coord[cluster1[i]]
        for j in range(nc1):
            D1[i][j] = D[cluster1[i]][cluster1[j]]
    for i in range(nc2):
        C2[i] = coord[cluster2[i]]
        for j in range(nc2):
            D2[i][j] = D[cluster2[i]][cluster2[j]]
    return D1, D2, C1, C2

def plot_splited_problem(coord: list, ncars: list):
    m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)

    for i in range(len(ncars)):
        tmp = list(coord[i].items())
        for j, x in enumerate(tmp):
            if j < ncars[i]:
                folium.Marker(
                    location=x[1][::-1],
                    popup=f"cluster{i}",
                    icon=folium.Icon(icon="car", prefix="fa", color=_colors[i]),
                ).add_to(m)
            else:
                folium.Marker(
                    location=x[1][::-1],
                    popup=f"person{i}",
                    icon=folium.Icon(
                        icon="user", prefix="fa", color="white", icon_color=_colors[i]
                    ),
                ).add_to(m)

    return m

D1, D2, cluster1, cluster2 = devide(q, D, index2coordinate, result)
plot_splited_problem([cluster1, cluster2], [ncars // 2, ncars - ncars // 2])

def splitProblem(q, ncars, npeople, D, k1=None):
    """問題を分割するためのモデルを作成する関数"""
    n = ncars + npeople
    if k1 is None:  # 係数が大きすぎると逆に問題がありそうなので、なるべく小さく設定
        k1 = 2 * int(max([sum(D[i]) for i in range(n)])) + 3
    half_cars = ncars // 2
    half_emp = npeople // 2
    cost = einsum("ij,i,j->", D, (2 * q - 1), (2 * q - 1))
    constraints = equal_to(sum_poly(q[:ncars]), half_cars) + equal_to(
        sum_poly(q[ncars:n]), half_emp
    )
    model = cost + k1 * constraints
    return model, half_cars, half_emp

def construct(ncars, npeople, D, C, k1=None, k2=None, alpha=1):
    """分割後の小規模な問題のためのモデルを作成する関数"""
    D = regularizeDistance(D)
    q = BinarySymbolGenerator().array(npeople, ncars)
    cost = setObjective(q, ncars, npeople, D, C, alpha=alpha)
    constraints = setConstraints(q, ncars, npeople, C, k1=k1, k2=k2, alpha=alpha)
    model = cost + constraints
    return model, q

def make_data(coord, q_values, C, last=0, data=None, nframe=1):
    if data is None:
        columns = ["latitude", "longitude", "size", "name", "time", "color"]
        data = {label: list() for label in columns}
    npeople = len(q_values)
    ncars = len(q_values[0])
    answer = dict()
    for i in range(npeople):
        car = np.where(np.array(q_values[i]) == 1)[0][-1]
        if car not in answer:
            answer[car] = []
        answer[car].append(i + ncars)

    loc = [[], []]
    for k in range(ncars):
        status = "active"
        car_loc = coord[k]
        if k in answer:
            tmp = answer[k]
            x = [coord[p][0] for p in tmp] + [car_loc[0]]
            y = [coord[p][1] for p in tmp] + [car_loc[1]]
        else:
            x = car_loc[:1]
            y = car_loc[1:]
            status = "empty"
        loc[0] += y
        loc[1] += x
        for i in range(nframe):
            data["latitude"] += list(
                map(lambda a: ((nframe - i) * a + y[-1] * i) / nframe, y)
            )
            data["longitude"] += list(
                map(lambda a: ((nframe - i) * a + x[-1] * i) / nframe, x)
            )
            data["size"] += [0.5] * (len(x) - 1) + [3]
            data["name"] += [f"group{k+last}({status})"] * len(x)
            data["time"] += [i] * len(x)
            data["color"] += [_colors[k + last]] * len(x)
    return data

def plot_final_result(data):
    m = folium.Map([sum(lat) / 2, sum(lon) / 2], tiles="OpenStreetMap", zoom_start=12)
    df = pd.DataFrame(data)
    for _name in data["name"]:
        tmp = df[df["name"] == _name]
        x = list(tmp["longitude"])
        y = list(tmp["latitude"])
        folium.Marker(
            location=[y[-1], x[-1]],
            popup=f"cluster_{_name}",
            icon=folium.Icon(icon="car", prefix="fa", color=list(tmp["color"])[-1]),
        ).add_to(m)
        for a, b, c in zip(y[:-1], x[:-1], list(tmp["color"])[:-1]):
            folium.Marker(
                location=[a, b],
                popup=f"person_{_name}",
                icon=folium.Icon(icon="user", prefix="fa", color="white", icon_color=c),
            ).add_to(m)
    return m

from collections import namedtuple

Parameter = namedtuple("Config", ("bit_size", "alpha"))

param = Parameter(bit_size=50, alpha=10)

from collections import deque
from IPython.display import display
import pandas as pd


def solve(ncars, npeople, D, C, coord, param, debug=False):
    print("問題設定")
    display(simple_plot(coord, ncars))
    client = FixstarsClient()
    client.parameters.timeout = 200  # 制限時間
    solver = Solver(client)
    print("問題を分割しています...", end="")
    queue = deque([(ncars, npeople, D, coord)])
    while (queue[0][1] + C) * queue[0][0] > param.bit_size:
        (ncars, npeople, D, coord) = queue.popleft()

        q = BinarySymbolGenerator().array(ncars + npeople)

        model, ncars1, npeople1 = splitProblem(q, ncars, npeople, D)
        result = solver.solve(model)
        if len(result) == 0:
            raise RuntimeError("No feasible solution was found.")
        D1, D2, C1, C2 = devide(q, D, coord, result)
        queue.append((ncars1, npeople1, D1, C1))
        queue.append((ncars - ncars1, npeople - npeople1, D2, C2))
    print("完了")
    print("結果を描写しています...", end="")
    m = plot_splited_problem([x[-1] for x in queue], [x[0] for x in queue])
    display(m)
    print("完了")

    print("分割後の問題を解いています...")
    client = FixstarsClient()
    client.parameters.timeout = 1000  # 制限時間
    solver = Solver(client)
    index = 0
    last = 0
    data = None
    while queue:
        index += 1
        (ncars, npeople, D, coord) = queue.pop()
        model, q = construct(ncars, npeople, D, C, alpha=param.alpha)
        result = solver.solve(model)
        if len(result) == 0:
            raise RuntimeError("No feasible solution was found.")
        print(f"子問題{index}が解けました")
        q_values = q.decode(result[0].values)
        data = make_data(coord, q_values, C, data=data, last=last)
        last += ncars
    print("結果を描写しています...")
    m = plot_final_result(data)
    display(m)

problem = generate_problem(lon, lat, parking, C=12, seed=0)

solve(*problem, param)