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,20 +271,19 @@ 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 if torch .cuda .is_available ():
278277 device = "cuda"
279278 elif torch .mps .is_available ():
280279 device = "mps"
281280 else :
282281 device = "cpu"
283-
282+
284283 self .transforms = [
285- DeferredLinearOperator ()
286- for _ in range (len (self .settings ))
284+ DeferredLinearOperator () for _ in range (len (self .settings ))
287285 ]
288-
286+
289287 D = X .shape [1 ]
290288
291289 freq_list = []
@@ -294,8 +292,7 @@ def __init__(
294292
295293 if len (s .bandwidths ) == 0 :
296294 full = np .zeros ((1 , D ), dtype = np .int32 )
297-
298-
295+
299296 local = np .atleast_2d (
300297 s .mode .index_set_without_zeros (
301298 np .array (s .bandwidths , dtype = np .int32 )
@@ -310,13 +307,9 @@ def __init__(
310307 freq_list .append (full )
311308
312309 freq = np .vstack (freq_list )
313-
310+
314311 X_torch = torch .tensor (X , dtype = torch .float64 , device = device )
315- I_torch = torch .tensor (
316- freq ,
317- dtype = torch .float64 ,
318- device = device
319- )
312+ I_torch = torch .tensor (freq , dtype = torch .float64 , device = device )
320313
321314 def trafo (fhat ):
322315 if self .system == "cos" or self .system == "cheb" :
@@ -325,8 +318,8 @@ def trafo(fhat):
325318
326319 kernel = 1.0
327320 for i in range (D ):
328- Xi = LazyTensor (X_torch [:, None , i : i + 1 ].contiguous ())
329- Ki = LazyTensor (I_torch [None , :, i : i + 1 ].contiguous ())
321+ Xi = LazyTensor (X_torch [:, None , i : i + 1 ].contiguous ())
322+ Ki = LazyTensor (I_torch [None , :, i : i + 1 ].contiguous ())
330323 kernel = kernel * (2 * torch .pi * Xi * Ki ).cos ()
331324
332325 fhat_torch = (
@@ -350,7 +343,9 @@ def trafo(fhat):
350343 two_pi_phase = - 2 * torch .pi * phase_fwd
351344 kernel = two_pi_phase .cos () + 1j * two_pi_phase .sin ()
352345
353- fhat_torch = torch .tensor (fhat , dtype = torch .complex128 , device = device )
346+ fhat_torch = torch .tensor (
347+ fhat , dtype = torch .complex128 , device = device
348+ )
354349 fhat_j = LazyTensor (fhat_torch [None , :, None ].contiguous ())
355350
356351 try :
@@ -365,15 +360,15 @@ def trafo(fhat):
365360 def adjoint (f ):
366361 if self .system == "cos" or self .system == "cheb" :
367362 mult = np .sqrt (2.0 ) ** np .count_nonzero (freq , axis = 1 )
368-
363+
369364 f_torch = torch .as_tensor (
370365 f , dtype = torch .float64 , device = device
371366 ).contiguous ()
372367
373368 kernel = 1.0
374369 for i in range (D ):
375- Ki = LazyTensor (I_torch [:, None , i : i + 1 ].contiguous ())
376- Xi = LazyTensor (X_torch [None , :, i : i + 1 ].contiguous ())
370+ Ki = LazyTensor (I_torch [:, None , i : i + 1 ].contiguous ())
371+ Xi = LazyTensor (X_torch [None , :, i : i + 1 ].contiguous ())
377372 kernel = kernel * (2 * torch .pi * Xi * Ki ).cos ()
378373
379374 f_i = LazyTensor (f_torch [None , :, None ].contiguous ())
@@ -387,7 +382,9 @@ def adjoint(f):
387382 print ("Error in KeOps adjoint:" , e )
388383 return None
389384 else :
390- f_torch = torch .tensor (f , dtype = torch .complex128 , device = device ).contiguous ()
385+ f_torch = torch .tensor (
386+ f , dtype = torch .complex128 , device = device
387+ ).contiguous ()
391388
392389 X_i = LazyTensor (X_torch [None , :, :].contiguous ())
393390 K_j = LazyTensor (I_torch [:, None , :].contiguous ())
@@ -409,10 +406,15 @@ def adjoint(f):
409406 return None
410407
411408 keops_dtype = np .float64 if self .system == "cos" else np .complex128
412- self .transforms = [DeferredLinearOperator (
413- dtype = keops_dtype , shape = (X .shape [0 ], len (freq )), mfunc = trafo , rmfunc = adjoint
414- )]
415-
409+ self .transforms = [
410+ DeferredLinearOperator (
411+ dtype = keops_dtype ,
412+ shape = (X .shape [0 ], len (freq )),
413+ mfunc = trafo ,
414+ rmfunc = adjoint ,
415+ )
416+ ]
417+
416418 else :
417419 self .transforms = []
418420 s1 = self .settings [0 ]
@@ -492,15 +494,15 @@ def adjoint_worker(i):
492494 adjoint_worker (i )
493495
494496 return fhat
495-
497+
496498 elif self .algorithm == "keops" :
497- return GroupedCoefficients (self .settings , self .transforms [0 ].H @ other )
498-
499+ return GroupedCoefficients (self .settings , self .transforms [0 ].H @ other )
500+
499501 elif self .algorithm == "direct" :
500502 return GroupedCoefficients (
501503 self .settings , (self .matrix .conj ()).T @ other
502504 )
503-
505+
504506 elif isinstance (other , GC ): # `f = F*fhat` (fhat = other)
505507 if self .settings != other .settings :
506508 raise ValueError (
@@ -528,10 +530,10 @@ def worker(i):
528530 worker (i )
529531
530532 return sum (results )
531-
533+
532534 elif self .algorithm == "keops" :
533535 return self .transforms [0 ] @ other .data
534-
536+
535537 elif self .algorithm == "direct" :
536538 return self .matrix @ other .data
537539 else :
0 commit comments