uafgeotools / mtuq

moment tensor uncertainty quantification
BSD 2-Clause "Simplified" License
65 stars 22 forks source link

CMA-ES implementation and Force example #236

Closed thurinj closed 8 months ago

thurinj commented 8 months ago

Here is a brief list of the main contributions of this pull request:

  1. Implementation of the CMA-ES (wiki) algorithm in MTUQ. This comes with its own set of utils and plotting. It integrates most of the pre-existing MTUQ classes and functions so as not to alter the main part of the code.
  2. Refinement of the matplotlib backend. I added Lune and Force Sphere visualization options based on matplotlib. The only piece missing to be 100% competitive with GMT is the beachball plots (which would be quite challenging to implement). Nevertheless, it is easier to work with to implement new options (plot type, scaling, combining plots, etc.). This is the default backend for most of the CMA-ES plotting functions, although the CMA-ES output is also compatible with the GMT backend. This also includes some default plots for polarity misfits (both on the Lune and Force Sphere).
  3. Added the Barry arm landslide data for the Force source type of MTUQ. This comes with bundled data and two examples for the Gridsearch and CMA-ES mode.
  4. Implemented point force polarity misfit. As it appears that the DPRK nuclear test of 2017-09-03 is very well fitted by a point force, we can also use the data polarity with a point force source. Implemented from R. Madariaga book chapter on Seismic Source Theory. (PDF available here at the time of this PR).

I have made a lot of tests and integrated feedback from @SeismoFelix and @ammcpherson who helped test an earlier version of the CMA-ES implementation.

The CMA-ES solver has been tested with several data groups and polarity misfits. It includes the option to balance the contribution of each waveform group that is given as input (and will otherwise default to normalize each misfit contribution by the data norm).

It is initialized with the following:

parameter_list = initialize_mt(Mw_range=[4,6], src_type=src_type)

# Creating list of important objects to be passed to solving and plotting functions later
DATA = [data_bw, data_sw]
MISFIT = [misfit_bw, misfit_sw]
PROCESS = [process_bw, process_sw]
GREENS = [greens_bw, greens_sw]

popsize = 48 # -- CMA-ES population size - number of mutants
CMA = CMA_ES(parameter_list , origin=origin, lmbda=popsize, event_id=event_id)
iter = 60

which uses most of the default MTUQ objects. The solve method will perform 60 CMA-ES iterations, with a population of 48 samples/mutants to evolve, and plot waveform misfit and misfit map every 10 steps along the way. The misfit map is updated by the new one rather than being kept on disk.

CMA.Solve(DATA, stations, MISFIT, PROCESS, GREENS, iter, plot_interval=10, misfit_weights=[2., 3.])

If you want to be more specific with your result processing, it is possible to return the whole CMA-ES history of samples with the following:

result = CMA.mutants_logger_list

which returns a similar object to the random grid search solutions (mtuq.grid_search.MTUQDataFrame) and is thus compatible with other base functions in the code.

(The fork was open for PR earlier today but an mpi4py error prevented the tests to pass and it was discarded for this updated version. I added a check in code that should help with the tests.)

rmodrak commented 8 months ago

Julien, It appears the request includes CMA-ES (very useful) as well as many other changes that I'm just starting to go through. Thanks for your patience as myself and others I work with get into the details. Looking forward to discussing at the next AFRL project meeting. -Ryan

thurinj commented 8 months ago

Will cancel this PR and restructure.