11import numpy as np
2+ import pykeops
23import torch
4+ from pykeops .torch import LazyTensor
35
46from pyGroupedTransforms import *
5- from .NFFTtools import index_set_without_zeros
6-
7- import pykeops
8- from pykeops .torch import LazyTensor
97
8+ from .NFFTtools import index_set_without_zeros
109
1110# All code that is linked to NFMTtools or to system = "mixed" is not tested yet....
1211
@@ -272,13 +271,12 @@ def __init__(
272271 bandwidths = s .bandwidths , X = np .copy (X [:, u ], order = "C" )
273272 )
274273 elif algorithm == "keops" :
275- self .matrix = np .empty ((0 ,0 ), dtype = object )
276-
274+ self .matrix = np .empty ((0 , 0 ), dtype = object )
275+
277276 self .transforms = [
278- DeferredLinearOperator ()
279- for _ in range (len (self .settings ))
277+ DeferredLinearOperator () for _ in range (len (self .settings ))
280278 ]
281-
279+
282280 D = X .shape [1 ]
283281
284282 freq_list = []
@@ -287,9 +285,8 @@ def __init__(
287285
288286 if len (s .bandwidths ) == 0 :
289287 full = np .zeros ((1 , D ), dtype = np .int32 )
290-
291-
292- else :
288+
289+ else :
293290 local = index_set_without_zeros (
294291 np .array (s .bandwidths , dtype = np .int32 )
295292 ).T
@@ -302,15 +299,11 @@ def __init__(
302299 freq_list .append (full )
303300
304301 freq = np .vstack (freq_list )
305-
302+
306303 X_torch = torch .tensor (X , dtype = torch .float64 , device = "cuda" )
307- I_torch = torch .tensor (
308- freq ,
309- dtype = torch .float64 ,
310- device = "cuda"
311- )
312-
313- def trafo (fhat ):
304+ I_torch = torch .tensor (freq , dtype = torch .float64 , device = "cuda" )
305+
306+ def trafo (fhat ):
314307 X_i = LazyTensor (X_torch [:, None , :].contiguous ())
315308 K_j = LazyTensor (I_torch [None , :, :].contiguous ())
316309 phase_fwd = (X_i * K_j ).sum (- 1 )
@@ -327,32 +320,41 @@ def trafo(fhat):
327320 except Exception as e :
328321 print ("Error in KeOps trafo:" , e )
329322 return None
330-
323+
331324 def adjoint (f ):
332- f_torch = torch .tensor (f , dtype = torch .complex128 , device = "cuda" ).contiguous ()
333-
325+ f_torch = torch .tensor (
326+ f , dtype = torch .complex128 , device = "cuda"
327+ ).contiguous ()
328+
334329 X_i = LazyTensor (X_torch [None , :, :].contiguous ())
335330 K_j = LazyTensor (I_torch [:, None , :].contiguous ())
336331
337332 phase = (K_j * X_i ).sum (- 1 )
338- kernel = (- 2 * torch .pi * phase ).cos () - 1j * (- 2 * torch .pi * phase ).sin ()
333+ kernel = (- 2 * torch .pi * phase ).cos () - 1j * (
334+ - 2 * torch .pi * phase
335+ ).sin ()
339336
340337 f_i = LazyTensor (f_torch [None , :, None ].contiguous ())
341338
342339 try :
343340 result = (kernel * f_i ).sum (dim = 1 )
344341 torch .cuda .synchronize ()
345342 result = result .squeeze ().cpu ().numpy ()
346-
343+
347344 return result
348345 except Exception as e :
349346 print ("Error in KeOps adjoint:" , e )
350347 return None
351348
352- self .transforms = [DeferredLinearOperator (
353- dtype = np .complex128 , shape = (X .shape [0 ], len (freq )), mfunc = trafo , rmfunc = adjoint
354- )]
355-
349+ self .transforms = [
350+ DeferredLinearOperator (
351+ dtype = np .complex128 ,
352+ shape = (X .shape [0 ], len (freq )),
353+ mfunc = trafo ,
354+ rmfunc = adjoint ,
355+ )
356+ ]
357+
356358 else :
357359 self .transforms = []
358360 s1 = self .settings [0 ]
@@ -432,15 +434,15 @@ def adjoint_worker(i):
432434 adjoint_worker (i )
433435
434436 return fhat
435-
437+
436438 elif self .algorithm == "keops" :
437- return GroupedCoefficients (self .settings , self .transforms [0 ].H @ other )
438-
439+ return GroupedCoefficients (self .settings , self .transforms [0 ].H @ other )
440+
439441 elif self .algorithm == "direct" :
440442 return GroupedCoefficients (
441443 self .settings , (self .matrix .conj ()).T @ other
442444 )
443-
445+
444446 elif isinstance (other , GC ): # `f = F*fhat` (fhat = other)
445447 if self .settings != other .settings :
446448 raise ValueError (
@@ -468,10 +470,10 @@ def worker(i):
468470 worker (i )
469471
470472 return sum (results )
471-
473+
472474 elif self .algorithm == "keops" :
473475 return self .transforms [0 ] @ other .data
474-
476+
475477 elif self .algorithm == "direct" :
476478 return self .matrix @ other .data
477479 else :
0 commit comments