Coverage for pygeodyn/shear/shear_algo.py: 93%

1import numpy as np 

2import h5py 

3import os 

4import sklearn.covariance as skcov 

5from pygeodyn.shear import blackBoxFormula2temp as bbf 

6from pygeodyn.shear import templib_temp as templib 

7from pygeodyn.shear.generic import Generic 

8 

9class ShearInversion(Generic): 

10 """ 

11 Compute the shear at CMB at analysis time 

12 """ 

13 

14 def __init__(self, cfg, test=""): 

15 """ 

16 Init instance of ShearInverion 

17 

18 :param cfg: variable containing configuration dictionary 

19 :type cfg: dict 

20 :param test: select test prior directory (default: "") 

21 :type test: str 

22 """ 

23 

24 Generic.__init__(self, cfg, test) 

25 

26 

27 def build_prior(self): 

28 

29 # Load prior 

30 u_prior = self.load_prior_flow() 

31 y1_error_of_representativeness, y2_error_of_representativeness = self.load_prior_err() 

32 

33 # Compute prior matrices 

34 Pss_inv_with_glasso = self.compute_prior_shear(u_prior) 

35 ( 

36 Pee_y1_tilda_with_glasso, 

37 Pee_y2_tilda_with_glasso 

38 ) = self.compute_prior_error_shear(y1_error_of_representativeness, y2_error_of_representativeness) 

39 

40 return Pee_y1_tilda_with_glasso, Pee_y2_tilda_with_glasso, Pss_inv_with_glasso 

41 

42 

43 def load_prior_err(self): 

44 """ 

45 Load prior error data 

46 :return: y1_error_of_representativeness, y2_error_of_representativeness 

47 :rtype: ndarray (Ntimes x Ly(Ly + 2)), ndarray (Ntimes x Ly(Ly + 2)) 

48 """ 

49 

50 # Load prior error 

51 if self.test == "": 

52 filename = self.prior_dir_shear + "/y12_err_repr_Ly_{}_Lu_{}_Lb_{}.npy".format(self.Ly, self.Lu, self.Lb) 

53 else: 

54 filename = self.test + "/test_shear/CED_y12_err_repr_Ly_18_Lu_18_Lb_13.npy" 

55 

56 if os.path.isfile(filename): 

57 y12_error_of_representativeness = np.load(filename) 

58 y1_error_of_representativeness = y12_error_of_representativeness[:,:self.Ny] 

59 y2_error_of_representativeness = y12_error_of_representativeness[:,self.Ny:] 

60 else: 

61 raise IOError(filename + " is not found you must provide the shear error of representativeness with the coefficients Ly = {}, Lu = {}, Lb = {} and for a prior = {}. To compute it, you can use the routine parallel code compute_error.py in pygeodyn/pygeodyn/shear/".format(self.Ly, self.Lu, self.Lb, self.prior_type_shear)) 

62 

63 return y1_error_of_representativeness, y2_error_of_representativeness 

64 

65 

66 def load_prior_flow(self): 

67 """ 

68 Load prior flow 

69 :return: u_prior 

70 :rtype: ndarray (Ntimes x 2Lu(Lu + 2)) 

71 """ 

72 

73 # Load prior flow 

74 if self.test == "": 

75 if self.prior_type_shear == "coupled_earth": 

76 #sign convention coupled earth U 

77 TNM = -np.loadtxt(self.prior_dir_shear + "/RealU.dat").T 

78 elif self.prior_type_shear == "geodynamo" or self.prior_type_shear == "midpath" or self.prior_type_shear == "71path": 

79 f = h5py.File(self.prior_dir_shear + "/Real.hdf5", "r") 

80 TNM = np.array(f['U']).T 

81 else: 

82 TNM = np.loadtxt(self.test + "/test_shear/tnmsnm").T 

83 

84 tnm = TNM[: TNM.shape[0]//2, :] 

85 snm = TNM[TNM.shape[0]//2 :, :] 

86 flowcoefs = np.concatenate((-tnm[0:self.Nu, :], snm[0:self.Nu, :]), axis=0) 

87 u_prior = flowcoefs.T 

88 

89 return u_prior 

90 

91 

92 def precompute_operators(self): 

93 """ 

94 Precompute operators 

95 """ 

96 

97 ######################################### 

98 # Calculating forward matrices to calculate  

99 # theta phi maps  

100 ######################################### 

101 

102 self.delta_plus_forward_matrix_complex = bbf.calculate_u_forward_complex( 

103 self.Lu, self.gauss_thetas, self.phis, bbf.generalised_Plm_plus, debug=True 

104 ) 

105 

106 self.delta_minus_forward_matrix_complex = bbf.calculate_u_forward_complex( 

107 self.Lu, self.gauss_thetas, self.phis, bbf.generalised_Plm_minus, debug=True 

108 ) 

109 

110 self.u_plus_forward_matrix_complex = np.array(self.delta_plus_forward_matrix_complex) 

111 self.u_minus_forward_matrix_complex = np.array(self.delta_minus_forward_matrix_complex) 

112 

113 

114 ######################################### 

115 # Defining parameters to rebuild the 

116 # forward matrices from the form 

117 # [c10 cos, c11 cos, c11 sin, ... ] 

118 # to 

119 # [c1-1 real, c1-1 imag, c10 real, c10 imag, ... ] 

120 ######################################### 

121 

122 ( 

123 self.u_in_real_imag_form_dict, 

124 self.Nu_real_imag_form, 

125 u_cos_sin_dict, 

126 Nu_cossin_form, 

127 ) = bbf.build_real_imag_and_cossin_dict_pos(self.Lu) 

128 

129 ######################################### 

130 # Transformation matrix is built this way 

131 # [c10 cos, c11 cos, c11 sin, ... ] 

132 # to 

133 # [c1-1 real, c1-1 imag, c10 real, c10 imag, ... ] 

134 ######################################### 

135 

136 #u_Matrix_from_cossin_to_R_I_form = ( 

137 # bbf.building_transformation_matrix_from_cossin_to_R_I_form( 

138 # self.Nu_real_imag_form, Nu_cossin_form, self.u_in_real_imag_form_dict, u_cos_sin_dict 

139 # ) 

140 #) 

141 # 

142 #u_Matrix_from_R_I_to_cossin = ( 

143 # bbf.building_transformation_matrix_from_R_I_to_cossin_form( 

144 # self.Nu_real_imag_form, Nu_cossin_form, self.u_in_real_imag_form_dict, u_cos_sin_dict 

145 # ) 

146 #) 

147 

148 ######################################### 

149 # Matrix, that describes multiplication by 1j 

150 ######################################### 

151 

152 #Matrix_multiply_by_i = bbf.building_transformation_matrix_multiply_by_i( 

153 # self.Nu_real_imag_form, self.u_in_real_imag_form_dict 

154 #) 

155 ######################################### 

156 # Defining parameters to rebuild the 

157 # forward matrices from the form 

158 # [g10 cos, g11 cos, g11 sin, ... ] 

159 # to 

160 # [g1-1 real, g1-1 imag, g10 real, g10 imag, ... ] 

161 # for Br 

162 ######################################### 

163 

164 ( 

165 b_in_real_imag_form_dict, 

166 Nb_real_imag_form, 

167 b_cos_sin_dict, 

168 Nb_cossin_form, 

169 ) = bbf.build_real_imag_and_cossin_dict_pos(self.Lb) 

170 

171 # Same for the flow coefficients 

172 

173 #( 

174 # u_in_complex_form_dict, 

175 # Nu_complex_form, 

176 # u_cos_sin_dict, 

177 # Nu_cossin_form, 

178 #) = bbf.build_complex_and_cossin_dict_pos(self.Lu) 

179 

180 # Same for coefficients of change of magnetic field in time 

181 

182 ( 

183 _, 

184 _, 

185 _, 

186 self.Nsv_cossin_form, 

187 ) = bbf.build_real_imag_and_cossin_dict_pos(self.Lsv) 

188 

189 ######################################### 

190 # Transformation matrix is built this way 

191 # [g10 cos, g11 cos, g11 sin, ... ] 

192 # to 

193 # [g1-1 real, g1-1 imag, g10 real, g10 imag, ... ] 

194 # for Br 

195 ######################################### 

196 

197 self.b_Matrix_from_cossin_to_R_I_form = ( 

198 bbf.building_transformation_matrix_from_cossin_to_R_I_form_B( 

199 Nb_real_imag_form, Nb_cossin_form, b_in_real_imag_form_dict, b_cos_sin_dict 

200 ) 

201 ) 

202 

203 

204 ######################################### 

205 # Pre calculating intermediate maps for  

206 # the flow 

207 ######################################### 

208 

209 

210 self.u_plus_data_matrices, self.u_minus_data_matrices = bbf.calculate_u_data_matrices( 

211 self.Nu_real_imag_form, 

212 self.u_in_real_imag_form_dict, 

213 self.gauss_thetas, 

214 self.phis, 

215 debug=True, 

216 u_plus_forward_matrix=self.u_plus_forward_matrix_complex, 

217 u_minus_forward_matrix=self.u_minus_forward_matrix_complex, 

218 forward_type="complex", 

219 ) 

220 

221 

222 ######################################### 

223 # The same but for "s" term in calculation  

224 # of the shear 

225 ######################################### 

226 

227 ( 

228 self.u_plus_forward_matrix_for_s, 

229 self.u_minus_forward_matrix_for_s, 

230 ) = bbf.calculate_u_plus_minus_forward_temp_for_s(self.Lu, self.gauss_thetas, self.phis, debug=True) 

231 

232 

233 ######################################### 

234 # The same but for "v0" term in calculation  

235 # of the shear 

236 ######################################### 

237 

238 

239 self.u_zero_data_matrices = bbf.calculate_u_zero_data_matrices( 

240 self.Nu_real_imag_form, self.u_in_real_imag_form_dict, self.gauss_thetas, self.phis, debug=True 

241 ) 

242 

243 ######################################### 

244 # Matricies that allow to recalculate  

245 # u+lm and u-lm coefficients from tlm and clm 

246 ######################################### 

247 

248 self.TCtUmUp = templib.calculate_TlmClm_to_UminusUplus_matrix(self.Lu) 

249 

250 self.UmUptTC = np.linalg.inv(self.TCtUmUp.T.dot(self.TCtUmUp)).dot(self.TCtUmUp.T) 

251 

252 ######################################### 

253 # Matricies to calculate u+lm u-lm  

254 # from u-lm and backwards 

255 ######################################### 

256 

257 self.u_minus_plus_from_u_minus = templib.calculate_u_minus_plus_from_u_minus(self.Lu) 

258 

259 self.u_minus_from_u_minus_plus = templib.calculate_u_minus_from_u_minus_plus(self.Lu) 

260 

261 # Calculation of some supporting information  

262 # for inverse problem for the shear 

263 # (Generally should be removed to some library) 

264 

265 ( 

266 Y_in_real_imag_form_dict, 

267 self.Ny_real_imag_form, 

268 _, 

269 _, 

270 ) = bbf.build_real_imag_and_cossin_dict_pos(self.Ly) 

271 

272 REMOVE_LINES_Y_Real = [] 

273 self.STAYED_LINES_Y_Real = [] 

274 n = 0 

275 for l_ in range(1, self.Ly + 1): 

276 for m_ in range(-l_, l_ + 1): 

277 if m_ == 0: 

278 self.STAYED_LINES_Y_Real.append(n) 

279 REMOVE_LINES_Y_Real.append(n + 1) 

280 else: 

281 self.STAYED_LINES_Y_Real.append(n) 

282 self.STAYED_LINES_Y_Real.append(n + 1) 

283 n += 2 

284 

285 Y_positive_to_full_real = np.zeros( 

286 ( 

287 self.Ny_real_imag_form, 

288 self.Ly * (self.Ly + 2) 

289 ) 

290 ) 

291 

292 ny_positive_real = 0 

293 for ly in range(1, self.Ly + 1): 

294 for my in range(0, ly + 1): 

295 temp = (-1) ** my 

296 

297 y_pos_reim_pos = list(Y_in_real_imag_form_dict.keys())[ 

298 list(Y_in_real_imag_form_dict.values()).index((ly, my, "r")) 

299 ] 

300 y_pos_reim_neg = list(Y_in_real_imag_form_dict.keys())[ 

301 list(Y_in_real_imag_form_dict.values()).index((ly, -my, "r")) 

302 ] 

303 

304 Y_positive_to_full_real[y_pos_reim_pos, ny_positive_real] = 1 

305 Y_positive_to_full_real[y_pos_reim_neg, ny_positive_real] = temp 

306 ny_positive_real += 1 

307 

308 if my != 0: 

309 Y_positive_to_full_real[y_pos_reim_pos+1, ny_positive_real] = 1 

310 Y_positive_to_full_real[y_pos_reim_neg+1, ny_positive_real] = -temp 

311 ny_positive_real += 1 

312 

313 Y_positive_to_full_real = Y_positive_to_full_real[self.STAYED_LINES_Y_Real] 

314 

315 self.Y_full_to_positive_real = np.linalg.inv(Y_positive_to_full_real.T.dot(Y_positive_to_full_real)).dot(Y_positive_to_full_real.T) 

316 

317 

318 REMOVE_LINES_Y_Imag = [] 

319 self.STAYED_LINES_Y_Imag = [] 

320 n = 0 

321 for l_ in range(1, self.Ly + 1): 

322 for m_ in range(-l_, l_ + 1): 

323 if m_ == 0: 

324 REMOVE_LINES_Y_Imag.append(n) 

325 self.STAYED_LINES_Y_Imag.append(n + 1) 

326 else: 

327 self.STAYED_LINES_Y_Imag.append(n) 

328 self.STAYED_LINES_Y_Imag.append(n + 1) 

329 n += 2 

330 

331 Y_positive_to_full_imag = np.zeros( 

332 ( 

333 self.Ny_real_imag_form, 

334 self.Ly * (self.Ly + 2) 

335 ) 

336 ) 

337 

338 ny_positive_imag = 0 

339 for ly in range(1, self.Ly + 1): 

340 for my in range(0, ly + 1): 

341 temp = (-1) ** my 

342 

343 if my != 0: 

344 y_pos_reim_pos = list(Y_in_real_imag_form_dict.keys())[ 

345 list(Y_in_real_imag_form_dict.values()).index((ly, my, "r")) 

346 ] 

347 y_pos_reim_neg = list(Y_in_real_imag_form_dict.keys())[ 

348 list(Y_in_real_imag_form_dict.values()).index((ly, -my, "r")) 

349 ] 

350 

351 Y_positive_to_full_imag[y_pos_reim_pos, ny_positive_imag] = 1 

352 Y_positive_to_full_imag[y_pos_reim_neg, ny_positive_imag] = -temp 

353 ny_positive_imag += 1 

354 

355 Y_positive_to_full_imag[y_pos_reim_pos+1, ny_positive_imag] = 1 

356 Y_positive_to_full_imag[y_pos_reim_neg+1, ny_positive_imag] = temp 

357 ny_positive_imag += 1 

358 

359 

360 else: 

361 y_pos_reim_pos = list(Y_in_real_imag_form_dict.keys())[ 

362 list(Y_in_real_imag_form_dict.values()).index((ly, my, "i")) 

363 ] 

364 y_pos_reim_neg = list(Y_in_real_imag_form_dict.keys())[ 

365 list(Y_in_real_imag_form_dict.values()).index((ly, -my, "i")) 

366 ] 

367 

368 Y_positive_to_full_imag[y_pos_reim_pos, ny_positive_imag] = 1 

369 Y_positive_to_full_imag[y_pos_reim_neg, ny_positive_imag] = temp 

370 ny_positive_imag += 1 

371 

372 Y_positive_to_full_imag = Y_positive_to_full_imag[self.STAYED_LINES_Y_Imag] 

373 

374 self.Y_full_to_positive_imag = np.linalg.inv(Y_positive_to_full_imag.T.dot(Y_positive_to_full_imag)).dot(Y_positive_to_full_imag.T) 

375 

376 ( 

377 delta_in_re_im_form_dict, 

378 self.N_delta_re_im_form, 

379 _, 

380 _, 

381 ) = bbf.build_real_imag_and_cossin_dict_pos(self.Lu, q=True) 

382 

383 s_plus_from_s_minus = np.zeros( 

384 ( 

385 self.N_delta_re_im_form, 

386 self.N_delta_re_im_form 

387 ) 

388 ) 

389 

390 for lu in range(0, self.Lu + 1): 

391 for mu in range(-lu, lu + 1): 

392 temp = (-1) ** mu 

393 

394 s_pos_reim_plus = list(delta_in_re_im_form_dict.keys())[ 

395 list(delta_in_re_im_form_dict.values()).index((lu, mu, "r")) 

396 ] 

397 s_pos_reim_minus = list(delta_in_re_im_form_dict.keys())[ 

398 list(delta_in_re_im_form_dict.values()).index((lu, -mu, "r")) 

399 ] 

400 

401 s_plus_from_s_minus[s_pos_reim_plus, s_pos_reim_minus] = temp 

402 s_plus_from_s_minus[s_pos_reim_plus + 1, s_pos_reim_minus + 1] = -temp 

403 

404 self.s_minus_plus_from_s_minus = np.concatenate((np.eye(self.N_delta_re_im_form), s_plus_from_s_minus), axis=0) 

405 

406 #s_minus_from_s_minus_plus = np.linalg.inv(self.s_minus_plus_from_s_minus.T.dot(self.s_minus_plus_from_s_minus)).dot(self.s_minus_plus_from_s_minus.T) 

407 

408 self.ulm_to_clm_matrix = np.zeros((self.Nu_real_imag_form, self.Nu_real_imag_form)) 

409 nu = 0 

410 for lu in range(1, self.Lu + 1): 

411 for mu in range(-lu, lu + 1): 

412 temp = np.sqrt(2) / np.sqrt(lu * (lu + 1)) 

413 self.ulm_to_clm_matrix[nu, nu] = temp 

414 nu += 1 

415 self.ulm_to_clm_matrix[nu, nu] = temp 

416 nu += 1 

417 

418 

419 def compute_prior_shear(self, u_prior): 

420 """ 

421 Compute a prior matrix for the shear  

422 

423 :param u_prior: Toroidal/poloidal flow Schmidt semi-normalized spherical harmonic coefficients 

424 :type u_prior: ndarray (Ntimes x 2Lu(Lu + 2)) 

425 :return Pss_inv_with_glasso: Covariance matrix Toroidal/poloidal flow Schmidt semi-normalized spherical harmonic coefficients 

426 :rtype Pss_inv_with_glasso: ndarray (2Lu(Lu + 2) x 2Lu(Lu + 2)) 

427 """ 

428 

429 # Calculating a priory matrix for the shear  

430 

431 u_covar_coefs = np.array(u_prior) 

432 u_m_covar_coefs = templib.flowToMMatrix(u_covar_coefs, self.Lu) 

433 Pss_from_dynamo_covar_cossin = np.matmul(u_m_covar_coefs.T, u_m_covar_coefs) / (u_m_covar_coefs.shape[0] - 1) 

434 

435 # Applying graphical lasso to shear prior information 

436 

437 Css_temp = np.linalg.inv(np.diag(np.sqrt(np.diag(Pss_from_dynamo_covar_cossin)))).dot(Pss_from_dynamo_covar_cossin).dot(np.linalg.inv(np.diag(np.sqrt(np.diag(Pss_from_dynamo_covar_cossin))))) 

438 Css_tilda = skcov.graphical_lasso(Css_temp, self.glasso_lambda_u, max_iter=100)[0] 

439 Pss_tilda = np.diag(np.sqrt(np.diag(Pss_from_dynamo_covar_cossin))).dot(Css_tilda).dot(np.diag(np.sqrt(np.diag(Pss_from_dynamo_covar_cossin)))) 

440 Pss_inv_with_glasso = np.linalg.inv(self.u_minus_from_u_minus_plus.dot(self.TCtUmUp).dot(Pss_tilda).dot(self.TCtUmUp.T).dot(self.u_minus_plus_from_u_minus)) 

441 

442 return Pss_inv_with_glasso 

443 

444 

445 def compute_prior_error_shear(self, y1_error_of_representativeness, y2_error_of_representativeness): 

446 """ 

447 Computing a priory error matrix for the shear product  

448 

449 :param y1_error_of_representativeness: Toroidal flow error Schmidt semi-normalized spherical harmonic coefficients 

450 :type y1_error_of_representativeness: ndarray (Ntimes x Lu(Lu + 2)) 

451 :param y2_error_of_representativeness: Poloidal flow error Schmidt semi-normalized spherical harmonic coefficients 

452 :type y2_error_of_representativeness: ndarray (Ntimes x Lu(Lu + 2)) 

453 :return Pee_y1_tilda_with_glasso: Covariance matrix Toroidal flow error Schmidt semi-normalized spherical harmonic coefficients 

454 :rtype Pee_y1_tilda_with_glasso: ndarray (Lu(Lu + 2) x Lu(Lu + 2)) 

455 :return Pee_y2_tilda_with_glasso: Covariance matrix Poloidal flow error Schmidt semi-normalized spherical harmonic coefficients 

456 :rtype Pee_y2_tilda_with_glasso: ndarray (Lu(Lu + 2) x Lu(Lu + 2)) 

457 """ 

458 

459 # Calculating a priory error matrix for the shear product  

460 Pee_y1_apriory_temp = y1_error_of_representativeness.T.dot(y1_error_of_representativeness) / y1_error_of_representativeness.shape[0] 

461 Pee_y2_apriory_temp = y2_error_of_representativeness.T.dot(y2_error_of_representativeness) / y2_error_of_representativeness.shape[0] 

462 

463 # Applying graphical lasso to shear prior information 

464 Cee_y1_temp = np.linalg.inv(np.diag(np.sqrt(np.diag(Pee_y1_apriory_temp)))).dot(Pee_y1_apriory_temp).dot(np.linalg.inv(np.diag(np.sqrt(np.diag(Pee_y1_apriory_temp))))) 

465 Cee_y2_temp = np.linalg.inv(np.diag(np.sqrt(np.diag(Pee_y2_apriory_temp)))).dot(Pee_y2_apriory_temp).dot(np.linalg.inv(np.diag(np.sqrt(np.diag(Pee_y2_apriory_temp))))) 

466 

467 Cee_y1_tilda = skcov.graphical_lasso(Cee_y1_temp, self.glasso_lambda_err, max_iter=100)[0] 

468 Cee_y2_tilda = skcov.graphical_lasso(Cee_y2_temp, self.glasso_lambda_err, max_iter=100)[0] 

469 

470 Pee_y1_tilda_with_glasso = np.diag(np.sqrt(np.diag(Pee_y1_apriory_temp))).dot(Cee_y1_tilda).dot(np.diag(np.sqrt(np.diag(Pee_y1_apriory_temp)))) 

471 Pee_y2_tilda_with_glasso = np.diag(np.sqrt(np.diag(Pee_y2_apriory_temp))).dot(Cee_y2_tilda).dot(np.diag(np.sqrt(np.diag(Pee_y2_apriory_temp)))) 

472 

473 return Pee_y1_tilda_with_glasso, Pee_y2_tilda_with_glasso 

474 

475 

476 def compute_shear(self, U, DU, MF, SV, Pee_y1_tilda_with_glasso, Pee_y2_tilda_with_glasso, Pss_inv_with_glasso): 

477 """ 

478 Compute the shear  

479 

480 :param U: Toroidal/poloidal flow in Schmidt semi-normalized spherical harmonic coefficients 

481 :type U: ndarray (2Lu(Lu + 2)) 

482 :param DU: Toroidal/poloidal rate of change flow in Schmidt semi-normalized spherical harmonic coefficients 

483 :type DU: ndarray (2Lu(Lu + 2)) 

484 :param MF: Radial Main Field in Schmidt semi-normalized spherical harmonic coefficients 

485 :type MF: ndarray (Lb(Lb + 2)) 

486 :param SV: Radial Secular Variation in Schmidt semi-normalized spherical harmonic coefficients 

487 :type SV: ndarray (Lsv(Lsv + 2)) 

488 :param Pee_y1_tilda_with_glasso: Covariance matrix Toroidal flow error Schmidt semi-normalized spherical harmonic coefficients 

489 :type Pee_y1_tilda_with_glasso: ndarray (Lu(Lu + 2) x Lu(Lu + 2)) 

490 :param Pee_y2_tilda_with_glasso: Covariance matrix Poloidal flow error Schmidt semi-normalized spherical harmonic coefficients 

491 :type Pee_y2_tilda_with_glasso: ndarray (Lu(Lu + 2) x Lu(Lu + 2)) 

492 :param Pss_inv_with_glasso: Covariance matrix Toroidal/poloidal flow Schmidt semi-normalized spherical harmonic coefficients 

493 :type Pss_inv_with_glasso: ndarray (2Lu(Lu + 2) x 2Lu(Lu + 2)) 

494 :return SH_coef: Toroidal/poloidal shear in Schmidt semi-normalized spherical harmonic coefficients 

495 :rtype SH_coef: ndarray (2Lu(Lu + 2)) 

496 """ 

497 

498 U, DU, MF, SV = self.load_data(U, DU, MF, SV) 

499 

500 ( 

501 Br_glm_Rc_r_i, 

502 Br_data_matrix, 

503 Br_dt_data_matrix, 

504 u_true, 

505 u_dt_true 

506 ) = self.prepare_data(U, DU, MF, SV) 

507 

508 self.build_forward_matrices(Br_glm_Rc_r_i, Br_data_matrix, Br_dt_data_matrix) 

509 

510 SH_coef = self.compute_inverse( 

511 Pee_y1_tilda_with_glasso, 

512 Pee_y2_tilda_with_glasso, 

513 Pss_inv_with_glasso, 

514 u_true, 

515 u_dt_true 

516 ) 

517 

518 return SH_coef 

519 

520 

521 def load_data(self, U, DU, MF, SV): 

522 """ 

523 Load U DU MF and SV coefficient vectors at required time 

524 """ 

525 

526 measures = ["U", "DU", "MF", "SV"] 

527 for i, data in enumerate([U, DU, MF, SV]): 

528 if np.max(data.shape) == data.reshape(-1,).shape[0]: 

529 data = data.reshape(-1,) 

530 else: 

531 raise ValueError(measures[i] + " must be 1D but has dimensions " + str(data.shape)) 

532 

533 return U, DU, MF, SV 

534 

535 

536 def prepare_data(self, U, DU, MF, SV): 

537 """ 

538 Perform U, DU, MF, SV transformations 

539 """ 

540 

541 ######################################### 

542 # Defining coefficients that are used  

543 # in calculation  

544 # glm - main field coefficients 

545 # dtglm - coefficients of change of  

546 # the magnetic field 

547 # tlmclm - flow coefficients 

548 ######################################### 

549 

550 glm = np.array(MF) 

551 dtglm = np.array(SV) 

552 

553 # multiplying magnetic field coefficients by -1^m 

554 # and recalcuting it to the Earth Core radius 

555 

556 (Br_glm_Rc_2, _, _) = templib.magnFieldFromRaToRaRcM(glm, self.Lb) 

557 (gdot_Rc_coeffs_2, _, _) = templib.magnFieldFromRaToRaRcM(dtglm, self.Lb) 

558 

559 # recalcilating to another representation and building 

560 # physical values on theta phi 

561 

562 Br_glm_Rc_r_i = self.b_Matrix_from_cossin_to_R_I_form.dot(Br_glm_Rc_2) 

563 Br_data_matrix = bbf.calculate_Br_from_ReIm(self.Lb, self.gauss_thetas, self.phis, Br_glm_Rc_r_i) 

564 

565 Br_glm_Rc_dt_r_i = self.b_Matrix_from_cossin_to_R_I_form.dot(gdot_Rc_coeffs_2) 

566 Br_dt_data_matrix = bbf.calculate_Br_from_ReIm(self.Lb, self.gauss_thetas, self.phis, Br_glm_Rc_dt_r_i) 

567 

568 ######################################### 

569 # multiplying all flow coefficients by -1^m 

570 

571 u_m = templib.flowToM(U, self.Lu) 

572 

573 # Doing the same for the time change of  

574 # the flow coefficients 

575 

576 u_m_dt = templib.flowToM(DU, self.Lu) 

577 

578 # caclculating u+lm and u-lm from tlm clm 

579 

580 u_true = self.TCtUmUp.dot(u_m) 

581 u_dt_true = self.TCtUmUp.dot(u_m_dt) 

582 

583 return Br_glm_Rc_r_i, Br_data_matrix, Br_dt_data_matrix, u_true, u_dt_true 

584 

585 

586 def build_forward_matrices(self, Br_glm_Rc_r_i, Br_data_matrix, Br_dt_data_matrix): 

587 """ 

588 Build forward matrices 

589 """ 

590 

591 ######################################### 

592 self.t_plus_matrix_from_ulm, self.t_minus_matrix_from_ulm = bbf.calculate_t_from_u( 

593 self.Ly, 

594 self.Lb, 

595 self.Nu_real_imag_form, 

596 self.gauss_thetas, 

597 self.gauss_weights, 

598 self.phis, 

599 np.array(Br_glm_Rc_r_i), 

600 self.u_plus_data_matrices, 

601 self.u_minus_data_matrices, 

602 tmax=self.tmax, 

603 debug=True, 

604 ) 

605 

606 ######################################### 

607 self.v_zero_matrix_from_ulm = bbf.calculate_v0_from_u_minus_plus( 

608 self.Ly, 

609 self.Lb, 

610 self.Ny_real_imag_form, 

611 self.gauss_thetas, 

612 self.gauss_weights, 

613 self.phis, 

614 np.array(Br_glm_Rc_r_i), 

615 self.u_plus_data_matrices, 

616 self.u_minus_data_matrices, 

617 tmax=self.tmax, 

618 debug=True, 

619 ) 

620 

621 ######################################### 

622 self.p_plus_matrix_from_clm, self.p_minus_matrix_from_clm = bbf.calculate_p_from_clm( 

623 self.Ly, 

624 self.Lb, 

625 self.Ny_real_imag_form, 

626 self.Nu_real_imag_form, 

627 np.array(Br_glm_Rc_r_i), 

628 self.gauss_thetas, 

629 self.gauss_weights, 

630 self.phis, 

631 self.u_zero_data_matrices, 

632 tmax=self.tmax, 

633 debug=True, 

634 ) 

635 

636 ######################################### 

637 self.s_plus_matrix_from_delta, self.s_minus_matrix_from_delta = bbf.calculate_s_from_delta( 

638 self.Ly, 

639 self.Lb, 

640 self.Lu, 

641 np.array(Br_glm_Rc_r_i), 

642 self.gauss_thetas, 

643 self.gauss_weights, 

644 self.phis, 

645 self.u_plus_forward_matrix_for_s, 

646 self.u_minus_forward_matrix_for_s, 

647 tmax=self.tmax, 

648 debug=True, 

649 ) 

650 

651 ############################# 

652 

653 #( 

654 # _, 

655 # _, 

656 # dBr_dt_matrix_from_u_minus_plus, 

657 #) = templib.openOrCalculateDBrDtForward( 

658 # self.Lb, 

659 # self.Lsv, 

660 # self.Lq, 

661 # self.tmax, 

662 # self.Nu_real_imag_form, 

663 # self.Br_data_matrix, 

664 # self.u_in_real_imag_form_dict, 

665 # self.gauss_thetas, 

666 # self.gauss_weights, 

667 # self.phis, 

668 # self.u_plus_forward_matrix_complex, 

669 # self.u_minus_forward_matrix_complex, 

670 # print_text="Opening or calculating dBr/dt forward matrix Br...", 

671 #) 

672 

673 ( 

674 self.v_plus_matrices, 

675 self.v_minus_matrices, 

676 _, 

677 ) = templib.openOrCalculateDBrDtForward( 

678 self.Lb, 

679 self.Ly, 

680 self.Lq, 

681 self.tmax, 

682 self.Nu_real_imag_form, 

683 Br_data_matrix, 

684 self.u_in_real_imag_form_dict, 

685 self.gauss_thetas, 

686 self.gauss_weights, 

687 self.phis, 

688 self.u_plus_forward_matrix_complex, 

689 self.u_minus_forward_matrix_complex 

690 ) 

691 

692 self.v_plus_matrices = templib.from_Ra_to_Rc(self.v_plus_matrices, self.Ly) 

693 self.v_minus_matrices = templib.from_Ra_to_Rc(self.v_minus_matrices, self.Ly) 

694 

695 ( 

696 self.v_plus_matrices_dBdt, 

697 self.v_minus_matrices_dBdt, 

698 _, 

699 ) = templib.openOrCalculateDBrDtForward( 

700 self.Lb, 

701 self.Ly, 

702 self.Lq, 

703 self.tmax, 

704 self.Nu_real_imag_form, 

705 Br_dt_data_matrix, 

706 self.u_in_real_imag_form_dict, 

707 self.gauss_thetas, 

708 self.gauss_weights, 

709 self.phis, 

710 self.u_plus_forward_matrix_complex, 

711 self.u_minus_forward_matrix_complex 

712 ) 

713 

714 self.v_plus_matrices_dBdt = templib.from_Ra_to_Rc(self.v_plus_matrices_dBdt, self.Ly) 

715 self.v_minus_matrices_dBdt = templib.from_Ra_to_Rc(self.v_minus_matrices_dBdt, self.Ly) 

716 

717 

718 ( 

719 v_plus_matrices_2, 

720 v_minus_matrices_2, 

721 _, 

722 ) = templib.openOrCalculateDBrDtForward( 

723 self.Lb, 

724 self.Ly, 

725 self.Lq, 

726 self.tmax, 

727 self.Nu_real_imag_form, 

728 Br_data_matrix, 

729 self.u_in_real_imag_form_dict, 

730 self.gauss_thetas, 

731 self.gauss_weights, 

732 self.phis, 

733 self.u_plus_forward_matrix_complex, 

734 self.u_minus_forward_matrix_complex 

735 ) 

736 

737 

738 ############################### 

739 self.lv_plus_matrix = np.array(v_plus_matrices_2) 

740 self.lv_minus_matrix = np.array(v_minus_matrices_2) 

741 n = 0 

742 for lv in range(1, self.Ly): 

743 for mv in range(-lv, lv + 1): 

744 self.lv_plus_matrix[n, :] = self.lv_plus_matrix[n, :] * -lv 

745 self.lv_minus_matrix[n, :] = self.lv_minus_matrix[n, :] * -lv 

746 n += 1 

747 self.lv_plus_matrix[n, :] = self.lv_plus_matrix[n, :] * -lv 

748 self.lv_minus_matrix[n, :] = self.lv_minus_matrix[n, :] * -lv 

749 n += 1 

750 

751 self.lv_plus_matrix = templib.from_Ra_to_Rc(self.lv_plus_matrix, self.Ly) 

752 self.lv_minus_matrix = templib.from_Ra_to_Rc(self.lv_minus_matrix, self.Ly) 

753 

754 

755 def compute_inverse(self, Pee_y1_tilda_with_glasso, Pee_y2_tilda_with_glasso, Pss_inv_with_glasso, u_true, u_dt_true): 

756 """ 

757 Solve the shear inverse problem 

758 """ 

759 

760 #Ny_real_imag_stayed = len(self.STAYED_LINES_Y_Real) 

761 

762 A1 = ( 

763 -np.concatenate((self.t_minus_matrix_from_ulm, self.t_plus_matrix_from_ulm), axis=1) 

764 + self.v_zero_matrix_from_ulm 

765 - np.concatenate( 

766 ( 

767 self.p_minus_matrix_from_clm.dot(self.ulm_to_clm_matrix), 

768 self.p_plus_matrix_from_clm.dot(self.ulm_to_clm_matrix), 

769 ), 

770 axis=1, 

771 ) 

772 ) 

773 A1 = np.array(A1[:self.Ny_real_imag_form, :])[self.STAYED_LINES_Y_Real] 

774 A1 = self.Y_full_to_positive_real.dot(A1) 

775 

776 B1 = np.concatenate((self.s_minus_matrix_from_delta, self.s_plus_matrix_from_delta), axis=1) 

777 B1 = np.array(B1[:self.Ny_real_imag_form, :])[self.STAYED_LINES_Y_Real] 

778 B1_minus = (B1.dot(self.s_minus_plus_from_s_minus))[:, 2:] 

779 B1_minus = self.Y_full_to_positive_real.dot(B1_minus) 

780 

781 A1_conducting_dBdt = ( 

782 self.tau1 * np.concatenate((self.v_minus_matrices_dBdt, self.v_plus_matrices_dBdt), axis=1) 

783 ) 

784 A1_conducting_dBdt = np.array(A1_conducting_dBdt[:self.Ny_real_imag_form, :])[self.STAYED_LINES_Y_Real] 

785 A1_conducting_dBdt = self.Y_full_to_positive_real.dot(A1_conducting_dBdt) 

786 

787 A1_conducting_Br = ( 

788 self.tau1 * np.concatenate((self.v_minus_matrices, self.v_plus_matrices), axis=1) 

789 ) 

790 A1_conducting_Br = np.array(A1_conducting_Br[:self.Ny_real_imag_form, :])[self.STAYED_LINES_Y_Real] 

791 A1_conducting_Br = self.Y_full_to_positive_real.dot(A1_conducting_Br) 

792 

793 #y1_true = A1.dot(u_true) 

794 y1_true_conducting = A1.dot(u_true) + A1_conducting_dBdt.dot(u_true) + A1_conducting_Br.dot(u_dt_true) 

795 

796 ############################### 

797 

798 A2 = ( 

799 -np.concatenate((-self.t_minus_matrix_from_ulm, self.t_plus_matrix_from_ulm), axis=1) 

800 - np.concatenate( 

801 ( 

802 -self.p_minus_matrix_from_clm.dot(self.ulm_to_clm_matrix), 

803 self.p_plus_matrix_from_clm.dot(self.ulm_to_clm_matrix), 

804 ), 

805 axis=1, 

806 ) 

807 - np.concatenate((-self.lv_minus_matrix, self.lv_plus_matrix), axis=1) 

808 ) 

809 

810 A2 = np.array(A2[:self.Ny_real_imag_form, :])[self.STAYED_LINES_Y_Imag] 

811 A2 = self.Y_full_to_positive_imag.dot(A2) 

812 

813 # A2_conducting_dBdt = ( 

814 # tau1 * np.concatenate((-v_minus_matrices_dBdt, v_plus_matrices_dBdt), axis=1) 

815 # ) 

816 # A2_conducting_dBdt = np.array(A2_conducting_dBdt[:Ny_real_imag_form, :])[STAYED_LINES_Y_Imag] 

817 # A2_conducting_dBdt = Y_full_to_positive_imag.dot(A2_conducting_dBdt) 

818 

819 # A2_conducting_Br = ( 

820 # tau1 * np.concatenate((-v_minus_matrices, v_plus_matrices), axis=1) 

821 # )