Learning SciPy for Numerical and Scientific Computing

Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 Sergio Rojas (srojas@usb.ve) and Erik A Christensen (erikcny@aol.com).

NOTE: This IPython notebook should be read alonside the corresponding chapter in the book, where each piece of code is fully explained.

Chapter 2. Working with Numpy Array as a first step to Scipy

Summary

In this chapter we explore in depth the creation and basic manipulation of the object array used by SciPy, as an overview of the NumPy libraries. In particular, we'll introduce the principles of slicing and masking, which simplify the coding of algorithms to the point of transforming an otherwise unreadable sequence of loops and primitive commands, into an intuitive and self-explanatory set of object calls and methods. We also will learn that the nonbasic modules in NumPy are replicated as modules in SciPy itself. The chapter roughly follows the same structure as the official NumPy reference (which the reader can access at the SciPy pages at docs.scipy.org/doc/numpy/reference ). There are other good sources that cover NumPy with rigor, and we refer you to any of that other material for a more detailed coverage of this topic.

References

Scipy Reference
http://docs.scipy.org/doc/scipy/reference/
Image manipulation and processing using Numpy and Scipy
http://scipy-lectures.github.io/advanced/image_processing/
Numpy and Scipy Documentation
http://docs.scipy.org/doc/
Numpy Tutorial
http://wiki.scipy.org/Tentative_NumPy_Tutorial
Numpy Data type objects
http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html

Working with Numpy Array as a first step to Scipy



In [1]:

    
%matplotlib inline
import numpy
import matplotlib.pyplot as plt
import scipy.misc



In [2]:

    
img=scipy.misc.lena()



In [3]:

    
plt.rcParams['figure.figsize'] = (12.0, 8.0)
plt.gray()
plt.imshow(img)
plt.show()



In [4]:

    
print(img[0:3,0:7])









    



[[162 162 162 161 162 157 163]
 [162 162 162 161 162 157 163]
 [162 162 162 161 162 157 163]]



In [5]:

    
img[1,1:6]=0
print(img[0:3,0:7])









    



[[162 162 162 161 162 157 163]
 [162   0   0   0   0   0 163]
 [162 162 162 161 162 157 163]]

Object essentials



In [6]:

    
print(img.dtype, img.shape, img.size)









    



(dtype('int64'), (512, 512), 262144)



In [7]:

    
print(img[32,67])

Datatype



In [8]:

    
img=scipy.misc.lena().astype('float32')
scores = numpy.array([101,103,84], dtype='float32')
print(scores)









    



[ 101.  103.   84.]



In [9]:

    
print(scores.dtype)









    



float32



In [10]:

    
scores = numpy.float32([101,103,84])
print(scores)









    



[ 101.  103.   84.]



In [11]:

    
print(scores.dtype)









    



float32



In [12]:

    
a=numpy.array(['Cleese', 'Idle', 'Gilliam'], dtype='str_')
print(a)









    



['Cleese' 'Idle' 'Gilliam']



In [13]:

    
print(a.dtype)

|S7



In [14]:

    
dt=numpy.dtype([ ('name', numpy.str_, 16), ('grades', numpy.float64, (2,)) ])



In [15]:

    
MA141=numpy.array([ ('Cleese', (7.0,8.0)), ('Gilliam', (9.0,10.0)) ], dtype=dt)
print(MA141)









    



[('Cleese', [7.0, 8.0]) ('Gilliam', [9.0, 10.0])]



In [16]:

    
MA141









    Out[16]:





array([('Cleese', [7.0, 8.0]), ('Gilliam', [9.0, 10.0])], 
      dtype=[('name', 'S16'), ('grades', '<f8', (2,))])



In [17]:

    
print(MA141.dtype)









    



[('name', 'S16'), ('grades', '<f8', (2,))]

Indexing



In [18]:

    
A=numpy.array([[1,2,3,4,5,6,7,8],[2,4,6,8,10,12,14,16]])
print(A[0:2, 0:8:2])









    



[[ 1  3  5  7]
 [ 2  6 10 14]]



In [19]:

    
print(A[0:2, 8:0:-2])









    



[[ 8  6  4  2]
 [16 12  8  4]]



In [20]:

    
print(A[ [0,0,1,1], [0,3,2,5] ])









    



[ 1  4  6 12]



In [21]:

    
print(A[ numpy.ix_( [0,1], [0,3] )])









    



[[1 4]
 [2 8]]

The array object



In [22]:

    
value=0; img=scipy.misc.lena()
for item in img.flat: value+=item
print(value)

32518120

The next instruction will create the file "lena.txt" in the current directory (refers to the textbook for additional details)



In [23]:

    
img.tofile("lena.txt",sep=" ",format="%i")



In [24]:

    
A=numpy.array([11,13,15,17,19,18,16,14,12,10])
print (type(A))









    



<type 'numpy.ndarray'>



In [25]:

    
print(A.argsort(kind='mergesort'))









    



[9 0 8 1 7 2 6 3 5 4]



In [26]:

    
A.argsort(kind='mergesort')









    Out[26]:





array([9, 0, 8, 1, 7, 2, 6, 3, 5, 4])



In [27]:

    
A.sort()
print (A)









    



[10 11 12 13 14 15 16 17 18 19]



In [28]:

    
A=numpy.array([[1,1,1],[2,2,2],[3,3,3]])
print(A.mean())

2.0



In [29]:

    
print(A.mean(axis=0))









    



[ 2.  2.  2.]



In [30]:

    
print(A.mean(axis=1))









    



[ 1.  2.  3.]



In [31]:

    
print(img.min(), img.max(), img.ptp())









    



(25, 245, 220)



In [32]:

    
A=img.clip(img.min(),img.min()+100)
print(A.min(), A.max(), A.ptp())









    



(25, 125, 100)

Routines for array creation



In [33]:

    
Z=numpy.zeros((5,5), dtype=int)
print(Z)









    



[[0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [0 0 0 0 0]]



In [34]:

    
U=numpy.ones((2,2), dtype=int)
print(U)









    



[[1 1]
 [1 1]]



In [35]:

    
I=numpy.identity(3, dtype=int)
print(I)









    



[[1 0 0]
 [0 1 0]
 [0 0 1]]



In [36]:

    
D=numpy.eye(4,k=1) + numpy.eye(4,k=-1)
print( D)









    



[[ 0.  1.  0.  0.]
 [ 1.  0.  1.  0.]
 [ 0.  1.  0.  1.]
 [ 0.  0.  1.  0.]]



In [37]:

    
A=3.0*(numpy.eye(5,k=4) + numpy.eye(5,k=-4))
A[0:2,0:2]=5*U; A[2:5,2:5]=6*I
print(A)









    



[[ 5.  5.  0.  0.  3.]
 [ 5.  5.  0.  0.  0.]
 [ 0.  0.  6.  0.  0.]
 [ 0.  0.  0.  6.  0.]
 [ 3.  0.  0.  0.  6.]]



In [38]:

    
B=numpy.fromfunction( (lambda i,j: i*j), (4,4), dtype=int)
print(B)



In [39]:

    
print(B==0)









    



[[ True  True  True  True]
 [ True False False False]
 [ True False False False]
 [ True False False False]]



In [40]:

    
B += numpy.fromfunction((lambda i,j:i*i+j*j), (4,4))*(B==0)
print(B)









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-40-74d67701c2f7> in <module>()
----> 1 B += numpy.fromfunction((lambda i,j:i*i+j*j), (4,4))*(B==0)
      2 print(B)

TypeError: Cannot cast ufunc add output from dtype('float64') to dtype('int64') with casting rule 'same_kind'



In [ ]:

    
print( B%2!=0)



In [ ]:

    
print(numpy.where(B%2!=0))



In [ ]:

    
numpy.putmask( B, B%2!=0, B**2+1)
print(B)



In [ ]:

    
L1=numpy.arange(-1,1,0.3)
print(L1)



In [ ]:

    
L2=numpy.linspace(-1,1,4)
print(L2)



In [ ]:

    
L3= numpy.logspace(-1,1,4)
print(L3)



In [ ]:

    
print(numpy.meshgrid(L2,L3))



In [ ]:

    
print(numpy.mgrid[0:5,0:5])



In [ ]:

    
print(numpy.ogrid[0:5,0:5])



In [ ]:

    
B=numpy.ones((3,3))
checker2by2=numpy.zeros((6,6))
checker2by2[0:3,0:3]=checker2by2[3:6,3:6]=B



In [ ]:

    
print(numpy.tile(checker2by2,(4,4)))

Routines for the combination of two or more arrays

Loading in memory necessary modules



In [ ]:

    
import scipy.ndimage 
import numpy as np 
import matplotlib.pyplot as plt

Loading in memory the image



In [ ]:

    
text = scipy.ndimage.imread('Chap_02_text_image.png')

Preprocessing the image to bring it as close as possible to the greyscale approximation



In [ ]:

    
text = np.mean(text.astype(float)/255,-1)*2-1

Let’s take a look to the text image



In [ ]:

    
plt.imshow(text)
plt.show()

Find the pattern in the image to be match



In [ ]:

    
letterE = text[37:53,275:291]

Let’s take a look to the text pattern



In [ ]:

    
plt.imshow(letterE)
plt.show()



In [ ]:

    
corr = scipy.ndimage.correlate(text,letterE)
eLocations = (corr >= 0.95 * corr.max())



In [ ]:

    
CorrLocIndex = np.where(eLocations==True)



In [ ]:

    
x=CorrLocIndex[1] 
x



In [ ]:

    
y=CorrLocIndex[0] 
y



In [ ]:

    
thefig=plt.figure() 
plt.subplot(211) 
plt.imshow(text, cmap=plt.cm.gray, interpolation='nearest') 
plt.axis('off') 
plt.subplot(212) 
plt.imshow(text, cmap=plt.cm.gray, interpolation='nearest') 
plt.autoscale(False) 
plt.plot(x,y,'wo',markersize=10) 
plt.axis('off') 
plt.show()

Routines for array manipulation

The next two Ipython inputs are repeated here for completeness (they were introduced earlier in this notebook)



In [ ]:

    
B=numpy.ones((3,3))
checker2by2=numpy.zeros((6,6))



In [ ]:

    
checker2by2[0:3,0:3]=checker2by2[3:6,3:6]=B
numpy.tile(checker2by2,(4,4))



In [ ]:

    
print(checker2by2)



In [ ]:

    
print(numpy.vsplit(checker2by2,3))



In [ ]:

    
a=numpy.array([-numpy.pi, numpy.pi])



In [ ]:

    
print(numpy.vstack((a, numpy.sin(a))))



In [ ]:

    
def max100(x):
    return max(x,100)

Refers to the textbook for why the next commands output a ValueError



In [ ]:

    
max100(a)



In [ ]:

    
numpy.vectorize(max100)(a)

Routines to extract information from arrays



In [ ]:

    
A=numpy.array([5,1,1,2,1,1,2,2,10,3,3,4,5])
numpy.unique(A)



In [ ]:

    
numpy.bincount(A)



In [ ]:

    
A=numpy.fromfunction((lambda i,j: (i+1)*(-1)**(i*j)), (4,4))
print(A)

As discussed in the book, the next command will generate a RuntimeWarning from python computational engine



In [ ]:

    
B=numpy.log2(A)
print(B)



In [ ]:

    
print(numpy.sum(B), numpy.nansum(B))

This is the end of the working codes shown and thoroughly discussed in Chapter 2 of the book Learning SciPy for Numerical and Scientific Computing

Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 Sergio Rojas (srojas@usb.ve) and Erik A Christensen (erikcny@aol.com).