From df79f4df6e727164e24d28d6fdd1a4c6471fde11 Mon Sep 17 00:00:00 2001 From: sfomel Date: Wed, 15 Jan 2025 16:39:11 -0600 Subject: [PATCH] cleanup --- book/gallery/marmousi/oway/SConstruct | 91 -- book/gallery/teapot/myangmig/SConstruct | 27 - book/jlu/apefsnsep/curvemod/SConstruct | 231 ----- .../apefsnsep/curvemod/curvelet_denoise2d.py | 109 --- book/jlu/apefsnsep/nmostack/SConstruct | 135 --- .../apefsnsep/nmostack/curvelet_denoise2d.py | 75 -- book/jlu/fxyspfdenoi/cmp/SConstruct | 356 -------- .../jlu/fxyspfdenoi/cmp/curvelet_denoise3d.py | 171 ---- book/jlu/fxyspfdenoi/real/SConstruct | 174 ---- book/rsf/tutorial2017/field/SConstruct | 125 --- book/rsf/tutorials/facies/SConstruct | 13 - book/rsf/tutorials/nn/SConstruct | 807 ------------------ book/xjtu/test/maxwell/SConstruct | 298 ------- 13 files changed, 2612 deletions(-) delete mode 100644 book/gallery/marmousi/oway/SConstruct delete mode 100644 book/gallery/teapot/myangmig/SConstruct delete mode 100755 book/jlu/apefsnsep/curvemod/SConstruct delete mode 100755 book/jlu/apefsnsep/curvemod/curvelet_denoise2d.py delete mode 100755 book/jlu/apefsnsep/nmostack/SConstruct delete mode 100755 book/jlu/apefsnsep/nmostack/curvelet_denoise2d.py delete mode 100755 book/jlu/fxyspfdenoi/cmp/SConstruct delete mode 100755 book/jlu/fxyspfdenoi/cmp/curvelet_denoise3d.py delete mode 100755 book/jlu/fxyspfdenoi/real/SConstruct delete mode 100644 book/rsf/tutorial2017/field/SConstruct delete mode 100644 book/rsf/tutorials/facies/SConstruct delete mode 100644 book/rsf/tutorials/nn/SConstruct delete mode 100644 book/xjtu/test/maxwell/SConstruct diff --git a/book/gallery/marmousi/oway/SConstruct b/book/gallery/marmousi/oway/SConstruct deleted file mode 100644 index 6c7a7edc0c..0000000000 --- a/book/gallery/marmousi/oway/SConstruct +++ /dev/null @@ -1,91 +0,0 @@ -from rsf.proj import * -from rsf.gallery import marmousi - -marmousi.getvel('vel') -Result('vel','grey min2=0 max2=9.2 min1=0 max1=3 scalebar=y color=j allpos=y title="Velocity Model" bias=1.5 barreverse=y') - -marmousi.get_zodata('exp') -Result('exp','grey title="Exploding Reflector Data" ') - -# From velocity to slowness -Flow('slowness','vel','transp | transp plane=23 | math output=1/input') - -# Fourier transform and transpose - -fmax = 50 # maximum frequency -nf = 206 # number of frequencies - -Flow('fft','exp', - 'fft1 | window max1=%g | transp plane=12 | transp plane=23' % fmax) - -# Extended split-step migration -Flow('mig','fft slowness', - ''' - zomig3 ompnth=1 mode=m --readwrite=y verb=y - nrmax=40 dtmax=0.001 slo=${SOURCES[1]} | - window | transp - ''',split=[3,nf,[0]],reduce='add') - -Result('mig', - ''' - grey title=Migration - label2=Depth unit2=km - label1=Distance unit1=km - ''') - - -nw=200 -jw=1 -ow=1 -nt=751 -dt=0.004 -ot=0 -nmx=2301 -dmx=0.004 -omx=0 - -marmousi.get_shots('shots') - -Result('shots','byte | grey3 flat=n frame1=250 frame2=80 frame3=100 title=Shots') -# Receiver shot gather wavefield FFT -Flow('rfft','shots', - ''' - fft1 | window squeeze=n n1=%d min1=%g j1=%d | - spray axis=3 n=1 o=0 d=1 label=hy | spray axis=5 n=1 o=0 d=1 label=sy - ''' % (nw,ow,jw), local=1) - -# Source wavefield FFT -Flow('sfft',None, - ''' - spike k1=1 n1=%d d1=%g | - fft1 | window squeeze=n n1=%d min1=%g j1=%d | put label1=w - ''' % (nt,dt,nw,ow,jw),local=1) - -# Interpolate wavefields on surface grid -Flow('rwav swav','rfft sfft','srsyn nx=%d dx=%g ox=%g wav=${SOURCES[1]} swf=${TARGETS[1]}' % (nmx,dmx,omx),local=1) -# w, x, y, s - -# Transpose and setting coordinates for 3-D migration -Flow('rtra','rwav','transp plane=12 | transp plane=23',local=1) -Flow('stra','swav','transp plane=12 | transp plane=23',local=1) -# x, y, w, s -# Prepare slowness on 3-D grid -Flow('slow','vel', - ''' - transp | window j1=1 | math output=1/input | - spray axis=2 n=1 d=1 o=0 | - window j3=1 squeeze=n - ''',local=1) - -Flow('img cig','stra rtra slow', - ''' - srmig3 nrmax=20 dtmax=5e-05 eps=0.01 --readwrite=y verb=y ompnth=1 - tmx=16 rwf=${SOURCES[1]} slo=${SOURCES[2]} cig=${TARGETS[1]} - itype=o - ''',split=[4,500,[0,1]],reduce='add') - -Result('image','img','grey title="PSPI wave equation image"',local=1) - -Result('agc','img','agc rect1=600 rect2=400 | grey title="PSPI wave equation image"',local=1) - -End() diff --git a/book/gallery/teapot/myangmig/SConstruct b/book/gallery/teapot/myangmig/SConstruct deleted file mode 100644 index da1c4a6346..0000000000 --- a/book/gallery/teapot/myangmig/SConstruct +++ /dev/null @@ -1,27 +0,0 @@ -from rsf.proj import * - -SConscript('../zomig/SConstruct') - -data3='../zomig/final_stack.rsf' -Flow ('data', data3,'put o2=0 d2=33.55 unit2=m o3=0 d3=33.55 unit3=m n4=1 o4=0 d4=100 unit4=m') - -Flow('vrms3d','../zomig/Vrms_interp.rsf', - ''' - math output="input*0.305" | - spray axis=2 n=189 d=33.55 o=0 - | spray axis=3 n=346 d=33.55 o=0 | put o2=0 d2=33.55 o3=0 d3=33.55 unit2=m unit3=m unit1=s | stack axis=4 | stack axis=4 | stack axis=4 - ''') - -# migration -Flow ('image', 'data vrms3d', - ''' - tmigda vel=${SOURCES[1]} - is3d=y isDipAz=n - iyo=0 iyn=1 iyd=33.55 - ixo=0 ixn=188 ixd=33.55 - dipn=81 dipo=-40 dipd=1 - sdipn=81 sdipo=-40 sdipd=1 - ''') -Result ('image','byte | grey3 wanttitle=n') - -End() diff --git a/book/jlu/apefsnsep/curvemod/SConstruct b/book/jlu/apefsnsep/curvemod/SConstruct deleted file mode 100755 index b0e8b1de7a..0000000000 --- a/book/jlu/apefsnsep/curvemod/SConstruct +++ /dev/null @@ -1,231 +0,0 @@ -from rsf.proj import * -import ssl - -ssl._create_default_https_context = ssl._create_unverified_context -### -# To finishing curvelet denoising part, it needs to install pylops and curvelops. -# pylops : pip install pylops -# curvelops : follow the instructions in the link: https://github.com/PyLops/curvelops -### - -### curve model -input = 'input.curve.segy' -Fetch(input,'yang') - -Flow("curve hcurve bcurve", - input, - """ - segyread tape=$SOURCE read=d hfile=${TARGETS[1]} bfile=${TARGETS[2]} - """, - stdin=0) -Flow("curve2", "curve", "window min1=1.2 max1=2 | bandpass fhi=60") - -# Nonstationary PEFs -Flow("cpad cmask", "curve2", "lpad jump=4 mask=${TARGETS[1]}") -Flow("cdmask", "cpad", "math output=1.") -Flow( - "capef", "cpad cdmask", """ - apef a=20,3 jump=4 rect1=20 rect2=3 niter=200 verb=y - maskin=${SOURCES[1]} - """) -Flow("data", "cpad capef cmask", - "miss4 filt=${SOURCES[1]} mask=${SOURCES[2]} verb=y") -Result( - "cm-data", "data", """ - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -# FK -Flow("fkdata", "data", 'fft1 | fft3 axis=2 | math output="abs(input)" | real') -Result( - "cm-fkdata", "fkdata", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz screenht=6. screenratio=0.5 - """) - -# Add noise -Flow( - "noiz", "data", """ - noise rep=n type=n seed=20106 range=0.0234818 | - smooth rect1=3 - """) -Result( - "cm-noiz", "noiz", """ - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -# FK -Flow("fknoiz", "noiz", "fft1 | fft3 axis=2 | math output='abs(input)' | real") -Result( - "cm-fknoiz", "fknoiz", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz screenht=6. screenratio=0.5 - """) -# snr -Flow("diff", "data noiz", "add scale=1,-1 ${SOURCES[1]}") -Flow("snr", "data diff", "snr2 noise=${SOURCES[1]}") - -############################################################################## -### PEF -# hpef -Flow("hnpef hnlag", "noiz", "hpef a=5,1 lag=${TARGETS[1]} ") -Flow("hspef hslag", "noiz", "hpef a=11,4 lag=${TARGETS[1]} ") -Flow( - "hsign", "noiz hspef hnpef", """ - signoi epsilon=1.0 sfilt=${SOURCES[1]} nfilt=${SOURCES[2]} - spitz=n niter=100 - """) -Result( - "cm-pef", "hsign", """ - window n3=1 f3=0 | - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -Result( - "cm-pefnoiz", "hsign", """ - window n3=1 f3=1 | - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -# FK -Flow( - "fkhsign", "hsign", - 'window n3=1 f3=0 | fft1 | fft3 axis=2 | math output="abs(input)" | real') -Result( - "cm-fkpef", "fkhsign", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz screenht=6. screenratio=0.5 - """) -# snr -Flow("hdiff", "data hsign", "add scale=1,-1 ${SOURCES[1]}") -Flow("hsnr", "data hdiff", "snr2 noise=${SOURCES[1]}") - -### APEF -Flow( - "anpef anpfpred", "noiz", """ - apef a=9,1 rect1=300 rect2=1 niter=200 - jump=1 pred=${TARGETS[1]} verb=y - """) -Flow( - "aspef aspefpred", "noiz", """ - apef a=11,4 rect1=30 rect2=15 niter=200 - jump=1 pred=${TARGETS[1]} verb=y - """) -Flow( - "asign", "noiz aspef anpef", """ - apefsignoi sfilt=${SOURCES[1]} nfilt=${SOURCES[2]} - niter=1000 eps=0.25 verb=y - """) -Flow("anoiz", "noiz asign", "add scale=1,-1 ${SOURCES[1]}") -Result( - "cm-apef", "asign", """ - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -Result( - "cm-apefnoiz", "anoiz", """ - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -# FK -Flow("fkasign", "asign", - " fft1 | fft3 axis=2 | math output='abs(input)' | real") -Result( - "cm-fkapef", "fkasign", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz screenht=6. screenratio=0.5 - """) -# snr -Flow("adiff", "data asign", "add scale=1,-1 ${SOURCES[1]}") -Flow("asnr", "data adiff", "snr2 noise=${SOURCES[1]}") - -### FX decon -Flow("patch", "noiz", "patch w=401,80") -Flow("wpatch", "patch", "window") -fxds = [] -mpas = [] -for nw in range(0, 5): - data = "data%d" % nw - fxd = "fx%d" % nw - Flow(data, "wpatch", "window n3=1 f3=%d" % nw) - Flow(fxd, data, "fxdecon lenf=8 n2w=10") - fxds.append(fxd) - -Flow( - "fxpatch", fxds, """ - cat ${SOURCES[1:%d]} axis=3 | transp plane=34 | - patch inv=y weight=y - """ % len(fxds)) -Flow("fxdif", "noiz fxpatch", "add scale=1,-1 ${SOURCES[1]}") - -Result( - "cm-fx", "fxpatch", """ - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -Result( - "cm-fxnoiz", "fxdif", """ - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) - -Flow("fkfxpatch", "fxpatch", - " fft1 | fft3 axis=2 | math output='abs(input)' | real") -Result( - "cm-fkfx", "fkfxpatch", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz screenht=6. screenratio=0.5 - """) -# snr -Flow("fxdiff", "data fxpatch", "add scale=1,-1 ${SOURCES[1]}") -Flow("fxsnr", "data fxdiff", "snr2 noise=${SOURCES[1]}") - -#! Need to install python package: pylops, curvelops -#! Follow "python curvelet_denoise.py" to generate "ct_denoise.rsf" -### curvelet -#Flow("ct_denoise", "ct_denoise.H", "dd form=native") - -import sys -python = sys.executable - -Command('ct_denoise.rsf','curvelet_denoise2d.py',python + ' $SOURCE') - -Result( - "cm-ct", "ct_denoise", """ - put d1=0.002 d2=0.25 o1=1.2 o2=0 label1=Time unit1=s label2=Trace | - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -Flow("ct_denoise_err", "noiz ct_denoise", "add scale=1,-1 ${SOURCES[1]}") -Result( - "cm-ctnoiz", "ct_denoise_err", """ - put d1=0.002 d2=0.25 o1=1.2 o2=0 label1=Time unit1=s label2=Trace | - grey yreverse=y transp=y poly=y label2=Position title= - screenratio=0.4 screenht=6. labelsz=5. titlesz=7 clip=0.02 - labelfat=2 font=2 titlefat=2 unit2=km - """) -# FK -Flow( - "fkct", "ct_denoise", """ - put d1=0.002 d2=0.25 o1=1.2 o2=0 | - fft1 | fft3 axis=2 | math output="abs(input)" | real - """) -Result( - "cm-fkct", "fkct", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz screenht=6. screenratio=0.5 - """) -# snr -Flow("ctdiff", "data ct_denoise", "add scale=1,-1 ${SOURCES[1]}") -Flow("ctsnr", "data ctdiff", "snr2 noise=${SOURCES[1]}") - -End() diff --git a/book/jlu/apefsnsep/curvemod/curvelet_denoise2d.py b/book/jlu/apefsnsep/curvemod/curvelet_denoise2d.py deleted file mode 100755 index 46a7b8057d..0000000000 --- a/book/jlu/apefsnsep/curvemod/curvelet_denoise2d.py +++ /dev/null @@ -1,109 +0,0 @@ -import numpy as np -import os -import pylops -from curvelops import FDCT2D -import m8r - - - -def byte2rsf(rsffile, inpfile, shape): - """ - convert numpy array to rsf file - """ - rsf = m8r.Output(rsffile) - inp = np.array(inpfile, dtype="float32") - rsf.put("n1", shape[1]) - rsf.put("n2", shape[0]) - rsf.put("o1", 0.0) - rsf.put("o2", 0.0) - - rsf.write(inp) - rsf.close() - - -# show data -os.system( - """ - sfgrey < noiz.rsf title= clip=0.02 | sfpen - """) - -# input data -inp = m8r.Input("noiz.rsf") -n1 = inp.int("n1") -n2 = inp.int("n2") -d = inp.read(shape=(n2, n1)) - -# print(d.shape) -nx, nt = d.shape -dx, dt = 0.04238, 0.004 -x, t = np.arange(nx) * dx, np.arange(nt) * dt - -Cop = FDCT2D(d.shape) -Wop = pylops.signalprocessing.DWT2D(d.shape, wavelet="db9", level=2) -# print(Cop.shape) - - -def reconstruct(data, op, perc=0.1, inv=False): - """ - Convenience function to calculate reconstruction using top - `perc` percent of coefficients of a given operator `op`. - """ - y = op * data.ravel() - denoise = np.zeros_like(y) - - # Order coefficients by strength - strong_idx = np.argsort(-np.abs(y)) - strong = np.abs(y)[strong_idx] - - # Select only top `perc`% coefficients - strong_idx = strong_idx[:int(np.rint(len(strong_idx) * perc))] - denoise[strong_idx] = y[strong_idx] - if inv: - data_denoise = op.inverse(denoise).reshape(data.shape) - else: - data_denoise = (op.H * denoise).reshape(data.shape) - return np.real(data_denoise), strong - - -#perc = 0.08 -cperc = 0.07 -d_dct, dct_strong = reconstruct(d, Cop, perc=cperc, inv=True) -# dct_err = d - d_dct -print(d_dct.shape) - - -byte2rsf("ct_denoise.rsf", d_dct, d_dct.shape) -os.system( - """ - sfput < ct_denoise.rsf d1=0.002 d2=0.25 o1=1.2 o2=0 | \ - sfgrey title= clip=0.02 font=2 titlefat=2 labelfat=3 gridfat=3 | sfpen - """) - -# byte2rsf("ct_denoise_err.rsf", dct_err, dct_err.shape) -# os.system( -# """ -# sfput < ct_denoise_err.rsf d1=0.002 d2=0.25 o1=1.2 o2=0 | \ -# sfgrey title= clip=0.02 font=2 titlefat=2 labelfat=3 gridfat=3 | sfpen -# """) - - -#perc = 0.08 -# wperc = 0.05 -# d_dwt, dwt_strong = reconstruct(d, Wop, perc=wperc) -# dwt_err = d - d_dwt -# print(d_dwt.shape) - - -# byte2rsf("wt_denoise.rsf", d_dwt, d_dwt.shape) -# os.system( -# """ -# sfput < wt_denoise.rsf d1=0.002 d2=0.25 o1=1.2 o2=0 | \ -# sfgrey title= clip=0.02 font=2 titlefat=2 labelfat=3 gridfat=3 | sfpen -# """) - -# byte2rsf("wt_denoise_err.rsf", dwt_err, dwt_err.shape) -# os.system( -# """ -# sfput < wt_denoise_err.rsf d1=0.002 d2=0.25 o1=1.2 o2=0 | \ -# sfgrey title= clip=0.02 font=2 titlefat=2 labelfat=3 gridfat=3 | sfpen -# """) diff --git a/book/jlu/apefsnsep/nmostack/SConstruct b/book/jlu/apefsnsep/nmostack/SConstruct deleted file mode 100755 index b56e0542bb..0000000000 --- a/book/jlu/apefsnsep/nmostack/SConstruct +++ /dev/null @@ -1,135 +0,0 @@ -from rsf.proj import * -import ssl - -ssl._create_default_https_context = ssl._create_unverified_context -### -# To finishing curvelet denoising part, it needs to install pylops and curvelops. -# pylops : pip install pylops -# curvelops : follow the instructions in the link: https://github.com/PyLops/curvelops -### - -### read data -sgy = "nmo_stack_agc.sgy" -Fetch(sgy,'yang') -Flow( - "oldstack toldstack oldstack.asc oldstack.bin", sgy, """ - segyread tfile=${TARGETS[1]} hfile=${TARGETS[2]} bfile=${TARGETS[3]} - """) -Result( - "stack1", "oldstack", """ - window n2=100 min1=2 max1=4 | put d2=0.04238 | - wiggle transp=y yreverse=y color=n wanttitle=n - wherexlabel=top poly=y parallel2=n n2tic=20 labelfat=3 gridfat=3 - label1=Time label2=Distance unit1=s unit2=km - """) -Flow("data", "oldstack", "window n2=100 min1=2.5 max1=5 | put d2=0.04238") -Result("nmo-data", "data", - "grey title= font=2 titlefat=2 labelfat=3 gridfat=3") -# FK -Flow("fkdata", "data", 'fft1 | fft3 axis=2 | math output="abs(input)" | real') -Result( - "nmo-fkdata", "fkdata", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz - """) - -# bandpass or not -Flow("noise", "data", "bandpass flo=30") -Result("noise", - "grey title=noise clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") - -# hpef -Flow("hnpef hnlag", "data", "hpef a=9,1 lag=${TARGETS[1]} ") -Flow("hspef hslag", "data", "hpef a=11,4 lag=${TARGETS[1]} ") -Flow( - "hsign", "data hspef hnpef", """ - signoi epsilon=1.0 sfilt=${SOURCES[1]} nfilt=${SOURCES[2]} - spitz=n niter=100 - """) -Flow("hss", "hsign", "window n3=1 f3=0") -Flow("hnn", "hsign", "window n3=1 f3=1") -Result("nmo-hpef", "hss", - "grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") -Result("nmo-hpefnoiz", "hnn", - "grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") -# FK -Flow("fkhpef", "hss", "fft1 | fft3 axis=2 | math output='abs(input)' | real") -Result( - "nmo-fkhpef", "fkhpef", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 label1=Frequency - label2=Wavenumber unit1=Hz - """) - -### apef -Flow( - "napef npred", "noise", """ - apef a=9,1 rect1=60 rect2=1 niter=500 - jump=1 pred=${TARGETS[1]} verb=y - """) -Flow( - "sapef spred", "data", """ - apef a=9,4 rect1=60 rect2=20 niter=500 - jump=1 pred=${TARGETS[1]} verb=y - """) -Flow( - "ass", "data sapef napef", """ - apefsignoi sfilt=${SOURCES[1]} nfilt=${SOURCES[2]} - niter=800 eps=15 verb=y - """) -Flow("ann", "data ass", "add scale=1,-1 ${SOURCES[1]}") -Result("nmo-apef", "ass", - "grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") -Result("nmo-apefnoiz", "ann", - "grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") -# FK -Flow("fkapef", "ass", "fft1 | fft3 axis=2 | math output='abs(input)' | real") -Result( - "nmo-fkapef", "fkapef", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz - """) - -### fxdecon -Flow("fxsign", "data", "fxdecon lenf=8 n2w=30") -Flow("fxnoiz", "data fxsign", "add scale=1,-1 ${SOURCES[1]}") -Result("nmo-fx", "fxsign", - "grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") -Result("nmo-fxnoiz", "fxnoiz", - "grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3") -# FK -Flow("fkfxdecon", "fxsign", - "fft1 | fft3 axis=2 | math output='abs(input)' | real") -Result( - "nmo-fkfx", "fkfxdecon", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz - """) - -#! Need to install python package: pylops, curvelops -#! Follow "python curvelet_denoise.py" to generate "ct_denoise.rsf" -### curvelet -#Flow("ct_denoise", "ct_denoise.H", "dd form=native") -Result( - "nmo-ct", "ct_denoise", """ - put d1=0.004 d2=0.04238 o1=2.5 o2=0 label1=Time unit1=s label2=Trace | - grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3 - """) -Flow("ct_denoise_err", "data ct_denoise", "add scale=1,-1 ${SOURCES[1]}") -Result( - "nmo-ctnoiz", "ct_denoise_err", """ - put d1=0.004 d2=0.04238 o1=2.5 o2=0 label1=Time unit1=s label2=Trace | - grey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3 - """) -# FK -Flow( - "fkct", "ct_denoise", """ - put d1=0.004 d2=0.04238 o1=2.5 o2=0| - fft1 | fft3 axis=2 | math output="abs(input)" | real - """) -Result( - "nmo-fkct", "fkct", """ - grey scalebar=n title= unit2= color=i max1=100 labelfat=2 font=2 - label1=Frequency label2=Wavenumber unit1=Hz - """) - -End() diff --git a/book/jlu/apefsnsep/nmostack/curvelet_denoise2d.py b/book/jlu/apefsnsep/nmostack/curvelet_denoise2d.py deleted file mode 100755 index 5d5138b569..0000000000 --- a/book/jlu/apefsnsep/nmostack/curvelet_denoise2d.py +++ /dev/null @@ -1,75 +0,0 @@ -import numpy as np -import os -import pylops -from curvelops import FDCT2D -import m8r - - -def byte2rsf(rsffile, inpfile, shape): - """ - convert numpy array to rsf file - """ - rsf = m8r.Output(rsffile) - inp = np.array(inpfile, dtype="float32") - rsf.put("n1", shape[1]) - rsf.put("n2", shape[0]) - rsf.put("o1", 0.0) - rsf.put("o2", 0.0) - - rsf.write(inp) - rsf.close() - - -# show data -os.system(""" - sfgrey < data.rsf title= | sfpen - """) - -inp = m8r.Input("data.rsf") -n1 = inp.int("n1") -n2 = inp.int("n2") -d = inp.read(shape=(n2, n1)) -# print(d.shape) - -nx, nt = d.shape -dx, dt = 0.04238, 0.004 -x, t = np.arange(nx) * dx, np.arange(nt) * dt - -Cop = FDCT2D(d.shape) -Wop = pylops.signalprocessing.DWT2D(d.shape, wavelet="db9", level=2) -# print(Cop.shape) - - -def reconstruct(data, op, perc=0.1, inv=False): - """ - Convenience function to calculate reconstruction using top - `perc` percent of coefficients of a given operator `op`. - """ - y = op * data.ravel() - denoise = np.zeros_like(y) - - # Order coefficients by strength - strong_idx = np.argsort(-np.abs(y)) - strong = np.abs(y)[strong_idx] - - # Select only top `perc`% coefficients - strong_idx = strong_idx[:int(np.rint(len(strong_idx) * perc))] - denoise[strong_idx] = y[strong_idx] - if inv: - data_denoise = op.inverse(denoise).reshape(data.shape) - else: - data_denoise = (op.H * denoise).reshape(data.shape) - return np.real(data_denoise), strong - - -#perc = 0.08 -cperc = 0.06 -d_dct, dct_strong = reconstruct(d, Cop, perc=cperc, inv=True) -# dct_err = d - d_dct -# print(d_dct.shape) - -byte2rsf("ct_denoise.rsf", d_dct, d_dct.shape) -os.system(""" - sfput < ct_denoise.rsf d1=0.008 d2=0.04238 o1=0 o2=0 | \ - sfgrey title= clip=3.10373 font=2 titlefat=2 labelfat=3 gridfat=3 | sfpen - """) \ No newline at end of file diff --git a/book/jlu/fxyspfdenoi/cmp/SConstruct b/book/jlu/fxyspfdenoi/cmp/SConstruct deleted file mode 100755 index dad82951ed..0000000000 --- a/book/jlu/fxyspfdenoi/cmp/SConstruct +++ /dev/null @@ -1,356 +0,0 @@ -import math -from rsf.proj import * - -### -# To finishing curvelet denoising part, need to install pylops and curvelops. -# pylops : pip install pylops -# curvelops : follow the instructions: https://github.com/PyLops/curvelops -### - -# Plot font and screen ratio, width, height for uniform ppt figs. -p1 = 0.7 -sr = 1.0 -sw = 7.0 -sh = 10.0 -xll = 2.0 -fat = 2 - -# Synthetic CMP Parameters -nx = 101 -delx = 0.0125 * 2 -ox = -1.25 -ny = nx -dely = delx -oy = ox -nt = 1001 -dt = 0.004 -fw = 10.0 - -# 4 Reflection Events -# 111111111 Make a layer: z = z0 + x*dx + y*dy - -dx = 0.4 -dy = 0.1 -dt = 0.004 -z0 = 0.8 -# z0=0.75 - -cg = 1 / math.sqrt(1 + dx * dx + dy * dy) -ca = -dx * cg -cb = -dy * cg -d = z0 * cg - -mx = 0 -my = 0 - -D = d - mx * ca - my * cb - -v0 = 2.5 -t0 = 2 * D / v0 -it0 = int(t0 / dt) -it1 = it0 -# print it1 -Flow("spike1", None, - """ - spike n1=%d k1=%d nsp=1 | ricker1 frequency=%g | - spray axis=2 n=%d o=%g d=%g | - spray axis=3 n=%d o=%g d=%g - """ % (nt, it1, fw, ny, oy, dely, nx, ox, delx)) - -wx = 1 - ca * ca -wy = 1 - cb * cb -wxy = -ca * cb - -Flow("vel1.asc", None, - "echo %g %g %g n1=3 data_format=ascii_float in=$TARGET" % (wx, wy, wxy)) -Flow("vel1","vel1.asc", - """ - dd form=native | scale dscale=%g | spray axis=1 n=%d - """ % (1 / (v0 * v0), nt), - local=1) -Flow("cmp1", "spike1 vel1", "inmo3 velocity=${SOURCES[1]}", local=1) - -# 222222222222 Make another layer: z = z0 + x*dx + y*dy -dx = 0.7 -dy = 0.41 -dt = 0.004 -z0 = 0.85 - -cg = 1 / math.sqrt(1 + dx * dx + dy * dy) -ca = -dx * cg -cb = -dy * cg -d = z0 * cg - -mx = 0 -my = 0 - -D = d - mx * ca - my * cb - -v0 = 1.7 -t0 = 2 * D / v0 -it0 = int(t0 / dt) -it2 = it0 -# print it2 -Flow("spike2", None, - """ - spike n1=%d k1=%d nsp=1 | ricker1 frequency=%g | - spray axis=2 n=%d o=%g d=%g | - spray axis=3 n=%d o=%g d=%g - """ % (nt, it2, fw, ny, oy, dely, nx, ox, delx)) - -wx = 1 - ca * ca -wy = 1 - cb * cb -wxy = -ca * cb - -Flow("vel2.asc", None, - "echo %g %g %g n1=3 data_format=ascii_float in=$TARGET" % (wx, wy, wxy)) -Flow("vel2","vel2.asc", - ''' - dd form=native | scale dscale=%g | - spray axis=1 n=%d - ''' % (1 / (v0 * v0), nt), local=1) -Flow("cmp2", "spike2 vel2", "inmo3 velocity=${SOURCES[1]}", local=1) - -# 3333333333333 Make yet another layer: z = z0 + x*dx + y*dy -dx = 0.1 -dy = 0.9 -dt = 0.004 -z0 = 1.1 - -cg = 1 / math.sqrt(1 + dx * dx + dy * dy) -ca = -dx * cg -cb = -dy * cg -d = z0 * cg - -mx = 0 -my = 0 - -D = d - mx * ca - my * cb - -v0 = 1.75 -t0 = 2 * D / v0 -it0 = int(t0 / dt) -it3 = it0 -# print it3 -Flow("spike3", None, - """ - spike n1=%d k1=%d nsp=1 | ricker1 frequency=%g | - spray axis=2 n=%d o=%g d=%g | - spray axis=3 n=%d o=%g d=%g - """ % (nt, it3, fw, ny, oy, dely, nx, ox, delx)) - -wx = 1 - ca * ca -wy = 1 - cb * cb -wxy = -ca * cb - -Flow("vel3.asc", None, - "echo %g %g %g n1=3 data_format=ascii_float in=$TARGET" % (wx, wy, wxy)) -Flow("vel3","vel3.asc", - """ - dd form=native | scale dscale=%g | spray axis=1 n=%d - """ % (1 / (v0 * v0), nt), local=1) -Flow("cmp3", "spike3 vel3", "inmo3 velocity=${SOURCES[1]}", local=1) - -# 444444444444 Make a 4th layer: z = z0 + x*dx + y*dy -dx = 0.2 -dy = 0.1 -dt = 0.004 -z0 = 1.3 - -cg = 1 / math.sqrt(1 + dx * dx + dy * dy) -ca = -dx * cg -cb = -dy * cg -d = z0 * cg - -mx = 0 -my = 0 -D = d - mx * ca - my * cb - -v0 = 2.0 -t0 = 2 * D / v0 -it0 = int(t0 / dt) -it4 = it0 -# print it4 -Flow("spike4", None, - """ - spike n1=%d k1=%d nsp=1 | ricker1 frequency=%g | - spray axis=2 n=%d o=%g d=%g | - spray axis=3 n=%d o=%g d=%g - """ % (nt, it4, fw, ny, oy, dely, nx, ox, delx)) - -wx = 1 - ca * ca -wy = 1 - cb * cb -wxy = -ca * cb - -Flow("vel4.asc", None, - "echo %g %g %g n1=3 data_format=ascii_float in=$TARGET" % (wx, wy, wxy)) -Flow("vel4","vel4.asc", - """ - dd form=native | scale dscale=%g | spray axis=1 n=%d - """ % (1 / (v0 * v0), nt), local=1) -Flow("cmp4", "spike4 vel4", "inmo3 velocity=${SOURCES[1]}", local=1) - -# Add events to create CMP. -Flow("cmp", "cmp1 cmp2 cmp3 cmp4", - """ - math cmp2=${SOURCES[1]} cmp3=${SOURCES[2]} cmp4=${SOURCES[3]} - output="input+cmp2+cmp3+cmp4" | - window max1=1.5 min1=0.5 j1=2 | put o1=0 - """, local=1) -Result("cmpmod", "cmp", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("noise", "cmp", "noise seed=1 range=0.3 ") -Result("cmpnoise", "noise", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("sig0", "cmp", 'math output="input*input" | stack axis=0') -Flow("noi0", "cmp noise", - """ - math A=${SOURCES[0]} B=${SOURCES[1]} output="(A-B)*(A-B)" | - stack axis=0 - """) -Flow("snr0", "sig0 noi0", - "math A=${SOURCES[0]} B=${SOURCES[1]} output='10*log(A/B)/log(10)'") - -### f-x SPF -Flow("spf2", "noise", - """ - fxspfdenoise2 lambda1=4.5 lambda3=1 - na1=4 na2=0 verb=y - """) -Result("cmpspf2", "spf2", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("errspf2", "noise spf2", "add ${SOURCES[1]} scale=1,-1 ") -Result("cmperrspf2", "errspf2", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("noi1", "cmp spf2", - """ - math A=${SOURCES[0]} B=${SOURCES[1]} output="(A-B)*(A-B)" | - stack axis=0 - """) -Flow("snr1", "sig0 noi1", - "math A=${SOURCES[0]} B=${SOURCES[1]} output='10*log(A/B)/log(10)'") - -### f-x-y SPF -# padding data -length_pad = 50 -Flow("noise_pad", "noise", - """ - window n3=%d | reverse which=4 | put o3=0 d3=0.025 - """ % length_pad) -Flow("noises", "noise_pad noise", "cat ${SOURCES[1:2]} axis=3") -Flow("spf3", "noises", - """ - fxyspfdenoise3 lambda1=4.5 lambda2=4.5 lambda3=1 na1=2 na2=2 verb=y | - window f3=%d | put o3=-1.25 - """ % length_pad) -Result("cmpspf3", "spf3", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("errspf3", "noise spf3", "add ${SOURCES[1]} scale=1,-1 ") -Result("cmperrspf3", "errspf3", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("noi2", "cmp spf3", - """ - math A=${SOURCES[0]} B=${SOURCES[1]} output="(A-B)*(A-B)" | - stack axis=0 - """) -Flow("snr2", "sig0 noi2", - "math A=${SOURCES[0]} B=${SOURCES[1]} output='10*log(A/B)/log(10)'") - -### f-x-y-RNA -Flow("fnoise", "noise", "fft1 | transp plane=13 | transp plane=12") -Flow("shifts", "fnoise", "cshifts2 ns1=2 ns2=2 | transp plane=34") -Flow("flt pre", "shifts fnoise", - """ - clpf match=${SOURCES[1]} pred=${TARGETS[1]} rect1=7 rect2=7 niter=30 - """) -Flow("rna", "pre", "transp plane=13 | transp plane=23 | fft1 inv=y") -Result("cmprna", "rna", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("errrna", "noise rna", "add ${SOURCES[1]} scale=1,-1 ") -Result("cmperrrna", "errrna", - """ - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("noi3", "cmp rna", - """ - math A=${SOURCES[0]} B=${SOURCES[1]} output="(A-B)*(A-B)" | - stack axis=0 - """) -Flow("snr3", "sig0 noi3", - "math A=${SOURCES[0]} B=${SOURCES[1]} output='10*log(A/B)/log(10)'") - -#! Need to install python package: pylops, curvelops -#! Follow "python curvelet_denoise.py" to generate "ct_denoise.rsf" - -import sys -python = sys.executable - -Command('ct_denoise.rsf','curvelet_denoise3d.py',python + ' $SOURCE') - -### curvelet -Result("cmpct", "ct_denoise", - """ - put d1=0.008 d2=0.025 d3=0.025 01=0 o2=-1.25 o3=-1.25 - label1="Time" unit1="s" | - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("ct_denoise_err", "noise ct_denoise", "add ${SOURCES[1]} scale=1,-1") -Result("cmperrct", "ct_denoise_err", - """ - put d1=0.008 d2=0.025 d3=0.025 01=0 o2=-1.25 o3=-1.25 - label1="Time" unit1="s" | - byte gainpanel=all clip=0.1 | - grey3 flat=n frame1=80 frame2=50 frame3=50 - point1=0.6 point2=0.7 title= - label2="X" label3="Y" unit2="km" unit3="km" - """) -Flow("noi4", "cmp ct_denoise", - """ - math A=${SOURCES[0]} B=${SOURCES[1]} output="(A-B)*(A-B)" | - stack axis=0 - """) -Flow("snr4", "sig0 noi4", - "math A=${SOURCES[0]} B=${SOURCES[1]} output='10*log(A/B)/log(10)'") - -End() diff --git a/book/jlu/fxyspfdenoi/cmp/curvelet_denoise3d.py b/book/jlu/fxyspfdenoi/cmp/curvelet_denoise3d.py deleted file mode 100755 index 4a828ea422..0000000000 --- a/book/jlu/fxyspfdenoi/cmp/curvelet_denoise3d.py +++ /dev/null @@ -1,171 +0,0 @@ -import numpy as np -import os -import pylops -from curvelops import FDCT3D -import m8r -# import pyvista as pv - - -def byte2rsf(rsffile, inpfile, shape): - """ - convert numpy array to rsf file - """ - rsf = m8r.Output(rsffile) - inp = np.array(inpfile, dtype="float32") - rsf.put("n1", shape[2]) - rsf.put("n2", shape[1]) - rsf.put("n3", shape[0]) - rsf.put("o1", 0.0) - rsf.put("o2", 0.0) - rsf.put("o3", 0.0) - - rsf.write(inp) - rsf.close() - - -# show data -os.system( - """ - sfbyte < noise.rsf gainpanel=all clip=0.1 | \ - sfgrey3 flat=n frame1=80 frame2=50 frame3=50 \ - point1=0.6 point2=0.7 title= label2="X" label3="Y" unit2="km" unit3="km" | \ - sfpen - """) - -# input data -# inputfile = '../testdata/sigmoid.npz' -inp = m8r.Input("noise.rsf") -n1 = inp.int("n1") -n2 = inp.int("n2") -n3 = inp.int("n3") -d = inp.read(shape=(n3, n2, n1)) - -print(d.shape) -# d = np.load(inputfile) -# d = d['sigmoid'] -ny, nx, nt = d.shape -dy, dx, dt = 8, 8, 0.004 -y, x, t = np.arange(ny) * dy, np.arange(nx) * dx, np.arange(nt) * dt - -Cop = FDCT3D(d.shape) -print(Cop.shape) - - -def reconstruct(data, op, perc=0.1): - """ - Convenience function to calculate reconstruction using top - `perc` percent of coefficients of a given operator `op`. - """ - y = op * data.ravel() - denoise = np.zeros_like(y) - - # Order coefficients by strength - strong_idx = np.argsort(-np.abs(y)) - strong = np.abs(y)[strong_idx] - - # Select only top `perc`% coefficients - strong_idx = strong_idx[:int(np.rint(len(strong_idx) * perc))] - denoise[strong_idx] = y[strong_idx] - - data_denoise = op.inverse(denoise).reshape(data.shape) - return np.real(data_denoise), strong - - -#perc = 0.08 -perc = 0.1 -d_dct, dct_strong = reconstruct(d, Cop, perc=perc) -dct_err = d - d_dct -print(d_dct.shape) - - -byte2rsf("ct_denoise.rsf", d_dct, d_dct.shape) -os.system( - """ - sfput < ct_denoise.rsf d1=0.008 d2=0.025 d3=0.025 01=0 o2=-1.25 o3=-1.25 | \ - sfbyte gainpanel=all clip=0.1 | \ - sfgrey3 flat=n frame1=80 frame2=50 frame3=50 \ - point1=0.6 point2=0.7 title= label2="X" label3="Y" unit2="km" unit3="km" | \ - sfpen - """) - -byte2rsf("ct_denoise_err.rsf", dct_err, dct_err.shape) -os.system( - """ - sfput < ct_denoise_err.rsf d1=0.008 d2=0.025 \ - d3=0.025 01=0 o2=-1.25 o3=-1.25 | \ - sfbyte gainpanel=all clip=0.1 | \ - sfgrey3 flat=n frame1=80 frame2=50 frame3=50 \ - point1=0.6 point2=0.7 title= label2="X" label3="Y" unit2="km" unit3="km" | \ - sfpen - """) - - -# pyvista -# grid = pv.UniformGrid() -# grid.dimensions = np.array(d_dct.shape) + 1 -# # Edit the spatial reference -# grid.origin = (0.0, 0.0, 0.0) # The bottom left corner of the data set -# grid.spacing = (1, 1, 1) # These are the cell sizes along each axis - -# # Add the data values to the cell data -# grid.cell_arrays["values"] = d_dct.flatten(order="F") # Flatten the array! -# p = pv.Plotter() -# # boxes -# bounds = [0, 50, 50, 101, 0, 126] -# clipped = grid.clip_box(bounds) -# # noninteractive -# # Set a custom position and size -# sargs = dict( -# height=0.25, -# vertical=True, -# position_x=0.05, -# position_y=0.05, -# title_font_size=20, -# label_font_size=16, -# shadow=True, -# n_labels=10, -# italic=True, -# fmt="%.1f", -# font_family="arial", -# color="black" -# ) - -# p = pv.Plotter() -# # p.add_mesh(grid, scalar_bar_args=sargs, cmap="seismic") -# p.add_mesh( -# # grid, -# clipped, -# # color="000066", -# # style="wireframe", -# style="surface", -# # show_scalar_bar=False, -# scalar_bar_args=sargs, -# stitle="Amplitude", -# cmap="gray", -# # cmap="seismic", -# # cmap="coolwarm", -# lighting=False, -# # show_edges=True, -# ) -# p.set_background(color="white") -# p.show_bounds( -# xlabel="Ytrace", -# ylabel="Xtrace", -# zlabel="Time", -# grid="front", -# location="outer", -# all_edges=True, -# color="black" -# ) -# cpos = [ -# (-296.8959883360211, 375.2608528067814, -159.3368341522591), -# (50.0, 75.0, 100.0), -# (0.3432324736098539, -0.3546432150194057, -0.8697238982000901), -# ] -# # p.view_vector((100, 150, 0), viewup=True) -# p.camera_position = cpos - -# cpos = p.show("cmpnoise") -# # cpos = p.show("Qdome", screenshot="qdome.png") - -# print(cpos) diff --git a/book/jlu/fxyspfdenoi/real/SConstruct b/book/jlu/fxyspfdenoi/real/SConstruct deleted file mode 100755 index 48511b47a7..0000000000 --- a/book/jlu/fxyspfdenoi/real/SConstruct +++ /dev/null @@ -1,174 +0,0 @@ -"""denoise test in 3d model""" -from rsf.prog import RSFROOT -from rsf.proj import * -from rsf.recipes.beg import server - -Fetch('image3d.rsf','cup',server) -### -# To finishing fx-emdpf denoising part, it needs to install Matlab, -# and rebuild Madagascar with "API=Matlab" option. -### - -### load data -Flow("data", "image3d", "dd form=native") -Result("realdata", "data", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) - -### f-x SPF -Flow("spf2", "data", - """ - fxspfdenoise2 lambda1=280 lambda3=55 na1=4 na2=0 verb=y - """) -Result("realspf2", "spf2", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) -Flow("errspf2", "data spf2", "add ${SOURCES[1]} scale=1,-1 ") -Result("realerrspf2", "errspf2", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) - -### f-x-y SPF -# padding data -length_pad = 100 -Flow("data_pad", "data", - """ - window n3=%d | reverse which=4 | put o3=0 d3=0.03 - """ % length_pad) -Flow("datas", "data_pad data", "cat ${SOURCES[1:2]} axis=3") -# denoise -Flow("spf3", "datas", - """ - fxyspfdenoise3 lambda1=600 lambda2=600 lambda3=90 na1=2 na2=2 verb=y | - window f3=%d | put o3=0 - """ % length_pad) -Result("realspf3", "spf3", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) -Flow("errspf3", "data spf3", "add ${SOURCES[1]} scale=1,-1 ") -Result("realerrspf3", "errspf3", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) - -### f-x EMD (emdpf) -matlab = WhereIs("matlab") -matROOT = "../Matfun" -matfun = "Real" -matlabpath = os.environ.get("MATLABPATH", os.path.join(RSFROOT, "lib")) - -if not matlab: - sys.stderr.write("\nCannot find Matlab.\n") - sys.exit(1) - -fxemds = [] -for n in range(0, 310): - fxdat = "fxdat-%s" % n - fxemd = "fxemd-%s" % n - Flow(fxdat, "data", "window n3=1 f3=%d" % n) - Flow(fxemd, [fxdat, os.path.join(matROOT, matfun + ".m")], - """ - MATLABPATH=%(matlabpath)s %(matlab)s -nosplash -nojvm -r - "addpath %(matROOT)s;%(matfun)s('${SOURCES[0]}','${TARGETS[0]}');quit" - """ % vars(), - stdin=0, - stdout=-1) - fxemds.append(fxemd) - -Flow("fxemd", fxemds, "cat ${SOURCES[1:%d]} axis=3 " % len(fxemds)) -Flow("errfxemd", "data fxemd", "add ${SOURCES[1]} scale=1,-1 ") -Result("realfxemd", "fxemd", - """ - put o1=0 o2=0 o3=0 d1=0.001 d2=0.03 d3=0.03 | - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) -Result("realerrfxemd", "errfxemd", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) - -### f-x-y RNA (consume large memory space) -# Patch f=73*2/5 (patching for reducing computational burden) -Flow('patch','data','patch w=140,266,310 p=10,1,1') -tpres = [] -tpre2ds = [] - -for nwt in range(0,10): - fd = 'fd-%d' % nwt - shiftsa= 'shiftsa-%d' % nwt - sh1 = 'sh1-%d' % nwt - shifts = 'shifts-%d' % nwt - flt = 'flt-%d' % nwt - pre = 'pre-%d' % nwt - tpre = 'tpre-%d' % nwt - Flow(fd,'patch', - ''' - window n4=1 f4=%d | fft1 | transp plane=13 memsize=1000 | - transp plane=12 memsize=1000 - ''' % nwt ) - - Flow(shifts,fd, - ''' - cshifts2 ns1=2 ns2=2 | transp plane=34 memsize=1000 - ''' ) - - Flow([flt, pre],[shifts, fd], - ''' - clpf match=${SOURCES[1]} pred=${TARGETS[1]} rect1=10 rect2=10 niter=10 - ''') - Flow(tpre,pre, - ''' - transp plane=13 memsize=1000 | transp plane=23 memsize=1000 | - fft1 inv=y - ''') - tpres.append(tpre) -Flow('rna',tpres, - 'cat ${SOURCES[1:%d]} axis=4 | patch inv=y weight=y dim=3' % len(tpres)) - -Result("realrna", "rna", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) -Flow("errrna", "data rna", "add ${SOURCES[1]} scale=1,-1 ") -Result("realerrrna", "errrna", - """ - byte gainpanel=e clip=1.5 | - grey3 flat=y frame1=250 frame2=47 frame3=254 - point1=0.7 point2=0.45 title= color=g - label2=X label3=Y unit2=km unit3=km - """) - - - - - - -End() diff --git a/book/rsf/tutorial2017/field/SConstruct b/book/rsf/tutorial2017/field/SConstruct deleted file mode 100644 index 296c922aec..0000000000 --- a/book/rsf/tutorial2017/field/SConstruct +++ /dev/null @@ -1,125 +0,0 @@ -from rsf.proj import * - -# Download one CMP gather from Viking Graben -Fetch('cmp.HH','viking') - -Flow('cmp','cmp.HH','dd form=native') -Plot('cmp','grey title="CMP Gather" clip=8.66') - -Flow('rad','cmp','hypradon adj=y np=101 dp=0.01 op=0') -Plot('rad', - ''' - grey title="Hyperbolic Radon Transform" - label2=Slowness unit2=s/km - ''') - -Flow('adj','rad', - 'hypradon adj=n nx=60 dx=0.05 ox=0.287') -Plot('adj','grey title="Adjoint Hyperbolic Radon Transform" ') - -Result('adj','cmp rad adj','SideBySideAniso') - -Flow('inv','cmp rad', - ''' - conjgrad hypradon mod=${SOURCES[1]} - np=101 dp=0.01 op=0 nx=60 dx=0.05 ox=0.287 niter=10 - ''') -Plot('inv', - ''' - grey title="Hyperbolic Radon Transform" - label2=Slowness unit2=s/km - ''') - -Flow('cmp2','inv', - 'hypradon adj=n nx=60 dx=0.05 ox=0.287') -Plot('cmp2', - 'grey title="Inverse Hyperbolic Radon Transform" clip=8.66') - -Result('inv','cmp inv cmp2','SideBySideAniso') - -# Velocity scan - -Flow('vscan','cmp', - 'vscan semblance=y half=n v0=1.2 nv=131 dv=0.02') -Plot('vscan', - 'grey color=j allpos=y title="Semblance Scan" unit2=km/s') - -# Automatic pick - -Flow('vpick','vscan', - ''' - mutter inner=y x0=1.4 half=n v0=0.45 t0=0.5 | - pick rect1=50 vel0=1.5 - ''') -Plot('vpick', - ''' - graph yreverse=y transp=y plotcol=7 plotfat=7 - pad=n min2=1.2 max2=3.8 wantaxis=n wanttitle=n - ''') - -Plot('vscan2','vscan vpick','Overlay') - -# NMO - -Flow('nmo','cmp vpick','nmo half=n velocity=${SOURCES[1]}') -Plot('nmo','grey title="NMO with Primary Velocity" clip=9.83') - -Result('vscan2','cmp vscan2 nmo','SideBySideAniso') - -Flow('nmorad0','cmp', - 'hypradon adj=y np=161 dp=0.005 op=0') -Flow('nmorad','nmo nmorad0', - ''' - conjgrad hypradon mod=${SOURCES[2]} - np=161 dp=0.005 op=0 nx=60 dx=0.05 ox=0.287 niter=10 - ''') -Plot('nmorad', - ''' - grey title="Hyperbolic Radon Transform" - label2=Slowness unit2=s/km - ''') - -Flow('nmo2','nmorad', - 'hypradon adj=n nx=60 dx=0.05 ox=0.287') -Plot('nmo2', - 'grey title="Inverse Hyperbolic Radon Transform" clip=9.83') - -Result('nmorad','nmo nmorad nmo2','SideBySideAniso') - -Flow('cut','nmorad','cut min2=0.2') -Plot('cut', - ''' - grey title="Hyperbolic Radon Transform" - label2=Slowness unit2=s/km - ''') - -Flow('signal','cut', - 'hypradon adj=n nx=60 dx=0.05 ox=0.287') -Plot('signal', - 'grey title="Inverse Hyperbolic Radon Transform" clip=9.83') - -Result('nmocut','nmo cut signal','SideBySideAniso') - -# Inverse NMO - -Flow('inmo','signal vpick', - 'inmo half=n velocity=${SOURCES[1]} | mutter half=n v0=1.2') -Plot('inmo','grey title=Signal clip=8.66') - -Flow('pvscan','inmo', - 'vscan semblance=y half=n v0=1.2 nv=131 dv=0.02') -Plot('pvscan', - 'grey color=j allpos=y title="Primary Semblance Scan" ') - -Flow('mult','cmp inmo','add scale=1,-1 ${SOURCES[1]}') -Plot('mult','grey title=Noise clip=8.66') - -Flow('mvscan','mult', - 'vscan semblance=y half=n v0=1.2 nv=131 dv=0.02') -Plot('mvscan', - 'grey color=j allpos=y title="Multiples Semblance Scan" ') - -Result('inmo','cmp inmo mult','SideBySideAniso') -Result('pvscan','vscan pvscan mvscan','SideBySideAniso') - -End() diff --git a/book/rsf/tutorials/facies/SConstruct b/book/rsf/tutorials/facies/SConstruct deleted file mode 100644 index df1f39c17c..0000000000 --- a/book/rsf/tutorials/facies/SConstruct +++ /dev/null @@ -1,13 +0,0 @@ -from rsf.proj import * - -# The dataset consists of 5 wireline log measurements, two indicator -# variables and a facies label at half foot intervals. -filename='facies_vectors.csv' - -Fetch(filename,'1610_Facies_classification', - server='https://raw.githubusercontent.com', - top='seg/tutorials-2016/master') - -Flow('vectors',filename,'csv2rsf') - -End() diff --git a/book/rsf/tutorials/nn/SConstruct b/book/rsf/tutorials/nn/SConstruct deleted file mode 100644 index 1d68f50167..0000000000 --- a/book/rsf/tutorials/nn/SConstruct +++ /dev/null @@ -1,807 +0,0 @@ -from __future__ import division -from rsf.proj import * - -# Sigmoid - -Flow('spike1',None,'spike n1=50 o1=-5 d1=%g | math output=x1'%(10/49)) -Flow('sigmaspike1','spike1','math output="1/(1+exp(-input))"') -Flow('spike2',None,'spike n1=80 o1=-6 d1=%g | math output=x1'%(12/78)) -Flow('sigmaspike2','spike2','math output="1/(1+exp(-input))"') -Plot('sigmaspike1','graph grid=y label1= unit1= title="Sigmoid function" max1=6.3 min1=-6.3') -Plot('sigmaspike2','graph wanttile=n wantaxis=n label1= unit1= dash=2 title="Sigmoid function" max1=6.3 min1=-6.3') -Result('sigma','sigmaspike1 sigmaspike2','Overlay') - -# Tanh - -Flow('tanhspike1','spike1','math output="(exp(input)-exp(-input))/(exp(input)+exp(-input))"') -Flow('tanhspike2','spike2','math output="(exp(input)-exp(-input))/(exp(input)+exp(-input))"') -Plot('tanhspike1','graph grid=y label1= unit1= title="Tanh function" max1=6.3 min1=-6.3') -Plot('tanhspike2','graph wanttile=n wantaxis=n label1= unit1= dash=2 title="Tanh function" max1=6.3 min1=-6.3') -Result('tanh','tanhspike1 tanhspike2','Overlay') - -# ReLU - -Flow('maskspike1','spike1','mask min=0 | dd type=float') -Flow('maskspike2','spike2','mask min=0 | dd type=float') -Flow('reluspike1','spike1 maskspike1','mul ${SOURCES[1]}') -Flow('reluspike2','spike2 maskspike2','mul ${SOURCES[1]}') -Plot('reluspike1','graph grid=y label1= unit1= title="ReLU function" max1=6.3 min1=-6.3') -Plot('reluspike2','graph wanttile=n wantaxis=n label1= unit1= dash=2 title="ReLU function" max1=6.3 min1=-6.3') -Result('relu','reluspike1 reluspike2','Overlay') - -# Load and process data - -# download well-log data -Fetch('R-39.las','1808_Neural_networks', - server='https://raw.githubusercontent.com', - top='seg/tutorials-2018/master') - -# Convert to RSF -Flow('R-39','R-39.las','las2rsf $SOURCE $TARGET',stdin=0,stdout=-1) - -# Examine with "< R-39.rsf sfheaderattr segy=n desc=y" - -# Skip the ends of the logs -for case in ('DT4P','DT4S','RHOB','DEPT'): - Flow(case,'R-39','headermath output=%s segy=n | window n2=7743' % case) - -Flow('VP1','DT4P','math output="1e6/input"') -Flow('VS1','DT4S','math output="1e6/input"') - -# Replace dodgy data in the VS log with mean -Flow('VSmask','VS1','mask max=0 | dd type=float') -Flow('VSmaskpos','VSmask','math output="1-input" | dd type=float') -Flow('VSpos','VS1 VSmaskpos','mul ${SOURCES[1]}') -Flow('VSneg','VSmask','math output="input*1973.8"') -Flow('VS2','VSpos VSneg','add ${SOURCES[1]}') - -# Take every 5th sample -Flow('VP','VP1','window j1=5 | window n1=501') -Flow('VS','VS2','window j1=5 | window n1=501') -Flow('RHO','RHOB','window j1=5 | window n1=501') -Flow('DEPTH','DEPT','window j1=5 | window n1=501') - -# Make upper layers and normalize - -Flow('VPupwn','VP','window n1=500 | rtoc') -Flow('VSupwn','VS','window n1=500 | rtoc') -Flow('RHOupwn','RHO','window n1=500 | rtoc') - -Flow('VPup','VP','window n1=500 | math output="(input-3461.27)/393.466"') -Flow('VSup','VS','window n1=500 | math output="(input-1905.71)/248.724"') -Flow('RHOup','RHO','window n1=500 | math output="(input-2455.2)/109.047"') -Flow('DEPTHup','DEPTH','window n1=500') - -# Make lower layers and normalize -Flow('VPlowwn','VP','window f1=1 n1=500 | rtoc') -Flow('VSlowwn','VS','window f1=1 n1=500 | rtoc') -Flow('RHOlowwn','RHO','window f1=1 n1=500 | rtoc') - -Flow('VPlow','VP','window f1=1 n1=500 | math output="(input-3461.43)/393.381"') -Flow('VSlow','VS','window f1=1 n1=500 | math output="(input-1905.76)/248.718"') -Flow('RHOlow','RHO','window f1=1 n1=500 | math output="(input-2455.9)/108.041"') -Flow('DEPTHlow','DEPTH','window f1=1 n1=500') - -# Angle of incidence - -def theta(target=None,source=None,env=None): - 'Angle of incidence' - import numpy, m8r - numpy.random.seed(42) - min_theta=0 - max_theta=20 - n=500 - theta = numpy.random.random(n)*(max_theta-min_theta)+min_theta - - rsf = m8r.Output(str(target[0])) - rsf.put("n1",n) - rsf.put("o1",0) - rsf.put("d1",1) - rsf.write(theta) - rsf.close() - return 0 - -Command('thetawn.rsf',None,action=Action(theta)) -Flow('theta','thetawn','math output="(input-9.97123)/5.97377"') - -#Flow('thetawn',None,'spike o1=0.10123167692437374 n1=500 d1=0.03959531912 | math output=x1 | put o1=0 d1=1') -#Flow('theta',None,'spike o1=0.10123167692437374 n1=500 d1=0.03959531912 | math output=x1 | put o1=0 d1=1 | math output="(input-9.98026)/5.7208"') - -# Create training data -Flow('X','VPup VSup RHOup VPlow VSlow RHOlow theta','cat ${SOURCES[1:7]} axis=2') - -# Create training label - -# Zoeppritz solution for P-P reflectivity -from math import pi -Flow('theta1','thetawn','math output="input*%g/180" | rtoc'%(pi)) -Flow('p','theta1 VPupwn','math x=${SOURCES[0]} y=${SOURCES[1]} output="sin(x)/y"') -Flow('theta2','p VPlowwn','math x=${SOURCES[0]} y=${SOURCES[1]} output="asin(x*y)"') -Flow('phi1','p VSupwn','math x=${SOURCES[0]} y=${SOURCES[1]} output="asin(x*y)"') -Flow('phi2','p VSlowwn','math x=${SOURCES[0]} y=${SOURCES[1]} output="asin(x*y)"') -Flow('a','RHOlowwn phi2 RHOupwn phi1','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*(1-2*sin(y)^2)-z*(1-2*sin(t)^2)"') -Flow('b','RHOlowwn phi2 RHOupwn phi1','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*(1-2*sin(y)^2)+2*z*sin(t)^2"') -Flow('c','RHOupwn phi1 RHOlowwn phi2','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*(1-2*sin(y)^2)+2*z*sin(t)^2"') -Flow('d','RHOlowwn VSlowwn RHOupwn VSupwn','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="2*(x*y^2-z*t^2)"') -Flow('E','theta1 VPupwn theta2 VPlowwn b c','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="(k*cos(x)/y)+(m*cos(z)/t)"') -Flow('F','phi1 VSupwn phi2 VSlowwn b c','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="(k*cos(x)/y)+(m*cos(z)/t)"') -Flow('G','a d theta1 VPupwn phi2 VSlowwn','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="x-y*cos(z)/t*cos(k)/m"') -Flow('H','a d theta2 VPlowwn phi1 VSupwn','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="x-y*cos(z)/t*cos(k)/m"') -Flow('D','E F G H p','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*y+z*t*k^2"') -Flow('rpp','D F b theta1 VPupwn c theta2 VPlowwn H p a d phi2 VSlowwn','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} l=${SOURCES[6]} n=${SOURCES[7]} o=${SOURCES[8]} p=${SOURCES[9]} q=${SOURCES[10]} r=${SOURCES[11]} s=${SOURCES[12]} f=${SOURCES[13]} output="(1/x) * (y*(z*(cos(t)/k) - m*(cos(l)/n)) - o*p^2 * (q + r*(cos(t)/k)*(cos(s)/f)))" | real') - -# Plot to QC - -for case in ('VP','VS','RHO'): - graph='' - if case=='VP': - graph+='plotcol=6 title="Vp" label2=Vp unit2="m/s"' - if case=='VS': - graph+='plotcol=5 title="Vs" label2=Vs unit2="m/s"' - if case=='RHO': - graph+='plotcol=4 title="Rho" label2=Rho unit2="kg/m3"' - if case=='rpp': - graph+='plotcol=7 title="Rpp" label2= unit2=' - Plot(case,['DEPTH',case], - ''' - cmplx ${SOURCES[1]} | window | - graph grid=y transp=y yreverse=y label1=Depth unit1=m labelsz=12 - ''' + graph) -Plot('rpp','put o1=2193.04 d1=0.762 | graph transp=y grid=y yreverse=y symbol=* symbolsz=10 labelsz=12 label2= unit2= title="Rpp" label1=Depth unit1=m min2=-0.4 max2=0.4') - -Result('QC','VP VS RHO rpp','SideBySideAniso') - -# Split into training and validation set - -Flow('Xtrain1 Ytrain1','X rpp','shuffle3 seed=42 axis=1 pi2=${SOURCES[1]} po2=${TARGETS[1]}') -Flow('Xtrain','Xtrain1','transp | window n1=400') -Flow('Ytrain','Ytrain1','transp | window n1=400') - -def numpy_load(target=None,source=None,env=None): - 'convert from numpy to RSF format' - import numpy, m8r - data = numpy.load(str(source[0])) - data = data.astype('float64') - #print(data.dtype) - rsf = m8r.Output(str(target[0])) - rsf.put("n2",data.shape[0]) - rsf.put("o2",2193.04) - rsf.put("d2",0.762) - rsf.put("n1",data.shape[1]) - rsf.put("o1",0) - rsf.put("d1",1) - rsf.put("label2","Time") - rsf.put("unit2","s") - rsf.write(data) - rsf.close() - return 0 - -def numpy_load2(target=None,source=None,env=None): - 'convert from numpy to RSF format' - import numpy, m8r - data = numpy.load(str(source[0])) - data=data.astype('float64') - rsf = m8r.Output(str(target[0])) - rsf.put("n1",data.shape[0]) - rsf.put("o1",2193.04) - rsf.put("d1",0.762) - rsf.put("n2",1) - rsf.put("o2",0) - rsf.put("d2",1) - rsf.put("label1","Time") - rsf.put("unit1","s") - rsf.write(data) - rsf.close() - return 0 - -# Load prepared data - -Command('Xtrain-1-1.rsf','X_train.npy',action=Action(numpy_load)) -Command('Ytrain-1.rsf','y_train.npy',action=Action(numpy_load2)) -Command('Yval-1.rsf','y_val.npy',action=Action(numpy_load2)) -Command('Xval-1-1.rsf','X_val.npy',action=Action(numpy_load)) -Command('Xblind-1-1.rsf','X_blind.npy',action=Action(numpy_load)) -Command('Yblind-1.rsf','y_blind.npy',action=Action(numpy_load2)) - -Flow('Xtrain-1','Xtrain-1-1','transp') -Flow('Xval-1','Xval-1-1','transp') -Flow('Xblind-1','Xblind-1-1','transp') - -# Create patch -Flow('Xtrainpatch','Xtrain','sfpatch w=1,7 p=400,1') - -Flow('Ytrainpatch','Ytrain','sfpatch w=1,1 p=400,1 | spray n=1 axis=2') - -Flow('Xtrainpatch2','Xtrain','sfpatch w=100,7') - -Flow('Ytrainpatch2','Ytrain','sfpatch w=100,1 | spray n=1 axis=2') - -Flow('Xtrainpatch3','Xtrain','sfpatch w=20,7') - -Flow('Ytrainpatch3','Ytrain','sfpatch w=20,1 | spray n=1 axis=2') - -# Validation data -Flow('Xval','Xtrain1','transp | window f1=400') -Flow('Yval','Ytrain1','transp | window f1=400') - -# Testing set - -# Take every 5th sample -Flow('VPtest','VP1','window j1=5 | window f1=1348') -Flow('VStest','VS2','window j1=5 | window f1=1348') -Flow('RHOtest','RHOB','window j1=5 | window f1=1348') -Flow('DEPTHtest','DEPT','window j1=5 | window f1=1348') - -# Make upper layers and normalize with training scaler - -Flow('VPupwntest','VPtest','window n1=200 | rtoc') -Flow('VSupwntest','VStest','window n1=200 | rtoc') -Flow('RHOupwntest','RHOtest','window n1=200 | rtoc') - -Flow('VPuptest','VPtest','window n1=200 | math output="(input-3461.27)/393.466"') -Flow('VSuptest','VStest','window n1=200 | math output="(input-1905.71)/248.724"') -Flow('RHOuptest','RHOtest','window n1=200 | math output="(input-2455.2)/109.047"') -Flow('DEPTHuptest','DEPTHtest','window n1=200') - -# Make lower layers and normalize with training scaler -Flow('VPlowwntest','VPtest','window f1=1 n1=200 | rtoc') -Flow('VSlowwntest','VStest','window f1=1 n1=200 | rtoc') -Flow('RHOlowwntest','RHOtest','window f1=1 n1=200 | rtoc') - -Flow('VPlowtest','VPtest','window f1=1 n1=200 | math output="(input-3461.43)/393.381"') -Flow('VSlowtest','VStest','window f1=1 n1=200 | math output="(input-1905.76)/248.718"') -Flow('RHOlowtest','RHOtest','window f1=1 n1=200 | math output="(input-2455.9)/108.041"') -Flow('DEPTHlowtest','DEPTHtest','window f1=1 n1=200') - -# Create testing data -Flow('Xtest','VPuptest VSuptest RHOuptest VPlowtest VSlowtest RHOlowtest thetatest','cat ${SOURCES[1:7]} axis=2') - -Flow('Xtestwn','VPupwntest VSupwntest RHOupwntest VPlowwntest VSlowwntest RHOlowwntest thetawntest','cat ${SOURCES[1:7]} axis=2') - -# Angle of incidence - -Flow('thetawntest',None,'spike o1=0 n1=200 d1=1 | math output=30 | put o1=0 d1=1') -Flow('thetatest',None,'spike o1=0 n1=200 d1=1 | math output=30 | put o1=0 d1=1 | math output="(input-9.98026)/5.7208"') - -# Create true testing label -Flow('theta1test','thetawntest','math output="input*%g/180" | rtoc'%(pi)) -Flow('ptest','theta1test VPupwntest','math x=${SOURCES[0]} y=${SOURCES[1]} output="sin(x)/y"') -Flow('theta2test','ptest VPlowwntest','math x=${SOURCES[0]} y=${SOURCES[1]} output="asin(x*y)"') -Flow('phi1test','ptest VSupwntest','math x=${SOURCES[0]} y=${SOURCES[1]} output="asin(x*y)"') -Flow('phi2test','ptest VSlowwntest','math x=${SOURCES[0]} y=${SOURCES[1]} output="asin(x*y)"') -Flow('atest','RHOlowwntest phi2test RHOupwntest phi1test','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*(1-2*sin(y)^2)-z*(1-2*sin(t)^2)"') -Flow('btest','RHOlowwntest phi2test RHOupwntest phi1test','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*(1-2*sin(y)^2)+2*z*sin(t)^2"') -Flow('ctest','RHOupwntest phi1test RHOlowwntest phi2test','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*(1-2*sin(y)^2)+2*z*sin(t)^2"') -Flow('dtest','RHOlowwntest VSlowwntest RHOupwntest VSupwntest','math x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="2*(x*y^2-z*t^2)"') -Flow('Etest','theta1test VPupwntest theta2test VPlowwntest btest ctest','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="(k*cos(x)/y)+(m*cos(z)/t)"') -Flow('Ftest','phi1test VSupwntest phi2test VSlowwntest btest ctest','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="(k*cos(x)/y)+(m*cos(z)/t)"') -Flow('Gtest','atest dtest theta1test VPupwntest phi2test VSlowwntest','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="x-y*cos(z)/t*cos(k)/m"') -Flow('Htest','atest dtest theta2test VPlowwntest phi1test VSupwntest','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} output="x-y*cos(z)/t*cos(k)/m"') -Flow('Dtest','Etest Ftest Gtest Htest ptest','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} output="x*y+z*t*k^2"') -Flow('rpptest','Dtest Ftest btest theta1test VPupwntest ctest theta2test VPlowwntest Htest ptest atest dtest phi2test VSlowwntest','math k=${SOURCES[4]} x=${SOURCES[0]} y=${SOURCES[1]} z=${SOURCES[2]} t=${SOURCES[3]} m=${SOURCES[5]} l=${SOURCES[6]} n=${SOURCES[7]} o=${SOURCES[8]} p=${SOURCES[9]} q=${SOURCES[10]} r=${SOURCES[11]} s=${SOURCES[12]} f=${SOURCES[13]} output="(1/x) * (y*(z*(cos(t)/k) - m*(cos(l)/n)) - o*p^2 * (q + r*(cos(t)/k)*(cos(s)/f)))" | real') - -# Plot to QC - -for case in ('VPtest','VStest','RHOtest'): - graph='' - if case=='VPtest': - graph+='plotcol=6 title="Vp" label2=Vp unit2="m/s"' - if case=='VStest': - graph+='plotcol=5 title="Vs" label2=Vs unit2="m/s"' - if case=='RHOtest': - graph+='plotcol=4 title="Rho" label2=Rho unit2="kg/m3"' - if case=='rpp': - graph+='plotcol=7 title="Rpp" label2= unit2=' - Plot(case,['DEPTHtest',case], - ''' - cmplx ${SOURCES[1]} | window | - graph grid=y transp=y yreverse=y label1=Depth unit1=m labelsz=12 - ''' + graph) -Plot('rpptest','put o1=3220.22 d1=0.762 | graph transp=y grid=y yreverse=y symbol=* symbolsz=10 labelsz=12 label2= unit2= title="Rpp" label1=Depth unit1=m min2=-0.12 max2=0.12') - -Result('QCtest','VPtest VStest RHOtest rpptest','SideBySideAniso') - -# Initialize parameters - -#Flow('W1-0-30',None,'spike n2=7 n1=300 mag=0 | noise seed=42 mean=0 var=1 | math output="input*sqrt(1/7)" | put d1=0.762 o1=2193.04 d2=1 o2=0') -Flow('W1-0-30-1 W2-0-30-1',None,'inittwolayer n1=300 n2=7 w2=${TARGETS[1]} mean=0 var=0.1 seed=42') -Flow('W1-0-30','W1-0-30-1','transp | put d1=0.762 o1=2193.04 d2=1 o2=0') -Flow('W2-0-30','W2-0-30-1','transp | put d1=0.762 o1=2193.04 d2=1 o2=0') - -Flow('b1-0-30',None,'spike n1=300 mag=0 n2=1 | put d1=0.762 o1=2193.04 d2=1 o2=0') -#Flow('W2-0-30',None,'spike n1=300 mag=0 n2=1 | noise seed=42 mean=0 var=1 | math output="input*sqrt(1/300)" | put d1=0.762 o1=2193.04 d2=1 o2=0') -#Flow('W2-0-30',None,'inittwolayer n1=300 n2=1 mean=0 var=0.1 seed=42') -Flow('b2-0-30',None,'spike n1=1 mag=0 n2=1 | put d1=0.762 o1=2193.04 d2=1 o2=0') -#i=1 - -Flow('vW1-0-30',None,'spike n2=7 n1=300 mag=0') -Flow('vb1-0-30',None,'spike n2=1 n1=300 mag=0') -Flow('vW2-0-30',None,'spike n2=300 n1=1 mag=0') -Flow('vb2-0-30',None,'spike n2=1 n1=1 mag=0') - -Flow('sW1-0-30',None,'spike n2=7 n1=300 mag=0') -Flow('sb1-0-30',None,'spike n2=1 n1=300 mag=0') -Flow('sW2-0-30',None,'spike n2=300 n1=1 mag=0') -Flow('sb2-0-30',None,'spike n2=1 n1=1 mag=0') - -Flow('lossm vlossm W1m b1m W2m b2m','Xtrainpatch Ytrainpatch Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=1 seed=42 stop=20 lossfunc=0 reg=0 alpha=0') - -Flow('lossa vlossa W1a b1a W2a b2a','Xtrainpatch Ytrainpatch Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.0001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=2 seed=42 stop=40 lossfunc=0 reg=0 alpha=0') - -Flow('lossa10 vlossa10 W1a10 b1a10 W2a10 b2a10','Xtrainpatch2 Ytrainpatch2 Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.0001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=2 seed=42 stop=40 lossfunc=0 reg=0 alpha=0') - -Flow('lossa10reg1 vlossa10reg1 W1a10reg1 b1a10reg1 W2a10reg1 b2a10reg1','Xtrainpatch3 Ytrainpatch3 Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.0001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=2 seed=42 stop=40 lossfunc=1 reg=1 alpha=0') - -Flow('loss-1a10 vloss-1a10 W1-1a10 b1-1a10 W2-1a10 b2-1a10','Xtrainpatch3 Ytrainpatch3 Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.0001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=2 seed=42 stop=40 lossfunc=1 reg=0 alpha=0') - -Flow('lossa20 vlossa20 W1a20 b1a20 W2a20 b2a20','Xtrainpatch3 Ytrainpatch3 Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.0001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=2 seed=42 stop=40 lossfunc=0 reg=0 alpha=0') - -Flow('loss vloss W1 b1 W2 b2','Xtrainpatch Ytrainpatch Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=0 opt=0 seed=42 stop=20 lossfunc=0 reg=0 alpha=0') - -Flow('losstanh vlosstanh W1tanh b1tanh W2tanh b2tanh','Xtrainpatch Ytrainpatch Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=1 opt=1 seed=42 stop=20 lossfunc=0 reg=0 alpha=0') - -Flow('lossrelu vlossrelu W1relu b1relu W2relu b2relu','Xtrainpatch Ytrainpatch Xval Yval W1-0-30 b1-0-30 W2-0-30 b2-0-30','twolayer label=${SOURCES[1]} valdata=${SOURCES[2]} vallabel=${SOURCES[3]} niter=100 lr=0.001 valloss=${TARGETS[1]} weight1out=${TARGETS[2]} bias1out=${TARGETS[3]} weight2out=${TARGETS[4]} bias2out=${TARGETS[5]} weight1=${SOURCES[4]} bias1=${SOURCES[5]} weight2=${SOURCES[6]} bias2=${SOURCES[7]} act=2 opt=1 seed=42 stop=20 lossfunc=0 reg=0 alpha=0') - -#lossistrain=[] -#lossisval=[] -#while i