Mailman 3 January 2014 - NumPy-Discussion

Deprecation of boolean substract and negative (the - operator)
by sebastian 30 Jan '14

30 Jan '14

Hey all, recently we had a small discussion about deprecating some of the operators for boolean arrays. This discussion seemed to have ended by large in the consense that while most boolean operators are well defined and should be kept, the `-` one is not very well defined on boolean arrays and has the problem of the inconsistency: - np.array(False) == True False - np.array(False) == False # leading to: False - (-np.arry(False)) != False + np.array(False) So that it is preferable to use one of the binary operators for this operation. For now this would only be a deprecation, but both operators are probably used out there. So if you have any serious doubt about starting this deprecation please note it here. The Pull request to implement such a deprecation is: https://github.com/numpy/numpy/pull/4105 Regards, Sebastian

1 0

ENH: Type specific binary search functions for `searchsorted`
by Jaime Fernández del Río 30 Jan '14

30 Jan '14

Hi, I have just added a new PR: https://github.com/numpy/numpy/pull/4244 >From the commit message: This PR replaces the generic binary search functions used by `searchsorted` with type specific ones for numeric types. This results in a speed-up of calls to `searchsorted` which is highly dependent on the size of the 'haystack' and the 'needle', with typical values for large enough needles in the 1.5x - 3.0x for direct searches (i.e. without a sorter argument) and 1.2x - 2.0x for indirect searches. A summary benchmark on float and int arrays can be found here. Furthermore, the type specific binary search functions can take strided inputs for all their arguments, which is a step in the right direction to eventually add an axis argument to`searchsorted`. Any comments and reviews are more than welcome! Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

1 0

Memory leak in numpy?
by Chris Laumann 29 Jan '14

29 Jan '14

Hi all- I think I just found a memory leak in numpy, or maybe I just don’t understand generators. Anyway, the following snippet will quickly eat a ton of RAM: P = randint(0,2, (20,13)) for i in range(50): for ai in ndindex((2,)*13): j = P.dot(ai) If you replace the last line with something like j = ai, the memory leak goes away. I’m not exactly sure what’s going on but the .dot seems to be causing the memory taken by the tuple ai to be held. This devours RAM in python 2.7.5 (OS X Mavericks default I believe), numpy version 1.8.0.dev-3084618. I’m upgrading to the latest Superpack (numpy 1.9) right now but I somehow doubt this behavior will change. Any thoughts? Best, Chris -- Chris Laumann Sent with Airmail

5 6

Re: [Numpy-discussion] scipy image processing memory leak in python 2.7?
by Joseph McGlinchy 29 Jan '14

29 Jan '14

Hi numpy list! I am trying to do some image processing on a number of images, 72 to be specific. I am seeing the python memory usage continually increase. I will say this is being done also using arcpy so that COULD be an issue, but want to rule out scipy first since I am only using arcpy to grab a set of polygons and convert them to a binary image, and to save images to disk. I pass the raster to my numpy/scipy processing. The code is organized as: 1) Main loop iterating over features 2) Function which does image processing using scipy 3) Class which implements some of the image processing features I placed 2 sleep calls in the code.. one at the end of the image processing function and one at the end of each iteration in the main loop. I am seeing a 20MB memory increase each iteration while looking at my processes in Task Manager: Iteration 1: 204,312K and 209,908K, respectively Iteration 2: 233,728K and 230,192K, respectively Iteration 3: 255,676K and 250,600K, respectively Any ideas? Thanks so much! Joe 1) Main loop #for each feature class, convert to raster by "Value" and reclass to 1 and 0 temp_raster = os.path.join(scratchWorkspace,"temp_1.tif") temp_raster_elev = os.path.join(scratchWorkspace,"temp_elev.tif") temp_reclass = os.path.join(scratchWorkspace,"temp_reclass.tif") remap = arcpy.sa.RemapValue([[1,1], ["NoData",0]]) for i,fc in enumerate(potWaterFeatNames): #try to delete the temp rasters in case they exist try: #delete temp rasters arcpy.Delete_management(temp_raster) arcpy.Delete_management(temp_raster_elev) arcpy.Delete_management(temp_reclass) except e: arcpy.AddMessage(e.message()) #clear in memory workspace" arcpy.Delete_management("in_memory/temp_table") #print "on feature {} of {}".format(i+1, len(potWaterFeatNames)) arcpy.AddMessage("on feature {} of {}".format(i+1, len(potWaterFeatNames))) arcpy.AddMessage("fc = {}".format(fc)) arcpy.PolygonToRaster_conversion(fc, "Value", temp_raster, "", "", 1) outReclass = arcpy.sa.Reclassify(temp_raster,"Value",remap) outReclass.save(temp_reclass) del outReclass #call function to process connected components try: proc_im = processBinaryImage(temp_reclass, scratchWorkspace, sm_pixels) except Exception as e: print "FAILED! scipy image processing" arcpy.AddMessage(e.message()) #convert processed raster to polygons try: out_fc = os.path.join(outWorkspace,fc + "_cleaned") arcpy.RasterToPolygon_conversion(proc_im, out_fc) arcpy.Delete_management(proc_im) except Exception as e: print "FAILED! converting cleaned binary raster to polygon" arcpy.AddMessage(e.message()) #delete zero value polyons, gridcode = 0 try: uc = arcpy.UpdateCursor(out_fc, "gridcode=0",arcpy.Describe(out_fc).spatialReference) #loop through 0 rows and delete for row in uc: uc.deleteRow(row) del uc except Exception as e: print "FAILED! deleting rows from table" arcpy.AddMessage(e.message()) #check that number of polygons with gridcode = 1 is greater than 0 count = arcpy.GetCount_management(out_fc) #add elevation field back in if (count > 0): arcpy.PolygonToRaster_conversion(fc, "z_val", temp_raster_elev, "", "", 1) arcpy.sa.ZonalStatisticsAsTable(out_fc, "Id", temp_raster_elev, "in_memory/temp_table","#","MEAN") arcpy.JoinField_management(out_fc, "Id", "in_memory/temp_table", "Id", ["MEAN"]) else: arcpy.Delete_management(out_fc) #delete temp rasters arcpy.Delete_management(temp_raster) arcpy.Delete_management(temp_raster_elev) arcpy.Delete_management(temp_reclass) #python garbage collection #collected = gc.collect() #print "Garbage collector: collected %d objects." % (collected) print "sleeping for 10 seconds" time.sleep(10) 2) Image processing function def processBinaryImage(imageName, save_dir, sm_pixels): fname = os.path.basename(imageName) imRaster = arcpy.Raster(imageName) #Grab AGS info from image #use describe module to grab info descData = arcpy.Describe(imRaster) cellSize = descData.meanCellHeight extent = descData.Extent spatialReference = descData.spatialReference pnt = arcpy.Point(extent.XMin, extent.YMin) del imRaster converter = ConvertRasterNumpy(); imArray = converter.rasterToNumpyArray(imageName) imMin = np.min(imArray) imMax = np.max(imArray) print imMin, imMax arcpy.AddMessage("image min: " + str(imMin)) arcpy.AddMessage("image max: " + str(imMax)); #other flags show_image = False #verbose save_flag = False #under verbose gray = False #threshold but keep gray levels make_binary = False #flag if grayscale is true, to binarize final image save_AGS_rasters = True #create filter object filter = thresholdMedianFilter() lowValue = 0 tImage = filter.thresholdImage(imArray, lowValue, gray) #median filter image filtImage1 = filter.medianFiltImage(tImage,1) #create structuring element for morphological operations sElem = np.array([[0., 1., 0.], [1., 1., 1.], [0., 1., 0.]]) #open filtered image gray = False mImage = filter.mOpenImage(filtImage1, sElem, gray) #set list to store info change num_its = 100 the_change = np.zeros(num_its) for it in range(100): prev = mImage filtImage = filter.medianFiltImage(mImage,1) mImage = filter.mOpenImage(filtImage, sElem, gray) #calculate difference (m_info, z_info) = filter.calcInfoChange(prev, mImage) the_change[it] = z_info del filtImage del prev #if the change is less than 5% of the initial change, exit cur_perc = the_change[it]/the_change[0]*100 if cur_perc < 5.0: print "exiting filter on iteration " + str(it) print "change is less than 5% (this iteration: " + str(cur_perc) + ")" break ############################################################################ # now we have a binary mask. Let's find and label the connected components # ############################################################################ #clear some space del tImage del filtImage1 del m_info label_im_init, nb_labels = ndimage.label(mImage) #Compute size, mean_value, etc. of each region: label_im = label_im_init.copy() del label_im_init sizes = ndimage.sum(mImage, label_im, range(nb_labels + 1)) mean_vals = ndimage.sum(imArray, label_im, range(1, nb_labels + 1)) #clean up small components mask_size = sizes < sm_pixels remove_pixel = mask_size[label_im] label_im[remove_pixel] = 0; labels = np.unique(label_im) label_im = np.searchsorted(labels, label_im) #make label image to a binary image, and convert to arcgis raster label_im[label_im > 0] = 1 label_im = np.array(label_im, dtype = 'float32') print label_im.dtype saveit = False if ~saveit: outRaster = save_dir + "\\" + fname[:-4] + "_CC_" + str(sm_pixels) + ".tif" temp = arcpy.NumPyArrayToRaster(label_im,pnt,cellSize,cellSize) arcpy.DefineProjection_management(temp, spatialReference) arcpy.CopyRaster_management(temp, outRaster, "DEFAULTS","0","","","","8_BIT_UNSIGNED") #clear more space del mImage del nb_labels del sizes del mean_vals del mask_size del remove_pixel del label_im del labels del temp del the_change del sElem del filter del imArray print 'sleeping' time.sleep(20) return outRaster 3) Image processing class class thresholdMedianFilter(): def thresholdImage(self, imArray, thresh, binary = True): """ threshold the image. values equal or below thresh are 0, above are 1""" tImage = imArray tImage[tImage <= thresh] = 0 if binary: tImage[tImage > thresh] = 1 return tImage def medianFiltImage(self,imArray,n=1): """median filter the image n amount of times. a single time is the default""" for n in range(n): prev = imArray imArray = signal.medfilt2d(prev) del prev return imArray def mOpenImage(self,imArray, sElem, gray = False): """ morphological opening """ #Mnp = np.array(morph.binary_erosion(imArray, sElem)) #Mnp = np.array(morph.binary_dilation(Mnp, sElem)) if gray: imArray1 = np.array(morph.grey_dilation(imArray, structure = sElem)) imArray2 = np.array(morph.grey_erosion(imArray1, structure = sElem), dtype = 'float32') del imArray1 return imArray2 else: Mnp1 = np.array(morph.binary_dilation(imArray, sElem)) Mnp2 = np.array(morph.binary_erosion(Mnp1, sElem), dtype = 'float32') del Mnp1 return Mnp2 def calcInfoChange(self,imArray1, imArray2): """calculate entropy of an image""" diff = imArray1 - imArray2 m_norm = sum(abs(diff)) #Manhattan Norm z_norm = norm(diff.ravel(), 0) #Zero norm del diff return (m_norm, z_norm) Joe McGlinchy | Imagery Scientist Database Services - 3D and Imagery Team ESRI || 380 New York St. || Redlands, CA 92373 || USA T 909-793-2853, ext. 4783

3 4

windows and C99 math
by Julian Taylor 27 Jan '14

27 Jan '14

hi, numpys no-C99 fallback keeps turning up issues in corner cases, e.g. hypot https://github.com/numpy/numpy/issues/2385 log1p https://github.com/numpy/numpy/issues/4225 these only seem to happen on windows, on linux and mac it seems to use the C99 math library just fine. Are our binary builds for windows not correct or does windows just not support C99 math? Hopefully it is the former. Any insight is appreciated (and patches to fix the build even more!) Cheers, Julian

4 3

bug in comparing object arrays to None (?)
by Charles G. Waldman 27 Jan '14

27 Jan '14

Hi Numpy folks. I just noticed that comparing an array of type 'object' to None does not behave as I expected. Is this a feature or a bug? (I can take a stab at fixing it if it's a bug, as I believe it is). >>> np.version.full_version '1.8.0' >>> a = np.array(['Frank', None, 'Nancy']) >>> a array(['Frank', None, 'Nancy'], dtype=object) >>> a == 'Frank' array([ True, False, False], dtype=bool) # Return value is an array >>> a == None False # Return value is scalar (BUG?)

2 1

another interesting high performance vector lib (yeppp)
by Neal Becker 27 Jan '14

27 Jan '14

http://www.yeppp.info/

2 1

ANN: numexpr 2.3 (final) released
by Francesc Alted 27 Jan '14

27 Jan '14

========================== Announcing Numexpr 2.3 ========================== Numexpr is a fast numerical expression evaluator for NumPy. With it, expressions that operate on arrays (like "3*a+4*b") are accelerated and use less memory than doing the same calculation in Python. It wears multi-threaded capabilities, as well as support for Intel's MKL (Math Kernel Library), which allows an extremely fast evaluation of transcendental functions (sin, cos, tan, exp, log...) while squeezing the last drop of performance out of your multi-core processors. Look here for a some benchmarks of numexpr using MKL: https://github.com/pydata/numexpr/wiki/NumexprMKL Its only dependency is NumPy (MKL is optional), so it works well as an easy-to-deploy, easy-to-use, computational engine for projects that don't want to adopt other solutions requiring more heavy dependencies. Numexpr is already being used in a series of packages (PyTables, pandas, BLZ...) for helping doing computations faster. What's new ========== The repository has been migrated to https://github.com/pydata/numexpr. All new tickets and PR should be directed there. Also, a `conj()` function for computing the conjugate of complex arrays has been added. Thanks to David Menéndez. See PR #125. Finallly, we fixed a DeprecationWarning derived of using ``oa_ndim == 0`` and ``op_axes == NULL`` with `NpyIter_AdvancedNew()` and NumPy 1.8. Thanks to Mark Wiebe for advise on how to fix this properly. Many thanks to Christoph Gohlke and Ilan Schnell for his help during the testing of this release in all kinds of possible combinations of platforms and MKL. In case you want to know more in detail what has changed in this version, see: https://github.com/pydata/numexpr/wiki/Release-Notes or have a look at RELEASE_NOTES.txt in the tarball. Where I can find Numexpr? ========================= The project is hosted at GitHub in: https://github.com/pydata/numexpr You can get the packages from PyPI as well (but not for RC releases): http://pypi.python.org/pypi/numexpr Share your experience ===================== Let us know of any bugs, suggestions, gripes, kudos, etc. you may have. Enjoy data! -- Francesc Alted

2 3

Numpy Enhancement Proposal: group_by functionality
by Eelco Hoogendoorn 26 Jan '14

26 Jan '14

Hi all, Please critique my draft exploring the possibilities of adding group_by support to numpy: http://pastebin.com/c5WLWPbp In nearly ever project I work on, I require group_by functionality of some sort. There are other libraries that provide this kind of functionality, such as pandas for instance, but I will try to make the case here that numpy ought to have a solid core of group_by functionality. Primarily, one may argue that the concept of grouping values by a key is far more general than a pandas dataframe. In particular, one often needs a simple one-line transient association between some keys and values, and trying to wrangle your problem into the more permanent and specialized datastructure that a dataframe is, is simply not called for. As a simple compact example: key1 = list('abaabb') key2 = np.random.randint(0,2,(6,2)) values = np.random.rand(6,3) print group_by((key1, key2)).median(values) Points of note; we can group by arbitrary combinations of keys, and subarrays can also act as keys. group_by has a rich set of reduction functionality, which performs efficient per-group reductions, as well as various ways to split your values per group. Also, the code here has a lot of overlap with np.unique and related arraysetops. functions like np.unique are easily reimplemented using the groundwork laid out here, and also may be extended to benefit from the generalizations made, allowing for a wider variety of objects to have their unique values taken; note the axis keyword here, meaning that what is unique here are the images found along the first axis; not the elements of shuffled. #create a stack of images images = np.random.rand(4,64,64) #shuffle the images; this is a giant mess now; how to find all the original ones? shuffled = images[np.random.randint(0,4,200)] #there you go print unique(shuffled, axis=0) Some more examples and unit tests can be found at the end of the module. Id love to hear your feedback on this. Specifically: - Do you agree numpy would benefit from group_by functionality? - Do you have suggestions for further generalizations/extensions? - Any commentary on design decisions / implementation? Regards, Eelco Hoogendoorn

3 7

ANN: BLZ 0.6.1 has been released
by Francesc Alted 26 Jan '14

26 Jan '14

Announcing BLZ 0.6 series ========================= What it is ---------- BLZ is a chunked container for numerical data. Chunking allows for efficient enlarging/shrinking of data container. In addition, it can also be compressed for reducing memory/disk needs. The compression process is carried out internally by Blosc, a high-performance compressor that is optimized for binary data. The main objects in BLZ are `barray` and `btable`. `barray` is meant for storing multidimensional homogeneous datasets efficiently. `barray` objects provide the foundations for building `btable` objects, where each column is made of a single `barray`. Facilities are provided for iterating, filtering and querying `btables` in an efficient way. You can find more info about `barray` and `btable` in the tutorial: http://blz.pydata.org/blz-manual/tutorial.html BLZ can use numexpr internally so as to accelerate many vector and query operations (although it can use pure NumPy for doing so too) either from memory or from disk. In the future, it is planned to use Numba as the computational kernel and to provide better Blaze (http://blaze.pydata.org) integration. What's new ---------- BLZ has been branched off from the Blaze project (http://blaze.pydata.org). BLZ was meant as a persistent format and library for I/O in Blaze. BLZ in Blaze is based on previous carray 0.5 and this is why this new version is labeled 0.6. BLZ supports completely transparent storage on-disk in addition to memory. That means that *everything* that can be done with the in-memory container can be done using the disk as well. The advantages of a disk-based container is that the addressable space is much larger than just your available memory. Also, as BLZ is based on a chunked and compressed data layout based on the super-fast Blosc compression library, the data access speed is very good. The format chosen for the persistence layer is based on the 'bloscpack' library and described in the "Persistent format for BLZ" chapter of the user manual ('docs/source/persistence-format.rst'). More about Bloscpack here: https://github.com/esc/bloscpack You may want to know more about BLZ in this blog entry: http://continuum.io/blog/blz-format In this version, support for Blosc 1.3 has been added, that meaning that a new `cname` parameter has been added to the `bparams` class, so that you can select you preferred compressor from 'blosclz', 'lz4', 'lz4hc', 'snappy' and 'zlib'. Also, many bugs have been fixed, providing a much smoother experience. CAVEAT: The BLZ/bloscpack format is still evolving, so don't trust on forward compatibility of the format, at least until 1.0, where the internal format will be declared frozen. Resources --------- Visit the main BLZ site repository at: http://github.com/ContinuumIO/blz Read the online docs at: http://blz.pydata.org/blz-manual/index.html Home of Blosc compressor: http://www.blosc.org User's mail list: blaze-dev(a)continuum.io ---- Enjoy! Francesc Alted Continuum Analytics, Inc.

3 2