- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Faculty Research and Publications /
- Seismic inversion through operator overloading
Open Collections
UBC Faculty Research and Publications
Seismic inversion through operator overloading 2008
pdf
Page Metadata
Item Metadata
Title | Seismic inversion through operator overloading |
Creator |
Herrmann, Felix J. |
Date Created | 2008-03-20 |
Date Issued | 2007 |
Description | Inverse problems in (exploration) seismology are known for their large to very large scale. For instance, certain sparsity-promoting inversion techniques involve vectors that easily exceed 230 unknowns while seismic imaging involves the construction and application of matrix-free discretized operators where single matrix-vector evaluations may require hours, days or even weeks on large compute clusters. For these reasons, software development in this field has remained the domain of highly technical codes programmed in low-level languages with little eye for easy development, code reuse and integration with (nonlinear) programs that solve inverse problems. Following ideas from the Symes’ Rice Vector Library and Bartlett’s C++ object-oriented interface, Thyra, and Reduction/Transformation operators (both part of the Trilinos software package), we developed a software-development environment based on overloading. This environment provides a pathway from in-core prototype development to out-of-core and MPI ’production’ code with a high level of code reuse. This code reuse is accomplished by integrating the out-of-core and MPI functionality into the dynamic object-oriented programming language Python. This integration is implemented through operator overloading and allows for the development of a coordinate-free solver framework that (i) promotes code reuse; (ii) analyses the statements in an abstract syntax tree and (iii) generates executable statements. In the current implementation, we developed an interface to generate executable statements for the out-of-core unix-pipe based (seismic) processing package RSF-Madagascar (rsf.sf.net). The modular design allows for interfaces to other seismic processing packages and to in-core Python packages such as numpy. So far, the implementation overloads linear operators and elementwise reduction/transformation operators. We are planning extensions towards nonlinear operators and integration with existing (parallel) solver frameworks such as Trilinos. |
Extent | 4719091 bytes |
Subject |
Software Python Inverse Problems Object Oriented SLIMpy Linear Operators Augmented Matrix |
Genre |
Other |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Seismic Laboratory for Imaging and Modeling |
Date Available | 2008-03-20 |
DOI | 10.14288/1.0107422 |
Affiliation |
Science, Faculty of Earth and Ocean Sciences, Department of |
Citation | Herrmann, Felix J. 2007. Seismic inversion through operator overloading. Conference on Applied Inverse Problems 2007. |
Peer Review Status | Unreviewed |
Scholarly Level | Faculty |
Copyright Holder | Felix J. Herrmann |
URI | http://hdl.handle.net/2429/604 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/340/items/1.0107422/source |
Download
- Media
- herrmann07AIP2.pdf
- herrmann07AIP2.pdf [ 4.5MB ]
- Metadata
- JSON: 1.0107422.json
- JSON-LD: 1.0107422+ld.json
- RDF/XML (Pretty): 1.0107422.xml
- RDF/JSON: 1.0107422+rdf.json
- Turtle: 1.0107422+rdf-turtle.txt
- N-Triples: 1.0107422+rdf-ntriples.txt
- Citation
- 1.0107422.ris
Full Text
Seismic inversion through operator overloading Felix J. Herrmann joint work with C. Brown, H. Modzelewski, G. Hennenfent, S. Ross Ross Seismic Laboratory for Imaging and Modeling http://slim.eos.ubc.ca AIP 2007, Vancouver, June 29 Synopsis Challenge: integrate & scale “Pipe”-based software with object-oriented abstraction fine-grained efficiency coarse-grained abstraction Opportunity: Create an object-oriented interface Implement algorithms modeled directly from the math Independent of the lower level software. (p)SLIMpy: A Collection of Python classes: vector, linear operator, R/T operation Benefits: Reusable End-Users Mathematicians Geophysicists Motivation Inverse problems in (exploration) seismology are large scale # of unknowns exceeds 230 matrix-free implementation of operators matrix-vector operations take hours, days, weeks Software development highly technical coding <=> reusable OO programming little code reuse emphasis on flows <=> iterations as part of nonlinear matrix-free optimization Motivation cont’d Interested in solving large-scale (nonlinear) optimization element-wise reduction/transformation operations linear matrix-free operators nonlinear operators (future) Design criteria code reuse & clarity seamless upscaling no overhead & possible speedup Opportunity Create an interface for the user to implement algorithms Object-oriented Interfaces the ANA (coordinate-free abstract numerical algorithms) with lower level software (Madagascar) Independent of the lower level software that can be in-core out-of-core (pipe-based) serial or MPI Our solution Operator overloading OO abstraction of vectors and operators coordinate free SLIMpy: a compiler for coarse grained “instruction set” for ANA’s concrete vector & operator class in Python Madagascar = instruction set Interactive scripting in Python math “beauty” at no performance sacrifice scalable (p)SLIMpy Provides a collection of Python classes and functions, to represent algebraic expressions Is a tool used to design and implement algorithms clean code; No more code than is necessary Designed to facilitate the transfer of knowledge between the end user and the algorithm designer Context SLIMpy is based on ideas list by: William Symes’ Rice Vector Library. http://www.trip.caam.rice.edu/txt/tripinfo/rvl.html http://www.trip.caam.rice.edu/txt/tripinfo/rvl2006.pdf Ross Bartlett's C++ object-oriented interface, Thyra. http://trilinos.sandia.gov/packages/thyra/index.html Reduction/Transformation operators (both part of the Trilinos software package). http://trilinos.sandia.gov/ PyTrilinos by William Spotz. http://trilinos.sandia.gov/packages/pytrilinos/ Who should use SLIMpy? Anyone working on large-to-extremely large scale optimization: NumPy, Matlab etc. etc. unix pipe-based (Madagascar, SU, SEPlib etc.) seamless migration from in-core to out-of-core to parallel Anyone who would like to produce code that is: readable & reusable deployable in different environments integrable with existing OO solver libraries Write solver once, deploy “everywhere” ... y = vector(‘data.rsf’) A1 = fdct2(domain=y.space).adj() A2 = fft2(domain=y.space).adj() A = aug_oper([A1, A2]) solver = GenThreshLandweber(10,5,thresh=None) x=solver.solve(A,y) Abstraction Let data be a vector y ∈ Rn. Let A1 := CT ∈ Cn×M be the inverse curvelet transform and A2 := FH ∈ Cn×n the inverse Fourier transform. Define A := [ A1 A2 ] and x = [ xT1 xT2 ]T Solve x̃ = arg min x ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! Vector & linear operator definition Math SLIMpy Matlab RSF y=data y=vector(‘data.rsf’) y=load(‘data’) <y.rsf A=CT C=linop(domain,range).adj() defined as function sffdct inv=y Instruction set for ANA’s Solvers consist of reduction/transformation operations that include element-wise addition, subtraction, multiplication vector inner products norms l1, l2 etc application of matrix-free matrix vector multiplications including adjoints composition of linear operators augmentation of linear operators Reduction transformation operations Math SLIMpy Matlab RSF y=a+b y=a+b y=a+b <a.rsf sfmath b=b.rsf output=input+b > y.rsf y=aTb y=inner(a,b) y=a’*b ? y=diag(a)*b y=a*b y=a.*b <a.rsf sfmath b=b.rsf output=input*b > y.rsf Linear operators Math SLIMpy Matlab RSF y=Ax y=A*x y=A(x) <x.rsf sffft2 > y.rsf z=AHy y=A.adj()*y z=A(y,’transp’) <y.rsf sffft2 inv=y > z.rsf A = [ A1 A2 ] A=aug_oper ([A1,A2]) new function complicated A = BCT A=comp ([B,C.adj()]) new function complicated Example - Define Linear Operator Define and apply a Linear Operator #User defined fft operator class fft_user(newlinop): command = "fft1" params = () def __init__(self,space,opt='y',inv='n'): self.inSpace = space self.kparams = dict(opt=opt,inv=inv) newlinop.__init__(self) #Initialize/Define the User created Linear Operator F = fft_user(vec1.getSpace()) final_create = F * vec1 Use predefined operators through a plug-in system. Create and import your own operators. Reuse linear operator def’s. Example - Compound Operator Define a compound operator: form new operators by combining predefined linear operators. adjoints automatically created. RM = weightoper(wvec=remove,inSpace=data.getSpace()) C = fdct3(data.getSpace(),cpxIn=True,*transparams) REMV=comp([C.adj(),RM,C]) PAD = Pad_wdl(data.getSpace(),padlen=padlen,winlen=winlen) F = fft3_wdl(PAD.range(),axis=1,sym='y',pad=1) T = transpoper(F.range(),memsize=2000,plane=[1,3]) PFT=comp([T,F,PAD]) P = Multpred_wdl(T.range(),filt=1,input=filterf) PEF = comp([PAD.adj(),F.adj(),T.adj(),P,PFT]) interp3d_remv.pyby Deli Wang Example - Augmented Matrix Define an augmented matrix: #Creating vector space of a 10 by 10 vecSpace = space(n1=10,n2=10, d... #Creating a vector of zeros from th... x = vecSpace.generateNoisyData() y = vecSpace.generateNoisyData() A = fft1(vecSpace) D = dwt(vecSpace) # Define the Augmented Matrix V = aug_vec( [[x], [y]]) L = aug_oper([[ D,A ], [ A,D ]]) # Multiple the two matrix’s RES = L*V Set up augmented linear systems. Close to visual representation of the matrix system. Produce testable code Automatic domain-range checking on linear operators. Automatic dot-test checking on all linear operators with optional flag. Allow dry-run for program testing. Dot-test check Use --dottest flag at command prompt. separate from the application at no overhead. Parses the script for the linear operators. Dot-test check Clean code domain = self.oper.domain() range = self.oper.range() domainNoise = self.generateNoisyData(domain) rangeNoise = self.generateNoisyData(range) self.operinv = self.oper.adj() trans = self.oper * domainNoise inv = self.operinv * rangeNoise tmp2 = rangeNoise * trans.conj() tmp1 = inv.conj() * domainNoise import numpy norm1 = numpy.sum(tmp1[:]) norm2 = numpy.sum(tmp2[:]) ratio = abs(norm1/norm2) self.assertAlmostEqual (norm1,norm2,places=11,msg=msg) print ratio if (0 == pid[4]) { /* reads from p[2], runs the adjoint, and writes to p[3] */ close(p[2][1]); close(STDIN_FILENO); dup(p[2][0]); close(p[3][0]); close(STDOUT_FILENO); dup(p[3][1]); argv[argc-1][4]='1'; execvp(argv[0],argv); _exit(4); } if (0 == pid[5]) { /* reads from p[3] and multiplies it with random model */ close(p[3][1]); close(STDIN_FILENO); dup(p[3][0]); pip = sf_input("in"); init_genrand(mseed); dp = 0.; for (msiz=nm, mbuf=nbuf; msiz > 0; msiz -= mbuf) { if (msiz < mbuf) mbuf=msiz; sf_floatread(buf,mbuf,pip); for (im=0; im < mbuf; im++) { dp += buf[im]*genrand_real1 (); } } sf_warning("L'[d]*m=%g",dp); SLIMpy Dot-Test Code Some of Madagascar’s 200 Line Dot-Test Code Python’s operators: +,-,*,/ are overloaded operator: * overloaded for linear op’s Linear operator constructor automatically generates adjoints & dot-tests keeps track of number type and domain & range Constructors for compound operators augmented systems Lurking problem w.r.t. efficiency ... y=A*x+A*z not as fast as y=A*(x+z) Observations Abstract Syntax Tree (AST) Operator overloading allows us to create an Abstract Syntax Tree (AST) Abstract Syntax Tree allows for analysis compound coarse-grained statements in ANA’s remove inefficiencies translate statements into a concrete instruction set do optimization Abstract Syntax Tree (AST) An AST is a finite, labeled, directed tree where: Internal nodes are labeled by operators Leaf nodes represent variables/vectors AST is used as an intermediate between a parse tree and a data structure. Executing single unix pipe-based commands is inefficient better to chain together reduce disk IO becomes complex when iterating Piped-based optimization /Volumes/Users/dwang/tools/python/2.5/bin/python ./interp3d_remv.py -n data=shot256.rsf filter=srme256.rsf remove=wac3D_256.rsf transparams=6,8,0 solverparams=2,1 eigenvalue=40.2 output=interp80_fcrsi.rsf mask=mask256_80.rsf flag=1 sfmath n1=256 n2=256 n3=256 output=0 | sffdct3 nbs=6 nbd=8 ac=0 maponly=y sizes=sizes256_256_256_6_8_0.rsf None < sfmath output="0" n2="256" n3="256" type="float" n1="256" | sfput head="tfile" o3="0" o2="0" d2="15" d3="15" o1="0" d1=".00" > tmp.walacs () shot256.rsf < sfcostaper nw1="15" nw2="15" > tmp.uNaFyB () None < sfmath output="0" n2="256" n3="256" type="float" n1="256" | sfput head="tfile.rsf" d2="15" o2="0" o3="0" d3="15" o1="0" d1=".00" > tmp.hyqCQa () None < sfcmplx tmp.uNaFyB tmp.hyqCQa | sfheadercut mask="mask256_80.rsf" | sffdct3 nbs="6" ac="0" nbd="8" inv="n" sizes="sizes256_256_256_6_8_0.rsf" > tmp.XiJ7Xs () wac3D_256.rsf < sfmath output="vec*input" vec="tmp.XiJ7Xs" | sffdct3 nbs="6" ac="0" sizes="sizes256_256_256_6_8_0.rsf" inv="y" nbd="8" > tmp.KdkI78 () srme256.rsf < sfcostaper nw1="15" nw2="15" > tmp.C1Qrmt () None < sfcmplx tmp.C1Qrmt tmp.walacs | sffdct3 nbs="6" ac="0" nbd="8" inv="n" sizes="sizes256_256_256_6_8_0.rsf" > tmp.TvVRiJ () wac3D_256.rsf < sfmath output="vec*input" vec="tmp.TvVRiJ" | sffdct3 nbs="6" ac="0" sizes="sizes256_256_256_6_8_0.rsf" inv="y" nbd="8" > tmp.DZjYCe () tmp.DZjYCe < sfmath output="0.0248756218905*input" | sfpad n1="512" | sffft3 inv="n" pad="1" sym="y" axis="1" | sftransp plane="13" memsize="2000" > tmp.F9Y0IZ () tmp.KdkI78 < sfheadercut mask="mask256_80.rsf" | sfpad n1="512" | sffft3 inv="n" pad="1" sym="y" axis="1" | sftransp plane="13" memsize="2000" | sfmultpred input="tmp.F9Y0IZ" adj="1" filt="1" | sftransp plane="13" memsize="2000" | sffft3 inv="y" pad="1" sym="y" axis="1" | sfwindow n1="256" | sffdct3 nbs="6" ac="0" nbd="8" inv="n" sizes="sizes256_256_256_6_8_0.rsf" > tmp.brima0 () tmp.brima0 < sfsort ascmode="False" memsize="500" > tmp.8rOLD8 () ${SCALAR01} = getitem(tmp.8rOLD8, 72299, ) ${SCALAR02} = getitem(tmp.8rOLD8, 71575815, ) Visualization dnoise.py OuterN=5 InnerN=5 --dot > dnoise_5x5.dot 3 iterations 25 iterations Optimization Currently implemented: Unix pipe-based optimization Unique to SLIMpy Assembles commands into longest possible pipe structures “Language” Specific Optimization Madagascar Goals within reach: Symbolic Optimization eg. A(x) + A(y) = A( x+y ) Parallel Optimization Load balancing Optimization Optimization of the AST is modular. Users can daisy chain each optimization function together specify which optimizations to perform #get the current AST Tree = getGraph() # perform Optimizations O1 = symbolicOptim( Tree ) # Perform Symbolic Optimizations O2 = pipe_Optim( O1 ) # Optimize for Unix pipe Structure O3 = language_rsf_Optim( O2 ) # Optimize for RSF Example Calculate tree Walkthrough Simple three step optimization # Create a Vector Space of a 10 x 10 image vecspace = VectorSpace(n1=10,n2=10,plugin=”rsf”) x = vecspace.zeros() # Vector of zeros y = vecspace.zeros() # Vector of zeros A = fft( vecspace ) # fft operator V = aug_vec( [[x],[y]] ) # augmented vector system M = aug_oper([[ A, 0 ], # augmented operator [ 0, A ] ) # system ans = M(V) # apply the operator to the vector O1 O2 O3Tree x y Zeros Zeros fft + fft Code AST Symbolic optimization Ax + Ay x y fft + fft O1 = symbolicOptim( Tree ) A(x + y) x y fft + O1 O2 O3Tree More efficient to add the vectors first. If we know A its a Linear Operator. only do one FFT computation Unix Pipe optimization O2 = pipe_Optim( O1 )x y fft + Zeros Zeros Zeros | add y | fft > ans Zeros > y O1 O2 O3Tree Compresses individual commands into minimal number of command pipes. Language specific O3 = language_rsf_Optim( O2 ) Zeros | add y | fft > ans Zeros > y sfmath output="0+y" | fft > ans Zeros > y O1 O2 O3Tree find shortcuts to reduce workload At What Cost What are the performance costs? linear time - with respect to the number of operations depends on the optimize functions used Performance cost 100 iterations of the solver dnoise.py OuterN=10 InnerN=10 --debug=display ... Display: Code ran in : 1.09 seconds Complexity : 1212 nodes Ran : 302 commands dnoise.py OuterN=30 InnerN=30 --debug=display ... Display: Code ran in : 6.51 seconds Complexity : 10812 nodes Ran : 2702 commands 900 iterations of the solver Pathway to parallel SLIMpy is scalable Includes parallellization Embarrassingly parallel Separate different branches of the AST Domain decomposition Slice the data-set into more manageable peaces Domain decomposition with partition of unity Embarrassingly parallel Separate branches of the AST can be easily distributed to different processors. cpu 1 ... cpu r ... cpu n xn Apply An V = aug_vec( [[x1], ..., [xn] ) # augmented vector system M = aug_oper([[ A1, ..., 0 ], # augmented operator [ 0, ..., Ar, ..., 0 ] [ 0, ..., An ] ) ans = M * V x1 Apply A1 xr Apply Ar Domain decomposition 4 windows along each dimension Taper functions shown on top Partition of unity 2D/3D Overlap arrays in-core with PETSc Dashed lines show window edges space in between lines is tapered overlap Domain decomposition dnoise.py serial parallel presto! mpirun [options] dnoise.py data=data.rsf output=res.rsf [pSLIMpy options] dnoise.py data=data.rsf output=res.rsf [pSLIMpy options] Linear operator A communicates overlap converts non-overlapping windows to overlapping ones Linear Operator B applies taper at edges of windows Combined operator T = BA satisfies: THT = I In Mathematical Terms... 〈Tx, y〉 = 〈x, THy〉 Performance Simple benchmark test 10 forward/transpose PWFDCT’s Decomposition Global Size Local Size Overlap Execution Time 2x2 1024x1024 544x544 16 51.97 630x630 64 66.94 2048x2048 1056x1056 16 201.85 1152x1152 64 237.83 4x4 2048x2048 544x544 16 52.67 630x630 64 67.26 4096x4096 1056x1056 16 202.53 1152x1152 64 240.25 8x8 4096x4096 544x544 16 61.89 630x630 64 87.84 8192x8192 1056x1056 16 220.68 1152x1152 64 257.50 Future Benefits of AST Integrate with the many software tools for AST analysis. Algorithm efficiency Advanced optimization techniques. Easily adapt AST optimizations to other platforms like Matlab. Observations End Users Mathematician: Implementing coordinate-free algorithms Easy to go from theory to practical applications Possible to compound and augment linear systems Geophysicist: Simplifying a large process flow to a single line No compromising of functionality Easier to control Parameter and domain-range tests Pathway to Reproducible Research Reuse SLIMpy code for a number of applications. Concrete implementation with overloading in Python. Focal Transform Interpolation SLIMpy App by Deli Wang SLIMpy code very compact. Interpolation around 25 lines. Easily readable. Wrapped in a SConstruct file for easy use of Reproducible Research. Set parameters in the SConstruct and run the application with one command. Uses reusable ANA. Picture by Deli Wang Primary operator Frequency slice from data matrix with dominant primaries. Receivers Shots Shots Receivers Frequency ∆P Solve Curvelet-based processing 3 SPARSITY-PROMOTING INVERSION Our solution strategy is built on the premise that seismic data and images have a sparse representation, x0, in the curvelet domain. To exploit this property, our forward model reads y = Ax0 + n (1) with y a vector with noisy and possibly incomplete mea- surements; A the modeling matrix that includes CT ; and n, a zero-centered white Gaussian noise. Because of the redundancy of C and/or the incompleteness of the data, the matrix A can not readily be inverted. However, as long as the data, y, permits a sparse vector, x0, the ma- trix, A, can be inverted by a sparsity-promoting program (Candès et al., 2006b; Donoho, 2006) of the following type: P! : { x̃ = argminx ‖x‖1 s.t. ‖Ax− y‖2 ≤ ! f̃ = ST x̃ (2) in which ! is a noise-dependent tolerance level, ST the inverse transform and f̃ the solution calculated from the vector x̃ (the symbol ˜ denotes a vector obtained by non- linear optimization) that minimizes P!. Nonlinear programs such as P! are not new to seismic data processing and imaging. Refer, for instance, to the extensive literature on spiky deconvolution (Taylor et al., 1979) and transform-based interpolation techniques such as Fourier-based reconstruction (Sacchi and Ulrych, 1996). By virtue of curvelets’ high compression rates, the non- linear program P! can be expected to perform well when CT is included in the modeling operator. Despite its large- scale and nonlinearity, the solution of the convex problem P! can effectively be approximated with a limited (< 250) number of iterations of a threshold-based cooling method derived from work by Figueiredo and Nowak (2003) and Elad et al. (2005). Each step involves a descent projection, followed by a soft thresholding. SEISMIC DATA RECOVERY The reconstruction of seismic wavefields from regularly- sampled data with missing traces is a setting where a curvelet-based method will perform well (see e.g. Herr- mann, 2005; Hennenfent and Herrmann, 2006a, 2007). As with other transform-based methods, sparsity is used to reconstruct the wavefield by solving P!. It is also shown that the recovery performance can be increased when in- formation on the major primary arrivals is included in the modeling operator. Curvelet-based recovery The reconstruction of seismic wavefields from incomplete data corresponds to the inversion of the picking operator R. This operator models missing data by inserting zero traces at source-receiver locations where the data is miss- ing. The task of the recovery is to undo this operation by filling in the zero traces. Since seismic data is sparse in the curvelet domain, the missing data can be recovered by compounding the picking operator with the curvelet modeling operator, i.e., A := RCT . With this defini- tion for the modeling operator, solving P! corresponds to seeking the sparsest curvelet vector whose inverse curvelet transform, followed by the picking, matches the data at the nonzero traces. Applying the inverse transform (with S := C in P!) gives the interpolated data. An example of curvelet based recovery is presented in Figure 1, where a real 3-D seismic data volume is recov- ered from data with 80% traces missing (see Figure 1(b)). The missing traces are selected at random according to a discrete distribution, which favors recovery (see e.g. Hen- nenfent and Herrmann, 2007), and corresponds to an av- erage sampling interval of 125m . Comparing the ’ground truth’ in Figure 1(a) with the recovered data in Figure 1(c) shows a successful recovery in case the high-frequencies are removed (compare the time slices in Figure 1(a) and 1(c)). Aside from sparsity in the curvelet domain, no prior information was used during the recovery, which is quite remarkable. Part of the explanation lies in the curvelet’s ability to locally exploit the 3-D structure of the data and this suggests why curvelets are successful for complex datasets where other methods may fail. Focused recovery In practice, additional information on the to-be-recovered wavefield is often available. For instance, one may have access to the predominant primary arrivals or to the ve- locity model. In that case, the recently introduced focal transform (Berkhout and Verschuur, 2006), which ’decon- volves’ the data with the primaries, incorporates this addi- tional information into the recovery process. Application of this primary operator,∆P, adds a wavefield interaction with the surface, mapping primaries to first-order surface- related multiples (see e.g. Verschuur and Berkhout, 1997; Herrmann, 2007). Inversion of this operator, strips the data off one interaction with the surface, focusing pri- maries to (directional) sources, which leads to a sparser curvelet representation. By compounding the non-adaptive curvelet transform with the data-adaptive focal transform, i.e., A := R∆PCT , the recovery can be improved by solving P!. The solution of P! now entails the inversion of ∆P, yielding the spars- est set of curvelet coefficients that matches the incomplete data when ’convolved’ with the primaries. Applying the inverse curvelet transform, followed by ’convolution’ with ∆P yields the interpolation, i.e. ST :=∆PCT. Compar- ing the curvelet recovery with the focused curvelet recov- ery (Fig ?? and ??) shows an overall improvement in the recovered details. SEISMIC SIGNAL SEPARATION Predictive multiple suppression involves two steps, namely multiple prediction and the primary-multiple separation. In practice, the second step appears difficult and adap- Recovery with focussing wit A := R∆PCT ST := ∆PCT y = RP(:) R = picking operator. Focal Transform Interpolation SConstruct file used to generated result. Efficient flow system only re-calculates what has changed. Flow('wac3D_256',None, ''' math n1=55521587 output="1" |sfput d1=0.004 o1=0| pad beg1=16777216>rwac3D.rsf && sfmath <rwac3D.rsf output="input*0" >iwac3D.rsf && sfcmplx rwac3D.rsf iwac3D.rsf|sfput d1=0.004 o1=0 o2=0 o3=0 n2=1 n3=1 d2=15 d3=15>$TARGET && sfrm rwac3D.rsf iwac3D.rsf ''',stdout=-1) Flow('interp80_fcrsi',['shot256','srme256','wac3D_256','mask256_80',interp_focal_remv], python+' '+interp_focal_remv + ' -n ' + """ data=${SOURCES[0]} filter=${SOURCES[1]} remove=${SOURCES[2]} transparams=6,8,0 solverparams=2,1 eigenvalue=40.2 output=$TARGET mask=${SOURCES[3]} flag=1 """,stdin=0,stdout=-1) Result('interp80','interp80_fcrsi.rsf',cube('Interp80% missing')) Result('cubshot','shot256.rsf',cube('shot')) End() Focal Transform Interpolation This small Focal Transform Interpolation application by Deli Wang produces... TAP = Taper_wdl(data.getSpace(),taplen=15) data = TAP*data filter = TAP*filter data = cmplx(data,data.getSpace().zeros()) filter = cmplx(filter,filter.getSpace().zeros()) #Remove fine scale wavelet operator RM = weightoper (wvec=remove,inSpace=data.getSpace()) #2D/3D, complex/real, curvelet transform C = fdct3(data.getSpace(),cpxIn=True,*transparams) PAD = Pad_wdl(data.getSpace (),padlen=padlen,winlen=winlen) #1D complex in/out FFT F = fft3_wdl(PAD.range(),axis=1,sym='y',pad=1) #Data cube transpose T = transpoper(F.range(),memsize=2000,plane=[1,3]) D = weightoper (wvec=eigenvalue,inSpace=filter.getSpace()) R = pickingoper(data.getSpace(),mask) REAL = Real_wdl(data.getSpace()) PFT = comp([T,F,PAD]) #Remove fine scale wavelet operator REMV = comp([C.adj(),RM,C]) filter = REMV*filter filterf = PFT*(D*filter) P = Multpred_wdl(T.range(),filt=1,input=filterf) #WCC in frequency domain PEF = comp([PAD.adj(),F.adj(),T.adj(),P,T,F,PAD]) #Define global linear operator if (flag==1): A = comp([R,PEF.adj(),C.adj()]) else: A = comp([R,C.adj()]) data = R*data data = REMV*data #Define and run the solver thresh = logcooling(thrparams[0],thrparams[1]) solver = GenThreshLandweber(solverparams [0],solverparams[1],thresh=thresh) x = solver.solve(A,data) interp3d_remv.py by Deli Wang Focal Transform Interpolation ...over 30 pipes utilizing more than 90 commands to obtain its result. /Volumes/Users/dwang/tools/python/2.5/bin/python ./interp3d_remv.py -n data=shot256.rsf filter=srme256.rsf remove=wac3D_256.rsf transparams=6,8,0 solverparams=2,1 eigenvalue=40.2 output=interp80_fcrsi.rsf mask=mask256_80.rsf flag=1 sfmath n1=256 n2=256 n3=256 output=0 | sffdct3 nbs=6 nbd=8 ac=0 maponly=y sizes=sizes256_256_256_6_8_0.rsf None < sfmath output="0" n2="256" n3="256" type="float" n1="256" | sfput head="tfile" o3="0" o2="0" d2="15" d3="15" o1="0" d1=".00" > tmp.walacs () shot256.rsf < sfcostaper nw1="15" nw2="15" > tmp.uNaFyB () None < sfmath output="0" n2="256" n3="256" type="float" n1="256" | sfput head="tfile.rsf" d2="15" o2="0" o3="0" d3="15" o1="0" d1=".00" > tmp.hyqCQa () None < sfcmplx tmp.uNaFyB tmp.hyqCQa | sfheadercut mask="mask256_80.rsf" | sffdct3 nbs="6" ac="0" nbd="8" inv="n" sizes="sizes256_256_256_6_8_0.rsf" > tmp.XiJ7Xs () wac3D_256.rsf < sfmath output="vec*input" vec="tmp.XiJ7Xs" | sffdct3 nbs="6" ac="0" sizes="sizes256_256_256_6_8_0.rsf" inv="y" nbd="8" > tmp.KdkI78 () srme256.rsf < sfcostaper nw1="15" nw2="15" > tmp.C1Qrmt () None < sfcmplx tmp.C1Qrmt tmp.walacs | sffdct3 nbs="6" ac="0" nbd="8" inv="n" sizes="sizes256_256_256_6_8_0.rsf" > tmp.TvVRiJ () wac3D_256.rsf < sfmath output="vec*input" vec="tmp.TvVRiJ" | sffdct3 nbs="6" ac="0" sizes="sizes256_256_256_6_8_0.rsf" inv="y" nbd="8" > tmp.DZjYCe () tmp.DZjYCe < sfmath output="0.0248756218905*input" | sfpad n1="512" | sffft3 inv="n" pad="1" sym="y" axis="1" | sftransp plane="13" memsize="2000" > tmp.F9Y0IZ () tmp.KdkI78 < sfheadercut mask="mask256_80.rsf" | sfpad n1="512" | sffft3 inv="n" pad="1" sym="y" axis="1" | sftransp plane="13" memsize="2000" | sfmultpred input="tmp.F9Y0IZ" adj="1" filt="1" | sftransp plane="13" memsize="2000" | sffft3 inv="y" pad="1" sym="y" axis="1" | sfwindow n1="256" | sffdct3 nbs="6" ac="0" nbd="8" inv="n" sizes="sizes256_256_256_6_8_0.rsf" > tmp.brima0 () tmp.brima0 < sfsort ascmode="False" memsize="500" > tmp.8rOLD8 () ${SCALAR01} = getitem(tmp.8rOLD8, 72299, ) ${SCALAR02} = getitem(tmp.8rOLD8, 71575815, ) Future goals MPI Interfacing through SLIMpy Allows for a user defined MPI definition. Will work with a number of different MPI platforms (Mpich, LAM). Easily call MPI features through SLIMpy script. Integration with existing OO solver frameworks such as Trilinos Automatic Differentiation Conclusions Use a scripting language to access low-level implementations of (linear) operators (seismic processing tools). Easy to use automatic checking tools such as domain- range checks and dot-test. Overloading and abstraction with small overhead. AST allows for optimization. Reusable ANAs and Applications. Is growing into a “compiler” for ANA’s .... SLIMpy Web Pages More information about SLIMpy can be found at the SLIM Homepage: Auto-books and tutorials can be found at the SLIMpy Generated Websites: http://slim.eos.ubc.ca http://slim.eos.ubc.ca/SLIMpy Acknowledgments Madagascar Development Team Sergey Fomel CurveLab 2.0.2 Developers Emmanuel Candes, Laurent Demanet, David Donoho and Lexing Ying SINBAD project with financial support BG Group, BP, Chevron, ExxonMobil, and Shell SINBAD is part of a collaborative research and development grant (CRD) number 334810-05 funded by the Natural Science and Engineering Research Council (NSERC)
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 6 | 0 |
China | 5 | 0 |
Germany | 3 | 1 |
France | 2 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 10 | 1 |
Beijing | 5 | 0 |
Ashburn | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: