結果

問題 No.5007 Steiner Space Travel
ユーザー r3yoheir3yohei
提出日時 2023-05-02 02:16:44
言語 Rust
(1.77.0 + proconio)
結果
AC  
実行時間 982 ms / 1,000 ms
コード長 21,393 bytes
コンパイル時間 1,730 ms
コンパイル使用メモリ 170,868 KB
実行使用メモリ 4,500 KB
スコア 8,830,879
最終ジャッジ日時 2023-05-02 02:17:20
合計ジャッジ時間 34,027 ms
ジャッジサーバーID
（参考情報） 		judge15 / judge12
純コード判定しない問題か言語
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 981 ms
4,368 KB
testcase_01 AC 981 ms
4,368 KB
testcase_02 AC 981 ms
4,372 KB
testcase_03 AC 982 ms
4,368 KB
testcase_04 AC 982 ms
4,368 KB
testcase_05 AC 981 ms
4,372 KB
testcase_06 AC 982 ms
4,372 KB
testcase_07 AC 982 ms
4,368 KB
testcase_08 AC 981 ms
4,368 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 982 ms
4,372 KB
testcase_13 AC 981 ms
4,372 KB
testcase_14 AC 982 ms
4,368 KB
testcase_15 AC 982 ms
4,372 KB
testcase_16 AC 981 ms
4,372 KB
testcase_17 AC 982 ms
4,372 KB
testcase_18 AC 982 ms
4,372 KB
testcase_19 AC 981 ms
4,368 KB
testcase_20 AC 981 ms
4,368 KB
testcase_21 AC 981 ms
4,368 KB
testcase_22 AC 982 ms
4,500 KB
testcase_23 AC 981 ms
4,368 KB
testcase_24 AC 982 ms
4,372 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 981 ms
4,372 KB
testcase_29 AC 982 ms
4,372 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#![allow(non_snake_case, unused)]

use std::{io::*, hash};
use std::{collections::*, fmt::format};
use std::{cmp::*, vec};
use crate::input::{*, marker::*};
use crate::rand::Xoshiro256;
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());
    }
}


mod rand {
    pub(crate) struct Xoshiro256 {
        s0: u64,
        s1: u64,
        s2: u64,
        s3: u64,
    }

    impl Xoshiro256 {
        pub(crate) fn new(mut seed: u64) -> Self {
            let s0 = split_mix_64(&mut seed);
            let s1 = split_mix_64(&mut seed);
            let s2 = split_mix_64(&mut seed);
            let s3 = split_mix_64(&mut seed);
            Self { s0, s1, s2, s3 }
        }

        fn next(&mut self) -> u64 {
            let result = (self.s1.wrapping_mul(5)).rotate_left(7).wrapping_mul(9);
            let t = self.s1 << 17;

            self.s2 ^= self.s0;
            self.s3 ^= self.s1;
            self.s1 ^= self.s2;
            self.s0 ^= self.s3;
            self.s2 ^= t;
            self.s3 = self.s3.rotate_left(45);

            result
        }

        pub(crate) fn gen_usize(&mut self, lower: usize, upper: usize) -> usize {
            assert!(lower < upper);
            let count = upper - lower;
            (self.next() % count as u64) as usize + lower
        }

        pub(crate) fn gen_i64(&mut self, lower: i64, upper: i64) -> i64 {
            assert!(lower < upper);
            let count = upper - lower;
            (self.next() % count as u64) as i64 + lower
        }

        pub(crate) fn gen_f64(&mut self) -> f64 {
            const UPPER_MASK: u64 = 0x3ff0000000000000;
            const LOWER_MASK: u64 = 0xfffffffffffff;
            let result = UPPER_MASK | (self.next() & LOWER_MASK);
            let result: f64 = unsafe { std::mem::transmute(result) };
            result - 1.0
        }

        pub(crate) fn gen_bool(&mut self, prob: f64) -> bool {
            self.gen_f64() < prob
        }
    }

    fn split_mix_64(x: &mut u64) -> u64 {
        *x = x.wrapping_add(0x9e3779b97f4a7c15);
        let mut z = *x;
        z = (z ^ z >> 30).wrapping_mul(0xbf58476d1ce4e5b9);
        z = (z ^ z >> 27).wrapping_mul(0x94d049bb133111eb);
        z ^ z >> 31
    }
}

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 MAP_SIZE: i64 = 1000;
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(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 }
    }

    fn dist_sq(&self, other: &Self) -> i64 {
        let dx = self.x - other.x;
        let dy = self.y - other.y;
        (dx * dx + dy * dy) as i64
    }
}

#[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![];
        for _ in 0..N {
            input! {
                ai: i64,
                bi: i64,
            }
            pos.push(Position::new(ai, bi));
        }
        let dist = warshall_floyd(&pos);
        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 != !0 {
        // toがstartになるまでprev[to]を辿っていく
        // 途中のエネルギーを保存しておく
        path.push(tt);
        tt = prev[tt];
    }
    path.reverse();
    path
}

#[derive(Clone)]
struct State {
    pos: Vec<Position>, // 先頭N個に惑星の場所，続くM個にステーションの場所が格納されている
    path: Vec<usize>, // 道順 (惑星orステーション，番号)
}
impl State {
    fn new(input: &Input) -> Self {
        // 惑星の位置をinputからもらってくる
        let mut pos = input.pos.clone();
        // ステーションの位置は雑に初期化する
        for _ in 0..M {
            pos.push(pos[0]);
        }
        // 道順は，とりあえず0~N,0の順に通るよう初期化する
        let mut path = vec![];
        for i in 0..N {
            path.push(i);
        }
        path.push(0);
        Self {pos, path}
    }
    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 {
            let dist_sq = squared_distance(self.pos[prev], prev, self.pos[next], next);
            dist_sq
        }
    }
    // 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(&input.pos, prev, next);
        for v in path {
            new_solution.path.push(v);
        }
    }

    new_solution
}


fn main() {
    get_time();
    let mut rng = Xoshiro256::new(42);
    let input = Input::new();
    let mut state = State::new(&input);

    let mut iter = 0;
    let mut accepted_count = 0;
    let mut update_best_count = 0;
    let mut current_score = state.calc_score_all(&input);
    let mut best_score = current_score;
    let mut best_state = state.clone();
    while get_time() < TL {
        iter += 1;
        // if (iter & ((1 << 6) - 1)) == 0 {
            
        // }
        let t = get_time() / TL;
        let T = T0.powf(1.0 - t) * T1.powf(t);
        // 乱数で生じさせた数値に応じて近傍を選択
        let neighbor_type = rng.gen_usize(0, 10);
        if neighbor_type < 5 {
            // 近傍3: ある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;
            const MAX_DELTA: f64 = 400.;
            const MIN_DELTA: f64 = 10.;
            let delta = (MAX_DELTA * (1. - get_time()) + MIN_DELTA * get_time()) 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::with_capacity(state.path.len());
            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 || 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;
            }
        } else {
            // 近傍4: 2-opt
            let from = rng.gen_usize(1, state.path.len() - 1);
            let to = rng.gen_usize(from + 1, state.path.len());

            let i0 = state.path[from - 1];
            let i1 = state.path[from];
            let i2 = state.path[to - 1];
            let i3 = state.path[to];

            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 || rng.gen_bool(f64::exp(-score_diff as f64 / T)) {
                // 解の更新
                current_score = new_score;
                accepted_count += 1;
                state.path[from..to].reverse();
            }
        }
        // どの近傍でもベスト更新は常に採用する
        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);
        }
    }

    best_state.print();
    eprintln!("iter: {}", iter);
    eprintln!("accepted: {}", accepted_count);
    eprintln!("update best: {}", update_best_count);
    eprintln!("energy: {}", best_score);
    eprintln!(
        "score: {}",
        (1e9 / (1e3 + (best_score as f64).sqrt())).round() as i64
    );
    eprintln!("time: {}", get_time());
}
0