Saturday, February 25, 2012

Defining functions for TikZ

The "easiest" way I've found to define "functions" in TikZ is to use the LaTeX \newcommand outside of the \tikzpicture environment.
\newcommand{\cir}[3]{ \def \x{ (#1,#2) circle (#3cm) } \draw \x ; }
And then use this inside the tikzpicture as:
\cir{0}{0}{1}
However, since the circle is drawn within the \newcommand, I can't rotate the figure inside the \tikzpicture environment.

Tuesday, December 20, 2011

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

Thursday, December 8, 2011

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







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

Wednesday, October 19, 2011

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


Friday, April 22, 2011

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/