Pre-stack modeling operator for multiple wavelets

未关闭
#648 2 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 2 天内回复

还没有人认领这个 Issue。

评估

难度
4/5
预计耗时
3-5 天
新手友好度
30/100
Issue 类型
缺陷
描述清晰度
基本清楚
活跃度
停滞
技术栈
python
领域
backend

调研方向

Start at the pre-stack inversion operator’s explicit=False path and read how it constructs Convolve1D operators and combines them with BlockDiag; compare this with the dense multi-wavelet path shown in the issue. First focus on the stated 1D case, leaving 2D and 3D extension for later. Done means sparse results for multiple wavelets are consistent with the dense implementation.

由索引模型根据 Issue 内容生成。

描述

enhancement

🛠 Multi-Wavelet Support in Pre-Stack Inversion Operator (explicit=False)

📖 Description

The convolution operator used by pylops for pre-stack inversion (explicit=False) is implemented using the Conv1D class. Currently, pylops supports only a single wavelet for all angles, meaning the inversion assumes the same wavelet for each angle.

I attempted to modify the code to allow a different wavelet for each angle (mid-angle in my case, as I’m using partial stacks). However, my current results are only consistent when using dense operators (explicit=True), whereas I would expect them to also work in the sparse case.


🔍 Observations

  • ✅ The dense implementation (explicit=True) behaves as expected.
  • ❌ The sparse implementation (explicit=False) does not produce comparable results.
  • ⚠️ The issue seems related to how the Conv1D operator is built when handling multiple wavelets.

The code using Conv1D for multiple wavelets is only working for spatdims=1 at the moment, the idea is to extend for 2D or 3D cases when i adress this issue.


💡 Code Snippet

✅ Dense Implementation (Working)
# Handling multiple wavelets
if isinstance(wav, list):
    print("Dense multi-wavelet operator")
    n = len(wav)
    C = np.zeros((n * nt0, n * nt0))

    wsub_list = []
    for w in wav:
        wsub = np.asarray(convmtx(w, nt0, len(w) // 2)[:nt0])
        wsub_list.append(wsub)

    C = get_block_diag(theta)(*wsub_list)

else:
    # Create wavelet operator
    C = ncp.asarray(convmtx(wav, nt0, len(wav) // 2)[:nt0])
    C = [C] * ntheta
    C = get_block_diag(theta)(*C)

# Combine operators
M = ncp.dot(C, ncp.dot(G, D))
Preop = MatrixMult(M, otherdims=spatdims, dtype=dtype)

❌ Sparse implementation (Not-working)
    else:
        if isinstance(wav, list):
            print("Dense multi-wavelet operator")
            # Criar uma lista de operadores de convolução
            Cop_list = []

            for w in wav:
                Cop = Convolve1D(
                    (nt0),
                    h=w,
                    offset=len(w) // 2,  # Define o centro da wavelet
                    axis=0,
                    dtype='float32',
                    method='fft'
                )
                Cop_list.append(Cop.tosparse())  # Adiciona à lista de operadores'
            

            Cop = BlockDiag(Cop_list)

        else:
            print("Sparse single wavelet operator")
            # Create wavelet operator
            Cop = Convolve1D(
                dims,
                h=wav,
                offset=len(wav) // 2,
                axis=0,
                dtype=dtype,
            )
主要语言
Python
星标
539
派生
126
平均合并
2 天 17 小时
30 天内合并 PR
21

环境准备

  • 没有 Dockerfile 或 Docker Compose 文件
  • 没有 Pull Request 模板
  • 阅读贡献指南

从这里开始

  1. 先读完整个 Issue,再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库,在一个分支上完成修改。
  4. 提交 Pull Request,并在描述里引用这个 Issue 编号。

PyLops/pylops 的其他 Issue

查看 PyLops/pylops 的全部 Issue

相似的 Issue

更多 Python Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。