diff --git a/README.md b/README.md index 76035621..3ec4c4d0 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@
MulensModel is package for modeling microlensing (or μ-lensing) events.
-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. diff --git a/docs/source/conf.py b/docs/source/conf.py index 14703cde..332b0cd1 100644 --- a/docs/source/conf.py +++ b/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. diff --git a/docs/source/install.rst b/docs/source/install.rst index 9666f118..70bc6865 100644 --- a/docs/source/install.rst +++ b/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. diff --git a/examples/use_cases/use_case_25_fluxes.py b/examples/use_cases/use_case_25_fluxes.py index a2f16c0a..f0f51f8d 100644 --- a/examples/use_cases/use_case_25_fluxes.py +++ b/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', diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 3fdabcb5..579ec29c 100644 --- a/source/MulensModel/event.py +++ b/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 diff --git a/source/MulensModel/model.py b/source/MulensModel/model.py index 25f3122b..e6368487 100644 --- a/source/MulensModel/model.py +++ b/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, diff --git a/source/MulensModel/mulensdata.py b/source/MulensModel/mulensdata.py index a93bd39a..8c755637 100644 --- a/source/MulensModel/mulensdata.py +++ b/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): """ diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index a93b2e43..78e02071 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/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) + diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index b4ee1ee6..973ca765 100644 --- a/source/MulensModel/trajectory.py +++ b/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) / diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py index c75be76c..bdc2ab81 100644 --- a/source/MulensModel/version.py +++ b/source/MulensModel/version.py @@ -1,2 +1,2 @@ -__version__ = "1.1.5" +__version__ = "1.2.0"