I'm just starting to understand the Python syntax and I created a module that does what I wanted, but really slow.Here are the stats of cProfile, top 10 ordered by internal time:   ncalls  tottime  percall  cumtime  percall filename:lineno(function)  1291716   45.576    0.000  171.672    0.000 geometry.py:10(barycentric_coords)  6460649   31.617    0.000   31.617    0.000 {numpy.core.multiarray.array}  2583432   15.979    0.000   15.979    0.000 {method 'reduce' of 'numpy.ufunc'objects}     2031   12.032    0.006  193.333    0.095 geometry.py:26(containing_tet)  1291716   10.944    0.000   58.323    0.000 linalg.py:244(solve)  1291716    7.075    0.000    7.075    0.000 {numpy.linalg.lapack_lite.dgesv}  1291716    5.750    0.000    9.865    0.000 linalg.py:99(_commonType)  2583432    5.659    0.000    5.659    0.000 {numpy.core.multiarray._fastCopyAndTranspose}  1291716    5.526    0.000    7.299    0.000 twodim_base.py:169(eye)  1291716    5.492    0.000   12.791    0.000 numeric.py:1884(identity)To process data and create a 24 bitmap of 300*300 pixels it would take about 1 1/2 hour on my laptop with the latest intel i5, a SSD disk and 12gb RAM!I am not sure if it's a good idea to post all the code, but if not how can you get the entire picture?The first thing colors_PPM.py does is to access a function called geometry.geometry.py itself calls tetgen, tetgen.py calls tetgen.exe.The module colors_PPM.py reads a list in targets.csv and a list in XYZcolorlist_D65.csv, it excludes some elements of XYZcolorlist_D65.csv and then reads one by one the rows of targets.csv, performs a delaunay triangulation via tetgen and returns 4 names[] and 4 bcoords[].Then random is used to choose one name by a series of if elif else tests.Finally the result is exported in a bitmap file epson_gamut.pbm.Do you see any way this could run faster? I know it should seem quite a mess..csv files structures examples:'XYZcolorlist_D65.csv'255 63 127,35.5344302104,21.380721966,20.3661095969255 95 127,40.2074945517,26.5282949405,22.7094284437255 127 127,43.6647438365,32.3482625492,23.6181801523255 159 127,47.1225628354,39.1780944388,22.9366615044255 223 159,61.7379149646,62.8387601708,32.3936200864255 255 159,70.7428790853,78.6134546144,29.5579371353255 0 127,32.0951763469,18.3503537537,19.0863164396255 31 127,32.281389139,18.5592317077,18.6802029444255 191 127,52.6108977261,48.5621713952,21.7645428218255 223 127,59.7600830083,60.9770436618,20.9338174593...'targets.csv'30,5,330.34,5,330.68,5,331.02,5,331.36,5,331.7,5,332.04,5,332.38,5,332.72,5,333.06,5,333.4,5,333.74,5,3...This is geometry.py:#geometry.pyimport numpy as npimport numpy.linalg as laimport tetgendef barycentric_coords(vertices, point):    T = (np.array(vertices[:-1])-vertices[-1]).T    v = np.dot(la.inv(T), np.array(point)-vertices[-1])    v.resize(len(vertices))    v[-1] = 1-v.sum()    return vdef tetgen_of_hull(points):    tg_all = tetgen.TetGen(points)    hull_i = set().union(*tg_all.hull)    hull_points = [points[i] for i in hull_i]    tg_hull = tetgen.TetGen(hull_points)    return tg_hull, hull_idef containing_tet(tg, point):    for i, tet in enumerate(tg.tets):        verts = [tg.points[j] for j in tet]        bcoords = barycentric_coords(verts, point)        if (bcoords >= 0).all():            return i, bcoords    return None, NoneThis is tetgen.py:#tetgenimport tempfile, subprocess, osclass TetGen:    def __init__(self, points):        self.points = points        node_f = tempfile.NamedTemporaryFile(suffix=.node, delete=False);        node_f.write(%i 3 0 0\n % len(points))        for i, point in enumerate(points):            node_f.write(%i %f %f %f\n % (i, point[0], point[1], point[2]))        node_f.close()        subprocess.call([tetgen, node_f.name], stdout=open(os.devnull, 'wb'))        #subprocess.call([C:\Users\gary\Documents\eclipse\Light Transformer Setup\tetgen, node_f.name], stdout=open(os.devnull, 'wb'))        ele_f_name = node_f.name[:-5] + .1.ele        face_f_name = node_f.name[:-5] + .1.face        ele_f_lines = [line.split() for line in open(ele_f_name)][1:-1]        face_f_lines = [line.split() for line in open(face_f_name)][1:-1]        self.tets = [map(int, line[1:]) for line in ele_f_lines]        self.hull = [map(int, line[1:]) for line in face_f_lines]if __name__ == '__main__':    from pprint import pprint    points = [(0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)]    tg = TetGen(points)    pprint(tg.tets)    pprint(tg.hull)And this is the module colors_PPM.py:import geometryimport csvimport numpy as npimport randomimport cv2S = 0img = cv2.imread(MAP.tif, -1)height, width = img.shapeppm = file(epson gamut.ppm, 'w')ppm.write(P3 + \n + str(width) +   + str(height) + \n + 255 + \n)# PPM file headerall_colors = [(name, float(X), float(Y), float(Z))              for name, X, Y, Z in csv.reader(open('XYZcolorlist_D65.csv'))]tg, hull_i = geometry.tetgen_of_hull([(X,Y,Z) for name, X, Y, Z in all_colors])colors = [all_colors[i] for i in hull_i]print (thrown out:        + , .join(set(zip(*all_colors)[0]).difference(zip(*colors)[0])))targets = ((float(X), float(Y), float(Z))           for X, Y, Z in csv.reader(open('targets.csv')))for target in targets:    X, Y, Z = target    target_point = (np.array([X,Y,Z]))    tet_i, bcoords = geometry.containing_tet(tg, target_point)    if tet_i == None:           ppm.write(str(255 255 255) + \n)        continue         # not in gamut    else:        A = bcoords[0]        B = bcoords[1]        C = bcoords[2]        D = bcoords[3]        R = random.uniform(0,1)        names = [colors[i][0] for i in tg.tets[tet_i]]        if R <= A:            S = names[0]         elif R <= A+B:            S = names[1]        elif R <= A+B+C:            S = names[2]        else:            S = names[3]        ppm.write(str(S) + \n)print doneppm.close()
Faster computation of barycentric coordinates for many points
python;performance;numpy;python 2.7;computational geometry
1. IntroductionThanks for running the profiler. As you can see from the output, most of the runtime is being spent in your containing_tet function.The first thing to say is that you have made this question unnecessarily difficult for us because your functions have no documentation. We have to read and reverse-engineer your code to try to figure out what the purpose of each function is. Python allows you to write a docstring for each function in which you explain what the function does, what arguments it takes, and what values it returns.It looks to me as though containing_tet tries to find a tetrahedron containing point, and it does so by computing the barycentric coordinates of point within each tetrahedron and checking that all the coordinates are non-negative.So you could have written:def containing_tet(tg, point):    If point is inside the i'th tetrahedron in tg, return the pair    (i, bcoords) where bcoords are the barycentric coordinates of    point within that tetrahedron. Otherwise, return (None, None).    and so on for your other functions. Writing this kind of documentationwill help your colleagues (and you in a few months when you haveforgotten all the details).2. DiagnosisWhy is containing_tet slow? Well, you say that you are calling itmany times, and each time it has to repeat some work. First, itcollects the vertices of each tetrahedron:verts = [tg.points[j] for j in tet]and then in barycentric_coords it computes the transform matrix:T = (np.array(vertices[:-1])-vertices[-1]).T... la.inv(T) ...but these computations will be the same every time (they do not depend on point). It would be best to compute them just once.Then you need to reorganize your code so that each operation isvectorized: that is, you should read all the targets into a NumPyarray, and then compute the barycentric coordinates for all thetargets at once.3. Using scipy.spatialHaving written the above, however, I am wondering why you did not use scipy.spatial.Delaunay? Is it because TetGen does a better job?Since I don't have a copy of TetGen handy, if I had to write this code I would generate the triangulations like this:>>> import numpy as np>>> import scipy.spatial>>> points = np.array([(0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1)])>>> tri = scipy.spatial.Delaunay(points)>>> tri.simplicesarray([[3, 2, 4, 0],       [3, 1, 4, 0],       [3, 6, 2, 4],       [3, 6, 7, 4],       [3, 5, 1, 4],       [3, 5, 7, 4]], dtype=int32)and then if I have an array of targets to query:>>> targets = np.array([[.1,.1,.1], [.9,.9,.9], [.1,.6,.7], [.4,.9,.1]])I can find which tetrahedron each point belongs to by calling scipy.spatial.Delaunay.find_simplex:>>> tetrahedra = tri.find_simplex(targets)>>> tetrahedraarray([0, 3, 1, 2], dtype=int32)And then I can find the barycentric coordinates of each point within its tetrahedron like this:>>> X = tri.transform[tetrahedra,:3]>>> Y = targets - tri.transform[tetrahedra,3]>>> b = np.einsum('ijk,ik->ij', X, Y)>>> bcoords = np.c_[b, 1 - b.sum(axis=1)]>>> bcoordsarray([[ 0.1,  0. ,  0.1,  0.8],       [ 0.1,  0. ,  0.8,  0.1],       [ 0.6,  0.1,  0.1,  0.2],       [ 0.1,  0.3,  0.5,  0.1]])This is essentially following the recipe in the scipy.spatial.Delaunay documentation, except that I transform each point using the affine transformation for the tetrahedron it was found in. Note that in the final result, all the barycentric coordinates are in the range 01 as you would expect.The bit of the computation that is tricky to figure out how to vectorize is the computation of b. numpy.dot computes the dot product of two arrays, but here I need an array of dot products. I could loop over the elements of X and Y in Python, like this:>>> b = np.array([x.dot(y) for x, y in zip(X, Y)])but that would be using slow Python iteration rather than fast NumPy vector operations. Hence the rather hairy use of numpy.einsum.Note that you'll have to do something about the points that were not found in any tetrahedron. scipy.spatial.Delaunay.find_simplex returns -1 for these points. You could mask out the points you need, as described here, or you could try using a masked array. (Or try both and see which is faster.)4. Answers to questionsIn comments you asked:What is vectorization? See the What is NumPy? section of the NumPy documentation, in particular the section starting:Vectorization describes the absence of any explicit looping, indexing, etc., in the code - these things are taking place, of course, just behind the scenes (in optimized, pre-compiled C code).Taking advantage of NumPy vector operations is the whole point for using NumPy. For example, this program computes 9 million multiplications and additions in Python:>>> v = np.arange(3000)>>> sum(i * j for i in v for j in v)20236502250000It takes about ten seconds on this computer. But the corresponding operation in NumPy:>>> np.sum(np.outer(v, v))20236502250000is about 100 times faster, because the loops involved in this computation run in fast machine code and not in slow Python code. So when you are working on a program in NumPy, you need to scrutinize every Python loop  that is, every for statement  to see if you can turn it into a single NumPy operation.How to structure your code? See 5 below. Notice how using NumPy efficiently requires us to turn the code inside out. Instead of looping over the targets at top level, and in each loop iteration performing a series of operations on that one target, we lift the series of operations up to top level and make each operation run over all elements of a NumPy array.How to read the colors from the colorlist CSV? In NumPy each array must have a single data type (that's one the requirements for NumPy to be able to process them quickly). So you must read the names into one array and the points into another array. An easy way to do this is to call numpy.loadtxt twice:>>> colors = np.loadtxt('XYZcolorlist_D65.csv', usecols=(0,), delimiter=',', ...                     converters={0:lambda s:s.split()}, dtype=np.uint8)>>> colors[:5]array([[255,  63, 127],       [255,  95, 127],       [255, 127, 127],       [255, 159, 127],       [255, 223, 159]], dtype=uint8)>>> points = np.loadtxt('XYZcolorlist_D65.csv', usecols=(1,2,3), delimiter=',')>>> points[:5]array([[ 35.53443021,  21.38072197,  20.3661096 ],       [ 40.20749455,  26.52829494,  22.70942844],       [ 43.66474384,  32.34826255,  23.61818015],       [ 47.12256284,  39.17809444,  22.9366615 ],       [ 61.73791496,  62.83876017,  32.39362009]])5. Revised code.import numpy as npimport scipy.spatial# Configuration.POINTS_FILENAME = 'XYZcolorlist_D65.csv'TARGETS_FILENAME = 'targets.csv'OUTPUT_FILENAME = 'gamut.ppm'DEFAULT_COLOR = np.array([[255, 255, 255]], dtype=np.uint8)# Load colorscolors = np.loadtxt(POINTS_FILENAME, usecols=(0,), delimiter=',',                     converters={0:lambda s:s.split()}, dtype=np.uint8)# Load pointspoints = np.loadtxt(POINTS_FILENAME, usecols=(1, 2, 3), delimiter=',')# Load targetstargets = np.loadtxt(TARGETS_FILENAME, delimiter=',')ntargets = len(targets)# Compute Delaunay triangulation of points.tri = scipy.spatial.Delaunay(points)# Find the tetrahedron containing each target (or -1 if not found)tetrahedra = tri.find_simplex(targets)# Affine transformation for tetrahedron containing each targetX = tri.transform[tetrahedra, :3]# Offset of each target from the origin of its containing tetrahedronY = targets - tri.transform[tetrahedra, 3]# First three barycentric coordinates of each target in its tetrahedron.# The fourth coordinate would be 1 - b.sum(axis=1), but we don't need it.b = np.einsum('...jk,...k->...j', X, Y)# Cumulative sum of barycentric coordinates of each target.bsum = np.c_[b.cumsum(axis=1), np.ones(ntargets)]# A uniform random number in [0, 1] for each target.R = np.random.uniform(0, 1, size=(ntargets, 1))# Randomly choose one of the tetrahedron vertices for each target,# weighted according to its barycentric coordinates, and get its# color.C = colors[tri.simplices[tetrahedra, np.argmax(R <= bsum, axis=1)]]# Mask out the targets where we failed to find a tetrahedron.C[tetrahedra == -1] = DEFAULT_COLOR# Determine width and height of image.# (Since I don't have your TIFF, this is the best I can do!)width, height = 1, ntargets# Write output as image in PPM format.with open(OUTPUT_FILENAME, 'wb') as ppm:    ppm.write(P3\n{} {}\n255\n.format(width, height).encode('ascii'))    np.savetxt(ppm, C, fmt='%d')