Calculate angle between planes

Question

I have written working code calculating the angle between the adjacent planes. I read subsequently from standart input: n - amount of triangles, m - amount of vertices ind[] - indices of the vertices given coord[] - coordinates of the vertices given

The output is supposed to be the maximum angle between the adjacent planes in rad.

Function calculate_angle() iterates over the amount of triangles so that I compare only new ones with the old ones (for z in range(k, n))

list_result = [i for i in index1 if i in index2] used for the adjacency detection: it asks, whether the indices of the first triangle coordinates == to the second. If there are at leat 2 of them - we start to calculate the normals to triangles (the angle between surfaces = to the angle between their normals)

However, if the number of triangles is more than 10⁴ it starts to work very slowly:

import numpy as np
import math
import sys
import cProfile, pstats, io
pr = cProfile.Profile()
pr.enable()

num_triangles, num_vertices = input().split()
n = int(num_triangles)  # amount of triangles
m = int(num_vertices)  # amount of vertices
ind = []
coord = []
angles_list =[]

for i in range(n):
    ind.append([int(j) for j in input().split()])  # indices of the vertices given
for j in range(m):
    coord.append([float(k) for k in input().split()])  # coordinates of the vertices given


def unit_vector(v):
    # Returns the unit vector of the vector.
    return v/ np.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2])


def angle_between(v1, v2):
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return math.acos(max(min(np.dot(v1_u, v2_u), 1), -1))   # (np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))



def calculate_angle():
    for k in range(0, n):
        for z in range(k, n):
            index1 = ind[k]
            index2 = ind[z]
            trignum =0
            list_result = [i for i in index1 if i in index2]
            if (ind[k] != ind[z])&(len(list_result) >= 2)&(trignum <= 3):
                trignum = trignum + 1
                n1 = angle_norm(index1)
                n2 = angle_norm(index2)
                ang = angle_between(n1, n2)
                angles_list.append(ang)
    return max(angles_list)


def cross(a, b):
    c = [a[1]*b[2] - a[2]*b[1],
         a[2]*b[0] - a[0]*b[2],
         a[0]*b[1] - a[1]*b[0]]
    return c

def angle_norm(triangle_index):
    p0 = coord[triangle_index[0]]
    p1 = coord[triangle_index[1]]
    p2 = coord[triangle_index[2]]
    V1 = np.array(p1) - np.array(p0)
    V2 = np.array(p2) - np.array(p0)
    an = cross(V1, V2)
    return an



print(calculate_angle())

pr.disable()
s = io.StringIO()
sortby = 'tottime'
ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
ps.print_stats()
print(s.getvalue())

The analysis results:

0.0
         12130580 function calls in 26.836 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  1048576    6.038    0.000    6.038    0.000 /Users/mac/PycharmProjects/untitled3/Course:21(unit_vector)
     1284    5.309    0.004    5.309    0.004 {built-in method builtins.input}
  4194304    4.645    0.000    4.645    0.000 {built-in method numpy.core.multiarray.array}
  1048576    3.304    0.000   10.566    0.000 /Users/mac/PycharmProjects/untitled3/Course:55(angle_norm)
  1048576    2.617    0.000    2.617    0.000 /Users/mac/PycharmProjects/untitled3/Course:49(cross)
        1    2.090    2.090   21.516   21.516 /Users/mac/PycharmProjects/untitled3/Course:33(calculate_angle)
   524288    0.996    0.000    8.224    0.000 /Users/mac/PycharmProjects/untitled3/Course:26(angle_between)
   524288    0.531    0.000    0.531    0.000 {built-in method numpy.core.multiarray.dot}
   819840    0.487    0.000    0.487    0.000 /Users/mac/PycharmProjects/untitled3/Course:39(<listcomp>)
   524288    0.329    0.000    0.329    0.000 {built-in method builtins.min}
   524289    0.249    0.000    0.249    0.000 {built-in method builtins.max}
   524288    0.095    0.000    0.095    0.000 {built-in method math.acos}
   819840    0.080    0.000    0.080    0.000 {built-in method builtins.len}
   525571    0.055    0.000    0.055    0.000 {method 'append' of 'list' objects}
        1    0.008    0.008    0.008    0.008 {built-in method builtins.print}
     1280    0.002    0.000    0.002    0.000 /Users/mac/PycharmProjects/untitled3/Course:16(<listcomp>)
     1284    0.000    0.000    0.000    0.000 {method 'split' of 'str' objects}
        3    0.000    0.000    0.000    0.000 /Users/mac/PycharmProjects/untitled3/Course:18(<listcomp>)
        1    0.000    0.000    0.000    0.000 /Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/codecs.py:318(decode)
        1    0.000    0.000    0.000    0.000 {built-in method _codecs.utf_8_decode}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
Here's what I already tried to optimise: 

1. I got rid of couple of np built-in functions, e.g. np.cross() and np.linalg.norm(), that gave me a couple of seconds.
2. It was for z in range(1, n), I changed 1 to k in order to not take into account already calculated triangles.

I also tried to make faster input, but to no avail. (Used map, used sys.std)

Please, can someone tell me how to make it significantly faster? I'm not well-acquainted with graphs, and I have a bad feeling about this...

Answer 1

One immediate advice is to precompute angle_norm for each triangle. That alone will give you some boost.

It also lets reformulate the problem as finding a point set diameter. It is a classical problem, and a generic case is quite hard. You may want to look here for an introduction and references.

However, in your case you deal with normalized vectors; their endpoints lay on the unitary sphere, and the problem reduces to a simpler 2D case, which admits the \$O(n\log{n})\$ solution.

PS

About precomputing normals:

        normals = []
        def compute_normals(triangles):
            for t in triangles:
                normals.append(unit_vector(angle_norm(t))

        def calculate_angle():
            ....
            if (....):
                ....
                # Already computed!
                # n1 = angle_norm(index1)
                # n2 = angle_norm(index2)
                ang = angle_between(normals[index1], normals[index2])
                angles_list.append(ang)