生成多個互不重疊的不同半徑圓

https://www.zhihu.com/question/53012468

固定區域，生成所有圓，統計圓的面積，最後面積最大的就是答案。

import numpy as np
from scipy.optimize import minimize
from scipy.spatial.distance import pdist
import matplotlib.pyplot as plt
import os


def obj_grad_round(x, radii, pair_index_to_index_pair, repulsion_scale, attraction_scale):
    n = len(radii)
    centroids = np.reshape(x, (n, 2))
    circles_pdist_sqrd = pdist(centroids, 'sqeuclidean')
    obj = 0
    grad = np.zeros((n, 2))
    for pidx, dist_sqrd in enumerate(circles_pdist_sqrd):
        i, j = pair_index_to_index_pair[pidx]
        touching_dist_sqrd = (radii[i] + radii[j]) ** 2
        grad_inc = 4 * (dist_sqrd - touching_dist_sqrd) * (centroids[i, :] - centroids[j, :])
        scale = repulsion_scale if dist_sqrd < touching_dist_sqrd else attraction_scale
        obj += scale * (dist_sqrd - touching_dist_sqrd) ** 2
        grad[i, :] += scale * grad_inc
        grad[j, :] -= scale * grad_inc

    return obj, grad.flatten()


def obj_grad_square(x, radii, pair_index_to_index_pair, repulsion_scale, attraction_scale):
    n = len(radii)
    centroids = np.reshape(x, (n, 2))
    circles_pdist_sqrd = pdist(centroids, 'sqeuclidean')
    obj = 0
    grad = np.zeros((n, 2))
    for pidx, dist_sqrd in enumerate(circles_pdist_sqrd):
        i, j = pair_index_to_index_pair[pidx]
        touching_dist_sqrd = (radii[i] + radii[j]) ** 2
        diff = centroids[i, :] - centroids[j, :]
        if dist_sqrd < touching_dist_sqrd:
            # repulsion
            obj += repulsion_scale * (dist_sqrd - touching_dist_sqrd) ** 2
            grad_inc = 4 * (dist_sqrd - touching_dist_sqrd) * diff
            grad[i, :] += repulsion_scale * grad_inc
            grad[j, :] -= repulsion_scale * grad_inc
        else:
            # attraction
            diff_sqrd = diff ** 2
            if max(diff_sqrd) >= touching_dist_sqrd:
                # do not overlap after projection, needs attraction
                if diff_sqrd[0] >= diff_sqrd[1]:
                    obj_inc = (diff_sqrd[0] - touching_dist_sqrd) ** 2
                    grad_inc = np.array([4 * (diff_sqrd[0] - touching_dist_sqrd) * diff[0], 0])
                else:
                    obj_inc = (diff_sqrd[1] - touching_dist_sqrd) ** 2
                    grad_inc = np.array([0, 4 * (diff_sqrd[1] - touching_dist_sqrd) * diff[1]])

                obj += attraction_scale * obj_inc
                grad[i, :] += attraction_scale * grad_inc
                grad[j, :] -= attraction_scale * grad_inc

    return obj, grad.flatten()


def find_layout(radii, shape='square', repulsion_scale=1e6, attraction_scale=1, attraction_decay=0.3, overlapping_tol=0., seed=None, verbose=False):
    assert shape == 'round' or shape == 'square'
    obj_grad = obj_grad_round if shape == 'round' else obj_grad_square
    n = len(radii)
    n_pairs = int(n * (n-1) / 2)
    pair_index_to_index_pair = [None] * n_pairs
    pidx = 0
    for i in range(n-1):
        for j in range(i+1, n):
            pair_index_to_index_pair[pidx] = i, j
            pidx += 1

    if seed:
        np.random.seed(seed)

    if verbose:
        options = {'disp': True}
        history_folder = 'history_%s' % shape
        if not os.path.exists(history_folder):
            os.mkdir(history_folder)
    else:
        options = None

    overlapped = True
    n_trials = 1
    x0 = np.random.randn(2 * n)

    while overlapped:
        if verbose:
            print('repulsion scale = %g\tattraction scale = %g' % (repulsion_scale, attraction_scale))
        opt_result = minimize(obj_grad, x0, args=(radii, pair_index_to_index_pair, repulsion_scale, attraction_scale), jac=True, options=options)
        centroids = np.reshape(opt_result.x, (n, 2))
        if verbose:
            if not opt_result.success:
                print(opt_result.message)
            draw_layout(centroids, radii, '%s/%02d_attraction=%g.png' % (history_folder, n_trials, attraction_scale))
            np.savetxt('%s/%02d_attraction=%g.dat' % (history_folder, n_trials, attraction_scale), centroids, fmt='%g')
        overlapped = check_overlapping(centroids, radii, overlapping_tol, verbose)
        if overlapped:
            attraction_scale *= attraction_decay
            n_trials += 1
            x0 = opt_result.x
        else:
            overlapped = False

    return centroids


def check_overlapping(centroids, radii, tol, verbose=False):
    n = len(radii)
    for i in range(n-1):
        for j in range(i+1, n):
            dij = np.linalg.norm(centroids[i, :] - centroids[j, :])
            rhs = radii[i] + radii[j] - tol
            if dij < rhs:
                if verbose:
                    print('Circle {i} and circle {j} overlap ({dij} = dist(C_{i}, C_{j}) < r_{i} + r_{j} - tol = {rhs}).'.format(i=i, j=j, dij=dij, rhs=rhs))
                return True
    return False


def draw_layout(centroids, radii, not_show_and_save_to=None):
    fig, ax = plt.subplots()
    for centroid, radius in zip(centroids, radii):
        circle = plt.Circle(centroid, radius)
        ax.add_artist(circle)
    x = centroids[:, 0]
    y = centroids[:, 1]
    xmin = min(x - radii)
    xmax = max(x + radii)
    ymin = min(y - radii)
    ymax = max(y + radii)
    plt.axis([xmin, xmax, ymin, ymax])
    ax.set_aspect('equal')
    if not_show_and_save_to:
        plt.savefig(not_show_and_save_to)
    else:
        plt.show()
    plt.close()


if __name__ == '__main__':
    seed = 1
    np.random.seed(seed)
    n_circles = 100
    radii = np.random.rand(n_circles) / 10 + 0.1    # [0.1, 0.2)
    overlapping_tol = 1e-4
    shape = 'square'  # 'round' or 'square'
    verbose = False
    layout = find_layout(radii, shape, overlapping_tol=overlapping_tol, verbose=verbose)
    draw_layout(layout, radii, 'final_%s.png' % shape)
    np.savetxt('radii.dat', radii, fmt='%g')
    np.savetxt('final_%s.dat' % shape, layout, fmt='%g')