Toolbox additions yt

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

pyLIMA/tests/test_toolbox.py

@@ -1,6 +1,9 @@
 import numpy as np
+import unittest
+from unittest.mock import MagicMock, Mock
 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 +40,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)

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 <input_lightcurve_file> [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()
+