src.acoustools.BEM.Gradients.H_Gradient
1import torch 2from torch import Tensor 3 4from vedo import Mesh 5 6import hashlib, pickle 7 8from acoustools.Utilities import device, DTYPE, forward_model_batched, forward_model_grad, forward_model_second_derivative_unmixed 9from acoustools.BEM.Forward_models import compute_bs, compute_A, get_cache_or_compute_H 10from acoustools.Mesh import get_centres_as_points, get_normals_as_points, board_name 11import acoustools.Constants as Constants 12 13from acoustools.BEM.Gradients.E_Gradient import get_G_partial 14 15 16 17def grad_H(points: Tensor, scatterer: Mesh, transducers: Tensor, return_components:bool = False, 18 path:str='', H:Tensor=None, use_cache_H:bool=True) ->tuple[Tensor,Tensor, Tensor] | tuple[Tensor,Tensor, Tensor, Tensor,Tensor, Tensor, Tensor]: 19 ''' 20 Computes the gradient of H wrt scatterer centres\n 21 Ignores `points` - for compatability with other gradient functions, takes centres of the scatterers 22 :param scatterer: The mesh used (as a `vedo` `mesh` object) 23 :param transducers: Transducers to use 24 :param return_components: if true will return the subparts used to compute the derivative 25 :return grad_H: The gradient of the H matrix wrt the position of the mesh 26 ''' 27 28 if H is None: 29 H = get_cache_or_compute_H(scatterer, transducers, use_cache_H, path) 30 31 32 # centres = torch.tensor(scatterer.cell_centers().points).to(device).T.unsqueeze_(0) 33 centres = get_centres_as_points(scatterer) 34 35 36 M = centres.shape[2] 37 38 B = compute_bs(scatterer,transducers) 39 A = compute_A(scatterer) 40 A_inv = torch.inverse(A).to(DTYPE) 41 42 43 Bx, By, Bz = forward_model_grad(centres, transducers) 44 Bx = Bx.to(DTYPE) 45 By = By.to(DTYPE) 46 Bz = Bz.to(DTYPE) 47 48 49 Ax, Ay, Az = get_G_partial(centres,scatterer,transducers) 50 # Ax *= -1 51 # Ay *= -1 52 # Az *= -1 53 54 Ax = (-1* Ax) 55 Ay = (-1* Ay) 56 Az = (-1* Az) 57 58 59 60 eye = torch.eye(M).to(bool) 61 Ax[:,eye] = 0 62 Ay[:,eye] = 0 63 Az[:,eye] = 0 64 65 # A_inv_x = (-1*A_inv @ Ax @ A_inv).to(DTYPE) 66 # A_inv_y = (-1*A_inv @ Ay @ A_inv).to(DTYPE) 67 # A_inv_z = (-1*A_inv @ Az @ A_inv).to(DTYPE) 68 69 # Hx_old = (A_inv_x@B) + (A_inv@Bx) 70 # Hy_old = (A_inv_y@B) + (A_inv@By) 71 # Hz_old = (A_inv_z@B) + (A_inv@Bz) 72 73 74 Hx = A_inv @ (Bx - Ax @ H) 75 Hy = A_inv @ (By - Ay @ H) 76 Hz = A_inv @ (Bz - Az @ H) 77 78 79 Hx = Hx.to(DTYPE) 80 Hy = Hy.to(DTYPE) 81 Hz = Hz.to(DTYPE) 82 83 if return_components: 84 return Hx, Hy, Hz, A, A_inv, Ax, Ay, Az 85 else: 86 return Hx, Hy, Hz 87 88 89def grad_2_H(points: Tensor, scatterer: Mesh, transducers: Tensor, A:Tensor|None = None, 90 A_inv:Tensor|None = None, Ax:Tensor|None = None, Ay:Tensor|None = None, Az:Tensor|None = None) -> Tensor: 91 ''' 92 Computes the second derivative of H wrt scatterer centres\n 93 Ignores `points` - for compatability with other gradient functions, takes centres of the scatterers 94 :param scatterer: The mesh used (as a `vedo` `mesh` object) 95 :param transducers: Transducers to use 96 :param A: The result of a call to `compute_A` 97 :param A_inv: The inverse of `A` 98 :param Ax: The gradient of A wrt the x position of scatterer centres 99 :param Ay: The gradient of A wrt the y position of scatterer centres 100 :param Az: The gradient of A wrt the z position of scatterer centres 101 :return Haa: second order unmixed gradient of H wrt scatterer positions 102 ''' 103 104 centres = get_centres_as_points(scatterer) 105 M = centres.shape[2] 106 107 B = compute_bs(scatterer,transducers) 108 109 Fx, Fy, Fz = forward_model_grad(centres, transducers) 110 Fx = Fx.to(DTYPE) 111 Fy = Fy.to(DTYPE) 112 Fz = Fz.to(DTYPE) 113 Fa = torch.stack([Fx,Fy,Fz],dim=3) 114 115 Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(centres, transducers) 116 Faa = torch.stack([Fxx,Fyy,Fzz],dim=3) 117 118 F = forward_model_batched(centres, transducers) 119 120 if A is None: 121 A = compute_A(scatterer) 122 123 if A_inv is None: 124 A_inv = torch.inverse(A) 125 126 if Ax is None or Ay is None or Az is None: 127 Ax, Ay, Az = get_G_partial(centres,scatterer,transducers) 128 eye = torch.eye(M).to(bool) 129 Ax[:,eye] = 0 130 Ay[:,eye] = 0 131 Az[:,eye] = 0 132 Ax = Ax.to(DTYPE) 133 Ay = Ay.to(DTYPE) 134 Az = Az.to(DTYPE) 135 Aa = torch.stack([Ax,Ay,Az],dim=3) 136 137 138 A_inv_x = (-1*A_inv @ Ax @ A_inv).to(DTYPE) 139 A_inv_y = (-1*A_inv @ Ay @ A_inv).to(DTYPE) 140 A_inv_z = (-1*A_inv @ Az @ A_inv).to(DTYPE) 141 142 143 A_inv_a = torch.stack([A_inv_x,A_inv_y,A_inv_z],dim=3) 144 145 m = centres.permute(0,2,1) 146 m = m.expand((M,M,3)) 147 148 m_prime = m.clone() 149 m_prime = m_prime.permute((1,0,2)) 150 151 vecs = m - m_prime 152 vecs = vecs.unsqueeze(0) 153 154 155 # norms = torch.tensor(scatterer.cell_normals).to(device) 156 norms = get_normals_as_points(scatterer,permute_to_points=False) 157 norms = norms.expand(1,M,-1,-1) 158 159 norm_norms = torch.norm(norms,2,dim=3) 160 vec_norms = torch.norm(vecs,2,dim=3) 161 vec_norms_cube = vec_norms**3 162 vec_norms_five = vec_norms**5 163 164 distance = torch.sqrt(torch.sum(vecs**2,dim=3)) 165 vecs_square = vecs **2 166 distance_exp = torch.unsqueeze(distance,3) 167 distance_exp = distance_exp.expand(-1,-1,-1,3) 168 169 distance_exp_cube = distance_exp**3 170 171 distaa = torch.zeros_like(distance_exp) 172 distaa[:,:,:,0] = (vecs_square[:,:,:,1] + vecs_square[:,:,:,2]) 173 distaa[:,:,:,1] = (vecs_square[:,:,:,0] + vecs_square[:,:,:,2]) 174 distaa[:,:,:,2] = (vecs_square[:,:,:,1] + vecs_square[:,:,:,0]) 175 distaa = distaa / distance_exp_cube 176 177 dista = vecs / distance_exp 178 179 Aaa = (-1 * torch.exp(1j*Constants.k * distance_exp) * (distance_exp*(1-1j*Constants.k*distance_exp))*distaa + dista*(Constants.k**2 * distance_exp**2 + 2*1j*Constants.k * distance_exp -2)) / (4*torch.pi * distance_exp_cube) 180 181 Baa = (distance_exp * distaa - 2*dista**2) / distance_exp_cube 182 183 Caa = torch.zeros_like(distance_exp).to(device) 184 185 vec_dot_norm = vecs[:,:,:,0]*norms[:,:,:,0]+vecs[:,:,:,1]*norms[:,:,:,1]+vecs[:,:,:,2]*norms[:,:,:,2] 186 187 Caa[:,:,:,0] = ((( (3 * vecs[:,:,:,0]**2) / (vec_norms_five) - (1)/(vec_norms_cube))*(vec_dot_norm)) / norm_norms) - ((2*vecs[:,:,:,0]*norms[:,:,:,0]) / (norm_norms*vec_norms_cube**3)) 188 Caa[:,:,:,1] = ((( (3 * vecs[:,:,:,1]**2) / (vec_norms_five) - (1)/(vec_norms_cube))*(vec_dot_norm)) / norm_norms) - ((2*vecs[:,:,:,1]*norms[:,:,:,1]) / (norm_norms*vec_norms_cube**3)) 189 Caa[:,:,:,2] = ((( (3 * vecs[:,:,:,2]**2) / (vec_norms_five) - (1)/(vec_norms_cube))*(vec_dot_norm)) / norm_norms) - ((2*vecs[:,:,:,2]*norms[:,:,:,2]) / (norm_norms*vec_norms_cube**3)) 190 191 Gx, Gy, Gz, A_green, B_green, C_green, Aa_green, Ba_green, Ca_green = get_G_partial(centres, scatterer, transducers, return_components=True) 192 193 Gaa = 2*Ca_green*(B_green*Aa_green + A_green*Ba_green) + C_green*(B_green*Aaa + 2*Aa_green*Ba_green + A_green*Baa)+ A_green*B_green*Caa 194 Gaa = Gaa.to(DTYPE) 195 196 areas = torch.Tensor(scatterer.celldata["Area"]).to(device) 197 areas = torch.unsqueeze(areas,0) 198 areas = torch.unsqueeze(areas,0) 199 areas = torch.unsqueeze(areas,3) 200 201 Gaa = Gaa * areas 202 # Gaa = torch.nan_to_num(Gaa) 203 eye = torch.eye(Gaa.shape[2]).to(bool) 204 Gaa[:,eye] = 0 205 206 207 A_inv_a = A_inv_a.permute(0,3,2,1) 208 Fa = Fa.permute(0,3,1,2) 209 210 A_inv = A_inv.unsqueeze(1).expand(-1,3,-1,-1) 211 Faa = Faa.permute(0,3,1,2) 212 213 Fa = Fa.to(DTYPE) 214 Faa = Faa.to(DTYPE) 215 216 Gaa = Gaa.permute(0,3,2,1) 217 Aa = Aa.permute(0,3,2,1) 218 Aa = Aa.to(DTYPE) 219 220 X1 = A_inv_a @ Fa + A_inv @ Faa 221 X2 = (A_inv @ (Aa @ A_inv @ Aa - Gaa)@A_inv) @ F 222 X3 = A_inv_a@Fa 223 224 225 Haa = X1 + X2 + X3 226 227 return Haa 228 229 230def get_cache_or_compute_H_2_gradients(scatterer:Mesh,board:Tensor,use_cache_H_grad:bool=True, path:str="Media", print_lines:bool=False) -> Tensor: 231 ''' 232 Get second derivatives of H using cache system. Expects a folder named BEMCache in `path`\n 233 :param scatterer: The mesh used (as a `vedo` `mesh` object) 234 :param board: Transducers to use 235 :param use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it 236 :param path: path to folder containing `BEMCache/ ` 237 :param print_lines: if true prints messages detaling progress 238 :return: second derivatives of H 239 ''' 240 if use_cache_H_grad: 241 242 f_name = scatterer.filename+"--"+ board_name(board) 243 f_name = hashlib.md5(f_name.encode()).hexdigest() 244 f_name = path+"/BEMCache/" + f_name +"_2grad"+ ".bin" 245 246 try: 247 if print_lines: print("Trying to load H 2 grads at", f_name ,"...") 248 Haa = pickle.load(open(f_name,"rb")) 249 Haa = Haa.to(device) 250 except FileNotFoundError: 251 if print_lines: print("Not found, computing H grad 2...") 252 Haa = grad_2_H(None, transducers=board, **{"scatterer":scatterer }) 253 f = open(f_name,"wb") 254 pickle.dump(Haa,f) 255 f.close() 256 else: 257 if print_lines: print("Computing H grad 2...") 258 Haa = grad_2_H(None, transducers=board, **{"scatterer":scatterer }) 259 260 return Haa 261 262 263def get_cache_or_compute_H_gradients(scatterer,board,use_cache_H_grad=True, path="Media", print_lines=False) -> tuple[Tensor, Tensor, Tensor]: 264 ''' 265 Get derivatives of H using cache system. Expects a folder named BEMCache in `path`\\ 266 :param scatterer: The mesh used (as a `vedo` `mesh` object)\\ 267 :param board: Transducers to use \\ 268 :param use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it\\ 269 :param path: path to folder containing BEMCache/ \\ 270 :param print_lines: if true prints messages detaling progress\\ 271 Returns derivatives of H 272 ''' 273 274 if use_cache_H_grad: 275 276 f_name = scatterer.filename +"--"+ board_name(board) 277 f_name = hashlib.md5(f_name.encode()).hexdigest() 278 f_name = path+"/BEMCache/" + f_name +"_grad"+ ".bin" 279 280 try: 281 if print_lines: print("Trying to load H grads at", f_name ,"...") 282 Hx, Hy, Hz = pickle.load(open(f_name,"rb")) 283 Hx = Hx.to(device) 284 Hy = Hy.to(device) 285 Hz = Hz.to(device) 286 except FileNotFoundError: 287 if print_lines: print("Not found, computing H Grads...") 288 Hx, Hy, Hz = grad_H(None, transducers=board, **{"scatterer":scatterer }, path=path) 289 f = open(f_name,"wb") 290 pickle.dump((Hx, Hy, Hz),f) 291 f.close() 292 else: 293 if print_lines: print("Computing H Grad...") 294 Hx, Hy, Hz = grad_H(None, transducers=board, **{"scatterer":scatterer }, path=path) 295 296 return Hx, Hy, Hz
def
grad_H( points: torch.Tensor, scatterer: vedo.mesh.Mesh, transducers: torch.Tensor, return_components: bool = False, path: str = '', H: torch.Tensor = None, use_cache_H: bool = True) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor] | tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]:
18def grad_H(points: Tensor, scatterer: Mesh, transducers: Tensor, return_components:bool = False, 19 path:str='', H:Tensor=None, use_cache_H:bool=True) ->tuple[Tensor,Tensor, Tensor] | tuple[Tensor,Tensor, Tensor, Tensor,Tensor, Tensor, Tensor]: 20 ''' 21 Computes the gradient of H wrt scatterer centres\n 22 Ignores `points` - for compatability with other gradient functions, takes centres of the scatterers 23 :param scatterer: The mesh used (as a `vedo` `mesh` object) 24 :param transducers: Transducers to use 25 :param return_components: if true will return the subparts used to compute the derivative 26 :return grad_H: The gradient of the H matrix wrt the position of the mesh 27 ''' 28 29 if H is None: 30 H = get_cache_or_compute_H(scatterer, transducers, use_cache_H, path) 31 32 33 # centres = torch.tensor(scatterer.cell_centers().points).to(device).T.unsqueeze_(0) 34 centres = get_centres_as_points(scatterer) 35 36 37 M = centres.shape[2] 38 39 B = compute_bs(scatterer,transducers) 40 A = compute_A(scatterer) 41 A_inv = torch.inverse(A).to(DTYPE) 42 43 44 Bx, By, Bz = forward_model_grad(centres, transducers) 45 Bx = Bx.to(DTYPE) 46 By = By.to(DTYPE) 47 Bz = Bz.to(DTYPE) 48 49 50 Ax, Ay, Az = get_G_partial(centres,scatterer,transducers) 51 # Ax *= -1 52 # Ay *= -1 53 # Az *= -1 54 55 Ax = (-1* Ax) 56 Ay = (-1* Ay) 57 Az = (-1* Az) 58 59 60 61 eye = torch.eye(M).to(bool) 62 Ax[:,eye] = 0 63 Ay[:,eye] = 0 64 Az[:,eye] = 0 65 66 # A_inv_x = (-1*A_inv @ Ax @ A_inv).to(DTYPE) 67 # A_inv_y = (-1*A_inv @ Ay @ A_inv).to(DTYPE) 68 # A_inv_z = (-1*A_inv @ Az @ A_inv).to(DTYPE) 69 70 # Hx_old = (A_inv_x@B) + (A_inv@Bx) 71 # Hy_old = (A_inv_y@B) + (A_inv@By) 72 # Hz_old = (A_inv_z@B) + (A_inv@Bz) 73 74 75 Hx = A_inv @ (Bx - Ax @ H) 76 Hy = A_inv @ (By - Ay @ H) 77 Hz = A_inv @ (Bz - Az @ H) 78 79 80 Hx = Hx.to(DTYPE) 81 Hy = Hy.to(DTYPE) 82 Hz = Hz.to(DTYPE) 83 84 if return_components: 85 return Hx, Hy, Hz, A, A_inv, Ax, Ay, Az 86 else: 87 return Hx, Hy, Hz
Computes the gradient of H wrt scatterer centres
Ignores points - for compatability with other gradient functions, takes centres of the scatterers
Parameters
- scatterer: The mesh used (as a
vedomeshobject) - transducers: Transducers to use
- return_components: if true will return the subparts used to compute the derivative
Returns
The gradient of the H matrix wrt the position of the mesh
def
grad_2_H( points: torch.Tensor, scatterer: vedo.mesh.Mesh, transducers: torch.Tensor, A: torch.Tensor | None = None, A_inv: torch.Tensor | None = None, Ax: torch.Tensor | None = None, Ay: torch.Tensor | None = None, Az: torch.Tensor | None = None) -> torch.Tensor:
90def grad_2_H(points: Tensor, scatterer: Mesh, transducers: Tensor, A:Tensor|None = None, 91 A_inv:Tensor|None = None, Ax:Tensor|None = None, Ay:Tensor|None = None, Az:Tensor|None = None) -> Tensor: 92 ''' 93 Computes the second derivative of H wrt scatterer centres\n 94 Ignores `points` - for compatability with other gradient functions, takes centres of the scatterers 95 :param scatterer: The mesh used (as a `vedo` `mesh` object) 96 :param transducers: Transducers to use 97 :param A: The result of a call to `compute_A` 98 :param A_inv: The inverse of `A` 99 :param Ax: The gradient of A wrt the x position of scatterer centres 100 :param Ay: The gradient of A wrt the y position of scatterer centres 101 :param Az: The gradient of A wrt the z position of scatterer centres 102 :return Haa: second order unmixed gradient of H wrt scatterer positions 103 ''' 104 105 centres = get_centres_as_points(scatterer) 106 M = centres.shape[2] 107 108 B = compute_bs(scatterer,transducers) 109 110 Fx, Fy, Fz = forward_model_grad(centres, transducers) 111 Fx = Fx.to(DTYPE) 112 Fy = Fy.to(DTYPE) 113 Fz = Fz.to(DTYPE) 114 Fa = torch.stack([Fx,Fy,Fz],dim=3) 115 116 Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(centres, transducers) 117 Faa = torch.stack([Fxx,Fyy,Fzz],dim=3) 118 119 F = forward_model_batched(centres, transducers) 120 121 if A is None: 122 A = compute_A(scatterer) 123 124 if A_inv is None: 125 A_inv = torch.inverse(A) 126 127 if Ax is None or Ay is None or Az is None: 128 Ax, Ay, Az = get_G_partial(centres,scatterer,transducers) 129 eye = torch.eye(M).to(bool) 130 Ax[:,eye] = 0 131 Ay[:,eye] = 0 132 Az[:,eye] = 0 133 Ax = Ax.to(DTYPE) 134 Ay = Ay.to(DTYPE) 135 Az = Az.to(DTYPE) 136 Aa = torch.stack([Ax,Ay,Az],dim=3) 137 138 139 A_inv_x = (-1*A_inv @ Ax @ A_inv).to(DTYPE) 140 A_inv_y = (-1*A_inv @ Ay @ A_inv).to(DTYPE) 141 A_inv_z = (-1*A_inv @ Az @ A_inv).to(DTYPE) 142 143 144 A_inv_a = torch.stack([A_inv_x,A_inv_y,A_inv_z],dim=3) 145 146 m = centres.permute(0,2,1) 147 m = m.expand((M,M,3)) 148 149 m_prime = m.clone() 150 m_prime = m_prime.permute((1,0,2)) 151 152 vecs = m - m_prime 153 vecs = vecs.unsqueeze(0) 154 155 156 # norms = torch.tensor(scatterer.cell_normals).to(device) 157 norms = get_normals_as_points(scatterer,permute_to_points=False) 158 norms = norms.expand(1,M,-1,-1) 159 160 norm_norms = torch.norm(norms,2,dim=3) 161 vec_norms = torch.norm(vecs,2,dim=3) 162 vec_norms_cube = vec_norms**3 163 vec_norms_five = vec_norms**5 164 165 distance = torch.sqrt(torch.sum(vecs**2,dim=3)) 166 vecs_square = vecs **2 167 distance_exp = torch.unsqueeze(distance,3) 168 distance_exp = distance_exp.expand(-1,-1,-1,3) 169 170 distance_exp_cube = distance_exp**3 171 172 distaa = torch.zeros_like(distance_exp) 173 distaa[:,:,:,0] = (vecs_square[:,:,:,1] + vecs_square[:,:,:,2]) 174 distaa[:,:,:,1] = (vecs_square[:,:,:,0] + vecs_square[:,:,:,2]) 175 distaa[:,:,:,2] = (vecs_square[:,:,:,1] + vecs_square[:,:,:,0]) 176 distaa = distaa / distance_exp_cube 177 178 dista = vecs / distance_exp 179 180 Aaa = (-1 * torch.exp(1j*Constants.k * distance_exp) * (distance_exp*(1-1j*Constants.k*distance_exp))*distaa + dista*(Constants.k**2 * distance_exp**2 + 2*1j*Constants.k * distance_exp -2)) / (4*torch.pi * distance_exp_cube) 181 182 Baa = (distance_exp * distaa - 2*dista**2) / distance_exp_cube 183 184 Caa = torch.zeros_like(distance_exp).to(device) 185 186 vec_dot_norm = vecs[:,:,:,0]*norms[:,:,:,0]+vecs[:,:,:,1]*norms[:,:,:,1]+vecs[:,:,:,2]*norms[:,:,:,2] 187 188 Caa[:,:,:,0] = ((( (3 * vecs[:,:,:,0]**2) / (vec_norms_five) - (1)/(vec_norms_cube))*(vec_dot_norm)) / norm_norms) - ((2*vecs[:,:,:,0]*norms[:,:,:,0]) / (norm_norms*vec_norms_cube**3)) 189 Caa[:,:,:,1] = ((( (3 * vecs[:,:,:,1]**2) / (vec_norms_five) - (1)/(vec_norms_cube))*(vec_dot_norm)) / norm_norms) - ((2*vecs[:,:,:,1]*norms[:,:,:,1]) / (norm_norms*vec_norms_cube**3)) 190 Caa[:,:,:,2] = ((( (3 * vecs[:,:,:,2]**2) / (vec_norms_five) - (1)/(vec_norms_cube))*(vec_dot_norm)) / norm_norms) - ((2*vecs[:,:,:,2]*norms[:,:,:,2]) / (norm_norms*vec_norms_cube**3)) 191 192 Gx, Gy, Gz, A_green, B_green, C_green, Aa_green, Ba_green, Ca_green = get_G_partial(centres, scatterer, transducers, return_components=True) 193 194 Gaa = 2*Ca_green*(B_green*Aa_green + A_green*Ba_green) + C_green*(B_green*Aaa + 2*Aa_green*Ba_green + A_green*Baa)+ A_green*B_green*Caa 195 Gaa = Gaa.to(DTYPE) 196 197 areas = torch.Tensor(scatterer.celldata["Area"]).to(device) 198 areas = torch.unsqueeze(areas,0) 199 areas = torch.unsqueeze(areas,0) 200 areas = torch.unsqueeze(areas,3) 201 202 Gaa = Gaa * areas 203 # Gaa = torch.nan_to_num(Gaa) 204 eye = torch.eye(Gaa.shape[2]).to(bool) 205 Gaa[:,eye] = 0 206 207 208 A_inv_a = A_inv_a.permute(0,3,2,1) 209 Fa = Fa.permute(0,3,1,2) 210 211 A_inv = A_inv.unsqueeze(1).expand(-1,3,-1,-1) 212 Faa = Faa.permute(0,3,1,2) 213 214 Fa = Fa.to(DTYPE) 215 Faa = Faa.to(DTYPE) 216 217 Gaa = Gaa.permute(0,3,2,1) 218 Aa = Aa.permute(0,3,2,1) 219 Aa = Aa.to(DTYPE) 220 221 X1 = A_inv_a @ Fa + A_inv @ Faa 222 X2 = (A_inv @ (Aa @ A_inv @ Aa - Gaa)@A_inv) @ F 223 X3 = A_inv_a@Fa 224 225 226 Haa = X1 + X2 + X3 227 228 return Haa
Computes the second derivative of H wrt scatterer centres
Ignores points - for compatability with other gradient functions, takes centres of the scatterers
Parameters
- scatterer: The mesh used (as a
vedomeshobject) - transducers: Transducers to use
- A: The result of a call to
compute_A - A_inv: The inverse of
A - Ax: The gradient of A wrt the x position of scatterer centres
- Ay: The gradient of A wrt the y position of scatterer centres
- Az: The gradient of A wrt the z position of scatterer centres
Returns
second order unmixed gradient of H wrt scatterer positions
def
get_cache_or_compute_H_2_gradients( scatterer: vedo.mesh.Mesh, board: torch.Tensor, use_cache_H_grad: bool = True, path: str = 'Media', print_lines: bool = False) -> torch.Tensor:
231def get_cache_or_compute_H_2_gradients(scatterer:Mesh,board:Tensor,use_cache_H_grad:bool=True, path:str="Media", print_lines:bool=False) -> Tensor: 232 ''' 233 Get second derivatives of H using cache system. Expects a folder named BEMCache in `path`\n 234 :param scatterer: The mesh used (as a `vedo` `mesh` object) 235 :param board: Transducers to use 236 :param use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it 237 :param path: path to folder containing `BEMCache/ ` 238 :param print_lines: if true prints messages detaling progress 239 :return: second derivatives of H 240 ''' 241 if use_cache_H_grad: 242 243 f_name = scatterer.filename+"--"+ board_name(board) 244 f_name = hashlib.md5(f_name.encode()).hexdigest() 245 f_name = path+"/BEMCache/" + f_name +"_2grad"+ ".bin" 246 247 try: 248 if print_lines: print("Trying to load H 2 grads at", f_name ,"...") 249 Haa = pickle.load(open(f_name,"rb")) 250 Haa = Haa.to(device) 251 except FileNotFoundError: 252 if print_lines: print("Not found, computing H grad 2...") 253 Haa = grad_2_H(None, transducers=board, **{"scatterer":scatterer }) 254 f = open(f_name,"wb") 255 pickle.dump(Haa,f) 256 f.close() 257 else: 258 if print_lines: print("Computing H grad 2...") 259 Haa = grad_2_H(None, transducers=board, **{"scatterer":scatterer }) 260 261 return Haa
Get second derivatives of H using cache system. Expects a folder named BEMCache in path
Parameters
- scatterer: The mesh used (as a
vedomeshobject) - board: Transducers to use
- use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it
- path: path to folder containing
BEMCache/ - print_lines: if true prints messages detaling progress
Returns
second derivatives of H
def
get_cache_or_compute_H_gradients( scatterer, board, use_cache_H_grad=True, path='Media', print_lines=False) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
264def get_cache_or_compute_H_gradients(scatterer,board,use_cache_H_grad=True, path="Media", print_lines=False) -> tuple[Tensor, Tensor, Tensor]: 265 ''' 266 Get derivatives of H using cache system. Expects a folder named BEMCache in `path`\\ 267 :param scatterer: The mesh used (as a `vedo` `mesh` object)\\ 268 :param board: Transducers to use \\ 269 :param use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it\\ 270 :param path: path to folder containing BEMCache/ \\ 271 :param print_lines: if true prints messages detaling progress\\ 272 Returns derivatives of H 273 ''' 274 275 if use_cache_H_grad: 276 277 f_name = scatterer.filename +"--"+ board_name(board) 278 f_name = hashlib.md5(f_name.encode()).hexdigest() 279 f_name = path+"/BEMCache/" + f_name +"_grad"+ ".bin" 280 281 try: 282 if print_lines: print("Trying to load H grads at", f_name ,"...") 283 Hx, Hy, Hz = pickle.load(open(f_name,"rb")) 284 Hx = Hx.to(device) 285 Hy = Hy.to(device) 286 Hz = Hz.to(device) 287 except FileNotFoundError: 288 if print_lines: print("Not found, computing H Grads...") 289 Hx, Hy, Hz = grad_H(None, transducers=board, **{"scatterer":scatterer }, path=path) 290 f = open(f_name,"wb") 291 pickle.dump((Hx, Hy, Hz),f) 292 f.close() 293 else: 294 if print_lines: print("Computing H Grad...") 295 Hx, Hy, Hz = grad_H(None, transducers=board, **{"scatterer":scatterer }, path=path) 296 297 return Hx, Hy, Hz
Get derivatives of H using cache system. Expects a folder named BEMCache in path\
Parameters
- scatterer: The mesh used (as a
vedomeshobject)\ - board: Transducers to use \
- use_cache_H_grad: If true uses the cache system, otherwise computes H and does not save it\
- path: path to folder containing BEMCache/ \
- print_lines: if true prints messages detaling progress\ Returns derivatives of H