from rsf.proj import*
from rsf.prog import RSFROOT

clip=0.8 #display percentage
nfw=7

def Grey(data,other):
        Result(data,'grey label2=Position unit2="km" label1=Time unit1="s" title="" wherexlabel=b scalebar=y wheretitle=t %s'%other)

##########################################
#    Make synthetic data:data1* & data2
##########################################
Flow('data1',None,
     '''
     spike n1=512 k1=130,150,180,250,300 nsp=5|
     ricker1 frequency=30 |
     spray axis=2 n=256 d=0.01 o=-1 label=Offset unit=km |
     nmostretch inv=y v0=4 half=n |
     scale axis=12
     ''')
Flow('data2',None,
     '''
     spike n1=512 k1=130,150,180,250,300 nsp=5 |
     ricker1 frequency=30 |
     spray axis=2 n=256 d=0.01 o=1 label=Offset unit=km |
     nmostretch inv=y v0=4 half=n |
     scale axis=12
     ''')
#############################################
#               Experiment
#############################################
## Apply dithering (can also use sfblend, 
#  they are equal in the case of small dithering range)
Flow('dither','data1',
     '''
     window n1=1 |
     noise rep=y seed=122011 var=0.01
     ''')
Flow('data2d','data2 dither',
     '''
     datstretch datum=${SOURCES[1]}
     ''')
## Blending 
Flow('datas','data1 data2d','add ${SOURCES[1]}')

# Apply shift and get the blended data2
Flow('udatas','datas dither',
     '''
     datstretch datum=${SOURCES[1]} inv=y
     ''')

#######################################################################
## filtering using median filter and space-varying median filter
#######################################################################
## using median filter
Flow('deblended1mf','datas','transp | mf nfw=%d | transp'%nfw)
Flow('deblended2mf','udatas','transp | mf nfw=%d | transp'%nfw)

## using svmf filter
Flow('deblended1tsmf L1-tsmf0','datas','transp  | tsmf nfw=%d L=${TARGETS[1]} l1=4 | transp'%nfw)
Flow('deblended2tsmf L2-tsmf','udatas','transp | tsmf nfw=%d L=${TARGETS[1]} l1=4  | transp'%nfw)

Flow('deblended1svmf-mf','datas','transp | mf nfw=%d | transp'%nfw)
Flow('deblended2svmf-mf','udatas','transp | mf nfw=%d | transp'%nfw)
Flow('simi1','deblended1svmf-mf datas','similarity other=${SOURCES[1]} rect1=5 rect2=5 | transp  | scale axis=2')
Flow('simi2','deblended2svmf-mf udatas','similarity other=${SOURCES[1]} rect1=5 rect2=5 | transp | scale axis=2')
Result('simi1','transp | grey color=j label2=Position unit2="km" label1=Time unit1="s" title="Signal reliability" wherexlabel=b scalebar=y wheretitle=t scalebar=y allpos=y')

Result('simi2','transp | grey color=j label2=Position unit2="km" label1=Time unit1="s" title="Signal reliability" wherexlabel=b scalebar=y wheretitle=t scalebar=y allpos=y')

Result('energy1-0','deblended1mf','scale axis=2 | math output="abs(input)" | grey color=j label2=Position unit2="km" label1=Time unit1="s" title="Signal energy" wherexlabel=b scalebar=y wheretitle=t scalebar=y allpos=y')

Result('energy2-0','deblended2mf','scale axis=2 | math output="abs(input)"  | grey color=j label2=Position unit2="km" label1=Time unit1="s" title="Signal energy" wherexlabel=b scalebar=y wheretitle=t scalebar=y allpos=y')

Plot('label-e-1',None,
        '''
        box x0=6.1 y0=5.6 label="" xt=-0.5 yt=0.5
        ''')
Plot('label-e-2',None,
        '''
        box x0=9.3 y0=5.6 label="" xt=-0.5 yt=0.5
        ''')
Plot('label-e-3',None,
        '''
        box x0=10.1 y0=4.1 label="" xt=0.5 yt=0.5
        ''')
Plot('label-e-4',None,
        '''
        box x0=4.4 y0=6.5 label="" xt=-0.5 yt=0.5
        ''')
Result('energy1','Fig/energy1-0.vpl label-e-1 label-e-2 label-e-3 label-e-4','Overlay')
Result('energy2','Fig/energy2-0.vpl label-e-1 label-e-2 label-e-3 label-e-4','Overlay')

Flow('deblended1svmf L1','datas simi1','transp |  svmf nfw=%d L=${TARGETS[1]} similarity=${SOURCES[1]} lambda1=0.50 lambda2=0.75 lambda3=0.85 lambda4=0.95 | transp'%nfw)
Flow('deblended2svmf L2','udatas simi2','transp | svmf nfw=%d L=${TARGETS[1]} similarity=${SOURCES[1]} lambda1=0.50 lambda2=0.75 lambda3=0.85 lambda4=0.95 | transp'%nfw)


#Flow('deblended1svmf','datas','transp  | svmf nfw=%d | transp'%nfw)
#Flow('deblended2svmf','udatas','transp | tsmf nfw=%d | transp'%nfw)
## noise section using mf
Flow('diff1mf-t','datas deblended1mf','add scale=1,-1 ${SOURCES[1]}')
Flow('diff2mf-t','udatas deblended2mf','add scale=1,-1 ${SOURCES[1]}')
Flow('diff1svmf-t','datas deblended1svmf','add scale=1,-1 ${SOURCES[1]}')
Flow('diff2svmf-t','udatas deblended2svmf','add scale=1,-1 ${SOURCES[1]}')

## noise section using svmf
Flow('diff1svmf','datas deblended1svmf','add scale=1,-1 ${SOURCES[1]}')
Flow('diff2svmf','udatas deblended2svmf','add scale=1,-1 ${SOURCES[1]}')

## frame
Plot('label1',None,
        '''
        box x0=11.5 y0=6.5 label="" xt=-0.5 yt=0.5
        ''')
Plot('label2',None,
        '''
        box x0=10.5 y0=5.5 label="Signal" xt=0.5 yt=0.5
        ''')

## Ploting
Grey('data1','scalebar=n clip=%f   '%clip)
Grey('data2','scalebar=n clip=%f   '%clip)
Grey('datas','scalebar=n clip=%f   '%clip)
Grey('udatas','scalebar=n clip=%f   '%clip)
Grey('deblended1mf','scalebar=n clip=%f   '%clip)
Grey('deblended2mf','scalebar=n clip=%f   '%clip)
Grey('deblended1tsmf','scalebar=n clip=%f   '%clip)
Grey('deblended2tsmf','scalebar=n clip=%f   '%clip)
Grey('deblended1svmf','scalebar=n clip=%f   '%clip)
Grey('deblended2svmf','scalebar=n clip=%f   '%clip)

Grey('diff1mf-t','scalebar=n clip=%f   '%clip)
Result('diff1mf','Fig/diff1mf-t.vpl label1','Overlay')
Grey('diff1svmf-t','scalebar=n clip=%f   '%clip)
Result('diff1svmf','Fig/diff1svmf-t.vpl label1','Overlay')
Result('L1','transp | grey label2=Position unit2="km" label1=Time unit1="s" title="Filter length" wherexlabel=b scalebar=y wheretitle=t scalebar=y color=j allpos=y')
Result('L2','transp | grey label2=Position unit2="km" label1=Time unit1="s" title="Filter length" wherexlabel=b scalebar=y wheretitle=t scalebar=y color=j allpos=y')

Result('L1-tsmf0','transp | grey label2=Position unit2="km" label1=Time unit1="s" title="Filter length" wherexlabel=b scalebar=y wheretitle=t scalebar=y color=j allpos=y')
Result('L1-tsmf','Fig/L1-tsmf0.vpl label-e-1 label-e-2 label-e-3 label-e-4','Overlay')

Result('L2-tsmf','transp | grey label2=Position unit2="km" label1=Time unit1="s" title="Filter length" wherexlabel=b scalebar=y wheretitle=t scalebar=y color=j allpos=y')

Grey('diff2mf-t','scalebar=n clip=%f   '%clip)
Result('diff2mf','Fig/diff2mf-t.vpl label2','Overlay')
Grey('diff2svmf-t','scalebar=n clip=%f   '%clip)
Result('diff2svmf','Fig/diff2svmf-t.vpl label2','Overlay')

Flow('fkdata1','data1','put d2=1 | fft1 | fft3 axis=2 pad=1')
Flow('fkdatas','datas','put d2=1 | fft1 | fft3 axis=2 pad1=1')
Flow('fkdeblended1mf','deblended1mf','put d2=1 | fft1 | fft3 axis=2 pad1=1')
Flow('fkdeblended1svmf','deblended1svmf','put d2=1 | fft1 | fft3 axis=2 pad1=1 ')
Result('fkdata1','cabs | grey title="" label2="Normalized wavenumber" unit2="" color=j ')
Result('fkdatas','cabs | grey title="" label2="Normalized wavenumber" unit2="" color=j ')
Result('fkdeblended1mf','cabs | grey title="" label2="Normalized wavenumber" unit2="" color=j ')
Result('fkdeblended1svmf','cabs | grey title="" label2="Normalized wavenumber" unit2="" color=j ')


End()