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,9 @@ def __init__(
287285
288286 if len (s .bandwidths ) == 0 :
289287 full = np .zeros ((1 , D ), dtype = np .int32 )
290-
288+
291289 elif self .system == "cos" :
292-
290+
293291 local = np .atleast_2d (
294292 s .mode .index_set_without_zeros (
295293 np .array (s .bandwidths , dtype = np .int32 )
@@ -300,7 +298,7 @@ def __init__(
300298
301299 for i , dim in enumerate (s .u ):
302300 full [:, dim ] = local [:, i ]
303-
301+
304302 else :
305303 local = index_set_without_zeros (
306304 np .array (s .bandwidths , dtype = np .int32 )
@@ -314,13 +312,9 @@ def __init__(
314312 freq_list .append (full )
315313
316314 freq = np .vstack (freq_list )
317-
315+
318316 X_torch = torch .tensor (X , dtype = torch .float64 , device = "cuda" )
319- I_torch = torch .tensor (
320- freq ,
321- dtype = torch .float64 ,
322- device = "cuda"
323- )
317+ I_torch = torch .tensor (freq , dtype = torch .float64 , device = "cuda" )
324318
325319 def trafo (fhat ):
326320 if self .system == "cos" or self .system == "cheb" :
@@ -329,8 +323,8 @@ def trafo(fhat):
329323
330324 kernel = 1.0
331325 for i in range (D ):
332- Xi = LazyTensor (X_torch [:, None , i : i + 1 ].contiguous ())
333- Ki = LazyTensor (I_torch [None , :, i : i + 1 ].contiguous ())
326+ Xi = LazyTensor (X_torch [:, None , i : i + 1 ].contiguous ())
327+ Ki = LazyTensor (I_torch [None , :, i : i + 1 ].contiguous ())
334328 kernel = kernel * (2 * torch .pi * Xi * Ki ).cos ()
335329
336330 fhat_torch = (
@@ -354,7 +348,9 @@ def trafo(fhat):
354348 two_pi_phase = - 2 * torch .pi * phase_fwd
355349 kernel = two_pi_phase .cos () + 1j * two_pi_phase .sin ()
356350
357- fhat_torch = torch .tensor (fhat , dtype = torch .complex128 , device = "cuda" )
351+ fhat_torch = torch .tensor (
352+ fhat , dtype = torch .complex128 , device = "cuda"
353+ )
358354 fhat_j = LazyTensor (fhat_torch [None , :, None ].contiguous ())
359355
360356 try :
@@ -370,16 +366,16 @@ def adjoint(f):
370366 if self .system == "cos" or self .system == "cheb" :
371367 mult = np .sqrt (2.0 ) ** np .count_nonzero (freq , axis = 1 )
372368 mult_torch = torch .tensor (mult , dtype = torch .float64 , device = "cuda" )
373-
369+
374370 f_torch = (
375371 torch .as_tensor (f , dtype = torch .float64 , device = "cuda" )
376372 * mult_torch
377373 ).contiguous ()
378374
379375 kernel = 1.0
380376 for i in range (D ):
381- Ki = LazyTensor (I_torch [:, None , i : i + 1 ].contiguous ())
382- Xi = LazyTensor (X_torch [None , :, i : i + 1 ].contiguous ())
377+ Ki = LazyTensor (I_torch [:, None , i : i + 1 ].contiguous ())
378+ Xi = LazyTensor (X_torch [None , :, i : i + 1 ].contiguous ())
383379 kernel = kernel * (2 * torch .pi * Xi * Ki ).cos ()
384380
385381 f_i = LazyTensor (f_torch [None , :, None ].contiguous ())
@@ -393,7 +389,9 @@ def adjoint(f):
393389 print ("Error in KeOps adjoint:" , e )
394390 return None
395391 else :
396- f_torch = torch .tensor (f , dtype = torch .complex128 , device = "cuda" ).contiguous ()
392+ f_torch = torch .tensor (
393+ f , dtype = torch .complex128 , device = "cuda"
394+ ).contiguous ()
397395
398396 X_i = LazyTensor (X_torch [None , :, :].contiguous ())
399397 K_j = LazyTensor (I_torch [:, None , :].contiguous ())
@@ -415,10 +413,15 @@ def adjoint(f):
415413 return None
416414
417415 keops_dtype = np .float64 if self .system == "cos" else np .complex128
418- self .transforms = [DeferredLinearOperator (
419- dtype = keops_dtype , shape = (X .shape [0 ], len (freq )), mfunc = trafo , rmfunc = adjoint
420- )]
421-
416+ self .transforms = [
417+ DeferredLinearOperator (
418+ dtype = keops_dtype ,
419+ shape = (X .shape [0 ], len (freq )),
420+ mfunc = trafo ,
421+ rmfunc = adjoint ,
422+ )
423+ ]
424+
422425 else :
423426 self .transforms = []
424427 s1 = self .settings [0 ]
@@ -498,15 +501,15 @@ def adjoint_worker(i):
498501 adjoint_worker (i )
499502
500503 return fhat
501-
504+
502505 elif self .algorithm == "keops" :
503- return GroupedCoefficients (self .settings , self .transforms [0 ].H @ other )
504-
506+ return GroupedCoefficients (self .settings , self .transforms [0 ].H @ other )
507+
505508 elif self .algorithm == "direct" :
506509 return GroupedCoefficients (
507510 self .settings , (self .matrix .conj ()).T @ other
508511 )
509-
512+
510513 elif isinstance (other , GC ): # `f = F*fhat` (fhat = other)
511514 if self .settings != other .settings :
512515 raise ValueError (
@@ -534,10 +537,10 @@ def worker(i):
534537 worker (i )
535538
536539 return sum (results )
537-
540+
538541 elif self .algorithm == "keops" :
539542 return self .transforms [0 ] @ other .data
540-
543+
541544 elif self .algorithm == "direct" :
542545 return self .matrix @ other .data
543546 else :
0 commit comments