Skip to content

Commit

Permalink
Merge pull request #349 from simpeg/fix_issue_96
Browse files Browse the repository at this point in the history
rename TRME to RME
  • Loading branch information
kkappler authored Aug 16, 2024
2 parents 63b3abb + a21b21b commit c5e42a6
Show file tree
Hide file tree
Showing 7 changed files with 20 additions and 21 deletions.
6 changes: 3 additions & 3 deletions aurora/pipelines/transfer_function_helpers.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@
from aurora.time_series.xarray_helpers import handle_nan
from aurora.transfer_function.regression.base import RegressionEstimator
from aurora.transfer_function.regression.iter_control import IterControl
from aurora.transfer_function.regression.TRME import TRME
from aurora.transfer_function.regression.TRME_RR import TRME_RR
from aurora.transfer_function.regression.RME import RME
from aurora.transfer_function.regression.RME_RR import RME_RR

# from aurora.transfer_function.weights.coherence_weights import compute_multiple_coherence_weights
from aurora.transfer_function.weights.edf_weights import (
Expand All @@ -24,7 +24,7 @@
import numpy as np
import xarray as xr

ESTIMATOR_LIBRARY = {"OLS": RegressionEstimator, "RME": TRME, "RME_RR": TRME_RR}
ESTIMATOR_LIBRARY = {"OLS": RegressionEstimator, "RME": RME, "RME_RR": RME_RR}


def get_estimator_class(estimation_engine: str) -> RegressionEstimator:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@
from scipy.linalg import solve_triangular


class TRME(MEstimator):
class RME(MEstimator):
def __init__(self, **kwargs):
"""
Constructor.
Expand All @@ -102,7 +102,7 @@ def __init__(self, **kwargs):
Its the number of points that didn't get weighted /total number of points
"""
super(TRME, self).__init__(**kwargs)
super(RME, self).__init__(**kwargs)
self.qr_input = "X"

def update_y_hat(self) -> None:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,15 @@
from loguru import logger


class TRME_RR(MEstimator):
class RME_RR(MEstimator):
def __init__(self, **kwargs):
"""
Constructor.
Robust remote reference estimator. Z is the reference station data,
and is the same size as X (see regression.base.RegressionEstimator(
"""
super(TRME_RR, self).__init__(**kwargs)
super(RME_RR, self).__init__(**kwargs)
self._Z = kwargs.get("Z", None)
self.Z = self._Z.to_array().data.T
self.qr_input = "Z"
Expand All @@ -38,7 +38,7 @@ def check_for_nan(self):
cond3 = np.isnan(self.Z).any()
nans_present = cond1 or cond2 or cond3
if nans_present:
logger.error("Missing data not allowed for TRME_RR class")
logger.error("Missing data not allowed for RME_RR class")
raise Exception

def check_for_enough_data_for_rr_estimate(self) -> None:
Expand Down
6 changes: 3 additions & 3 deletions aurora/transfer_function/regression/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -342,8 +342,8 @@ def qr_decomposition(
Parameters
----------
X: numpy array
In TRME this is the Input channels X
In TRME_RR this is the RR channels Z
In RME this is the Input channels X
In RME_RR this is the RR channels Z
sanity_check: boolean
check QR decomposition is working correctly. Set to True for debugging.
Can probably be deprecated.
Expand Down Expand Up @@ -505,7 +505,7 @@ def _get_channel_names(
Returns list of channel names.
If X is a numpy array, names will be created.
These are needed by TRME.estimate() to return xarrays.
These are needed by RME.estimate() to return xarrays.
Parameters
----------
Expand Down
4 changes: 2 additions & 2 deletions aurora/transfer_function/regression/helper_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@ def rme_beta(r0: float) -> float:
Returns a normalization factor to correct for bias in residual variance.
Details:
- This is an RME specific property. It represents a bias in the calculation of
residual_variance which we correct for in TRME and TRME_RR. The implemented
- This is a Robust M-estimator specific property. It represents a bias in the calculation of
residual_variance which we correct for in RME and RME_RR. The implemented
formula is an approximation. This is approximately equal to 1/beta where beta is
defined by Equation A3 in Egbert & Booker 1986.
Expand Down
13 changes: 6 additions & 7 deletions aurora/transfer_function/regression/m_estimator.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""
This module contains the MEstimator class - an extension of RegressionEstimator.
MEstimator has the class methods that are common to both TRME and TRME_RR.
MEstimator has the class methods that are common to both RME and RME_RR.
See Notes in RME, RME_RR for more details.
"""
Expand Down Expand Up @@ -33,8 +33,7 @@ def __init__(self, **kwargs):
then it increases linearly to 1 (r0)
Psi' is something like 1 between -r0, and r0
Psi' is zero outside
So the expectiaon value of psi' is the number of points outside
its the number of points that didnt get weighted /total number of points
So the expectation value of psi' is the number of points outside, (that didn't get weighted) divided by the total number of points
"""
super(MEstimator, self).__init__(**kwargs)
Expand Down Expand Up @@ -117,8 +116,8 @@ def residual_variance_method1(self) -> np.ndarray:
"""
returns the residual variance of the output channels.
This is the method that was originally in TRME_RR.m. It seems more correct
than the one in TRME, but also has more computational overhead.
This is the method that was originally in RME_RR.m. It seems more correct
than the one in RME, but also has more computational overhead.
"""
res = self.Yc - self.Y_hat # initial estimate of error variance
residual_variance = np.sum(np.abs(res * np.conj(res)), axis=0) / self.n_data
Expand Down Expand Up @@ -196,7 +195,7 @@ def update_y_cleaned_via_huber_weights(self) -> None:
w = np.minimum(r0s / residuals, 1.0)
self.Yc[:, k] = w * self.Y[:, k] + (1 - w) * self.Y_hat[:, k]
self.expectation_psi_prime[k] = 1.0 * np.sum(w == 1) / self.n_data
self.update_QHYc() # note the QH is different in TRME_RR vs TRME
self.update_QHYc() # note the QH is different in RME_RR vs RME
return

def initial_estimate(self) -> None:
Expand All @@ -210,7 +209,7 @@ def initial_estimate(self) -> None:

def apply_huber_regression(self) -> None:
"""
This is the 'convergence loop' from TRME, TRME_RR
This is the 'convergence loop' from RME, RME_RR
TODO: Consider not setting iter_control.number_of_iterations
- Instead, Initialize a new iter_control object
Expand Down
2 changes: 1 addition & 1 deletion tests/synthetic/test_processing.py
Original file line number Diff line number Diff line change
Expand Up @@ -318,7 +318,7 @@ def main():
# Just like process_synthetic_1, but the window is ridiculously long so that we
# encounter the underdetermined problem. We actually pass that test but in testing
# I found that at the next band over, which has more data because there are multipe
# FCs the sigma in TRME comes out as negative. see issue #4 and issue #55.
# FCs the sigma in RME comes out as negative. see issue #4 and issue #55.
# Returns
# -------
#
Expand Down

0 comments on commit c5e42a6

Please sign in to comment.