diff --git a/examples/example_16/ob03235_2_full.yaml b/examples/example_16/ob03235_2_full.yaml index 2a8212e2..6ae84dba 100644 --- a/examples/example_16/ob03235_2_full.yaml +++ b/examples/example_16/ob03235_2_full.yaml @@ -37,9 +37,18 @@ starting_parameters: # parameters: t_0 u_0 t_E # See also ob08092-o4_starting_file.yaml fit_constraints: - negative_blending_flux_sigma_mag: 20. + negative_blending_flux_sigma_mag: 20. + #one can specify for which dataset soft constrain on blending flux should be applied: sigma dataset_label(s) + #negative_blending_flux_sigma_mag: 20. "OGLE I-band" "OB03235_MOA.txt" # Alternative sharp constraint: # no_negative_blending_flux: True + #color constrains, where color=source_mag_from_dataset_k-source_mag_from_dataset_m: gauss, mu ,sigma, k, m + color : gauss 0.3 0.01 "OGLE I-band" "OB03235_MOA.txt" + #Alternative for binary source models: + #color source 1 : gauss 0.3 0.01 "OGLE I-band" "OB03235_MOA.txt" + #and/or + #color source 2 : gauss 0.3 0.01 "OGLE I-band" "OB03235_MOA.txt" + prior: t_E: Mroz et al. 2017 # Other possibility: diff --git a/examples/example_16/ulens_model_fit.py b/examples/example_16/ulens_model_fit.py index f8bd511b..571e59db 100644 --- a/examples/example_16/ulens_model_fit.py +++ b/examples/example_16/ulens_model_fit.py @@ -9,6 +9,7 @@ import warnings import math import numpy as np +import shlex from scipy.interpolate import interp1d from matplotlib import pyplot as plt from matplotlib import gridspec, rcParams, rcParamsDefault @@ -38,7 +39,7 @@ except Exception: raise ImportError('\nYou have to install MulensModel first!\n') -__version__ = '0.36.3' +__version__ = '0.37.0' class UlensModelFit(object): @@ -197,14 +198,14 @@ class UlensModelFit(object): saved to temporary files and deleted at the end. ``multimodal`` (*bool*) - do you want multiple modes in - the prosterior to be detected and reported separately? + the posterior to be detected and reported separately? ``n_live_points`` (*int*) - number of live points, default value is 400. ``sampling efficiency`` (*float*) - requested sampling efficiency. MultiNest documentation suggests 0.8 (default value) for parameter - estimatrion and 0.3 for evidence evalutation. + estimation and 0.3 for evidence evaluation. ``evidence tolerance`` (*float*) - requested tolerance of ln(Z) calculation; default is 0.5 and should work well in most cases. @@ -218,10 +219,27 @@ class UlensModelFit(object): blending flux if *True* ``'negative_blending_flux_sigma_mag'`` - impose a prior that - disfavours models with negative blending flux using gaussian prior + disfavors models with negative blending flux using gaussian prior for negative values; the value provided should be on the order of *20.* + ``'color'`` - specify gaussian prior for colors of the sources. + Parameters: + *mean* and *sigma* are floats in magnitudes, *dataset_label* + are str defined in MulensData.plot_properties['label'] + + ``'color source 1'`` - specify gaussian prior for color of + the primary source in binary source model. + Parameters: + *mean* and *sigma* are floats in magnitudes, *dataset_label* + are str defined in MulensData.plot_properties['label'] + + ``'color source 2'`` - specify gaussian prior for color of + the secondary source in binary source model. + Parameters: + *mean* and *sigma* are floats in magnitudes, *dataset_label* + are str defined in MulensData.plot_properties['label'] + ``'prior'`` - specifies the priors for quantities. It's also a *dict*. Possible key-value pairs: @@ -306,13 +324,14 @@ class UlensModelFit(object): These files will have three columns: time, residuals, and uncertainties. """ + def __init__( self, photometry_files, starting_parameters=None, prior_limits=None, model=None, fixed_parameters=None, min_values=None, max_values=None, fitting_parameters=None, fit_constraints=None, plots=None, other_output=None - ): + ): self._check_MM_version() self._photometry_files = photometry_files self._starting_parameters_input = starting_parameters @@ -462,7 +481,7 @@ def _set_default_parameters(self): xi_argument_of_latitude_reference='\\xi_u', xi_eccentricity='\\xi_e', xi_omega_periapsis='\\xi_{\\omege}', - ) + ) self._latex_conversion_other = dict() def _guess_fitting_method(self): @@ -478,7 +497,7 @@ def _guess_fitting_method(self): "Both starting_parameters and prior_limits were defined " "which makes impossible to choose the fitting method. " "These settings indicate EMCEE and pyMultiNest " - "rescpectively, and cannot be both set.") + "respectively, and cannot be both set.") method = "MultiNest" if method is None: raise ValueError( @@ -697,7 +716,7 @@ def _check_plots_parameters_best_model(self): for key in ['legend', 'rcParams', 'second Y scale']: if key in self._plots['best model']: if not isinstance(self._plots['best model'][key], dict): - msg = ('The value of {:} (in best model setttings)' + msg = ('The value of {:} (in best model settings)' 'must be a dictionary, but you provided {:}') args = [key, type(self._plots['best model'][key])] raise TypeError(msg.format(*args)) @@ -868,7 +887,7 @@ def _check_model_parameters(self): def _check_other_fit_parameters(self): """ - Check if there aren't any other inconsistenties between settings + Check if there aren't any other inconsistencies between settings """ if self._fit_method == "MultiNest": if self._min_values is not None or self._max_values is not None: @@ -989,7 +1008,7 @@ def _check_residual_files(self): """ Check if provided names of output files with residuals make sense. We do not check here if the number of files provided is the same - as the number of input datests. + as the number of input datasets. """ existing = [] names = [] @@ -1031,7 +1050,7 @@ def _get_datasets(self): out = '{:} vs {:}'.format( len(self._datasets), len(self._residuals_files)) raise ValueError('The number of datasets and files for ' - 'residuals ouptut do not match: ' + out) + 'residuals output do not match: ' + out) def _get_1_dataset(self, file_, kwargs): """ @@ -1299,7 +1318,7 @@ def _set_dict_safetly(self, target, source, keys_mapping): def _check_output_files_MultiNest(self): """ - Check if output files exist and warn about overwrtting them. + Check if output files exist and warn about overwriting them. If they directory doesn't exist then raise error. """ @@ -1434,7 +1453,7 @@ def _parse_fit_constraints(self): self._priors = None if self._fit_constraints is None: - self._fit_constraints = {"no_negative_blending_flux": False} + self._set_default_fit_constraints() return if isinstance(self._fit_constraints, list): @@ -1444,9 +1463,26 @@ def _parse_fit_constraints(self): "the code. Most probably what you need is:\n" + "fit_constraints = {'no_negative_blending_flux': True}") + self._parse_fit_constraints_keys() + self._parse_fit_constraints_fluxes() + + if 'prior' in self._fit_constraints: + self._parse_fit_constraints_prior() + + def _parse_fit_constraints_keys(self): + """ + Validate the keys in the provided fit_constraints. + """ allowed_keys_flux = { "no_negative_blending_flux", "negative_blending_flux_sigma_mag"} - allowed_keys = {*allowed_keys_flux, "prior"} + + allowed_keys_color = {'color', + 'color source 1', + 'color source 2', } + + allowed_keys = {*allowed_keys_flux, + *allowed_keys_color, + "prior"} used_keys = set(self._fit_constraints.keys()) if len(used_keys - allowed_keys) > 0: raise ValueError('unrecognized constraint: {:}'.format( @@ -1458,13 +1494,106 @@ def _parse_fit_constraints(self): if "no_negative_blending_flux" not in self._fit_constraints: self._fit_constraints["no_negative_blending_flux"] = False - key = "negative_blending_flux_sigma_mag" - if key in used_keys: - self._fit_constraints[key] = mm.Utils.get_flux_from_mag( - self._fit_constraints[key]) + self._check_color_constraints_conflict(allowed_keys_color) - if 'prior' in self._fit_constraints: - self._parse_fit_constraints_prior() + def _set_default_fit_constraints(self): + """ + Set default fitting constraints if none are provided. + + """ + self._fit_constraints = {"no_negative_blending_flux": False} + + def _check_color_constraints_conflict(self, allowed_keys_color): + """ + Check for conflicts among color constraints. + """ + used_keys = set(self._fit_constraints.keys()) + + if len(used_keys.intersection(allowed_keys_color)) >= 2: + if 'color' in used_keys: + raise ValueError( + 'You cannot specify both color and ' + + str(used_keys.intersection(allowed_keys_color)-{'color'})) + + def _parse_fit_constraints_fluxes(self): + """msg += + Process each constraint fit_constraints. + """ + for key, value in self._fit_constraints.items(): + if key == "negative_blending_flux_sigma_mag": + self._parse_fit_constraints_soft_blending(key, value) + elif key in ['color', 'color source 1', 'color source 2']: + self._parse_fit_constraints_color(key, value) + + def _parse_fit_constraints_soft_blending(self, key, value): + """ + Check if soft fit constraint on blending flux are correctly defined. + """ + if isinstance(value, float): + sigma = float(value) + sets = list(range(len(self._datasets))) + + else: + + sigma = float(value.split()[0]) + sets = list(map(self._get_no_of_dataset, + shlex.split(value, posix=False)[1:])) + if len(sets) > len(self._datasets): + raise ValueError( + "dataset number specified in" + + "negative_blending_flux_sigma_mag" + + "do not match with provided datasets") + + self._fit_constraints[key] = [ + mm.Utils.get_flux_from_mag(sigma), sets] + + def _parse_fit_constraints_color(self, key, value): + """ + Check if fit constraint on color are correctly defined. + """ + self._check_unique_datasets_labels() + words = shlex.split(value, posix=False) + + if len(words) != 5 or words[0] != 'gauss': + msg = "Something went wrong in parsing prior for " + msg += "{:}: {:}" + if len(words) == 3 and words[0] == 'gauss': + msg += "color priors require the specification" + msg += "of datasets that should be used for color calculation" + raise ValueError(msg.format(key, value)) + try: + settings = [words[0], float(words[1]), float( + words[2]), words[3], words[4]] + except Exception: + raise ValueError('error in parsing: ' + words[1] + " " + + words[2] + " " + words[3] + " " + words[4]) + if settings[2] < 0.: + raise ValueError('sigma cannot be negative: ' + words[2]) + + settings[3] = self._get_no_of_dataset(settings[3]) + settings[4] = self._get_no_of_dataset(settings[4]) + + if settings[3] == settings[4]: + raise ValueError( + "in " + key + " color have to be from different datasets") + n = len(self._datasets)-1 + if (0 >= settings[3] >= n) or (0 >= settings[4] >= n): + raise ValueError( + "label specified in color prior" + + "do not match with provided datasets") + + self._fit_constraints[key] = settings + + def _check_unique_datasets_labels(self): + """ + Check if the labels of datasets are unique. + """ + labels = [ + dataset.plot_properties['label'] + for dataset in self._datasets + ] + if len(labels) != len(set(labels)): + raise ValueError("Declared labels of datasets must be unique.") def _parse_fit_constraints_prior(self): """ @@ -1494,6 +1623,7 @@ def _parse_fit_constraints_prior(self): if settings[2] < 0.: raise ValueError('sigma cannot be negative: ' + words[2]) priors[key] = settings + else: raise KeyError( "Unrecognized key in fit_constraints/prior: " + key) @@ -1502,6 +1632,30 @@ def _parse_fit_constraints_prior(self): if len(priors) > 0: self._priors = priors + def _get_no_of_dataset(self, label): + """ + Returns the index of a dataset with a specific label. + Parameters : + label: *str* ,*int* + Label of the dataset defined by + MulensData.plot_properties['label'], + or name of the data file if label is not specified, + or a sequential index of the dataset. + + Returns : + index: *int* + Sequential index of the dataset from [0,1,...,n_datasets-1] + + """ + + if '"' in label: + label = label.strip('"') + for (i, dataset) in enumerate(self._datasets): + if dataset.plot_properties['label'] == label: + return i + raise KeyError( + "Unrecognized dataset label in fit_constraints/prior: " + label) + def _read_prior_t_E_data(self): """ read data that specify t_E prior and parse them appropriately @@ -2001,6 +2155,7 @@ def _ln_prior(self, theta): value = self._model.parameters.parameters[parameter] ln_prior += self._get_ln_prior_for_1_parameter( value, prior_settings) + else: raise ValueError('prior not handled: ' + parameter) @@ -2013,8 +2168,8 @@ def _check_valid_Cassan08_trajectory(self): """ sampling = self._model.parameters.uniform_caustic_sampling valid = sampling.check_valid_trajectory( - self._model.parameters.x_caustic_in, - self._model.parameters.x_caustic_out) + self._model.parameters.x_caustic_in, + self._model.parameters.x_caustic_out) return valid def _get_ln_prior_for_1_parameter(self, value, settings): @@ -2114,6 +2269,33 @@ def _get_fluxes(self): return fluxes + def _sumup_inside_prior(self, fluxes, key, inside, index_plus): + """ + Calculates the contribution to the ln_prior + from specified color constraints + Parameters : + fluxes: *array* + Array with fluxes of the current model. + key: *str* + constrain key. + inside: *float* + ln_prior contribution + index_plus: *int* + For a single source, index_plus=0; + for a binary source, index_plus=0 or 1. + Returns : + inside: *float* + Evaluated ln_prior contribution + """ + settings = self._fit_constraints[key] + index1 = (settings[3])*self._n_fluxes_per_dataset + index_plus + index2 = (settings[4])*self._n_fluxes_per_dataset + index_plus + value = mm.Utils.get_mag_from_flux( + fluxes[index1])-mm.Utils.get_mag_from_flux(fluxes[index2]) + inside += self._get_ln_prior_for_1_parameter(value, settings[:-2]) + + return inside + def _run_flux_checks_ln_prior(self, fluxes): """ Run the checks on fluxes - are they in the prior? @@ -2126,15 +2308,47 @@ def _run_flux_checks_ln_prior(self, fluxes): if fluxes[blend_index] < 0.: return outside + inside += self._apply_negative_blending_flux_sigma_mag_prior(fluxes) + inside += self._apply_color_prior(fluxes) + + return inside + + def _apply_negative_blending_flux_sigma_mag_prior(self, fluxes): + """ + Apply the negative blending flux sigma magnitude priotr. + """ + inside = 0.0 key = "negative_blending_flux_sigma_mag" + if key in self._fit_constraints: - blend_index = self._n_fluxes_per_dataset - 1 - if fluxes[blend_index] < 0.: - sigma = self._fit_constraints[key] - inside += -0.5 * (fluxes[blend_index] / sigma)**2 + sigma, datasets = self._fit_constraints[key] + for i, dataset in enumerate(self._datasets): + if i in datasets: + blend_index = ((i + 1) * self._n_fluxes_per_dataset) - 1 + if fluxes[blend_index] < 0.0: + inside += -0.5 * (fluxes[blend_index] / sigma) ** 2 return inside + def _apply_color_prior(self, fluxes): + """ + Apply the color constraints. + """ + inside = 0.0 + key = 'color' + if key in self._fit_constraints: + for i in range(self._n_fluxes_per_dataset - 1): + inside += self._sumup_inside_prior(fluxes, key, inside, i) + + key = 'color source 1' + if key in self._fit_constraints: + inside += self._sumup_inside_prior(fluxes, key, inside, 0) + + key = 'color source 2' + if key in self._fit_constraints: + inside += self._sumup_inside_prior(fluxes, key, inside, 1) + return inside + def _update_best_model_EMCEE(self, ln_prob, theta, fluxes): """ Check if the current model is the best one and save information. @@ -2562,6 +2776,9 @@ def _shift_t_0_in_samples(self): self._samples[:, :, index] = ( self._samples[:, :, index] - self._shift_t_0_val) + if name in self._fixed_parameters.keys(): + self._shift_t_0_val = int(self._fixed_parameters[name]) + def _get_fluxes_to_print_EMCEE(self): """ prepare values to be printed for EMCEE fitting diff --git a/source/MulensModel/data b/source/MulensModel/data deleted file mode 120000 index e67b4559..00000000 --- a/source/MulensModel/data +++ /dev/null @@ -1 +0,0 @@ -../../data \ No newline at end of file diff --git a/source/MulensModel/data b/source/MulensModel/data new file mode 100644 index 00000000..e67b4559 --- /dev/null +++ b/source/MulensModel/data @@ -0,0 +1 @@ +../../data \ No newline at end of file