up [pdf]
from rsf.proj import *

# Download data
Fetch('bay.h','bay')

# Convert format
Flow('bay','bay.h',
     '''
     dd form=native |
     window f2=500 n2=1500 
     ''')

# Display
def plot(title):
    return '''
    grey allpos=y title="%s" yreverse=n 
    label1=South-North label2=West-East 
    screenratio=0.8
    ''' % title

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

# Program for running average
run = Program('run.c')

# COMMENT ABOVE AND UNCOMMENT BELOW IF YOU WANT TO USE PYTHON 
# run = Command('run.exe','run.py','cp $SOURCE $TARGET')
# AddPostAction(run,Chmod(run,0o755))

w = 30

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

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

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

import sys

if sys.platform=='darwin':
     gtime = WhereIs('gtime')
     if not gtime:
          print("For computing CPU time, install gtime.")
else:
     gtime = WhereIs('gtime') or WhereIs('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,'bay %s' % run[0],
             '''
             ( (%s -f "%%S %%U"
             ./${SOURCES[1]} < ${SOURCES[0]} 
              w1=%d w2=%d how=%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
sfwindow
sfgrey
sfadd
sfspike
sfstack
sfcat
sftransp
sfgraph

data/bay/bay.h