Tuesday, November 22, 2011

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

Monday, November 7, 2011

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