Rescaling

This is an efficient function for rescaling n-dimensional arrays.

def rescale( data, newmin, newmax ):
    oldmin, oldmax = np.min( data ), np.max( data )
    newscale = float( newmax - newmin )
    oldscale = float( oldmax - oldmin )
    if oldscale == 0: return None
    ratio = float( newscale / oldscale )
    c = data.copy()
    c -= oldmin
    c *= ratio
    c += newmin
    return c

Merging Segments in an Image

I'll walk though some code for merging segments in an image using a two classes, Segment and Segments. The first class is for 


class Segment:
    def __init__( self, elements, neighbors, class ):
        self.elements = [ elements ] # tuple addresses
        self.neighbors = neighbors # class labels
        self.class = class # personal class
        self.quality = None # either textured or smooth
        self.mean = None # mean reflectance
        self.texture = None # texture vector


<< OTHER CODE >>

    def merge( self, a, b ):

        # enhance readabilityness

        newclass = a.class

        oldclass = b.class

        # a subsumes b's elements 

        a.elements += b.elements 

        # set of class labels

        # the class labels neighboring a and the class label of a

        mine = set( a.neighbors + [ newclass ] )

        # set of class labels

        # the class labels neighboring b and the class lebel of b

        theirs = set( b.neighbors + [ oldclass ] )

        # set of class labels for elements of a (after adding b to a)

        ours = set( [ self.element_segment[i] for i in a.elements ] )

        # update a's new neighbors

        a.neighbors = list( mine.union( theirs ) - ours ) 

        # update class elementship of b's elements

        for m in b.elements:

            self.element_segment[ m ] = newclass 

            self.segment_element[ newclass ].append( m ) 

        # remove the oldclass from the dict of classes

        del self.segment_element[ oldclass ] 

        # remove the oldclass from the dict of segments

        del self.segments[ oldclass ]

        # update old neighbors of the oldclass

        for s in self.segments:

            if oldclass in self.segments[s].neighbors:

                while oldclass in self.segments[s].neighbors:

                    self.segments[s].neighbors.remove( oldclass )

                if newclass not in self.segments[s].neighbors:

                    self.segments[s].neighbors.append( newclass )

        

What the merge() function does, is figure out the new elements of the new segment (easy), figure out the new neighbors of the new segment (harder), and then do some cleaning operations.

In the merge() function, newclass = a.class and oldclass = b.class take the class labels of the segments to be merged, a and b. In this implementation, a is dominant, and subsumes b's class label.

Then, with a.elements += b.elements, we add the elements of b to the elements of a. These elements are tuple addresses for image chips, or individual pixels.

Imagine the following setup:

|---|---|---|---|

| 3 | 8 | 8 | 8 |

|---|---|---|---|

| 3 | 1 | 2 | 7 |

|---|---|---|---|

| 4 | 1 | 2 | 7 |

|---|---|---|---|

| 4 | 4 | 5 | 6 |

|---|---|---|---|

Here, segment a is denoted by the ones, and segment b is denoted by the twos. The line mine = set( a.neighbors + [ newclass ] ) is just the set {3,4,2,8,1}, the class labels for the n,s,e,w neighbors of segment a, which are {3,4,2,8}, and the class label for a, which is {1}. Similarly for theirs we have the class labels for the neighbors of b, and the class label for b, and this ends up being the set {1,8,7,5,2}. Then ours is the class labels in the new segment, which will be the set {1,2}, since 1 and 2 are the class labels of segment a and segment b.

Then the new neighbors of the new segment will be found by taking the union of mine and theirs, which is {3,8,4,7,5,1,2}, and then tossing ours, giving us the set {3,4,8,7,5}, and these are the tiles ordinally neighboring the tiles labeled 1 and 2. Note that tile 6 is not considered a neighbor because it is to the south-east of tile 2; only tiles that share a whole side are counted as neighbors.

After cleaning up the requisite dictionaries that keep track of these things, we should have the following diagram,

|---|---|---|---|

| 3 | 8 | 8 | 8 |

|---|---|---|---|

| 3 | 1 | 1 | 7 |

|---|---|---|---|

| 4 | 1 | 1 | 7 |

|---|---|---|---|

| 4 | 4 | 5 | 6 |

|---|---|---|---|

pywt versus PyWT

To be specific, I am talking about the difference between pywt (PyWavelets) and PyWT (Python Web Toolkit).

I thought I had installed the wavelets package, as I had used pywt freely for several years; however, when I executed,

$ pip install pywt

I got the web toolkit package instead, even though it goes by PyWF now. Executing the inverese:

$ pip uninstall pywt

I got the reassuring,

$ pip uninstall pywtUninstalling PyWT:Proceed (y/n)? y  Successfully uninstalled PyWT

But when I went to use the wavelet package, it was still accessing the ununinstalled web toolkit, which does absolutely nothing for wavelets. So... I imported pywt at the terminal again and looked at the path:

>>> pywt.__path__
['/usr/local/lib/python2.6/dist-packages/pywt']

And then I went in like gangbangers and shot up that joint with,

$ sudo rm -r /usr/local/lib/python2.6/dist-packages/pywt/

Then everything worked out okay.

    >>> import pywt
    >>> from pywt import dwt
    >>> dwt([1,2,3,4,5,6,7,8],'db4')
    (array([  7.06453146,  4.23073611,  1.41360717,  2.83605428,
              5.6633906 ,  8.49718595, 11.3143149  ]),
     array([  2.37131306e-02,  4.09620864e-02,  -6.46752170e-02,
             -9.90107302e-13, -2.37131306e-02,  -4.09620864e-02,
              6.46752170e-02]))
    >>> 

Matrix Inversion with HDF5 files

#!/usr/bin/env python

import sys
import h5py
import numpy as np



def doolittle_h5( a_filename, side_size ):

import h5py
import numpy as np

fa = h5py.File( a_filename, 'r' )
a = fa['chunked']

fl = h5py.File( 'l.h5', 'r+' )
l = fl[u'data']

fu = h5py.File( 'u.h5', 'r+' )
u = fu[u'data']

n = side_size

for i in range( n ):
    for j in range( i, n ):
        u[i,j] = a[i,j]
        for k in range( i ):
            u[i,j] = u[i,j] - l[i,k] * u[k,j]
    for j in range( i, n ):
        l[j,i] = a[j,i]
        for k in range( i ):
            l[j,i] = l[j,i] - l[j,k] * u[k,i]
        l[j,i] = l[j,i] / float( u[i,i] )
    print 'd: '+i

fa.close()
fl.close()
fu.close()

return None
 

def solve_L_system_h5( a_filename, side_size, b ):

import h5py
import numpy as np

f = h5py.File( a_filename, 'r' )
a = f[u'data'][ :side_size, :side_size ]

x = np.zeros_like( b ).astype( np.float32 )

for i in range( b.size ):
    x[i] = b[i]
    for j in range( i ):
        x[i] -= a[i,j] * x[j]

f.close()

return x
 

def solve_U_system_h5( a_filename, side_size, b ):

import h5py
import numpy as np

f = h5py.File( a_filename, 'r' )
a = f[u'data'][ :side_size, :side_size ]

x = np.zeros_like( b ).astype( np.float32 )

for i in range( b.size-1,-1,-1 ):
    x[i] = b[i]
    for j in range( i+1, b.size ):
        x[i] -= a[i,j] * x[j]
    x[i] /= float( a[i,i] )

f.close()

return x


def inverse_h5( a_filename, side_size ):

import h5py
import numpy as np

doolittle_h5( a_filename, side_size )

print "Done doing little."

finv = h5py.File( 'ainv.h5', 'w' )
dinv = finv.create_dataset( u'data', ( side_size, side_size ), np.float32, chunks=(side_size,1) )
inv = finv[u'data']

for i in range( side_size ):
    b = np.zeros((side_size,)) ; b[i] = 1.0
    b.astype( np.float32 )
    y = solve_L_system_h5( 'l.h5', side_size, b )
    x = solve_U_system_h5( 'u.h5', side_size, y )
    inv[:,i] = x
    print "inverting: "+str(i/float(side_size))+"% done"

finv.close()

return None


#fa = h5py.File( 'a.h5', 'w' )
#d = fa.create_dataset( 'chunked', ( 3,3 ), np.float32 )
#a = fa['chunked']

#da = [ [ 2,-1,-2 ], [ -4,6,3 ], [ -4,-2,8 ] ]
#for i in range( 3 ):
#    for j in range( 3 ):
#        a[i,j] = da[i][j]

#fa.close()

#inverse_h5('a.h5',3)

# at command line: d0=date; nice -n -20 python doolittle.py xx.h5 250000; d1=date;
inverse_h5( sys.argv[1], int( sys.argv[2] ) )

Deleting tons of files at once

Okay, suppose you have tons of files you need to delete.. so many that you can't open up the folder, or use ls to delete them in one fell swoop. Suppose further that your files are labeled like, "x_0_0_1_2.dat". Then you can write a loop in you bash terminal as follows:

$ for i in 0 1 2 ; do rm x_${i}_* ; done

The key thing to note here is the curly braces around the iterating variable. If you leave off the braces, it won't work.

$ for i in 3 4 5 ; do rm x_$i_* ; done
bash: /bin/rm: Argument list too long
bash: /bin/rm: Argument list too long
bash: /bin/rm: Argument list too long

An Online Physics Textbook!

It claims to be a work in progress, but it seems pretty complete to me. You can find it here at http://physics.info/

PyCUDA Documentation

Here is the PyCUDA documentation.