pca removal schemes

src/scope/input.txt

@@ -67,6 +67,7 @@ tell_type              data-driven   # type of telluric simulation. supported mo
 time_dep_tell          False         # whether the tellurics are time-dependent or not.
 
 # Analysis Parameters
+pca_removal            subtract      # how to remove the PCA fit. options are ``subtract`` and ``divide``; the former weights absolute errors, and the latter relative errors.
 n_princ_comp           4             # number of principal components to remove from the simulated data before cross-correlation.
 divide_out_of_transit  False         # *experimental*. Whether to divide the in-transit data by median out-of-transit or not. only used if observation is set to transmission.
 out_of_transit_dur     1.          # *experimental*. Duration of the out-of-transit data, in units of transit duration. only used if observation is set to transmission and divide_out_of_transit is set to True.

src/scope/input_output.py

@@ -378,6 +378,9 @@ def parse_arguments():
     parser.add_argument("--v_rot", type=float, default=4.5, help="Rotation velocity")
     parser.add_argument("--scale", type=float, default=1.0, help="Scale factor")
     parser.add_argument("--v_sys", type=float, default=0.0, help="Systemic velocity")
+    parser.add_argument(
+        "--pca_rmeove", type=str, default="subtract", help="PCA removal scheme"
+    )
     parser.add_argument(
         "--modelname", type=str, default="yourfirstsimulation", help="Model name"
     )

src/scope/run_simulation.py

@@ -59,6 +59,7 @@ def make_data(
     v_sys_measured=0.0,
     vary_throughput=True,
     instrument="IGRINS",
+    pca_removal="subtract",
 ):
     """
     Creates a simulated HRCCS dataset. Main function.
@@ -258,10 +259,9 @@ def make_data(
             flux_cube[j] -= np.mean(flux_cube[j])
             flux_cube[j] /= np.std(flux_cube[j])
             flux_cube[j], pca_noise_matrix[j] = perform_pca(
-                flux_cube[j], n_princ_comp, return_noplanet=True
+                flux_cube[j], n_princ_comp, pca_removal=pca_removal
             )
-            # todo: think about the svd
-            # todo: redo all analysis centering on 0?
+
     else:
         for j in range(n_order):
             for i in range(n_exposure):
@@ -309,6 +309,7 @@ def calc_log_likelihood(
     LD=True,
     b=0.0,  # impact parameter
     v_sys_measured=0.0,
+    pca_removal="subtract",
 ):
     """
     Calculates the log likelihood and cross-correlation function of the data given the model parameters.
@@ -396,7 +397,9 @@ def calc_log_likelihood(
         if do_pca:
             # this is the "model reprocessing" step.
             model_flux_cube *= A_noplanet[order]
-            model_flux_cube, _ = perform_pca(model_flux_cube, n_princ_comp, False)
+            model_flux_cube, _ = perform_pca(
+                model_flux_cube, n_princ_comp, pca_removal=pca_removal
+            )
 
         logl, ccf = calc_ccf(model_flux_cube, flux_cube[order], n_pixel)
         CCF += ccf
@@ -445,6 +448,7 @@ def simulate_observation(
     planet_name="yourfirstplanet",
     n_kp=200,
     n_vsys=200,
+    pca_removal="subtract",
     **kwargs,
 ):
     """
@@ -560,6 +564,7 @@ def simulate_observation(
         v_sys_measured=v_sys,
         LD=LD,
         instrument=instrument,
+        pca_removal=pca_removal,
     )
 
     run_name = f"{n_princ_comp}_NPC_{blaze}_blaze_{star}_star_{telluric}_telluric_{SNR}_SNR_{tell_type}_{time_dep_tell}_{wav_error}_{order_dep_throughput}_{seed}"
@@ -593,6 +598,7 @@ def simulate_observation(
                 observation=observation,
                 v_sys_measured=v_sys,
                 LD=LD,
+                pca_removal=pca_removal,
             )
             lls[l, k], ccfs[l, k] = res
 

src/scope/scrape_igrins_etc.py

@@ -0,0 +1,71 @@
+"""
+pings the IGRINS exposure time calculator.
+"""
+import re
+
+from selenium import webdriver
+from selenium.common.exceptions import ElementNotInteractableException
+from selenium.webdriver.common.by import By
+
+from scope.logger import *
+
+logger = get_logger()
+
+
+def scrape_igrins_etc(kmag, exposure_time):
+    """
+    pings the IGRINS exposure time calculator.
+
+    Parameters
+    ----------
+    kmag : float
+        K-band magnitude of the target.
+    exptime : float
+        Exposure time in seconds.
+
+    Returns
+    -------
+    float
+        Exposure time in seconds.
+    """
+
+    driver = webdriver.Chrome()
+    web_address = "https://igrins-jj.firebaseapp.com/etc/simple"
+
+    driver.get(web_address)
+
+    try:
+        input_fields = driver.find_elements(By.TAG_NAME, "input")
+        # first is kmag, second is SNR, third is kmag, fourth is integration time
+        for i, input_field in enumerate(input_fields):
+            input_field.clear()
+            if i in [0, 2]:
+                input_field.send_keys(str(kmag))
+            if i == 3:
+                input_field.send_keys(str(exposure_time))
+        response_page = driver.page_source
+
+    # now need to grab the exposure time at the end
+
+    except ElementNotInteractableException:
+        logger.error("Could not interact with IGRINS SNR website.")
+
+    snr_value = extract_snr_from_html(response_page)
+
+    return snr_value
+
+
+def extract_snr_from_html(html_content):
+    """
+    Extracts the SNR value from the HTML content of the IGRINS SNR calculator.
+
+    """
+    # Regular expression to match the number following the comment structure for SNR
+    pattern = r"<!-- react-text: 43 -->\s*(\d+)\s*<!-- /react-text -->"
+    match = re.search(pattern, html_content)
+    if match:
+        snr_value = match.group(1)
+        logger.info(f"Extracted SNR value: {snr_value}")
+    else:
+        logger.warning("SNR value not found.")
+    return snr_value

src/scope/tests/test_utils.py

@@ -77,7 +77,7 @@ def test_perform_pca_shapes(self):
         """
 
         n_princ_comp = 5
-        scaled_cube, pca = perform_pca(self.cube, n_princ_comp, return_noplanet=True)
+        scaled_cube, pca = perform_pca(self.cube, n_princ_comp)
         assert scaled_cube.shape == self.cube.shape
         assert pca.shape == self.cube.shape
 
@@ -87,7 +87,7 @@ def test_perform_pca_variance_removed(self):
         """
 
         n_princ_comp = 5
-        scaled_cube, pca = perform_pca(self.cube, n_princ_comp, return_noplanet=True)
+        scaled_cube, pca = perform_pca(self.cube, n_princ_comp)
         assert np.std(scaled_cube) < np.std(
             pca
         )  # the variance should be contained in the main components
@@ -98,7 +98,7 @@ def test_perform_pca_variance_removed_compared_input(self):
         """
 
         n_princ_comp = 5
-        scaled_cube, pca = perform_pca(self.cube, n_princ_comp, return_noplanet=True)
+        scaled_cube, pca = perform_pca(self.cube, n_princ_comp)
         assert np.std(scaled_cube) < np.std(
             self.cube
         )  # the variance should be contained in the main components
@@ -109,10 +109,11 @@ def test_perform_pca_variance_removed_compared_input_diffcompoents(self):
         """
 
         n_princ_comp = 5
-        scaled_cube5, pca5 = perform_pca(self.cube, n_princ_comp, return_noplanet=True)
+        scaled_cube5, pca5 = perform_pca(self.cube, n_princ_comp)
         n_princ_comp = 10
         scaled_cube10, pca10 = perform_pca(
-            self.cube, n_princ_comp, return_noplanet=True
+            self.cube,
+            n_princ_comp,
         )
         assert np.std(scaled_cube10) < np.std(
             scaled_cube5

src/scope/utils.py

@@ -234,33 +234,45 @@ def calc_limb_darkening(u1, u2, a, b, Rstar, ph, LD):
 
 
 @njit
-def perform_pca(input_matrix, n_princ_comp, return_noplanet=False):
+def perform_pca(input_matrix, n_princ_comp, pca_removal="subtract"):
     """
-    Perform PCA using SVD.
+    Perform PCA using SVD and remove principal components via subtraction or division.
 
-    SVD is written as A = USV^H, where ^H is the Hermitian operator.
+    Parameters
+    ----------
+    input_matrix : ndarray
+        2D array where rows are observations and columns are features.
+    n_princ_comp : int
+        Number of principal components to retain.
+    mode : str, optional
+        Method for removing PCA fit, either "subtract" or "divide".
+        Default is "subtract".
 
-    Inputs
-    ------
-        :input_matrix:
-        :n_princ_comp: number of principal components to keep
+    Returns
+    -------
+    corrected_data : ndarray
+        The input matrix with PCA components removed.
+    A_pca_fit : ndarray
+        The PCA reconstruction used for removal.
     """
-    u, singular_values, vh = np.linalg.svd(
-        input_matrix, full_matrices=False
-    )  # decompose
-    if return_noplanet:
-        s_high_variance = singular_values.copy()
-        s_high_variance[n_princ_comp:] = 0.0  # keeping the high-variance terms here
-        s_matrix = np.diag(s_high_variance)
-        A_noplanet = np.dot(u, np.dot(s_matrix, vh))
-
-    singular_values[:n_princ_comp] = 0.0  # zero out the high-variance terms here
-    s_matrix = np.diag(singular_values)
-    arr_planet = np.dot(u, np.dot(s_matrix, vh))
-
-    if return_noplanet:
-        return arr_planet, A_noplanet
-    return arr_planet, arr_planet
+    # Perform SVD decomposition
+    u, singular_values, vh = np.linalg.svd(input_matrix, full_matrices=False)
+
+    # Construct PCA model of systematics
+    s_high_variance = singular_values.copy()
+    s_high_variance[n_princ_comp:] = 0.0  # Keep the high-variance trends
+    s_matrix = np.diag(s_high_variance)
+    A_pca_fit = np.dot(u, np.dot(s_matrix, vh))  # Systematics model
+
+    # Apply chosen removal method
+    if pca_removal == "subtract":
+        corrected_data = input_matrix - A_pca_fit
+    elif pca_removal == "divide":
+        corrected_data = input_matrix / (A_pca_fit + 1e-8)  # Avoid division by zero
+    else:
+        raise ValueError("Mode must be 'subtract' or 'divide'.")
+
+    return corrected_data, A_pca_fit
 
 
 @njit