#864786 (C++17(gcc12)) No.5007 Steiner Space Travel

• 問題: No.5007 Steiner Space Travel
• ユーザー: bird01
• 提出日時: 2023-04-29 18:41:43
• 言語: C++17(gcc12) (gcc 12.3.0 + boost 1.87.0)
• 結果: TLE
• 実行時間: -
• コード長: 17,395 bytes
• コンパイル時間: 2,418 ms
• コンパイル使用メモリ: 144,940 KB
• 実行使用メモリ: 4,500 KB
• スコア: 0
• 最終ジャッジ日時: 2023-04-29 18:42:29
• 合計ジャッジ時間: 45,418 ms
• ジャッジサーバーID （参考情報）: judge11 / judge16

ファイルパターン	結果
other	TLE * 30

ソースコード

// 1.惑星だけでNearstNeighbor法
// 2.k-means++法でクラスタリングして宇宙ステーションを配置
// 2.calc_dist_planetの距離計算が間違ていたので修正
// 3.各辺の間に宇宙ステーションの挿入を試みる(cuの所属クラスタの宇宙ステーションの経由のみを考慮する簡易版)
// 4.2-opt
// 5.焼きなまし,calc_score,2-opt
// 6.Or-opt
// 7.挿入

#include <algorithm>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <iomanip>
#include <assert.h>

#define rep(i,a,b) for(int i=a;i<b;i++)
#define fore(i,a) for(auto &i:a)
#define all(x) (x).begin(),(x).end()
#define sz(x) (int)((x).size())
#define fo(x) cout<<x<<endl
#define vi vector<int>
#define vvi vector<vi>
#define vvvi vector<vvi>
#define vvvvi vector<vvvi>
#define vp vector<P>
template <typename T> bool chmin(T &a, const T& b){if(a>b){a=b; return true;}return false;}
template <typename T> bool chmax(T &a, const T& b){if(a<b){a=b; return true;}return false;}
#define dump(x)  cerr << #x << " = " << (x) << endl;
#define debug(x) cerr << #x << " = " << (x) << " (L" << __LINE__ << ")" << " " << __FILE__ << endl;
using namespace std;
using ll = long long;
using P = pair<int,int>;
template<typename T> using min_priority_queue = priority_queue<T, vector<T>, greater<T>>;
//const long long INF=1e18;
const int INF=1000000000;

const int64_t CYCLES_PER_SEC = 3000000000;
const double TIMELIMIT = 0.95;
constexpr int RANDOM_SEED = 42; // 乱数のシード値

struct Timer {
  int64_t start;
  Timer() { reset(); }
  void reset() { start = getCycle(); }
  void plus(double a) { start -= (a * CYCLES_PER_SEC); }
  inline double get() { return (double)(getCycle() - start) / CYCLES_PER_SEC; }
  inline int64_t getCycle() {
    uint32_t low, high;
    __asm__ volatile ("rdtsc" : "=a" (low), "=d" (high));
    return ((int64_t)low) | ((int64_t)high << 32);
  }
};

mt19937 mt(RANDOM_SEED);
uniform_int_distribution<> rand24(0,23);

struct Input{
  int N,M;
  vector<P> coord;
  Input(int N=0,int M=0, vector<P> coord={}): N(N), M(M), coord(coord) {}
};

struct Solver {
  Timer timer;

  // グローバル変数==============================
  Input input;
  ll best_score;
  vector<P> best_ans; // 経由順
  vector<P> best_station; // 宇宙ステーションの座標
  vvi dist_planet; // 全惑星間の距離

  // 便利関数==============================
  // 乱数
  double Randouble(){
    return 1.0 * rand() / RAND_MAX;
  }
  // 距離の2乗を求める
  template<typename T>
  T dist2(T a,T b){return a*a+b*b;}
  // ans上の2点間の距離を求める
  int ans_to_dist2(Input &input, vector<P> &ans, vector<P> &station, int i, int j){
    int t1,idx1,t2,idx2;
    tie(t1,idx1)=ans[i];
    tie(t2,idx2)=ans[j];
    // 両方惑星なら早期return
    if(t1==1 and t2==1) return dist_planet[idx1-1][idx2-1];

    int x1,y1,x2,y2;
    if(t1==1) tie(x1,y1)=input.coord[idx1-1];
    else tie(x1,y1)=station[idx1-1];
    if(t2==1) tie(x2,y2)=input.coord[idx2-1];
    else tie(x2,y2)=station[idx2-1];
    int d=dist2(x2-x1,y2-y1);
    if(t1==1) d*=5; // 惑星は5倍
    if(t2==1) d*=5; // 惑星は5倍
    return d;
  }
  // 区間を逆回りにする
  void reverse_section(vector<P> &ans, int l, int r){
    while(l<r) swap(ans[l++],ans[r--]);
  }
  // 2-opt
  void two_opt(Input &input, vector<P> &ans, vector<P> &station, int i, int j){
    // 距離を求める
    int d1 = ans_to_dist2(input,ans,station,i,i+1) + ans_to_dist2(input,ans,station,j,j+1);
    int d2 = ans_to_dist2(input,ans,station,i,j) + ans_to_dist2(input,ans,station,i+1,j+1);
    // 改善できるなら入れ替える
    if(d2<d1) reverse_section(ans,i+1,j);
  }
  // すべてを2-opt
  void two_opt_all(Input &input, vector<P> &ans, vector<P> &station){
    rep(i,0,sz(ans)-3){
      rep(j,i+2,sz(ans)-1){ // 始点と終点は同一なので-1
        two_opt(input,ans,station,i,j);
      }
    }
  }
  // k-means++法でM個のクラスタ中心を求める
  vi select_cluster_centers(Input &input){
    vi centers; // クラスタ中心の集合
    vi dist(input.N,INF); // 最寄のクラスタ中心までの距離

    // データ点からランダムに1つ選びクラスタ中心とする
    uniform_int_distribution<> randN(0,input.N-1);
    int center=randN(mt);
    centers.emplace_back(center); // クラスタ中心を追加

    // 残りのM-1個のクラスタ中心を求める
    rep(i,0,input.M-1){
      int cx,cy;
      tie(cx,cy)=input.coord[center];

      // 各点について最寄のクラスタ中心までの距離を求める
      rep(j,0,input.N){
        int x,y;
        tie(x,y)=input.coord[j];
        chmin(dist[j],dist2(x-cx,y-cy));
      }

      // 最寄のクラスタ中心までの距離の内,最長のものを持つ点を求める
      int d=0;
      rep(j,0,input.N){
        if(chmax(d,dist[j])){
          center=j;
        }
      }

      centers.emplace_back(center); // クラスタ中心を追加
    }

    return centers;
  }
  // k-means法でクラスタリング
  vi k_means(Input &input, vector<pair<double,double>> &center_pos){
    vi cluster_num(input.N,-1); // 各点が属するクラスタ番号
    while(1){
      vi cluster_num_buf(input.N,-1);
      vector<double> mn_dist(input.N,INF); // 各点の最寄クラスタまでの距離
      // 各点を最寄のクラスタ中心に分類する
      rep(i,0,input.N){ // 点をまわす
        int x,y;
        tie(x,y)=input.coord[i];
        rep(j,0,input.M){ // クラスタ中心をまわす
          double cx,cy;
          tie(cx,cy)=center_pos[j];
          // 距離を更新
          if(chmin(mn_dist[i],dist2(x-cx,y-cy))){
            cluster_num_buf[i]=j;
          }
        }
      }

      // 終了判定,更新
      if(cluster_num==cluster_num_buf) break; // 変化がなくなったら抜ける
      swap(cluster_num,cluster_num_buf);

      // 新しいクラスタ中心を求める
      rep(i,0,input.M) center_pos[i]=make_pair(0,0); // 初期化
      vi point_cnt(input.M,0); // そのクラスタに所属する点の個数
      // 座標を加算
      rep(i,0,input.N){
        int c_idx=cluster_num[i];
        int x,y;
        tie(x,y)=input.coord[i];
        center_pos[c_idx].first+=x;
        center_pos[c_idx].second+=y;
        point_cnt[c_idx]++; // クラスタに所属する点の個数をインクリメント
      }
      // 点の数で除算
      rep(i,0,input.M){
        center_pos[i].first=round((double)center_pos[i].first/point_cnt[i]);
        center_pos[i].second=round((double)center_pos[i].second/point_cnt[i]);
      }
    }
    return cluster_num;
  }
  // 宇宙ステーションを配置,各惑星のクラスター番号も返す
  vector<P> locate_station(Input &input,vi &cluster_num){
    // k-means++
    vi centers=select_cluster_centers(input); // M個のクラスタ中心を求める
    vector<pair<double,double>> center_pos(input.M); // 各クラスタ中心の座標 // 0初期化されるか
    rep(i,0,input.M){
      center_pos[i]=input.coord[centers[i]];
    }
    // k-means
    cluster_num=k_means(input,center_pos); // cluster_numは各惑星のクラスター番号,center_posも変わる
    // center_posをdoubleからintにする
    vector<P> station_pos(input.M);
    rep(i,0,input.M){
      station_pos[i].first=round(center_pos[i].first);
      station_pos[i].second=round(center_pos[i].second);
    }

    return station_pos;
  }
  // 各辺の間に宇宙ステーションの挿入を試みる
  vector<P> insert_station(Input &input, vector<P> &ans, vector<P> station, vi cluster_num){
    vector<P> res;
    res.emplace_back(ans[0]); // スタートはそのまま追加
    // 各径路間にステーションを挿入するか判定
    rep(i,0,sz(ans)-1){
      int planet_idx_cu=ans[i].second-1;
      int planet_idx_nx=ans[i+1].second-1;
      // 距離が有利なら挿入する
      int station_num=cluster_num[planet_idx_cu]; // 宇宙ステーションの番号
      int px1,py1,px2,py2,sx,sy;
      tie(px1,py1)=input.coord[planet_idx_cu]; // 惑星の座標
      tie(px2,py2)=input.coord[planet_idx_nx]; // 惑星の座標
      tie(sx,sy)=station[station_num];// 宇宙ステーションの座標
      int d=dist2(px1-sx,py1-sy)*5 + dist2(px2-sx,py2-sy)*5; // 宇宙ステーションを経由する場合のコスト
      if(dist_planet[planet_idx_cu][planet_idx_nx]>d){
        res.emplace_back(P(2,station_num+1));
      }
      // 次の惑星に行く
      res.emplace_back(ans[i+1]);
    }
    return res;
  }
  // スコア計算
  ll calc_score(Input &input, vector<P> &ans, vector<P> &station){
    int sum=0;
    rep(i,0,sz(ans)-1){
      sum+=ans_to_dist2(input,ans,station,i,i+1);
    }
    return round(1000000000/(1000+sqrt(sum)));
  }
  // 全惑星間の距離(a^2*D^2)を求める
  vvi calc_dist_planet(Input &input){
    vvi res(input.N,vi(input.N,0));
    rep(i,0,input.N-1){
      int x1,y1;
      tie(x1,y1)=input.coord[i];
      rep(j,i+1,input.N){
        int x2,y2;
        tie(x2,y2)=input.coord[j];
        int d=dist2(x2-x1,y2-y1)*25;
        res[i][j]=d;
        res[j][i]=d;
      }
    }
    return res;
  }
  // 初期解を生成する
  vector<P> get_initial_solution(Input &input, vector<P> &station, vvi &dist_planet){
    vector<P> res;
    // 宇宙ステーションの位置を決める,各惑星のクラスタ番号を得る=====
    vi cluster_num; // 各惑星のクラスタ番号
    station=locate_station(input,cluster_num);
    // 惑星だけでNearstNeighbor法=====
    vi visited(input.N,0);
    visited[0]=1;
    int cu=0;
    res.emplace_back(P(1,1)); // 最初の惑星からスタート
    rep(i,0,input.N-1){
      // 最短距離の惑星を求める
      int nx=-1;
      int mn_d=INF;
      rep(i,0,input.N){
        if(visited[i]) continue;
        if(chmin(mn_d,dist_planet[cu][i])){
          nx=i;
        }
      }
      res.emplace_back(P(1,nx+1));
      visited[nx]=1;
      swap(cu,nx);
    }
    res.emplace_back(P(1,1)); // 最初の惑星に戻る
    // 各辺の間に宇宙ステーションを挿入
    res=insert_station(input,res,station,cluster_num);

    // 2-opt
    /*
    uniform_int_distribution<> rand_ans(0,sz(res)-2);
    rep(i,0,100000){
      int a=rand_ans(mt);
      int b=rand_ans(mt);
      while(abs(a-b)<=1) b=rand_ans(mt);
      two_opt(input,res,station,a,b);
    }
    */
    /*
    rep(i,0,10){
      two_opt_all(input,res,station);
    }
    */

    // 各辺の間に宇宙ステーションを挿入
    //res=insert_station(input,res,station,cluster_num);

    return res;
  }
  // 2-opt(make_new_ans用)
  void mna_two_opt(vector<P> &ans, vector<P> &station){
    uniform_int_distribution<> rand_ans(0,sz(ans)-2);
    int a=rand_ans(mt);
    int b=rand_ans(mt);
    while(abs(a-b)<=1) b=rand_ans(mt);
    two_opt(input,ans,station,a,b);
  }
  // Or-opt(make_new_ans用)
  void mna_or_opt(Input &input, vector<P> &ans, vector<P> &station){
    uniform_int_distribution<> rand_ans(0,sz(ans)-2);
    // 移動対象の区間を選択(区間長は1に固定)
    // 0は選ばない,sz(ans)-2は選ばない
    int tar_section=rand_ans(mt);
    while(tar_section==0 or tar_section==sz(ans)-2) tar_section=rand_ans(mt);

    // 移動先を選択
    // tar_section-1,tar_section,tar_section+1は選ばない
    int saki=rand_ans(mt);
    while(tar_section-1<=saki and saki<=tar_section+1) tar_section=rand_ans(mt);

    int a=tar_section-1;
    int b=tar_section;
    int c=tar_section+1;
    int d=tar_section+2;
    int e=saki;
    int f=saki+1;
    int cu_d=ans_to_dist2(input,ans,station,a,b)+ans_to_dist2(input,ans,station,c,d)+ans_to_dist2(input,ans,station,e,f);
    int nx_d1=ans_to_dist2(input,ans,station,a,d)+ans_to_dist2(input,ans,station,e,b)+ans_to_dist2(input,ans,station,c,f);
    int nx_d2=ans_to_dist2(input,ans,station,a,d)+ans_to_dist2(input,ans,station,e,c)+ans_to_dist2(input,ans,station,b,f); // 反転あり

    bool rev=nx_d1>nx_d2; // 区間反転した方が良い場合
    int mn_nx_d=min(nx_d1,nx_d2);
    // 更新可能な場合
    if(mn_nx_d < cu_d){
      P tar1,tar2;
      // 反転を伴う場合
      if(rev){
        tar1=ans[c];
        tar2=ans[b];
      // 反転を伴わない場合
      }else{
        tar1=ans[b];
        tar2=ans[c];
      }

      // 更新
      vector<P> res;
      rep(i,0,sz(ans)){
        if(i==b or i==c) continue; // 移動対象はスキップ
        res.emplace_back(ans[i]);
        if(i==e){
          res.emplace_back(tar1);
          res.emplace_back(tar2);
        }
      }
      swap(ans,res);
    }
  }
  // 挿入
  void mna_insert(Input &input, vector<P> &ans, vector<P> &station){
    // 挿入対象の点を選ぶ
    // 始点,終点以外の点をランダムに指定
    uniform_int_distribution<> rand_ans(0,sz(ans)-2);
    int tar=rand_ans(mt);
    while(tar==0) tar=rand_ans(mt);
    
    // 被挿入対象の点を選ぶ
    int insert_tar=rand_ans(mt);
    while(insert_tar==tar-1 or insert_tar==tar) insert_tar=rand_ans(mt);

    // 1.その点をなくして別の場所に入れる
    // 2.その点はそのままに別の場所に入れる
    int a=tar-1;
    int b=tar;
    int c=tar+1;
    int d=insert_tar;
    int e=insert_tar+1;
    int cu_d=ans_to_dist2(input,ans,station,a,b)+ans_to_dist2(input,ans,station,b,c)+ans_to_dist2(input,ans,station,d,e);
    int nx_d1=ans_to_dist2(input,ans,station,a,c)+ans_to_dist2(input,ans,station,d,b)+ans_to_dist2(input,ans,station,b,e); // 削除あり
    int nx_d2=ans_to_dist2(input,ans,station,a,b)+ans_to_dist2(input,ans,station,b,c)+ans_to_dist2(input,ans,station,d,b)+ans_to_dist2(input,ans,station,b,e); // 削除なし
    
    bool is_remove_b=false;
    // 1.何もしない
    if(cu_d<=nx_d1 and cu_d<=nx_d2){
      return;
    // 2.点をなくして入れる
    }else if(nx_d1<cu_d and nx_d1<nx_d2){
      is_remove_b=true;
    // 3.点はそのままで入れる
    }else if(nx_d2<cu_d and nx_d2<nx_d1){
      // 何もしない
    }
    
    // 更新
    vector<P> res;
    rep(i,0,sz(ans)){
      if(is_remove_b and i==b) continue; // 挿入対象はスキップ
      res.emplace_back(ans[i]);
      if(i==d){
        res.emplace_back(ans[b]);
      }
    }
    swap(ans,res);
  }
  // 交換
  // Or-opt
  // λ-opt
  // pre_ansからnew_ansを作る
  void make_new_ans(Input &input, vector<P> &ans, vector<P> &station){
    uniform_int_distribution<> rand_mna(0,2);
    int sel=rand_mna(mt);
    if(sel==0) mna_two_opt(ans,station);
    else if(sel==1) mna_or_opt(input,ans,station);
    else if(sel==2) mna_insert(input,ans,station);
  }
  
  
  // 初期解生成
  void init(){
    dist_planet=calc_dist_planet(input); // 全惑星間の距離を求める
    vector<P> station;
    vector<P> ans=get_initial_solution(input,station,dist_planet);
    best_ans=ans;
    best_station=station;
  }
  
  void SA(){
    double ti = timer.get();
    double start_time=ti;

    best_score = calc_score(input,best_ans,best_station);
    vector<P> pre_ans=best_ans;
    vector<P> pre_station=best_station;
    ll pre_score = best_score;
    cerr << best_score << endl;
    double kankaku = TIMELIMIT / 10;
    double ne = ti;
    int cnt = 0;
    int cnt2 = 0;

    double start_temp = 10;
    double end_temp = 1;
    while(1){
      ti = timer.get();
      if(TIMELIMIT < ti-start_time) break;
      
      double temp = start_temp + (end_temp - start_temp) * (ti-start_time) / TIMELIMIT;

      if (ti > ne) {
        ne += kankaku;
        cerr << temp << " " << pre_score << " " << best_score << endl;
      }
      
      // 新たな解の生成==============================
      vector<P> new_ans=pre_ans;
      vector<P> new_station=pre_station;
      make_new_ans(input,new_ans,new_station);
      
      // new_scoreを計算
      ll new_score=calc_score(input,new_ans,new_station); // スコア計算
      
      // 採用するか否かを決定
      double prob = exp((new_score-pre_score)/temp);
      //prob/=1000000000;
      if (prob*10000 > 10000 * (rand() % 10000)) {
        pre_ans=new_ans;
        pre_station=new_station;
        pre_score=new_score;

        if(chmax(best_score,new_score)){
          best_ans=new_ans;
          best_station=new_station;
          best_score=new_score;
        }
        cnt2++;
      }
      cnt++;
    }
    
    dump(cnt);
		dump(cnt2);
    //ti = timer.get();
    //cerr << ti << endl;
  }
  
  void input_func() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    int N,M;
    cin >> N >> M;
    vector<P> coord(N);
    rep(i,0,N){
      int a,b;
      cin >> a >> b;
      coord[i]=P(a,b);
    }
    input={N,M,coord};
  }

  int solve() {
    init();
    SA();
    return 0;
  }
  
  void output(){
    dump(calc_score(input,best_ans,best_station));
    rep(i,0,input.M){
      int x,y;
      tie(x,y)=best_station[i];
      cout << x << ' ' << y << endl;
    }
    fo(sz(best_ans));
    rep(i,0,sz(best_ans)){
      int t,r;
      tie(t,r)=best_ans[i];
      cout << t << ' ' << r << endl;
    }
  }
};

Solver solver;
int main(int argc, char *argv[]) {
  solver.input_func();
  solver.solve(); // 入力に合うグラフを作成
  solver.output(); // クエリに答える
  return 0;
}