From 24feecc322509377507cd917170fb822114806c7 Mon Sep 17 00:00:00 2001 From: "E. Joshua Rigler" Date: Tue, 11 May 2021 15:48:09 -0600 Subject: [PATCH 1/3] Fix calculation of SECS error variance The secs.fit() function originally had an option called obs_var, which was intended to hold independent/uncorrelated error variances of each input channel as a function of time. For better consistency with Pulkkinen et al. (EPS-2003), and with the USGS Geomag-IMP package, this was modified to be obs_std, or the error standard deviation of each input channel. Practically speaking, this changes little, especially since the default values were all unity, and the square root difference between obs_var and obs_std didn't come into play. But obs_std is arguably more mathematically "correct" for weighted least-squares. The secs.fit() function still generates error variances for each SEC in the form of secs.sec_amps_var, becuase this is what is typically expected when one propogates the observation error through to the predictions. Note, however, that secs.sec_amps_var is not complete. Even though we assume uncorrelated errors in the input observations, the errors associated with the estimated SECs are not, in general, uncorrelated. But if we stored the full error covariance matrix for each time step, which would be required to fully and properly map observation errors (as specified with obs_std) to the actual SECS predictions, memory requirements would increase by a factor equal to the number of SEC basis functions, which could quickly become computationally prohibitive. --- pysecs/secs.py | 44 +++++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/pysecs/secs.py b/pysecs/secs.py index 11d2dd9..a3a800a 100644 --- a/pysecs/secs.py +++ b/pysecs/secs.py @@ -71,7 +71,7 @@ def nsec(self): nsec += len(self.sec_cf_loc) return nsec - def fit(self, obs_loc, obs_B, obs_var=None, epsilon=0.05): + def fit(self, obs_loc, obs_B, obs_std=None, epsilon=0.05): """Fits the SECS to the given observations. Given a number of observation locations and measurements, @@ -88,10 +88,10 @@ def fit(self, obs_loc, obs_B, obs_var=None, epsilon=0.05): obs_B: ndarray (ntimes x nobs x 3 [Bx, By, Bz]) An array containing the measured/observed B-fields. - obs_var : ndarray (ntimes x nobs x 3 [varX, varY, varZ]), optional - Variances in the components at each observation location. Can be used to - weight different observation locations more/less heavily. Infinite variance - effectively eliminates the observation from the fit. + obs_std : ndarray (ntimes x nobs x 3 [varX, varY, varZ]), optional + Standard error of vector components at each observation location. + This can be used to weight different observations more/less heavily. + An infinite value eliminates the observation from the fit. Default: ones(nobs x 3) equal weights epsilon : float @@ -106,9 +106,9 @@ def fit(self, obs_loc, obs_B, obs_var=None, epsilon=0.05): # Just a single snapshot given, so expand the dimensionality obs_B = obs_B[np.newaxis, ...] - # Assume unit variance of all measurements - if obs_var is None: - obs_var = np.ones(obs_B.shape) + # Assume unit standard error of all measurements + if obs_std is None: + obs_std = np.ones(obs_B.shape) ntimes = len(obs_B) @@ -121,39 +121,41 @@ def fit(self, obs_loc, obs_B, obs_var=None, epsilon=0.05): # Calculate the singular value decomposition (SVD) # NOTE: T_obs has shape (nobs, 3, nsec), we reshape it - # to (nobs*3, nsec); obs_var has shape (ntimes, nobs, 3), + # to (nobs*3, nsec); obs_std has shape (ntimes, nobs, 3), # we reshape it to (ntimes, nobs*3), then loop over ntimes # to solve using (potentially) time-dependent observation - # error variances to weight the observations + # standard errors to weight the observations for i in range(ntimes): # Only (re-)calculate SVD when necessary - if i == 0 or not np.all(obs_var[i] == obs_var[i-1]): + if i == 0 or not np.all(obs_std[i] == obs_std[i-1]): - # Weight T_obs with obs_var + # Weight T_obs with obs_std svd_in = (T_obs.reshape(-1, self.nsec) / - obs_var[i].ravel()[:, np.newaxis]) + obs_std[i].ravel()[:, np.newaxis]) # Find singular value decompostion U, S, Vh = np.linalg.svd(svd_in, full_matrices=False) - # Divide by infinity (1/S) gives zero weights + # Eliminate singular values less than epsilon by setting their + # reciprocal to zero (setting S to infinity firsts avoids + # divide-by-zero warings) + S[S < epsilon * S.max()] = np.inf W = 1./S - # Eliminate the small singular values (less than epsilon) - # by giving them zero weight - W[S < epsilon*S.max()] = 0. - - # Update VWU if obs_var changed + # Update VWU if obs_std changed VWU = Vh.T @ (np.diag(W) @ U.T) # Solve for SEC amplitudes and error variances # shape: ntimes x nsec - self.sec_amps[i, :] = (VWU @ (obs_B[i]/obs_var[i]).reshape(-1).T).T + self.sec_amps[i, :] = (VWU @ (obs_B[i] / obs_std[i]).reshape(-1).T).T # Maybe we want the variance of the predictions sometime later...? # shape: ntimes x nsec - self.sec_amps_var[i, :] = np.sum((Vh.T * W)**2, axis=1) + valid = np.isfinite(obs_std[i].reshape(-1)) + self.sec_amps_var[i, :] = np.sum( + (VWU[:,valid] * obs_std[i].reshape(-1)[valid])**2, + axis=1) return self From 6b307318b9c2f0f25357ab7523af14c213ae9c58 Mon Sep 17 00:00:00 2001 From: "E. Joshua Rigler" Date: Tue, 11 May 2021 16:18:12 -0600 Subject: [PATCH 2/3] Fix test functions I neglected to run pytest before committing recent changes to secs.py related to obs_var -> obs_std changes in the fit() method. Everything should pass now. --- tests/test_secs.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/test_secs.py b/tests/test_secs.py index 673dd76..8a3b973 100644 --- a/tests/test_secs.py +++ b/tests/test_secs.py @@ -330,17 +330,17 @@ def test_fit_multi_time(): assert_allclose(expected, secs.sec_amps) -def test_fit_obs_var(): +def test_fit_obs_std(): """Test that variance on observations changes the results.""" secs = pysecs.SECS(sec_df_loc=[[1., 0., R_EARTH + 1e6], [-1., 0., R_EARTH + 1e6]]) obs_loc = np.array([[0, 0, R_EARTH]]) obs_B = np.ones((2, 1, 3)) obs_B[1, :, :] *= 2 - obs_var = np.ones(obs_B.shape) + obs_std = np.ones(obs_B.shape) # Remove the z component from the fit of the second timestep - obs_var[1, :, 2] = np.inf - secs.fit(obs_loc, obs_B, obs_var=obs_var) + obs_std[1, :, 2] = np.inf + secs.fit(obs_loc, obs_B, obs_std=obs_std) expected = np.array([[6.40594202e+13, -7.41421248e+13], [1.382015e+14, -1.382015e+14]]) assert_allclose(expected, secs.sec_amps, rtol=1e-6) @@ -353,10 +353,10 @@ def test_fit_epsilon(): obs_loc = np.array([[0, 0, R_EARTH]]) obs_B = np.ones((2, 1, 3)) obs_B[1, :, :] *= 2 - obs_var = np.ones(obs_B.shape) + obs_std = np.ones(obs_B.shape) # Remove the z component from the fit of the second timestep - obs_var[1, :, 2] = np.inf - secs.fit(obs_loc, obs_B, obs_var=obs_var, epsilon=0.8) + obs_std[1, :, 2] = np.inf + secs.fit(obs_loc, obs_B, obs_std=obs_std, epsilon=0.8) expected = np.array([[-5.041352e+12, -5.041352e+12], [1.382015e+14, -1.382015e+14]]) assert_allclose(expected, secs.sec_amps, rtol=1e-6) From 58bdd928fbc24fe6b8bc14ba1f122e700fc6e373 Mon Sep 17 00:00:00 2001 From: "E. Joshua Rigler" Date: Tue, 11 May 2021 16:23:28 -0600 Subject: [PATCH 3/3] Fix secs.fit() usage in plot_animation.py This addresses recent changes to secs.fit() where obs_var (error variances) was changed to obs_std (error standard deviations). --- examples/plot_animation.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/plot_animation.py b/examples/plot_animation.py index 36ed916..75ca977 100644 --- a/examples/plot_animation.py +++ b/examples/plot_animation.py @@ -57,15 +57,15 @@ B_obs[:, :, 0] *= 2*np.sin(np.deg2rad(obs_lat_lon_r[:, 0])) B_obs[:, :, 1] *= 2*np.sin(np.deg2rad(obs_lat_lon_r[:, 1])) -B_var = np.ones(B_obs.shape) +B_std = np.ones(B_obs.shape) # Ignore the Z component -B_var[..., 2] = np.inf -# Can modify the variance as a function of time to +B_std[..., 2] = np.inf +# Can modify the standard error as a function of time to # see how that changes the fits too -# B_var[:, 0, 1] = 1 + ts +# B_std[:, 0, 1] = 1 + ts # Fit the data, requires observation locations and data -secs.fit(obs_loc=obs_lat_lon_r, obs_B=B_obs, obs_var=B_var) +secs.fit(obs_loc=obs_lat_lon_r, obs_B=B_obs, obs_std=B_std) # Create prediction points # Extend it a little beyond the observation points (-11, 11)