diff --git a/src/scope/input.txt b/src/scope/input.txt index 6bb11c6..75d0376 100644 --- a/src/scope/input.txt +++ b/src/scope/input.txt @@ -16,6 +16,7 @@ Planet name: WASP-166 b modelname wasp_166_150seed_ld # name of the model. used for saving the output files. seed 150 # seed for the random number generator when generating noisy data. set to None for a random seed. +log_level INFO # logging level above which statements are printed. Supported levels are ``DEBUG``, ``INFO``, ``WARNING``, ``ERROR``, and ``CRITICAL``. # Filepaths planet_spectrum_path /home/asavel/scratch/retrieval_stuff/fm_wasp_77ab_noisy/fm_wasp_77ab_noisy/spectrum_tran_WASP-166b_0.3_c_to_o.pic diff --git a/src/scope/input_output.py b/src/scope/input_output.py index f357f51..649d514 100644 --- a/src/scope/input_output.py +++ b/src/scope/input_output.py @@ -2,14 +2,18 @@ Module to handle the input files. """ +import argparse import io import os -import argparse +import warnings from datetime import datetime import pandas as pd from scope.calc_quantities import * +from scope.logger import * + +logger = get_logger() class ScopeConfigError(Exception): @@ -42,14 +46,14 @@ def query_database( # check whether planet name is in database if planet_name not in df["pl_name"].values: - print(f"Planet name {planet_name} not found in database.") + logger.warning(f"Planet name {planet_name} not found in database.") try: # Use the mapped parameter name if it exists, otherwise use the original db_parameter = parameter_mapping.get(parameter, parameter) value = df.loc[df["pl_name"] == planet_name, db_parameter].values[0] return float(value) except Exception as e: - print(f"Error querying database for {planet_name}, {parameter}: {e}") + logger.warning(f"Error querying database for {planet_name}, {parameter}: {e}") return np.nan @@ -88,7 +92,7 @@ def unpack_lines(content): def coerce_nulls(data, key, value): if value == "NULL": data[key] = np.nan - + logger.warning(f" {key} is null on input") return data @@ -294,51 +298,120 @@ def write_input_file(data, output_file_path="input.txt"): f.write(f"{param:<23} {value}\n") f.write("\n") - print(f"Input file written to {output_file_path}") - + logger.info(f"Input file written to {output_file_path}") def parse_arguments(): - parser = argparse.ArgumentParser(description='Simulate observation') - - # Input file argument (required) - parser.add_argument('--input_file', type=str, default="input.txt", - help='Input file with parameters (default: input.txt)') - - # All other parameters without default values - parser.add_argument('--planet_spectrum_path', type=str, help='Path to planet spectrum') - parser.add_argument('--star_spectrum_path', type=str, help='Path to star spectrum') - parser.add_argument('--data_cube_path', type=str, help='Path to data cube') - parser.add_argument('--phase_start', type=float, help='Start phase') - parser.add_argument('--phase_end', type=float, help='End phase') - parser.add_argument('--n_exposures', type=int, help='Number of exposures') - parser.add_argument('--observation', type=str, help='Observation type') - parser.add_argument('--blaze', type=lambda x: x.lower() == 'true', help='Blaze flag') - parser.add_argument('--n_princ_comp', type=int, help='Number of principal components') - parser.add_argument('--star', type=lambda x: x.lower() == 'true', help='Star flag') - parser.add_argument('--SNR', type=float, help='Signal to noise ratio') - parser.add_argument('--telluric', type=lambda x: x.lower() == 'true', help='Telluric flag') - parser.add_argument('--tell_type', type=str, help='Telluric type') - parser.add_argument('--time_dep_tell', type=lambda x: x.lower() == 'true', help='Time dependent telluric') - parser.add_argument('--wav_error', type=lambda x: x.lower() == 'true', help='Wavelength error flag') - parser.add_argument('--rv_semiamp_orbit', type=float, help='RV semi-amplitude orbit') - parser.add_argument('--order_dep_throughput', type=lambda x: x.lower() == 'true', help='Order dependent throughput') - parser.add_argument('--Rp', type=float, help='Planet radius (Jupiter radii)') - parser.add_argument('--Rstar', type=float, help='Star radius (solar radii)') - parser.add_argument('--kp', type=float, help='Planetary orbital velocity (km/s)') - parser.add_argument('--v_rot', type=float, help='Rotation velocity') - parser.add_argument('--scale', type=float, help='Scale factor') - parser.add_argument('--v_sys', type=float, help='Systemic velocity') - parser.add_argument('--modelname', type=str, help='Model name') - parser.add_argument('--divide_out_of_transit', type=lambda x: x.lower() == 'true', help='Divide out of transit') - parser.add_argument('--out_of_transit_dur', type=float, help='Out of transit duration') - parser.add_argument('--include_rm', type=lambda x: x.lower() == 'true', help='Include RM effect') - parser.add_argument('--v_rot_star', type=float, help='Star rotation velocity') - parser.add_argument('--a', type=float, help='Semi-major axis') - parser.add_argument('--lambda_misalign', type=float, help='Misalignment angle') - parser.add_argument('--inc', type=float, help='Inclination') - parser.add_argument('--seed', type=int, help='Random seed') - parser.add_argument('--LD', type=lambda x: x.lower() == 'true', help='Limb darkening') - parser.add_argument('--vary_throughput', type=lambda x: x.lower() == 'true', help='Vary throughput') - + + parser = argparse.ArgumentParser(description="Simulate observation") + + # Required parameters + parser.add_argument( + "--planet_spectrum_path", type=str, default=".", help="Path to planet spectrum" + ) + parser.add_argument( + "--star_spectrum_path", type=str, default=".", help="Path to star spectrum" + ) + parser.add_argument( + "--data_cube_path", type=str, default=".", help="Path to data cube" + ) + + # Optional parameters with their defaults matching your function + parser.add_argument( + "--phase_start", + type=float, + default=0, + help="Start phase of the simulated observations", + ) + parser.add_argument( + "--phase_end", + type=float, + default=1, + help="End phase of the simulated observations", + ) + parser.add_argument( + "--n_exposures", type=int, default=10, help="Number of exposures" + ) + parser.add_argument( + "--observation", type=str, default="emission", help="Observation type" + ) + parser.add_argument("--blaze", type=bool, default=True, help="Blaze flag") + parser.add_argument( + "--n_princ_comp", type=int, default=4, help="Number of principal components" + ) + parser.add_argument("--star", type=bool, default=True, help="Star flag") + parser.add_argument("--SNR", type=float, default=250, help="Signal to noise ratio") + parser.add_argument("--telluric", type=bool, default=True, help="Telluric flag") + parser.add_argument( + "--tell_type", type=str, default="data-driven", help="Telluric type" + ) + parser.add_argument( + "--time_dep_tell", type=bool, default=False, help="Time dependent telluric" + ) + parser.add_argument( + "--wav_error", type=bool, default=False, help="Wavelength error flag" + ) + parser.add_argument( + "--rv_semiamp_orbit", type=float, default=0.3229, help="RV semi-amplitude orbit" + ) + parser.add_argument( + "--order_dep_throughput", + type=bool, + default=True, + help="Order dependent throughput", + ) + parser.add_argument( + "--Rp", type=float, default=1.21, help="Planet radius (Jupiter radii)" + ) + parser.add_argument( + "--Rstar", type=float, default=0.955, help="Star radius (solar radii)" + ) + parser.add_argument( + "--kp", type=float, default=192.02, help="Planetary orbital velocity (km/s)" + ) + 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( + "--modelname", type=str, default="yourfirstsimulation", help="Model name" + ) + parser.add_argument( + "--divide_out_of_transit", + type=bool, + default=False, + help="Divide out of transit", + ) + parser.add_argument( + "--out_of_transit_dur", type=float, default=0.1, help="Out of transit duration" + ) + parser.add_argument( + "--include_rm", type=bool, default=False, help="Include RM effect" + ) + parser.add_argument( + "--v_rot_star", type=float, default=3.0, help="Star rotation velocity" + ) + parser.add_argument("--a", type=float, default=0.033, help="Semi-major axis") + parser.add_argument( + "--lambda_misalign", type=float, default=0.0, help="Misalignment angle" + ) + parser.add_argument("--inc", type=float, default=90.0, help="Inclination") + parser.add_argument("--seed", type=int, default=42, help="Random seed") + parser.add_argument("--LD", type=bool, default=True, help="Limb darkening") + parser.add_argument( + "--vary_throughput", type=bool, default=True, help="Vary throughput" + ) + parser.add_argument( + "--log_level", + type=str, + default="INFO", + choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"], + help="Logging level (default: INFO)", + ) + + # For file input option + parser.add_argument( + "--input_file", type=str, default="input.txt", help="Input file with parameters" + ) + + return parser.parse_args() diff --git a/src/scope/logger.py b/src/scope/logger.py new file mode 100644 index 0000000..806be3f --- /dev/null +++ b/src/scope/logger.py @@ -0,0 +1,61 @@ +import logging + +# Global logger instance +_logger = None + + +def setup_logging(log_file=None, log_level=logging.INFO): + """ + Configure logging for the entire application + """ + global _logger + + # Create logger + _logger = logging.getLogger("simulation") + _logger.setLevel(log_level) + + # Clear any existing handlers + if _logger.hasHandlers(): + _logger.handlers.clear() + + # Create formatter + formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) + + # Create console handler + console_handler = logging.StreamHandler() + console_handler.setLevel(log_level) + console_handler.setFormatter(formatter) + _logger.addHandler(console_handler) + + # Create file handler if log_file is specified + if log_file: + file_handler = logging.FileHandler(log_file) + file_handler.setLevel(log_level) + file_handler.setFormatter(formatter) + _logger.addHandler(file_handler) + + return _logger + + +def get_logger(): + """ + Get the configured logger instance + """ + global _logger + if _logger is None: + # Configure with defaults if not already configured + _logger = logging.getLogger("simulation") + _logger.setLevel(logging.INFO) + + # Add a default console handler if no handlers exist + if not _logger.handlers: + handler = logging.StreamHandler() + formatter = logging.Formatter( + "%(asctime)s - %(name)s - %(levelname)s - %(message)s" + ) + handler.setFormatter(formatter) + _logger.addHandler(handler) + + return _logger diff --git a/src/scope/run_simulation.py b/src/scope/run_simulation.py index 85b037c..2fcabdc 100644 --- a/src/scope/run_simulation.py +++ b/src/scope/run_simulation.py @@ -4,9 +4,19 @@ calculating the log likelihood and cross-correlation function of the data given the model parameters. """ + +import logging +import os + +import numpy as np +from tqdm import tqdm + from scope.broadening import * from scope.ccf import * from scope.input_output import * +from scope.logger import * + + from scope.noise import * from scope.tellurics import * from scope.utils import * @@ -81,13 +91,17 @@ def make_data( :just_tellurics: (array) the telluric model that's multiplied to the dataset. """ - print(seed) + np.random.seed(seed) + logger.info(f"Seed set to {seed}") rv_planet, rv_star = calc_rvs( v_sys, v_sys_measured, Kp, rv_semiamp_orbit, phases ) # measured in m/s now + logger.debug(f"RV planet: {rv_planet}, RV star: {rv_star}") + + flux_cube = np.zeros( (n_order, n_exposure, n_pixel) ) # will store planet and star signal @@ -589,7 +603,7 @@ def simulate_observation( if __name__ == "__main__": args = parse_arguments() - + # First, parse the input file to get base parameters inputs = parse_input_file(args.input_file) @@ -604,5 +618,12 @@ def simulate_observation( if value is not None: inputs[key] = value + logger = setup_logging(log_level=inputs["log_level"]) + logger.debug(f"Parsed inputs: {inputs}") + + # Call the simulation function with the merged parameters - simulate_observation(**inputs) + try: + simulate_observation(**inputs) + except Exception as e: + logger.error(f"An error occurred: {e}", exc_info=True)