6
\$\begingroup\$

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 104 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...

\$\endgroup\$
  • \$\begingroup\$ "I also tried to make faster input, but to no avail." What do you mean in particular, could you extend a bit about that. Also did you already use a decent profiling tool to discover the bottlenecks in your code? \$\endgroup\$ – πάντα ῥεῖ May 5 '17 at 17:59
  • 1
    \$\begingroup\$ @πάνταῥεῖ The code clearly uses cProfile. \$\endgroup\$ – 200_success May 5 '17 at 18:04
  • \$\begingroup\$ My bad. I should have pointed that out before directing him here. \$\endgroup\$ – Carcigenicate May 5 '17 at 18:04
  • 1
    \$\begingroup\$ It would be nice if you included a small sample of the input in the question. \$\endgroup\$ – 200_success May 5 '17 at 18:07
  • 1
    \$\begingroup\$ @πάνταῥεῖ, I added the results into the question. The bottleneck is the calculate_angle() function. I have no better idea than to think i's because of the 2 for-cycles. I think the only faster algorithm would be some kind of a graph -thing, but I'm not an expert, as I pointed out, so I need some help. \$\endgroup\$ – TheDoctor May 5 '17 at 18:22
4
\$\begingroup\$

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)
\$\endgroup\$
  • 1
    \$\begingroup\$ sounds very interesting, but could you please elaborate about precomputing? I cannot understand how to achieve that. Thank you for the article. \$\endgroup\$ – TheDoctor May 5 '17 at 22:43
  • 1
    \$\begingroup\$ @TheDoctor See edit. \$\endgroup\$ – vnp May 5 '17 at 22:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.