src.acoustools.Solvers

  1from acoustools.Utilities import *
  2from acoustools.Optimise.Constraints import constrain_phase_only
  3from acoustools.Constraints import constrain_amplitude, constrain_field, constrain_field_weighted
  4from acoustools.Optimise.Objectives import target_gorkov_BEM_mse_sine_objective
  5from acoustools.Optimise.Constraints import sine_amplitude
  6import torch
  7
  8from torch import Tensor
  9from types import FunctionType
 10
 11from vedo import Mesh
 12
 13
 14
 15def wgs_solver_unbatched(A, b, K):
 16    '''
 17    @private
 18    unbatched WGS solver for transducer phases, better to use `wgs_solver_batch` \\
 19    `A` Forward model matrix to use \\ 
 20    `b` initial guess - normally use `torch.ones(N,1).to(device)+0j`\\
 21    `k` number of iterations to run for \\
 22    returns (hologram image, point phases, hologram)
 23    '''
 24    #Written by Giorgos Christopoulos 2022
 25    AT = torch.conj(A).T.to(device)
 26    b0 = b.to(device)
 27    x = torch.ones(A.shape[1],1).to(device) + 0j
 28    for kk in range(K):
 29        y = torch.matmul(A,x)                                   # forward propagate
 30        y = y/torch.max(torch.abs(y))                           # normalize forward propagated field (useful for next step's division)
 31        b = torch.multiply(b0,torch.divide(b,torch.abs(y)))     # update target - current target over normalized field
 32        b = b/torch.max(torch.abs(b))                           # normalize target
 33        p = torch.multiply(b,torch.divide(y,torch.abs(y)))      # keep phase, apply target amplitude
 34        r = torch.matmul(AT,p)                                  # backward propagate
 35        x = torch.divide(r,torch.abs(r))                        # keep phase for hologram  
 36                    
 37    return y, p, x
 38
 39def wgs_solver_batch(A, b, iterations):
 40    '''
 41    @private
 42    batched WGS solver for transducer phases\\
 43    `A` Forward model matrix to use \\ 
 44    `b` initial guess - normally use `torch.ones(self.N,1).to(device)+0j`\\
 45    `iterations` number of iterations to run for \\
 46    returns (point pressure ,point phases, hologram)
 47    '''
 48    AT = torch.conj(A).mT.to(device).to(DTYPE)
 49    
 50    b = b.expand(A.shape[0],-1,-1)
 51    b0 = b.to(device).expand(A.shape[0],-1,-1).to(DTYPE)
 52    
 53    x = torch.ones(A.shape[2],1).to(device).to(DTYPE) + 0j
 54    for kk in range(iterations):
 55        p = A@x
 56        p,b = constrain_field_weighted(p,b0,b)
 57        x = AT@p
 58        x = constrain_amplitude(x)
 59    y =  torch.abs(A@x) 
 60    return y, p, x
 61
 62def wgs(points:Tensor,iterations:int = 200, board:Tensor|None = None, A:Tensor|None = None, b:Tensor|None=None, 
 63        return_components:bool=False, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None) -> Tensor:
 64    '''
 65    Weighted-GS algorithm\n
 66    :param points: Points to use
 67    :param iter: Number of iterations for WGS, default:`200`
 68    :param board: The Transducer array, If `None` uses default two 16x16 arrays
 69    :param A: Forward model matrix to use 
 70    :param b: initial guess - If none will use `torch.ones(N,1).to(device)+0j`
 71    :param return_components: IF True will return `hologram image, point phases, hologram` else will return `hologram`, default False
 72    :return: hologram
 73
 74    ```Python
 75    from acoustools.Solvers import wgs
 76    from acoustools.Utilities import create_points, propagate_abs
 77
 78    p = create_points(2,4)
 79    print(p)
 80    x = wgs(p)
 81    print(propagate_abs(x,p))
 82    
 83    p = p.squeeze(0)
 84    print(p)
 85    x = wgs(p)
 86    print(propagate_abs(x,p))
 87    ```
 88    '''
 89    if board is None:
 90        board = TRANSDUCERS
 91
 92    if len(points.shape) > 2:
 93        N = points.shape[2]
 94        batch=True
 95    else:
 96        N = points.shape[1]
 97        batch=False
 98
 99    if A is None:
100        A = forward_model(points, board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms).to(DTYPE)
101    if b is None:
102        b = torch.ones(N,1).to(device).to(DTYPE)+0j
103
104    if batch:
105        img,pha,act = wgs_solver_batch(A,b,iterations)
106    else:
107        img,pha,act = wgs_solver_unbatched(A,b,iterations)
108
109    if return_components:
110        return img,pha,act
111    return act
112
113
114def gspat_solver(R,forward, backward, target, iterations,return_components=False):
115    '''
116    @private
117    GS-PAT Solver for transducer phases\\
118    `R` R Matrix\\
119    `forward` forward propagation matrix\\
120    `backward` backward propagation matrix\\
121    `target` initial guess - can use `torch.ones(N,1).to(device)+0j`
122    `iterations` Number of iterations to use\\
123    returns (hologram, point activations)
124    '''
125    #Written by Giorgos Christopoulos 2022
126    field = target 
127
128    for _ in range(iterations):
129        
130#     amplitude constraint, keeps phase imposes desired amplitude ratio among points     
131        target_field = constrain_field(field, target)
132#     backward and forward propagation at once
133        field = torch.matmul(R,target_field)
134#     AFTER THE LOOP
135#     impose amplitude constraint and keep phase, after the iterative part this step is different following Dieg
136    target_field = torch.multiply(target**2,torch.divide(field,torch.abs(field)**2))
137#     back propagate 
138    complex_hologram = torch.matmul(backward,target_field)
139#     keep phase 
140    phase_hologram = torch.divide(complex_hologram,torch.abs(complex_hologram))
141    if return_components:
142        points = torch.matmul(forward,phase_hologram)
143
144        return phase_hologram, points
145    else:
146        return phase_hologram
147
148
149def gspat(points:Tensor|None=None, board:Tensor|None=None,A:Tensor|None=None,B:Tensor|None=None, 
150          R:Tensor|None=None ,b:Tensor|None = None, iterations:int=200, return_components:bool=False, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None) -> Tensor:
151    '''
152    GSPAT Solver\n
153    :param points: Target point positions
154    :param board: The Transducer array, if None uses default two 16x16 arrays
155    :param A: The Forward propagation matrix, if `None` will be computed 
156    :param B: The backwards propagation matrix, if `None` will be computed 
157    :param R: The R propagation matrix, if `None` will be computed 
158    :param b: initial guess - If None will use `torch.ones(N,1).to(device)+0j`
159    :param iterations: Number of iterations to use
160    :param return_components: IF True will return `hologram, pressure` else will return `hologram`, default True
161    :return: Hologram
162    '''
163
164    if board is None:
165        board = TRANSDUCERS
166
167    if A is None:
168        A = forward_model(points,board,p_ref=p_ref,norms=norms, k=k, transducer_radius=transducer_radius)
169    if B is None:
170        B = torch.conj(A).mT.to(DTYPE)
171    if R is None:
172        R = A@B
173
174    if b is None:
175        if is_batched_points(points):
176            b = torch.ones(points.shape[2],1).to(device).to(DTYPE)
177        else:
178            b = torch.ones(points.shape[1],1).to(device).to(DTYPE)
179    
180    
181    if return_components:
182        phase_hologram,pres = gspat_solver(R,A,B,b, iterations,return_components)
183        return phase_hologram,pres
184    
185    phase_hologram = gspat_solver(R,A,B,b, iterations,return_components)
186    return phase_hologram
187
188
189def naive_solver_batched(points,board=TRANSDUCERS, activation=None, A=None, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None):
190    '''
191    @private
192    Batched naive (backpropagation) algorithm for phase retrieval\\
193    `points` Target point positions\\
194    `board` The Transducer array, default two 16x16 arrays\\
195    `activation` Initial starting point activation \\
196    returns (point activations, hologram)
197    '''
198    if activation is None:
199        activation = torch.ones(points.shape[2],1, device=device, dtype=DTYPE) +0j
200    
201    if A is None:
202        A = forward_model_batched(points,board,p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
203    back = torch.conj(A).mT
204    trans = back@activation
205    trans_phase=  constrain_amplitude(trans)
206    out = A@trans_phase
207
208
209    return out, trans_phase
210
211def naive_solver_unbatched(points,board=TRANSDUCERS, activation=None,A=None, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None):
212    '''
213    @private
214    Unbatched naive (backpropagation) algorithm for phase retrieval\\
215    `points` Target point positions\\
216    `board` The Transducer array, default two 16x16 arrays\\
217    `activation` Initial starting point activation \\
218    returns (point activations, hologram)
219    '''
220    if activation is None:
221        activation = torch.ones(points.shape[1]) +0j
222        activation = activation.to(device)
223    if A is None:
224        A = forward_model(points,board,p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
225    back = torch.conj(A).T
226    trans = back@activation
227    trans_phase=  constrain_amplitude(trans)
228    out = A@trans_phase
229
230
231    return out, trans_phase
232
233def naive(points:Tensor, board:Tensor|None = None, return_components:bool=False, activation:Tensor|None=None, A=None, 
234          p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None, iterations=None) -> Tensor:
235    '''
236    Naive solver\n
237    :param points: Target point positions
238    :param board: The Transducer array, default two 16x16 arrays
239    :param return_components: If `True` will return `hologram, pr
240    :param iterations: Ignroed - for compatability
241    :param A: propagator to use
242    :return: hologram
243    '''
244    if board is None:
245        board = TRANSDUCERS
246    if is_batched_points(points):
247        out,act = naive_solver_batched(points,board=board, activation=activation, A=A, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
248    else:
249        out,act = naive_solver_unbatched(points,board=board, activation=activation, A=A, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
250    if return_components:
251        return act, out
252    return act
253
254def ph_thresh(z_last,z,threshold):
255    '''
256    @private
257    Phase threshhold between two timesteps point phases, clamps phase changes above `threshold` to be `threshold`\\
258    `z_last` point activation at timestep t-1\\
259    `z` point activation at timestep t\\
260    `threshold` maximum allowed phase change\\
261    returns constrained point activations
262    '''
263
264    ph1 = torch.angle(z_last)
265    ph2 = torch.angle(z)
266    dph = ph2 - ph1
267    
268    dph = torch.atan2(torch.sin(dph),torch.cos(dph)) 
269    
270    dph[dph>threshold] = threshold
271    dph[dph<-1*threshold] = -1*threshold
272    
273    ph2 = ph1 + dph
274    z = abs(z)*torch.exp(1j*ph2)
275    
276    return z
277
278def soft(x,threshold):
279    '''
280    @private
281    Soft threshold for a set of phase changes, will return the change - threshold if change > threshold else 0\\
282    `x` phase changes\\
283    `threshold` Maximum allowed hologram phase change\\
284    returns new phase changes
285    '''
286    y = torch.max(torch.abs(x) - threshold,0).values
287    y = y * torch.sign(x)
288    return y
289
290def ph_soft(x_last,x,threshold):
291    '''
292    @private
293    Soft thresholding for holograms \\
294    `x_last` Hologram from timestep t-1\\
295    `x` Hologram from timestep t \\
296    `threshold` Maximum allowed phase change\\
297    returns constrained hologram
298    '''
299    pi = torch.pi
300    ph1 = torch.angle(x_last)
301    ph2 = torch.angle(x)
302    dph = ph2 - ph1
303
304    dph[dph>pi] = dph[dph>pi] - 2*pi
305    dph[dph<-1*pi] = dph[dph<-1*pi] + 2*pi
306
307    dph = soft(dph,threshold)
308    ph2 = ph1 + dph
309    x = abs(x)*torch.exp(1j*ph2)
310    return x
311
312def temporal_wgs(A:Tensor, y:Tensor, K:int,ref_in:Tensor, ref_out:Tensor,T_in:float,T_out:float) -> Tensor:
313    '''
314    Based off `
315    Giorgos Christopoulos, Lei Gao, Diego Martinez Plasencia, Marta Betcke, 
316    Ryuji Hirayama, and Sriram Subramanian. 2023. 
317    Temporal acoustic point holography. (2024) https://doi.org/10.1145/3641519.3657443 \n
318    WGS solver for hologram where the phase change between frames is constrained \n
319    :param A: Forward model  to use
320    :param y: initial guess to use normally use `torch.ones(self.N,1).to(device)+0j`
321    :param K: Number of iterations to use
322    :param ref_in: Previous timesteps hologram
323    :param ref_out: Previous timesteps point activations
324    :param T_in: Hologram phase change threshold
325    :param T_out: Point activations phase change threshold
326    :return: (hologram image, point phases, hologram)
327    
328
329    '''
330    #ref_out -> points
331    #ref_in-> transducers
332    AT = torch.conj(A).mT.to(device)
333    y0 = y.to(device)
334    x = torch.ones(A.shape[2],1).to(device) + 0j
335    for kk in range(K):
336        z = torch.matmul(A,x)                                   # forward propagate
337        z = z/torch.max(torch.abs(z),dim=1,keepdim=True).values # normalize forward propagated field (useful for next step's division)
338        z = ph_thresh(ref_out,z,T_out); 
339        
340        y = torch.multiply(y0,torch.divide(y,torch.abs(z)))     # update target - current target over normalized field
341        y = y/torch.max(torch.abs(y),dim=1,keepdim=True).values # normalize target
342        p = torch.multiply(y,torch.divide(z,torch.abs(z)))      # keep phase, apply target amplitude
343        r = torch.matmul(AT,p)                                  # backward propagate
344        x = torch.divide(r,torch.abs(r))                        # keep phase for hologram    
345        x = ph_thresh(ref_in,x,T_in);    
346    return y, p, x
347
348
349
350
351
352
353def gradient_descent_solver(points: Tensor, objective: FunctionType, board:Tensor|None=None, optimiser:torch.optim.Optimizer=torch.optim.Adam, lr: float=0.01, 
354                            objective_params:dict={}, start:Tensor|None=None, iters:int=200, 
355                            maximise:bool=False, targets:Tensor=None, constrains:FunctionType=constrain_phase_only, log:bool=False, return_loss:bool=False,
356                            scheduler:torch.optim.lr_scheduler.LRScheduler=None, scheduler_args:dict=None, save_each_n:int = 0, save_set_n:list[int] = None,
357                            init_type:Literal['rand', 'ones','focal','trap']|Tensor='rand',
358                            p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None) -> Tensor:
359    '''
360    Solves phases using gradient descent\n
361    :param points: Target point positions 
362    :param objective: Objective function - must take have an input of (`transducer_phases, points, board, targets, **objective_params`), `targets` may be `None` for unsupervised
363    :param board: The Transducer array, default two 16x16 arrays
364    :param optimiser: Optimiser to use (should be compatable with the interface from from `torch.optim`). Default `torch.optim.Adam`
365    :param lr: Learning Rate to use. Default `0.01`
366    :param objective_params: Any parameters to be passed to `objective` as a dictionary of `{parameter_name:parameter_value}` pairs. Default `{}`
367    :param start: Initial guess. If None will default to a random initilsation of phases 
368    :param iters: Number of optimisation Iterations. Default 200
369    :param maximise: Set to `True` to maximise the objective, else minimise. Default `False`
370    :param targets: Targets to optimise towards for supervised optimisation, unsupervised if set to `None`. Default `None`
371    :param constrains: Constraints to apply to result 
372    :param log: If `True` prints the objective values at each step. Default `False`
373    :param return_loss: If `True` save and return objective values for each step as well as the optimised result 
374    :param scheduler: Learning rate scheduler to use, if `None` no scheduler is used. Default `None` 
375    :param scheduler_args: Parameters to pass to `scheduler`
376    :param save_each_n: For n>0 will save the optimiser results at every n steps. Set either `save_each_n` or `save_set_iters`
377    :param save_set_iters: List containing exact iterations to save optimiser results at. Set either `save_each_n` or `save_set_iters`
378    :param init_type: type of initialisation to use. rand:random, ones:tensor of 1, focal:naive focal point, trap:two focal points offset in the z-axis
379    :return: optimised result and optionally the objective values and results (see `return_loss`, `save_each_n` and `save_set_iters`). If either are returned both will be returned but maybe empty if not asked for
380    
381    ```Python
382    from acoustools.Optimise.Objectives import propagate_abs_sum_objective
383    from acoustools.Solvers import gradient_descent_solver
384    from acoustools.Optimise.Constraints import constrain_phase_only
385
386    p = create_points(4,2)
387    x = gradient_descent_solver(p,propagate_abs_sum_objective, 
388                                maximise=True, constrains=constrain_phase_only, 
389                                log=False, lr=1e-1)
390
391    print(propagate_abs(x,p))
392
393    ```
394    ''' 
395
396
397    if board is None:
398        board = TRANSDUCERS
399
400    losses = []
401    results = {}
402    B = points.shape[0] if points is not None else 1
403    N = points.shape[2] if points is not None else 1
404    M = board.shape[0]
405    if start is None:
406        # start = torch.ones((B,M,1)).to(device) +0j
407        if type(init_type) == Tensor:
408            start = init_type
409        else: 
410            if init_type == 'ones':
411                start = torch.ones((B,M,1))
412            elif init_type == 'focal':
413    
414                start = naive(points, board=board,return_components=False, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
415            elif init_type == 'trap':
416                new_points = points.expand(B,3,2*N).clone()
417                SCALE = 2
418                new_points[:,2,:N] -= Constants.wavelength / SCALE
419                new_points[:,2,N:] += Constants.wavelength / SCALE
420                target_phases = torch.zeros(B,2*N)
421                # target_phases[:,N:] = Constants.pi
422                activation = torch.exp(1j * target_phases).unsqueeze(2).to(device)
423                
424            
425                start = naive(new_points, board, return_components=False, activation=activation, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
426                
427            else: #rand is default
428                start = torch.exp(1j*torch.rand((B,M,1))*torch.pi)
429
430        start=start.to(device).to(DTYPE)
431    
432    
433    # param = torch.nn.Parameter(start).to(device)
434    param = start.requires_grad_()
435    optim = optimiser([param],lr)
436    if scheduler is not None:
437        scheduler = scheduler(optim,**scheduler_args)
438
439
440    for epoch in range(iters):
441        optim.zero_grad()       
442
443        loss = objective(param, points, board, targets, **objective_params)
444
445        if log:
446            print(epoch, loss.data.item())
447
448        if maximise:
449            loss *= -1
450        
451        if return_loss:
452            losses.append(loss)
453        
454        if save_each_n > 0 and epoch % save_each_n == 0:
455            results[epoch] = param.clone().detach()
456        elif save_set_n is not None and epoch in save_set_n:
457            results[epoch] = param.clone().detach()
458
459        
460        loss.backward(torch.tensor([1]*B).to(device))
461        optim.step()
462        if scheduler is not None:
463            scheduler.step()
464        
465        if constrains is not None:
466            param.data = constrains(param)
467
468    if return_loss or save_each_n > 0:
469        return param, losses, results
470    
471
472    return param
473    
474
475def iterative_backpropagation(points:Tensor,iterations:int = 200, board:Tensor|None = None, A:Tensor|None = None, 
476                              b:Tensor|None=None, return_components: bool=False, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, norms=None) -> list[Tensor]:
477    '''
478    IB solver for transducer phases\n
479    :param points: Points to use
480    :param iterations: Number of iterations for WGS, default`200`
481    :param board: The Transducer array, default two 16x16 arrays
482    :param A: Forward model matrix to use 
483    :param b: initial guess - If none will use `torch.ones(N,1).to(device)+0j`
484    :param return_components: IF True will return `hologram image, point phases, hologram` else will return `hologram`, default False
485    :return: (point pressure ,point phases, hologram)
486
487    ```Python
488    from acoustools.Solvers import iterative_backpropagation
489    from acoustools.Utilities import create_points, propagate_abs
490
491    p = create_points(2,1)
492    print(p)
493    x = iterative_backpropagation(p)
494    print(propagate_abs(x,p))
495    
496    p = p.squeeze(0)
497    print(p)
498    x = iterative_backpropagation(p)
499    print(propagate_abs(x,p))
500
501    ```
502    '''
503
504    if board is None:
505        board  = TRANSDUCERS
506    
507    if len(points.shape) > 2:
508        N = points.shape[2]
509        batch=True
510    else:
511        N = points.shape[1]
512        batch=False
513
514    if A is None:
515        A = forward_model(points, board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=norms)
516    
517    if batch:
518        M = A.shape[2]
519    else:
520        M = A.shape[1]
521
522
523    if b is None:
524        b = torch.ones(N,1).to(device).to(DTYPE) +0j
525    
526    AT =  torch.conj(A).mT.to(device)
527    x = torch.ones(M,1).to(device).to(DTYPE) 
528    for kk in range(iterations):
529        p = A@x
530        p = constrain_field(p,b)
531        x = AT@p
532        x = constrain_amplitude(x)
533    y =  torch.abs(A@x) 
534    if return_components:
535        return y, p, x
536    else:
537        return x
538    
539
540
541def gorkov_target(points:Tensor, objective:FunctionType = target_gorkov_BEM_mse_sine_objective,
542                  board:Tensor=None, U_targets:Tensor=None, iterations:int=100, lr:int=1e9,
543                  constraint:FunctionType=sine_amplitude, reflector:Mesh|None=None, path:str|None=None) -> Tensor:
544    '''
545    Phase retrieval to generate target gorkov values at points via `acoustools.Solvers.gradient_descent_solver` \n
546    :param points: points of interest
547    :param objective: Objevtive function to minimise, default `acoustools.Optimise.Objectives.target_gorkov_BEM_mse_sine_objective`
548    :param board: Board to use. Default `acoustools.Utilities.TOP_BOARD`
549    :param U_targets: Target Gorkov values
550    :param iterations: Iterations to use for the solver
551    :param lr: learning rate
552    :param constraint: constraint function to use in the optimiser. default `acoustools.Optimise.Constraints.sine_amplitude`
553    :param reflector: Mesh to use as reflector or None
554    :param path: BEM path
555    :returns hologram:
556    '''
557
558    if board is None:
559        board = TOP_BOARD
560
561    if type(U_targets) == float:
562        U_targets = torch.tensor([U_targets])
563    
564    if U_targets is None:
565        U_targets = torch.tensor([-1e-5])
566
567    x = gradient_descent_solver(points, objective, 
568                                board, log=True, targets=U_targets, iters=iterations, 
569                                lr=lr, constrains=constraint, objective_params={'root':path,'reflector':reflector})
570    
571    return x
572
573
574def kd_solver(points:Tensor, board:Tensor|None = None,k=c.k):
575    '''
576    Solver for one point by setting phases to be equal at target point \\
577    see `A volumetric display for visual, tactile and audio presentation using acoustic trapping`\\
578    :param points: point to use - must be only one point
579    :param board: Transducers, if None will use two 16x16 arrays
580    :param k: wavenumber 
581    '''
582    B = points.shape[0]
583    M = board.shape[0]
584
585    b = board.unsqueeze(0).permute((0,2,1))
586    p = points.expand(B,3,M)
587
588    distance = torch.sqrt(torch.sum((p - b)**2,dim=1)).unsqueeze_(0).mT
589    distance = distance.to(device).to(DTYPE)
590    return torch.exp(1j * -1* distance*k)
591
592
593def translate_hologram(x:Tensor,board:Tensor|None=None, dx:float=0, dy:float=0, dz:float=0, k:float=c.k):
594
595    '''
596    Translates an existing hologram by (dx,dy,dz)\\
597    :param x: Hologram
598    :param board: Transducers, if None will use two 16x16 arrays
599    :param dx: x translation
600    :param dy: y translation
601    :param dz: z translation
602    :param k: wavenumber 
603    '''
604    
605    if board is None:
606        board = TRANSDUCERS
607
608    p2 =  create_points(1,1,dx,dy,dz)
609
610    B = p2.shape[0]
611    M = board.shape[0]
612
613    b = board.unsqueeze(0).permute((0,2,1))
614    p2 = p2.expand(B,3,M)
615
616    distance_p2 = torch.sqrt(torch.sum((p2 - b)**2,dim=1)).unsqueeze_(0).mT.to(device).to(DTYPE)
617    distance_p = torch.sqrt(torch.sum((b)**2,dim=1)).unsqueeze_(0).mT.to(device).to(DTYPE)
618
619    distance = distance_p-distance_p2
620
621    kd = (k*distance)
622
623    phase = torch.angle(x)
624    phase_2 = phase + kd
625    x2 = torch.exp(1j * phase_2)
626
627    return x2