diff --git a/developers_board.md b/developers_board.md index e230b5e80..6dd5f8049 100644 --- a/developers_board.md +++ b/developers_board.md @@ -1,5 +1,4 @@ # Next big tasks - think about the order: -1. xallarap 2. full Keplerian motion of the lens (because 2-parameter approximation is unphysical and is already implemented) 3. triple sources 4. triple lenses @@ -8,17 +7,12 @@ 7. terrestial parallax # Next smaller tasks: -1. gradient for parallax model issue 2. ModelParameters - note that t_0_1 etc. are for "binary source" models -2. update VBBL 3. make directory with reference plots from examples 4. remove unused branches 5. UC40 - finish it and implement Event.get\_lc and plot\_lc -6. example 16 7. UC38 - finish it -8. chi2 gradient for rho or t_star - first move code from FitData to Trajectory and PointLens (which means new minor version) 9. faster FFP calculations -10. print Model/ModelParameters: add t\_0\_par somewhere - maybe only when it is != t\_0 11. "import MulensModel as mm" everywhere: use_cases/->30,31,34 examples/run_time_tests/check_system.py examples/checks.py TUTORIALS 12. satellite positions read from a file similar to data challenge files 13. check open issues on github @@ -42,7 +36,7 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * more on setup.py: [link](https://github.com/kennethreitz/setup.py) * compile with "pedantic" flags for compilers * Documentation - * magnification\_methods.pdf - add full references there and link this file in the code + * magnification\_methods.pdf - add full references there * Sagan workshop hands-on activity in MM * examples as ipython notebooks * Add \_\_repr\_\_ functions to Lens and Source @@ -61,8 +55,6 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * Fit.fit\_fluxes docstring to be updated * which\_parameters() - note that it doesnt work for binary source parameters, but the parameters work properly; just BSPL and rho\_2 etc. are optional * parallax models - * binary source - there should be one Fit less in Event.get\_chi2xxx functions - if there are 2 sources, then calculate magnification and use F\_S = F\_S\_1+F\_S\_2 and get it from self.model.fit; most probably adding some function to Fit would help - * test binary-lens binary-source * different t\_E for each source (correct Model.set\_times) * test binary source with exactly one rho\_X defined * add t\_eff\_1, t\_eff\_2 @@ -83,7 +75,6 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * get gamma/u LD coeffs from Claret papers etc. * [Lee+09](https://ui.adsabs.harvard.edu/abs/2009ApJ...695..200L/abstract) - gradient calculations for uniform source, also faster calculations - profile * FSPL for large sources using [Agol 2003](https://ui.adsabs.harvard.edu/abs/2003ApJ...594..449A/abstract) - * Xallarap (see below for references) * Quadratic limb darkening * Multi-lens ray shooting: * mapmaking version which adds new rays as needed (but remember that it runs for fixed (s,q) only!) @@ -132,7 +123,6 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * Event should sync information on which of the 3 types of parallax are used, so that if it is specified for event, then there will be exception if one dataset is missing earth\_coords etc. In general there should be some way to make sure which parallax types are used in which calculation of magnification. * Class Event should have not only set\_datasets() methods but also add\_datasets(), i.e. a similar method that appends datasets to self.\_datasets. * reduce calls to Fit.fit\_fluxes() - * add finite source in chi2\_gradient() * chi2\_gradient() should cope NaN values in a way similar to get\_chi2() * **check all functions that should pass** fit\_blending parameter - Event.plot\_model, what else??? Already done: Event.get\_ref\_fluxes() * chi2 with maximum value provided - if the chi2 for point-source gives chi2 larger than specified limit, then finite source calculations are not undertaken (this should significantly speed-up MultiNest) @@ -214,7 +204,7 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * when plotting data, make sure that max/min limits on Y axis include errorbars, if the errorbars are shown * export/save given data file in scale of other dataset and model * this line may be wrong for some values of char: kwargs['fmt'] = kwargs['fmt'].replace(char, "") - * for plt.scatter() the color can be set as 'facecolor', 'facecolors', or 'edgecolors' and this should be dealt with in _set_plot_properties() + * for plt.scatter() the color can be set as 'facecolor', 'facecolors', or 'edgecolors' and this should be dealt with in \_set\_plot\_properties() * for plotting X for bad data use large size and/or thinner line * separate colors (or even kwargs) for X-es as an option (to get contrasting colors see https://itsphbytes.wordpress.com/2016/08/29/complementary-colors-python-code/) * PointLens class: @@ -230,7 +220,6 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * _warning when too many annual parallax calculations are conducted_ * _\_get\_delta\_satellite() should be using self.times_ * annual parallax caching - if moved to MulensData, then would be faster because hashing of coords and time vector takes time - * maybe Trajectory should be able to plot itself, and Model.plot\_trajectory() should call it - it would be easier for binary sources etc. * colorscale time or magnification (see Fig 2 in [Ranc+19](https://arxiv.org/abs/1810.00014)) * plot in units of theta\_star (or even days) instead of theta\_E * UniformCausticSampling class: @@ -263,7 +252,6 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * _Hamiltonian MCMC [link 1](http://arogozhnikov.github.io/2016/12/19/markov_chain_monte_carlo.html) and [link 2](https://theclevermachine.wordpress.com/2012/11/18/mcmc-hamiltonian-monte-carlo-a-k-a-hybrid-monte-carlo/) and [link 3](https://colindcarroll.com/2019/04/11/hamiltonian-monte-carlo-from-scratch/)_ * example with [parallel EMCEE](https://emcee.readthedocs.io/en/stable/tutorials/parallel/) * _plot many models from posterior_ - * **chi2 per dataset** * **scipy.curve\_fit() and print parameter uncertainties** * **corner plots; they require [corner](https://github.com/dfm/corner.py), [pyhdust](https://github.com/danmoser/pyhdust), or [pyGTC](https://pypi.org/project/pyGTC/)** * _F\_s for MOA data for MB08310 differs from Janczak paper - is it caused by FSPL vs. FSBL models?_ @@ -277,8 +265,6 @@ Changes for planned v3 are here: [documents/MM\_v3.md](documents/MM_v3.md) * note in PSPL tutorial about plotting data in MulensData * add example that after calling Event.get\_chi2() use Event.fit to get e.g. magnifications so that the magnification is not calculated twice * **satellite data fitted and plotted - what is missing now?** - * use ast.literal\_eval() for .cfg files to read dict or list, e.g., for MulensData options. - * _high-level fitting example where we dont care how complicated it is, we just want to make it simple and useful for the user_ * some cfg files use "../data/..." - change it to MM.DATA\_PATH somehow * in emcee we should check if all the starting points are in prior * check if MM correctly interacts with scipy.optimize.leastsq and maybe add an example diff --git a/examples/example_16/ob03235_2_full.yaml b/examples/example_16/ob03235_2_full.yaml index 283140f09..b660b92e4 100644 --- a/examples/example_16/ob03235_2_full.yaml +++ b/examples/example_16/ob03235_2_full.yaml @@ -82,13 +82,14 @@ plots: # makes legend in 2 columns loc: lower center second Y scale: - # This adds second Y axis on the right-hand side. Only first line below is required. - magnifications: [2, 3, 4, 5, 6, 7, 8, 9] - # If you don't know what will be the range of magnifications on your plot, then make - # a test run with some very small and very large values and a warning will tell you - # what is exact range on the plot. - labels: [a, b, c, d, e, f, g, h] - # If you remove the line above, then magnification values will be printed. + # This adds second Y axis to the right side. Only magnifications key is required. + magnifications: optimal + # magnifications: [2, 3, 4, 5, 6, 7, 8, 9] + # If you want to provide magnification values but don't know what will be the range + # of magnifications on your plot, then make a test with very small and large numbers + # and a warning will tell you the exact range. + # labels: [a, b, c, d, e, f, g, h] + # The list of labels above can not be given if magnifications = "optimal" label: What is shown on this axis? color: magenta trajectory: diff --git a/examples/example_16/ulens_model_fit.py b/examples/example_16/ulens_model_fit.py index 5ea663177..614f29c2a 100644 --- a/examples/example_16/ulens_model_fit.py +++ b/examples/example_16/ulens_model_fit.py @@ -11,8 +11,8 @@ import numpy as np from scipy.interpolate import interp1d from matplotlib import pyplot as plt -from matplotlib import gridspec, rc, rcParams, rcParamsDefault -from matplotlib.backends.backend_pdf import PdfPages +from matplotlib import gridspec, rcParams, rcParamsDefault +# from matplotlib.backends.backend_pdf import PdfPages import_failed = set() try: @@ -38,7 +38,7 @@ except Exception: raise ImportError('\nYou have to install MulensModel first!\n') -__version__ = '0.33.0' +__version__ = '0.34.2' class UlensModelFit(object): @@ -334,6 +334,8 @@ def _check_MM_version(self): """ Check if MulensModel is new enough """ + # code_version = "{:} and {:}".format(mm.__version__, __version__) + # print('\nMulensModel and script versions:', code_version, end='\n\n') if int(mm.__version__.split('.')[0]) < 2: raise RuntimeError( "ulens_model_fit.py requires MulensModel in version " @@ -717,17 +719,27 @@ def _check_plots_parameters_best_model_Y_scale(self): 'Unknown settings for "second Y scale" in ' '"best model": {:}'.format(unknown)) if not isinstance(settings['magnifications'], list): - raise TypeError( - '"best model" -> "second Y scale" -> "magnifications" has to ' - 'be a list, not ' + str(type(settings['magnifications']))) - for value in settings['magnifications']: - if not isinstance(value, (int, float)): + if settings['magnifications'] != 'optimal': raise TypeError( - 'Wrong value in magnifications: ' + str(value)) + '"best model" -> "second Y scale" -> "magnifications" has ' + 'to be a list or "optimal", not ' + + str(type(settings['magnifications']))) + else: + for value in settings['magnifications']: + if not isinstance(value, (int, float)): + raise TypeError( + 'Wrong value in magnifications: ' + str(value)) if 'labels' not in settings: - settings['labels'] = [ - str(x) for x in settings['magnifications']] + if settings['magnifications'] != 'optimal': + settings['labels'] = [ + str(x) for x in settings['magnifications']] + else: + settings['labels'] = [] else: + if settings['magnifications'] == 'optimal': + raise ValueError( + 'In "best model" -> "second Y scale", labels can not be ' + 'provided if "magnifications" is defined as "optimal"') if not isinstance(settings['labels'], list): raise TypeError( '"best model" -> "second Y scale" -> "labels" has to be ' @@ -1113,11 +1125,11 @@ def _parse_fitting_parameters_EMCEE(self): 'got: ' + name) if path.exists(name): if path.isfile(name): - msg = "Exisiting file " + name + " will be overwritten" + msg = "Existing file " + name + " will be overwritten" warnings.warn(msg) else: raise ValueError("The path provided for posterior (" + - name + ") exsists and is a directory") + name + ") exists and is a directory") self._posterior_file_name = name[:-4] self._posterior_file_fluxes = None @@ -1284,7 +1296,7 @@ def _check_output_files_MultiNest(self): existing.append(file_name) if len(existing) > 0: - message = "\n\n Exisiting files will be overwritten " + message = "\n\n Existing files will be overwritten " message += "(unless you kill this process)!!!\n" warnings.warn(message + str(existing) + "\n") @@ -2309,6 +2321,11 @@ def _parse_results_EMCEE(self): This version works with EMCEE version 2.X and 3.0. """ + if self._yaml_results: + lst = [mm.__version__, __version__] + code_version = "MulensModel and script versions: {:}".format(lst) + print(code_version, **self._yaml_kwargs) + accept_rate = np.mean(self._sampler.acceptance_fraction) out = "Mean acceptance fraction: {0:.3f}".format(accept_rate) print(out) @@ -2323,6 +2340,11 @@ def _parse_results_EMCEE(self): print(out, **self._yaml_kwargs) self._extract_posterior_samples_EMCEE() + if self._yaml_results and isinstance(self._fixed_parameters, dict): + print("Fixed parameters:", **self._yaml_kwargs) + for (key, value) in self._fixed_parameters.items(): + print(" {:} : {:}".format(key, value), **self._yaml_kwargs) + print("Fitted parameters:") self._print_results(self._samples_flat) if self._yaml_results: @@ -2343,6 +2365,10 @@ def _parse_results_EMCEE(self): if self._yaml_results: self._print_yaml_best_model() + if self._shift_t_0 and self._yaml_results: + print("Plots shift_t_0 : {:}".format(self._shift_t_0_val), + **self._yaml_kwargs) + def _extract_posterior_samples_EMCEE(self): """ set self._samples_flat and self._samples for EMCEE @@ -2486,20 +2512,20 @@ def _shift_t_0_in_samples(self): if name in self._fit_parameters: index = self._fit_parameters.index(name) values = self._samples_flat[:, index] - mean = np.mean(values) + self._shift_t_0_val = int(np.mean(values)) try: - self._samples_flat[:, index] -= int(mean) + self._samples_flat[:, index] -= self._shift_t_0_val if 'trace' in self._plots: - self._samples[:, :, index] -= int(mean) + self._samples[:, :, index] -= self._shift_t_0_val except TypeError: fmt = ("Warning: extremely wide range of posterior {:}: " "from {:} to {:}") warnings.warn( fmt.format(name, np.min(values), np.max(values))) - self._samples_flat[:, index] = values - int(mean) + self._samples_flat[:, index] = values - self._shift_t_0_val if 'trace' in self._plots: self._samples[:, :, index] = ( - self._samples[:, :, index] - int(mean)) + self._samples[:, :, index] - self._shift_t_0_val) def _get_fluxes_to_print_EMCEE(self): """ @@ -2823,6 +2849,9 @@ def _save_figure(self, file_name, figure=None, dpi=None): kwargs = dict() if dpi is not None: kwargs = {'dpi': dpi} + if path.isfile(file_name): + msg = "Existing file " + file_name + " will be overwritten" + warnings.warn(msg) caller.savefig(file_name, **kwargs) plt.close() @@ -3071,18 +3100,71 @@ def _mark_second_Y_axis_in_best_plot(self): magnifications = settings['magnifications'] color = settings.get("color", "red") label = settings.get("label", "magnification") - labels = settings['labels'] + labels = settings.get("labels") ylim = plt.ylim() + ax2 = plt.gca().twinx() + (A_min, A_max, sb_fluxes) = self._second_Y_axis_get_fluxes(ylim) + out1, out2 = False, False + if magnifications == "optimal": + (magnifications, labels, out1) = self._second_Y_axis_optimal( + ax2, A_min, A_max) + flux = sb_fluxes[0] * magnifications + sb_fluxes[1] + out2 = self._second_Y_axis_warnings(flux, labels, magnifications, + A_min, A_max) + if out1 or out2: + ax2.get_yaxis().set_visible(False) + return + + ticks = mm.Utils.get_mag_from_flux(flux) + ax2.set_ylabel(label).set_color(color) + ax2.spines['right'].set_color(color) + ax2.set_ylim(ylim[0], ylim[1]) + ax2.tick_params(axis='y', colors=color) + plt.yticks(ticks, labels, color=color) + + def _second_Y_axis_get_fluxes(self, ylim): + """ + Get fluxes and limiting magnification values for the second Y axis + """ flux_min = mm.Utils.get_flux_from_mag(ylim[0]) flux_max = mm.Utils.get_flux_from_mag(ylim[1]) - (source_flux, blend_flux) = self._event.get_ref_fluxes() - if self._model.n_sources == 1: - total_source_flux = source_flux - else: - total_source_flux = sum(source_flux) - flux = total_source_flux * magnifications + blend_flux + + total_source_flux = sum(source_flux) + A_min = (flux_min - blend_flux) / total_source_flux + A_max = (flux_max - blend_flux) / total_source_flux + + return (A_min, A_max, [total_source_flux, blend_flux]) + + def _second_Y_axis_optimal(self, ax2, A_min, A_max): + """ + Get optimal values of magnifications and labels + """ + ax2.set_ylim(A_min, A_max) + A_values = ax2.yaxis.get_ticklocs().round(7) + A_values = A_values[(A_values >= max(1, A_min)) & (A_values < A_max)] + is_integer = [mag.is_integer() for mag in A_values] + if all(is_integer): + labels = [f"{int(x):d}" for x in A_values] + return (A_values, labels, False) + + fnum = np.array([str(x)[::-1].find(".") for x in A_values]) + labels = np.array([f"%0.{max(fnum)}f" % x for x in A_values]) + if max(fnum) >= 4 and len(fnum[fnum < 4]) < 3: + msg = ("The computed magnifications for the second Y scale cover" + " a range too small to be shown: {:}") + warnings.warn(msg.format(A_values)) + return (A_values, labels.tolist(), True) + if max(fnum) >= 4: + labels = np.array([f"{x:0.3f}" for x in A_values]) + + return (A_values[fnum < 4], labels[fnum < 4].tolist(), False) + + def _second_Y_axis_warnings(self, flux, labels, A_values, A_min, A_max): + """ + Issue warnings for negative flux or bad range of magnificaitons + """ if np.any(flux < 0.): mask = (flux > 0.) flux = flux[mask] @@ -3091,28 +3173,18 @@ def _mark_second_Y_axis_in_best_plot(self): "because they correspond to negative flux which cannot " "be translated to magnitudes.") warnings.warn(msg.format(np.sum(np.logical_not(mask)))) - A_min = (flux_min - blend_flux) / total_source_flux - A_max = (flux_max - blend_flux) / total_source_flux - if (np.min(magnifications) < A_min or np.max(magnifications) > A_max or + if (np.min(A_values) < A_min or np.max(A_values) > A_max or np.any(flux < 0.)): msg = ("Provided magnifications for the second (i.e., right-hand " "side) Y-axis scale are from {:} to {:},\nbut the range " "of plotted magnifications is from {:} to {:}, hence, " "the second scale is not plotted") - args = [min(magnifications), max(magnifications), - A_min[0], A_max[0]] + args = [np.min(A_values), np.max(A_values), A_min, A_max] warnings.warn(msg.format(*args)) - return - - ticks = mm.Utils.get_mag_from_flux(flux) + return True - ax2 = plt.gca().twinx() - ax2.set_ylabel(label).set_color(color) - ax2.spines['right'].set_color(color) - ax2.set_ylim(ylim[0], ylim[1]) - ax2.tick_params(axis='y', colors=color) - plt.yticks(ticks, labels, color=color) + return False def _make_trajectory_plot(self): """