diff --git a/src/scope/input.txt b/src/scope/input.txt index c29fb9e..553e195 100644 --- a/src/scope/input.txt +++ b/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. diff --git a/src/scope/input_output.py b/src/scope/input_output.py index 44ed3e7..4ee4549 100644 --- a/src/scope/input_output.py +++ b/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" ) diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py index 1278613..6bec9a9 100644 --- a/src/scope/run_simulation.py +++ b/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 diff --git a/src/scope/scrape_igrins_etc.py b/src/scope/scrape_igrins_etc.py new file mode 100644 index 0000000..2dee690 --- /dev/null +++ b/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"\s*(\d+)\s*" + 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 diff --git a/src/scope/tests/test_utils.py b/src/scope/tests/test_utils.py index e2080cf..a566e1b 100644 --- a/src/scope/tests/test_utils.py +++ b/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 diff --git a/src/scope/utils.py b/src/scope/utils.py index 9df8df3..81d1dbd 100644 --- a/src/scope/utils.py +++ b/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