Iterate over coordinates and correct constraint violations

Question

I am creating an algorithm to take a series of labelled points placed randomly and move them until they fulfil a set of constraints.

The constraints have been pre-computed and are each a lower and upper bound of the distance between two points. There is definitely a solution that satisfies all constraints.

The procedure is as follows:

• Every pair of points is checked. If the distance violates the constraint, the points are moved an equal amount apart/together to be at a uniformly-random distance between the bounds.
• After each pair has been checked, the number of violations is recomputed (as some pairs will be moved having been fixed) and if it is greater than zero, every pair is iterated over again.
• This process is repeated up to a max number of times. If violations are still present, the coordinates are re-randomised and the whole process starts again.

The code below works but it is slow when I have, say, 2000 points and each has a constraint with each other.

I am looking for input on how to speed up the code below. The procedure outlined above must be followed, but any algorithm that implements it is allowed.

I have a few ideas for improving it, including:

• Don't take the square roots: compare square distances instead.
• Have a list of pairs that need to be 'checked' rather than checking the distances of all pairs.
• Be clever with NumPy arrays.

But I would appreciate input on which of these is fastest, works best in Python etc. Does anyone have a feel for how good the 'best' performance of Python will be compared to, for example, C++?

import numpy as np
import random


# inters is a list of lists of lists: the first index is point_1_no, the second index is point_2_no and
# the third index is either 0 for the lower bound or 1 for the upper bound.
# See below for sample data and code for reading it into the inters list.

random_size = 20
locations = [np.array([random.random() * random_size for j in range(3)]) for i in range(len(inters))]
while not run_iters(locations, inters):
    locations = [np.array([random.random() * random_size for j in range(3)]) for i in range(len(inters))]


def run_iters(locations, inters, max_runs=600):
    """Runs the iterative procedure of fulfilling constraints until the sum of violations is zero.
    Returns True if a valid structure is generated and False if the maximum number of runs is reached.
    """
    reached = False
    no_points = len(locations)
    for i in range(max_runs):
        if not check_no_viols(locations, inters, no_points):
            run_constraints(locations, inters, no_points)
        else:
            reached = True
            break
    return reached


def run_constraints(locations, inters, no_points):
    """Cycles over each point pair with constraints and, if all constraints present are violated, moves the pair to
    fulfil one of the constraints."""
    for point_1_no in range(no_points):
        for point_2_no in range(point_1_no):
            diff = locations[point_1_no] - locations[point_2_no]
            distance = np.sqrt(diff.dot(diff))
            if check_viol(distance, inters[point_1_no][point_2_no]) == 1:
                move_pair(locations, point_1_no, point_2_no, distance, inters[point_1_no][point_2_no])


def check_no_viols(locations, inters, no_points):
    """Checks if the structure has any point pairs where the constraint is violated.
    Returns True if the sum of violations is zero and False otherwise.
    """
    no_viols = True
    for point_1_no in range(no_points):
        for point_2_no in range(point_1_no):
            diff = locations[point_1_no] - locations[point_2_no]
            distance = np.sqrt(diff.dot(diff))
            if check_viol(distance, inters[point_1_no][point_2_no]) == 1:
                no_viols = False
                break
    return no_viols


def check_viol(distance, constraint):
    """Returns 1 if the distance between two points violates the constraints present, otherwise returns 0."""
    if constraint is None:
        score = 0
    elif constraint[0] <= distance < constraint[1]:
        score = 0
    else:
        score = 1
    return score


def move_pair(locations, point_1_id, point_2_id, distance, constraint):
    """Move two points to fulfil their distance constraint. The points are moved apart or together equal amounts
    to be at a random distance between the lower and upper bounds of the constraint.
    """
    new_dist = constraint[0] + (constraint[1] - constraint[0]) * random.random()
    shift_dist = (new_dist - distance) / 2
    vec_2_to_1 = locations[point_1_id] - locations[point_2_id]
    unit_2_to_1 = vec_2_to_1 / np.sqrt(vec_2_to_1.dot(vec_2_to_1))
    locations[point_1_id] += unit_2_to_1 * shift_dist
    locations[point_2_id] -= unit_2_to_1 * shift_dist

Sample data can be found here for testing:

https://www.dropbox.com/s/9bb0v4lpwwdfyn7/constraint_example.zip?dl=0

And read in using something like:

import re
tot_atoms = 606
inters = [[None for j in range(tot_atoms)] for i in range(tot_atoms)]
in_file = open('int_example.txt')
for line in in_file:
    split = re.split(r'\s+', line.rstrip())
    inters[int(split[0])][int(split[1])] = [float(split[2]), float(split[3])]
    inters[int(split[1])][int(split[0])] = [float(split[2]), float(split[3])]
in_file.close()

Answer 1

\$ \newcommand{\pt}[2]{(#1 \!\times\! #2)} \$Here's some commentary about potential algorithms particularly suited to Python. Note that a general optimal algorithm might not be suitable for Python because Numpy almost requires you to be able to vectorize your calculations.

The aim, basically, would be to make maximal usage of scipy.spatial.distance. This will require moving the data into homogeneous pure-Numpy arrays. I suggest

1. \$\pt n3\$ array of points
2. \$\pt nn\$ array of lower limits
3. \$\pt nn\$ array of upper limits

You should then utilize

from scipy.spatial import distance

distances = distance.squareform(distance.pdist(points))

to find the \$\pt nn\$ matrix of distances.

The problem with the current technique is the quick invalidation of these generated distances. When you call move_pair or anything that emulates it, you can not further involve it in movements until a new call to distance.pdist(points).

This means that \$\mathcal{O}(n^2)\$ work must be done to move at most \$n\$ points once each. The alternative is what you have done: calculate distances and directions independently with each move_pair.

A faster algorithm (in Python) would require vectorized calls to move the points. The sensible thing to me would be a pseudo-physical simulation; each point has two forces it applies to every other point according to some potential. One should fall off exponentially (or at \$1/d\$) whereas the other should rise exponentially. Where the distance, \$d\$, is between the limits the forces should (approximately) cancel.

The equations can be alike to

import numpy
force_out = numpy.exp(distances - upper_limits)
force_in = numpy.reciprocal(lower_limits - distances)

I have seen something like this used effectively for an online visualisation, although cannot remember what. I suggest doing something similar to simulated annealing where the sharpness of the curve starts pretty flat (say, approx. linearithmic) and ends up close to an infinite potential well.

from numpy import linalg, newaxis

# Normalized vectors
differences = points - points[:, newaxis]
differences /= linalg.norm(differences, axis=-1)[..., newaxis]

# Apply "simulated annealing"... sort of
total_force = (force_out + force_in) ** time

# Get displacements
displacements = differences * total_force

Be careful not to ramp up the energy as well; annealing works as the states move slower towards the solution.

When reaching a non-solution, instead of restarting one should consider finding the most displaced fraction of points and permuting their positions with their close neighbours. This will allow near-solutions to complete more quickly. The radius of permuted points should probably increase with the number of failed attempts.

---

Note that these are very hand-wavy thoughts and I haven't tested them for efficiency or even effectivity. It's possible this algorithm would be rather slow at reaching a solution or would get to a close solution and then stop. It, however, seems like a sensible option.

The main problem with fast convergence with the original solution was its convergence; imagine something like

        A                           C
                    B

With the natural state A B C. This could iterate

A-B

            A                       C
                B

A-C

                   A         C
                B

B-C

                   A     C
                    B

A-B

                   A
                         C

                    B

A-C

                  A
                          C

                    B

...and so on.

With a more physical, simultaneous simulation, I'd expect convergence more like

        A→→→                     ←←←C
                    B→

              A→→           ←←C
                    B→

                 A   ↑     C
                     B

                 A    B    C

Namely, prior movements won't cause unwanted collateral damage since no move will be independent of other point's requirements.

Again, just rambling. I've not based this off of any paper nor any trials of my own. It just feels appropriate.