如何在 python 中打包球体?

How to pack spheres in python?

我正在尝试使用 python 对正方形中大小均匀的随机封闭堆积球体进行建模。 球体不应重叠但我不知道该怎么做

我到目前为止:

代码:

import random, math, pylab

def show_conf(L, sigma, title, fname):
    pylab.axes()
    for [x, y] in L:
        for ix in range(-1, 2):
            for iy in range(-1, 2):
                cir = pylab.Circle((x + ix, y + iy), radius=sigma,  fc='r')
                pylab.gca().add_patch(cir)
    pylab.axis('scaled')
    pylab.xlabel('eixo x')
    pylab.ylabel('eixo y')
    pylab.title(title)
    pylab.axis([0.0, 1.0, 0.0, 1.0])
    pylab.savefig(fname)
    pylab.close()

L = []
N = 8 ** 2

for i in range(N):
    posx = float(random.uniform(0, 1))
    posy = float(random.uniform(0, 1))
    L.append([posx, posy])

print L

N = 8 ** 2
eta = 0.3
sigma = math.sqrt(eta / (N * math.pi))
Q = 20
ltilde = 5*sigma

N_sqrt = int(math.sqrt(N) + 0.5)


titulo1 = '$N=$'+str(N)+', $\eta =$'+str(eta)
nome1 = 'inicial'+'_N_'+str(N) + '_eta_'+str(eta) + '.png'
show_conf(L, sigma, titulo1, nome1)

这是一个非常难的问题(可能 np-hard)。应该有很多可用的资源。

在我介绍一些更通用的方法之前,请查看 this wikipedia-site 以了解当前最著名的 packing-patterns 的概述,其中 N 个(正方形中的 N 个圆圈)。

你很幸运,在 python 中有一个现有的 circle-packing 实现(heuristic!),它在很大程度上基于现代 optimization-theory (difference of convex-functions + Concave-convex-procedure).

  • 使用的方法描述here(学术论文&link软件;2016!)
  • 使用的软件包是here
    • 有一个带有 circle_packing.py 的示例目录(与输出一起发布在下面)
  • 以下示例也适用于不同形状的圆

上面的例子software-package(沈新月的例子)

__author__ = 'Xinyue'
from cvxpy import *
import numpy as np
import matplotlib.pyplot as plt
import dccp

n = 10
r = np.linspace(1,5,n)

c = Variable(n,2)
constr = []
for i in range(n-1):
    for j in range(i+1,n):
        constr.append(norm(c[i,:]-c[j,:])>=r[i]+r[j])
prob = Problem(Minimize(max_entries(max_entries(abs(c),axis=1)+r)), constr)
#prob = Problem(Minimize(max_entries(normInf(c,axis=1)+r)), constr)
prob.solve(method = 'dccp', ccp_times = 1)

l = max_entries(max_entries(abs(c),axis=1)+r).value*2
pi = np.pi
ratio = pi*sum_entries(square(r)).value/square(l).value
print "ratio =", ratio
print prob.status

# plot
plt.figure(figsize=(5,5))
circ = np.linspace(0,2*pi)
x_border = [-l/2, l/2, l/2, -l/2, -l/2]
y_border = [-l/2, -l/2, l/2, l/2, -l/2]
for i in xrange(n):
    plt.plot(c[i,0].value+r[i]*np.cos(circ),c[i,1].value+r[i]*np.sin(circ),'b')
plt.plot(x_border,y_border,'g')
plt.axes().set_aspect('equal')
plt.xlim([-l/2,l/2])
plt.ylim([-l/2,l/2])
plt.show()

输出

您的任务修改:equal-sized 个圈子

只需替换:

r = np.linspace(1,5,n)

与:

r = [1 for i in range(n)]

输出

64 个圆圈的有趣示例(这需要一些时间!)

如果您想要更新版本的@leopold.talirz 解决方案,我建议您使用以下方法:

from cvxpy import *
import numpy as np
import matplotlib.pyplot as plt
import dccp

n = 10
r = np.linspace(1,5,n)

c = Variable(shape=(n,2))
constr = []
for i in range(n-1):
    for j in range(i+1,n):
        constr.append(norm(c[i,:]-c[j,:])>=r[i]+r[j])
prob = Problem(Minimize(max(max(abs(c),axis=1)+r)), constr)
#prob = Problem(Minimize(max_entries(normInf(c,axis=1)+r)), constr)
prob.solve(method = 'dccp', ccp_times = 1)

l = max(max(abs(c),axis=1)+r).value*2
pi = np.pi
ratio = pi*sum(square(r)).value/square(l).value
print("ratio =", ratio)
print(prob.status)

# plot
plt.figure(figsize=(5,5))
circ = np.linspace(0,2*pi)
x_border = [-l/2, l/2, l/2, -l/2, -l/2]
y_border = [-l/2, -l/2, l/2, l/2, -l/2]
for i in range(n):
    plt.plot(c[i,0].value+r[i]*np.cos(circ),c[i,1].value+r[i]*np.sin(circ),'b')
plt.plot(x_border,y_border,'g')
plt.axes().set_aspect('equal')
plt.xlim([-l/2,l/2])
plt.ylim([-l/2,l/2])
plt.show()