二维 & 三维费马点问题【没有理论研究,只有代码研究】

模拟退火什么的我只知道一个大概思想,然后会很矬的去实现它。近似算法真心伤不起啊唉。
例题:poj 2420 【二维多边形费马点】 , uva 5102 【福州现场赛E题,四边形费马点,有结论,所以模拟退火应该不算正解】, ZOJ 2928 【三维费马点】。

模拟退火神马的就是一个迭代过程,在每次迭代中向所有可能的方向试探,然后选取其中最优的一个作为下一次的迭代起点;为了让迭代能够终止,设定了一个步长,每次迭代之后步长乘以一个小于1的系数来减小步长。
迭代方式是相同的,就是向所有可能的方向试探,再按规则模拟计算结果,比如二维点就是向8个方向试探之后再计算到所有点的距离和,而三维点就是26个方向了。
初始步长作为上界,这个该怎么设定我一直不太清楚,我反正是用最大可能步长作为上界的;初始迭代点,随手设置一个吧,我用的是所有点坐标的算术平均数,也可以用原点之类,但是我觉得无论选择啥点都不能保证绝对的正确性的吧(否则,“图灵奖”。。。);迭代步长系数,好像我的写法比较矬?导致我必须把系数设到0.993以上才能过题,真是弱爆了;至于精度下界就设为1e-10就好了,反正不要钱,哈。

模拟退火的一大弊端在于,如果进入了某个次优解点的邻域,且此时步长已经很小的话,那很有可能就会陷入这个次优解邻域的漩涡里(俗话说叫掉坑里?好像有专业名词叫鞍点或者涡点什么的?我数学不太好不好意思。。),然后再也跑不出来。所以近似算法毕竟是“近似”算法,想避免掉入坑里,只能去优化写法或者找出答案的某种规律性,来分段或者剪枝迭代什么的。。我是玩不转这个了,如上所述的三题,我加起来交了五六十次了吧?总之是没有找到一组又能ac又快的好参数。唉。
我若爆了。。

附代码吧,以2928题为例:

#include <iostream>
#include 
<cstdio>
#include 
<cstring>
#include 
<cmath>
#define maxn 105
const double eps = 1e-12;
int comp(double x)
{
    
if(fabs(x) < eps)
        
return 0;
    
else if(x < -eps)
        
return -1;
    
else
        
return 1;
}
using namespace std;
struct point
{
    
double x,y,z;
    point(
double a = 0,double b = 0,double c = 0):x(a),y(b),z(c){}
    point 
operator - (const point p)
    {
        
return point(x - p.x,y - p.y,z - p.z);
    }
    point 
operator + (const point p)
    {
        
return point(x + p.x,y + p.y,z + p.z);
    }
    
double norm()
    {
        
return sqrt(x * x + y * y + z * z);
    }
    point 
operator *(double a)
    {
        
return point(x * a,y * a,z * a);
    }
}p[maxn],fermat,dir[
30];
int n;
int main()
{
    
while(scanf("%d",&n) == 1)
    {
        
for(int i = 0;i < n;i++)
            scanf(
"%lf %lf %lf",&p[i].x,&p[i].y,&p[i].z);
        
int cnt = 0;
        
for(int i=-1;i<=1;i++)
            
for(int j=-1;j<=1;j++)
                
for(int k = -1;k <= 1;k++)
                    dir[cnt
++= point(i,j,k);
        
double step = 100;
        point now,ans;
        
double temp = 1e10,calc;
        point sel;
        
for(int i = 0;i < n;i++)
            ans 
= ans + p[i];
        ans 
= ans * (1.0 / (double)n);
        
while(step > eps)
        {
            sel 
= ans;
            
for(int i = 0;i < cnt;i++)
            {
                calc 
= 0.0;
                now 
= ans + dir[i] * step;
                
for(int idx = 0;idx < n;idx++)
                    calc 
+= (p[idx] - now).norm();
                
if(comp(calc - temp) < 0)
                {
                    sel 
= now;
                    temp 
= calc;
                }
            }
            ans 
= sel;
            step 
*= 0.993;
        }
        printf(
"%.3f %.3f %.3f\n",ans.x,ans.y,ans.z);
    }
}

posted on 2011-09-09 20:23 BUPT-[aswmtjdsj] @ Penalty 阅读(542) 评论(0)  编辑 收藏 引用 所属分类: 计算几何ZOJ Solution ReportPOJ Solution Report UVAlive Solution Report


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理


<2011年9月>
28293031123
45678910
11121314151617
18192021222324
2526272829301
2345678

导航

统计

常用链接

留言簿(1)

随笔分类(150)

随笔档案(71)

搜索

积分与排名

最新评论

阅读排行榜

评论排行榜