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 <= distance < constraint: 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 + (constraint - constraint) * 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:
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)][int(split)] = [float(split), float(split)] inters[int(split)][int(split)] = [float(split), float(split)] in_file.close()