#863522 (Rust) No.5007 Steiner Space Travel

• 問題: No.5007 Steiner Space Travel
• ユーザー: r3yohei
• 提出日時: 2023-04-28 18:58:58
• 言語: Rust (1.77.0 + proconio)
• 結果: AC
• 実行時間: 982 ms / 1,000 ms
• コード長: 26,633 bytes
• コンパイル時間: 7,732 ms
• コンパイル使用メモリ: 174,596 KB
• 実行使用メモリ: 4,372 KB
• スコア: 8,240,065
• 最終ジャッジ日時: 2023-04-28 18:59:39
• 合計ジャッジ時間: 36,808 ms
• ジャッジサーバーID （参考情報）: judge11 / judge15

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	981 ms; 4,372 KB
testcase_01	AC	982 ms; 4,368 KB
testcase_02	AC	982 ms; 4,372 KB
testcase_03	AC	982 ms; 4,372 KB
testcase_04	AC	982 ms; 4,372 KB
testcase_05	AC	981 ms; 4,368 KB
testcase_06	AC	982 ms; 4,368 KB
testcase_07	AC	982 ms; 4,368 KB
testcase_08	AC	982 ms; 4,372 KB
testcase_09	AC	982 ms; 4,372 KB
testcase_10	AC	981 ms; 4,368 KB
testcase_11	AC	981 ms; 4,368 KB
testcase_12	AC	981 ms; 4,372 KB
testcase_13	AC	982 ms; 4,372 KB
testcase_14	AC	981 ms; 4,368 KB
testcase_15	AC	982 ms; 4,368 KB
testcase_16	AC	982 ms; 4,368 KB
testcase_17	AC	981 ms; 4,372 KB
testcase_18	AC	982 ms; 4,372 KB
testcase_19	AC	982 ms; 4,372 KB
testcase_20	AC	982 ms; 4,368 KB
testcase_21	AC	981 ms; 4,372 KB
testcase_22	AC	981 ms; 4,368 KB
testcase_23	AC	982 ms; 4,368 KB
testcase_24	AC	982 ms; 4,368 KB
testcase_25	AC	982 ms; 4,372 KB
testcase_26	AC	982 ms; 4,368 KB
testcase_27	AC	982 ms; 4,368 KB
testcase_28	AC	982 ms; 4,368 KB
testcase_29	AC	982 ms; 4,372 KB

ソースコード

#![allow(non_snake_case, unused)]

use std::{io::*, hash};
use std::{collections::*, fmt::format};
use std::{cmp::*, vec};
use std::time::SystemTime;
use std::collections::hash_map::DefaultHasher;
use std::hash::{Hash, Hasher};
use std::ops::{Add, Sub, Rem, Mul};
use crate::input::{*, marker::*};
use self::i64::kmeans;

// コピペで使える proconio もどき
// cunitacさんからお借りしています
// https://gist.github.com/cunitac/b00be62bf68c9fb6055d22eb77c17e14
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());
    }
}

// yukicoder用乱数生成器
// 現在時刻のハッシュ値を適当な範囲に変換する
// 確率値に応じてboolの値を返す
fn rand_gen_bool(probability: f64) -> bool {
    assert!((0.0..=1.0).contains(&probability), "Probability must be between 0 and 1");
    let mut rng = Xorshift::new();
    let random_float = rng.rand_f64();
    random_float <= probability
}

// 半開区間 [min, max) の範囲の乱数を得る
fn rand_gen_range<T>(min: T, max: T) -> T
where
    T: UnsignedInteger,
{
    let mut rng = Xorshift::new();
    T::generate_in_range(&mut rng, min, max)
}

trait UnsignedInteger:
    Copy + PartialEq + PartialOrd + Add<Output = Self> + Sub<Output = Self> + Rem<Output = Self> + Mul<Output = Self>
{
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self;
}

impl UnsignedInteger for usize {
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self {
        (rng.rand_u64() as usize % (max - min)) + min
    }
}

impl UnsignedInteger for u8 {
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self {
        (rng.rand_u64() as u8 % (max - min)) + min
    }
}

impl UnsignedInteger for u16 {
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self {
        (rng.rand_u64() as u16 % (max - min)) + min
    }
}

impl UnsignedInteger for u32 {
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self {
        (rng.rand_u64() as u32 % (max - min)) + min
    }
}

impl UnsignedInteger for u64 {
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self {
        (rng.rand_u64() % (max - min)) + min
    }
}

impl UnsignedInteger for u128 {
    fn generate_in_range(rng: &mut Xorshift, min: Self, max: Self) -> Self {
        (rng.rand_u64() as u128 % (max - min)) + min
    }
}

struct Xorshift {
    state: u64,
}

impl Xorshift {
    fn new() -> Self {
        let seed = SystemTime::now()
            .duration_since(SystemTime::UNIX_EPOCH)
            .expect("Failed to get the duration since UNIX_EPOCH")
            .as_micros() as u64;

        Xorshift { state: seed }
    }

    fn rand_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }
    
    fn rand_f64(&mut self) -> f64 {
        let u64_max_as_f64 = u64::MAX as f64;
        self.rand_u64() as f64 / u64_max_as_f64
    }
}

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, &centroids);

                loop {
                    let centroids = recompute_centroids(&xs, &clustering, k);
                    let new_clustering = nearest_centroids(&xs, &centroids);

                    // 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);

const N: usize = 100;
const M: usize = 8;
const ALPHA: i64 = 5;
const INF: i64 = 1_000_000_000;
const TL: f64 = 0.98;
// const TL: f64 = 10.0;
// この問題は解の改善幅が10^3オーダーくらい
// 仮に1000悪くなりT=1000のとき，e^(-1) = 1/2.7くらいの確率で採用される
// T=500ならe^(-2) = 1/(2.7)^2
const T0: f64 = 10000.; // 焼きなまし初期温度
const T1: f64 = 20.; // 焼きなまし終温度

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(xi: i64, y1: i64, i: usize, xj: i64, y2: i64, j: usize, discount: i64) -> i64 {
    if i < N && j < N {
        ((xi - xj) * (xi - xj) + (y1 - y2) * (y1 - y2)) * ALPHA * ALPHA / discount / discount
    } else if i < N || j < N {
        ((xi - xj) * (xi - xj) + (y1 - y2) * (y1 - y2)) * ALPHA / discount
    } else {
        ((xi - xj) * (xi - xj) + (y1 - y2) * (y1 - y2))
    }
}

#[derive(Clone)]
pub struct Input {
    pos: Vec<Vec<i64>>,
}
impl Input {
    fn new() -> Self {
        input! {
            _n: usize,
            _m: usize,
        }
        let mut pos = vec![];
        for _ in 0..N {
            input! {
                ai: i64,
                bi: i64,
            }
            pos.push(vec![ai, bi]);
        }
        Self { pos }
    }
}

// ワーシャルフロイドで全点対間最短経路を求める
pub fn warshall_floyd(pos: &Vec<Vec<i64>>) -> Vec<Vec<i64>> {
    let mut dist = vec![vec![INF; N+M]; N+M];
    // 自身への距離は0
    for i in 0..N+M {
        dist[i][i] = 0;    
    }
    // (i,j)間の距離を入れる
    for i in 0..N+M {
        for j in i+1..N+M {
            dist[i][j] = squared_distance(pos[i][0], pos[i][1], i, pos[j][0], pos[j][1], j, 7);
            dist[j][i] = dist[i][j];
        }
    }
    // ワーシャルフロイド法
    for k in 0..N+M {
        for i in 0..N+M {
            for j in 0..N+M {
                // (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<Vec<i64>>, s: usize, t: usize) -> (Vec<(usize, usize)>, Vec<i64>, i64) {
    // 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][0], pos[frm][1], frm, pos[to][0], pos[to][1], to, 1);
            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 energy = vec![];
    let mut sum_energy = 0;
    let mut tt = t;
    while tt != !0 {
        // toがstartになるまでprev[to]を辿っていく
        // 途中のエネルギーを保存しておく
        if !path.is_empty() {
            let &(last_type, mut last_index) = path.last().unwrap();
            // pathにはちゃんとした回答用のものが入っているので，惑星かステーションかによってindexをずらすかどうか見ないといけない
            last_index = if last_type == 1 {
                // よくわからないが型推論できないと怒られる
                last_index as usize
            } else {
                last_index + N
            };
            // 実際のスコアが必要な個所では割引かない
            let tmp_energy = squared_distance(pos[last_index][0], pos[last_index][1], last_index, pos[tt][0], pos[tt][1], tt, 1);
            energy.push(tmp_energy);
            sum_energy += tmp_energy;
        }
        if tt < N {
            // 惑星
            path.push((1, tt));
        } else {
            // ステーション
            path.push((2, tt - N));
        }
        tt = prev[tt];
    }
    path.reverse();
    energy.reverse();
    (path, energy, sum_energy)
}

// 始点sから開始して，未訪問の惑星で一番近いものへのパスを経路復元つきダイクストラで求めていく
pub fn update_path(mut s: usize, destination: usize, pos: &Vec<Vec<i64>>, dist_warshall_floyd: &Vec<Vec<i64>>, mut visited: &mut Vec<bool>) -> (Vec<(usize, usize)>, Vec<i64>, i64) {
    let mut path = vec![];
    let mut energy = vec![];
    let mut sum_energy = 0;
    // スタートがステーションである可能性を考慮
    if s < N {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 s == j || visited[j] || j == destination {continue;}
                if dist_warshall_floyd[s][j] < min_dist {
                    min_dist = dist_warshall_floyd[s][j];
                    t = j;
                }
            }
            // ここまで来てもt = !0のままなら，すべて訪問済みなのでcontinueしてよい
            if t == !0 {continue;}
        } else {
            // 最後は目的地に帰ってきたいので，visitedを無視して設定
            t = destination;
        }
        // eprintln!("s {}, t {}", s, t);
        // s->tパスをダイクストラで見つけ，全体に追加
        let (mut small_path, small_energy, small_sum_energy) = dijkstra(&pos, s, t);
        if t != destination {
            // 目的地以外では，tが次のsになるので，pathとしては含めなくてよい
            small_path.pop();
        }
        path.extend(small_path);
        energy.extend(small_energy);
        sum_energy += small_sum_energy;
        // 次のイテレーションのための更新
        s = t;
        // ゴールがステーションである可能性を考慮
        if t < N {visited[t] = true;}
    }
    (path, energy, sum_energy)
}

#[derive(Clone)]
struct State {
    turn: i64,
    is_done: bool,
    game_score: i64,
    evaluated_score: i64,
    pos: Vec<Vec<i64>>, // 先頭N個に惑星の場所，続くM個にステーションの場所が格納されている
    dist_warshall_floyd: Vec<Vec<i64>>, // 惑星と宇宙ステーション間の距離
    path: Vec<(usize, usize)>, // 道順 (惑星orステーション，番号)
    energy: Vec<i64>, // 上記パスを通った場合に消費するエネルギー
    sum_energy: i64, // 上記の合計
}
impl State {
    fn new(input: &Input) -> Self {
        // 各惑星をkmeansでクラスタリングする
        // どのクラスタに所属するかと，各クラスタのセントロイドを保持し，セントロイドは宇宙ステーションとする
        let (_, mut station) = kmeans(input.pos.clone(), 8);
        let mut pos = input.pos.clone();
        pos.extend(station);
        // 惑星と宇宙ステーションの位置をもとに，(i, j)間の距離をワーシャルフロイドで求めておく
        let dist_warshall_floyd = warshall_floyd(&pos);
        let mut visited = vec![false; N];
        let (path, tmp_energy, sum_energy) = update_path(0, 0, &pos, &dist_warshall_floyd, &mut visited);
        // pathとenergyの関係は，星-(エネルギー)-星のように，階差の関係であるからpath.len() - energy.len() == 1となる
        // これだと扱いがめんどうなので，energeyの先頭にdummyを入れておく
        let mut energy = vec![0];
        energy.extend(tmp_energy);
        Self {
            turn: 0,
            is_done: false,
            game_score: 0,
            evaluated_score: 0,
            pos,
            dist_warshall_floyd,
            path,
            energy,
            sum_energy,
        }
    }
    // 任意のステーションの位置を少しずらす
    pub fn shift_station(&mut self) {
        let target_station = rand_gen_range(0, M);
        let diff = 25_u64;
        let mut diff_x = rand_gen_range(1_u64, diff) as i64;
        if rand_gen_range(1_u8, 3_u8) == 1 {
            diff_x *= -1;
        }
        let mut diff_y = rand_gen_range(1_u64, diff) as i64;
        if rand_gen_range(1_u8, 3_u8) == 1 {
            diff_y *= -1;
        }
        // 0~1000の範囲を守ってずらす
        self.pos[N + target_station][0] = min(max(0, self.pos[N + target_station][0] + diff_x), 1000);
        self.pos[N + target_station][1] = min(max(0, self.pos[N + target_station][1] + diff_y), 1000);
        // 惑星と宇宙ステーションの位置をもとに，(i, j)間の距離をワーシャルフロイドで求めておく
        let dist_warshall_floyd = warshall_floyd(&self.pos);
        let mut visited = vec![false; N];
        let (path, energy, sum_energy) = update_path(0, 0, &self.pos, &dist_warshall_floyd, &mut visited);
        self.dist_warshall_floyd = dist_warshall_floyd;
        self.path = path;
        self.energy = energy;
        self.sum_energy = sum_energy;
    }
    // パス上の適当な両端点p1, p2を選択肢，p1~p2の経路を部分破壊&再構築する
    pub fn scrap_and_build(&mut self) {
        let mut tmp_energy = self.sum_energy;
        // 適当な惑星1以外のp1,p2を選ぶ
        let width = 20;
        let mut p1 = rand_gen_range(1, self.path.len() - 2 - width);
        let mut p2 = rand_gen_range(p1 + 1, min(p1 + width, self.path.len() - 1));
        // eprintln!("p1 {:?}, p2 {:?}", self.path[p1], self.path[p2]);
        // 惑星1~p1, p1~p2, p2~惑星1に分割する
        let path_s_p1 = &self.path[..p1];
        let crt_path_p1_p2 = &self.path[p1..=p2];
        // eprintln!("crt_path_p1_p2: {:?}", crt_path_p1_p2);
        let path_p2_t = &self.path[p2+1..];
        // エネルギーの方は，p1に至るまでに必要なエネルギーは保持しておきたいので，等号あり
        let energy_s_p1 = &self.energy[..=p1];
        // p1~p2に至るためのエネルギーは全て破棄してよい
        let crt_energy_p1_p2 = &self.energy[p1+1..=p2];
        let energy_p2_t = &self.energy[p2+1..];
        // 現在のp1~p2のエネルギーを調べて，引く
        for &e in crt_energy_p1_p2 {
            tmp_energy -= e;
        }
        // p1~p2の最短経路をダイクストラで再構築
        // p1~p2の経路内の惑星だけを未訪問とする
        // 厳密にp1~p2経路上の惑星だけを未訪問にするのは，経路が更新されづらく，厳しい．
        let mut visited = vec![true; N];
        for &(p_type, p_index) in crt_path_p1_p2 {
            if p_type == 1 {
                visited[p_index] = false;
                // p_indexだけでなく，p_indexの近くの惑星も未訪問にする
                // [ToDo] 前計算する
                // for i in 0..N {
                //     if !visited[i] {continue;}
                //     if squared_distance(self.pos[p_index][0], self.pos[p_index][1], p_index, self.pos[i][0], self.pos[i][1], i, ALPHA) <= 10000 {
                //         visited[i] = false
                //     }
                // }
            }
        }
        // p1,p2がステーションなら，インデックスにNを足さなければならない
        let p1_index = if self.path[p1].0 == 1 {
            self.path[p1].1
        } else {
            self.path[p1].1 + N
        };
        let p2_index = if self.path[p2].0 == 1 {
            self.path[p2].1
        } else {
            self.path[p2].1 + N
        };
        let (next_path_p1_p2, next_energy_p1_p2, sum_energy_p1_p2) = update_path(p1_index, p2_index, &self.pos, &self.dist_warshall_floyd, &mut visited);
        // 新しいパス･エネルギーを結合
        let next_path = [path_s_p1, &next_path_p1_p2, path_p2_t].concat();
        let next_energy = [energy_s_p1, &next_energy_p1_p2, energy_p2_t].concat();
        tmp_energy += sum_energy_p1_p2;
        self.path = next_path;
        self.energy = next_energy;
        self.sum_energy = tmp_energy;
    }
    // スコア計算
    pub fn compute_score(&mut self) {
        self.game_score = (1_000_000_000_f64 / (1_000_f64 + (self.sum_energy as f64).sqrt())).round() as i64;
    }
    // 1-indexに直しつつ回答を出力
    pub fn print(&self) {
        for i in 0..M {
            println!("{} {}", self.pos[i + N][0], self.pos[i + N][1]);
        }
        let v = self.path.len();
        println!("{}", v);
        for i in 0..v {
            let (ti, ri) = self.path[i];
            println!("{} {}", ti, ri + 1);
        }
    }
}

fn main() {
    // 2-optに見えるが，再訪ありなので単純には適応できない
    // ただ，宇宙ステーションに再訪してもいいが，惑星は再訪しないほうがよい？
    // 近傍: 
    // - 2-opt (ただし，最初と最後が惑星1なのは固定する)
    // - 少なくとも一方が宇宙ステーションであるパスの追加
    // - いずれかの宇宙ステーションの座標をずらす (あるいは，これはあらかじめk-meansで求めておく)
    get_time();
    let input = Input::new();
    let mut state = State::new(&input);
    state.compute_score();

    let mut iter = 0;
    let mut best_score = state.game_score;
    let mut best_state = state.clone();
    while get_time() < TL {
        iter += 1;
        let t = get_time() / TL;
        let T = T0.powf(1.0 - t) * T1.powf(t);
        let mut next_state = state.clone();
        if rand_gen_bool(0.2) {
            // eprintln!("shift");
            next_state.shift_station();
        } else {
            // eprintln!("scrap_and_build");
            next_state.scrap_and_build();
        }
        // next_state.scrap_and_build();
        next_state.compute_score();

        let crt_score = state.game_score;
        let next_score = next_state.game_score;
        // ベスト更新は常に採用する
        if next_score > best_score {
            best_score = next_score;
            best_state = next_state.clone();
            // eprintln!("update best: {}", best_score);
            // eprintln!("update iter: {}", iter);
        }
        // 現時刻と次時刻の更新については，スコアが高いか，焼きなましの許容範囲なら更新
        if next_score >= crt_score || rand_gen_bool(((next_score - crt_score) as f64 / T).exp()) {
            state = next_state;
            // if next_score - crt_score < 0 {eprintln!("annealing");}
        }
    }

    best_state.print();
    eprintln!("iter: {}", iter);
    eprintln!("score: {}", best_state.game_score);
    eprintln!("time: {}", get_time());
}