From 276d0d5844a08f53f8a827d1518972128de05397 Mon Sep 17 00:00:00 2001 From: "huluvu38@gmail.com" Date: Fri, 9 Feb 2024 13:06:57 +0100 Subject: [PATCH 1/5] Added a tool to examine lightcurves and select problematic data points for masking. --- pyLIMA/toolbox/examine_lightcurve.py | 141 +++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 pyLIMA/toolbox/examine_lightcurve.py diff --git a/pyLIMA/toolbox/examine_lightcurve.py b/pyLIMA/toolbox/examine_lightcurve.py new file mode 100644 index 0000000..95b4f24 --- /dev/null +++ b/pyLIMA/toolbox/examine_lightcurve.py @@ -0,0 +1,141 @@ +import matplotlib +matplotlib.use('TKAgg') # Set the backend before importing pyplot +import numpy as np +import matplotlib.pyplot as plt +from sys import argv, exit + +# Examine a lightcurve to get the indexes of a list of points that can be used for masking + +class PointBrowser: + """ + A tool for interactive point selection and navigation in a matplotlib + lightcurve plot. + + Click on a point to highlight it. Navigate through points using 'n'/'right' + and 'p'/'left' keys for next and previous points. Add a new point to the list + by pressing 'a' or remove it by pressing 'r'. Display the list of indexes of the + selected points in the terminal by pressing 't'. + """ + + def __init__(self, fig, ax, hjd, mag, idx, line): + """ + Initializes the PointBrowser instance with plot axis, data, and index. + + Parameters: + - fig: The matplotlib figure object for the plot. + - ax: The matplotlib axis object for plotting. + - hjd, mag: Arrays of data points for plotting. + - idx: Array of indices corresponding to data points. + - line: The matplotlib line object for interaction. + """ + self.fig = fig + self.ax = ax + self.hjd = hjd + self.mag = mag + self.idx = idx + self.line = line + self.lastind = 0 + self.text = ax.text(0.05, 0.95, 'selected: none', + transform=ax.transAxes, va='top') + self.selected, = ax.plot([], [], 'o', ms=8, alpha=0.7, color='teal', visible=False) + self.added_points, = ax.plot([], [], 'rx', ms=6, mew=2, alpha=0.7) # For visualizing added points + self.output = [] + + def onpress(self, event): + """Handles keyboard press events for navigating and modifying points.""" + if self.lastind is None or event.key not in ('n', 'p', 'a', 'r', 't', 'left', 'right'): + return + inc = 0 + if event.key in ('n', 'right'): + inc = 1 + elif event.key in ('p', 'left'): + inc = -1 + elif event.key == 'a' and self.idx[self.lastind] not in self.output: + self.output.append(self.idx[self.lastind]) + self.update_added_points() # Update visualization of added points + return + elif event.key == 'r' and self.idx[self.lastind] in self.output: + self.output.remove(self.idx[self.lastind]) + self.update_added_points() # Update visualization of added points + return + elif event.key == 't': + self.output.sort() + print(self.output) + return + + self.lastind = (self.lastind + inc) % len(self.hjd) # Cycle through points + self.update() + + def update_added_points(self): + """Updates the visualization of added points.""" + if not self.output: + self.added_points.set_data([], []) # No points to display + else: + # Get the hjd and mag values for added points + added_hjd = [self.hjd[i] for i in self.output] + added_mag = [self.mag[i] for i in self.output] + self.added_points.set_data(added_hjd, added_mag) + self.fig.canvas.draw() + + def onpick(self, event): + """Handles mouse pick events to select points.""" + if event.artist != self.line or not len(event.ind): + return True + + distances = np.hypot(event.mouseevent.xdata - self.hjd[event.ind], + event.mouseevent.ydata - self.mag[event.ind]) + indmin = distances.argmin() + self.lastind = event.ind[indmin] + self.update() + + def update(self): + """Updates the plot with the selected point highlighted.""" + if self.lastind is None: + return + dataind = self.lastind + self.selected.set_visible(True) + self.selected.set_data(self.hjd[dataind], self.mag[dataind]) + self.text.set_text(f'selected: {self.idx[dataind]}, HJD: {self.hjd[dataind]:.2f}') + self.update_added_points() # Ensure added points are updated alongside + self.fig.canvas.draw() + +if __name__ == '__main__': + if len(argv) != 2: + print('Usage: python examine_lightcurve.py .dat') + exit() + + datafile = argv[1] + try: + try: + hjd, mag, merr = np.loadtxt(datafile, delimiter=',', skiprows=1, usecols=(0, 1, 2), + dtype=float, unpack=True) + except: + hjd, mag, merr = np.loadtxt(datafile, delimiter=',', skiprows=1, usecols=(0, 2, 3), + dtype=float, unpack=True) + idx = np.arange(len(hjd)) + except Exception as e: + print(f'Error reading {datafile}: {e}') + exit() + + fig, ax = plt.subplots() + ax.set_title(datafile) + + # Set the labels for the x-axis and y-axis + ax.set_xlabel("Julian Date") + ax.set_ylabel("Magnitude") + line, = ax.plot(hjd, mag, 'k.', picker=5) # 5 points tolerance + plt.ylim(max(mag) + 0.1, min(mag) - 0.1) + ax.grid(True, which='both', color='gray', linestyle='-', linewidth=0.5, alpha=0.5) + + # Customize the x-axis with exact Julian dates + #tick_interval = 30 # Adjust this based on your dataset's range and density + #ticks = np.arange(min(hjd), max(hjd), tick_interval) + #ax.set_xticks(ticks) + #ax.set_xticklabels([f'{tick:.0f}' for tick in ticks], rotation=25, ha='right') # Format for clarity + + browser = PointBrowser(fig, ax, hjd, mag, idx, line) + fig.canvas.mpl_connect('pick_event', browser.onpick) + fig.canvas.mpl_connect('key_press_event', browser.onpress) + + plt.show() + From a7087c1ef82f94a1d671cb1f0ff1d3491fac9000 Mon Sep 17 00:00:00 2001 From: "huluvu38@gmail.com" Date: Tue, 13 Feb 2024 14:45:09 +0100 Subject: [PATCH 2/5] added /toolbox/bin_lightcurve.py and toolbox/examine_lightcurve.py with appropriate unit tests in tests/test_toolbox.py --- pyLIMA/tests/test_toolbox.py | 99 ++++++++++++++++++++++- pyLIMA/toolbox/bin_lightcurve.py | 114 +++++++++++++++++++++++++++ pyLIMA/toolbox/examine_lightcurve.py | 3 +- 3 files changed, 214 insertions(+), 2 deletions(-) create mode 100644 pyLIMA/toolbox/bin_lightcurve.py diff --git a/pyLIMA/tests/test_toolbox.py b/pyLIMA/tests/test_toolbox.py index 13c8a19..0ae8d16 100644 --- a/pyLIMA/tests/test_toolbox.py +++ b/pyLIMA/tests/test_toolbox.py @@ -1,6 +1,10 @@ import numpy as np +import unittest +from unittest.mock import MagicMock, Mock from pyLIMA.toolbox import brightness_transformation - +from pyLIMA.toolbox import brightness_transformation +from pyLIMA.toolbox.examine_lightcurve import PointBrowser # Adjust the import path as needed +from pyLIMA.toolbox.bin_lightcurve import weighted_mean, get_sigma, bin_data # Adjust the import path as needed def test_magnitude_to_flux(): flux = brightness_transformation.magnitude_to_flux(18.76) @@ -37,3 +41,96 @@ def test_noisy_observations(): flux_obs = brightness_transformation.noisy_observations(flux, exp_time=None) assert flux_obs != flux + +class TestPointBrowser(unittest.TestCase): + + def setUp(self): + # Mocking the matplotlib figure, axis, and line objects + self.fig = MagicMock() + self.ax = MagicMock() + self.line = MagicMock() + + # Example data + self.hjd = np.array([1, 2, 3, 4, 5]) + self.mag = np.array([10, 11, 9, 12, 8]) + self.idx = np.arange(len(self.hjd)) + + # Configure the ax mock to return a mock object when plot is called + mock_plot_return = Mock() + self.ax.plot.return_value = (mock_plot_return,) + + # Instantiate PointBrowser with mocked data + self.browser = PointBrowser(self.fig, self.ax, self.hjd, self.mag, self.idx, self.line) + + def test_initialization(self): + """Test that PointBrowser initializes with expected attributes.""" + self.assertEqual(self.browser.lastind, 0) + self.assertEqual(len(self.browser.output), 0) + self.assertIsNotNone(self.browser.fig) + self.assertIsNotNone(self.browser.ax) + self.assertIsNotNone(self.browser.line) + self.assertTrue(np.array_equal(self.browser.hjd, self.hjd)) + self.assertTrue(np.array_equal(self.browser.mag, self.mag)) + self.assertTrue(np.array_equal(self.browser.idx, self.idx)) + + def test_add_point(self): + # Simulate adding a point and test the outcome + self.browser.output.append(2) # Simulate action + self.assertIn(2, self.browser.output) # Test condition + + def test_remove_point(self): + # First, add a point to ensure there's something to remove + self.browser.output.append(2) + self.assertIn(2, self.browser.output, "Point was not added correctly.") + + # Now, simulate removing the point + self.browser.output.remove(2) + self.assertNotIn(2, self.browser.output, "Point was not removed correctly.") + +def test_examine_lightcurve(): + """Function to run all tests.""" + suite = unittest.TestLoader().loadTestsFromTestCase(TestPointBrowser) + unittest.TextTestRunner().run(suite) + +class TestBinLightcurve(unittest.TestCase): + + def setUp(self): + self.hjd = np.array([1, 2, 3, 4, 5]) + self.mags = np.array([10, 12, 11, 13, 14]) + self.merrs = np.array([0.5, 0.2, 0.3, 0.4, 0.1]) + + def test_weighted_mean_types(self): + """Test that weighted_mean function accepts correct types.""" + self.assertIsInstance(self.mags, np.ndarray, "Magnitudes must be a NumPy array") + self.assertIsInstance(self.merrs, np.ndarray, "Errors must be a NumPy array") + + # Perform a basic operation to ensure the function works with correct types + result = weighted_mean(self.mags, self.merrs) + self.assertIsInstance(result, float, "Weighted mean should be a float") + + def test_get_sigma_types(self): + """Test that get_sigma function accepts correct types.""" + self.assertIsInstance(self.merrs, np.ndarray, "Errors must be a NumPy array") + + # Perform a basic operation to ensure the function works with correct types + result = get_sigma(self.merrs) + self.assertIsInstance(result, float, "Sigma should be a float") + + def test_bin_data_types(self): + """Test that bin_data function accepts correct types.""" + self.assertIsInstance(self.hjd, np.ndarray, "HJDs must be a NumPy array") + self.assertIsInstance(self.mags, np.ndarray, "Magnitudes must be a NumPy array") + self.assertIsInstance(self.merrs, np.ndarray, "Errors must be a NumPy array") + self.assertIsInstance(1, (int, float), "tstart must be numeric") + self.assertIsInstance(5, (int, float), "tend must be numeric") + self.assertIsInstance(2, (int, float), "bin_size must be numeric") + + # Perform a basic operation to ensure the function works with correct types + binned_data = bin_data(self.hjd, self.mags, self.merrs, 1, 5, 2) + self.assertIsInstance(binned_data, np.ndarray, "Binned data should be a NumPy array") + self.assertEqual(binned_data.shape[1], 3, "Binned data should have three columns") + +def test_bin_lightcurve(): + """Function to run all tests in the lightcurve processing test suite.""" + suite = unittest.TestLoader().loadTestsFromTestCase(TestBinLightcurve) + unittest.TextTestRunner().run(suite) diff --git a/pyLIMA/toolbox/bin_lightcurve.py b/pyLIMA/toolbox/bin_lightcurve.py new file mode 100644 index 0000000..871919e --- /dev/null +++ b/pyLIMA/toolbox/bin_lightcurve.py @@ -0,0 +1,114 @@ +import numpy as np +from astropy.timeseries import TimeSeries +import os +import sys + +def weighted_mean(my_mags, my_errs): + ''' Function to return the weighted mean for the magnitude array + (inverse-variance weighting) + my_mags: magnitude measurements + my_merrs: uncertainties in the measurements + ''' + weights = 1.0/my_errs**2 + weighted_mean_mag = np.sum(my_mags * weights) / np.sum(weights) + + return weighted_mean_mag + + +def get_sigma(my_errs): + ''' Function to return a single value for the error array + merrs: uncertainties in the measurements + ''' + weights = 1.0/my_errs**2 + new_sigma = np.sqrt( np.sum(weights**2 * my_errs**2) / np.sum(weights)**2) + + return new_sigma + +def bin_data(hjds, mags, merrs, tstart, tend, bin_size): + ''' Function to bin a timeseries data set + hjds: heliocentric julian dates corresponding to the magnitude measurements + mags: magnitude measurements at those julian dates + merrs: uncertainties in the measurements at those julian dates + tstart: julian day to start from + tend: julian day to stop at + bin_size: the bin size in days + ''' + new_hjds = [] + new_mags = [] + new_merrs = [] + timestamp = tstart + while timestamp < tend: + next_timestamp = timestamp + bin_size + timeclips = np.where( (timestamp <= hjds) & (hjds < next_timestamp) ) + hjds_temp = hjds[timeclips] + mags_temp = mags[timeclips] + merrs_temp = merrs[timeclips] + + if hjds_temp.size != 0: + new_hjds.append(np.mean(hjds_temp)) + new_merrs.append(get_sigma(merrs_temp)) + new_mags.append(weighted_mean(mags_temp,merrs_temp)) + + timestamp = next_timestamp + + binned_data = np.empty([len(new_hjds),3]) + binned_data[:,0] = new_hjds + binned_data[:,1] = new_mags + binned_data[:,2] = new_merrs + + return binned_data + +if __name__ == '__main__': + run_check = True # Set to False if you don't want to get a plotting check at the end + # Check for command-line arguments + # Assumes the first line in the light curve file includes the column descriptions: time, magnitude, error + if len(sys.argv) < 2: + print("Usage: python bin_lightcurve.py [output_file]") + sys.exit(1) + + input_file = sys.argv[1] + if len(sys.argv) == 3: + output_file = sys.argv[2] + else: + # Create a default output filename by appending '_binned' before the file extension + base, ext = os.path.splitext(input_file) + output_file = f"{base}_binned{ext}" + + # Data points to be masked by index (if any), default is no masking + bad_data = [] + + # Read data as astropy TimeSeries Table and remove any masked data points + ts = TimeSeries.read(input_file, time_column='time', time_format='mjd') + ts.remove_rows(bad_data) + + # Split up data for processing and define binning ranges + hjds = ts['time'].value + mags = ts['magnitude'].value + merrs = ts['error'].value + + tstart = hjds[0] + tend = hjds[-1]+0.5 + + bin_size = 1.0 # Default value is daily binning when run from the command line + + binned_data = bin_data(hjds, mags, merrs, tstart, tend, bin_size) + + # Save the binned data + np.savetxt(output_file, binned_data, delimiter=',') + + # Plotting check (if needed) + if run_check: + import matplotlib.pyplot as plt + + new_data = bin_data(hjds, mags, merrs, tstart, tend, bin_size) + + plt.errorbar(new_data[:,0],new_data[:,1],yerr=new_data[:,2], fmt='ro', label="binned") + plt.errorbar(hjds,mags,yerr=merrs, fmt='k.',alpha=0.2, label="unbinned") + plt.gca().invert_yaxis() + # Add labels to the x-axis and y-axis + plt.xlabel("Julian Date") + plt.ylabel("Magnitude") + plt.grid(True, which='both', color='gray', linestyle='-', linewidth=0.5, alpha=0.5) + plt.legend() + plt.show() + diff --git a/pyLIMA/toolbox/examine_lightcurve.py b/pyLIMA/toolbox/examine_lightcurve.py index 95b4f24..26ecf97 100644 --- a/pyLIMA/toolbox/examine_lightcurve.py +++ b/pyLIMA/toolbox/examine_lightcurve.py @@ -101,11 +101,12 @@ def update(self): if __name__ == '__main__': if len(argv) != 2: - print('Usage: python examine_lightcurve.py .dat') + print('Usage: python examine_lightcurve.py /path/to/.csv') exit() datafile = argv[1] try: + # Skip the column descriptions in the first row try: hjd, mag, merr = np.loadtxt(datafile, delimiter=',', skiprows=1, usecols=(0, 1, 2), dtype=float, unpack=True) From 6306e9b203891a460089e87ca750904ef3673b15 Mon Sep 17 00:00:00 2001 From: ytsapras Date: Thu, 15 Feb 2024 15:46:24 +0100 Subject: [PATCH 3/5] Functions to save and load pickled pyLIMA model and fit objects. --- pyLIMA/outputs/share_results.py | 40 +++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 pyLIMA/outputs/share_results.py diff --git a/pyLIMA/outputs/share_results.py b/pyLIMA/outputs/share_results.py new file mode 100644 index 0000000..cdfa05f --- /dev/null +++ b/pyLIMA/outputs/share_results.py @@ -0,0 +1,40 @@ +import pickle + +def save_results(model_object, fit_object, filename, run_silent=True): + """ + Serializes and saves a model object and a fit object to a file. + + Parameters: + - model_object: The pyLIMA model object to be saved. + - fit_object: The pyLIMA fit object related to the model object to be saved. + - filename: String specifying the path to save the pickle file. + + Returns: + - None + """ + with open(filename, 'wb') as outfile: + pickle.dump((model_object, fit_object), outfile) + if not run_silent: + print ('Saved', model_object.model_type,'model fit results for event', + model_object.event.name) + +def load_results(filename, run_silent=True): + """ + Deserializes model and fit objects from a file and prints model information. + + Parameters: + - filename: String specifying the path to the pickle file containing the + pyLIMA model and fit serialized objects. + + Returns: + - model_object: The deserialized model object. + - fit_object: The deserialized fit object associated with the model object. + + Prints the model type and the event name of the model object. + """ + with open(filename, 'rb') as infile: + model_object, fit_object = pickle.load(infile) + if not run_silent: + print (model_object.model_type,'model fit results loaded for event', + model_object.event.name) + return model_object, fit_object From 15abc6a50ab418ebf8f519d8fb9137f81dce45cd Mon Sep 17 00:00:00 2001 From: ytsapras Date: Thu, 15 Feb 2024 16:05:06 +0100 Subject: [PATCH 4/5] fixed pyLIMA/tests/test_toolbox.py:5:28: F811 Redefinition of unused from line 4 --- pyLIMA/tests/test_toolbox.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyLIMA/tests/test_toolbox.py b/pyLIMA/tests/test_toolbox.py index 0ae8d16..756103a 100644 --- a/pyLIMA/tests/test_toolbox.py +++ b/pyLIMA/tests/test_toolbox.py @@ -2,7 +2,6 @@ import unittest from unittest.mock import MagicMock, Mock from pyLIMA.toolbox import brightness_transformation -from pyLIMA.toolbox import brightness_transformation from pyLIMA.toolbox.examine_lightcurve import PointBrowser # Adjust the import path as needed from pyLIMA.toolbox.bin_lightcurve import weighted_mean, get_sigma, bin_data # Adjust the import path as needed From 842ade328cd8e8a3c757b84641adc0ca6b3c9e4b Mon Sep 17 00:00:00 2001 From: ytsapras Date: Thu, 15 Feb 2024 16:14:36 +0100 Subject: [PATCH 5/5] fixed exception warning --- pyLIMA/toolbox/examine_lightcurve.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyLIMA/toolbox/examine_lightcurve.py b/pyLIMA/toolbox/examine_lightcurve.py index 26ecf97..004018f 100644 --- a/pyLIMA/toolbox/examine_lightcurve.py +++ b/pyLIMA/toolbox/examine_lightcurve.py @@ -100,6 +100,7 @@ def update(self): self.fig.canvas.draw() if __name__ == '__main__': + print('Assumes the first line of the .csv file contains the column headings.') if len(argv) != 2: print('Usage: python examine_lightcurve.py /path/to/.csv') exit() @@ -110,7 +111,7 @@ def update(self): try: hjd, mag, merr = np.loadtxt(datafile, delimiter=',', skiprows=1, usecols=(0, 1, 2), dtype=float, unpack=True) - except: + except ValueError: hjd, mag, merr = np.loadtxt(datafile, delimiter=',', skiprows=1, usecols=(0, 2, 3), dtype=float, unpack=True) idx = np.arange(len(hjd))