Fast loop to create an array of values [closed]

Question

I have a code that creates a 3D array of values from a minimum to a maximum in X and Z with constant Y.

Right now I make it in normal Python, and then I transform it in a np.array. Is there a way to make it directly a NumPy array? How can I translate the code in NumPy to make it faster?

targets = []

X = Y = 0
STEP_X = 0.1
STEP_Y = 0.2
MIN_X = X

for i in range (1000):
    for j in range (1000):

        targets = targets + [[X, Y, Z]]
        X += STEP_X

    X = MIN_X
    Z += STEP_Z

# Make it a np.array
targets = np.array(targets)

EDIT:

Now that this part is optimized (I didn't change the code yet but I will in a few hours), I have more code with another problem of vectorization.

In my code, I used the symbol #___________PART TO BE VECTORIZED to highlight the slow parts.

This code is an already optimized version of something 5 times slower that I wrote before:

• It takes a list of colors in the CIE XYZ colorspace
• calculates the sRGB color values of the XYZ colors
• generates a convex hull of the colors by Delaunay tetrahedralization
• then creates a list of targets, checks if the targets are inside the hull
• and finally generates a dithered SVG file containing a grid of circles (sRGB colors, random dither)

The two parts I could'nt find a way to vectorize with numpy are:

• the loop where I generate the targets

• and the loop where I associate RGB colors with a random dither algorithm at the end

Could you please have a look and orient me towards a better optimization?

PS, the import GEO is a module of mine to write a SVG file, it is documented after the main code, at the bottom of this page.

'''
Created on 12 juin 2014

@author: gary
thanks to a lot of help by Gareth Rees:
http://codereview.stackexchange.com/questions/41024/faster-computation-of-barycentric-coordinates-for-many-points
http://codereview.stackexchange.com/questions/41316/python-numpy-optimize-module-to-handle-big-file
'''
import GEO
import numpy as np
import scipy.spatial

#___________________________________________________________
#__________FUNCTIONS________________________________________

def iround(x):
    """iround(number) -> integer
    Round a number to the nearest integer.
    http://www.daniweb.com/software-development/python/threads/299459/round-to-nearest-integer"""
    return int(round(x) - .5) + (x > 0)


def XYZ2sRGB(X,Y,Z):
    """transforms CIE XYZ tristimulus values
    into sRGB values with gamma = 2.4"""

    X = float(X)
    Y = float(Y)
    Z = float(Z)

    var_X = X / 100        #X from 0 to  95.047      (Observer = 2deg, Illuminant = D65)
    var_Y = Y / 100        #Y from 0 to 100.000
    var_Z = Z / 100        #Z from 0 to 108.883

    var_R = var_X *  3.2406 + var_Y * -1.5372 + var_Z * -0.4986
    var_G = var_X * -0.9689 + var_Y *  1.8758 + var_Z *  0.0415
    var_B = var_X *  0.0557 + var_Y * -0.2040 + var_Z *  1.0570

    if ( var_R > 0.0031308 ): 
        var_R = 1.055 * ( var_R ** ( 1 / 2.4 ) ) - 0.055
    else:                     
        var_R = 12.92 * var_R
    if ( var_G > 0.0031308 ): 
        var_G = 1.055 * ( var_G ** ( 1 / 2.4 ) ) - 0.055
    else:                     
        var_G = 12.92 * var_G
    if ( var_B > 0.0031308 ): 
        var_B = 1.055 * ( var_B ** ( 1 / 2.4 ) ) - 0.055
    else:                     
        var_B = 12.92 * var_B

    R = var_R * 255
    if (R > 255):
        R = 255
    if (R < 0):
        R = 0

    G = var_G * 255
    if (G > 255):
        G = 255
    if (G < 0):
        G = 0

    B = var_B * 255 
    if (B > 255):
        B = 255
    if (B < 0):
        B = 0

    return iround(R), iround(G), iround(B)

#________________________________________________________________
#__________CODE__________________________________________________


# Configuration
POINTS_FILENAME = 'colors.csv'


# Load XYZ tristimulus colors values
XYZ = np.loadtxt(POINTS_FILENAME, usecols=(2,3,4), delimiter=',')
print "XYZ colors loaded"

# Load color names
colornames = np.loadtxt(POINTS_FILENAME, usecols=(1,), delimiter=',',
                        converters={0:lambda s:s.split()}, dtype=np.str)
print "colornames loaded"

# Make sRGB values of XYZ tristimulus colors values
sRGB = []

for i in range(len(XYZ)):

    rgb = XYZ2sRGB(XYZ[i][0],XYZ[i][1],XYZ[i][2])
    sRGB = sRGB + [rgb]

# Make it a np.array
sRGB = np.array(sRGB)
print "sRGB colors computed"

# Encode XYZ color of the support
SUPPORT = np.array([86.83449926, 90.41826972, 101.2739682])


# Average XYZ colors with a weighted amount of support 
# Parameters
SUPPORT_AMOUNT   = .3333333333333
REST             = 1 - SUPPORT_AMOUNT
SUPPORT_WEIGHTED = np.multiply(SUPPORT, SUPPORT_AMOUNT)
XYZ_WEIGHTED     = np.multiply(XYZ, REST)

# Resulting list of points
XYZplusSUPPORT   = np.add(SUPPORT_WEIGHTED, XYZ_WEIGHTED)


# Compute Delaunay tetrahedralization of the new points
tri = scipy.spatial.Delaunay(XYZplusSUPPORT, furthest_site=False) 

# indices of vertices
indices = tri.simplices

# vertices for each tetrahedron
vertices = XYZplusSUPPORT[indices]
print "tetrahedralization OK"

# Make XYZ target values

# Limits of the cube containing XYZ+SUPPORT values
MIN_X, MAX_X   = np.min(XYZplusSUPPORT[:,0]), np.max(XYZplusSUPPORT[:,0])
MIN_Z, MAX_Z   = np.min(XYZplusSUPPORT[:,2]), np.max(XYZplusSUPPORT[:,2])

# custom limits
print "custom limits for X/Z? (Y+ENTER)"

INFO = raw_input()

if INFO == 'Y':
    print 'MIN_X'
    MIN_X = float(raw_input())
    print 'MAX_X'
    MAX_X = float(raw_input())
    print 'MIN_Y'
    MIN_Y = float(raw_input())
    print 'MAX_Y'
    MAX_Y = float(raw_input())


# Target Y
#87.618, 76.303, 66, 56.681, 48.278, 40.749, 34
Y = 34
X, Z = MIN_X, MIN_Z

# Size of the canvas to project targets
SIZE_X, SIZE_Z = 48, 54

# Diameter of points of color, and frequency of the grid
DIAM = .15
FREQ = DIAM + (np.sqrt(np.pi)*(DIAM+((DIAM*np.sqrt(REST)-DIAM*REST)/REST))-2*DIAM)/2 

# Amount of steps on the canvas
STEPS_X = SIZE_X/FREQ
STEPS_Z = SIZE_Z/FREQ

# Range of axis X and axis Z
RANGE_X = MAX_X - MIN_X
RANGE_Z = MAX_Z - MIN_Z

# Size of a step in the colorspace
XYZ_STEP_X = RANGE_X/STEPS_X
XYZ_STEP_Z = RANGE_Z/STEPS_Z

# integer rounded amount of steps
ROUND_X = iround(STEPS_X+1)
ROUND_Z = iround(STEPS_Z+1)

# Targets container
targets = []

# Make targets
# _________PART TO BE VECTORIZED______
for i in range (ROUND_Z+1):
    for j in range (ROUND_X+1):

        targets = targets + [[X, Y, Z]]
        X += XYZ_STEP_X

    X = MIN_X
    Z += XYZ_STEP_Z

# Make it a np.array
targets = np.array(targets)
print "targets OK"

# Find the tetrahedron containing each target (or -1 if not found)
tet = tri.find_simplex(targets)

# Affine transformation for tetrahedron containing each target
U = tri.transform[tet, :3]

# Offset of each target from the origin of its containing tetrahedron
V = targets - tri.transform[tet, 3] 

# Barycentric coordinates of each target in its tetrahedron.
b = np.einsum('ijk,ik->ij', U, V)
bcoords = np.c_[b, 1 - b.sum(axis=1)]
print "bcoords OK"


# Get the sRGB color corresponding to each vertex 
C = sRGB[tri.simplices]



# A uniform random number in [0, 1] for each target.
RAND = np.random.uniform(0, 1, size=(len(targets)))
print "random OK"

# SVG file header
FILENAME = str(Y)+'.svg'
GEO.header(FILENAME, SIZE_X, SIZE_Z)


# Transpose the targets in Centimeters
TARGETS_CM_X = np.subtract(targets[:,0], MIN_X)
TARGETS_CM_X = np.divide(TARGETS_CM_X, RANGE_X)
TARGETS_CM_X = np.multiply(TARGETS_CM_X, SIZE_X)

TARGETS_CM_Z = np.subtract(targets[:,2], MIN_Z)
TARGETS_CM_Z = np.divide(TARGETS_CM_Z, RANGE_Z)
TARGETS_CM_Z = np.multiply(TARGETS_CM_Z, SIZE_Z)

print "target transposed to cm"

#_________________________________________
# PART TO BE VECTORIZED

for i in range(len(tet)):

    if(tet[i] != -1):

        R = RAND[i]

        x = TARGETS_CM_X[i]
        z = TARGETS_CM_Z[i]

        if R <= bcoords[i][0]:
            R,G,B = C[tet][i][0][0], C[tet][i][0][1], C[tet][i][0][2]

        elif R <= bcoords[i][0]+bcoords[i][1]:
            R,G,B = C[tet][i][1][0], C[tet][i][1][1], C[tet][i][1][2]

        elif R <= bcoords[i][0]+bcoords[i][1]+bcoords[i][2]:
            R,G,B = C[tet][i][2][0], C[tet][i][2][1], C[tet][i][2][2]

        else:
            R,G,B = C[tet][i][3][0], C[tet][i][3][1], C[tet][i][3][2]


        GEO.DISC(FILENAME, x, z, DIAM/2, R, G, B)

GEO.END(FILENAME)

print "file written"

The GEO module:

def header(filename,Xmax,Ymax):
""" header of a SVG file"""

        Xmax = float(Xmax)
        Ymax = float(Ymax)

        # SVG Header 
        f = open(str(filename), "w")
        f.write('<svg version="1.1"'+'\n')
        f.write('    baseProfile="full"'+'\n')
        f.write('    width="'+ str(Xmax*1/2.54*72) + '" '+'height="'+ str(Ymax*1/2.54*72) +'"'+'\n')
        f.write('    xmlns="http://www.w3.org/2000/svg">'+'\n')
        f.close()
        print "header written, filename is:", filename
        return None

def DISC(filename,x,y,radius, R, G, B):
            '''
            circle path in a SVG file
            x,y = center of the disk
            http://stackoverflow.com/questions/5737975/circle-drawing-with-svgs-arc-path
            '''

            radius = radius*1/2.54*72
            f = open(str(filename), "a")
            f.write('<path d="M'+str(x*1/2.54*72)+" "+str(y*1/2.54*72)+" \n") #moveto
            f.write('    m '+str(-radius)+ ',0 \n')
            f.write('    a '+str(radius)+','+str(radius)+ ' 0 1,0 ' + str(radius*2)+',0 \n')
            f.write('    a '+str(radius)+','+str(radius)+ ' 0 1,0 ' + str(-radius*2)+',0 \n')      
            f.write('    " fill = "rgb('+str(R)+','+str(G)+','+str(B)+')"/> \n')
            f.close()
            return None

def END(filename):
""" closing the svg file"""

            f = open(str(filename), "a")
            f.write('</svg>')
            f.close()
            return None

Answer 1

You could use the np.fromiter function and Python's built in itertools.product to create the array you need:

Note: I'm assuming you're using Python 2.x based on your print statements.

import itertools
import numpy as np

product = itertools.product(xrange(X, X + 1000*STEP_X, STEP_X),
                            [Y],
                            xrange(Z, Z + 1000*STEP_Z, STEP_Z))

targets = np.fromiter(product)

This should be faster because it uses iterators instead of creating and allocating an entire list.

---

UPDATE

Here are some style pointers and other minor improvements that I could see. Most of these recommendations stem from PEP8, the official Python style guide, so if you need a reference for my suggestions, you can head over there.

1. ALWAYS USE with. Whenever you deal with file access, use a with block as it is significantly less prone to user errors than using open() and close(). Luckily, you're code doesn't show the typical bug of not calling close() after an open(). However, its best to get into the habit of using with:

   with open('some_file.txt', 'r') as file: # Do stuff

2. Use underscores_in_names when naming variables and functions. For the most part you do this. However, your function names could be updated.

3. Function names should be verb-based as this style helps show that the function does something:

   # Currently...
   def XYZ2sRGB(...):
   
   # Better...
   def convert_to_RGB(...)
   

   A quick note: Typically I don't like using upper-case letters in anything except constants. However, because RGB is basically an acronym, capital letters seem appropriate.

4. Speaking about upper-case letters, convention says that only constants should be capitalized in Python. This is relatively significant because convention is the only way we can 'define' constants in Python as there is no syntactic way to do so.

5. Whitespace is your friend, however be careful not to overdo it. PEP8 actually calls extraneous whitespace a pet peeve. A few of the points mentioned in that section of PEP8 that are applicable are:

   # Bad                  # Good
   foo            = 1  |  foo = 1
   some_long_name = 0  |  some_long_name = 0
   --------------------+---------------------
   range (1000)        |  range(1000)
   --------------------+---------------------
   foo = ( x + 1 * 2 ) |  foo = (x + 1*2)
   

   The last example is really based on preference: simply use whitespace to group operations and operands together so that the calculation reads well.

6. Parenetheses aren't required in if statements (unless they group conditionals together). You can remove almost all of yours.

7. Use if ... elif ... when applicable. Take this group of statements:

   G = var_G * 255
   if (G > 255):
       G = 255
   if (G < 0):
       G = 0
   

   The second if will always be evaluated even if the first evaluated to True which means the second will evaluate to False. Because the two conditional are mutually exclusive, use and if-elif structure. Also, instead of basing your conditionals off of G (which requires a calculation beforehand) base your conditionals off of var_G:

   if var_G > 1:
       G = 255
   elif var_G < 0:
       G = 0
   else:
       G = var_G * 255
   

   This code only does the calculation if necessary and has the same number of possible comparisions (in the worst case).

8. Use str.format instead of string concatenation. While whether string formatting performs better than string concatenation is up in the air, its more conventional (and, in my opinion, MUCH cleaner) to use str.format:

   with open(str(filename), "a") as f:
       f.write('<path d="M{} {} \n'.format(x*1/2.54*72, y*1/2.54*72)) #moveto
       f.write('    m {},0 \n'.format(-radius))
       f.write('    a {} 0 1,0 {},0 \n'.format(radius, radius, radius*2))
       f.write('    a {} 0 1,0 {},0 \n'.format(radius, radius, -radius*2))
       f.write('    " fill = "rgb({},{},{})"/> \n'.format(R, G, B))