Merge branch 'rpoleski:master' into fix-interp2d-import

.gitignore

@@ -68,3 +68,4 @@ source/MulensModel/*.so
 
 # Coverage report:
 coverage.xml
+/source/MulensModel/htmlcov/

README.md

@@ -6,7 +6,7 @@
 
 [**Detailed documentation: https://rpoleski.github.io/MulensModel/**](https://rpoleski.github.io/MulensModel/)
 
-[Latest release: 2.18.0](https://github.com/rpoleski/MulensModel/releases/latest) and we're working on further developing the code.
+[Latest release: 2.19.0](https://github.com/rpoleski/MulensModel/releases/latest) and we're working on further developing the code.
 
 MulensModel can generate a microlensing light curve for a given set of microlensing parameters, fit that light curve to some data, and return a chi2 value. That chi2 (and its gradient in some cases) can then be input into an arbitrary likelihood function to find the best-fit parameters.
 
@@ -62,4 +62,4 @@ If you want to contribute to MulensModel, then please see [this file](CONTRIBUTI
 ![example workflow](https://github.com/rpoleski/MulensModel/actions/workflows/python-app.yml/badge.svg)
 [![Coverage Status](https://codecov.io/gh/rpoleski/MulensModel/branch/master/graph/badge.svg)](https://codecov.io/gh/rpoleski/MulensModel)
 
-file revised Sep 2023
+file revised Dec 2023

data/unit_test_files/ob140939_Spitzer_ref_v2.dat

@@ -0,0 +1,31 @@
+6814.07844 15.042 0.001    -0.023   -19.345 1.296867   -0.7731    0.7371
+6815.03086 15.039 0.001    -0.005    -4.165 1.319117   -0.7277    0.7329
+6815.82782 15.016 0.001    -0.009    -7.327 1.338840   -0.6898    0.7294
+6816.77402 15.028 0.001     0.026    22.175 1.363561   -0.6449    0.7253
+6817.58341 14.974 0.001    -0.009    -7.510 1.385806   -0.6064    0.7219
+6818.73166 14.932 0.001    -0.021   -18.234 1.419008   -0.5520    0.7171
+6819.53932 14.909 0.001    -0.023   -22.244 1.443406   -0.5138    0.7138
+6820.52634 14.889 0.001    -0.018   -17.159 1.474198   -0.4671    0.7098
+6821.48308 14.870 0.001    -0.011   -10.516 1.504799   -0.4220    0.7060
+6822.65891 14.856 0.001     0.005     4.743 1.542871   -0.3666    0.7014
+6823.46861 14.810 0.001    -0.020   -19.788 1.568958   -0.3285    0.6983
+6824.70223 14.790 0.001    -0.011   -10.841 1.607594   -0.2705    0.6937
+6825.49204 14.774 0.001    -0.009    -9.173 1.631038   -0.2335    0.6908
+6827.09946 14.745 0.001    -0.008    -7.411 1.673595   -0.1582    0.6850
+6827.80323 14.740 0.001    -0.001    -0.745 1.689287   -0.1253    0.6825
+6828.52669 14.723 0.001    -0.008    -7.725 1.703095   -0.0915    0.6800
+6830.38779 14.720 0.001     0.005     4.732 1.725832   -0.0047    0.6737
+6831.42212 14.721 0.001     0.009     8.424 1.729582    0.0434    0.6703
+6832.59984 14.703 0.001    -0.012   -11.978 1.725780    0.0982    0.6666
+6834.21402 14.729 0.001     0.001     0.721 1.707234    0.1730    0.6616
+6835.29855 14.737 0.001    -0.006    -5.597 1.687140    0.2232    0.6583
+6836.28071 14.762 0.001     0.003     2.840 1.664599    0.2686    0.6554
+6836.97713 14.774 0.001     0.002     2.114 1.646585    0.3007    0.6534
+6838.16978 14.809 0.001     0.012    11.890 1.612805    0.3557    0.6501
+6839.03530 14.844 0.001     0.027    26.868 1.586669    0.3956    0.6477
+6839.94137 14.841 0.001     0.002     2.256 1.558475    0.4372    0.6453
+6841.10244 14.877 0.001     0.010     9.571 1.521894    0.4905    0.6423
+6842.17330 14.924 0.001     0.030    25.773 1.488399    0.5396    0.6396
+6843.58489 14.948 0.001     0.018    15.543 1.445512    0.6042    0.6361
+6844.44113 14.973 0.001     0.021    17.904 1.420527    0.6434    0.6341
+6845.70596 15.007 0.001     0.024    21.116 1.385364    0.7011    0.6311

examples/example_17_1L2S_plotting.py

@@ -94,4 +94,9 @@ def generate_dataset(time, flux_1, flux_2, blend_flux, flux_err,
 my_event.plot_model(zorder=10, color='black')
 my_event.plot_data()
 
+# Plot the source trajectories
+plt.figure()
+plt.title('Model Source Trajectories')
+my_model.plot_trajectory()
+
 plt.show()

source/MulensModel/fitdata.py

@@ -647,16 +647,11 @@ def get_dataset_trajectory(self):
             trajectory: :py:class:`~MulensModel.trajectory.Trajectory`
                 Trajectory for given dataset.
         """
-        if self.dataset.ephemerides_file is None:
-            return self.model.get_trajectory(self.dataset.time)
-        else:
-            kwargs_ = {
-                'times': self.dataset.time,
-                'parallax': self.model.get_parallax(),
-                'coords': self.model.coords,
-                'satellite_skycoord': self.dataset.satellite_skycoord}
+        kwargs = dict()
+        if self.dataset.ephemerides_file is not None:
+            kwargs['satellite_skycoord'] = self.dataset.satellite_skycoord
 
-            return mm.Trajectory(parameters=self.model.parameters, **kwargs_)
+        return self.model.get_trajectory(self.dataset.time, **kwargs)
 
     def get_d_A_d_u_for_PSPL_model(self):
         """

source/MulensModel/model.py

@@ -599,6 +599,7 @@ def plot_trajectory(
             times = self.set_times(
                 t_range=t_range, t_start=t_start, t_stop=t_stop, dt=dt,
                 n_epochs=n_epochs)
+
         if satellite_skycoord is None:
             satellite_skycoord = self.get_satellite_coords(times)
         else:
@@ -755,21 +756,42 @@ def _plot_source_for_trajectory(self, trajectory, **kwargs):
             for (x, y) in zip(trajectory.x, trajectory.y):
                 axis.add_artist(plt.Circle((x, y), **kwargs))
 
-    def get_trajectory(self, times):
+    def get_trajectory(self, times, satellite_skycoord=None):
         """
         Get the source trajectory for the given set of times.
 
         Parameters :
             times:  *np.ndarray*, *list of floats*, or *float*
                 Epochs for which source positions are requested.
 
-        Returns : A `:py:class:`~MulensModel.trajectory.Trajectory` object.
+            satellite_skycoord: *astropy.SkyCoord*
+                Allows the user to specify that the trajectory is calculated
+                for a satellite. If *astropy.SkyCoord* object is provided,
+                then these are satellite positions for all epochs.
+                See also :py:func:`get_satellite_coords()`
+
+        Returns : A `:py:class:`~MulensModel.trajectory.Trajectory` object. If
+            n_sources > 1, returns a tuple of
+            `:py:class:`~MulensModel.trajectory.Trajectory`s
 
         """
+        if satellite_skycoord is None:
+            satellite_skycoord = self.get_satellite_coords(times)
+
         kwargs_ = {
             'times': times, 'parallax': self._parallax, 'coords': self._coords,
-            'satellite_skycoord': self.get_satellite_coords(times)}
-        return Trajectory(parameters=self.parameters, **kwargs_)
+            'satellite_skycoord': satellite_skycoord}
+        if self.n_sources == 1:
+            return Trajectory(parameters=self.parameters, **kwargs_)
+        elif self.n_sources == 2:
+            trajectory_1 = Trajectory(
+                parameters=self.parameters.source_1_parameters, **kwargs_)
+            trajectory_2 = Trajectory(
+                parameters=self.parameters.source_2_parameters, **kwargs_)
+            return (trajectory_1, trajectory_2)
+        else:
+            raise NotImplementedError(
+                "only 1 or 2 sources allowed here at this point")
 
     def set_times(
             self, t_range=None, t_start=None, t_stop=None, dt=None,

source/MulensModel/mulensdata.py

@@ -450,17 +450,26 @@ def _plot_datapoints(
 
         if show_errorbars:
             if np.any(y_err[self.good] < 0.):
-                warnings.warn("Cannot plot errorbars with negative values. "
-                              "Skipping dataset: " + self._get_name())
-                return
+                ind_neg_err = np.where((y_err < 0.) & self.good)
+                warnings.warn(
+                    "Some points have errorbars with negative values. " +
+                    "Setting to zero. \n" +
+                    "Dataset: " + self._get_name() +
+                    "\nEpochs: {0}".format(self.time[ind_neg_err]))
+                y_err[ind_neg_err] = 0.
+
             container = self._plt_errorbar(time_good, y_good,
                                            y_err[self.good], properties)
             if show_bad:
                 if np.any(y_err[self.bad] < 0.):
+                    ind_neg_err = np.where((y_err < 0.) & self.bad)
                     warnings.warn(
-                        "Cannot plot errorbars with negative values (bad "
-                        "data). Skipping dataset: " + self._get_name())
-                    return
+                        "Some (bad data) points have errorbars with " +
+                        "negative values. Setting to zero. \n" +
+                        "Dataset: " + self._get_name() +
+                        "\nEpochs: {0}".format(self.time[ind_neg_err]))
+                    y_err[ind_neg_err] = 0.
+
                 if not ('color' in properties_bad or 'c' in properties_bad):
                     properties_bad['color'] = container[0].get_color()
 

source/MulensModel/tests/test_FitData.py

@@ -12,10 +12,10 @@
 dir_4 = join(dir_2, 'fspl_derivs')
 
 SAMPLE_FILE_02 = join(dir_1, 'ob140939_OGLE.dat')  # HJD'
-SAMPLE_FILE_02_REF = join(dir_2, 'ob140939_OGLE_ref_v1.dat')  # HJD'
+SAMPLE_FILE_02_REF = join(dir_2, 'ob140939_OGLE_ref_v2.dat')  # HJD'
 SAMPLE_FILE_03 = join(dir_1, 'ob140939_Spitzer.dat')  # HJD'
 SAMPLE_FILE_03_EPH = join(dir_3, 'Spitzer_ephemeris_01.dat')  # UTC
-SAMPLE_FILE_03_REF = join(dir_2, 'ob140939_Spitzer_ref_v1.dat')  # HJD'
+SAMPLE_FILE_03_REF = join(dir_2, 'ob140939_Spitzer_ref_v2.dat')  # HJD'
 SAMPLE_FILE_04_WF = join(mm.DATA_PATH, 'WFIRST_1827.dat')
 SAMPLE_FILE_FSPL_51 = join(dir_4, 'fort.51')
 SAMPLE_FILE_FSPL_61 = join(dir_4, 'fort.61')
@@ -1006,6 +1006,31 @@ def test_FSPLDerivs_get_satellite_coords():
     np.testing.assert_almost_equal(result_2.dec.value, dec_2, decimal=3)
 
 
+def test_get_trajectory_1L2S_satellite_parallax():
+    """test parallax calculation with Spitzer data"""
+    model_parameters = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87,
+                        'pi_E_N': -0.248, 'pi_E_E': 0.234,
+                        't_0_par': 2456836.2}
+    coords = "17:47:12.25 -21:22:58.2"
+
+    model_with_par = mm.Model(model_parameters, coords=coords)
+    model_with_par.parallax(satellite=True, earth_orbital=True,
+                            topocentric=False)
+
+    data_Spitzer = mm.MulensData(
+        file_name=SAMPLE_FILE_03, ephemerides_file=SAMPLE_FILE_03_EPH)
+
+    fit = mm.FitData(model=model_with_par, dataset=data_Spitzer)
+
+    ref_Spitzer = np.loadtxt(SAMPLE_FILE_03_REF, unpack=True)
+
+    trajectory = fit.get_dataset_trajectory()
+
+    ratio_x = trajectory.x / ref_Spitzer[6]
+    ratio_y = trajectory.y / ref_Spitzer[7]
+    np.testing.assert_almost_equal(ratio_x, 1., decimal=2)
+    np.testing.assert_almost_equal(ratio_y, 1., decimal=3)
+
 # Tests to add:
 #
 # test get_chi2_gradient(), chi2_gradient:

source/MulensModel/tests/test_Model.py

@@ -3,10 +3,17 @@
 from math import isclose
 import unittest
 from astropy import units as u
+import os.path
 import warnings
 
 import MulensModel as mm
 
+dir_2 = os.path.join(mm.DATA_PATH, 'unit_test_files')
+dir_3 = os.path.join(mm.DATA_PATH, 'ephemeris_files')
+SAMPLE_FILE_02_REF = os.path.join(dir_2, 'ob140939_OGLE_ref_v2.dat')  # HJD'
+SAMPLE_FILE_03_EPH = os.path.join(dir_3, 'Spitzer_ephemeris_01.dat')  # UTC
+SAMPLE_FILE_03_REF = os.path.join(dir_2, 'ob140939_Spitzer_ref_v2.dat')  # HJD'
+
 
 def test_n_lenses():
     """check n_lenses property"""
@@ -765,14 +772,99 @@ def test_xallarap_at_t_0_plus_half_of_period_7_eccentric():
     almost(expected, model.get_magnification(t_0+d_time))
 
 
+class TestGetTrajectory(unittest.TestCase):
+
+    def setUp(self):
+        # Parallax model parameters
+        self.model_parameters_par = {
+            't_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87,
+            'pi_E_N': -0.248, 'pi_E_E': 0.234, 't_0_par': 2456836.2}
+        self.coords = "17:47:12.25 -21:22:58.2"
+
+        self.model_with_par = mm.Model(
+            self.model_parameters_par, coords=self.coords)
+        self.model_with_par.parallax(satellite=True, earth_orbital=True,
+                                     topocentric=False)
+
+        self.ref_OGLE = np.loadtxt(SAMPLE_FILE_02_REF, unpack=True)
+        self.times_OGLE = self.ref_OGLE[0] + 2450000.
+        self.ref_Spitzer = np.loadtxt(SAMPLE_FILE_03_REF, unpack=True)
+        self.ephemerides_file = SAMPLE_FILE_03_EPH
+        self.times_spz = self.ref_Spitzer[0] + 2450000.
+
+    def test_1L1S(self):
+        # straight-up trajectory for static point-lens model
+        t_0 = self.model_parameters_par['t_0']
+        u_0 = self.model_parameters_par['u_0']
+        t_E = self.model_parameters_par['t_E']
+        times = np.arange(t_0 - t_E, t_0 + t_E, 1.)
+
+        model = mm.Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E})
+        trajectory = model.get_trajectory(times)
+        y = np.ones(len(times)) * u_0
+        x = (times - t_0) / t_E
+        ratio_x = trajectory.x / x
+        ratio_y = trajectory.y / y
+        np.testing.assert_almost_equal(ratio_x, 1.)
+        np.testing.assert_almost_equal(ratio_y, 1.)
+
+    def test_1L1S_annual_parallax(self):
+        """case with annual parallax"""
+        trajectory = self.model_with_par.get_trajectory(
+            self.times_OGLE)
+
+        ratio_x = trajectory.x / self.ref_OGLE[6]
+        ratio_y = trajectory.y / self.ref_OGLE[7]
+        np.testing.assert_almost_equal(ratio_x, [1.] * len(ratio_x), decimal=3)
+        np.testing.assert_almost_equal(ratio_y, [1.] * len(ratio_y), decimal=3)
+
+    def test_1L1S_satellite_parallax_1(self):
+        """Case with satellite parallax (check test_Model_Parallax.py)"""
+        satellite_skycoord_obj = mm.SatelliteSkyCoord(
+            ephemerides_file=self.ephemerides_file)
+        satellite_skycoord = satellite_skycoord_obj.get_satellite_coords(
+            self.times_spz)
+        trajectory = self.model_with_par.get_trajectory(
+            self.times_spz, satellite_skycoord=satellite_skycoord)
+
+        ratio_x = trajectory.x / self.ref_Spitzer[6]
+        ratio_y = trajectory.y / self.ref_Spitzer[7]
+        np.testing.assert_almost_equal(ratio_x, [1.] * len(ratio_x), decimal=2)
+        np.testing.assert_almost_equal(ratio_y, [1.] * len(ratio_y), decimal=3)
+
+    def test_1L1S_satellite_parallax_2(self):
+        """Case with satellite parallax (check test_Model_Parallax.py)"""
+        model_with_sat_par = mm.Model(
+            self.model_parameters_par, ephemerides_file=self.ephemerides_file,
+            coords=self.coords)
+        trajectory = model_with_sat_par.get_trajectory(self.times_spz)
+
+        ratio_x = trajectory.x / self.ref_Spitzer[6]
+        ratio_y = trajectory.y / self.ref_Spitzer[7]
+        np.testing.assert_almost_equal(ratio_x, [1.] * len(ratio_x), decimal=2)
+        np.testing.assert_almost_equal(ratio_y, [1.] * len(ratio_y), decimal=3)
+
+    def test_1L2S(self):
+        """Binary source trajectories"""
+        model = mm.Model({
+            't_0_1': 5000., 'u_0_1': 0.005, 'rho_1': 0.001,
+            't_0_2': 5100., 'u_0_2': 0.0003, 't_star_2': 0.03, 't_E': 25.})
+        model_1 = mm.Model(model.parameters.source_1_parameters)
+        model_2 = mm.Model(model.parameters.source_2_parameters)
+
+        time = np.linspace(4900., 5200., 4200)
+
+        (traj_1_1L2S, traj_2_1L2S) = model.get_trajectory(time)
+        traj_1_1L1S = model_1.get_trajectory(time)
+        np.testing.assert_equal(traj_1_1L2S.x, traj_1_1L1S.x)
+        np.testing.assert_equal(traj_1_1L2S.y, traj_1_1L1S.y)
+
+        traj_2_1L1S = model_2.get_trajectory(time)
+        np.testing.assert_equal(traj_2_1L2S.x, traj_2_1L1S.x)
+        np.testing.assert_equal(traj_2_1L2S.y, traj_2_1L1S.y)
+
 # Tests to Add:
 #
-# test get_trajectory:
-#   straight-up trajectory
-#   case with annual parallax (check test_Model_Parallax.py)
-#   case with satellite parallax (check test_Model_Parallax.py)
-#   coords is propagating correctly (check test_Model_Parallax.py)
-#
 # test set_times:
 #   keywords to test:
 #     t_range=None, t_start=None, t_stop=None, dt=None, n_epochs=1000

source/MulensModel/tests/test_Model_Parallax.py

@@ -12,10 +12,10 @@
 dir_3 = os.path.join(mm.DATA_PATH, 'ephemeris_files')
 
 SAMPLE_FILE_02 = os.path.join(dir_1, 'ob140939_OGLE.dat')  # HJD'
-SAMPLE_FILE_02_REF = os.path.join(dir_2, 'ob140939_OGLE_ref_v1.dat')  # HJD'
+SAMPLE_FILE_02_REF = os.path.join(dir_2, 'ob140939_OGLE_ref_v2.dat')  # HJD'
 SAMPLE_FILE_03 = os.path.join(dir_1, 'ob140939_Spitzer.dat')  # HJD'
 SAMPLE_FILE_03_EPH = os.path.join(dir_3, 'Spitzer_ephemeris_01.dat')  # UTC
-SAMPLE_FILE_03_REF = os.path.join(dir_2, 'ob140939_Spitzer_ref_v1.dat')  # HJD'
+SAMPLE_FILE_03_REF = os.path.join(dir_2, 'ob140939_Spitzer_ref_v2.dat')  # HJD'
 SAMPLE_FILE_04_EPH = os.path.join(dir_3, 'Spitzer_ephemeris_03.dat')  # UTC
 SAMPLE_FILE_04_REF = os.path.join(dir_2, 'gets_70.dat')