from rsfproj import *

Fetch(['border.hh','elevation.HH',
       'alldata.hh','obsdata.hh',
       'coord.hh','predict.hh'],'rain')

box = 'min1=-185.556 max1=193.18275 min2=-127.262 max2=127.25044'
f2 = 0
def border(name,n2):
    global f2
    Flow([name,name+'x',name+'y'],'border.hh',
         '''
         window n1=1 f1=0 n2=%(n2)d f2=%(f2)d > ${TARGETS[1]} &&
         window n1=1 f1=1 n2=%(n2)d f2=%(f2)d > ${TARGETS[2]} < $SOURCE &&
         cmplx ${TARGETS[1]} ${TARGETS[2]}
         ''' % {'n2':n2,'f2':f2})
    Plot(name,'graph %s title=" " plotcol=6 plotfat=4' % box)
    f2 = f2 + n2

border('border1',338)
border('border2',234)
border('border3',717)
Plot('border',['border1','border2','border3'],'Overlay')

Plot('elev','elevation.HH',
          '''
          igrad |
          grey title=Elevation transp=n yreverse=n
          wantaxis=n wantlabel=n wheretitle=t wherexlabel=b
          ''')
Result('elev',['elev','border'],'Overlay')

select = '''
window n1=1 f1=0 > ${TARGETS[1]} &&
window n1=1 f1=1 > ${TARGETS[2]} < $SOURCE &&
cmplx ${TARGETS[1]} ${TARGETS[2]}
'''

Flow(['alldata','alldatax','alldatay'],'alldata.hh',select)
Result('alldata',
       'graph symbol=x title="All data locations" plotcol=7 ' + box)

Flow(['obs','obsx','obsy'],'obsdata.hh',select)
Plot('obs','graph symbol=o title="Observed data locations" plotcol=5 ' + box)
Plot('obsdata',['obs','border'],'Overlay')

Flow('coord','coord.hh','dd form=native')
Flow(['coord2','coordx','coordy'],'coord',select)
Plot('coord','coord2','graph symbol=x wanttitle=n plotcol=7 ' + box)
Plot('data',['coord','border'],'Overlay')

Result('raindata',['obsdata','data'],'SideBySideIso')

Flow('obsdata','obsdata.hh','dd form=native')
Flow('raindata','obsdata','window f1=2 | put head=$SOURCE')

def contour(cont):
    return '''
    contour nc=6 c0=-0.001 dc=100
    minval=0 maxval=525
    scalebar=y yreverse=n transp=n barlabel=' '
    wanttitle=n plotcol=%d plotfat=%d wantframe=y wantaxis=n
    ''' % (((0,10),(7,1))[cont])

cases = [10,100,1000,10000,-10,-20,-40,0]
for niter in cases:
    lapinter = 'lapinter%d' % niter
    pout = 'pout%d' % niter
    if (niter == 0):
        niter = 100
        name = "Gauss Estimated"
        extra = 'shape=y gauss=y pattern=${SOURCES[1]} eps=1.e-10'
        Flow(lapinter,'raindata lapinter-40',
             '''shapebin niter=%d xkey=0 ykey=1 filt1=25
             nx=376 ny=253 dx=1.00997 dy=1.00997 x0=-185.556 y0=-127.262
             %s''' % (niter,extra))
    else:
        if (niter < 0):
            niter = -niter
            name = "Shaping Regularization"
            extra = 'shape=y pout=${TARGETS[1]}'
            files = [lapinter,pout]
        else:
            prec = 0
            name = "Laplacian Regularization"
            extra = 'shape=n'
            files = lapinter
        Flow(files,'raindata',
             '''
             shapebin niter=%d xkey=0 ykey=1 filt1=25
             nx=376 ny=253 dx=1.00997 dy=1.00997 x0=-185.556 y0=-127.262
             %s
             ''' % (niter,extra))
    Plot('g'+lapinter,lapinter,
              '''grey scalebar=y yreverse=n transp=n allpos=y
              minval=0 maxval=525
              color=j clip=500 title="%s %d" ''' % (name,niter))
    for cont in range(2):
        Plot('c'+str(cont)+lapinter,lapinter,contour(cont))
    Plot(lapinter,['g'+lapinter,'c0'+lapinter,'c1'+lapinter],'Overlay')

reals = []
for seed in range(2000,2005):
    rand = 'rand%d' % seed
    real = 'real%d' % seed
    reals.append(real)
    Flow(rand,'pout-40',
         '''
         noise mean=152.776 var=1e6 rep=y seed=%d
          ''' % seed)
    Flow(real,['raindata',rand],
         '''
         shapebin shape=y pin=${SOURCES[1]}
         niter=100 xkey=0 ykey=1 filt1=25
         nx=376 ny=253 dx=1.00997 dy=1.00997 x0=-185.556 y0=-127.262
         ''')
    Plot('g'+real,real,
         '''grey scalebar=y yreverse=n transp=n allpos=y
         minval=0 maxval=525
         color=j clip=500 title="Realization %d" ''' % (seed-1999))
    for cont in range(2):
        Plot('c'+str(cont)+real,real,contour(cont))
    Plot(real,['g'+real,'c0'+real,'c1'+real],'Overlay')
Result('real',map(lambda seed: 'real%d' % seed,range(2000,2005)),'TwoRows')

Result('lapinter',map(lambda x: 'lapinter%d' % x, cases[:4]),'TwoRows')
Result('gaussinter',map(lambda x: 'lapinter%d' % x, [-10,-40]),
       'SideBySideIso')

Flow('predict','predict.hh','dd form=native')
for niter in [10000,-40,0]:
    lapinter = 'lapinter%d' % niter
    lapdat = 'lapdat%d' % niter
    lapstat = 'lapstat%d' % niter
    Flow(lapdat,[lapinter,'coord'],
         'extract head=${SOURCES[1]} xkey=0 ykey=1')
    if (niter < 0):
        name = "Shaping Regularization"
    else:
        name = "Laplacian Regularization"
    Plot(lapstat,['predict',lapdat],
         '''cmplx ${SOURCES[0:2]} |
         graph symbol="*" title="%s"
         screenratio=1 min1=0 max1=500 min2=0 max2=500
         ''' % name,stdin=0)
Result('stat','lapstat10000 lapstat-40','SideBySideIso')

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

End()