from rsf.proj import *
from math import *
import os

wf = 2*pi
nt = 501
dt = 0.004
ot = 0
nx = 501
dx = 0.01
ox = 0
nw = 200
dw = 0.0005
ow = 0


for eve in (1,2,3,4):
    spike='spike%d' % eve
    tpara='tpara%d'   % eve
    para='para%d'     % eve
    Flow(spike,None,
        '''
        spike n1=%d d1=%g o1=%g n2=%d d2=%g o2=%g nsp=1 k1=%d mag=1  p2=0| 
        ricker1 frequency=15 | put unit2=km label2=Distance
        ''' % (nt,dt,ot,nx,dx,ox,eve*80-30))

    Flow(tpara,spike,
        '''
        window n1=1 | math output="-sqrt(%g*%g+(x1-2.5)*(x1-2.5)/%g/%g)+%g"
        ''' % (0.004*(eve*80-30),0.004*(eve*80-30),2,2,0.004*(eve*80-30)))

    Flow(para,[spike, tpara],
        'datstretch datum=${SOURCES[1]} ')

Flow('para','para1 para2 para3 para4','add ${SOURCES[1]} ${SOURCES[2]} ${SOURCES[3]}')

Result('para','para','window j2=4 | wiggle label2=Distance unit2=km transp=y yreverse=y poly=y  title="Signal" ')


Flow('npara','para','noise var=1e-2')
Flow('nft','npara','fft1 ')

Result('npara','npara','window j2=5 | wiggle label2=Distance unit2=km transp=y yreverse=y poly=y  title="Noisy Signal" ')

nshifts = []
for s in range(1,5):

    nshift = 'nshift-%d' % s
    Flow(nshift,'nft','window f2=%d | pad end2=%d' % (s,s))
    nshifts.append(nshift)

    nshift = 'nshift+%d' % s
    Flow(nshift,'nft','window n2=%d | pad beg2=%d ' % (nx-s,s))
    nshifts.append(nshift)


Flow('nshifts',nshifts,'cat ${SOURCES[1:%d]} axis=3 | put o2=0 ' % len(nshifts))


Flow('nflt npref','nshifts nft',
     'clpf match=${SOURCES[1]} pred=${TARGETS[1]} rect2=20 rect1=3 niter=10 verb=n')

Flow('npre','npref','fft1 inv=y ')

Result('npre','npre','window j2=5 | wiggle label2=Distance unit2=km transp=y yreverse=y poly=y  title="FX RNA"')

Flow('ndiff','npara npre','math x=${SOURCES[1]} output="x-input"')

Result('ndiff','ndiff','grey pclip=98 color=j maxval=0.1 minval=-0.1  scalebar=y title="FX RNA" wherexlabel=b  wheretitle=t ')

Result('npar','nflt','window n3=1 f3=0 | real | grey label2=Distance unit2=km color=j title="Coef" wherexlabel=b  wheretitle=t  ')




########## Patch

Flow('patch0','npara','patch w=501,20 p=1,50 | patch inv=y weight=y ')

Flow('patch','npara','patch w=501,20 p=1,50  ')
Flow('wpatch','patch','window')
fxds = []
mpas = []
for nw in range(0,50):
    data = 'data%d' % nw
    fxd  = 'fx%d'   % nw
    Flow(data,'wpatch','window n3=1 f3=%d' % nw)
    Flow(fxd,data,'fxdecon lenf=4 n2w=20')
    fxds.append(fxd)

    lom = 'lom%d' %nw
    lag = 'lag%d' %nw
    mpa = 'mpa%d' %nw
    Flow([lom, lag],data,'hpef niter=200 a=4,5 lag=${TARGETS[1]}')
    Flow(mpa,[data,lom],'helicon filt=${SOURCES[1]}')
    mpas.append(mpa)
Flow('fxpatch',fxds,'cat ${SOURCES[1:%d]} axis=3 | transp plane=34 | patch inv=y weight=y' % len(fxds))
Flow('mpapatch',mpas,'cat ${SOURCES[1:%d]} axis=3 | transp plane=34 | patch inv=y weight=y' % len(mpas))

Result('fxpatch','fxpatch','window j2=5 | wiggle label2=Distance unit2=km transp=y yreverse=y poly=y  title="FX Decon"')
Flow('fxdiff','npara fxpatch','math x=${SOURCES[1]} output="x-input"')
Result('fxdiff','fxdiff',
       '''
       grey color=j pclip=98 label2=Distance unit2=km maxval=0.07 minval=-0.07 scalebar=y title="FX Decon" wheretitle=t wherexlabel=b
       ''')

Result('mpapatch','mpapatch',
       '''
       grey pclip=98 label2=Distance unit2=km color=j 
       maxval=0.07 minval=-0.07 title="TX PEF" wherexlabel=b  wheretitle=t scalebar=y
       ''' )
Result('tpefpatch','npara mpapatch',
       '''
       math x=${SOURCES[1]} output="input-x" | 
       window j2=5 | 
       wiggle label2=Distance unit2=km transp=y yreverse=y poly=y  title="TX PEF"
       ''')







End()