from rsf.proj import *

zmax = 171
xmin = 315
xmax = 930

minv = 1.5
maxv = 3.3
clip = 1.85

# fetch model
Fetch('marmvel.hh',"marm")

Flow('marmvel','marmvel.hh',
     '''
     dd form=native | 
     scale rscale=.001 | put
     label1=Depth label2=Position unit1=km unit2=km
     d1=0.004 d2=0.004
     ''')

# decimate model
Flow('marm','marmvel','window j1=2 j2=2')

Plot('marm',
     '''
     window n1=%d n2=%d f2=%d |
     grey color=j scalebar=y allpos=y title="Marmousi Model"
     barlabel=Velocity barunit=km/s barreverse=y
     minval=%g maxval=%g clip=%g bias=%g labelsz=10 titlesz=12 labelfat=6 titlefat=6
     ''' % (zmax,xmax-xmin+1,xmin,minv,maxv,clip,minv))

# prior model
Flow('init','marm','math output="1.5+x1"')

Plot('init',
     '''
     window n1=%d n2=%d f2=%d |
     grey color=j scalebar=y allpos=y title="Initial Model"
     barlabel=Velocity barunit=km/s barreverse=y
     minval=%g maxval=%g clip=%g bias=%g labelsz=10 titlesz=12 labelfat=6 titlefat=6
     ''' % (zmax,xmax-xmin+1,xmin,minv,maxv,clip,minv))

cgiter = 10
niter  = 4
eps    = 5

# DSR tomography
################

# data
Flow('data','marm','dsreiko')

Result('data',
       '''
       byte gainpanel=all bar=bar.rsf allpos=y pclip=100 | 
       grey3 color=j scalebar=y flat=n point1=0.35 point2=0.7
       title="DSR Prestack Modeling" barlabel=Traveltime barunit=s 
       label2=Receiver unit2=km label3=Source unit3=km
       screenratio=0.7 labelsz=6 titlesz=8 labelfat=3 titlefat=3
       ''')

# mask
Flow('mask','data',
     '''
     window n1=1 | math output="abs(x1-x2)" | mask max=6.004 | 
     window j2=4 | lpad jump=4 | window n2=1151
     ''')

# noise
Flow('nmask','data',
     '''
     window n1=1 | math output="abs(x1-x2)" | mask min=1.396 max=6.004 | 
     window j2=4 | lpad jump=4 | window n2=1151 | dd type=float
     ''')
Flow('noise','nmask',
     '''
     noise seed=20130516 range=0.6 rep=y |
     thr mode=hard thr=0.25 | add mode=p ${SOURCES[0]}
     ''')

Flow('ndata','noise data',
     '''
     spray axis=1 n=376 d=0.008 o=0 |
     add scale=1,1 ${SOURCES[1]} | math output="abs(input)"
     ''')

# DSR Tickhnov
Flow('dsrt gdsrt','init data mask',
     '''
     dsrtomo reco=${SOURCES[1]} mask=${SOURCES[2]} grad=${TARGETS[1]}
     verb=y shape=n cgiter=%d niter=%d eps=%g
     ''' % (cgiter,niter,eps))
Flow('ndsrt gndsrt','init ndata mask',
     '''
     dsrtomo reco=${SOURCES[1]} mask=${SOURCES[2]} grad=${TARGETS[1]}
     verb=y shape=n cgiter=%d niter=%d eps=%g
     ''' % (cgiter,niter,eps))

Plot('dsrt',
     '''
     window n1=%d n2=%d f2=%d |
     grey color=j allpos=y scalebar=y title="DSR Tomography"
     barlabel=Velocity barunit=km\/s barreverse=y
     minval=%g maxval=%g clip=%g bias=%g labelsz=10 titlesz=12 labelfat=6 titlefat=6
     ''' % (zmax,xmax-xmin+1,xmin,minv,maxv,clip,minv))
Plot('ndsrt',
     '''
     window n1=%d n2=%d f2=%d |
     grey color=j allpos=y scalebar=y title="DSR Tomography"
     barlabel=Velocity barunit=km\/s barreverse=y
     minval=%g maxval=%g clip=%g bias=%g labelsz=10 titlesz=12 labelfat=6 titlefat=6
     ''' % (zmax,xmax-xmin+1,xmin,minv,maxv,clip,minv))

Flow('cdsrt',None,
     '''
     echo 1. 0.370197 0.262643 0.149201 0.078508
     n1=5 o1=0 d1=1 data_format=ascii_float in=$TARGET
     ''')
Plot('cdsrt',
     '''
     graph wanttitle=n
     plotfat=4 labelsz=7 titlesz=9 labelfat=3 titlefat=3
     label1="Number of Linearization Updates" unit1= label2="Normalized l2 Misfit" unit2=
     screenratio=0.8 screenht=7.5
     ''')

Flow('cndsrt',None,
     '''
     echo 1. 0.742392 0.715866 0.700967 0.694760
     n1=5 o1=0 d1=1 data_format=ascii_float in=$TARGET
     ''')
Plot('cndsrt',
     '''
     graph wanttitle=n
     plotfat=4 labelsz=7 titlesz=9 labelfat=3 titlefat=3
     label1="Number of Linearization Updates" unit1= label2="Normalized l2 Misfit" unit2=
     screenratio=0.8 screenht=7.5
     ''')

# STD tomography
################

# shots
Flow('stemp2',None,
     'math n1=288 d1=0.032 o1=0. output=x1 | transp')
Flow('stemp1','stemp2','math output=0.')
Flow('shots','stemp1 stemp2',
     'cat axis=1 ${SOURCES[1]} ${SOURCES[0]}')

# data and mask
Flow('receiver toporec','marm shots',
     '''
     mkrcv shot=${SOURCES[1]} reco=${TARGETS[1]} offset1=750
     ''')

# noise
Flow('ntemp','noise','window j2=4')

ntoporec=[]
for n in range(288):
    ntoporec.append('ntoporec_'+str(n))
    start=max(n*0.032-6.004,0.)

    Flow('ntoporec_'+str(n),'ntemp',
         '''
         window n2=1 f2=%d | window min1=%g |
         pad n1out=1151 | put o1=0.
         ''' % (n,start))

Flow('ntoporec0',ntoporec,'cat axis=2 ${SOURCES[1:%d]}' % len(ntoporec))
Flow('ntoporec','ntoporec0 toporec','add scale=1,1 ${SOURCES[1]} | math output="abs(input)"')

# STD Tickhnov
Flow('stdt gstdt','init shots toporec receiver',
     '''
     fatomoomp shot=${SOURCES[1]} reco=${SOURCES[2]} recv=${SOURCES[3]}
     grad=${TARGETS[1]} shape=n stiter=%d niter=%d eps=%g verb=y
     ''' % (cgiter,niter,eps))
Flow('nstdt ngstdt','init shots ntoporec receiver',
     '''
     fatomoomp shot=${SOURCES[1]} reco=${SOURCES[2]} recv=${SOURCES[3]}
     grad=${TARGETS[1]} shape=n stiter=%d niter=%d eps=%g verb=y
     ''' % (cgiter,niter,eps))

Plot('stdt',
     '''
     window n1=%d n2=%d f2=%d |
     grey color=j allpos=y scalebar=y title="Standard Tomography"
     barlabel=Velocity barunit=km\/s barreverse=y
     minval=%g maxval=%g clip=%g bias=%g labelsz=10 titlesz=12 labelfat=6 titlefat=6
     ''' % (zmax,xmax-xmin+1,xmin,minv,maxv,clip,minv))
Plot('nstdt',
     '''
     window n1=%d n2=%d f2=%d |
     grey color=j allpos=y scalebar=y title="Standard Tomography"
     barlabel=Velocity barunit=km\/s barreverse=y
     minval=%g maxval=%g clip=%g bias=%g labelsz=10 titlesz=12 labelfat=6 titlefat=6
     ''' % (zmax,xmax-xmin+1,xmin,minv,maxv,clip,minv))

Flow('cstdt',None,
     '''
     echo 1. 0.375924 0.342480 0.234603 0.067037
     n1=5 o1=0 d1=1 data_format=ascii_float in=$TARGET
     ''')
Plot('cstdt','graph wanttitle=n wantaxis=n dash=1 plotfat=4 screenratio=0.8 screenht=7.5')

Flow('cnstdt',None,
     '''
     echo 1. 0.726534 0.701836 0.689851 0.671497
     n1=5 o1=0 d1=1 data_format=ascii_float in=$TARGET
     ''')
Plot('cnstdt','graph wanttitle=n wantaxis=n dash=1 plotfat=4 screenratio=0.8 screenht=7.5')

# plot
######
Result('marm','marm init','OverUnderAniso')
Result('tomo','stdt dsrt','OverUnderAniso')
Result('ntomo','nstdt ndsrt','OverUnderAniso')
Result('conv','cstdt cdsrt','Overlay')
Result('nconv','cnstdt cndsrt','Overlay')

#Flow('cdsrt0','cdsrt','dd form=native')
#Flow('cstdt0','cstdt','dd form=native')

#Result('conv','cdsrt0 cstdt0',
#       '''
#       cat axis=2 ${SOURCES[1]} | pygraph
#       label1="Normalized l2 Misfit" unit1=
#       label2="Number of Linear Updates" unit2= 
#       legend="DSR tomography:standard tomography" wherelegend=1
#       plotcol=0 dash=0,1 min1=0 min2=0 plotfat=3 format="eps"
#       ''', suffix='.eps')

End()