from rsf.proj import *

nr = 10000 # number of realizations
ns = 100   # number of samples

Flow('rand',None,
     '''
     spike n1=%d n2=%d d2=1 d1=1 o1=1 o2=1|
     noise seed=2006 rep=y
     ''' % (nr,ns))

means = []
ameans = []
varis = []
avaris = []
for i in range(ns):
    rand = 'rand%d' % i
    Flow(rand,'rand','window n2=%d' % (i+1))

    s  = 's%d'    % i
    ss = 'ss%d'   % i
    s2 = 's2-%d'  % i
    sem = 'sem%d' % i

    Flow(s,rand,'stack norm=n')
    Flow(ss,s,'math output="input*input" ')
    Flow(s2,rand,'math output="input*input" | stack norm=n')
    Flow(sem,[ss,s2],'add ${SOURCES[1]} mode=d scale=1,%d' % (i+1))

    mean = 'mean%d' % i
    Flow(mean,sem,'stack axis=1 norm=n | scale dscale=%g' % (1.0/nr))

    vari = 'vari%d' % i
    Flow(vari,[mean,sem],
         '''
         spray axis=1 n=%d o=1 d=1 |
         add scale=1,-1 ${SOURCES[1]} |
         math output="input*input" | stack axis=1 norm=n |
         math output="sqrt(input/%d)"
         ''' % (nr,nr))

    means.append(mean)
    varis.append(vari)

    if i > 0:
        h  = 'h%d'    % i
        h2 = 'h2-%d'  % i
        sh = 'sh%d'   % i
        avo = 'avo%d' % i

        Flow(h,rand,'math output=x2 | stack norm=n')
        Flow(h2,rand,'math output=x2*x2 | stack norm=n')
        Flow(sh,rand,'math output=input*x2 | stack norm=n')
        Flow(avo,[s,s2,h,h2,sh],
             '''
             math
             s=${SOURCES[0]} s2=${SOURCES[1]}
             h=${SOURCES[2]} h2=${SOURCES[3]} sh=${SOURCES[4]}
             output="(2*s*h*sh - s*s*h2 - %d*sh*sh)/(s2*(h*h-%d*h2))"
             ''' % (i+1,i+1))


        amean = 'amean%d' % i
        Flow(amean,avo,'stack axis=1 norm=n | scale dscale=%g' % (1.0/nr))

        avari = 'avari%d' % i
        Flow(avari,[amean,avo],
             '''
             spray axis=1 n=%d o=1 d=1 |
             add scale=1,-1 ${SOURCES[1]} |
             math output="input*input" | stack axis=1 norm=n |
             math output="sqrt(input/%d)"
             ''' % (nr,nr))

        ameans.append(amean)
        avaris.append(avari)

Flow('mean',means,'cat axis=1 ${SOURCES[1:%d]}' % ns)
Plot('mean',
     '''
     graph symbol=o title="Semblance Expectation"
     label1=N label2= unit1= unit2=
     min2=-0.01 max2=1.01 symbolsz=5
     ''')

Plot('tmean','mean',
     '''
     math output=1/x1 |
     graph plotcol=5 wanttitle=n wantaxis=n
     min2=-0.01 max2=1.01
     ''')

Result('mean','mean tmean','Overlay')

Flow('amean',ameans,'cat axis=1 ${SOURCES[1:%d]} | put o1=2' % (ns-1))
Plot('amean',
     '''
     graph symbol=o title="AVO Semblance Expectation"
     label1=N label2= unit1= unit2=
     min2=-0.01 max2=1.01 symbolsz=5
     ''')

Plot('tamean','amean',
     '''
     math output="2/x1" |
     graph plotcol=5 wanttitle=n wantaxis=n
     min2=-0.01 max2=1.01
     ''')

Result('amean','amean tamean','Overlay')

Flow('vari',varis,'cat axis=1 ${SOURCES[1:%d]}' % ns)
Plot('vari',
     '''
     graph symbol=o title="Semblance Standard Deviation"
     label1=N label2= unit1= unit2=
     min2=-0.01 max2=0.37 symbolsz=5
     ''')

Plot('tvari','vari',
     '''
     math output="sqrt(2*(x1-1)/(x1+2))/x1" |
     graph plotcol=5 wanttitle=n wantaxis=n
     min2=-0.01 max2=0.37
     ''')

Result('vari','vari tvari','Overlay')

Flow('avari',avaris,'cat axis=1 ${SOURCES[1:%d]} | put o1=2' % (ns-1))
Plot('avari',
     '''
     graph symbol=o title="AVO Semblance Standard Deviation"
     label1=N label2= unit1= unit2=
     min2=-0.01 max2=0.37 symbolsz=5
     ''')

Plot('tavari','avari',
     '''
     math output="sqrt(4*(x1-2)/(x1+2))/x1" |
     graph plotcol=5 wanttitle=n wantaxis=n
     min2=-0.01 max2=0.37
     ''')

Result('avari','avari tavari','Overlay')


End()