From dc2e58e8c3e61fa7589c78ca2a8c19a52d623733 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Sun, 1 Apr 2018 06:42:55 -0400 Subject: [PATCH 01/17] starting Event.chi2_gradient() --- source/MulensModel/event.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 4ddabe48..f0622fc0 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -272,6 +272,36 @@ def get_chi2_per_point(self, fit_blending=None): return chi2_per_point + def chi2_gradient(self, parameters): + """ + 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`. + + Returns : + gradient: *float* or *np.ndarray* + chi^2 gradient + """ + if not isinstance(parameters, list): + parameters = [parameters] + if self.model.n_lenses != 1: + raise NotImplementedError('Event.chi2_gradient() works only ' + + 'single lens models currently') + +# check which parameters are in self.model.parameters.as_dict() +# check which subset of (u_0, t_E, t_eff) is used +# run get_chi2_per_point() +# mag vs flux +# dA/du +# parallax treated separately because doesnt depend on parametrization + raise NotImplementedError('NOT FINISHED YET') + def get_ref_fluxes(self, data_ref=None): """ Get source and blending fluxes for the reference dataset. See From d8fca256507a729f4e47462f1d371d5f53a477ad Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Sun, 1 Apr 2018 16:21:36 -0400 Subject: [PATCH 02/17] data_and_err_in_input_fmt() added to MulensData --- source/MulensModel/mulensdata.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/source/MulensModel/mulensdata.py b/source/MulensModel/mulensdata.py index 9ddedf0a..db46b0da 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): """ From e2e692f1a7e5897466e2d1694f396eca159226b7 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Sun, 1 Apr 2018 16:22:51 -0400 Subject: [PATCH 03/17] data_and_err_in_input_fmt() used in Event and some work on chi2 gradient --- source/MulensModel/event.py | 50 +++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index f0622fc0..4479f39c 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 @@ -215,15 +215,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() diff = data - self.fit.get_input_format(data=dataset) chi2 = (diff / err_data)**2 @@ -258,21 +250,13 @@ def get_chi2_per_point(self, fit_blending=None): # Calculate chi^2 given the fit chi2_per_point = [] for (i, dataset) in enumerate(self.datasets): - 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() diff = data - self.fit.get_input_format(data=dataset) chi2_per_point.append((diff/err_data)**2) return chi2_per_point - def chi2_gradient(self, parameters): + def chi2_gradient(self, parameters, fit_blending=None): """ Calculate chi^2 gradient (also called Jacobian), i.e., :math:`d chi^2/d parameter`. @@ -284,6 +268,11 @@ def chi2_gradient(self, parameters): ``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 @@ -293,8 +282,27 @@ def chi2_gradient(self, 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') -# check which parameters are in self.model.parameters.as_dict() + # 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): + (data, err_data) = dataset.data_and_err_in_input_fmt() + factor_1 = data - self.fit.get_input_format(data=dataset) + factor_1 /= err_data**2 + if dataset.input_fmt == 'mag': + factor_1 *= -2.5 / (log(10.) * Utils.get_flux_from_mag(mag)) + factor_1 *= F_S # check which subset of (u_0, t_E, t_eff) is used # run get_chi2_per_point() # mag vs flux From 0b8ccd3c35dc682763221b4c06fa0a9978e11450 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Sun, 1 Apr 2018 18:52:58 -0400 Subject: [PATCH 04/17] Model.get_parallax() --- source/MulensModel/model.py | 12 ++++++++++++ source/MulensModel/trajectory.py | 6 +++--- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/source/MulensModel/model.py b/source/MulensModel/model.py index 3d5d3daa..bb4a0b14 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/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) / From 5564f806cb3881ff994225ccd5ede28e0d046b5d Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Sun, 1 Apr 2018 19:05:55 -0400 Subject: [PATCH 05/17] more work on chi2_gradient() --- source/MulensModel/event.py | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 4479f39c..29bfde65 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -279,6 +279,8 @@ def chi2_gradient(self, parameters, fit_blending=None): """ 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') @@ -302,11 +304,28 @@ def chi2_gradient(self, parameters, fit_blending=None): factor_1 /= err_data**2 if dataset.input_fmt == 'mag': factor_1 *= -2.5 / (log(10.) * Utils.get_flux_from_mag(mag)) - factor_1 *= F_S -# check which subset of (u_0, t_E, t_eff) is used -# run get_chi2_per_point() -# mag vs flux -# dA/du + factor_1 *= self.fit.flux_of_sources(dataset)[0] + + trajectory = Trajectory(dataset.time, self.model.parameters, + self.model.get_parallax, self.coords, + dataset.satellite_skycoord) + 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)) + d_x_d_u = trajectory.x / u + d_y_d_u = trajectory.y / u + + # 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: + pass + elif 't_E' not in as_dict: + pass + elif 'u_0' not in as_dict: + pass + else: + raise KeyError('') + # parallax treated separately because doesnt depend on parametrization raise NotImplementedError('NOT FINISHED YET') From 32d2c221109c925c3445386ec1847f7d5c4dee32 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Sun, 1 Apr 2018 21:38:40 -0400 Subject: [PATCH 06/17] chi2_gradient() mostly finished --- source/MulensModel/event.py | 76 ++++++++++++++++++++++++++++++------- 1 file changed, 63 insertions(+), 13 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 29bfde65..0b275e43 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -300,11 +300,11 @@ def chi2_gradient(self, parameters, fit_blending=None): # for (i, dataset) in enumerate(self.datasets): (data, err_data) = dataset.data_and_err_in_input_fmt() - factor_1 = data - self.fit.get_input_format(data=dataset) - factor_1 /= err_data**2 + factor = data - self.fit.get_input_format(data=dataset) + factor *= -2. / err_data**2 if dataset.input_fmt == 'mag': - factor_1 *= -2.5 / (log(10.) * Utils.get_flux_from_mag(mag)) - factor_1 *= self.fit.flux_of_sources(dataset)[0] + factor *= -2.5 / (log(10.) * Utils.get_flux_from_mag(mag)) + factor *= self.fit.flux_of_sources(dataset)[0] trajectory = Trajectory(dataset.time, self.model.parameters, self.model.get_parallax, self.coords, @@ -312,22 +312,72 @@ def chi2_gradient(self, parameters, fit_blending=None): 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)) - d_x_d_u = trajectory.x / u - d_y_d_u = trajectory.y / u + factor *= d_A_d_u + + factor_d_x_d_u = factor * trajectory.x / u + sum_d_x_d_u = np.sum(factor_d_x_d_u) + factor_d_y_d_u = factor * trajectory.y / u + sum_d_y_d_u = np.sum(factor_d_y_d_u) + dt = dataset.time - 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: - pass + if 't_0' in parameters: + gradient['t_0'] += -sum_d_x_d_u / as_dict['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 / as_dict['t_E']**2) elif 't_E' not in as_dict: - pass + if 't_0' in parameters: + gradient['t_0'] += ( + -sum_d_x_d_u * as_dict['u_0'] / as_dict['t_eff']) + if 'u_0' in parameters: + gradient['u_0'] += sum_d_y_d_u + np.sum( + factor_d_x_d_u * dt / as_dict['t_eff']) + if 't_eff' in parameters: + gradient['t_eff'] += np.sum(factor_d_x_d_u * -dt * + as_dict['u_0'] / as_dict['t_eff']**2) elif 'u_0' not in as_dict: - pass + if 't_0' in parameters: + gradient['t_0'] += -sum_d_x_d_u / as_dict['t_E'] + if 't_E' in parameters: + gradient['t_E'] += (np.sum(factor_d_x_d_u * dt) - + sum_d_y_d_u * as_dict['t_eff']) / as_dict['t_E']**2 + if 't_eff' in parameters: + gradient['t_eff'] += sum_d_y_d_u / as_dict['t_E'] else: - raise KeyError('') - -# parallax treated separately because doesnt depend on parametrization - raise NotImplementedError('NOT FINISHED YET') + 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, + dataset.satellite_skycoord) + dx = trajectory.x - trajectory_no_piE.x + dy = trajectory.y - trajectory_no_piE.y + pi_E_2 = as_dict['pi_E_N']**2 + as_dict['pi_E_E']**2 + E_normalized = as_dict['pi_E_E'] / pi_E_2 + N_normalized = as_dict['pi_E_N'] / pi_E_2 + delta_E = dx * E_normalized + dy * N_normalized + delta_N = dx * N_normalized - dy * E_normalized + 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) + 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) + + 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): """ From 1886deb122b47d400e1f9226e0f4b49a50cdecfd Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Mon, 2 Apr 2018 01:56:03 -0400 Subject: [PATCH 07/17] additional test of Event & MulensData --- source/MulensModel/event.py | 2 +- source/MulensModel/tests/test_Event.py | 27 ++++++++++++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 0b275e43..6f315883 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -303,7 +303,7 @@ def chi2_gradient(self, parameters, fit_blending=None): 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(mag)) + factor *= -2.5 / (log(10.) * Utils.get_flux_from_mag(data)) factor *= self.fit.flux_of_sources(dataset)[0] trajectory = Trajectory(dataset.time, self.model.parameters, diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index a93b2e43..1e177b7d 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -72,6 +72,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 From fff855818b536b28269fd7903466dab9eae3d68b Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Wed, 4 Apr 2018 15:24:00 -0400 Subject: [PATCH 08/17] adding dataset.good in chi2_gradient() --- source/MulensModel/event.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 6f315883..f7c9af34 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -297,7 +297,6 @@ def chi2_gradient(self, parameters, fit_blending=None): else: self.fit.fit_fluxes() - # for (i, dataset) in enumerate(self.datasets): (data, err_data) = dataset.data_and_err_in_input_fmt() factor = data - self.fit.get_input_format(data=dataset) @@ -314,11 +313,11 @@ def chi2_gradient(self, parameters, fit_blending=None): 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 + 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 + 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 - as_dict['t_0'] + 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. @@ -359,20 +358,21 @@ def chi2_gradient(self, parameters, fit_blending=None): trajectory_no_piE = Trajectory(dataset.time, self.model.parameters, parallax, self.coords, dataset.satellite_skycoord) - dx = trajectory.x - trajectory_no_piE.x - dy = trajectory.y - trajectory_no_piE.y - pi_E_2 = as_dict['pi_E_N']**2 + as_dict['pi_E_E']**2 - E_normalized = as_dict['pi_E_E'] / pi_E_2 - N_normalized = as_dict['pi_E_N'] / pi_E_2 - delta_E = dx * E_normalized + dy * N_normalized - delta_N = dx * N_normalized - dy * E_normalized + 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: From 60877d13b0f46b75b133aab272a44b2d1a3f8a1b Mon Sep 17 00:00:00 2001 From: Jennifer Yee Date: Wed, 11 Apr 2018 16:39:06 -0400 Subject: [PATCH 09/17] Created unit test for chi2_gradient(). Added import statement for Trajectory. --- source/MulensModel/event.py | 1 + source/MulensModel/tests/test_Event.py | 14 ++++++++++++++ 2 files changed, 15 insertions(+) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index f7c9af34..ee3e6379 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -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): diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index 1e177b7d..1d32b710 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""" @@ -183,3 +185,15 @@ 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(): + parameters = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87} + gradient = {'t_0': -0.0348918609, 'u_0': -0.0933733096, 't_E': 1.52778388} + + event = Event( + datasets=MulensData(file_name=SAMPLE_FILE_02), model=Model(parameters)) + result = event.chi2_gradient( + [key for key in parameters.keys()], fit_blending=False) + + for key in gradient.keys(): + np.testing.assert_almost_equal(gradient[key], result[key]) From 2e7c4a1be23e57ca5c0101a19a5a522a1d1fcdc5 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Thu, 12 Apr 2018 07:33:11 -0400 Subject: [PATCH 10/17] correcting test_event_chi2_gradient() --- source/MulensModel/tests/test_Event.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index 1d32b710..1be72694 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -188,12 +188,12 @@ def test_event_init_2(self): def test_event_chi2_gradient(): parameters = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87} + params = ['t_0', 'u_0', 't_E'] gradient = {'t_0': -0.0348918609, 'u_0': -0.0933733096, 't_E': 1.52778388} event = Event( datasets=MulensData(file_name=SAMPLE_FILE_02), model=Model(parameters)) - result = event.chi2_gradient( - [key for key in parameters.keys()], fit_blending=False) - - for key in gradient.keys(): - np.testing.assert_almost_equal(gradient[key], result[key]) + result = event.chi2_gradient(params, fit_blending=False) + + np.testing.assert_almost_equal([gradient[key] for key in params], result) + From ddb446c9e9a007caadb41e231b004b48bb7b87fc Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Thu, 12 Apr 2018 07:33:46 -0400 Subject: [PATCH 11/17] a few bugs in Event.chi2_gradient() --- source/MulensModel/event.py | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index ee3e6379..90aede5a 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -306,48 +306,52 @@ def chi2_gradient(self, parameters, fit_blending=None): factor *= -2.5 / (log(10.) * Utils.get_flux_from_mag(data)) factor *= self.fit.flux_of_sources(dataset)[0] + 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, - dataset.satellite_skycoord) + self.model.get_parallax(), self.coords, **kwargs) u_2 = trajectory.x**2 + trajectory.y**2 - u = np.sqrt(u_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] + 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] + 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 / as_dict['t_E'] + 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 / as_dict['t_E']**2) + 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'] / as_dict['t_eff']) + 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 / as_dict['t_eff']) + 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'] / as_dict['t_eff']**2) + 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 / as_dict['t_E'] + 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 * as_dict['t_eff']) / as_dict['t_E']**2 + sum_d_y_d_u * t_eff) / t_E**2 if 't_eff' in parameters: - gradient['t_eff'] += sum_d_y_d_u / as_dict['t_E'] + 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) From eabc5e7bc9d49176aa30d36b37f6a9dcd5b9e93a Mon Sep 17 00:00:00 2001 From: Jennifer Yee Date: Wed, 18 Apr 2018 15:51:08 -0400 Subject: [PATCH 12/17] Updated unit test for chi2_gradient() --- source/MulensModel/event.py | 29 ++++++++++++++++++++------ source/MulensModel/tests/test_Event.py | 8 ++++--- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 90aede5a..1196a247 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -299,12 +299,29 @@ def chi2_gradient(self, parameters, fit_blending=None): self.fit.fit_fluxes() for (i, dataset) in enumerate(self.datasets): - (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] + ## 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: diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index 1be72694..40ddcdfc 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -187,13 +187,15 @@ def test_event_init_2(self): ev = Event(datasets='some_string') def test_event_chi2_gradient(): + # fs = 11.0415734, fb = 0.0 parameters = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87} params = ['t_0', 'u_0', 't_E'] - gradient = {'t_0': -0.0348918609, 'u_0': -0.0933733096, 't_E': 1.52778388} + gradient = {'t_0': 236.206598, 'u_0': 101940.249, + 't_E': -1006.88678} + data = MulensData(file_name=SAMPLE_FILE_02) event = Event( - datasets=MulensData(file_name=SAMPLE_FILE_02), model=Model(parameters)) + datasets=[data], model=Model(parameters)) result = event.chi2_gradient(params, fit_blending=False) np.testing.assert_almost_equal([gradient[key] for key in params], result) - From 7fbeaf192369b3109100e2b05f707b6259ab2907 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Wed, 18 Apr 2018 17:00:51 -0400 Subject: [PATCH 13/17] going back to magnitudes in chi2 gradient --- source/MulensModel/event.py | 24 ++++++++++++------------ source/MulensModel/tests/test_Event.py | 4 +++- 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index 1196a247..a81fd4dd 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -300,12 +300,12 @@ def chi2_gradient(self, parameters, fit_blending=None): 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] + (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() @@ -316,12 +316,12 @@ def chi2_gradient(self, parameters, fit_blending=None): #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) + #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: diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index 40ddcdfc..7d5669b3 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -198,4 +198,6 @@ def test_event_chi2_gradient(): datasets=[data], model=Model(parameters)) result = event.chi2_gradient(params, fit_blending=False) - np.testing.assert_almost_equal([gradient[key] for key in params], result) + reference = np.array([gradient[key] for key in params]) + np.testing.assert_almost_equal(reference/result, 1., decimal=1) + From 28b0688c55de7c714c8017c16ea64702d55efd6c Mon Sep 17 00:00:00 2001 From: Jennifer Yee Date: Wed, 25 Apr 2018 16:42:42 -0400 Subject: [PATCH 14/17] Added parallax and fluxes to gradient test. --- examples/use_cases/use_case_25_fluxes.py | 57 ++++++++++++++++++++++++ source/MulensModel/tests/test_Event.py | 29 ++++++++---- 2 files changed, 78 insertions(+), 8 deletions(-) create mode 100644 examples/use_cases/use_case_25_fluxes.py diff --git a/examples/use_cases/use_case_25_fluxes.py b/examples/use_cases/use_case_25_fluxes.py new file mode 100644 index 00000000..e3b47de5 --- /dev/null +++ b/examples/use_cases/use_case_25_fluxes.py @@ -0,0 +1,57 @@ +""" +Use cases showing how to implement a flux constraint. +""" +import numpy as np +import glob + +import MulensModel + +def chi2_fun(theta, event, parameters_to_fit): + """ + for given event set attributes from parameters_to_fit (list of + str) to values from theta list + """ + for (key, val) in enumerate(parameters_to_fit): + setattr(event.model.parameters, val, theta[key]) + + chi2 = event.get_chi2() + constraint = get_color_constraint(event) + + return chi2 + constraint + +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=1 + + # Color constraint for OB161195 (I_KMT - L_Spitzer) + (source_color, sigma) = (0.78 0.03) + + f_s_ogle = event.get_source_flux(dataset=KMT) + f_s_spitz = event.get_source_flux(dataset=Spitzer) + + color = -2.5 * np.log10(f_s_ogle / f_s_spitz) + + return (color - source_color)**2 / sigma**2 + +# Add data from OB161195 +datasets = [] +filenames = [] # Add list of file names here +for file in filenames: + datasets.append(MulensModel.MulensData( + file_name=os.path.join( + MulensModel.MODULE_PATH, "data", "OB161195", file))) + +# Close-- model +model = MulensModel.Model( + {'t_0': 2457568.7692, 'u_0': -0.05321, 't_E': 9.96, 'rho': 0.00290, + 'pi_E_N': -0.2154, 'pi_E_E': -0.380, + 'alpha': np.rad2deg(-0.9684), 's': 0.9842, 'q': 0.0000543}) + +event = MulensModel.Event(datasets=datasets, model=model) +print(chi2_fun(event)) diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index 7d5669b3..52fc2f1e 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -188,16 +188,29 @@ def test_event_init_2(self): def test_event_chi2_gradient(): # fs = 11.0415734, fb = 0.0 - parameters = {'t_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87} - params = ['t_0', 'u_0', 't_E'] - gradient = {'t_0': 236.206598, 'u_0': 101940.249, + 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} + 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) data = MulensData(file_name=SAMPLE_FILE_02) - event = Event( - datasets=[data], model=Model(parameters)) - 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) + coords = '17:47:12.25 −21:22:58.7' + for test in [test_1, test_2]: + (parameters, params, gradient) = test + event = Event( + datasets=[data], model=Model(parameters), coords=coords) + 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) From febd48d624b84d6f9369d6339047605587b0b002 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Thu, 26 Apr 2018 14:22:56 -0400 Subject: [PATCH 15/17] bug in Event.chi2_gradient() --- source/MulensModel/event.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/MulensModel/event.py b/source/MulensModel/event.py index a81fd4dd..1aa89a11 100644 --- a/source/MulensModel/event.py +++ b/source/MulensModel/event.py @@ -379,7 +379,7 @@ def chi2_gradient(self, parameters, fit_blending=None): 'topocentric': False} trajectory_no_piE = Trajectory(dataset.time, self.model.parameters, parallax, self.coords, - dataset.satellite_skycoord) + **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'] From 4e5b22cd2852d60b696466c87edeab24409a0943 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Thu, 26 Apr 2018 14:23:45 -0400 Subject: [PATCH 16/17] commenting parallax and flux unit test --- source/MulensModel/tests/test_Event.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index 52fc2f1e..a8dff0bc 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -201,16 +201,15 @@ def test_event_chi2_gradient(): '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) - - coords = '17:47:12.25 −21:22:58.7' - for test in [test_1, test_2]: + 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( - datasets=[data], model=Model(parameters), coords=coords) + 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) - From 5d7d8b2f2f0b572de7c94c2f8cc9d1b6cc187946 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Wed, 2 May 2018 15:28:14 -0400 Subject: [PATCH 17/17] version update --- README.md | 2 +- docs/source/conf.py | 4 ++-- docs/source/install.rst | 4 ++-- source/MulensModel/tests/test_Event.py | 6 ++++-- source/MulensModel/version.py | 2 +- 5 files changed, 10 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index af3bb4c9..8f118250 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/source/MulensModel/tests/test_Event.py b/source/MulensModel/tests/test_Event.py index a8dff0bc..78e02071 100644 --- a/source/MulensModel/tests/test_Event.py +++ b/source/MulensModel/tests/test_Event.py @@ -195,14 +195,15 @@ def test_event_chi2_gradient(): 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} + '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'} @@ -213,3 +214,4 @@ def test_event_chi2_gradient(): reference = np.array([gradient[key] for key in params]) np.testing.assert_almost_equal(reference/result, 1., decimal=1) + diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py index 57b9f058..bdc2ab81 100644 --- a/source/MulensModel/version.py +++ b/source/MulensModel/version.py @@ -1,2 +1,2 @@ -__version__ = "1.1.0" +__version__ = "1.2.0"