Resolving conflicts after merging pspl branch

README.md

@@ -3,7 +3,7 @@
 <dl>MulensModel is package for modeling microlensing (or &mu;-lensing) 
 events. </dl>
 
-It is still under development. [Latest release: 1.1.0](https://github.com/rpoleski/MulensModel/releases/latest)
+It is still under development. [Latest release: 1.2.0](https://github.com/rpoleski/MulensModel/releases/latest)
 
 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 can then be input into an arbitrary likelihood function to find the best fit parameters.
 

docs/source/conf.py

@@ -58,9 +58,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u'1.1'
+version = u'1.2'
 # The full version, including alpha/beta/rc tags.
-release = u'1.1.0'
+release = u'1.2.0'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/source/install.rst

@@ -6,11 +6,11 @@ How to install?
 3. Unpack the archive.
 4. Add the path to the unpack directory to the ``PYTHONPATH``, e.g., if you've extracted the archive in your home directory (``/home/USER_NAME/``) in tcsh::
 
-    setenv PYTHONPATH /home/USER_NAME/MulensModel-1.1.0/source\:$PYTHONPATH
+    setenv PYTHONPATH /home/USER_NAME/MulensModel-1.2.0/source\:$PYTHONPATH
 
 in bash::
 
-    export PYTHONPATH=/home/USER_NAME/MulensModel-1.1.0/source:$PYTHONPATH
+    export PYTHONPATH=/home/USER_NAME/MulensModel-1.2.0/source:$PYTHONPATH
 
 In order to have this command invoked every time you open a terminal, please add this command to your startup file (``~/.cshrc``, ``~/.bashrc``, ``~/.profile``, or similar). If you didn't have ``PYTHONPATH`` defined before, then skip the last part of the above commands.
 

examples/use_cases/use_case_25_fluxes.py

@@ -22,11 +22,12 @@ def chi2_fun(theta, event, parameters_to_fit):
 def get_color_constraint(event):
     """
     Calculate the color constraint chi2 penalty
+
     KMT = *int*, dataset number for KMTC I-band
     Spitzer = *int*, dataset number for Spitzer
     """
-    KMT=0
-    Spitzer=9
+    KMT = 0
+    Spitzer = 9
 
     # Color constraint for OB161195 (I_KMT - L_Spitzer)
     (source_color, sigma) = (0.78 0.03)
@@ -38,6 +39,7 @@ def get_color_constraint(event):
 
     return (color - source_color)**2 / sigma**2
 
+# Add data from OB161195
 datasets = []
 filenames = ['KCT01I.dat', 'KCT41I.dat', 'KCT42I.dat', 'KSA01I.dat', 'KSA41I.dat',
              'KSA42I.dat', 'KSS01I.dat', 'KSS41I.dat', 'KSS42I.dat', 

source/MulensModel/event.py

@@ -1,5 +1,5 @@
 import numpy as np
-from math import fsum
+from math import fsum, log
 from astropy.coordinates import SkyCoord
 import astropy.units as u
 
@@ -8,6 +8,7 @@
 from MulensModel.mulensdata import MulensData
 from MulensModel.model import Model
 from MulensModel.coordinates import Coordinates
+from MulensModel.trajectory import Trajectory
 
 
 class Event(object):
@@ -217,15 +218,7 @@ def get_chi2_for_dataset(self, index_dataset, fit_blending=None):
         else:
             self.fit.fit_fluxes()
 
-        if dataset.input_fmt == "mag":
-            data = dataset.mag
-            err_data = dataset.err_mag
-        elif dataset.input_fmt == "flux":
-            data = dataset.flux
-            err_data = dataset.err_flux
-        else:
-            raise ValueError('Unrecognized data format: {:}'.format(
-                    dataset.input_fmt))
+        (data, err_data) = dataset.data_and_err_in_input_fmt()
 
         model = self.fit.get_input_format(data=dataset)
         diff = data - model
@@ -284,10 +277,155 @@ def get_chi2_per_point(self, fit_blending=None):
                 masked_model = self.fit.get_flux(data=dataset)[mask]
                 diff[mask] = dataset.flux[mask] - masked_model
                 err_data[mask] = dataset.err_flux[mask]
+
             chi2_per_point.append((diff/err_data)**2)
 
         return chi2_per_point
 
+    def chi2_gradient(self, parameters, fit_blending=None):
+        """
+        Calculate chi^2 gradient (also called Jacobian), i.e.,
+        :math:`d chi^2/d parameter`.
+
+        Parameters :
+            parameters: *str* or *list*, required
+                Parameters with respect to which gradient is calculated.
+                Currently accepted parameters are: ``t_0``, ``u_0``, ``t_eff``,
+                ``t_E``, ``pi_E_N``, and ``pi_E_E``. The parameters for
+                which you request gradient must be defined in py:attr:`~model`.
+
+            fit_blending: *boolean*, optional
+                Are we fitting for blending flux? If not then blending flux is
+                fixed to 0.  Default is the same as
+                :py:func:`MulensModel.fit.Fit.fit_fluxes()`.
+
+        Returns :
+            gradient: *float* or *np.ndarray*
+                chi^2 gradient
+        """
+        if not isinstance(parameters, list):
+            parameters = [parameters]
+        gradient = {param: 0 for param in parameters}
+
+        if self.model.n_lenses != 1:
+            raise NotImplementedError('Event.chi2_gradient() works only ' +
+                'single lens models currently')
+        as_dict = self.model.parameters.as_dict()
+        if 'rho' in as_dict or 't_star' in as_dict:
+            raise NotImplementedError('Event.chi2_gradient() is not working ' +
+                'for finite source models yet')
+
+        # Define a Fit given the model and perform linear fit for fs and fb
+        self.fit = Fit(
+            data=self.datasets, magnification=self.model.data_magnification)
+        if fit_blending is not None:
+            self.fit.fit_fluxes(fit_blending=fit_blending)
+        else:
+            self.fit.fit_fluxes()
+
+        for (i, dataset) in enumerate(self.datasets):
+            ## Original
+            (data, err_data) = dataset.data_and_err_in_input_fmt()
+            factor = data - self.fit.get_input_format(data=dataset)
+            factor *= -2. / err_data**2
+            if dataset.input_fmt == 'mag':
+                factor *= -2.5 / (log(10.) * Utils.get_flux_from_mag(data))
+            factor *= self.fit.flux_of_sources(dataset)[0]
+
+            ## np.log
+            #(data, err_data) = dataset.data_and_err_in_input_fmt()
+            #factor = data - self.fit.get_input_format(data=dataset)
+            #factor *= -2. / err_data**2
+            #if dataset.input_fmt == 'mag':
+            #    factor *= -2.5 / (np.log(10.) * Utils.get_flux_from_mag(data))
+            #factor *= self.fit.flux_of_sources(dataset)[0]
+
+            ## Fluxes
+            #f_source = self.fit.flux_of_sources(dataset)[0]
+            #f_blend = self.fit.blending_flux(dataset)
+            #model_flux = (f_source *
+            #              self.model.get_data_magnification(dataset) + f_blend)
+            #factor = (-2. * f_source * 
+            #           (dataset.flux - model_flux) / dataset.err_flux**2)
+
+            kwargs = {}
+            if dataset.ephemerides_file is not None:
+                kwargs['satellite_skycoord'] = dataset.satellite_skycoord
+            trajectory = Trajectory(dataset.time, self.model.parameters, 
+                    self.model.get_parallax(), self.coords, **kwargs)
+            u_2 = trajectory.x**2 + trajectory.y**2
+            u_ = np.sqrt(u_2)
+            d_A_d_u = -8. / (u_2 * (u_2 + 4) * np.sqrt(u_2 + 4))
+            factor *= d_A_d_u
+
+            factor_d_x_d_u = (factor * trajectory.x / u_)[dataset.good]
+            sum_d_x_d_u = np.sum(factor_d_x_d_u)
+            factor_d_y_d_u = (factor * trajectory.y / u_)[dataset.good]
+            sum_d_y_d_u = np.sum(factor_d_y_d_u)
+            dt = dataset.time[dataset.good] - as_dict['t_0']
+
+            # Exactly 2 out of (u_0, t_E, t_eff) must be defined and
+            # gradient depends on which ones are defined. 
+            if 't_eff' not in as_dict:
+                t_E = as_dict['t_E'].to(u.day).value
+                if 't_0' in parameters:
+                    gradient['t_0'] += -sum_d_x_d_u / t_E
+                if 'u_0' in parameters:
+                    gradient['u_0'] += sum_d_y_d_u
+                if 't_E' in parameters:
+                    gradient['t_E'] += np.sum(factor_d_x_d_u * -dt / t_E**2)
+            elif 't_E' not in as_dict:
+                t_eff = as_dict['t_eff'].to(u.day).value
+                if 't_0' in parameters:
+                    gradient['t_0'] += -sum_d_x_d_u * as_dict['u_0'] / t_eff
+                if 'u_0' in parameters:
+                    gradient['u_0'] += sum_d_y_d_u + np.sum(
+                            factor_d_x_d_u * dt / t_eff)
+                if 't_eff' in parameters:
+                    gradient['t_eff'] += np.sum(factor_d_x_d_u * -dt *
+                            as_dict['u_0'] / t_eff**2)
+            elif 'u_0' not in as_dict:
+                t_E = as_dict['t_E'].to(u.day).value
+                t_eff = as_dict['t_eff'].to(u.day).value
+                if 't_0' in parameters:
+                    gradient['t_0'] += -sum_d_x_d_u / t_E
+                if 't_E' in parameters:
+                    gradient['t_E'] += (np.sum(factor_d_x_d_u * dt) -
+                            sum_d_y_d_u * t_eff) / t_E**2
+                if 't_eff' in parameters:
+                    gradient['t_eff'] += sum_d_y_d_u / t_E
+            else:
+                raise KeyError('Something is wrong with ModelParameters in ' +
+                    'Event.chi2_gradient():\n', as_dict)
+
+            # Below we deal with parallax only.
+            if 'pi_E_N' in parameters or 'pi_E_E' in parameters:
+                parallax = {'earth_orbital': False, 'satellite': False,
+                        'topocentric': False}
+                trajectory_no_piE = Trajectory(dataset.time, 
+                    self.model.parameters, parallax, self.coords,
+                    **kwargs)
+                dx = (trajectory.x - trajectory_no_piE.x)[dataset.good]
+                dy = (trajectory.y - trajectory_no_piE.y)[dataset.good]
+                delta_E = dx * as_dict['pi_E_E'] + dy * as_dict['pi_E_N']
+                delta_N = dx * as_dict['pi_E_N'] - dy * as_dict['pi_E_E']
+                det = as_dict['pi_E_N']**2 + as_dict['pi_E_E']**2
+
+                if 'pi_E_N' in parameters:
+                    gradient['pi_E_N'] += np.sum(
+                        factor_d_x_d_u * delta_N + factor_d_y_d_u * delta_E)
+                    gradient['pi_E_N'] /= det
+                if 'pi_E_E' in parameters:
+                    gradient['pi_E_E'] += np.sum(
+                        factor_d_x_d_u * delta_E - factor_d_y_d_u * delta_N)
+                    gradient['pi_E_E'] /= det
+
+        if len(parameters) == 1:
+            out = gradient[parameters[0]]
+        else:
+            out = np.array([gradient[p] for p in parameters])
+        return out
+
     def get_ref_fluxes(self, data_ref=None):
         """
         Get source and blending fluxes for the reference dataset. See

source/MulensModel/model.py

@@ -337,6 +337,18 @@ def parallax(
             if topocentric is not None:
                 self._parallax['topocentric'] = topocentric
 
+    def get_parallax(self):
+        """
+        Returns *dict* that specifies the types of the microlensing parallax
+        that are included in calculations.
+
+        Returns :
+            parallax: *dict*
+                For keys ``'earth_orbital'``, ``'satellite'``,
+                and ``'topocentric'`` returns *bool*.
+        """
+        return self._parallax
+
     def plot_magnification(
             self, times=None, t_range=None, t_start=None, t_stop=None, dt=None,
             n_epochs=None, subtract_2450000=False, subtract_2460000=False,

source/MulensModel/mulensdata.py

@@ -284,6 +284,30 @@ def n_epochs(self):
         """
         return self._n_epochs
 
+    def data_and_err_in_input_fmt(self):
+        """
+        Gives photometry in input format (mag or flux).
+
+        Returns :
+            data: *np.ndarray*
+                Magnitudes or fluxes
+
+            data_err: *np.ndarray*
+                Uncertainties of magnitudes or of fluxes
+
+        """
+        if self.input_fmt == "mag":
+            data = self.mag
+            err_data = self.err_mag
+        elif dataset.input_fmt == "flux":
+            data = self.flux
+            err_data = self.err_flux
+        else:
+            raise ValueError('Unrecognized data format: {:}'.format(
+                    self.input_fmt))
+
+        return (data, err_data)
+
     @property
     def satellite_skycoord(self):
         """

source/MulensModel/tests/test_Event.py

@@ -13,6 +13,8 @@
 
 SAMPLE_FILE_01 = os.path.join(MulensModel.MODULE_PATH, 
                                     "data", "phot_ob08092_O4.dat")
+SAMPLE_FILE_02 = os.path.join(MulensModel.MODULE_PATH, 
+                                    "data", "ob140939_OGLE.dat")
 
 def test_event_get_chi2_1():
     """basic unit test on ob08092 OGLE-IV data"""
@@ -72,6 +74,33 @@ def test_event_get_chi2_2():
     chi2_3 = ev.get_chi2_for_dataset(1)
     np.testing.assert_almost_equal(chi2_3, answer)
 
+def test_event_get_chi2_3():
+    """test on ob08092 OGLE-IV data - MulensData.good & MulensData.bad"""
+    t_0 = 5379.57091
+    u_0 = 0.52298
+    t_E = 17.94002
+
+    bad = np.zeros(383, dtype='bool')
+    bad[300:350] = True
+    data_1 = MulensData(file_name=SAMPLE_FILE_01, bad=bad)
+    data_2 = MulensData(file_name=SAMPLE_FILE_01, good=~bad)
+
+    ev = Event()
+    mod = Model({'t_0': t_0, 'u_0': u_0, 't_E': t_E})
+    mod.set_datasets([data_1])
+    ev.model = mod
+    ev.datasets = [data_1]
+    chi2 = ev.get_chi2()
+    np.testing.assert_almost_equal(float(chi2), 343.46567, decimal=4, 
+                                   err_msg='problem in resulting chi2')
+
+    mod.set_datasets([data_2])
+    ev.model = mod
+    ev.datasets = [data_2]
+    chi2 = ev.get_chi2()
+    np.testing.assert_almost_equal(float(chi2), 343.46567, decimal=4,
+                                    err_msg='problem in resulting chi2')
+
 def test_event_get_chi2_double_source_simple():
     """basic test on ob08092 OGLE-IV data with added second source
     Note that currently this test hacks into internal functions of 
@@ -156,3 +185,33 @@ def test_event_init_1(self):
     def test_event_init_2(self):
         with self.assertRaises(TypeError):
             ev = Event(datasets='some_string')
+
+def test_event_chi2_gradient():
+    # fs = 11.0415734, fb = 0.0 
+    parameters_1 = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87}
+    params_1 = ['t_0', 'u_0', 't_E']
+    gradient_1 = {'t_0': 236.206598, 'u_0': 101940.249,
+                't_E': -1006.88678}
+    test_1 = (parameters_1, params_1, gradient_1)
+
+    parameters_2 = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87,
+                    'pi_E_N': -0.248, 'pi_E_E': 0.234} # This model also
+                    # used fluxes given above.
+    params_2 = ['t_0', 'u_0', 't_E', 'pi_E_N', 'pi_E_E', 'f_source', 'f_blend']
+    gradient_2 = {'t_0': 568.781786, 'u_0': 65235.3513, 't_E': -491.782005,
+                  'pi_E_N': -187878.357, 'pi_E_E': 129162.927,
+                  'f_source': -83124.5869, 'f_blend': -78653.242}
+    test_2 = (parameters_2, params_2, gradient_2)
+    # We're not applying the test above, yet. See 'for' loop below.
+
+    data = MulensData(file_name=SAMPLE_FILE_02)
+    kwargs = {'datasets': [data], 'coords': '17:47:12.25 −21:22:58.7'}
+
+    for test in [test_1]:#, test_2]:
+        (parameters, params, gradient) = test
+        event = Event(model=Model(parameters), **kwargs)
+        result = event.chi2_gradient(params, fit_blending=False)
+
+        reference = np.array([gradient[key] for key in  params])
+        np.testing.assert_almost_equal(reference/result, 1., decimal=1)
+

source/MulensModel/trajectory.py

@@ -87,9 +87,9 @@ def __init__(self, times, parameters, parallax=None,
 
     def get_xy(self):
         """
-        For a given set of parameters (a
-        :py:class:`~MulensModel.modelparameters.ModelParameters`
-        object), calculate the xy position of the source.
+        For a given set of parameters
+        (a :py:class:`~MulensModel.modelparameters.ModelParameters` object),
+        calculate the xy position of the source.
         """
         # Calculate the position of the source
         vector_tau = ((self.times - self.parameters.t_0) /

source/MulensModel/version.py

@@ -1,2 +1,2 @@
 
-__version__ = "1.1.5"
+__version__ = "1.2.0"