11import logging
2- import warnings
32
43import numpy as np
54
6- from pylops import LinearOperator
75from pylops .basicoperators import Diagonal , MatrixMult , Restriction , Transpose
6+ from pylops .LinearOperator import aslinearoperator
87from pylops .utils ._internal import _value_or_list_like_to_tuple
98from pylops .utils .backend import get_array_module
109
@@ -58,7 +57,7 @@ def _linearinterp(dims, iava, axis=-1, dtype="float64"):
5857 ) + Diagonal (weights , dims = dimsd , axis = axis , dtype = dtype ) * Restriction (
5958 dims , iva_r , axis = axis , dtype = dtype
6059 )
61- return Op , iava
60+ return Op , iava , dims , dimsd
6261
6362
6463def _sincinterp (dims , iava , axis = 0 , dtype = "float64" ):
@@ -76,18 +75,18 @@ def _sincinterp(dims, iava, axis=0, dtype="float64"):
7675 otherdims = ncp .array (dims )
7776 otherdims = ncp .roll (otherdims , - axis )
7877 otherdims = otherdims [1 :]
79- Op = MatrixMult (sinc , dims = otherdims , dtype = dtype )
78+ Op = MatrixMult (sinc , otherdims = otherdims , dtype = dtype )
8079
8180 # create Transpose operator that brings axis to first dimension
81+ dimsd = list (dims )
82+ dimsd [axis ] = len (iava )
8283 if axis > 0 :
8384 axes = np .arange (len (dims ), dtype = int )
8485 axes = np .roll (axes , - axis )
85- dimsd = list (dims )
86- dimsd [axis ] = len (iava )
8786 Top = Transpose (dims , axes = axes , dtype = dtype )
8887 T1op = Transpose (dimsd , axes = axes , dtype = dtype )
8988 Op = T1op .H * Op * Top
90- return Op
89+ return Op , dims , dimsd
9190
9291
9392def Interp (dims , iava , axis = - 1 , kind = "linear" , dtype = "float64" ):
@@ -193,9 +192,15 @@ def Interp(dims, iava, axis=-1, kind="linear", dtype="float64"):
193192 if kind == "nearest" :
194193 interpop , iava = _nearestinterp (dims , iava , axis = axis , dtype = dtype )
195194 elif kind == "linear" :
196- interpop , iava = _linearinterp (dims , iava , axis = axis , dtype = dtype )
195+ interpop , iava , dims , dimsd = _linearinterp (dims , iava , axis = axis , dtype = dtype )
197196 elif kind == "sinc" :
198- interpop = _sincinterp (dims , iava , axis = axis , dtype = dtype )
197+ interpop , dims , dimsd = _sincinterp (dims , iava , axis = axis , dtype = dtype )
199198 else :
200199 raise NotImplementedError ("kind is not correct..." )
200+ # add dims and dimsd to composite operators (not needed for neareast as
201+ # interpop is a Restriction operator already
202+ if kind != "nearest" :
203+ interpop = aslinearoperator (interpop )
204+ interpop .dims = dims
205+ interpop .dimsd = dimsd
201206 return interpop , iava
0 commit comments