From 4d2ebee7e4cf3d7cb47ab3f1a0f0d3b2e5863443 Mon Sep 17 00:00:00 2001 From: jenniferyee Date: Tue, 22 Feb 2022 19:15:11 -0500 Subject: [PATCH 1/7] Simplified parallax calculations, now stores delta NE. --- source/MulensModel/fitdata.py | 3 ++ source/MulensModel/trajectory.py | 65 +++++++++++++++++++++----------- 2 files changed, 46 insertions(+), 22 deletions(-) diff --git a/source/MulensModel/fitdata.py b/source/MulensModel/fitdata.py index 814a2a72..1da96278 100644 --- a/source/MulensModel/fitdata.py +++ b/source/MulensModel/fitdata.py @@ -686,6 +686,8 @@ def _get_d_u_d_params(self, parameters): parameters_no_piE.pop('pi_E_N') parameters_no_piE.pop('pi_E_E') + # Isn't this effectively delta tau, delta_beta? --> store that + # information in Trajectory rather than reconstructing it? trajectory_no_piE = Trajectory( self.dataset.time, ModelParameters(parameters_no_piE), **kwargs) @@ -693,6 +695,7 @@ def _get_d_u_d_params(self, parameters): dy = trajectory.y - trajectory_no_piE.y 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 gradient['pi_E_N'] = (d_u_d_x * delta_N + d_u_d_y * delta_E) / det gradient['pi_E_E'] = (d_u_d_x * delta_E - d_u_d_y * delta_N) / det diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index 7a0e8d39..2cf7979a 100644 --- a/source/MulensModel/trajectory.py +++ b/source/MulensModel/trajectory.py @@ -150,10 +150,12 @@ def get_xy(self): # If parallax is non-zero, apply parallax effects: if self.parameters.pi_E is not None: + # Warnings & Checks if self.coords is None: raise ValueError("You're trying to calculate trajectory in " + "a parallax model, but event sky " + "coordinates were not provided.") + keys = ['earth_orbital', 'satellite', 'topocentric'] if set([self.parallax[k] for k in keys]) == set([False]): raise ValueError( @@ -161,27 +163,12 @@ def get_xy(self): 'of parallax dict has to be True ' + '(earth_orbital, satellite, or topocentric)') - # Apply Earth Orbital parallax effect - if self.parallax['earth_orbital']: - [delta_tau, delta_u] = self._annual_parallax_trajectory() - vector_tau += delta_tau - vector_u += delta_u - - # Apply satellite parallax effect - if (self.parallax['satellite'] and - self.satellite_skycoord is not None): - [delta_tau, delta_u] = self._satellite_parallax_trajectory() - vector_tau += delta_tau - vector_u += delta_u - - # Apply topocentric parallax effect - if self.parallax['topocentric'] and self._earth_coords is not None: - # When you implement it, make sure the behavior - # depends on the access to the observatory location - # information as the satellite parallax depends on the - # access to satellite_skycoord. - raise NotImplementedError( - "The topocentric parallax effect not implemented yet") + # Calculate parallax offsets + + self._calculate_delta_NE() + [delta_tau, delta_u] = self._project_delta() + vector_tau += delta_tau + vector_u += delta_u # If 2 lenses, rotate trajectory relative to binary lens axis if self.parameters.n_lenses == 1: @@ -208,10 +195,44 @@ def get_xy(self): self._x = vector_x self._y = vector_y - def _project_delta(self, delta): + def _calculate_delta_NE(self): + self.delta_NE = {} + + if self.parallax['earth_orbital']: + delta_eo = self._get_delta_annual() + for direction in ['N', 'E']: + if direction in self.delta_NE.keys(): + self.delta_NE[direction] += delta_eo[direction] + else: + self.delta_NE[direction] = np.copy(delta_eo[direction]) + + # Apply satellite parallax effect + if (self.parallax['satellite'] and + self.satellite_skycoord is not None): + delta_sat = self._get_delta_satellite() + for direction in ['N', 'E']: + if direction in self.delta_NE.keys(): + self.delta_NE[direction] += delta_sat[direction] + else: + self.delta_NE[direction] = np.copy(delta_sat[direction]) + + # Apply topocentric parallax effect + if self.parallax['topocentric'] and self._earth_coords is not None: + # When you implement it, make sure the behavior + # depends on the access to the observatory location + # information as the satellite parallax depends on the + # access to satellite_skycoord. + raise NotImplementedError( + "The topocentric parallax effect not implemented yet") + + + def _project_delta(self, delta=None): """ Project N and E parallax offset vector onto the tau, beta plane. """ + if delta is None: + delta = self.delta_NE + delta_tau = (delta['N'] * self.parameters.pi_E_N + delta['E'] * self.parameters.pi_E_E) delta_beta = (-delta['N'] * self.parameters.pi_E_E + From 2f68d3851352dac31942b4702ac7b4f84c3465f6 Mon Sep 17 00:00:00 2001 From: jenniferyee Date: Thu, 24 Feb 2022 08:41:12 -0500 Subject: [PATCH 2/7] Made delta_NE a property of Trajectory. Cleaned up some commented code. --- source/MulensModel/fitdata.py | 38 +++++++++----------------------- source/MulensModel/trajectory.py | 23 ++++++++++++++----- 2 files changed, 27 insertions(+), 34 deletions(-) diff --git a/source/MulensModel/fitdata.py b/source/MulensModel/fitdata.py index 1da96278..b5d0d4a0 100644 --- a/source/MulensModel/fitdata.py +++ b/source/MulensModel/fitdata.py @@ -671,34 +671,16 @@ def _get_d_u_d_params(self, parameters): # Below we deal with parallax only. if 'pi_E_N' in parameters or 'pi_E_E' in parameters: - warnings.warn( - "\n\nTests indicate that chi2 gradient for models with " - "parallax has BUGS!!!\n It's better not to use it or contact " - "code authors.\n") - # JCY Not happy about this as it requires importing from other - # modules. It is inelegant, which in my experience often means it - # needs to be refactored. - kwargs = dict() - if self.dataset.ephemerides_file is not None: - kwargs['satellite_skycoord'] = self.dataset.satellite_skycoord - - parameters_no_piE = {**self.model.parameters.as_dict()} - parameters_no_piE.pop('pi_E_N') - parameters_no_piE.pop('pi_E_E') - - # Isn't this effectively delta tau, delta_beta? --> store that - # information in Trajectory rather than reconstructing it? - trajectory_no_piE = Trajectory( - self.dataset.time, ModelParameters(parameters_no_piE), - **kwargs) - dx = trajectory.x - trajectory_no_piE.x - dy = trajectory.y - trajectory_no_piE.y - 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 - gradient['pi_E_N'] = (d_u_d_x * delta_N + d_u_d_y * delta_E) / det - gradient['pi_E_E'] = (d_u_d_x * delta_E - d_u_d_y * delta_N) / det + # warnings.warn( + # "\n\nTests indicate that chi2 gradient for models with " + # "parallax has BUGS!!!\n It's better not to use it or contact " + # "code authors.\n") + + delta_N = trajectory.delta_NE['N'] + delta_E = trajectory.delta_NE['E'] + + gradient['pi_E_N'] = (d_u_d_x * delta_N + d_u_d_y * delta_E) + gradient['pi_E_E'] = (d_u_d_x * delta_E - d_u_d_y * delta_N) return gradient diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index 2cf7979a..12332381 100644 --- a/source/MulensModel/trajectory.py +++ b/source/MulensModel/trajectory.py @@ -134,6 +134,17 @@ def y(self): """ return self._y + @property + def delta_NE(self): + """ + *dict* + + Net North (key='N') and East (key='E') components of the parallax + offset calculated for each time stamp (so sum of the offsets from all + parallax types). + """ + return self._delta_NE + def get_xy(self): """ For a given set of parameters @@ -196,25 +207,25 @@ def get_xy(self): self._y = vector_y def _calculate_delta_NE(self): - self.delta_NE = {} + self._delta_NE = {} if self.parallax['earth_orbital']: delta_eo = self._get_delta_annual() for direction in ['N', 'E']: if direction in self.delta_NE.keys(): - self.delta_NE[direction] += delta_eo[direction] + self._delta_NE[direction] += delta_eo[direction] else: - self.delta_NE[direction] = np.copy(delta_eo[direction]) + self._delta_NE[direction] = np.copy(delta_eo[direction]) # Apply satellite parallax effect if (self.parallax['satellite'] and self.satellite_skycoord is not None): delta_sat = self._get_delta_satellite() for direction in ['N', 'E']: - if direction in self.delta_NE.keys(): - self.delta_NE[direction] += delta_sat[direction] + if direction in self._delta_NE.keys(): + self._delta_NE[direction] += delta_sat[direction] else: - self.delta_NE[direction] = np.copy(delta_sat[direction]) + self._delta_NE[direction] = np.copy(delta_sat[direction]) # Apply topocentric parallax effect if self.parallax['topocentric'] and self._earth_coords is not None: From 97897551cee2f6b277c7e4b337eb85f78344ec15 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Mon, 7 Mar 2022 19:54:46 +0100 Subject: [PATCH 3/7] PEP8 and cleanup --- source/MulensModel/fitdata.py | 9 ++------- source/MulensModel/trajectory.py | 17 ++++++----------- 2 files changed, 8 insertions(+), 18 deletions(-) diff --git a/source/MulensModel/fitdata.py b/source/MulensModel/fitdata.py index b5d0d4a0..bc877d20 100644 --- a/source/MulensModel/fitdata.py +++ b/source/MulensModel/fitdata.py @@ -671,16 +671,11 @@ def _get_d_u_d_params(self, parameters): # Below we deal with parallax only. if 'pi_E_N' in parameters or 'pi_E_E' in parameters: - # warnings.warn( - # "\n\nTests indicate that chi2 gradient for models with " - # "parallax has BUGS!!!\n It's better not to use it or contact " - # "code authors.\n") - delta_N = trajectory.delta_NE['N'] delta_E = trajectory.delta_NE['E'] - gradient['pi_E_N'] = (d_u_d_x * delta_N + d_u_d_y * delta_E) - gradient['pi_E_E'] = (d_u_d_x * delta_E - d_u_d_y * delta_N) + gradient['pi_E_N'] = d_u_d_x * delta_N + d_u_d_y * delta_E + gradient['pi_E_E'] = d_u_d_x * delta_E - d_u_d_y * delta_N return gradient diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index 12332381..90161e1e 100644 --- a/source/MulensModel/trajectory.py +++ b/source/MulensModel/trajectory.py @@ -161,7 +161,6 @@ def get_xy(self): # If parallax is non-zero, apply parallax effects: if self.parameters.pi_E is not None: - # Warnings & Checks if self.coords is None: raise ValueError("You're trying to calculate trajectory in " + "a parallax model, but event sky " + @@ -174,8 +173,6 @@ def get_xy(self): 'of parallax dict has to be True ' + '(earth_orbital, satellite, or topocentric)') - # Calculate parallax offsets - self._calculate_delta_NE() [delta_tau, delta_u] = self._project_delta() vector_tau += delta_tau @@ -202,11 +199,13 @@ def get_xy(self): raise NotImplementedError( "trajectory for more than 2 lenses not handled yet") - # Store trajectory self._x = vector_x self._y = vector_y def _calculate_delta_NE(self): + """ + XXX + """ self._delta_NE = {} if self.parallax['earth_orbital']: @@ -217,7 +216,6 @@ def _calculate_delta_NE(self): else: self._delta_NE[direction] = np.copy(delta_eo[direction]) - # Apply satellite parallax effect if (self.parallax['satellite'] and self.satellite_skycoord is not None): delta_sat = self._get_delta_satellite() @@ -227,16 +225,13 @@ def _calculate_delta_NE(self): else: self._delta_NE[direction] = np.copy(delta_sat[direction]) - # Apply topocentric parallax effect if self.parallax['topocentric'] and self._earth_coords is not None: - # When you implement it, make sure the behavior - # depends on the access to the observatory location - # information as the satellite parallax depends on the - # access to satellite_skycoord. + # When you implement it, make sure the behavior depends on the + # access to the observatory location information as the satellite + # parallax depends on the access to satellite_skycoord. raise NotImplementedError( "The topocentric parallax effect not implemented yet") - def _project_delta(self, delta=None): """ Project N and E parallax offset vector onto the tau, beta plane. From 5da1d2076c2308f5bdb9dc886b5a2f9ada17df4f Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Mon, 7 Mar 2022 20:04:20 +0100 Subject: [PATCH 4/7] Trajectory._calculate_delta_NE(): docstrings, simplifying it, and NE -> N_E --- source/MulensModel/trajectory.py | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index 90161e1e..75a017f3 100644 --- a/source/MulensModel/trajectory.py +++ b/source/MulensModel/trajectory.py @@ -143,7 +143,7 @@ def delta_NE(self): offset calculated for each time stamp (so sum of the offsets from all parallax types). """ - return self._delta_NE + return self._delta_N_E def get_xy(self): """ @@ -173,7 +173,7 @@ def get_xy(self): 'of parallax dict has to be True ' + '(earth_orbital, satellite, or topocentric)') - self._calculate_delta_NE() + self._calculate_delta_N_E() [delta_tau, delta_u] = self._project_delta() vector_tau += delta_tau vector_u += delta_u @@ -202,28 +202,22 @@ def get_xy(self): self._x = vector_x self._y = vector_y - def _calculate_delta_NE(self): + def _calculate_delta_N_E(self): """ - XXX + Calculate shifts caused by microlensing parallax effect. """ - self._delta_NE = {} + self._delta_N_E = {'N': 0., 'E': 0.} if self.parallax['earth_orbital']: - delta_eo = self._get_delta_annual() - for direction in ['N', 'E']: - if direction in self.delta_NE.keys(): - self._delta_NE[direction] += delta_eo[direction] - else: - self._delta_NE[direction] = np.copy(delta_eo[direction]) + delta_annual = self._get_delta_annual() + self._delta_N_E['N'] += delta_annual['N'] + self._delta_N_E['E'] += delta_annual['E'] if (self.parallax['satellite'] and self.satellite_skycoord is not None): - delta_sat = self._get_delta_satellite() - for direction in ['N', 'E']: - if direction in self._delta_NE.keys(): - self._delta_NE[direction] += delta_sat[direction] - else: - self._delta_NE[direction] = np.copy(delta_sat[direction]) + delta_satellite = self._get_delta_satellite() + self._delta_N_E['N'] += delta_satellite['N'] + self._delta_N_E['E'] += delta_satellite['E'] if self.parallax['topocentric'] and self._earth_coords is not None: # When you implement it, make sure the behavior depends on the From 8593e1dfb047084f3b0caceff79dd53bc8e77e33 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Mon, 7 Mar 2022 20:05:28 +0100 Subject: [PATCH 5/7] updating version number --- README.md | 4 ++-- source/MulensModel/version.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index bf1505ae..40094d01 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ [**Detailed documentation: https://rpoleski.github.io/MulensModel/**](https://rpoleski.github.io/MulensModel/) -[Latest release: 2.4.5](https://github.com/rpoleski/MulensModel/releases/latest) and we're working on further developing the code. +[Latest release: 2.5.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 can then be input into an arbitrary likelihood function to find the best-fit parameters. @@ -52,5 +52,5 @@ If you want to contribute to MulensModel, then please see [this file](CONTRIBUTI [![GitHub stars](https://badgen.net/github/stars/rpoleski/MulensModel)](https://GitHub.com/rpoleski/MulensModel/stargazers/) [![MIT license](https://img.shields.io/badge/License-MIT-blue.svg)](https://lbesson.mit-license.org/) -file revised Feb 2022 +file revised Mar 2022 diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py index c9e914fc..50062f87 100644 --- a/source/MulensModel/version.py +++ b/source/MulensModel/version.py @@ -1 +1 @@ -__version__ = "2.4.6" +__version__ = "2.5.0" From 6fd79c934d7078bfc9aaedee5f08b76e94b7d5ff Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Mon, 7 Mar 2022 23:40:23 +0100 Subject: [PATCH 6/7] Trajectory.parallax_delta_N_E instead of .delta_NE --- source/MulensModel/trajectory.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index 75a017f3..8e6bf873 100644 --- a/source/MulensModel/trajectory.py +++ b/source/MulensModel/trajectory.py @@ -135,7 +135,7 @@ def y(self): return self._y @property - def delta_NE(self): + def parallax_delta_N_E(self): """ *dict* @@ -231,7 +231,7 @@ def _project_delta(self, delta=None): Project N and E parallax offset vector onto the tau, beta plane. """ if delta is None: - delta = self.delta_NE + delta = self.parallax_delta_N_E delta_tau = (delta['N'] * self.parameters.pi_E_N + delta['E'] * self.parameters.pi_E_E) From 752e7774f2f838d7484ab617d9ec794f9e0e75a4 Mon Sep 17 00:00:00 2001 From: Radek Poleski Date: Tue, 8 Mar 2022 19:01:20 +0100 Subject: [PATCH 7/7] updating version to be consistent with master --- README.md | 2 +- source/MulensModel/version.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 40094d01..b101ed6e 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ [**Detailed documentation: https://rpoleski.github.io/MulensModel/**](https://rpoleski.github.io/MulensModel/) -[Latest release: 2.5.0](https://github.com/rpoleski/MulensModel/releases/latest) and we're working on further developing the code. +[Latest release: 2.7.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 can then be input into an arbitrary likelihood function to find the best-fit parameters. diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py index 50062f87..2614ce9d 100644 --- a/source/MulensModel/version.py +++ b/source/MulensModel/version.py @@ -1 +1 @@ -__version__ = "2.5.0" +__version__ = "2.7.0"