from matplotlib import pyplot as plt
import numpy as np
import pulp

# 問題を作る（最小化問題）
prob = pulp.LpProblem('example', pulp.LpMinimize)

# 0-1変数を3つ作る
a = pulp.LpVariable('a', cat='Binary')
b = pulp.LpVariable('b', cat='Binary')
c = pulp.LpVariable('c', cat='Binary')

# 目的関数を追加する（最初の += は目的関数と解釈される）
prob += 3 * a + 2 * b + 5 * c

# 制約条件を追加する（2回目以降の += は制約条件と解釈される）
prob += a + b + c == 2

# 解く
prob.solve(pulp.PULP_CBC_CMD(msg=False))
print('求解結果:', pulp.LpStatus[prob.status])

# 変数の値と目的関数の値を取り出す
print('a =', pulp.value(a))
print('b =', pulp.value(b))
print('c =', pulp.value(c))
print('最小値 =', pulp.value(prob.objective))

# 変数 v_1, v_2, ..., v_5 を作ってリストに入れる
v = []
for j in range(5):
    v.append(pulp.LpVariable(f'v_{j+1}', cat='Binary'))

print(v)

# 需要点（街）のデータ：各行が [x座標, y座標, 人口]
demands = [
    [2, 9, 500],
    [5, 9, 300],
    [8, 9, 200],
    [1, 6, 400],
    [4, 6, 700],
    [7, 6, 300],
    [9, 7, 100],
    [2, 3, 600],
    [5, 3, 200],
    [8, 3, 400],
    [3, 1, 100],
    [7, 1, 300],
]

# 配置点の座標：[x座標, y座標]（Part B ではこのまま使う）
pfls = [
    [2, 8],
    [5, 3],
    [7, 7],
    [9, 2],
]

# 配置する店舗の数 k
n_stores = 2

n_demands = len(demands)  # 需要点の数 n
n_pfls = len(pfls)        # 配置点の数 m

def plot_region(demands, pfls, x=None, y=None):
    plt.figure(figsize=(8, 8))
    plt.xlim(0, 11)
    plt.ylim(0, 11)

    # 需要点と店舗を結ぶ線（x が与えられたときだけ描く）
    if x is not None:
        for i in range(len(demands)):
            for j in range(len(pfls)):
                if x[i][j] == 1:
                    plt.plot([demands[i][0], pfls[j][0]],
                             [demands[i][1], pfls[j][1]], color='gray')

    # 需要点（青い丸）と、その番号と人口
    for i in range(len(demands)):
        plt.scatter(demands[i][0], demands[i][1], color='royalblue', s=400)
        plt.text(demands[i][0], demands[i][1], str(i + 1),
                 color='white', ha='center', va='center')
        plt.text(demands[i][0] + 0.25, demands[i][1] + 0.25,
                 str(demands[i][2]), color='royalblue')

    # 配置点（赤い丸）と、その番号。店舗が置かれた点（y[j]=1）は黄色にする
    for j in range(len(pfls)):
        if y is not None and y[j] == 1:
            plt.scatter(pfls[j][0], pfls[j][1], color='gold', s=400)
        else:
            plt.scatter(pfls[j][0], pfls[j][1], color='tomato', s=400)
        plt.text(pfls[j][0], pfls[j][1], str(j + 1),
                 color='white', ha='center', va='center')

    plt.show()


plot_region(demands, pfls)

# 距離 d_ij の計算：distance[i][j] = 需要点i と配置点j の距離
distance = np.zeros((n_demands, n_pfls))
for i in range(n_demands):
    for j in range(n_pfls):
        dx = ...  # x座標の差
        dy = ...  # y座標の差
        distance[i][j] = ...

print(distance)

# 変数の定義
# x[i][j]: 需要点i の住民が配置点j の店舗を利用するなら1
x = []
for i in range(n_demands):
    row = []
    for j in range(n_pfls):
        row.append(pulp.LpVariable(f'x_{i+1}_{j+1}', cat=...))
    x.append(row)

# y[j]: 配置点j に店舗を置くなら1
y = []
for j in range(n_pfls):
    y.append(...)

# 最小化問題として問題を作る
prob = pulp.LpProblem('p-median', pulp.LpMinimize)

# 目的関数：重み付き総移動距離（人口 × 距離 × x_ij の総和）
obj = 0
for i in range(n_demands):
    for j in range(n_pfls):
        obj += ...
prob += obj

# 制約条件1: 各需要点の住民はちょうど1つの店舗を利用する
for i in range(n_demands):
    prob += ...

# 制約条件2: 店舗が置かれていない配置点は利用できない（x_ij <= y_j）
for i in range(n_demands):
    for j in range(n_pfls):
        prob += ...

# 制約条件3: 配置する店舗の数はちょうど n_stores 個
prob += ...

# 求解
prob.solve(pulp.PULP_CBC_CMD(msg=False))
print('求解結果:', pulp.LpStatus[prob.status])

# 最適値（重み付き総移動距離）
f = pulp.value(prob.objective)
print(f'重み付き総移動距離（最適値）: {f:.2f} [人/km]')

# 変数の値（0 か 1）を取り出してリストにする
x_val = []
for i in range(n_demands):
    row = []
    for j in range(n_pfls):
        row.append(int(pulp.value(x[i][j])))
    x_val.append(row)

y_val = []
for j in range(n_pfls):
    y_val.append(int(pulp.value(y[j])))

print('店舗が置かれた配置点 (y):', y_val)

# 各需要点が利用する配置点の一覧（レポート設問4 の表に対応）
print()
print('需要点 ID | 利用する配置点 ID')
for i in range(n_demands):
    for j in range(n_pfls):
        if x_val[i][j] == 1:
            print(f'{i + 1:>8} | {j + 1}')

# 結果の図示（黄色 = 店舗が置かれた配置点、灰色の線 = 誰がどの店舗を使うか）
plot_region(demands, pfls, x=x_val, y=y_val)

レポートの設問	対応する場所
設問1　作成した地域（配置点の座標と店舗数）	Part C
設問2　数理モデル	第1回の解答を自分の設定（$n, m, k$）に合わせて書き直す
設問3　最適値（重み付き総移動距離）	Part C
設問4　需要点の割り当て	Part C
設問5　結果の図示と考察	Part C

記号	意味
$n$	需要点の数
$m$	配置点の数
$k$	配置する店舗の数
$p_i$	需要点 $i$ の人口
$d_{ij}$	需要点 $i$ と配置点 $j$ の距離
$x_{ij}$	需要点 $i$ の住民が配置点 $j$ の店舗を利用するなら1
$y_j$	配置点 $j$ に店舗を配置するなら1

やること	書き方
問題を作る	`prob = pulp.LpProblem('名前', pulp.LpMinimize)`
変数を作る	`a = pulp.LpVariable('a', cat='Binary')`（0か1しか取らない変数）
目的関数を追加する	`prob += 式` （最初の1回目の `+=`）
制約条件を追加する	`prob += 式 == 値` など（2回目以降の `+=`）
解く／値を取り出す	`prob.solve()` ／ `pulp.value(a)`

課題（第2回）：p-median 問題を pulp で実装する¶

ルール¶

0. p-median 問題の数理モデル（第1回の答え合わせ）¶

Part A. pulp 入門（10分）¶

Part B. p-median 問題を実装する¶

B-1. データ¶

B-2. 距離 $d_{ij}$ を計算する【実装1】¶

B-3. 変数 $x_{ij}, y_j$ を作る【実装2】¶

B-4. 目的関数と制約条件を書く【実装3】¶

B-5. 解いて検証する¶

Part C. 自分の地域で解いてレポートを作る¶

考察のヒント（設問5）¶

提出前チェックリスト¶