Implementing JUDI to work with custom datatype instead of SEGY

[kerim371]

Hi,

I'm trying to make JUDI work with custom data types.
As I know to process big data JUDI uses SegyIO and the latter has such objects as SeisCon (holds blocks) and SeisBlock (may retrieve data from SEGY).

I prepared the same objects (structs) called H5SeisCon and H5Seis (H5SeisCon's operator[] retrieves H5Seis and the latter's [] operator retrieves the trace data) and prepared functions to create judiVector from H5SeisCon. Thus judiVector{Float32, SeisCon} and judiVector{Float32, H5SeisCon} has pretty similar API. But the geometry is based on In Core Geometry struct.

Basically all my modifications are located in JUDI.jl/src/TimeModeling/h5geo/ dir and they don't overwrite anything within original JUDI API.

Now I'm trying to run the modeling_basic_2D.l example.
But I'm getting the error. The point of the error is that time_modeling receives the Vector{H5SeisCon} but it there should be H5SeisCon object. I can't understand which function is responsible to retrieve the H5SeisCon object from Vector before passing it to time_modeling function.

Also is it possible to disable distributed? It makes harder to debug the code

P.S.: I'm trying to implement this as JUDI works only with XY axes aligned models, seismic data and the model should be in the same length units and in the same coordinate reference system. I hope I could extend this with h5geo

[mloubout]

I can't understand which function is responsible to retrieve the H5SeisCon object from Vector before passing it to time_modeling function.

So it will depend on the operator. For the standard PDEFull operator, it will fo first through process_input_data then through subsample. THe second one is probably the one you need to add for your own HDF5 type:

https://github.com/kerim371/JUDI.jl/blob/e2bd96816592c8f5ffc7bd096ae4c0c54a63de21/src/TimeModeling/Utils/auxiliaryFunctions.jl#L755

You will also need to make sure it is then converted into an Array in time_modeling_serial in a similar way than SeisCon being read from disk.

Your addition look quite nice!!

[kerim371]

Thank you very much!

You will also need to make sure it is then converted into an Array in time_modeling_serial in a similar way than SeisCon being read from disk.

I think you are talking about places like this: if typeof(recData) == SegyIO.SeisCon then...

After I added subsample(a::Array{H5SeisCon,1}, i::Int64) = a[i] and overwritten the above line of code to something like if (typeof(recData) == SegyIO.SeisCon || typeof(recData) == H5SeisCon) ... then the modeling_basic_2D_example.jl ran fine.

How do you think is it possible to modify the above type checking typeof(recData) == SegyIO.SeisCon in JUDI library to something that maybe expanded to custom types like H5SeisCon?
It seems that it would be possible to define custom container like my H5SeisCon (and all the corresponding functions) outside the JUDI library and JUDI would work with it if typeof(recData) == SegyIO.SeisCon places were overwritten to take into account the idea that JUDI may work not only with SegyIO.SeisCon but also with external containers that privide the similar API.
In that case I would be able to not fork JUDI and explore JUDI without source code modification. It is easier then to keep my staff up to date with JUDI.

[mloubout]

I completely agree!! Making the API general so that you can define your own container and have JUDI take care of the rest would be perfect. And as you said much easier to maintain and work with on your or other user sides.

The main change in JUDI would be to create a "convert_to_array" of sort instead of that line and document it (in combination with subsample/process_input_data if needed) to have a nice "Build you own judiVector" WIKI. Your case would be a perfect example of it as well.

I'm happy to discuss the details on how to implement it and we'll be happy if you want to contribute to JUDI to help implement it.

[kerim371]

I would love to.

I can see two aspects:

1. it would be great to provide such 'abstract' API to not only judiVector but also generalize Out Of Core Geometry. I expect that when working with big data (hundred Gb or few Tb) In Core Geometry may take pretty much space on the RAM. I'm not sure but I believe JUDI GeometryOOC is more suited for large data than GeometryIC. In my case I added few members to GeometryOOC:

mutable struct GeometryOOC{T} <: Geometry{T}
    container::Array{SegyIO.SeisCon,1}
    dt::Array{T,1}
    nt::Array{Integer,1}
    t::Array{T,1}
    nsamples::Array{Integer,1}
    key::String
    segy_depth_key::String

    # h5geo members
    h5geo::PyCall.PyObject
    crs::String    # coordinate reference system
    lengthUnits::String
    temporalUnits::String
    model_origin_x::Float64
    model_origin_y::Float64
    model_orientation::Float64
end

and modified the function that reads data and coverts GeometryOOC to GeometryIC (function Geometry(geometry::GeometryOOC)) so that I transform units and coordinate system during the reading.

2. for now h5geo library doesn't provide any pip or conda installation and as it has several dependencies it is uncomfortable to built it from source for a usual user. Thus I think to faciliate our cooperation I should start with this task i.e. provide pip install h5geo. I have little experience of that and I have a feeling that it would take few days.

[mloubout]

That sounds like a good plan (point 2). I'll have a quick look at your code as well to see what can be used as inspiration to generalize the API.

Concerning GeometryOOC we have used it for large 3D dataset without issue and I don't think it has any stopper and is definitely much smaller than the IC version that stores all the corrdinates. But yes it would need to be generalized to allow different input format.

[kerim371]

@mloubout I'm back

It seems I prepared pip installation for h5geo project. I created git tag v0.1.0

To install it one should have C++ compiler supportind 17 standard available.
The command pip install git+https://github.com/tierra-colada/h5geo.git@v0.1.0 --verbose should build and install h5geo and all its dependencies (without GDAL yet) using CMake Superbuild technique. It usually takes few minutes to build and install it.

All dependencies will be stored in site-packages/h5geopy and in site-packages-h5gtpy (C++ wrapper to work with hdf5).
I tested it on Ubuntu 20.04 (compiler GCC 9.3) and on Windows 11 (MSVC 2019).

Now I think I'm going to start preparing Julia project that integrates h5geo with JUDI

[kerim371]

@mloubout hi,

I think I should define custom Geometry struct

mutable struct GeometryH5OOC{T} <: Geometry{T}
	h5geo::PyCall.PyObject			# retrieve coordinate ref system (CRS), length units, temporal units
	container::Array{H5SeisCon,1}	        # retrieve nSamp, dt, coordinates in the CRS and in units defined in 'h5geo' settings
	key::String					# source/receiver
	xkey:String					# x trace header name
	ykey:String					# y trace header name
	zkey:String					# z trace header name
	model_orientation::Float32		# if model is rotated around Z-axis then I should also rotate XY coordinates
end

Then in the same way JUDI loads GeometryIIC from GeometryOOC I will overload the function Geometry(geometry::GeometryH5OOC)

How do you think?

[mloubout]

Would work perfectly I think yes

[kerim371]

Just in case I think it should also be possible to slightly modify JUDI so that it was able to write data to custom formats.
For now I can see there are two functions that controls writing data: judiVector_to_SeisBlock and write_shot_record
Maybe it is worth to modify them (explicitely set the accepting types for arguments), is it possible?

[mloubout]

I have to think about this one, may need a bit of code re-write. But should be doable

[kerim371]

I just tested with custom H5GeometryOOC and it is working! Thank you very much for working on this!

The only thing is that get_computational_nt(H5GeometryOOC, ...) doesnt work:

ERROR: LoadError: type H5GeometryOOC has no field xloc
Stacktrace:
 [1] getproperty(x::H5GeometryOOC, f::Symbol)
   @ Base ./Base.jl:33
 [2] get_computational_nt(srcGeometry::H5GeometryOOC, recGeometry::H5GeometryOOC, model::Model; dt::Nothing)
   @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Utils/auxiliaryFunctions.jl:354
 [3] get_computational_nt(srcGeometry::H5GeometryOOC, recGeometry::H5GeometryOOC, model::Model)
   @ JUDI ~/.julia/dev/JUDI/src/TimeModeling/Utils/auxiliaryFunctions.jl:351

And one question related to perfomance issue: does perfomance of function n_samples(g::H5GeometryOOC, nsrc::Integer) is critical? For example for now to get n_samples I sum the size of N HDF5 datasets where N is equal to the number of sources.

[kerim371]

@mloubout I'm still about writing data in custom format.

What is the problem:
As I mentioned my model has some orientation and I read data from custom format. My OOC Geometry struct look:

mutable struct H5GeometryOOC <: Geometry{Float32}
	h5geo::PyCall.PyObject
	container::H5SeisCon    # I need to write the data in it
	key::String
	xkey::String
	ykey::String
	zkey::String
        do_coord_transform::Bool
        model_origin_x::Number
        model_origin_y::Number
	model_orientation::Number
end

At the time when I convert OOC Geometry to IC Geometry I rotate the coordinates.
But when JUDI writes data to segy function save_to_disk that calls write_shot_record accepts IC geometry:

function save_to_disk(shot::judiVector, srcGeometry::GeometryIC, srcData::Array, options::JUDIOptions, ::Val{true}) 
function write_shot_record(srcGeometry::GeometryIC, srcData, recGeometry::GeometryIC, recData, options)

thus even if I overwrite these functions I'm still unable to get all information needed (model_orientation, container) to correctly write the data.

As I can see save_to_disk is called only from time_modeling function

Maybe we could change the save_to_disk and write_shot_record so that it accept GeomOrNot or Geometry types?
Then in time_modeling instead of overwriting recGeometry and srcGeometry introduce new variables recGeometryIC and srcGeometryIC and pass recGeometry and srcGeometry to save data functions.

How do you think?
If that is possible then I could try it on my PC and if I succeded I will create PR.