search

slimgroup / JUDI.jl

Julia Devito inversion.
https://slimgroup.github.io/JUDI.jl
MIT License
96 stars 29 forks source link

Added Viking Graben line12 example #181

Closed kerim371 closed 9 months ago

kerim371 commented 1 year ago

Viking Grabel line 12 is 2D seismic marine field data.

The added example contains preprocessing in Madagascar and FWI/RTM. LSRTM is in the process

codecov[bot] commented 1 year ago

Codecov Report

Attention: 232 lines in your changes are missing coverage. Please review.

Comparison is base (783c142) 81.70% compared to head (b5421fc) 66.29%. Report is 56 commits behind head on master.

Files Patch % Lines
...imeModeling/Preconditioners/DataPreconditioners.jl 0.00% 47 Missing :warning:
src/TimeModeling/Types/ModelStructure.jl 77.10% 38 Missing :warning:
src/TimeModeling/Types/GeometryStructure.jl 59.75% 33 Missing :warning:
src/TimeModeling/Utils/auxiliaryFunctions.jl 78.22% 27 Missing :warning:
src/TimeModeling/LinearOperators/callable.jl 35.48% 20 Missing :warning:
src/TimeModeling/Types/judiComposites.jl 81.66% 11 Missing :warning:
...meModeling/Preconditioners/ModelPreconditioners.jl 0.00% 10 Missing :warning:
src/JUDI.jl 76.31% 9 Missing :warning:
src/TimeModeling/Modeling/propagation.jl 33.33% 8 Missing :warning:
src/TimeModeling/LinearOperators/lazy.jl 50.00% 6 Missing :warning:
... and 8 more
Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #181 +/- ## =========================================== - Coverage 81.70% 66.29% -15.41% =========================================== Files 28 28 Lines 2192 2433 +241 =========================================== - Hits 1791 1613 -178 - Misses 401 820 +419 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

mloubout commented 1 year ago

Thanks a lot! Will look at it when back from holiday

rafaelorozco commented 1 year ago

IM a student at SLIM lab. this PR and the PDF are so cool and informative. Thank you for sharing!

kerim371 commented 1 year ago

@mloubout generally I understood what and how to do. In the nearest days I will try to fix that.

kerim371 commented 1 year ago

@mloubout hi,

I'm trying to use judiFilter in that way:

Ml_freq = judiFilter(d_obs.geometry, 0.01, frq*1000f0)
Mr_freq = judiFilter(q.geometry, 0.01, frq*1000f0)

fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_freq[indsrc]*d_obs[indsrc], options=jopt)

but I get exception:

ERROR: MethodError: no method matching copy(::SeisCon)
Closest candidates are:
  copy(::Union{MathOptInterface.DualPowerCone, MathOptInterface.EqualTo, MathOptInterface.GreaterThan, MathOptInterface.Interval, MathOptInterface.LessThan, MathOptInterface.Parameter, MathOptInterface.PowerCone, MathOptInterface.Semicontinuous, MathOptInterface.Semiinteger}) at /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/MathOptInterface/4g9vU/src/sets.jl:2584
  copy(::Union{TransparentColor{C, T, N} where N, C} where {T, C<:Union{AbstractRGB{T}, AbstractGray{T}}}) at /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/ColorVectorSpace/QI5vM/src/ColorVectorSpace.jl:250
  copy(::DataStructures.Accumulator) at /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/DataStructures/59MD0/src/accumulator.jl:41
  ...
Stacktrace:
  [1] *(n::Float32, s::SeisCon)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Types/judiVector.jl:156
  [2] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [3] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [4] getindex
    @ ./broadcast.jl:575 [inlined]
  [5] copy
    @ ./broadcast.jl:922 [inlined]
  [6] materialize
    @ ./broadcast.jl:883 [inlined]
  [7] zero(::Type{Float32}, v::judiVector{Float32, SeisCon}; nsrc::Int64)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Types/judiVector.jl:137
  [8] zero(::Type{Float32}, v::judiVector{Float32, SeisCon})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Types/judiVector.jl:137
  [9] copy
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Types/abstract.jl:36 [inlined]
 [10] deepcopy
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Types/abstract.jl:37 [inlined]
 [11] matvec(#unused#::FrequencyFilter{Float32, 0.01f0, 5.0f0}, Din::judiVector{Float32, SeisCon})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Preconditioners/DataPreconditioners.jl:158
 [12] *(J::FrequencyFilter{Float32, 0.01f0, 5.0f0}, ms::judiVector{Float32, SeisCon})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Preconditioners/base.jl:28
 [13] top-level scope
    @ REPL[2]:1

It seems SeisCon is not supported. Probably SeisBlock is supported but I prefer not to load all the SEGY to RAM with SeisBlock. May we fix that?

Also does judiFilter works with Hz or with kHz? Is it possible to use LowPass filter instead of BandPass (I prefer not to touch low frequencies that may suffer from slope of BandPass filter I think).

JUDI v3.2.3

kerim371 commented 1 year ago

@mloubout hi,

Should I wait for judiFilter exception to be fixed so I could continue working on the PR?

mloubout commented 1 year ago

Should I wait for judiFilter exception to be fixed so I could continue working on the PR?

You can do F*get_data(dobs) to convert it to in-core for the time being. For the frequencies, just put something very low for fmin (like fmin=.001) to do low-pass

mloubout commented 1 year ago

Also not sure how much of it you can use for the processing but I made an initial Julia wrapper for RSF/Madagascar

https://github.com/mloubout/Madagascar.jl

kerim371 commented 1 year ago

Tested with FWI and RTM should work too.

mloubout commented 1 year ago

Ok I think it's mostly there. Can you

The final step would be to add some of the notebooks/md files to the documentation. But I can do that myself once this is merged

kerim371 commented 1 year ago

@mloubout hi,

It seems I cleaned it all up.

It seems github is able visualize matplotlib images without problems. But plotly wasn't working that is why in deghosting notebook we could not see some pictures. I modified proc notebooks so they work with matplotlib only.

Most likely the same problem we encounter when looking at initial_model notebook which uses Julia's Plots library for plotting. I guess if I use PyCall's matplotlib it will work fine. How do you think should I try?

  • I still think should only keep the Madagascar notebook used for this example, and you can add a links to that other repo for more information.

If you don't mind I suggest to ask opinion/permission the author of that repo @rmorel.

Basically two steps from preprocessing (geometry check and deghosting) of Viking Graben line 12 are based on @rmorel work.

If I simply use this package then it will be unsufficient. If we recommend the user to preprocess the data using original repo and then come back to this processing flow then it will be confusing. Moreover the original repo is a little old and may not work properly with newer version of pandas.

I think the processing should not be divided in two packages. So I suggest either leave it as is or improve the original repo probably with my additions.

rmorel commented 1 year ago

@kerim371, all the code and work in my repository is open source (GPL). Feel free to reuse or modify the code in any way you want.

I'm happy to hear that it's useful =).

mloubout commented 1 year ago

I think the processing should not be divided in two packages.

Sorry wasn't what I mean. I only meant to remove the notebooks that do processing steps not sued for the FWI/RTM so that this example is self-contained.

Thanks for the plots. Sorry I'm a bit busy these days but not forgetting about this great addition.

Would you mind If I pushed some changes to this PR myself. I won't change any of the actual content but I want to set it up with DrWatson.jl so that it is safely fully reproducible at least for the julia side of it.

kerim371 commented 1 year ago

Sorry wasn't what I mean. I only meant to remove the notebooks that do processing steps not sued for the FWI/RTM so that this example is self-contained.

If I understood you then processing notebooks are required for RTM/LSRTM. FWI doesn't require these notebooks. So you mean completely remove these processing notebooks from PR (for example create another repository on github with these notebooks) or restructure the folders so the user could understand that FWI doesn't require processing?

Would you mind If I pushed some changes to this PR myself. I won't change any of the actual content but I want to set it up with DrWatson.jl so that it is safely fully reproducible at least for the julia side of it.

Yes of course.

Also there are few things that I'm working on:

  1. try to replace julia's Plots by PyPlot in initial_model notebook so that the images were visible on github
  2. use ImageGathers for QC and further try to perform AVO
  3. add LSRTM script
  4. replace get_data(d_obs[indsrc]) by d_obs[indsrc] when SegyIO PR is merged

I don't know when I finish these so you can merge it whenever you want. If I succeeded then I will reopen/crete new PR with these changes.

mloubout commented 9 months ago

@kerim371 would you mind checking if this still runs properly with the new changes. I set up this project with DrWatson.jl so that all packages used are fixed at a given version making it fully reproducible.

kerim371 commented 9 months ago

@mloubout I made some small fixes.

I haven't worked with DrWatson until now and this is how I tried to use this field example:

  1. ] dev JUDI
  2. follow instructions in README.md like: preinstall Julia deps (along with DrWatson.jl)
  3. run script fw_lbfgs.jl

If that is the way how users are going to use it now then it should work fine.