USACO5.2.3-fence5-
我也覺得好像算法是錯的。
但就是過掉了，算法就是逐漸加大精度的搜索。
先找出最優解的大致范圍，然后繼續搜索。
但的確是一種解決問題的思路吧！加上卡時間，應該能拿很多分的。
兩次遞降精度掃描
/**//*
USER:zyn19961
TASK:fence3
LANG:C++
*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<fstream>
#include<iostream>
#include<algorithm>
using namespace std;
//
#define MM(a,i) memset(a,i,sizeof(a))
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define DFOR(i,b,a) for(int i=(b);i>=(a);i--)
//
const int maxn=151;
const int INF=~0U>>2;
//
typedef long long int64;
//
ifstream fin("fence3.in");
ofstream fout("fence3.out");
class line
{public:int x1,x2,y1,y2;};
line L[maxn];
int minx=INF,maxx=0,miny=INF,maxy=0;
int xx,yy;
int main(){
int n;
double ans=double(INF);
fin>>n;
FOR(i,1,n){
fin>>L[i].x1>>L[i].y1>>L[i].x2>>L[i].y2;
L[i].x1*=10,L[i].y1*=10,L[i].x2*=10,L[i].y2*=10;
if(L[i].y1==L[i].y2)if(L[i].x1>L[i].x2)swap(L[i].x1,L[i].x2),swap(L[i].y1,L[i].y2);
if(L[i].x1==L[i].x2)if(L[i].y1>L[i].y2)swap(L[i].x1,L[i].x2),swap(L[i].y1,L[i].y2);
minx=min(minx,L[i].x1),miny=min(miny,L[i].y1),maxx=max(maxx,L[i].x2),maxy=max(maxy,L[i].y2);
}
for(int x=minx;x<=maxx;x+=16)
for(int y=miny;y<=maxy;y+=16){
double t=0.0;
FOR(i,1,n){
if(L[i].x1==L[i].x2){
if(L[i].y1>y)t+=sqrt((y-L[i].y1)*(y-L[i].y1)+(x-L[i].x1)*(x-L[i].x1));
else if(L[i].y2<y)t+=sqrt((y-L[i].y2)*(y-L[i].y2)+(x-L[i].x1)*(x-L[i].x1));
else t+=abs(x-L[i].x1);
}
else{
if(L[i].x1>x)t+=sqrt((x-L[i].x1)*(x-L[i].x1)+(y-L[i].y1)*(y-L[i].y1));
else if(L[i].x2<x)t+=sqrt((x-L[i].x2)*(x-L[i].x2)+(y-L[i].y1)*(y-L[i].y1));
else t+=abs(y-L[i].y1);
}
}
if(t<ans)ans=t,xx=x,yy=y;
}
minx=xx-16,maxx=xx+16,miny=yy-16,maxy=yy+16;
FOR(x,minx,maxx)
FOR(y,miny,maxy){
double t=0.0;
FOR(i,1,n){
if(L[i].x1==L[i].x2){
if(L[i].y1>y)t+=sqrt((y-L[i].y1)*(y-L[i].y1)+(x-L[i].x1)*(x-L[i].x1));
else if(L[i].y2<y)t+=sqrt((y-L[i].y2)*(y-L[i].y2)+(x-L[i].x1)*(x-L[i].x1));
else t+=abs(x-L[i].x1);
}
else{
if(L[i].x1>x)t+=sqrt((x-L[i].x1)*(x-L[i].x1)+(y-L[i].y1)*(y-L[i].y1));
else if(L[i].x2<x)t+=sqrt((x-L[i].x2)*(x-L[i].x2)+(y-L[i].y1)*(y-L[i].y1));
else t+=abs(y-L[i].y1);
}
}
if(t<ans)ans=t,xx=x,yy=y;
}
int q;
if(ans-floor(ans)>0.5)q=int(ceil(ans));else q=int(floor(ans));
fout<<xx/10<<"."<<xx%10<<" ";
fout<<yy/10<<"."<<yy%10<<" ";
fout<<q/10<<"."<<q%10<<"\n";
fin.close();
fout.close();
return 0;
}
據說正解是模擬退火，等待學習。
其實這道題，本來就是考察近似算法，所以一開始的那種方法也是可以接受的。（不過效果更差罷了）
學習了一下模擬退火，主要懂了思想，就是有一定的概率接受較次的解，以防止陷入局部最優之中。
然后溫度T，單調遞減，目的是使這個概率逐漸降低。
但實現時非常不成功，需要大量的調試+枚舉算法對拍。
我的經歷：1.一次只向一個方向移動，x或y
2.設計移動距離關于T的函數比較困難，多項式顯然錯誤，因為T是遞減的冪函數，
概率函數是指數函數，我嘗試地試了對數函數，結果錯解無數。。。
3.T不一定要降到一個指定的值啊，只要周圍沒有更優解就結束。
4.其實我寫的東西還是精度遞降的形式（常數遞降），每個精度也不過嘗試一個常數次數罷了，只是加上了接受較劣解的函數而已。
——可以見得，設計各種函數、常數都是有一定技巧的——
不正版的模擬退火
/**//*
USER:zyn19961
TASK:fence3
LANG:C++
*/
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<fstream>
#include<iostream>
#include<algorithm>
using namespace std;
//
#define MM(a,i) memset(a,i,sizeof(a))
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define DFOR(i,b,a) for(int i=(b);i>=(a);i--)
//
const int maxn=151;
const int INF=~0U>>2;
//
typedef long long int64;
//
ifstream fin("fence3.in");
ofstream fout("fence3.out");
class line{public:int x1,x2,y1,y2;};
line L[maxn];
int minx=INF,maxx=0,miny=INF,maxy=0;
int xx,yy;
int n;
double come(int x,int y);
int main(){
srand(unsigned(time(NULL)));
double ans=double(INF);
fin>>n;
FOR(i,1,n){
fin>>L[i].x1>>L[i].y1>>L[i].x2>>L[i].y2;
L[i].x1*=10,L[i].y1*=10,L[i].x2*=10,L[i].y2*=10;
if(L[i].y1==L[i].y2)if(L[i].x1>L[i].x2)swap(L[i].x1,L[i].x2),swap(L[i].y1,L[i].y2);
if(L[i].x1==L[i].x2)if(L[i].y1>L[i].y2)swap(L[i].x1,L[i].x2),swap(L[i].y1,L[i].y2);
minx=min(minx,L[i].x1),miny=min(miny,L[i].y1),maxx=max(maxx,L[i].x2),maxy=max(maxy,L[i].y2);
}
int x,y;
int nowx=(maxx+minx)/2,nowy=(maxy+miny)/2;
double tt=come(nowx,nowy);
double const INIT_T=1000,RATE=0.95,FINAL_T=1E-10,loga=10;const int tim=5;
double T=INIT_T;
for(int i=10;i;i--){
for(int k=1;k<=200;k++){
if(tt<ans){ans=tt,xx=nowx,yy=nowy;}
int xxxx=rand()%4;
if(xxxx==1)x=nowx+i,y=nowy;
if(xxxx==2)x=nowx-i,y=nowy;
if(xxxx==3)x=nowx,y=nowy+i;
if(xxxx==0)x=nowx,y=nowy-i;
if(x>=minx&&x<=maxx&&y>=miny&&y<=maxy){
double t=come(x,y);
if(t<tt){tt=t,nowx=x,nowy=y,T*=RATE;}
else {
double rr=rand()%100000/100000.0;
double p=exp((tt-t)/T);
if(p>rr){tt=t,nowx=x,nowy=y,T*=RATE;}
}
}
}
}
int q;
if(ans-floor(ans)>0.5)q=int(ceil(ans));else q=int(floor(ans));
fout<<xx/10<<"."<<xx%10<<" ";
fout<<yy/10<<"."<<yy%10<<" ";
fout<<q/10<<"."<<q%10<<"\n";
fin.close();
fout.close();
return 0;
}
double come(int x,int y){
double t=0.0;
FOR(i,1,n){
if(L[i].x1==L[i].x2){
if(L[i].y1>y)t+=sqrt((y-L[i].y1)*(y-L[i].y1)+(x-L[i].x1)*(x-L[i].x1));
else if(L[i].y2<y)t+=sqrt((y-L[i].y2)*(y-L[i].y2)+(x-L[i].x1)*(x-L[i].x1));
else t+=abs(x-L[i].x1);
}
else{
if(L[i].x1>x)t+=sqrt((x-L[i].x1)*(x-L[i].x1)+(y-L[i].y1)*(y-L[i].y1));
else if(L[i].x2<x)t+=sqrt((x-L[i].x2)*(x-L[i].x2)+(y-L[i].y1)*(y-L[i].y1));
else t+=abs(y-L[i].y1);
}
}
return t;
}
據說正版代碼，部分代碼求解釋
#include <fstream>
#include <ctime>
#include <cmath>
#include <iomanip>
#define V 150
#define sqr(q) ((q)*(q))
#define ex(a,b) ({long c;c=(a);(a)=(b);(b)=c;})
#define ou(x,y,x1,y1) (sqrt(sqr((x)-(x1))+sqr((y)-(y1))))
using namespace std;
struct re{
long x1,y1,x2,y2;
}l[V];
ifstream fin("fence3.in");ofstream fout("fence3.out");
long n=0;
double dis(double x,double y,long k){
if(l[k].x1==l[k].x2){
if(y<l[k].y1)return ou(x,y,l[k].x1,l[k].y1);
if(y>l[k].y2)return ou(x,y,l[k].x2,l[k].y2);
return fabs(x-l[k].x1);
}else{
if(x<l[k].x1)return ou(x,y,l[k].x1,l[k].y1);
if(x>l[k].x2)return ou(x,y,l[k].x2,l[k].y2);
return fabs(y-l[k].y1);
}
}
int main(int argc,char** argv){
srand(size_t(time(NULL)));
fin >>n;
double x=rand()%100,y=rand()%100;double t=100;double best=0;
for(long i=1;i<=n;i++){
fin >>l[i].x1>>l[i].y1>>l[i].x2>>l[i].y2;
if(l[i].x1>l[i].x2)ex(l[i].x1,l[i].x2);
if(l[i].y1>l[i].y2)ex(l[i].y1,l[i].y2);
best+=dis(x,y,i);
}
long d=31;double temp1,temp2;
while(t>10e-3){
for(long i=1;i<=500;i++){
double x1,y1;
x1=t*(double(rand())/double(RAND_MAX))*(2*(rand()%2)-1);
y1=sqrt(sqr(t)-sqr(x1))*(2*(rand()%2)-1)+y;//這是神馬東西？
x1+=x;
double temp=0;
for(long k=1;k<=n;k++)temp+=dis(x1,y1,k);
double p=temp-best,yy=0;
if(p<0){
yy=1;
temp1=x1;temp2=y1;best=temp;
}else
yy=exp(-p/t);
double q=double(rand())/double(RAND_MAX);//這個概率咋么計算出來的？
if(q<yy)
{x=x1;
y=y1;}
}
d++;
t=t/(log10(d)/log10(20));
}
fin.close();
fout <<fixed<<setprecision(1);
fout<<temp1<<' '<<temp2<<' '<<best<<endl;
return 0;
}