Iterate over coordinates and correct constraint violations

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()

• I've been thinking about this, but there's no obvious way to give large speed improvements whilst retaining the sequential nature of the computation (even if emulated) without some really clever pairing of vertices. There are a lot of algorithms that could get the desired result, but all will be distinct from the given method. Jan 19, 2015 at 15:10
• As such, rewriting this with Numba or jitpy seems the best bet. jitpy would be easier but Numba might be faster. Jan 19, 2015 at 15:12
• Thanks @Veedrac, I will give those a look. Out of interest what other algorithms would you recommend to fulfil the constraints using another procedure? I was limited to this procedure as I was attempting to recreate and build on someone's data, but knowing other approaches would be useful in the long run. Jan 19, 2015 at 16:03
• I have now asked the more general question on the Computer Science Stack Exchange: cs.stackexchange.com/questions/37389/… Jan 19, 2015 at 17:26

$\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.

• This seems very similar to the "spring" model that D.W. mentioned in the CS crosspost. Lots of other interesting ideas there, too, and the Wikipedia link he gave is pretty neat too. Jan 19, 2015 at 21:22