workflow for inverse problems: make the inversion tools more standard and more user-friendly

[mpbl]

Update from Dimitri, Oct 2012: after talking with many developers and users, we should definitely make the inversion tools in SPECFEM3D more standard and more user friendly. Even if there are many options to adjust when solving inverse problems, and different users currently have different tools (in Matlab, Perl, Fortran, Python and so on), sooner rather than later we should probably try to define a single tool in which all these different options could be implemented, but at least we would have a unified framework, which we could document in the manual.

Currently many users complain that computing sensitivity kernels is easy with SPECFEM3D but that then it is not so easy to use them (and different people give different answers when they are asked). So far, at least the following groups have local tools: the Princeton group, the Toronto group, Carl in Alaska, Emanuele and his group at INGV, Alessia (at least for Flexwin, there are now two versions: in Fortran and in Python), and the Marseille / CNRS group. We should try to address this issue relatively quickly, because the more we wait the more these sets of tools will diverge.

One thing we should do quickly for sure is make sure we maintain the current tools of SPECFEM3D inside the same directory as SPECFEM3D itself; it seems that currently everything is outside, in an "ADJOINT_TOMO" directory, which is thus not downloaded by users when they get SPECFEM3D through SVN or from the CIG tar file; which is problematic.

[QuLogic]

Since the move to GitHub, most of the ADJOINT_TOMO directory has not also been exported. The fact that it's on SVN instead of here makes it one step farther for potential users.

[komatits]

Hi,

I agree, we should move it to Git as well.

Thanks, Dimitri.

On 11/15/2013 01:51 AM, Elliott Sales de Andrade wrote:

Since the move to GitHub, most of the |ADJOINT_TOMO| directory has not also been exported. The fact that it's on SVN instead of here makes it one step farther for potential users.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/23#issuecomment-28538671.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

[komatits]

ADJOINT_TOMO is still not moved to GitHub... that's bad. It seems to still be at https://geodynamics.org/svn/cig/seismo/3D/ADJOINT_TOMO/ADJOINT_TOMOGRAPHY_TOOLKIT

[komatits]

Vadim Monteiller tells me that Ludovic Métivier from Grenoble (France) has also developed a very nice package for that.

[komatits]

There is also OpenTOAST from Germany: http://www.opentoast.de/ and http://www.opentoast.de/Forward-modelling-code_SPECFEM3D.php

[komatits]

More details from Vadim:

"C'est Ludovic Métivier qui a développé les outils d'inversion pour SEISCOPE. Il a écrit une bibliothèque Fortran avec des entrées-sorties standardisées, du coup son package s'implémente facilement pour n'importe quel problème, et il a mis plein de méthodes dedans (CG, BFGS...) que l'on choisit avec des flags. L'utilisateur n'utilise qu'un seul call à une subroutine, le reste est géré par sa bibliothèque. Je crois que pour l'instant il ne distribue pas encore son outil car il est encore sous SEISCOPE, mais à terme (je crois qu'il faut attendre un an) il va mettre son code sous une licence type GPL. "

[komatits]

From Daniel @danielpeter Date: Mon, 12 May 2014 21:12:37 +0200

Haskell http://en.wikipedia.org/wiki/Haskell_%28programming_language%29 might be useful for developing an inverse framework, since it is easily taking care of dependencies. Max has written some small inversion toolbox for a toy problem with Haskell and became quite fond of it.

With Andreas Fichtner, we also just started a little project here to see how writing codes in chapel would work out performance-wise: http://en.wikipedia.org/wiki/Chapel_(programming_language) although chapel is at the moment not focussed on performance and it will probably need some tricks to have it run at speed.

It seems to me that especially the computer science people here like newer languages. they tend to say that computer scientists are more experimental than probably we as domain scientists.

[komatits]

From @komatits May 13, 2014

Regarding inversion frameworks, let me cc Max, Vadim, Hejun and Ludovic Métivier because it seems different people are currently developing different frameworks for that (the "ADJOINTTOMOGRAPHY" package by Hejun, Ryan et al, Vadim's Fortran routines, Ludovic's very nice package (according to Vadim that package is very flexible), and Max's new routines programmed in Haskell language ( http://en.wikipedia.org/wiki/Haskell%28programming_language%29 ). Sébastien Chevrot also has some processing tools written in MATLAB.

I was wondering if it would make sense at some point to discuss and see if all of us could try to come up with a unified and open-source package, to avoid duplicating the efforts. In a few months we could maybe organize a Skype call about this.

(Ludovic's package is developed within the SEISCOPE consortium of Jean Virieux http://seiscope2.osug.fr/ thus I do not think it is currently open source; Ludovic, do you think it could become open source in the future?).

[komatits]

From @komatits May 13, 2014

Nocedal himself has released his optimized implementation of L-BFGS-B version 3.0 open source:

http://www.ece.northwestern.edu/~nocedal/lbfgsb.html

[komatits]

From Jeroen Tromp, May 13, 2014:

I very much agree that we need to organize our efforts with regards to an inversion framework. There are lots of uncoordinated developments, and, like the SPECFEM source code, we should bring these under one umbrella that we all share.

In this context I should mention:

1. Matthieu Lefebvre has been developing a workflow for seismic imaging & inversion. He is using NetworkX and other graph tools. This workflow is basically the same as the adjoint tomography workflow.
2. James Smith, Lion Krischer, Wenjie Lei and Ebru Bozdag have been working on the pre-processing workflow using ADIOS/HDF5 as a new data format, called the Adaptable Seismic Data Format (ASDF).

Others actively involved in adjoint tomography are Min Chen and Christina Morency.

[komatits]

From Matthieu Lefebvre @mpbl May 13, 2014:

I agree that a functional language, in particular Haskell, is the most natural way for such a framework. However it is not commonly available on HPC clusters and the learning curve is probably too steep for a geophysicist.

The workflow Jeroen is mentioning has two (Python) packages, for now:

1) FlowJobs -- uses NetworkX to define direct acyclic graphs. Nodes are thought to be SAGA jobs and will be extended to allow BigJobs jobs (that is SAGA + a pilot framework). These APIs take care of machine / file transfer protocol dependant operations. When a job fails, all its successors are marked as failed, and therefore not executed. It is possible to dump the workflow graph at the end. An improvement will be to reload the sub-graph of failed jobs to re-run it.

2) The workflow itself defines nodes and directed edges between them. The whole set of operations is then run, according to job dependencies.

For these two packages, I have private bitbucket repositories that I can give access to.

James, Wenjie and Ryan have their own sets of processing tools that could be chained using FlowJobs. We might end-up defining a library of nodes that users can pick in order to define their own workflow graphs.

[komatits]

From Ludovic Métivier May 13, 2014:

I confirm a Fortran 90 package for large scale nonlinear optimization with bound constraints has been developed within the SEISCOPE consortium. The set of routines includes

• steepest descent
• nonlinear conjugate gradient
• L-BFGS
• truncated Newton methods

All the methods are implemented using a reverse communication protocole. A communication FLAG is used to tell the user whenever the following actions are required:

• computation of the misfit function and its gradient
• application of a preconditioner
• computation of a Hessian-vector product (only for truncated Newton methods)

The methods all use the same linesearch strategy (which satisfies the Wolfe criterion).

The access to this package is currently restricted to the sponsors of the SEISCOPE consortium. However, the code should become open-source in the next months, and available to be included into the development of any inversion codes.

[komatits]

OK, great. Very useful. I have updated this old thread/issue and cut and pasted the emails there.

In a few weeks/months let us probably have a Skype call about this.

[komatits]

Dear Ludovic,

Thank you very much for your answer. Very useful. We will be glad to try it and cite the related papers when it becomes open source; please keep us informed when you release it in a few months.

[komatits]

List of people who are potentially interested (in no particular order):

Dimitri Komatitsch, Matthieu Lefebvre, Jeroen Tromp, Daniel Peter, Ryan Modrak, Max Rietmann, Vadim Monteiller, Hejun Zhu, Sébastien Chevrot, Yi Wang, Qinya Liu, Elliott Sales de Andrade, Thierry Scotti, Ludovic Métivier, Paul Cristini, James A. Smith, Ebru BOZDAG, Wenjie Lei, Lion Krischer, Christina Morency, Min Chen, Emanuele Casarotti, Federica Magnoni, Clément Durochat, Alexis Bottero, Zheng Guangying 烟消ゞ云散, Mathilde Griveaux, Cédric Bellis, Jed Brown, Andreas Fichtner.

[komatits]

If we use functional languages such as Haskell we need to make sure that we will have people who can maintain that in the future (i.e. computer science people; I am not sure geophysicists could maintain that efficiently)

[jedbrown]

I'm interested in the requirements and rationale. Though I have written a fair amount of Haskell, it's not clear to me why Haskell would be considered for this purpose.

[rietmann]

Regarding the programming language discussion: In my (admittedly limited) experience Haskell’s beauty really takes a bump once one has non-side effect free operations like I/O and error handling. When dealing with real data, where stuff tends to go wrong, those two things are vital. I would like to proven wrong though. Aside from that I agree with the statement that it is just too different and adoption by domain scientists will most likely not happen.

I would argue that Haskell is ideally suited for handling errors from I/O operations, due to the expressive type system. I would say that Haskell's strengths are:

• Easy to write inversion in a "declarative" manner, where gradient calculations, updates, and simulations are tied together via a dependency system. This helps keep the code a bit simpler and makes it easier to split the code.
• Excellent error handling, logging, etc. can be expressed as part of the type system, allowing the software to be surprisingly elegant and robust.
• Excellent high-performance vector libraries and C-Interop (Haskell is statically compiled)
• Excellent parsing libraries make it simpler to deal with the many text-based config fies.

Downsides are of course:

• Lack of many libraries suited for seismology.
• Relatively steep learning curve

Python is the obvious choice for such a library, something I use myself, especially given the excellent scientific libraries like ObsPy. However, the lack of a type system makes larger Python libraries/packages, in my opinion, less robust. A sweet spot between Python and Haskell however is still missing. (there is an effort called MyPy which is bringing a type system to Python, but is still very new). Most companies I hear of that deploy large Python codes have a significant amount of testing to avoid many of these problems, but I find that building and maintaining a strong test suite is often difficult for academic projects (such as SPECFEM3D, which is lacking in this regard)

From my perspective the final solution will likely be a hybrid solution. I chose Haskell to develop a dependency-based system for computing inversions, but I haven't yet decided if it's the best solution. The most important thing for me is to be able to encode the inversion processes in a dependency graph, such that most steps are automated, eliminating simple mistakes.

[komatits]

From Lion Krischer @krischer
(one of the main developers of ObsPy (http://www.obspy.org) ) May 14, 2014:

Regarding the programming language discussion: In my (admittedly limited) experience Haskell’s beauty really takes a bump once one has non-side effect free operations like I/O and error handling. When dealing with real data, where stuff tends to go wrong, those two things are vital. I would like to proven wrong though. Aside from that I agree with the statement that it is just too different and adoption by domain scientists will most likely not happen.

I also tend to agree with Dimitri: C/C++/Fortran for HPC applications and most likely Python for the rest. With ObsPy (http://www.obspy.org) we solved most of our domain specific processing issues and the large Python ecosystem takes care of the rest. The most likely future contender in my opinion is the Julia language which has been designed for technical computing but the libraries are not there just yet.

We also developed a framework for dealing with adjoint inversions with seismological data. It mostly takes care of processing and data organisation and attempts to give some structure to what goes in and what comes out of each iteration. It still has some rough edges but should be done in the next couple of months. I think it could integrate nicely with the developments of Matthieu and Ludovic as each takes care of a different aspect of the whole workflow.

http://krischer.github.io/LASIF

Cheers!

Lion

[komatits]

From Ludovic Métivier, May 14, 2014: (translation of the first sentence: "To be more precise, my inversion tools will be made available open-source sometime in July on my web page"):

Pour être plus précis, le code sera dispo courant juillet sur ma page web. La politique SEISCOPE est de laisser le code un an à disposition exclusive des sponsors, ensuite de quoi on peut les proposer en open-source.

C'est super si ces codes vous intéressent et qu'ils peuvent vous être utiles. Je mettrai les deux papiers à citer en référence: c'est des applis FWI à l'échelle crustale où on regarde l'intérêt du Newton tronqué pour l'inversion.

[komatits]

Ludovic Métivier's tools are now open source.

[komatits]

ADJOINT_TOMO is now moved to GitHub, in a sub-directory of SPECFEM3D called ADJOINT_TOMOGRAPHY_TOOLS. Please maintain it there (only) from now on. Later today Eric @eheien is going to close the old SVN repository https://geodynamics.org/svn/cig/seismo/3D/ADJOINT_TOMO/ADJOINT_TOMOGRAPHY_TOOLKIT . That old svn repo was not used by anyone and thus the routines were quickly becoming obsolete and unused. By maintaining it on GitHub jointly with the rest things should quickly go back to normal and that great set of routines will be actively used and maintained again! We successfully did that for GEOCUBIT a few months ago, when we were in the same situation for that package (each group had its own local version, drifting from others of course).

[krischer]

Any reason for this not having its own repository? I am not sure how tight this package is tied to a specific SPECFEM version but it sounds like something that should live on its own.

[komatits]

Hi Lion,

In principle yes, in practice no, we tried this several times for several packages (including this one) in the past and it always failed, people just kept using obsolete or modified local versions and not committing anything. Having a single repo is the only thing that works in practice otherwise people just do not download the rest of the stuff.

Best wishes, Dimitri.

On 31/07/2014 17:53, Lion Krischer wrote:

Any reason for this not having its own repository? I am not sure how tight this package is tied to a specific SPECFEM version but it sounds like something that should live on its own.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/23#issuecomment-50778572.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

[komatits]

From Carl Tape @carltape to Emanuele @casarotti : Following some other discussions about the inversion workflow, there is quite a bit of work to do to get the workflow in order prior to Sept 1st. ADJOINT_TOMOGRAPHY_TOOLS is far from production mode right now.

Even with all adjoint sources made for the first iteration, there will still be much to sort out between iterations.

[komatits]

From Carl Tape @carltape : I followed a few different conversations on this and see that you brought ADJOINT_TOMOGRAPHY_TOOLS into SPECFEM3D.

https://github.com/geodynamics/specfem3d/tree/devel/ADJOINT_TOMOGRAPHY_TOOLS

These tools are definitely not at the standard of the rest of the package. The main chunk represented my 2007-8 scripts for the 2009 Science paper. Everything I needed was there, but it was written in a grad-student-get-the-job-done manner. The directories I used were

flexwin iterate_adj measure_adj

As I recall, Hejun adapted some of these tools into the directory ADJOINT_TOMOGRAPHY_TOOLKIT. That was intended to be stand-alone, too, but it was also implemented for the global code, so some changes would be needed. Also, I think there was an adapted program for the measurement/adjoint source code (measure_adj) above.

mtadj was a somewhat simpler version that Qinya Liu wrote. But, again, these were being stored for not-so-wide distribution. Obviously the redundancies are not a great thing, especially under the SPECFEM3D hood.

One thing that immediately seems to be an issue is that flexwin is already here as a stand-alone package: https://github.com/geodynamics/flexwin

Really it does not make much sense as a "CIG code," since it is a measurement tool within a larger workflow. But that's where it resides now. I cced Alessia here. She was working on a python version of flexwin. She also may have an opinion on whether to put it under SPECFEM3D or outside.

My understanding is that there are various "workflows" out there for time-domain (and freq-domain) adjoint-based inversions. (Quotes apply to ADJOINT_TOMOGRAPHY_TOOLS, too, since it is more a collection of scripts.) I think Lion is working on this: https://github.com/krischer/LASIF Probably that is close to what Andreas Fichtner is using. Po Chen's group has a different set of tools, not open-source as far as I know.

It has clearly been a challenge to converge on a clean workflow for inversions, though I know there are lots of efforts to tackle pieces of this (or all of it).

Sorry, this is a lot of text without direct recommendations. I agree that this is a problem. With the recent move, we have something inside SPECFEM3D, but it is probably not the best path forward. But it is a path and one that worked for me in 2008 -- or for Hejun in 2011.

[komatits]

Hi Carl @carltape ,

Thanks a lot. That is very useful.

Three other things that are worth checking (in particular the first one):

• Ludovic Métivier's very nice and efficient open-source set of tools: http://www-ljk.imag.fr/membres/Ludovic.Metivier/node2.html
• Nocedal himself has released his optimized implementation of L-BFGS-B version 3.0 open source: http://www.ece.northwestern.edu/~nocedal/lbfgsb.html
• There is also OpenTOAST from Germany: http://www.opentoast.de/ and http://www.opentoast.de/Forward-modelling-code_SPECFEM3D.php

[komatits]

I think at least the following groups have their own set of tools based in part on the ADJOINT_TOMOGRAPHY_TOOLS scripts, often combined with some calls to ObsPy: Vadim Monteiller and Yi Wang here in my group and in Sébastien Chevrot's group, people at Princeton, Elliott and Qinya at Toronto, Carl in Alaska, Hejun in Texas, Federica and Emanuele at INGV, maybe Min and Ying as well (the list is probably incomplete?).

We should also check with Ebru if she has another local version of the tools or if she uses the same tools as people at Princeton.

Vadim has also developed his own version of L-BFGS (he had done that before joining my group, while working on other projects without SPECFEM).

At some point let us try to talk about what is best to do in terms of trying to merge all of this in a single set of clean tools. We would all benefit from that.

[komatits]

From Dimitri, Sept 20, 2014:

I talked to Emanuele, Federica, Sébastien Chevrot, Vadim Monteiller, Yi Wang, Ebru, and to people in our group here, in particular Clément, Paul and Cédric. We are going to develop (more precisely: to finalize / combine / merge, since most of the bricks already exist) a set of common tools that will implement L-BFGS + the Wolfe condition of Nocedal (2006) that Vadim uses successfully.

These tools will run 100% in parallel and will avoid I/Os as much as possible by reading the input from disk only once at the beginning, then calling all the tools (summing, smoothing, preconditioning, model_update and so on) from the same parallel code, and writing the final result only once at the end

(reducing I/Os is critical on big machines, as observed e.g. on Titan in the US and on PRACE machines in Europe)

We will develop them in the current Git repository, in specfem3d/ADJOINT_TOMOGRAPHY_TOOLS.

Everybody is welcome to join the effort; the current situation (each group having its 'more or less OK but not so great' tools and scripts is not good; the sooner we switch to a better way of developing joint tools the better)

PS: Emanuele and Federica will probably organize a workshop of developers and users (of GEOCUBIT, SPECFEM, and ADJOINT_TOMOGRAPHY_TOOLS) in a couple of years. Excellent idea. I will be glad to attend and help organize it as well. But let us make sure the joint tools are ready long before that, ideally in the next few months. In my group Clément Durochat is in charge of that full time thus do not hesitate to contact him if you want (and cc me).

[komatits]

From Daniel Peter, Sept 20, 2014:

I just had some colleagues from Japan visiting, who were interested in doing tomography with both, global and local specfem3d versions. they would be happy to work with the tools available now, but would prefer to have identical ones for both packages. for that purpose i thought to create a new source directory, src/tomography/ in both packages, and have the corresponding tools (summing, preconditionig, smoothing, updates) compiled there (e.g. by 'make tomography'). As they prefer having the same tools for both packages, this would allow to easily switch between them. since the tools for the global version depend anyway on static compilation (values_from_mesher.h), i thought it would make sense to separate them and put them into the corresponding packages for now.

Of course in the long run it would be nicer to have a better framework built around and some more new/merged/combined features, but for now it would provide an easy quick start to them. just let me know what you think about.

[komatits]

From Ryan Modrak, Sept 22, 2014:

To follow up on your email, myself and others in Jeroen's group have been testing a python workflow package, which we are calling SEISFLOWS, capable of running various kinds of inversions and interfacing with SPECFEM2D, 3D, and 3D_Globe. In case anyone is interested, some more information about SEISFLOWS is given below. Please don't hesitate to contact me for more details.

In developing an adjoint tomography package, a major challenge for our group has been in trying to meet a rather diverse set of user needs. Exploration seismology, for example, deals with different objective functions, signal processing tools, inversion strategies, and file formats than regional and global tomography. Even within a purely regional and global context, differences between users in terms of programming preferences and underlying usage requirements have been a major reason behind the current fragmentation of our inversion tools.

To provide the flexibility needed to solve this problem, SEISFLOWS is designed to be highly modular, with separate 'workflow', 'system', 'solver', 'optimization', 'preprocessing', and 'postprocessing' components so that, for example, an NLCG optimization class can be exchanged for an LBFGS class, a SAC preprocessing module can be interchanged with an ADIOS module, or a PBS-TORQUE system interface can be swapped for a SLURM system interface. Crucial to this approach is the ability, offered by python, to derive new components from existing ones through inheritance and overloading. Thus, in case one needs to depart from the standard inversion workflow, one can draw selectively from the main code while continuing to work within and contribute to the overall package, without having to begin a separate line of development or introduce narrow or specialized modifications to the standard case.

Going back to the email thread, it sounds like the summing, smoothing, and preconditioning steps to be included in ADJOINT_TOMOGRAPHY_TOOLS fall within the scope of what have been calling postprocessing. I agree that refactoring and standardization of postprocessing tools across SPECFEM packages would be very helpful and work in this area would be greatly appreciated.

With respect to the larger nonlinear optimization worfklow, I was wondering, did you have a particular approach in mind? In my experience as someone looking at workflow issues, the combination of Fortran programs and shell scripts used in many SPECFEM applications seems somewhat too rigid for the task, but I would be eager to hear from anyone having ideas on the topic.

[komatits]

From Matthieu Lefebvre, Sept 22, 2014:

As Ryan mentioned, we are actively developing a couple of options to run seismic workflows.

Doing that we came across several challenges, mainly simplicity, scalability and portability.

Portability — As we are targeting diverse environments, it is a good option to rely on specific tools that have been developed with this aspect in mind. http://radical-cybertools.github.io/ is a good example of such a tool. It allows to transparently manage executions on local / remotes computers / clusters / grids as well as file transfers.

Scalability — Seismic workflows can involve a large number of processing tasks, some depending on the completion of others. Maximizing the workflow throughput is done in keeping tracks of the links between tasks in order to schedule them as soon as possible.

Simplicity — Simplicity along flexibility arises from a modular approach with several levels of granularity. Low level operations (make, setup, execution, …) are grouped into higher operations (forward, adjoint, update, …) that are in turn grouped into complete workflows. Standard operations and workflows are predefined and users can rearranges bits and pieces to suit particular needs.

[komatits]

From Jean-Paul Ampuero @jpampuero , Sept 22, 2014:

Will these tomography tools also include the moment tensor inversion tools? I am wondering if in the future some components of the tomography tools and workflows could be reused by my team for adjoint finite-source inversion. Somala et al implemented the linear source inversion problem with a conjugate gradient solver, but other algorithms would be needed for non-linear formulations of the problem and to impose positivity constraints.

S. N. Somala, J.-P. Ampuero and N. Lapusta Finite-fault source inversion using adjoint methods in 3D heterogeneous media Submitted to Geophys. J. Int. https://dl.dropboxusercontent.com/u/82498341/Somala%20et%20al%202013a%20GJI.pdf

[komatits]

From Dimitri, Sept 22, 2014:

Hi Jean-Paul,

Yes, they will. As you say that aspect is also very important.

BFGS for instance could be better than conjugate gradients.

[komatits]

From Dimitri, Sept 22, 2014:

Matthieu's two projects are developed at https://bitbucket.org/mpbl/flowjobs and https://bitbucket.org/mpbl/specfem3d_workflow

[komatits]

From Dimitri, Sept 22, 2014:

Hi Daniel,

Excellent idea I think. It will also make it much easier for Clément to make a parallel version of a calling tool that will call all of them with minimal I/Os.

Please let Clément know when you are done.

Let us call the new directory "tomography" rather than "tomo".

I am not sure what will remain in the "auxiliaries" directory after the move, if only stuff related to movies remains we could consider renaming "auxiliaries" to "movies" or "visualization" or similar ("auxiliaries" is indeed not very clear).

[komatits]

From Dimitri, Sept 22, 2014:

Dear Matthieu, Ryan and Daniel, dear all,

Excellent! Thanks for your three emails. I think that is exactly the tools we will need in the long run. Here is what I would suggest, following your emails:

1/ very good idea that Daniel moves the tomography tools (sum_kernels, smooth_kernels, model_updates etc) from the "auxiliaries" directory to a new directory called "tomography" (let us call it "tomography" rather than "tomo") for both packages (SPECFEM3D and SPECFEM3D_GLOBE). This will make it much clearer what they correspond to

2/ Clément can create a simple parallel version of Vadim's code for L-BFGS and the Wolfe condition; that tool will be very simple to use and very efficient (almost no I/Os)

3/ once ready, the great tools you are currently developing (SEISFLOWS, specfem3d_workflow, and FLOWJOBS) will be the best solution (combined in a single tool maybe?). They could call the parallel L-BFGS + Wolfe as one of their options. I think these tools will be really great and very efficient

4/ I have left a note in specfem3d/ADJOINT_TOMOGRAPHY_TOOLS/README to mention these new tools as well as the URL where they are developed.

[komatits]

See also https://github.com/geodynamics/specfem3d/issues/290

[komatits]

Done by Ryan and the Princeton group, here in our group Vadim Monteiller is also working on that and is almost done.

[komatits]

Reopened this issue, since the topic is currently being discussed by developers.

[komatits]

From Youyi Ruan, April 28, 2015:

The Seisstar website is built in order to host all the tools we are currently using for our global adjoint tomography and seismic imaging projects, and serves as an online library. We hope bringing all the resources together in one place will help other users to access these tools easily. Each package is hosted and maintained independently in their own Github repositories or websites.

Wenjie is in charge of the Seisstar website, so he and individual authors may know more about the current status of these packages.

I can speak part of these projects to my best knowledge, James, Lion and Matthieu are the main developers of ASDF toolboxes, they have finished the pyasdf package and implementing ASDF into SPECFEM3D.

Wenjie, Lion, Yanhua and myself are involved in the development of pre-processing workflows within a python framework, we have finished the pycmt3D (ready for beta release) and pyflex, and testing the multi-taper part of pyadjoint.

I think Seisstar is an open place and welcomes any contributors, and all we want is encouraging more people to try it.

[komatits]

From Ryan Modrak, April 29, 2015:

To add to Youyi's descriptions, let me try to summarize the Python inversion packages. Any additions, comments, or corrections are most welcome.

---

PySIT provides time- and frequency-domain waveform inversion capabilities. Rather than a set of preassembled workflows, PySIT is toolkit through which users can script their own workflows. Putting together migrations or inversions with PySIT is easy and flexible, in our experience.

As a flexible Python inversion toolkit, PySIT seems well suited to research and development. It is hard to judge how adaptable it is to production runs or to high performance environments. Because of an absence of portability/scalability features, it probably works best for small or medium scale inversions.

Forward and adjoint simulations are carried out with built-in 2D or 3D TDFD or FDFD solvers. Written in C with Python bindings, these solvers are very well integrated with other inversion components. Use of built-in solvers rather than external solvers avoids the need for complicated wrappers and makes the package easy to maintain, and hence robust.

Because of fundamental design issues, it is not possible, or at least not practical, to use PySIT and SPECFEM in the same workflow. This was the main factor behind our decision to develop SeisFlows.

PySIT was written by Russell Hewitt, who has a background in applied math and computer science. Through a mutual connection with Total, our group has visits and discussions with him.

---

LASIF is an interactive tool for processing data and launching simulations. Postprocessing and nonlinear optimization are not currently included. The work of a model update must be performed manually.

To a greater degree than other packages, LASIF is specialized for regional and global inversion. A major aspect of the interface is a fixed directory tree, which acts as a data structure containing the entire state of the inversion. With some tradeoffs when it comes to flexibility and possibility for abstraction, this approach has the advantage of being clear and simple, with directory names that make sense to regional and global seismologists.

As an interactive tool, adaptability to high performance environments is presumably limited. Because nonlinear optimization and workflow management capabilities are not currently included, LASIF cannot be used be as the basis for an automated inversion, though it is worth checking with Lion Krischner as new features may be in planning or development.

Besides developing LASIF, Lion Krischner is one of the main developers of obspy and our collaborator on data format issues. We have had extensive discussion with him about workflow issues, among other things.

---

SeisFlows is a Python waveform inversion package with both research and production capabilities, suitable for inexpensive 2D inversions through expensive 3D global 1 chunk inversions.

SeisFlows contains ready-to-go inversion and migration workflows, which can be adapted to individual user needs. To provide this flexibility, there is a plugin system through which users are offered choices in each of the following categories: workflow, system, solver, optimization, preprocessing, and postprocessing. In general, the tradeoff is a somewhat higher level of abstraction in exchange for greater flexibility and robustness.

While base classes are provided for all six categories above, special attention is paid to the optimization base class. Relative to NLCG, computational savings of roughly 50 percent can be realized with L-BFGS and an efficient line search in the majority of regional, global, and exploration problems. In the preprocessing base class, obspy, pyadjoint, and pyflex integration is currently in progress.

The easiest way to learn more about SeisFlows is to run some examples. If you have a Princeton account, simply follow the instructions in the SeisFlows documentation for running inversions on Tiger. If you do not have Princeton account, feel free to contact me to learn how to install and run it on your own cluster.

---

ASKI is described on its webpage as a waveform inversion package. At first glance though, it seems disconnected from mainstream waveform inversion.

Among other issues, the documentation states that, in carrying out a model update, a "sensitivity matrix is explicitly assembled" and a "linear system is solved using LAPACK libraries". In other words, model updates are computed using a full Newton/Gauss Newton method, which can be very costly.

[schumf]

Replying to the above comment (30 Apr 2015) on the Gauss Newton method of ASKI being very costly, I would like to mention the 2007 GJI paper by Chen et al.. They compare scattering-integral-based Gauss Newton FWI (which constitutes the methodology used by ASKI) with the adjoint conjugate-gradient method, with regard to computational expenses. They find that, in terms of the number of simulations to be done, scattering-integral-based inversion is more efficient than the adjoint conjugate-gradient method except for the case that the number of receiver components strongly exceeds the number of sources (threshold that can be deduced from the paper: number of receiver components being 23 to 41 times higher than the number of sources).

Obviously you pay this efficiency in the number of simulations by requiring to store complete wavefield information to hard disk. ASKI addresses this problem by using the frequency domain (storing spectral wavefields at only a small number of frequencies per iteration) and by pre-integration of the kernels onto an independent (in general much coarser) inversion grid that can be chosen according to the resolutional power of the inverted data. In our 2016 GJI paper on the methodology of ASKI we discuss these aspects in detail, give examples and argue that scattering-integral-based FWI thus becomes feasible today, even on regional scale.

[rmodrak]

I see that ASKI has been made to interface with both time- and frequency-domain forward solvers. While it is true that full Newton or Gauss Newton in the frequency domain is more efficient than in the time domain, it is simply not correct that full Newton or Gauss Newton is competitive with conjugate gradient or quasi-Newton methods either for seismic waveform inversion or general large-scale nonlinear optimization (see for example the chapter “Large-scale unconstrained optimization” in Nocedal and Wright or the review article “Numerical methods for large-scale nonlinear optimization” by Gould et al). One may be able to get around the huge cost in time, memory, or disk space of working with the full Hessian or Jacobian by projecting (preintegrating) onto a coarse basis, but only by working at a much coarser spacing than the resolution limit of seismic data at commonly used bandwidths.

It is problematic I think that Chen et al or Schumacher et al do not connect with the well-developed numerical optimization literature reviewed by Nocedal and Wright or Gould et al. Are you aware of truncated Newton methods? This is the approach actually taken by Akcelik et al and Epanomeritakis et al, whom you cite. The idea is to avoid full matrix assembly by working only with the action of the Hessian or Jacobian. Given than you submitted to GJI, I’m surprised, didn’t editors Metivier or Virieux have anything to say? Their work on applying truncated Newton methods to waveform inversion is becoming well known.

[schumf]

My motivation of giving my last comment was to make aware of the improvements we suggest in our GJI paper in order to make FWI based on waveform Fréchet kernels derived from Born scattering (in the context of Parker's "Geophysical Inverse Theory", 1994) more feasible for practical application, thus I intended to relax the very general statement "which can be very costly" at the end of the third last comment. In our paper we discuss in detail advantages as well as challenges arising from such a modified approach to what Chen et al. call "scattering-integral method". It is not apparent to us, why Chen et al.s' comparison of computational costs, to which we refer (Table 1 in their paper), is wrong.

It is certainly true that it depends on the specific seismic inversion problem whether the scattering-integral or the adjoint conjugate-gradient method, as compared by Chen et al., is more computationally expensive. ASKI simply offers an alternative approach, which eventually has to be chosen by the applicant, of course, depending on the problem to tackle.

[rmodrak]

At a glance, it is difficult to tell what exactly Chen et al are comparing, because they ignore the the standard terminology of numerical optimization, as well as the work of other seismologists, particularly those in exploration geophysics, who routinely use this terminology.

Also at glance, however, it is clear that the main work of ASKI does involve fully assembling and solving a linear system, which is very costly.

One way to put it is that any approach that involves fully assembling and solving a linear system is uncompetitive for large problems. In waveform inversion, the size of the model vector is on the order of millions to billions. The size of the Hessian is the square of that, so working with the fully assembled matrix is extremely inefficient.

More detailed calculations bear this out. A typical mesh in a regional, global or near surface seismology is roughly on the order of 500 x 500 x 500 grid points. (Since we are on a spectral element solver forum, I believe the coarsest mesh that handles surface wave propagation without artifacts is something like 128 x 128 elements at Earth’s surface, 640 x 640 integration points or so?) At this coarseness, you resolve periods as low as 20 s. With a dispersion condition for the mesh of about 5 grid points per wavelength, the resolution limit for waveform inversion is about twice the grid spacing. The memory or disk space requirements for such an inversion translate to mid terabytes to mid petabytes. By working with the action of the Hessian rather than the fully assembled matrix, this storage burden can be almost entirely avoided. Projecting to a coarser mesh may make full assembly feasible, but it provides no substitute for instead working with the action of the matrix operator rather than the assembled matrix.

On Wed, Mar 9, 2016 at 4:42 AM, schumf notifications@github.com wrote:

My motivation of giving my last comment https://github.com/geodynamics/specfem3d/issues/23#issuecomment-193656321 was to make aware of the improvements we suggest in our GJI paper in order to make FWI based on waveform Fréchet kernels derived from Born scattering (in the context of Parker's "Geophysical Inverse Theory", 1994) more feasible for practical application, thus I intended to relax the very general statement "which can be very costly" at the end of the third last comment https://github.com/geodynamics/specfem3d/issues/23#issuecomment-97622479. In our paper https://doi.org/10.1093/gji/ggv505 we discuss in detail advantages as well as challenges arising from such a modified approach to what Chen et al. call "scattering-integral method". It is not apparent to us, why Chen et al.s' https://doi.org/10.1111/j.1365-246X.2007.03429.x comparison of computational costs, to which we refer (Table 1 in their paper), is wrong.

It is certainly true that it depends on the specific seismic inversion problem whether the scattering-integral or the adjoint conjugate-gradient method, as compared by Chen et al., is more computationally expensive. ASKI simply offers an alternative approach, which eventually has to be chosen by the applicant, of course, depending on the problem to tackle.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/23#issuecomment-194211998 .

[schumf]

ASKI is not yet capable of global-scale FWI at high frequencies, and nobody ever claimed it was. Again, it certainly depends on the problems which method is favourable.

So far, we have applied ASKI e.g. to the Aegean, inverting data with up to 0.1 Hz frequency content (our 2016 GJI paper, section 4.2). Setting up and solving the full linear system (67.320 unknowns, 140.140 sensitivity equations, plus 134.640 sparse regularization equations) by a parallelized CG algorithm is negligible in terms of computation time (less than 0.1 percent) compared with solving the forward problem (regularly, not by path-specific approach) and computing the kernels.

ASKI, furthermore, does not perform coarsening just for the sake of computational feasibility, but instead aims at choosing an unstructured spatial model description that explicitly tries to exploit the locally varying resolving power of the inverted data, trying to avoid setting up insignificant parts of the sensitivity matrix, hence trying to reduce the Nullspace of the inverse problem. I suggest having a closer "glance" at our paper.

[rmodrak]

I looked closely at your package not long ago. I believe you should have a glance at the numerical optimization literature before claiming in the title of your paper that your package is computationally efficient.

On Thu, Mar 10, 2016 at 2:03 AM, schumf notifications@github.com wrote:

ASKI is not yet capable of global-scale FWI at high frequencies, and nobody ever claimed it was. Again, it certainly depends on the problems which method is favourable.

So far, we have applied ASKI e.g. to the Aegean, inverting data with up to 0.1 Hz frequency content (our 2016 GJI paper, section 4.2). Setting up and solving the full linear system (67.320 unknowns, 140.140 sensitivity equations, plus 134.640 sparse regularization equations) by a parallelized CG algorithm is negligible in terms of computation time (less than 0.1 percent) compared with solving the forward problem (regularly, not by path-specific approach) and computing the kernels.

ASKI, furthermore, does not perform coarsening just for the sake of computational feasibility, but instead aims at choosing an unstructured spatial model description that explicitly tries to exploit the locally varying resolving power of the inverted data, trying to avoid setting up insignificant parts of the sensitivity matrix, hence trying to reduce the Nullspace of the inverse problem. I suggest having a closer "glance" at our paper.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/23#issuecomment-194706810 .

[komatits]

Hi all,

If possible let us avoid making strong statements on the mailing list, let us try to keep the atmosphere and debate friendly. We should encourage people to contribute and support people who do.

Thanks, Best wishes,

Dimitri.

On 10/03/2016 20:36, rmodrak wrote:

I looked closely at your package not long ago. I believe you should have a glance at the numerical optimization literature before claiming in the title of your paper that your package is computationally efficient.

On Thu, Mar 10, 2016 at 2:03 AM, schumf notifications@github.com wrote:

ASKI is not yet capable of global-scale FWI at high frequencies, and nobody ever claimed it was. Again, it certainly depends on the problems which method is favourable.

So far, we have applied ASKI e.g. to the Aegean, inverting data with up to 0.1 Hz frequency content (our 2016 GJI paper, section 4.2). Setting up and solving the full linear system (67.320 unknowns, 140.140 sensitivity equations, plus 134.640 sparse regularization equations) by a parallelized CG algorithm is negligible in terms of computation time (less than 0.1 percent) compared with solving the forward problem (regularly, not by path-specific approach) and computing the kernels.

ASKI, furthermore, does not perform coarsening just for the sake of computational feasibility, but instead aims at choosing an unstructured spatial model description that explicitly tries to exploit the locally varying resolving power of the inverted data, trying to avoid setting up insignificant parts of the sensitivity matrix, hence trying to reduce the Nullspace of the inverse problem. I suggest having a closer "glance" at our paper.

— Reply to this email directly or view it on GitHub

https://github.com/geodynamics/specfem3d/issues/23#issuecomment-194706810 .

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/specfem3d/issues/23#issuecomment-195012675.

Dimitri Komatitsch CNRS Research Director (DR CNRS), Laboratory of Mechanics and Acoustics, UPR 7051, Marseille, France http://komatitsch.free.fr

[komatits]

See also https://github.com/geodynamics/specfem3d/issues/292