結果
| 問題 |
No.5007 Steiner Space Travel
|
| コンテスト | |
| ユーザー |
 ra5anchor
|
| 提出日時 | 2025-04-20 23:27:40 |
| 言語 | PyPy3 (7.3.15) |
| 結果 |
AC
|
| 実行時間 | 933 ms / 1,000 ms |
| コード長 | 15,801 bytes |
| コンパイル時間 | 587 ms |
| コンパイル使用メモリ | 82,172 KB |
| 実行使用メモリ | 87,996 KB |
| スコア | 8,510,552 |
| 最終ジャッジ日時 | 2025-04-20 23:28:11 |
| 合計ジャッジ時間 | 30,534 ms |
|
ジャッジサーバーID （参考情報） |
judge1 / judge5 |
| 純コード判定しない問題か言語 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 30 |
ソースコード
import sys
import math
from time import perf_counter
import random
import heapq
import argparse
N = 100
DISTS = [[0] * N for i in range(N)]
DSIZE = 0
class TimeKeeper:
def __init__(self):
self.start_time = perf_counter()
def is_time_over(self, LIMIT):
return (perf_counter() - self.start_time) >= LIMIT
def time_now(self):
return (perf_counter() - self.start_time)
class UnionFind():
def __init__(self, n):
self.n = n
self.parents = [-1] * n
def find(self, x):
if self.parents[x] < 0:
return x
else:
self.parents[x] = self.find(self.parents[x])
return self.parents[x]
def union(self, x, y):
x = self.find(x)
y = self.find(y)
if x == y:
return
if self.parents[x] > self.parents[y]:
x, y = y, x
self.parents[x] += self.parents[y]
self.parents[y] = x
def same(self, x, y):
return self.find(x) == self.find(y)
class Reverse:
def __init__(self, val):
self.val = val
def __lt__(self, other):
return self.val > other.val
def __repr__(self):
return repr(self.val)
class Heapq:
def __init__(self, arr, desc = False):
if desc:
for i in range(len(arr)):
arr[i] = Reverse(arr[i])
self.desc = desc
self.hq = arr
heapq.heapify(self.hq)
def pop(self):
if self.desc: return heapq.heappop(self.hq).val
else: return heapq.heappop(self.hq)
def push(self, a):
if self.desc: heapq.heappush(self.hq, Reverse(a))
else: heapq.heappush(self.hq, a)
def top(self):
if self.desc: return self.hq[0].val
else: return self.hq[0]
def size(self):
return len(self.hq)
def calc_total_dist(tour):
res = DISTS[tour[len(DISTS) - 1]][tour[0]]
for i in range(len(DISTS) - 1):
res += DISTS[tour[i]][tour[i + 1]]
return res
def calscore(stations, tour, ALPHA, id_to_xy):
s = 0
for (t0,i0),(t1,i1) in zip(tour, tour[1:]):
# i0 -= 1; i1 -= 1
if t0 == 1:
x0,y0 = id_to_xy[i0]
else:
x0,y0 = stations[i0]
if t1 == 1:
x1,y1 = id_to_xy[i1]
else:
x1,y1 = stations[i1]
d2 = (x0-x1)**2 + (y0-y1)**2
coef = {2:ALPHA**2, 3:ALPHA, 4:1}[t0+t1]
s += d2 * coef
return int(10**9 / (1000 + math.sqrt(s)) + 0.5)
def christofides():
global DISTS, DSIZE
# 最小全域木ベースとマッチングにより初期解の構築
start_time = perf_counter()
N = len(DISTS)
# 考慮する辺の枝刈り
coef = 200 #? 得られる辺の数=(coef * N)/2
if N <= coef:
threshold = DSIZE
else:
threshold = max(DSIZE * (1 - math.sqrt(1 - coef / N)), 0.05 * DSIZE)
# create MST by Kruskal
priority_edges = []
for i in range(N):
for j in range(i+1, N):
if DISTS[i][j] <= threshold:
priority_edges.append((DISTS[i][j], i, j))
priority_edges.sort()
G = [[] for _ in range(N)]
# selected_edges = []
deg = [0] * N
uf = UnionFind(N)
edge_cnt = 0
for w, a, b in priority_edges:
if not uf.same(a, b):
uf.union(a, b)
G[a].append(b)
G[b].append(a)
# selected_edges.append((a, b, w)) #! 使ってない
deg[a] += 1
deg[b] += 1
edge_cnt += 1
if edge_cnt == N-1:
break
# find odd-degree nodes from MST
odd_nodes = set()
for i in range(N):
if deg[i] % 2 == 1:
odd_nodes.add(i)
# print("number of odd node:", len(odd_nodes), file=sys.stderr)
# add matching edges to MST
# １つずつ貪欲にマッチング #!最適マッチングとは限らない
NN = len(odd_nodes)
min_length = 10**18
for iloop in range(NN):
not_used_nodes = list(odd_nodes)
random.shuffle(not_used_nodes)
matched_edges = []
tot_length = 0
while len(not_used_nodes) > 0:
v = not_used_nodes.pop()
length = 10**18
for u in not_used_nodes:
if v != u and DISTS[v][u] < length:
length = DISTS[v][u]
closest = u
u = closest
matched_edges.append((v, u, length))
tot_length += length
not_used_nodes.remove(u)
# print(tot_length, file=sys.stderr)
if tot_length < min_length:
min_length = tot_length
best_edges = matched_edges[:]
for v, u, length in best_edges:
G[u].append(v)
G[v].append(u)
# print("loop", iloop, file=sys.stderr)
def create_eulerian_tour(G0, start, reverse):
# オイラーツアー
# find an eulerian tour ([v0, v1, v2, ..., vn-1, v0])
tour = []
stack = [start]
G = [[] for _ in range(N)]
if reverse == False:
for i in range(N):
G[i] = G0[i][:]
else:
for i in range(N):
G[i] = G0[i][::-1]
while stack:
v = stack[-1]
if G[v]:
u = G[v].pop()
# remove the reverse direction as well
G[u].remove(v)
stack.append(u)
else:
# No more edges, add to tour
stack.pop()
tour.append(v)
return tour[::-1] # Reverse to get the correct order
def create_hamilton_tour(tour_in, reverse):
# 重複する頂点があれば飛ばしてハミルトン閉路を作る
visited = [False] * N
tour = []
tour_in = tour_in[::-1] if reverse else tour_in
for v in tour_in:
if not visited[v]:
visited[v] = True
tour.append(v)
return tour
min_dist = 10**18
start_node = 0
for rev1 in [False, True]:
for rev2 in [False, True]:
eulerian_tour = create_eulerian_tour(G, start_node, reverse=rev1)
tour = create_hamilton_tour(eulerian_tour, reverse=rev2)
tot_dist = calc_total_dist(tour)
if tot_dist < min_dist:
min_dist = tot_dist
best_tour = tour
return min_dist, best_tour
def opt_mix(tour_in, LIMIT, temp0, temp1):
global DISTS, DSIZE
N = len(DISTS)
start_time = perf_counter()
temp = temp0
tour = tour_in[:]
loop = 0
cnt2, cnt3 = 0, 0
exp_ = [0] * (10**5+1)
for i in range(10**5+1):
exp_[i] = math.exp(-i / 10**4) # exp(x/10**4)
exp_[10**5] = 0
X = 4
while True:
loop += 1
if loop % 10000 == 0:
tnow = perf_counter()
progress = min((tnow - start_time) / LIMIT, 1.0)
temp = temp0 + (temp1 - temp0) * progress ** 1.0
if tnow - start_time > LIMIT:
break
if loop % X >= 1:
# 2-opt
cnt2 += 1
a = (cnt2 // N) % N
b = cnt2 % N
if a > b:
a, b = b, a
c = (a + 1) % N
d = (b + 1) % N
if a == b or a == d:
continue
# 焼きなまし
sc0 = DISTS[tour[a]][tour[c]] + DISTS[tour[b]][tour[d]]
sc1 = DISTS[tour[a]][tour[b]] + DISTS[tour[c]][tour[d]]
if sc1 < sc0:
probability = 1.0
# print("opt2", loop, int(sc0-sc1), file=sys.stderr)
else:
# probability_ = math.exp((sc0 - sc1) / temp)
x = int(-(sc0 - sc1) * 10**4 / temp)
if x > 10**5:
continue
probability = exp_[x]
if probability == 1.0 or probability > random.random():
# 変更
x, y = c, b
if b - c > a + N - d:
x, y = d, a + N
while x < y:
tour[x%N], tour[y%N] = tour[y%N], tour[x%N]
x += 1
y -= 1
else:
# 3-opt
cnt3 += 1
a = (cnt3 // N) % N
b = (a + 1) % N
c = (a + 2) % N
d = cnt3 % N
e = (d + 1) % N
if d in {a, b, c}:
continue
ab = DISTS[tour[a]][tour[b]]
bc = DISTS[tour[b]][tour[c]]
de = DISTS[tour[d]][tour[e]]
ac = DISTS[tour[a]][tour[c]]
db = DISTS[tour[d]][tour[b]]
be = DISTS[tour[b]][tour[e]]
# 焼きなまし
sc0 = ab + bc + de
sc1 = ac + db + be
if sc1 < sc0:
probability = 1.0
# print("opt3", loop, int(sc0-sc1), file=sys.stderr)
else:
# probability_ = math.exp((sc0 - sc1) / temp)
x = int(-(sc0 - sc1) * 10**4 / temp)
if x > 10**5:
continue
probability = exp_[x]
if probability == 1.0 or probability > random.random():
# 変更
city = tour.pop(b)
if b < e:
tour.insert(d, city)
else:
tour.insert(e, city)
tot_dist = calc_total_dist(tour)
print(" SA loop:", loop, file=sys.stderr)
return tot_dist, tour
# K-meansクラスタリング
def euclidean_distance(p1, p2):
return sum((a - b) ** 2 for a, b in zip(p1, p2)) ** 0.5
def initialize_centroids(points, k):
return random.sample(points, k)
def assign_clusters(points, centroids):
clusters = []
for p in points:
ds = [euclidean_distance(p, c) for c in centroids]
clusters.append(ds.index(min(ds)))
return clusters
def update_centroids(points, clusters, k):
new_centroids = []
for i in range(k):
pts = [points[j] for j in range(len(points)) if clusters[j] == i]
if pts:
avg = [sum(dim)/len(pts) for dim in zip(*pts)]
new_centroids.append(tuple(avg))
return new_centroids
def kmeans(points, k, max_iters=100):
centroids = initialize_centroids(points, k)
for _ in range(max_iters):
clusters = assign_clusters(points, centroids)
new_centroids = update_centroids(points, clusters, k)
if centroids == new_centroids:
break
centroids = new_centroids
return clusters, centroids
def main(DEBUG):
global DISTS, DSIZE
def distance(city1, city2):
return math.sqrt((city1[0] - city2[0])**2 + (city1[1] - city2[1])**2)
tk = TimeKeeper()
LIMIT = 0.8 # 制限時間(s)
LIMIT1 = 0.4
random.seed(0)
ALPHA = 5
# 入力
N, M = map(int, input().split())
points = [tuple(map(int, input().split())) for _ in range(N)]
for i in range(N):
for j in range(i+1, N):
DISTS[i][j] = distance(points[i], points[j])
DISTS[j][i] = DISTS[i][j]
DSIZE = 1000 * 1.41
id_to_xy = dict()
for i,(x,y) in enumerate(points):
id_to_xy[i] = (x,y)
d1, tour = christofides() # MSTベース
# temp:焼きなましの温度
temp0 = (DSIZE / N**0.5) / 10
temp1 = temp0 / 10
t_remain = LIMIT1 - tk.time_now()
d2, tour = opt_mix(tour, t_remain, temp0, temp1) # 焼きなまし
# （フォーマット）訪問開始を1にしてlabel付け
stt_idx = tour.index(0)
tour = tour[stt_idx:] + tour[:stt_idx] + [0]
tour = [(1,x) for x in tour]
tour0 = tour[:]
# tour0 : ステーションなしの
##############
# ステーション初期配置
stations = []
labels, cents = kmeans(points, M)
for i in range(M):
x, y = cents[i]
stations.append((int(x), int(y)))
found = True
while found == True:
found = False
for i in range(len(tour)-1):
label0, id0 = tour[i]
label1, id1 = tour[i+1]
if label0 == 2 or label1 == 2:
continue
# print(id0, id1, file=sys.stderr)
d_org = DISTS[id0][id1]
d_min = 10**18
st_min = -1
for j in range(M):
d = (distance(points[id0], stations[j]) +
distance(points[id1], stations[j]))
if d < d_min:
d_min = d
st_min = j
if d_min < d_org * ALPHA:
tour.insert(i+1, (2, st_min))
found = True
break
sc = calscore(stations, tour, ALPHA, id_to_xy)
current_sc = sc
best_sc = sc
best_stations = stations[:] # 最良解
best_tour = tour[:] # 最良ツアー
# ステーション配置を変えていく
loop = 0
DD0 = 500
time0 = tk.time_now()
temp0 = max(500, 1)
while not tk.is_time_over(LIMIT):
loop += 1
elapsed = tk.time_now() - time0
progress = min(elapsed / (LIMIT - time0), 1.0)
# 線形クーリング
temp = temp0 * (1 - progress)
if temp < 1e-6:
temp = 1e-6
DD = max(10, min(500, int(DD0 * (1 - progress))))
# --- 近傍解の生成 ---
nst = random.randrange(M)
x0, y0 = stations[nst]
dx = random.randint(-DD, DD)
dy = random.randint(-DD, DD)
nx = max(0, min(1000, x0 + dx))
ny = max(0, min(1000, y0 + dy))
stations[nst] = (nx, ny) # 仮に更新
# stations_org = stations[:]
# for i in range(M):
# x0, y0 = stations[i]
# dx = random.randint(-DD, DD)
# dy = random.randint(-DD, DD)
# nx = max(0, min(1000, x0 + dx))
# ny = max(0, min(1000, y0 + dy))
# stations[i] = (nx, ny) # 仮に更新
tour = tour0[:] # ステーションなしの経路
found = True
while found == True:
found = False
for i in range(len(tour)-1):
label0, id0 = tour[i]
label1, id1 = tour[i+1]
if label0 == 2 or label1 == 2:
continue
# print(id0, id1, file=sys.stderr)
d_org = DISTS[id0][id1]
d_min = 10**18
st_min = -1
for j in range(M):
d = (distance(points[id0], stations[j]) +
distance(points[id1], stations[j]))
if d < d_min:
d_min = d
st_min = j
if d_min < d_org * ALPHA:
tour.insert(i+1, (2, st_min))
found = True
break
new_sc = calscore(stations, tour, ALPHA, id_to_xy)
delta = new_sc - current_sc
if delta >= 0 or math.exp(delta / temp) > random.random():
# 受け入れ
current_sc = new_sc
if new_sc > best_sc:
best_sc = new_sc
best_stations = stations[:]
best_tour = tour[:]
print(f" best loop={loop} sc={best_sc}", file=sys.stderr)
else:
# 棄却 → 元に戻す
# stations = stations_org[:]
stations[nst] = (x0, y0)
# 出力
print('SC', best_sc, file=sys.stderr)
for x,y in best_stations:
print(x, y)
print(len(best_tour))
for label, id in best_tour:
print(label, id+1)
if __name__ == '__main__':
parser = argparse.ArgumentParser(description='Debug mode')
parser.add_argument('--debug', action='store_true', help='Enable debug mode')
args = parser.parse_args()
main(args.debug)