rpoleski / MulensModel

Microlensing Modelling package
https://rpoleski.github.io/MulensModel/
Other
57 stars 15 forks source link

Version 3 - main thread #116

Closed rpoleski closed 1 day ago

rpoleski commented 10 months ago

What is decided to be done for Version 3:

Added later (only high-level ones):

Also check MM_v3.md file for less changes that were not accepted yet.

Not decided yet and important:

jenniferyee commented 10 months ago

For shifting alpha by 180 deg, I think we should make this change for V3 because it's wrong and we should fix it. However, is it possible to set up the fix in a similar way to MAG_ZEROPOINT so the user can set it to the V2 convention if desired.

jenniferyee commented 10 months ago

@rpoleski

Why are there two ways to calculate point source, binary lens magnification? We have both WittMao95 and VBBL. In the new architecture, this would be two separation magnification classes. Except the user never accesses these. Anyway, seems very complicated, so there must be a reason...

rpoleski commented 10 months ago

Long time ago VBBL had a feature/bug: it checked if source-lens distance is larger than 10 and if so, then point lens approximation was used. I don't know if it's still true. At some point we decided to leave both versions since they were already coded.

rapoliveira commented 10 months ago

There is a FutureWarning in the function _set_coords() of MulensData class: "coords will be deprecated in future. There is no reason to tie this to a given dataset"

Should this variable be removed at this point? This is the only deprecated missing in MulensData now.

rpoleski commented 10 months ago

@jenniferyee You have 2 defintions of FitData.get_d_A_d_rho(). Issue was spotted by @rapoliveira

jenniferyee commented 9 months ago

@rpoleski Is there an example for the exception that is raised from VBBL PSBL calculation (or even a unit test)? Just having the code check for an exception makes debugging hard. I introduced some bugs with my refactor, and they weren't caught because the calculation just switched to WM95. Relatedly, if VBBL returns a magnification vector, why do we have the check magnification < 1?

see BinaryLens.point_source_magnification() method.

jenniferyee commented 9 months ago

@rpoleski I don't think a try/except is the right way to handle VBBL PSBL failures.

  1. The except clause is too broad, so we can miss fixable errors. In my tests, VBBL PSBL is failing because it is receiving the wrong number of arguments ("this function takes at least 7 arguments (4 given)"). This is also a problem in the master branch. So either this never worked or something changed with VBBL.

  2. It would be better to model this after the way we deal with problems with point_source method --> point_source_point_lens method, so make it the responsibility of the user to change the calculation if it fails for mysterious reasons. Then, we would have "point_source" and "point_source_VBBL" --> VBBL PSBL calculation and "point_source_WittMao95" --> WM95 calculation. This can be easily added to the current architecture. It would also prevent the complete suppression of errors.

(I also need help fixing the VBBL PSBL error.)

jenniferyee commented 9 months ago

Note on "use_planet_frame": these should just be implemented as child classes of the standard frame.

So there would be BinaryLensPointSourceWM95Magnification() AND BinaryLensPointSourceWM95PlanetFrameMagnification().

These names are getting long, so we might also consider abbreviations, where BLPS_WM95Magnification inherits BinaryLensPointSourceWM95Magnification.

rpoleski commented 9 months ago
  1. The except clause is too broad, so we can miss fixable errors. In my tests, VBBL PSBL is failing because it is
    receiving the wrong number of arguments ("this function takes at least 7 arguments (4 given)"). This is also a problem in the master branch. So either this never worked or something changed with VBBL.

I've checked that. It turned out that the bug was introduced while merging shear & convergence branch. I've made a quick fix in 3e31b57 . Proper solution requires more work on importing VBBL :(

rpoleski commented 9 months ago

First, I think it would be best to fully solve issue with VBBL 4 vs. 7 parameters in master branch, then merge it to Version3_1 branch.

Second, maybe this error was the caused somebody (Keto Zhang?) to report that our approach to importing C++ is too slow and we should try to pass vectors once instead of passing floats in a loop. Based on that comment, I've started vbbl_vectors branch, which is still not finished.

jenniferyee commented 9 months ago

Finished updating documentation. Still need help with fixing binarylens calculations. This is a prerequisite for finishing the binarylenswithshear refactor.

Next task: fixing binary source implementation. On deck: replacing the set_magnification_objects with "factories" to make it easier for the user to add their own magnification calculations. Side quest: create a parent MagnificationObject class of which the children will be PointSourcePointLensMagnification and BinaryLensPointSourceMagnification, see comments in binarylens.py

jenniferyee commented 9 months ago

@rpoleski documentation for BinaryLensPointSourceWM95Magnification has missing information. See boldface.

rpoleski commented 9 months ago

@jenniferyee Do you mean the reference to WM95 paper? If not, then please be specific.

jenniferyee commented 7 months ago

Status Update: Implemented a SourceParameters class and the ModelParams unit tests now pass. However, Example 11 and several other binary source unit tests don't work because it checks for a "ModelParameters" object rather than a "SourceParameters" object. In addition, there is a bunch of code in "SourceParameters" that is copy-pasted directly from ModelParameters. This means these two objects are related. The question is whether they are both children of a "MulensParameters" class or whether "ModelParameters" is a child of "SourceParameters." (I don't think SourceParameters is a child of ModelParameters because SourceParameters is more simple.) Regardless, the checks for "ModelParameters" objects need to be updated.

Relatedly: there is no real reason to store the source information differently for 1 source rather than multiple sources.

rpoleski commented 6 months ago

I was supposed to make notes on classes called factories (wiki link). Here is an example code:

model_1 = mm.GetModel({'t_0': 0, 'u_0': 0.5, 't_E': 10})
print(type(model_1))
# returns ModelPointSourcePointLens

model_2 = mm.GetModel({'t_0': 0, 'u_0': 0.5, 't_E': 10, 's': 1.5, 'q': 0.1, 'alpha': 123.})
print(type(model_2))
# returns ModelPointSourceBinaryLens

Note - a single class GetModel produced objects of different types based on the input it got.

jenniferyee commented 6 months ago

Current thoughts on refactoring ModelParameters to support multiple sources: Because we can add various components to a model independently (e.g. additional sources or lenses, parallax), this is similar to how astropy.modeling allows you to add various models together. Ergo, the ModelParameters refactor should/could follow the basic architecture of astropy.modeling.CompoundModel. (still figuring out how that would work in practice) https://docs.astropy.org/en/stable/api/astropy.modeling.CompoundModel.html#astropy.modeling.CompoundModel

https://docs.astropy.org/en/stable/_modules/astropy/modeling/core.html#CompoundModel

For more on simplifying architectures using inheritance, see also: https://realpython.com/inheritance-composition-python/

jenniferyee commented 6 months ago

Spent some time thinking about class diagrams, inheritance, and the order of adding effects. Although we might think finite source effects, extra lenses, parallax, etc. can be added in any order, the code needs a specific structure for adding these components. I propose:

  1. Basic Model
  2. N x Extra Lenses (Including Orbital Motion)
  3. Finite (Primary Source)
  4. M x Extra Sources (Including Finite Source, then xallarap)
  5. Parallax

Hence, a Binary Lens model with Binary Sources (both exhibiting finite source effects) including parallax would be built as follows:

Basic Model (t0, u0, tE)

A PointSource model with a Binary Lens exhibiting full keplerian orbital motion and parallax would be:

Basic Model (t0, u0, tE)

A triple source, point lens model, with finite source effects only for the third source: Basic Model (t0, u0, tE)

rpoleski commented 6 months ago

I'm still not getting why you want to change ModelParameters. For me the current master branch is quite good. Ok, one obvious change is to keep multiple sources as a list, not hard-coded _1 and _2.

Minor comment:

A triple source, point lens model, with finite source effects only for the third source: Basic Model (t0, u0, tE) ...

I will keep our old decision to call these parameters (t_0_1, u_0_1, t_E) if there is more than one source. The idea was that parameters ending in _N are for specific source and all other (most importantly t_E but also pi_E_E/N) are shared by all sources. This approach seems clear for user and relatively easy to code.

jenniferyee commented 6 months ago

Multisource problems all resolved. The thoughts listed above were too complicated.

jenniferyee commented 6 months ago

merged master into Version 3_1, but there are a lot of problems with incorporating xallarap. The next task will be to resolve the unit test failures. commit 721c4fc5e4587ff96565fa30b144ba3ad58d3e33

jenniferyee commented 6 months ago

Merge complete: 28284c7406307a24fc30784b67e8ef7f2c95ec02

All the unit tests that fail also failed pre-merge.

rpoleski commented 6 months ago

After discussion on May 20, 2024:

rpoleski commented 6 months ago

@jenniferyee, Where do you save plots from examples?

jenniferyee commented 5 months ago

Unit Tests needed in test_FitData.py

I found that satellite_skycoord was not being correctly masked when creating MagnificationCurve objects in FitData. _set_data_magnification_curves(). I fixed this, but it should have a unit test.

Likewise, I noticed that behavior for n_sources == 1 and 2 was explicitly defined, which probably means FitData still needs to be updated to handle arbitrary numbers of sources. Therefore, we should write a unit test for this first.

rpoleski commented 5 months ago

@jenniferyee - please make the tests faster. I've found that 3 tests take half of the total tests time:

I've checked that the first one tests small number of values, but a large number is calculated earlier. Whole trajectories (defined in TestPointSourcePointLensMagnification:setUp()) have magnifications calculated, which takes a lot of time. The second test seems to be identical to the first one.

rpoleski commented 5 months ago

@jenniferyee File binarylens.py has your comment in BinaryLensPointSourceWM95Magnification.get_magnification():

        # JCY: I think this is wrong:
        # The origin of the coordinate system is at the center of mass and
        # both masses are on X axis with higher mass at negative X; this
        # means that the higher mass is at (X, Y)=(-s*q/(1+q), 0) and
        # the lower mass is at (s/(1+q), 0).

What is supposed to be wrong? Did you mean (20 lines later):

        self._position_z1 = -self.separation + 0.j
        self._position_z2 = 0. + 0.j

? I think the above is correct - the change of coordinate system is done sooner and the above lines are consistent with the polynomial coefficients.

rpoleski commented 5 months ago

@jenniferyee , in FiniteSourceUniformGould94Magnification in the method get_d_A_d_u() there is no multiplication by self.b0. When I added it, 5 more tests failed. I think the tests have wrong values. Can you comment on that?

jenniferyee commented 5 months ago

@jenniferyee - please make the tests faster. I've found that 3 tests take half of the total tests time:

  • tests/test_PointLens.py::TestFiniteSourceLDYoo04DirectMagnification::test_magnification
  • tests/test_PointLens.py::TestFiniteSourceLDYoo04DirectMagnification::test_get_magnification
  • tests/test_PointLens.py::TestFiniteSourceLDYoo04DirectMagnification::test_b1

I've checked that the first one tests small number of values, but a large number is calculated earlier. Whole trajectories (defined in TestPointSourcePointLensMagnification:setUp()) have magnifications calculated, which takes a lot of time. The second test seems to be identical to the first one.

Fixed in commit e336a45

jenniferyee commented 5 months ago

@jenniferyee File binarylens.py has your comment in BinaryLensPointSourceWM95Magnification.get_magnification():

        # JCY: I think this is wrong:
        # The origin of the coordinate system is at the center of mass and
        # both masses are on X axis with higher mass at negative X; this
        # means that the higher mass is at (X, Y)=(-s*q/(1+q), 0) and
        # the lower mass is at (s/(1+q), 0).

What is supposed to be wrong? Did you mean (20 lines later):

        self._position_z1 = -self.separation + 0.j
        self._position_z2 = 0. + 0.j

? I think the above is correct - the change of coordinate system is done sooner and the above lines are consistent with the polynomial coefficients.

I added this comment in commit daf1afe

Prior to this commit, the statement

"Calculate point source magnification for given position. The origin of the coordinate system is at the center of mass and both masses are on X axis with higher mass at negative X; this means that the higher mass is at (X, Y)=(-s*q/(1+q), 0) and the lower mass is at (s/(1+q), 0)."

was part of the documentation for BinaryLensPointSourceWM95Magnification.get_magnification(). However, I do not think this statement is correct for WM95.

Anyway, the main action item is that the documentation for BinaryLensPointSourceWM95Magnification.source_x needs the origin of the WM95 calculations to be added. Then, the above comment can be deleted.

rpoleski commented 5 months ago

Do we need .source_x/y to be public in all the binary classes? It's based on the trajectory that is public and convention that is required by WM95, VBBL, or AC methods. Hence, .source_x/y are the implementation detail. I've checked that they're not directly tested. Do you agree to remove public access to these properties and keep them as private?

EDIT: done in 8ee4d2ee

jenniferyee commented 5 months ago

@jenniferyee , in FiniteSourceUniformGould94Magnification in the method get_d_A_d_u() there is no multiplication by self.b0. When I added it, 5 more tests failed. I think the tests have wrong values. Can you comment on that?

yes. This is wrong. But there is a b0 factor in get_d_A_d_params, so maybe the tests fail because there are too many factors of b0. Also, the current get_d_A_d_u() is identical to get_d_A_d_u_PSPL(), so the d_A_d_u tests may be doing the same comparison as for the PSPL case, which is wrong and would fail if you fixed get_d_A_d_u. Anyway,

Next Steps:

  1. rederive d_A_d_u and d_A_d_params
  2. Update tests
  3. Update code.
rpoleski commented 4 months ago

I've added the class for 2L1S magnification that calls VBBL first and switches to WM95 if the former fails: commit 267cc0f1. @jenniferyee you may want to review it.

jenniferyee commented 4 months ago

issues with get_d_A_d_u resolved. @rpoleski you may want to review it. commit d5bbcae874de4809e4675db6d97ffd385faba923

jenniferyee commented 4 months ago

Unit Tests needed in test_FitData.py

  • [ ] test_bad_data_w_ephemerides
  • [ ] test_multiple_sources

I found that satellite_skycoord was not being correctly masked when creating MagnificationCurve objects in FitData. _set_data_magnification_curves(). I fixed this, but it should have a unit test.

Likewise, I noticed that behavior for n_sources == 1 and 2 was explicitly defined, which probably means FitData still needs to be updated to handle arbitrary numbers of sources. Therefore, we should write a unit test for this first.

Other n_source tests that are needed:

Generally:

jenniferyee commented 4 months ago

ation that calls VBBL first and switches to WM95 if the former fails: commit 267cc0f1. @jenniferyee you may want to review it.

Looks fine to me.

rpoleski commented 4 months ago

I'm guessing more than 2 source and xallarap fails currently, so I suggest to forbid that.

jenniferyee commented 4 months ago

@rpoleski for multi-source models, having a single line of parameters in repr is very long + we have the keys ordered by type, but I would prefer them to be ordered by source. see example_25_multiple_sources/ex_25_1L4S_model.py

Should we change repr for multi-source models to

  1. group parameters by source?
  2. output source information on separate lines/in separate blocks?
jenniferyee commented 4 months ago

Note to self: Apparently, N sources was not implemented in model.py.

Next step: write a unit test (or several UTs), e.g., calculate the flux from a model as would be needed for FitData.

Then: implement N sources in model.py and retry ex_25_1L4S_model.py

jenniferyee commented 3 months ago

Implemented unit tests for multi-source models in Model and FitData. They fail (as expected). So, the next step is to implement multiple sources in model.py

6381f246a8f71a443ee5a03f5f3692af78696c3e 70cb82784da1aa6b985a3e2ac9b1e9b318c84696 c84e666973ddc679f0cf6cd855833ef566d8fca7

jenniferyee commented 3 months ago

Implemented multi-source Models. Example 25 1L4S plotting works.

@rpoleski One significant API change was to change default for separate keyword to None. For N sources, the default behavior is now separate = True. See 01dcaf32db48eaec3020edb569ab2471cf26fc64 6a6df2c6064a8bca4cb9ba09866bf18af3d8ccb4

jenniferyee commented 3 months ago

It is now possible to plot models with arbitrary numbers of sources.

@rpoleski I updated the definitions of t_0_par and t_0_kep so that if they are not set by the user, they are initialized to t_0_1 for multi-source models.

See: 4c5c8340deaf5d2c5f7c5c01aad340f13d651126 and prior commits

rpoleski commented 3 months ago

I'm behind with @jenniferyee posts for 2 weeks, so let's try to get me on track:

  1. Demonstrate maximally complex model --> 2LNS w/ parallax and orbital motion

Note that we already have 2L2S model in example 19.

  1. __repr__ for multi-source models - I agree that it's getting very long and can be confusing. We have parameters that affect all sources (e.g., t_E or pi_E_E) and specific ones (e.g., t_0_1). I'm guessing it would be best to print joint parameters first (2 lines) and then specific ones in the following lines (2*N lines). For 3 sources there will be 8 lines total.
  2. t_0_par/kep - OK
  3. I've seen you added code like: self.__getattr__('_source_{0}_parameters'.format(source_num)), which I understand but at the same time I find a bit complicated. Simpler approach: self._source_i_parameters is a list and above is equivalent to self._source_i_parameters[source_num-1] (we can also add None as 0-th element of that list, so that we don't have to type -1 in these indexes. What do you think about that?
rpoleski commented 3 months ago
  1. I've just realized that if t_0_kep is not provided, but t_0_par is provided, then the former should default to the latter. Honestly, I thought we already have that coded.
jenniferyee commented 3 months ago

as of ac630e4856eed7b101fee7bca3a46e4a8f4ad6b9, multiple sources is nearly fully implemented. Remaining action item: separate gammas for the different sources, e.g., self._check_gamma_for_2_sources(gamma) and unit tests.

jenniferyee commented 3 months ago
  1. I've just realized that if t_0_kep is not provided, but t_0_par is provided, then the former should default to the latter. Honestly, I thought we already have that coded.

@rpoleski Is the inverse also true? As in, if t_0_par is not provided, but t_0_kep is provided, then t_0_par = t_0_kep?

jenniferyee commented 3 months ago
  1. I've seen you added code like: self.__getattr__('_source_{0}_parameters'.format(source_num)), which I understand but at the same time I find a bit complicated. Simpler approach: self._source_i_parameters is a list and above is equivalent to self._source_i_parameters[source_num-1] (we can also add None as 0-th element of that list, so that we don't have to type -1 in these indexes. What do you think about that?

@rpoleski The problem with the self._source_i_parameters as a list approach is that source_1_parameters and source_2_parameters were already implemented. It would be strange to change the API just because we now have a source 3. So, the reason for introducing self.__getattr__() is to maintain the existing system for referencing source parameters.

rpoleski commented 3 months ago
  1. Ok, I see there are ._source... and .source... attributes and the __getattr__ is called on the former only 3 times. Changing that one doesn't change API so we could do that easily. Probably not worth effort.

  2. As in, if t_0_par is not provided, but t_0_kep is provided, then t_0_par = t_0_kep?

I think it should be that way.

  1. Minor simplification - can I remove the first 3 lines of ModelParameters.__getattr__()? Doing that doesn't make any test to fail.
  2. What is the purpose of Model._N_source_attr?
rpoleski commented 2 months ago

alpha shift by 180 deg done in 6eef71c. Examples that still need corrections: 1, 8, 13, 16, and 22. I also have to consider uniformcausticsampling.py file.

EDIT - currently remaining: only 8.

rpoleski commented 1 month ago

Relatively small tasks:

rpoleski commented 1 month ago

@jenniferyee Function which_parameters() was just removed: 0e587b2