kaggle Santa 2022 赛后总结
比赛表述:
https://www.kaggle.com/competitions/santa-2022/overview/description
2023年春节前参加了kaggle santa-2022 比赛,最后获得了铜牌成绩。这次比赛我几乎全程参与了整个过程,所以希望分享给大家。
初衷
从比赛描述来看,数据量小,参与门槛不高;时间跨度为1个月,占用时间不会太多;奖金丰厚,吸引不少高手参与。
过程
Follow了讨论区,每天都会查看讨论区里大家发布的想法和评论,也参与了一些讨论。同时在code区里download代码,将有用的高效率的function 合并到自己的代码里,我认为discuss和code 是获得奖牌的必要条件。
查看同类型的往年解决方案,发现各年来的解决方案中确实存在些共性,那就是使用LKH作为提升分数的关键工具,这也是我最终使用的方法。
我的解决方案
将原问题转换为近似的TSP问题
图片中一个像素点具备x,y 坐标,以及颜色信息rgb,input 作为 城市的X,Y,Z坐标
cost 近似为 城市间的权重 即 sqrt(abs(ax-bx)+abs(ay-by))+3*color_cost
生成TSP 文件
NAME: distances
TYPE: TSP
COMMENT: TSP
DIMENSION: 66049
EDGE_WEIGHT_TYPE: SPECIAL
EDGE_WEIGHT_FORMAT: FUNCTION
NODE_COORD_TYPE: THREED_COORDS
NODE_COORD_SECTION
添加约束
两类:由初始位置引起的约束 和 以机械臂配置引起的约束
具体实现
1.起始位置:64 0;-32 0;-16 0;-8 0;-4 0;-2 0;-1 0;-1 0
仅能在y轴方向移动
if((ax==0&&ay==0)||(bx==0&&by==0))
{
if(((ax==0&&abs(ay)==1))||(bx==0&&abs(by)==1))
return 0;
else
return 800*scal+dx+dy;
}
每次移动是每个机械臂仅能移动一个单位,近似为在x,y方向只能移动一个单位
if(dx>1||dy>1)
return 1000*scal+dx+dy;
特殊约束
仅能在(+-128,+-128)穿越四个象限
int get_xq(int x,int y)
{
if(x>=0&&y>0)
return 1;
if(x<0&&y>=0)
return 2;
if(x<0&&&y<0)
return 3;
if(x>=0&&y<=0)
return 4;
return 0;
}
if(ax==0&&ay==128&&abs(bx)==1&&by==128)
return 0;
if(ax==0&&ay==-128&&abs(bx)==1&&by==-128)
return 0;
if(ax==128&&ay==0&&bx==128&&abs(by)==1)
return 0;
if(ax==-128&&ay==0&&ax==-128&&abs(by)==1)
return 0;
if(bx==0&&by==128&&abs(ax)==1&&ay==128)
return 0;
if(bx==0&&by==-128&&abs(ax)==1&&by==-128)
return 0;
if(bx==128&&by==0&&bx==128&&abs(ay)==1)
return 0;
if(bx==-128&&by==0&&ax==-128&&abs(ay)==1)
return 0;
if(get_xq(ax,ay)!=get_xq(bx,by))
return 400*scal;
先访问abs(max(x,y))<64 的内部再访问外部,仅能通过(+-64,0),(0,+-64)穿越区域
#define MAX(a, b) (((a) > (b)) ? (a) : (b))
if(ax==0&&abs(ay)==64&&bx==0&&abs(ay-by)==1)
return 0;
if(bx==0&&abs(by)==64&&ax==0&&abs(ay-by)==1)
return 0;
if(ay==0&&abs(ax)==64&&by==0&&abs(bx-ax)==1)
return 0;
if(by==0&&abs(bx)==64&&bx==0&&abs(ax-bx)==1)
return 0;
int d_a = MAX(abs(ax),abs(ay));
int d_b = MAX(abs(bx),abs(by));
if(d_a<64&&d_b>=64)return 200*scal+dx;
if(d_b<64&&d_a>=64) return 200*scal+dy;
其他分界区域
int bound[] = {8,16,24,32,40,48,64,72,80,96,112,120};
int limit = 128;
for (int i = 0;i <sizeof(bound)/4;i++)
{
if(
(abs(ay)>bound[i]&&abs(by)<=bound[i])
|| (abs(by)>bound[i]&&abs(ay)<=bound[i])||(abs(ax)>bound[i]&&abs(bx)<=bound[i])
|| (abs(bx)>bound[i]&&abs(ax)<=bound[i])
){
dx+=0.05;dy+=0.05;
}
}
禁止访问的正方形边界
if(d_a==d_b&&(d_a == 1||d_a==4||d_a==8))
{
return 300*scal+dx+dy;
}
获取结果
运行LKH 并获取结果,tour 文件,LKH会单线程运行,完整的一个run大概需要32h
从tour文件中获取path的方法:
采用贪心方法,关键函数如下:
-
get_neighbors(config,p),上个点的机械臂配置为config,下个点的坐标为p,可以从该函数获取p的可能配置,返回为list结构,list的每个单元为np.array()
这个函数是从discuss 复制过来,通过hash大大节约了查询时间,将原有的寻路时间从30分钟降低到了5分钟左右
L = 8
lengths = np.array([int(2**max(L-2-i,0)) for i in range(L)])
# pre-computes the "squares"
squares = np.zeros((len(lengths), (lengths[0] * 8), 2), dtype = 'int')
for i, length in enumerate(lengths):
r = np.arange((-length),(length + 1))
h = np.array(np.meshgrid(r,r)).reshape((2,-1)).T # full square
h = h[np.maximum(np.abs(h[:,0]),np.abs(h[:,1])) == length] # boundary of square
squares[i, :len(h)] = h
# pre-computes all combinations
nbhd_index = np.array(list(product(range(3),repeat=L)))
nbhd_index = nbhd_index[1:,:] # deleted neighborhood
@nb.njit()
def get_neighbors(config,p):
nbhds= np.zeros(((3**8 - 1), len(config), 2), dtype = 'int')
for i in range(L):
x = config[i,:].reshape(1, -1)
h = squares[i, :(lengths[i] * 8)]
x_n = np.concatenate((x,h[np.abs(h-x).sum(axis=1)==1]))
if(len(x_n)!=3):
print (config)
assert len(x_n)==3 # exactly 3 neighbors for each link
nbhds[:,i] = x_n[nbhd_index[:,i]]
res = []
for i in nbhds:
if((get_position(i) == p).all() and (i !=config).any()):
res.append(i)
return res
# return nbhds
reverse path
当到达下一点时,发现1次移动到达不了,则修改上个点配置试试看,还找不到,就修改上上个点的配置试试,一直往前去找。
该方法是个笨方法,也是贪心算法的一点改进。
结果优化
这部分主要来自讨论区的方法,
通过2 opt search ,最后实践结果大概降低了cost 5~10分左右
这段程序使用numba 的加速及并行大大降低了搜索时间
from numba import njit, prange
@nb.njit(fastmath=True,nopython=True,parallel=True)
def opt2(opt2path):
for j in prange(1,len(opt2path)):
for i in prange(len(opt2path)-2-j):
if(j==1):
current_cost= total_cost([opt2path[i], opt2path[i+1],opt2path[i+2],opt2path[i+3]], image)
new_cost= total_cost([opt2path[i], opt2path[i+2],opt2path[i+1],opt2path[i+3]], image)
# if(current_cost<5):
# continue
if(new_cost<current_cost):
opt2path[i+1],opt2path[i+j+1] = opt2path[i+j+1].copy(),opt2path[i+1].copy()
break
else:
distance1 = np.abs(get_position(opt2path[i])-get_position(opt2path[i+j+1])).max()
distance2 = np.abs(get_position(opt2path[i+1])-get_position(opt2path[i+j+2])).max()
distance3 = np.abs(get_position(opt2path[i+2])-get_position(opt2path[i+j+1])).max()
distance4 = np.abs(get_position(opt2path[i+1])-get_position(opt2path[i+j])).max()
if(distance1>2 or distance2>2 or distance3>2 or distance4>2):
continue
current_cost = total_cost([opt2path[i], opt2path[i+1],opt2path[i+2]], image)
current_cost += total_cost([opt2path[i+j],opt2path[i+j+1], opt2path[i+j+2]], image)
new_cost = total_cost([opt2path[i], opt2path[i+j+1],opt2path[i+2]], image)
new_cost+= total_cost([opt2path[i+j],opt2path[i+1],opt2path[i+j+2]],image)
# if(current_cost <6):
# continue
if(new_cost<current_cost):
opt2path[i+1],opt2path[i+j+1] = opt2path[i+j+1].copy(),opt2path[i+1].copy()
break
else:
e1 = get_neighbors(opt2path[i],get_position(opt2path[i+j+1]))
e2 = get_neighbors(opt2path[i+j+2],get_position(opt2path[i+1]))
if(len(e1)>0 and len(e2)>0):
min_cost = current_cost
for t1 in e1:
diffs = np.abs(np.asarray(t1) - np.asarray(opt2path[i+2])).sum(axis=1)
if(diffs.max()>1):
continue
for t2 in e2:
if(j==1):
diffs = np.abs(np.asarray(t2) - np.asarray(t1)).sum(axis=1)
if(diffs.max()>1):
continue
new_cost = total_cost([opt2path[i],t1,t2,opt2path[i+j+2]], image)
else:
diffs = np.abs(np.asarray(t2) - np.asarray(opt2path[i+j])).sum(axis=1)
if(diffs.max()>1):
continue
new_cost = total_cost([opt2path[i],t1,opt2path[i+2]],image)
new_cost+= total_cost([opt2path[i+j],t2,opt2path[i+j+2]],image)
if(new_cost<min_cost):
min_cost = new_cost
new_t1 = t1
new_t2 = t2
if(min_cost <current_cost):
opt2path[i+1],opt2path[i+j+1] = new_t1.copy(),new_t2.copy()
print(total_cost(opt2path,image))
break
结果
74473.520276204200
获得铜牌,排名 87/875
教训:
-
添加约束过多导致 low bound 偏高,下界和讨论区中描述不符,导致从根本上无法获取最优结果。
-
获取路径时greedy search 可以改进为beam search 进行优化
方法优点:
-
根本方法 LKH+ 约束 正确
-
从修改LKH distance_special.c 入手,大大降低了调试难度。