Submission #44786610

---

Source Code Expand

#include <bits/stdc++.h>
#ifdef MY_DEBUG
  #include "dprint.cpp"
  # define dprint(a, ...);\
      do {\
          printf("line : %d, func : %s\n",\
                __LINE__, __func__);\
          dprint(#a,a,##__VA_ARGS__ );\
      } while (0)
#else
  #define dprint(a,...) 
#endif
using namespace std;
using P =pair<int,int>;

int L,N,S;
vector<P> EXIT_COR;
int measure_count=0;
int n_sample=1;
struct parameta{
  int pw1=20;int pw2=90;
  int h1=50;int h2=(int)min(1000.0,23.635*sqrt(S)+78);
  int baseline=min(0,1000-h2);
  int h00=(int)(49*sqrt(S))+50;
  double acc_prob=0.8;
  double acc_prob2=0.9;
  parameta(){
    h00=(int)(49*sqrt(S)+50);
    if(h00>=1000){
      h00=1000;
      h2=(-6.9*sqrt(S))+634;
      h1=10.7*sqrt(S)-177.19;
    }
    if(S==900||S==841||S==784){
      h1=950;
      h2=1000;
    }
  }
};
parameta param;

// 正規分布の累積密度関数を計算する関数
vector<double>normalCDF_memo(2002,0);
void normalCDF_init(int mean,int stddev){
  for (int t = -1000; t <= 1000; t++){
    normalCDF_memo[t+1000]=0.5 * (1.0 + std::erf(((double)t - mean) / (stddev * std::sqrt(2.0))));
  }
}
double normalCDF(int x,int mean) {
  return normalCDF_memo[x-mean+1000];
    // return 0.5 * (1.0 + std::erf((x - mean) / (stddev * std::sqrt(2.0))));
}


struct temperature{
  vector<vector<int>>T;

  void prickle(int y,int x,int h,int w){
    assert (y-w>=0);
    assert (y+w<L);
    assert (x-w>=0);
    assert (x+w<L);
    if(w==0){
      T[y][x]+=h;
      return;
    }

    for (int i = -w; i <= w; i++){
      for (int j = -w; j <= w; j++){
        int dy=abs(i);
        int dx=abs(j);
        T[y+i][x+j]+=(w-max(dx,dy)+1)*h/w;
      }
    }
  }

  void init(){
    int pw1=param.pw1;int pw2=param.pw2;
    int h1=param.h1;int h2=param.h2;
    int baseline=param.baseline;

    //temp
    T=vector<vector<int>>(L,vector<int>(L,baseline));
    int c1=L/3/2;
    int c2=L/2;//もうっか所あるので注意
    int c3=L-1-L/3/2;

    int w1=c1*pw1/100; 
    int w2=c2*pw2/100-1;
    //prickle(c1,c1,h1,w1);
    prickle(c1,c2,h1,w1);
    //prickle(c1,c3,h1,w1);

    prickle(c2,c1,h1,w1);
    prickle(c2,c2,h2,w2);
    prickle(c2,c3,h1,w1);
    
    //prickle(c3,c1,h1,w1);
    prickle(c3,c2,h1,w1);
    //prickle(c3,c3,h1,w1);
    // T[0][0]=1000;
    T[c2][c2]=param.h00;

    //clip
    for (int i = 0; i < L; i++){
      for (int j = 0; j < L; j++){
       T[i][j]=max(0,min(1000,T[i][j])); 
      }
    }
  }
  
  int value(int y,int x){
    return T[y][x];
  }

  double get_p(int y,int x,int low,int high){//座標y,xで、温度lowからhighを出力する確率。
    // clip影響
    double res=normalCDF(high, T[y][x]) - normalCDF(low ,T[y][x]);
    if(low==0) res=normalCDF(high, T[y][x]) - 0;
    if(high==1000)res=1 - normalCDF(low ,T[y][x]);
    return res;
  }

  double get_p(int y,int x,int measure){//座標y,xで、温度measureを出力する確率。
    //T[y][x]
    // clip影響
    double res= normalCDF(measure+1, T[y][x]) - normalCDF(measure ,T[y][x]);
    if(measure==0) res=normalCDF(measure+1, T[y][x]) -0;
    if(measure==1000)res=1-normalCDF(measure ,T[y][x]);
    return res;
  }

};

temperature Temp;

struct Judge {
    void print_temperature(temperature& T) {
        for (int i = 0; i < L; i++){
          for (int j = 0; j < L; j++){
            cout << T.value(i,j) << " \n"[j==L-1];
          } 
        }
        cout.flush();
    }
 
    int measure(int i, int y, int x) {
      int ave=0;
      for (int t = 0; t < n_sample; t++){
          measure_count++;
          cout << i << " " << y << " " << x << endl; // endl does flush
          int v;
          cin >> v;
          ave+=v;
          if (v == -1) {
              cerr << "#something went wrong. i=" << i << " y=" << y << " x=" << x << endl;
              exit(1);
          }
      }

        return ave/n_sample;
    }
 
    void answer(const vector<pair<int,double>>& estimate) {
        cout << "-1 -1 -1" << endl;
        for (int i = 0; i < N; i++){
            cout<<"#entrance i:"<<i<<"probablity:"<<estimate[i].second<<endl;
            cout << estimate[i].first << endl;
            
        } 

    }
};

vector<vector<double>>KLmemo;
double KL(int  m1,int m2){
  if(KLmemo[m1][m2]!=-1)return KLmemo[m1][m2];
  //ex1とex2の分布の距離を測る
  double ans=0;
  //https://sucrose.hatenablog.com/entry/2013/07/20/190146

  ans=(double)(m1-m2)/S;
  ans=ans*ans/2;
  KLmemo[m1][m2]=ans;
  return ans;

}

void update_ll(vector<vector<double>>& ll,int m,int index,int y, int x){
  for (int ex = 0; ex < N; ex++){
    double p=Temp.get_p((EXIT_COR[ex].first+y+L)%L,(EXIT_COR[ex].second+x+L)%L,m);
    ll[index][ex]+=log10(p);
  }
  //最大の数で割り算。
  double maxi=-DBL_MAX;
  for (int i = 0; i < N; i++){
   maxi=max(maxi,ll[index][i]); 
  }
  for (int i = 0; i < N; i++){
   ll[index][i]-=maxi; 
  }

  ll[index][N]=0;
  for (int i = 0; i < N; i++){//N番目には合計値を保存しておく。(高速化用) ここは最大値
   ll[index][N]+=pow(10,ll[index][i]); 
  }

}

double ll_to_p(const vector<double> &ll,int ex){
  return pow(10,ll[ex])/ll[N];//後半はupdateの時に取得済み
}

P global_cor(int ex,int y,int x){
  return make_pair((EXIT_COR[ex].first+y+L)%L,(EXIT_COR[ex].second+x+L)%L);
}

double calc_IAB(const vector<pair<double,int>> &Pa,int y,int x){
  double ans=0;//よくわからない値 分布の距離の期待値
  for (int ti = 0; ti <4; ti++){//N<=10 本とっはmin
   for (int tj = ti+1; tj <4; tj++){
    double pi=(Pa[ti].first);
    double  pj=(Pa[tj].first);

    P x1=global_cor(Pa[ti].second,y,x);
    P x2=global_cor(Pa[tj].second,y,x);
    // cout << "#"<<ti<<","<<tj << endl;
    ans+=pi*pj*KL(Temp.value(x1.first,x1.second),Temp.value(x2.first,x2.second));
   }   
  }
  // double scale=0;
  // for (int i = 0; i < 10; i++){
  //   scale+=(double)(Pa[i]/100)/1e6;
  // }
  return ans;
}
void select_yx2(double& maxi,int & y,int& x,const  vector<pair<double,int>>&Pa,int cy,int cx){
  //本当はいらないんだが、dx+=2,dy+=2のために必要
    P pointmax= EXIT_COR[Pa[0].second];
    int dy=-(pointmax.first-cy);
    int dx=-(pointmax.second-cx);
    double ret= calc_IAB(Pa,dy,dx);//llの情報を持っているとき、dydxを叩くことで得られる情報量。（減少するエントロピーの量）
    if(ret>maxi){
      // cout << "#kokoko" << endl;
      maxi=ret;
      y=dy;
      x=dx;
    }
}
P select_yx(const vector<double> &ll){
  int y=0;
  int x=0;
  double maxi=0;
  //calc_IAB内で使うが、使いまわせるので先に計算しておく。
  vector<pair<double,int>>Pa(N);
  for (int i = 0; i < N; i++){
    // int t=ll_to_p(ll,i)*1e7;
    // Pa[i]=t*100+i;
    Pa[i]=make_pair(ll_to_p(ll,i),i);
  }
  sort(Pa.rbegin(),Pa.rend());

  int FORMIN= -L/2-1;
  int FORMAX=L/2+1;
  
  for (int dx =FORMIN; dx <=FORMAX ; dx+=1){//TODO//本当は+1にしたいんだが、計算量的に間に合わん
   for (int dy = FORMIN; dy <= FORMAX; dy+=1){
    double ret= calc_IAB(Pa,dy,dx)-0.00000005*(10+abs(dy)+abs(dx));//llの情報を持っているとき、dydxを叩くことで得られる情報量。（減少するエントロピーの量）
    if(ret>maxi){
      maxi=ret;
      y=dy;
      x=dx;
    }
   } 
  }

  // int c1=L/3/2;
  // int c2=L/2;//もうっか所あるので注意
  // int c3=L-1-L/3/2;
  // select_yx2(maxi,y,x,Pa,c1,c2);
  // select_yx2(maxi,y,x,Pa,c2,c1);
  // select_yx2(maxi,y,x,Pa,c2,c2);
  // select_yx2(maxi,y,x,Pa,c2,c3);
  // select_yx2(maxi,y,x,Pa,c3,c1);
    


  return make_pair(y,x);
}

void debug_print_ll(vector<vector<double>>& ll,int index){
  std::cout << std::fixed<<"#";
  for (int i = 0; i < N; i++){
    std::cout << std::setprecision(2) << ll_to_p(ll[index],i)<<" ";
  }
  std::cout << std::defaultfloat;
  cout<<endl;
}


bool solve(vector<vector<double>>& ll,int index,double fast_cut_rate){
  // cout<<"#-----------------------------"<<endl;
  double maxi=-1.0;
  // debug_print_ll(ll,index);
  for (int i = 0; i < N; i++){
    double t=ll_to_p(ll[index],i);
    if(t>maxi){
      maxi=t;
    }
  }

  if(maxi>=fast_cut_rate) {
    cout << "#"<<"end because maxi>=fast_cut_rate:"<<fast_cut_rate << endl;
    return true;
  }

  if(isnan(maxi)||isinf(maxi)||maxi<0.0){
    cout << "#"<<"isnan(maxi)" << endl;
    return true;
  }

  //yxの選択(max)
  P ret = select_yx(ll[index]);////3,125,000 //890,625

  int y=ret.first;
  int x=ret.second;

  //計測実行
  int m=Judge().measure(index,y,x);

  update_ll(ll,m,index,y,x);
  // debug_print_ll(ll,index);
  return false;
}

void select_pi_greedy(const vector<vector<double>>& ll,vector<pair<int,double>>& pi){

  vector<tuple<double,int,int>>pq;//選択確率、entrance,exit

  for (int entrance = 0; entrance < N; entrance++){
   for (int exit = 0; exit < N; exit++){
    double t=ll_to_p(ll[entrance],exit);
    pq.push_back(make_tuple(t,entrance,exit));
   } 
  }
  sort(pq.rbegin(),pq.rend());

  vector<bool>used_ent(N,false);
  vector<bool>used(N,false);
  for(auto v:pq){
    double p = get<0>(v);
    int ent= get<1>(v);
    int exi= get<2>(v);
    if(used[exi]==true)continue;
    if(used_ent[ent]==true)continue;
    used_ent[ent]=true;
    used[exi]=true;
    pi[ent]=make_pair(exi,p);

  }

}

int main(int argc, char *argv[]) {
  clock_t start_time = clock(); // 開始時刻
  clock_t TIME_LIMIT = 3800000;//TODO

  ios::sync_with_stdio(false);
  cin.tie(nullptr);
  // 10≤L≤50 60≤N≤100 S∈{i^2∣1≤i≤30} 
  cin >>L>>N>>S;
  // if(S>=400){
  //   S/=4;
  //   n_sample=4;
  // }
  EXIT_COR=vector<P> (N);
  KLmemo=vector<vector<double>>(1001,vector<double>(1001,-1));
  normalCDF_init(0,S);
  for (int i = 0; i < N; i++){
    int y,x;
    cin >>y>>x;
    EXIT_COR[i]=make_pair(y,x);
  }

    //パラメタ生成
  param = parameta();
  if(argc==3){
    // param.pw1=atoi(argv[1]);
    param.h1 =atoi(argv[1]);

    // param.pw2=atoi(argv[3]);
    param.h2 =atoi(argv[2]);
    // param.h00 =atoi(argv[5]);

    // param.acc_prob = atof(argv[6]);
  }else{
    // cerr<<"#not param. is it ok?"<<endl;
    // dprint("not param");
  }
  //配置
  Temp = temperature();
  Temp.init();
  Judge().print_temperature(Temp);


  //計測
  vector<vector<double>> ll(N, vector<double>(N+1,0));
  for (int i = 0; i < N; i++){//N番目には合計値を保存しておく。(高速化用)
   ll[i][N]=N;//*10^0 
  }


  int index=0;
  //////一週目////
  while (true) { // 時間の許す限り回す
    
    clock_t now_time = clock();
    if(now_time-start_time>TIME_LIMIT){
      cout<<"#end because timelimit"<<endl;
      break;
    }
    if(measure_count>=10000){
      cout<<"# end because count>=10000"<<endl;
      break;
    }
    bool ret = solve(ll,index,param.acc_prob);
    if(ret==true){
      index++;
      if(index>=N)break;
    }
  }

cout << "#"<<"----微調整----" << endl;
  //////微調整////
  vector<pair<int,double>>pi(N,make_pair(0,0.0));
  while(true){
    //系列得る
    select_pi_greedy(ll,pi);
    
    clock_t now_time = clock();
    if(now_time-start_time>TIME_LIMIT){
      cout<<"#end because timelimit"<<endl;
      break;
    }
    if(measure_count>=10000){
      cout<<"# end because count>=10000"<<endl;
      break;
    }
    //１　piの中で最も自身がないやつを探す、->大体外れなので、
    int index_weak=0;
    double mini=1.0;
    for (int i = 0; i < N; i++){
      if(mini>pi[i].second){
        mini=pi[i].second;
        index_weak=i;
      }
    }
    if(mini>=param.acc_prob2) {//もし受理できるならそのまま抜ける
      cout << "#"<<"end because mini>=param.acc_prob2!!!!!:" <<param.acc_prob2<< endl;
      break;
    }
    //2　それについて再調査して系列アップデート(fast cutなし（してしまうので）)
    solve(ll,index_weak,1.1);
    }

  //回答
  Judge().answer(pi);

  return 0;
}

Submission Info

Judge Result