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