結果
| 問題 | No.5007 Steiner Space Travel |
| ユーザー |
 r3yohei
|
| 提出日時 | 2023-07-11 22:17:40 |
| 言語 | Rust (1.77.0 + proconio) |
| 結果 |
AC
|
| 実行時間 | 999 ms / 1,000 ms |
| コード長 | 39,303 bytes |
| コンパイル時間 | 13,554 ms |
| コンパイル使用メモリ | 185,652 KB |
| 実行使用メモリ | 4,504 KB |
| スコア | 9,049,879 |
| 最終ジャッジ日時 | 2023-07-11 22:18:27 |
| 合計ジャッジ時間 | 37,302 ms |
|
ジャッジサーバーID （参考情報） |
judge11 / judge13 |
| 純コード判定しない問題か言語 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 998 ms
4,380 KB |
| testcase_01 | AC | 998 ms
4,384 KB |
| testcase_02 | AC | 998 ms
4,376 KB |
| testcase_03 | AC | 999 ms
4,380 KB |
| testcase_04 | AC | 998 ms
4,376 KB |
| testcase_05 | AC | 998 ms
4,376 KB |
| testcase_06 | AC | 998 ms
4,380 KB |
| testcase_07 | AC | 998 ms
4,380 KB |
| testcase_08 | AC | 997 ms
4,376 KB |
| testcase_09 | AC | 998 ms
4,380 KB |
| testcase_10 | AC | 998 ms
4,380 KB |
| testcase_11 | AC | 998 ms
4,376 KB |
| testcase_12 | AC | 998 ms
4,376 KB |
| testcase_13 | AC | 998 ms
4,380 KB |
| testcase_14 | AC | 998 ms
4,380 KB |
| testcase_15 | AC | 998 ms
4,376 KB |
| testcase_16 | AC | 998 ms
4,380 KB |
| testcase_17 | AC | 999 ms
4,376 KB |
| testcase_18 | AC | 998 ms
4,380 KB |
| testcase_19 | AC | 999 ms
4,376 KB |
| testcase_20 | AC | 998 ms
4,380 KB |
| testcase_21 | AC | 998 ms
4,376 KB |
| testcase_22 | AC | 999 ms
4,376 KB |
| testcase_23 | AC | 999 ms
4,376 KB |
| testcase_24 | AC | 998 ms
4,376 KB |
| testcase_25 | AC | 999 ms
4,376 KB |
| testcase_26 | AC | 998 ms
4,376 KB |
| testcase_27 | AC | 998 ms
4,504 KB |
| testcase_28 | AC | 998 ms
4,376 KB |
| testcase_29 | AC | 999 ms
4,376 KB |
ソースコード
#![allow(non_snake_case, unused)]
use self::i64::kmeans;
use crate::input::{marker::*, *};
use std::{cmp::*, vec};
use std::{collections::*, fmt::format};
use std::{hash, io::*};
use std::env;
// proconio
pub mod input {
use std::{
cell::RefCell,
fmt::Debug,
io::{stdin, BufRead, BufReader, Stdin},
str::{FromStr, SplitWhitespace},
};
thread_local!(
pub static STDIN_SOURCE: RefCell<Source> = RefCell::new(Source {
stdin: BufReader::new(stdin()),
tokens: "".split_whitespace(),
});
);
pub struct Source {
stdin: BufReader<Stdin>,
tokens: SplitWhitespace<'static>,
}
impl Source {
pub fn next_token(&mut self) -> Option<&str> {
self.tokens.next().or_else(|| {
let mut input = String::new();
self.stdin.read_line(&mut input).unwrap();
self.tokens = Box::leak(input.into_boxed_str()).split_whitespace();
self.tokens.next()
})
}
}
#[macro_export]
macro_rules! read_value {
(from $s:expr, [$t:tt; $n:expr]) => {
(0..$n).map(|_| $crate::read_value!(from $s, $t)).collect::<Vec<_>>()
};
(from $s:expr, [$t:tt]) => {
let n = $crate::read_value!(from $s, usize);
$crate::read_value!(from $s, [$t; n])
};
(from $s:expr, $t:ty) => {
<$t as $crate::input::Readable>::read(&mut $s)
};
(from $s:expr, $($t:tt),* $(,)?) => {
($($crate::read_value!(from $s, $t)),*)
};
($($r:tt)*) => {
$crate::input::STDIN_SOURCE.with(|s| {
let mut s = s.borrow_mut();
$crate::read_value!(from s, $($r)*)
})
}
}
#[macro_export]
macro_rules! input {
() => {
};
($x:tt: $t:tt, $($r:tt)*) => {
let $x = $crate::read_value!($t);
$crate::input!($($r)*);
};
(mut $x:tt: $t:tt, $($r:tt)*) => {
let mut $x = $crate::read_value!($t);
$crate::input!($($r)*);
};
(from $s:expr, $x:tt, $t:tt, $($r:tt)*) => {
let $x = $crate::read_value!(from $s, $t);
$crate::input!(from $s, $($r)*);
};
(from $s:expr, mut $x:tt, $t:tt, $($r:tt)*) => {
let mut $x = $crate::read_value!(from $s, $t);
$crate::input!(from $s, $($r)*);
};
($($r:tt)*) => {
$crate::input!($($r)*,);
};
}
pub trait Readable {
type Output;
fn read(source: &mut Source) -> Self::Output;
}
impl<T: FromStr<Err = E>, E: Debug> Readable for T {
type Output = T;
fn read(source: &mut Source) -> T {
source.next_token().unwrap().parse().unwrap()
}
}
pub mod marker {
macro_rules! impl_readable {
($e:ident, $t:ty, $u:ty, $f:expr) => {
pub enum $e {}
impl $crate::input::Readable for $e {
type Output = $t;
fn read(mut source: &mut $crate::input::Source) -> $t {
$f($crate::read_value!(from source, $u))
}
}
};
}
impl_readable!(Usize1, usize, usize, |x| x - 1);
impl_readable!(Isize1, isize, isize, |x| x - 1);
impl_readable!(Chars, Vec<char>, String, |x: String| x.chars().collect());
impl_readable!(Bytes, Vec<u8>, String, |x: String| x.bytes().collect());
}
}
/// seed値を4つの初期状態値に分割するためのsplit mix 64
struct SplitMix64 {
s: u64,
}
impl SplitMix64 {
fn new(seed: u64) -> Self {
Self { s: seed }
}
fn next_u64(&mut self) -> u64 {
self.s = self.s.wrapping_add(0x9e3779b97f4a7c15);
let mut z = self.s;
z = (z ^ z >> 30).wrapping_mul(0xbf58476d1ce4e5b9);
z = (z ^ z >> 27).wrapping_mul(0x94d049bb133111eb);
z ^ z >> 31
}
}
/// Xoshiro256による乱数生成器
struct Xoshiro256 {
s: [u64; 4],
}
impl Xoshiro256 {
fn new(seed: u64) -> Self {
let mut split_mix_64 = SplitMix64::new(seed);
let mut s = [0; 4];
for si in &mut s {
*si = split_mix_64.next_u64();
}
Self { s }
}
fn next_u64(&mut self) -> u64 {
let result = (self.s[1].wrapping_mul(5)).rotate_left(7).wrapping_mul(9);
let t = self.s[1] << 17;
self.s[2] ^= self.s[0];
self.s[3] ^= self.s[1];
self.s[1] ^= self.s[2];
self.s[0] ^= self.s[3];
self.s[2] ^= t;
self.s[3] = self.s[3].rotate_left(45);
result
}
fn gen_usize(&mut self, lower: usize, upper: usize) -> usize {
assert!(lower < upper);
let count = upper - lower;
(self.next_u64() % count as u64) as usize + lower
}
fn gen_i64(&mut self, lower: i64, upper: i64) -> i64 {
assert!(lower < upper);
let count = upper - lower;
(self.next_u64() % count as u64) as i64 + lower
}
fn gen_f64(&mut self) -> f64 {
const UPPER_MASK: u64 = 0x3ff0000000000000;
const LOWER_MASK: u64 = 0xfffffffffffff;
let result = UPPER_MASK | (self.next_u64() & LOWER_MASK);
let result: f64 = unsafe { std::mem::transmute(result) };
result - 1.0
}
fn gen_bool(&mut self, prob: f64) -> bool {
self.gen_f64() < prob
}
fn fisher_yates_shuffle<T>(&mut self, items: &mut [T]) {
for i in (1..items.len()).rev() {
let j = (self.next_u64() as usize) % (i + 1);
items.swap(j, i);
}
}
}
macro_rules! impl_kmeans {
($kind: ty, $modname: ident) => {
// Since we can't overload methods in rust, we have to use namespacing
pub mod $modname {
use std::$modname::MAX;
/// computes sum of squared deviation between two identically sized vectors
/// `x`, and `y`.
fn distance(x: &[$kind], y: &[$kind]) -> $kind {
x.iter()
.zip(y.iter())
.fold(0, |dist, (&xi, &yi)| dist + (xi - yi) * (xi - yi))
}
/// Returns a vector containing the indices z<sub>i</sub> in {0, ..., K-1} of
/// the centroid nearest to each datum.
fn nearest_centroids(xs: &[Vec<$kind>], centroids: &[Vec<$kind>]) -> Vec<usize> {
xs.iter()
.map(|xi| {
// Find the argmin by folding using a tuple containing the argmin
// and the minimum distance.
let (argmin, _) = centroids.iter().enumerate().fold(
(0_usize, MAX),
|(min_ix, min_dist), (ix, ci)| {
let dist = distance(xi, ci);
if dist < min_dist {
(ix, dist)
} else {
(min_ix, min_dist)
}
},
);
argmin
})
.collect()
}
/// Recompute the centroids given the current clustering
fn recompute_centroids(
xs: &[Vec<$kind>],
clustering: &[usize],
k: usize,
) -> Vec<Vec<$kind>> {
let ndims = xs[0].len();
// NOTE: Kind of inefficient because we sweep all the data from each of the
// k centroids.
(0..k)
.map(|cluster_ix| {
let mut centroid: Vec<$kind> = vec![0; ndims];
let mut n_cluster: $kind = 0;
xs.iter().zip(clustering.iter()).for_each(|(xi, &zi)| {
if zi == cluster_ix {
n_cluster += 1;
xi.iter().enumerate().for_each(|(j, &x_ij)| {
centroid[j] += x_ij;
});
}
});
centroid.iter().map(|&c_j| c_j / n_cluster).collect()
})
.collect()
}
/// Assign the N D-dimensional data, `xs`, to `k` clusters using K-Means clustering
pub fn kmeans(xs: Vec<Vec<$kind>>, k: usize) -> (Vec<usize>, Vec<Vec<$kind>>) {
assert!(xs.len() >= k);
// Rather than pulling in a dependency to radomly select the staring
// points for the centroids, we're going to deterministally choose them by
// slecting evenly spaced points in `xs`
let n_per_cluster: usize = xs.len() / k;
let centroids: Vec<Vec<$kind>> =
(0..k).map(|j| xs[j * n_per_cluster].clone()).collect();
let mut clustering = nearest_centroids(&xs, ¢roids);
loop {
let centroids = recompute_centroids(&xs, &clustering, k);
let new_clustering = nearest_centroids(&xs, ¢roids);
// loop until the clustering doesn't change after the new centroids are computed
if new_clustering
.iter()
.zip(clustering.iter())
.all(|(&za, &zb)| za == zb)
{
// We need to use `return` to break out of the `loop`
return (clustering, centroids);
} else {
clustering = new_clustering;
}
}
}
}
};
}
impl_kmeans!(i64, i64);
/*
[I 2023-06-25 10:41:46,398] Trial 2230 finished with value: 301121.87 and parameters:
{'CALC_TEMP_FREQUENCY': 2851, 'KICK_THRESHOLD': 50621, 'MAX_DELTA': 257.236083625224, 'MIN_DELTA': 14.080818963193781,
'MULTI_START': 1, 'NUM_STATION': 5, 'PROB_2OPT': 64, 'PROB_3OPT': 66, 'PROB_DOUBLE_BRIDGE': 0.756448734745971, 'PROB_INSERT': 77,
'PROB_INSERT_STATION': 83, 'PROB_MOVE_STATION': 81, 'T1': 164.0285193747435, 'TO': 14096.678262676927}.
Best is trial 2230 with value: 301121.87.
*/
const N: usize = 100;
const M: usize = 8;
const ALPHA: i64 = 5;
const MAP_SIZE: i64 = 1000;
const INF: i64 = 1_000_000_000;
const TL: f64 = 0.996;
const IS_ANNEALING: bool = true;
const T0: f64 = 14096.678262676927;
const T1: f64 = 164.0285193747435;
const MULTI_START: f64 = 1.0;
const PROB_2OPT: usize = 64;
const PROB_3OPT: usize = 66;
const PROB_INSERT: usize = 77;
const PROB_MOVE_STATION: usize = 81;
const PROB_INSERT_STATION: usize = 83;
const MAX_DELTA: f64 = 257.236083625224;
const MIN_DELTA: f64 = 14.080818963193781;
const KICK_THRESHOLD: usize = 50621;
const PROB_DOUBLE_BRIDGE: f64 = 0.756448734745971;
const CALC_TEMP_FREQUENCY: usize = 2851;
const NUM_STATION: usize = 5;
pub trait ChangeMinMax {
fn change_min(&mut self, v: Self) -> bool;
fn change_max(&mut self, v: Self) -> bool;
}
impl<T: PartialOrd> ChangeMinMax for T {
fn change_min(&mut self, v: T) -> bool {
*self > v && {
*self = v;
true
}
}
fn change_max(&mut self, v: T) -> bool {
*self < v && {
*self = v;
true
}
}
}
pub fn get_time() -> f64 {
static mut STIME: f64 = -1.0;
let t = std::time::SystemTime::now()
.duration_since(std::time::UNIX_EPOCH)
.unwrap();
let ms = t.as_secs() as f64 + t.subsec_nanos() as f64 * 1e-9;
unsafe {
if STIME < 0.0 {
STIME = ms;
}
// ローカル環境とジャッジ環境の実行速度差はget_timeで吸収しておくと便利
#[cfg(feature = "local")]
{
(ms - STIME) * 1.5
}
#[cfg(not(feature = "local"))]
{
(ms - STIME)
}
}
}
// 頂点の種類によって距離の二乗を返す
// あまり宇宙ステーションを贔屓し過ぎないようなパラメータも考慮
pub fn squared_distance(point_i: Position, i: usize, point_j: Position, j: usize) -> i64 {
if i < N && j < N {
((point_i.x - point_j.x) * (point_i.x - point_j.x)
+ (point_i.y - point_j.y) * (point_i.y - point_j.y))
* ALPHA
* ALPHA
} else if i < N || j < N {
((point_i.x - point_j.x) * (point_i.x - point_j.x)
+ (point_i.y - point_j.y) * (point_i.y - point_j.y))
* ALPHA
} else {
((point_i.x - point_j.x) * (point_i.x - point_j.x)
+ (point_i.y - point_j.y) * (point_i.y - point_j.y))
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub struct Position {
x: i64,
y: i64,
}
impl Position {
fn new(x: i64, y: i64) -> Self {
Self { x, y }
}
}
#[derive(Clone)]
pub struct Input {
pos: Vec<Position>, // 惑星の位置
dist: Vec<Vec<i64>>, // ステーションを無視した惑星同士の距離
}
impl Input {
fn new() -> Self {
input! {
_n: usize,
_m: usize,
}
let mut pos = vec![];
// kmeansでステーションの初期位置を定めるため，posを変換
let mut pos_kmeans = vec![];
for _ in 0..N {
input! {
ai: i64,
bi: i64,
}
pos.push(Position::new(ai, bi));
pos_kmeans.push(vec![ai, bi]);
}
// ステーションの効果を無視して惑星間でワーシャルフロイド
let dist = warshall_floyd(&pos);
// kmeansでステーションの初期位置を決める
let (_, station_kmeans) = kmeans(pos_kmeans, NUM_STATION);
let mut station = vec![];
// kmeansした箇所
for i in 0..NUM_STATION {
station.push(Position::new(station_kmeans[i][0], station_kmeans[i][1]));
}
// kmeansしてない箇所はランダムに入れる
let mut rng = Xoshiro256::new(16384);
for i in NUM_STATION..M {
station.push(Position::new(rng.gen_i64(0, 1000), rng.gen_i64(0, 1000)));
}
pos.extend(station);
Self { pos, dist }
}
}
// ワーシャルフロイドで全点対間最短経路を求める
pub fn warshall_floyd(pos: &Vec<Position>) -> Vec<Vec<i64>> {
let mut dist = vec![vec![0; pos.len()]; pos.len()];
// (i,j)間の距離を入れる
for i in 0..pos.len() {
for j in i + 1..pos.len() {
dist[i][j] = squared_distance(pos[i], i, pos[j], j);
dist[j][i] = dist[i][j];
}
}
// ワーシャルフロイド法
for k in 0..pos.len() {
for i in 0..pos.len() {
for j in 0..pos.len() {
// (i,j)間はkを経由したほうが短くなるか調べる
if dist[i][j] > dist[i][k] + dist[k][j] {
dist[i][j] = dist[i][k] + dist[k][j];
}
}
}
}
dist
}
// (i, j)間の最短パスをダイクストラで見つける
pub fn dijkstra(pos: &Vec<Position>, s: usize, t: usize) -> Vec<usize> {
// s->tパスをダイクストラで見つける
let mut dist = vec![INF; N + M];
dist[s] = 0;
let mut prev = vec![!0; N + M]; // 経路復元用
let mut bh = BinaryHeap::new();
bh.push((Reverse(0), s));
while let Some((Reverse(d), frm)) = bh.pop() {
// 目的地についているなら，終了する
if frm == t {
break;
}
// 見ようとしているものより既にいい経路が見つかっていれば飛ばす
if dist[frm] < d {
continue;
}
for to in 0..N + M {
let energy = squared_distance(pos[frm], frm, pos[to], to);
if d + energy < dist[to] {
// コストを更新したほうがいいなら，更新しつつ優先度付きキューに入れる
dist[to] = d + energy;
prev[to] = frm; // 頂点 frm -> toにたどり着いた
bh.push((Reverse(dist[to]), to));
}
}
}
// 経路復元
let mut path = vec![];
let mut tt = t;
while tt != s {
// toがstartになるまでprev[to]を辿っていく
// 途中のエネルギーを保存しておく
path.push(tt);
tt = prev[tt];
}
path.reverse();
path
}
// 惑星1から開始して，未訪問の惑星で一番近いものへのパスを経路復元つきダイクストラで求めていく
pub fn init_path(input: &Input) -> Vec<usize> {
let mut path = vec![];
let mut visited = vec![false; N];
let mut s = 0;
visited[s] = true;
for i in 0..N {
// sから未訪問の惑星のうち，一番近いものを見つける
let mut min_dist = INF;
let mut t = !0;
if i < N - 1 {
// 目的地以外の最大N-1個の探索
for j in 0..N {
if visited[j] || j == 0 {
continue;
}
if input.dist[s][j] < min_dist {
min_dist = input.dist[s][j];
t = j;
}
}
// ここまで来てもt = !0のままなら，すべて訪問済みなのでcontinueしてよい
if t == !0 {
continue;
}
} else {
// 最後は目的地に帰ってきたいので，visitedを無視して設定
t = 0;
}
// s->tパスをダイクストラで見つけ，全体に追加
let mut small_path = dijkstra(&input.pos, s, t);
path.extend(small_path);
// 次のイテレーションのための更新
s = t;
visited[t] = true;
}
path
}
#[derive(Clone)]
struct State {
pos: Vec<Position>, // 先頭N個に惑星の場所，続くM個にステーションの場所が格納されている
path: Vec<usize>, // 道順 (惑星orステーション，番号)
}
impl State {
fn new(input: &Input, greedy_path: &Vec<usize>) -> Self {
Self {
pos: input.pos.clone(),
path: greedy_path.clone(),
}
}
fn calc_score_all(&self, input: &Input) -> i64 {
let mut score = 0;
for w in self.path.windows(2) {
let (prev, next) = (w[0], w[1]);
score += self.calc_score(input, prev, next);
}
score
}
#[inline]
fn calc_score(&self, input: &Input, prev: usize, next: usize) -> i64 {
if prev < N && next < N {
input.dist[prev][next]
} else {
squared_distance(self.pos[prev], prev, self.pos[next], next)
}
}
// 1-indexに直しつつ回答を出力
pub fn print(&self) {
for i in 0..M {
println!("{} {}", self.pos[i + N].x, self.pos[i + N].y);
}
let v = self.path.len();
println!("{}", v);
for &i in self.path.iter() {
if i < N {
println!("1 {}", i + 1);
} else {
println!("2 {}", i + 1 - N);
}
}
}
}
fn restore_solution(input: &Input, solution: &State) -> State {
let mut new_solution = solution.clone();
new_solution.path.clear();
new_solution.path.push(0);
for w in solution.path.windows(2) {
let (prev, next) = (w[0], w[1]);
let path = dijkstra(&solution.pos, prev, next);
for v in path {
new_solution.path.push(v);
}
}
new_solution
}
fn main() {
get_time();
let input = Input::new();
let mut total_iter = 0;
let mut each_iter = 0.; // 多点スタートの各イテレーション
let mut accepted_count = 0;
let mut update_best_count = 0;
let mut double_bridge_count = 0;
let mut rin_sakamichi_count = 0;
// 初期解を貪欲に決定
let mut greedy_path = vec![0]; // このdijkstraの経路復元はスタート位置を保持しないのでここで代入
greedy_path.extend(init_path(&input));
let mut best_state = State::new(&input, &greedy_path);
let mut best_score = best_state.calc_score_all(&input);
while get_time() < TL {
// 多点スタートで振り出しに戻る箇所
each_iter += 1.;
// 乱数シードも変更する
let mut rng = Xoshiro256::new(8192 + each_iter as u64);
// stateを初期のものに変更する
let mut state = State::new(&input, &greedy_path);
let mut current_score = state.calc_score_all(&input);
// 解を更新しなかった回数(閾値を超えたらdouble-bridgeでkickする)
let mut not_accepted_count = 0;
// 温度
let mut t = get_time() / (TL * each_iter / MULTI_START);
let mut T = T0.powf(1.0 - t) * T1.powf(t);
// 特定の近傍の次は必ず2-optするようにする
let mut is_select_2opt = false;
// 各多点のイテレーションは，TL / MULTI_START秒ずつ割り当てられる
while get_time() < TL * each_iter / MULTI_START {
total_iter += 1;
if total_iter % CALC_TEMP_FREQUENCY == 0 {
t = get_time() / (TL * each_iter / MULTI_START);
T = T0.powf(1.0 - t) * T1.powf(t);
}
// 乱数で生じさせた数値に応じて近傍を選択
let neighbor_type = rng.gen_usize(0, 100);
if is_select_2opt || neighbor_type < PROB_2OPT || get_time() > (TL * each_iter / MULTI_START) * 0.95 {
// 2-opt
// 最後の方は交差路をなくすために2-optばかりする
// ある(u,next_u), (v,next_v)に対し，(u,v), (next_u, next_v) と繋ぎ変えたほうが嬉しいのか調べる
let mut u = !0;
let mut v = !0;
while u == v {
u = rng.gen_usize(0, state.path.len() - 1);
v = rng.gen_usize(0, state.path.len() - 1);
}
if u > v {std::mem::swap(&mut u, &mut v);}
let i0 = state.path[u];
let i1 = state.path[u + 1];
let i2 = state.path[v];
let i3 = state.path[v + 1];
let d01 = state.calc_score(&input, i0, i1);
let d23 = state.calc_score(&input, i2, i3);
let d02 = state.calc_score(&input, i0, i2);
let d13 = state.calc_score(&input, i1, i3);
// スコア計算
let score_diff = d02 + d13 - d01 - d23;
let new_score = current_score + score_diff;
if score_diff < 0 || IS_ANNEALING && rng.gen_bool(f64::exp(-score_diff as f64 / T))
{
// 解の更新
current_score = new_score;
accepted_count += 1;
state.path[u + 1..v + 1].reverse();
is_select_2opt = false;
} else {
not_accepted_count += 1;
}
} else if neighbor_type < PROB_3OPT {
// 3-opt
// 20230119 梅谷 メタヒューリスティクスの設計と実装 を参照
// (u,next_u),(v,next_v),(w,next_w)のつなぎ替えのうちで最もいいものを採用する
// type1(Or-opt): (u,next_v),(w,next_u),(v,next_w)
// type2: (u,w),(next_v,next_u),(v,next_w)
// type3: (u,next_v),(w,v),(next_u,next_w)
// type4: (u,v),(next_u,w),(next_v,next_w)
// 3点の選び方
// (0..state.path.len()-1)をフィッシャー･イェーツ シャッフルし，先頭3点を選ぶ
let mut points: Vec<usize> = (0..state.path.len()-1).collect();
rng.fisher_yates_shuffle(&mut points);
let (first_three, _) = points.split_at_mut(3);
first_three.sort();
let (u, v, w) = (first_three[0], first_three[1], first_three[2]);
let i0 = state.path[u];
let i1 = state.path[u + 1];
let i2 = state.path[v];
let i3 = state.path[v + 1];
let i4 = state.path[w];
let i5 = state.path[w + 1];
let mut original_score = state.calc_score(&input, i0, i1);
original_score += state.calc_score(&input, i2, i3);
original_score += state.calc_score(&input, i4, i5);
// type1
// u -> v+1
let mut type1_score = state.calc_score(&input, i0, i3);
// w -> u+1
type1_score += state.calc_score(&input, i4, i1);
// v -> w+1
type1_score += state.calc_score(&input, i2, i5);
// type2
// u -> w
let mut type2_score = state.calc_score(&input, i0, i4);
// v+1 -> u+1
type2_score += state.calc_score(&input, i3, i1);
// v -> w+1
type2_score += state.calc_score(&input, i2, i5);
// type3
// u -> v+1
let mut type3_score = state.calc_score(&input, i0, i3);
// w -> v
type3_score += state.calc_score(&input, i4, i2);
// u+1 -> w+1
type3_score += state.calc_score(&input, i1, i5);
// type4
// u -> v
let mut type4_score = state.calc_score(&input, i0, i2);
// u+1 -> w
type4_score += state.calc_score(&input, i1, i4);
// v+1 -> w+1
type4_score += state.calc_score(&input, i3, i5);
// スコア計算
let each_type_score = vec![type1_score, type2_score, type3_score, type4_score];
let mut score_diff = INF;
let mut best_type = !0;
for (idx, &type_score) in each_type_score.iter().enumerate() {
if score_diff.change_min(type_score - original_score) {
best_type = idx;
}
}
let new_score = current_score + score_diff;
if score_diff < 0 || IS_ANNEALING && rng.gen_bool(f64::exp(-score_diff as f64 / T))
{
// 解の更新
current_score = new_score;
accepted_count += 1;
// 採用された交換のtypeによって，pathの変え方を分岐する
let part1 = state.path[u + 1..v + 1].to_vec();
let part2 = state.path[v + 1..w + 1].to_vec();
if best_type == 0 {
// path[0..u,u+1..v,v+1..w,w+1..0]をpath[0..u,v+1..w,u+1..v,w+1..0]としたい
state.path.splice(
u + 1..w + 1,
part2.iter().cloned().chain(part1.iter().cloned()),
);
} else if best_type == 1 {
state.path.splice(
u + 1..w + 1,
part2.iter().rev().cloned().chain(part1.iter().cloned()),
);
} else if best_type == 2 {
state.path.splice(
u + 1..w + 1,
part2.iter().cloned().chain(part1.iter().rev().cloned()),
);
} else {
state.path.splice(
u + 1..w + 1,
part1
.iter()
.rev()
.cloned()
.chain(part2.iter().rev().cloned()),
);
}
is_select_2opt = true;
} else {
not_accepted_count += 1;
}
} else if neighbor_type < PROB_INSERT {
// 挿入近傍
// あるu, vを選び，uをvの後ろへ入れる
let mut u = !0;
let mut v = !0;
while u == v {
u = rng.gen_usize(1, state.path.len() - 2);
v = rng.gen_usize(1, state.path.len() - 2);
}
if u > v {std::mem::swap(&mut u, &mut v);}
let i0 = state.path[u - 1];
let i1 = state.path[u];
let i2 = state.path[u + 1];
let i3 = state.path[v];
let i4 = state.path[v + 1];
let mut original_score = state.calc_score(&input, i0, i1);
original_score += state.calc_score(&input, i1, i2);
original_score += state.calc_score(&input, i3, i4);
// u-1 -> u+1
let mut insert_score = state.calc_score(&input, i0, i2);
// v -> u
insert_score += state.calc_score(&input, i3, i1);
// u -> v+1
insert_score += state.calc_score(&input, i1, i4);
// スコア計算
let score_diff = insert_score - original_score;
let new_score = current_score + score_diff;
if score_diff < 0 || IS_ANNEALING && rng.gen_bool(f64::exp(-score_diff as f64 / T))
{
// 解の更新
current_score = new_score;
accepted_count += 1;
// path[0..u-1,u,u+1..v,v+1..0]をpath[0..u-1,u+1..v,u,v+1..0]としたい
let tmp = state.path[u];
state.path.remove(u);
state.path.insert(v, tmp);
is_select_2opt = true;
} else {
not_accepted_count += 1;
}
} else if neighbor_type < PROB_MOVE_STATION {
// あるstationを一旦削除し、ランダムにずらした上で、各辺でstationを使う/使わないを貪欲に決め直す
let station_id = N + rng.gen_usize(0, M);
let mut tmp_path = state.path.clone();
tmp_path.retain(|&v| v != station_id);
let old_p = state.pos[station_id];
let mut p = old_p;
let delta = (MAX_DELTA * (1. - t) + MIN_DELTA * t) as i64;
p.x = rng.gen_i64((p.x - delta).max(0), (p.x + delta).min(MAP_SIZE) + 1);
p.y = rng.gen_i64((p.y - delta).max(0), (p.y + delta).min(MAP_SIZE) + 1);
state.pos[station_id] = p;
let mut new_path = vec![];
let mut new_score = 0;
// path上の点で，新たにずらしたステーションを使った方がよいかどうかを判断する
// warshall floydに近い気持ち
for w in tmp_path.windows(2) {
let (prev, next) = (w[0], w[1]);
new_path.push(prev);
let old_dist = state.calc_score(&input, prev, next);
let new_dist = state.calc_score(&input, prev, station_id)
+ state.calc_score(&input, station_id, next);
if new_dist < old_dist {
new_score += new_dist;
new_path.push(station_id);
} else {
new_score += old_dist;
}
}
// 惑星1に帰ってくる
new_path.push(0);
let score_diff = new_score - current_score;
if score_diff < 0 || IS_ANNEALING && rng.gen_bool(f64::exp(-score_diff as f64 / T))
{
// 解の更新
current_score = new_score;
accepted_count += 1;
state.path = new_path;
} else {
// ロールバック
state.pos[station_id] = old_p;
not_accepted_count += 1;
}
} else if neighbor_type < PROB_INSERT_STATION {
// stationを適当な位置に挿入する
let station_id = N + rng.gen_usize(0, M);
let index = rng.gen_usize(1, state.path.len());
let prev = state.path[index - 1];
let next = state.path[index];
let old_score = state.calc_score(&input, prev, next);
let new_score = state.calc_score(&input, prev, station_id)
+ state.calc_score(&input, station_id, next);
let score_diff = new_score - old_score;
if score_diff <= 0 || rng.gen_bool(f64::exp(-score_diff as f64 / T)) {
// 解の更新
current_score += score_diff;
accepted_count += 1;
state.path.insert(index, station_id);
} else {
not_accepted_count += 1;
}
} else {
// 適当な位置のstationを削除する
let mut index = 0;
let mut trial = 0;
while state.path[index] < N && trial < 10 {
index = rng.gen_usize(0, state.path.len());
trial += 1;
}
if state.path[index] < N {
continue;
}
let station_id = state.path[index];
let prev = state.path[index - 1];
let next: usize = state.path[index + 1];
let old_score = state.calc_score(&input, prev, station_id)
+ state.calc_score(&input, station_id, next);
let new_score = state.calc_score(&input, prev, next);
let score_diff = new_score - old_score;
if score_diff <= 0 || rng.gen_bool(f64::exp(-score_diff as f64 / T)) {
// 解の更新
current_score += score_diff;
accepted_count += 1;
state.path.remove(index);
} else {
not_accepted_count += 1;
}
}
// どの近傍でもベスト更新は常に採用する
if current_score < best_score {
best_score = current_score;
best_state = state.clone();
update_best_count += 1;
// eprintln!("update best: {}", best_score);
// eprintln!("update iter: {}", iter);
}
// 近傍遷移回数が一定回行われなければkickする
if not_accepted_count == KICK_THRESHOLD {
not_accepted_count = 0;
if rng.gen_bool(PROB_DOUBLE_BRIDGE) {
// double-bridge (4-optの1つ)
double_bridge_count += 1;
// (u,next_u),(v,next_v),(w,next_w),(x,next_x)を以下のように繋ぎ変える
// (u,next_w),(x,next_v),(w,next_u),(v,next_x)
let mut points: Vec<usize> = (0..state.path.len()-1).collect();
rng.fisher_yates_shuffle(&mut points);
let (first_four, _) = points.split_at_mut(4);
first_four.sort();
let (u, v, w, x) = (first_four[0], first_four[1], first_four[2], first_four[3]);
let i0 = state.path[u];
let i1 = state.path[u + 1];
let i2 = state.path[v];
let i3 = state.path[v + 1];
let i4 = state.path[w];
let i5 = state.path[w + 1];
let i6 = state.path[x];
let i7 = state.path[x + 1];
let mut original_score = state.calc_score(&input, i0, i1);
original_score += state.calc_score(&input, i2, i3);
original_score += state.calc_score(&input, i4, i5);
original_score += state.calc_score(&input, i6, i7);
// double-bridge
// u -> w+1
let mut kick_score = state.calc_score(&input, i0, i5);
// x -> v+1
kick_score += state.calc_score(&input, i6, i3);
// w -> u+1
kick_score += state.calc_score(&input, i4, i1);
// v -> x+1
kick_score += state.calc_score(&input, i2, i7);
let new_score = current_score + (kick_score - original_score);
// kick近傍は必ず更新
current_score = new_score;
// 採用された交換のtypeによって，pathの変え方を分岐する
let part1 = state.path[u + 1..v + 1].to_vec();
let part2 = state.path[v + 1..w + 1].to_vec();
let part3 = state.path[w + 1..x + 1].to_vec();
// path[0..u,u+1..v,v+1..w,w+1..x,x+1..0]をpath[0..u,w+1..x,v+1..w,u+1..v,x+1..0]としたい
state.path.splice(
u + 1..x + 1,
part3
.iter()
.cloned()
.chain(part2.iter().cloned().chain(part1.iter().cloned())),
);
} else {
// 坂道輪サーチ
rin_sakamichi_count += 1;
state = best_state.clone();
current_score = best_score;
}
is_select_2opt = true;
}
}
eprintln!("=== annealing ===");
eprintln!("total iter: {}", total_iter);
eprintln!("each iter: {}", each_iter);
eprintln!("accepted: {}", accepted_count);
eprintln!("update best: {}", update_best_count);
eprintln!("kicked by double bridge: {}", double_bridge_count);
eprintln!("kicked by rin sakamichi search: {}", rin_sakamichi_count);
eprintln!("energy: {}", best_score);
eprintln!(
"score: {}",
(1e9 / (1e3 + (best_score as f64).sqrt())).round() as i64
);
eprintln!("time: {}", get_time());
eprintln!();
}
let solution = restore_solution(&input, &best_state);
solution.print();
let final_score = solution.calc_score_all(&input);
eprintln!(
"final score: {}",
(1e9 / (1e3 + (final_score as f64).sqrt())).round() as i64
);
eprintln!("time: {}", get_time());
}