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  ?

EDIT:

Now that this part is optimized (I didn't chnage the code yet but I will in a few hours), I have a larger code with another problem of vectorization, any help would be appreciated!

Thanks

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

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?

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.

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.

EDIT:

Now that this part is optimized (I didn't chnage the code yet but I will in a few hours), I have a larger code with another problem of vectorization, any help would be appreciated!

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?

Thanks

PS, the import GEO is a module of mine to write a SVG file, it is documenter after the main code, 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