diff --git a/analytical_correct.py b/analytical_correct.py index 53cac66..e8e9dac 100644 --- a/analytical_correct.py +++ b/analytical_correct.py @@ -4,6 +4,7 @@ import parameters as params import argparse +import warnings def V(n): @@ -103,8 +104,10 @@ def adjust_phi_angle(p): rxy = ele_pos - proj_rxyz_rz x = np.cross(p, dp_loc) cos_phi = np.dot(rxy, x.T) / np.dot(np.linalg.norm(rxy, axis=1).reshape(len(rxy),1), np.linalg.norm(x, axis=1).reshape(1, len(x))) + if abs(cos_phi).max() - 1 > 1e-10: + warnings.warn("cos_phi out of [-1 - 1e-10, 1 + 1e-10]", RuntimeWarning) cos_phi = np.nan_to_num(cos_phi) - phi_temp = np.arccos(cos_phi) + phi_temp = np.arccos(np.maximum(-1, np.minimum(1, cos_phi))) phi = phi_temp range_test = np.dot(rxy, p.T) for i in range(len(r_ele)): @@ -187,6 +190,12 @@ def compute_phi(s12, s23, s34, I): I = 1. n = np.arange(1, 100) +# test for NaN bug +src_pos = [-3.850000e+00, 6.668396e+00, 0.05] +snk_pos = [-3.850000e+00, 6.668396e+00, -0.05] +s12, s23, s34 = conductivity(params.sigma_skull20) +assert not np.isnan(compute_phi(s12, s23, s34, I)).any() + for dipole in params.dipole_list: print('Now computing for dipole: ', dipole['name']) src_pos = dipole['src_pos'] diff --git a/figure3.py b/figure3.py index de1ba2e..75b1526 100644 --- a/figure3.py +++ b/figure3.py @@ -17,7 +17,7 @@ const='Sri98_no_bn1', default='Sri98', dest='sri98', - help='a path to the result directory') + help='use Analytical_Sri98_no_bn1*.npz files') args = parser.parse_args() diff --git a/figure4.py b/figure4.py index 989b639..542dea1 100644 --- a/figure4.py +++ b/figure4.py @@ -22,7 +22,7 @@ const='Sri98_no_bn1', default='Sri98', dest='sri98', - help='a path to the result directory') + help='use Analytical_Sri98_no_bn1*.npz files') args = parser.parse_args()