up [pdf]
from rsf.proj import *

# Download data
Fetch('austin-w.HH','bay')

# Convert format
Flow('dem','austin-w.HH','dd form=native')

# Display
def plot(title):
    return '''
    grey clip=250 allpos=y title="%s" 
    screenratio=1
    ''' % title

Result('dem',plot('Digital Elevation Map'))

# Running median program
run = Program('running.c')

w = 30

# !!! CHANGE BELOW !!!
Flow('ave','dem %s' % run[0],
     './${SOURCES[1]} w1=%d w2=%d what=fast' % (w,w))
Result('ave',plot('Signal'))

# Difference
Flow('res','dem ave','add scale=1,-1 ${SOURCES[1]}')
Result('res',plot('Noise') + ' allpos=n')

#############################################################

import sys

if sys.platform=='darwin':
     gtime = WhereIs('gtime')
else:
     gtime = WhereIs('gtime') or WhereIs('time')

if not gtime:
    print '''
For computing CPU time, please install GNU time!
    '''

# slow or fast
for case in ('fast','slow'):

    ts = []
    ws = []

    time = 'time-' + case
    wind = 'wind-' + case



    # loop over window size
    for w in range(3,16,2):
        itime = '%s-%d' % (time,w)
        ts.append(itime)

        iwind = '%s-%d' % (wind,w)
        ws.append(iwind)

        # measure CPU time



        Flow(iwind,None,'spike n1=1 mag=%d' % (w*w))
        Flow(itime,'dem %s' % run[0],
             '''
             ( (%s -f "%%S %%U"
             ./${SOURCES[1]} < ${SOURCES[0]} 
              w1=%d w2=%d what=%s > /dev/null ) 2>&1 )
              > time.out &&
             (tail -1 time.out; 
              echo in=time0.asc n1=2 data_format=ascii_float)
              > time0.asc &&
             dd form=native < time0.asc | stack axis=1 norm=n
             > $TARGET &&
             /bin/rm time0.asc time.out
             ''' % (gtime,w,w,case),stdin=0,stdout=-1)

    Flow(time,ts,'cat axis=1 ${SOURCES[1:%d]}' % len(ts))
    Flow(wind,ws,'cat axis=1 ${SOURCES[1:%d]}' % len(ws))

    # complex numbers for plotting
    Flow('c'+time,[wind,time],
         '''
         cat axis=2 ${SOURCES[1]} |
         transp |
         dd type=complex
         ''')

# Display CPU time
Plot  ('time','ctime-fast ctime-slow',
       '''
       cat axis=1 ${SOURCES[1]} | transp |
       graph dash=0,1 wanttitle=n
       label2="CPU Time" unit2=s
       label1="Window Size" unit1=
       ''',view=1)

End()

sfdd
sfgrey
sfadd
sfspike
sfstack
sfcat
sftransp
sfgraph

data/bay/austin-w.HH