纯模拟退火题解
jia_shengyuan
2021-12-01 21:46:01
## 前言
考场上T1脑抽花了1h，导致这题没时间观察性质，没发现交换方差，只能写了个朴素退火，$100\to40$。
这是一篇****纯模拟退火****的题解，相较于其他题解应该是很简洁的。
#### $40\ pts$ 做法
按题目要求每次随机一个数进行操作，模拟退火即可。
#### $80\ pts$ 做法
首先观察性质，$a_i=a_{i-1}+a_{i+1}-a_i$ 的实质就是交换相邻两项的差分，而通过若干次操作我们可以将原差分序列重新排列成任意我们想要的顺序。
观察到交换差分的性质后，无论是推样例还是数感，都不难感觉到方差最小时，差分应该先减后增（你可以认为是让尽可能多的数集中，从而减小方差）。
因此，我们可以先构造出一种差分先减后增的初始解，之后每次随机两个位置并交换差分，在此基础上再模拟退火，即可获得 $80\ pts$ 的好成绩。
#### $100\ pts$ 做法
我们可以把交换差分改为随机把一个差分插入到另一个地方。观察到许多插入操作的最优方案是****一定****的。例如对于差分序列 `3 1 0 2 4`，若要移动最左测的 $3$，则放到 $2,4$ 中间一定是最优的。
因此，我们每次可以随机一个值，并把它放到差分拐点的另一侧，使得差分序列仍然满足先减后增的性质，在此基础上模拟退火，并调出较好的参数，即可通过本题。
由于涉及到插入与删除，且有吸氧，我们不妨使用STL的 `deque` 简化代码。
我们可以用样例 $4$ 进行测试，以获得较好的参数。实际上我的代码在样例 $4$ 上的出解率也只有 $80\%$，然而交了几次都能过掉这题，希望能出hack数据。（刚交的题解发现有问题……然后就再交了一次，希望不会导致管理看到好几篇吧）
```cpp
#include <algorithm>
#include <climits>
#include <cstdio>
#include <random>
#include <deque>
#include <ctime>
#include <cmath>
#include <cassert>
#define ll long long
// #define err printf
#define err(...)
const int maxn = 1e4+9;
using namespace std;
typedef deque<int>::iterator iter;
int a[maxn],b[maxn],mini=INT_MAX,n;
deque<int> cf;//差分
iter lpt,rpt;
ll ans;
mt19937 rnd(time(0));
void Read(){
scanf("%d",&n);
for(int i=0; i<n; i++){
scanf("%d",&a[i]);
a[i]*=n; // 为防止出现小数，开始时把每一项乘以n即可
}
for(int i=0; i<n-1; i++){
b[i]=a[i+1]-a[i];
}
}
void Init(){
sort(b,b+n-1);
bool f=0;
// 构造初始解
for(int i=0; i<n-1; i++){
mini = min(mini,b[i]);
if(f) cf.push_front(b[i]);
else cf.push_back(b[i]);
f = !f;
}
// 计算最小值有几个，用于退火时判断当前数在拐点左侧还是右侧
f=0;
for(int i=0; i<n-1; i++){
if(cf[i]==mini){
if(!f){
lpt = cf.begin()+i;
rpt = cf.begin()+i;
f=1;
} else {
rpt = cf.begin()+i;
}
}
}
}
ll pf(ll x) { return x*x; }
// 计算方差。由于整体加减不影响方差，我们默认第一项是0
ll calc(){
ll sum=0,cur=0,pj,fc;
for(int solo : cf){
cur += solo;
sum += cur;
}
pj = sum/n;
fc = pf(pj);//a[0]的贡献
cur = 0;
for(int solo : cf){
cur += solo;
fc += pf(cur-pj);
}
return fc/n;
}
void SA(){
double T = 2e5;//实测这个初温有很好的出解率
ll now=ans,newans,de;
int rndsize = n-1-(rpt-lpt+1),val,len=rpt-lpt+1;
iter p1;
while(T>1e-5){
// 处理拐点区域
for(int i=0; i<n-1; i++){
if(cf[i] == mini){
lpt = cf.begin()+i;
rpt = lpt+len-1;
break;
}
}
deque<int> tmp = cf;
// 随机选一个数插入到别的地方
p1 = cf.begin()+(rnd()%rndsize);
if(p1 >= lpt){ //从右侧插入到左侧
p1 += len;
val = *p1;
cf.erase(p1);
iter now = cf.begin();
while(*now > val) now++;
cf.insert(now,val);
} else { // 从左侧插入到右侧
val = *p1;
cf.erase(p1);
iter now = cf.end();
while(*(now-1) > val) {
now--;
}
cf.insert(now,val);
}
// 更新答案
newans = calc(); de=newans-now;
if(de<0 || exp(-de/T)*UINT_MAX>rnd()) {
ans = min(ans,newans);
now = newans;
} else cf = tmp;
T*=0.97;
}
}
int main(){
Read();
Init();
ans = calc();
if(rpt-lpt+1 == n-1){ //如果差分全一样，直接输出防止诡异错误
printf("%lld",ans);
} else {
// 退火时尽可能把时间用满，增加出解率
clock_t s = clock();
SA();
clock_t len = clock()-s, limit = CLOCKS_PER_SEC*0.975;
while(clock()+len < limit) SA();
printf("%lld",ans);
}
return 0;
}
```