Skip to content

Commit

Permalink
Final improvements to magnifications range, finished issue rpoleski#93
Browse files Browse the repository at this point in the history
  • Loading branch information
@rapoliveira
rapoliveira committed Oct 27, 2023
1 parent 5747f7f commit a1d7d0b
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 31 deletions.
5 changes: 2 additions & 3 deletions examples/example_16/ob03235_2_full.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,9 +86,8 @@ plots:
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. Another option is to set the following option to
# True, so that the magnification range is calculated automatically:
recalculate magnification ticks: True
# what is exact range. Another option is to set it to "optimal", where this range
# will be obtained automatically, but labels cannot be given in this case.
labels: [a, b, c, d, e, f, g, h]
# If you remove the line above, then magnification values will be printed.
label: What is shown on this axis?
Expand Down
75 changes: 47 additions & 28 deletions examples/example_16/ulens_model_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,8 +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')
# 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 "
Expand Down Expand Up @@ -712,25 +712,34 @@ def _check_plots_parameters_best_model_Y_scale(self):
This function assumes that the second Y scale will be plotted.
"""
settings = self._plots['best model']['second Y scale']
allowed = set(['color', 'label', 'labels', 'magnifications',
'recalculate magnification ticks'])
allowed = set(['color', 'label', 'labels', 'magnifications'])
unknown = set(settings.keys()) - allowed
if len(unknown) > 0:
raise ValueError(
'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 '
Expand Down Expand Up @@ -2355,7 +2364,7 @@ def _parse_results_EMCEE(self):
self._print_best_model()
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_value),
**self._yaml_kwargs)
Expand Down Expand Up @@ -3095,14 +3104,34 @@ def _mark_second_Y_axis_in_best_plot(self):
labels = settings['labels']

ylim = plt.ylim()
ax2 = plt.gca().twinx()
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)
A_min = (flux_min - blend_flux) / total_source_flux
A_max = (flux_max - blend_flux) / total_source_flux

if magnifications == "optimal":
ax2.set_ylim(A_min, A_max)
ticks = ax2.yaxis.get_ticklocs()[1:-1]
magnifications = ticks.tolist()
is_integer = [mag.is_integer for mag in magnifications]
if all(is_integer):
labels = [f"%d" % int(x) for x in ticks]
else:
max_n = max([len(str(x))-str(x).find('.')-1 for x in ticks])
labels = [f"%0.{max_n}f" % x for x in ticks]
if max_n > 3:
msg = ("The computed magnifications for the second Y scale"
" cover a range too small to be shown: {:}")
warnings.warn(msg.format(magnifications))
ax2.get_yaxis().set_visible(False)
return

flux = total_source_flux * magnifications + blend_flux
if np.any(flux < 0.):
mask = (flux > 0.)
Expand All @@ -3112,10 +3141,7 @@ 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
ax2 = plt.gca().twinx()


if (np.min(magnifications) < A_min or np.max(magnifications) > A_max or
np.any(flux < 0.)):
msg = ("Provided magnifications for the second (i.e., right-hand "
Expand All @@ -3124,16 +3150,9 @@ def _mark_second_Y_axis_in_best_plot(self):
"the second scale is not plotted")
args = [min(magnifications), max(magnifications),
A_min[0], A_max[0]]
if settings.get("recalculate magnification ticks") is not True:
warnings.warn(msg.format(*args))
ax2.get_yaxis().set_visible(False)
return
else:
ax2.set_ylim(A_min, A_max)
ticks = ax2.yaxis.get_ticklocs()[1:-1]
flux = total_source_flux * ticks.tolist() + blend_flux
max_n = max([len(str(x))-str(x).find('.')-1 for x in ticks])
labels = [f"%0.{max_n}f" % x for x in ticks]
warnings.warn(msg.format(*args))
ax2.get_yaxis().set_visible(False)
return

ticks = mm.Utils.get_mag_from_flux(flux)

Expand Down

0 comments on commit a1d7d0b

Please sign in to comment.