diff --git a/examples/example_16/ob03235_2_full.yaml b/examples/example_16/ob03235_2_full.yaml index 3be65d64..23780c8e 100644 --- a/examples/example_16/ob03235_2_full.yaml +++ b/examples/example_16/ob03235_2_full.yaml @@ -5,6 +5,14 @@ photometry_files: # - {file_name: some_K2_data.txt, phot_fmt: flux, bandpass: 'Kp', ephemerides_file: K2_ephemeris_01.dat} # For files with time starting with 245...: # - {..., add_2450000: False} +fit_method: EMCEE +# Options: EMCEE, MultiNest, UltraNest. +# If fit_method is not given, EMCEE is used if starting_parameters are given; +# and MultiNest or UltraNest is used if prior_limits are given, depending on +# the combination of fitting_parameters. +# If MultiNest/UltraNest is used, this file cannot have starting_parameters, +# fit_constraints (except negative_blending_flux_sigma_mag and prior->t_E for +# UltraNest), min_values, max_values and plots->trace. model: methods: 2452800. point_source 2452833. VBBL 2452845. point_source 2452860. default method: point_source_point_lens @@ -36,6 +44,17 @@ starting_parameters: # file: ob08092-o4_starting_file_input.txt # parameters: t_0 u_0 t_E # See also ob08092-o4_starting_file.yaml +# prior_limits: +# t_0: [2452847.5, 2452848.5] +# u_0: [0.08, 0.18] +# t_E: 58.7 64.7 +# rho: [0.0004, 0.0018] +# pi_E_N: [-0.05, 0.05] +# pi_E_E: [-0.05, 0.05] +# s: [1.077, 1.125] +# q: [0.001, 0.02] +# alpha: [219.4, 229.4] +# # Only MultiNest and UltraNest, cannot be given with starting_parameters fit_constraints: negative_blending_flux_sigma_mag: 20. #one can specify for which dataset soft constrain on blending flux should be applied: sigma dataset_label(s) @@ -76,12 +95,27 @@ max_values: s: 2. alpha: 360. fitting_parameters: + ## EMCEE only n_walkers: 20 n_steps: 4 n_burn: 2 progress: True posterior file: ob03235_2_models.npy posterior file fluxes: all + ## MultiNest only (basename, multimodal, evidence tolerance=0.5) + # basename: out_ob08092_O4_MN- + # multimodal: True + # evidence tolerance: 0.5 + ## UltraNest only (log directory, derived parameter names, show_status, + ## dlogz=0.5, frac_remain=0.01, max_num_improvement_loops=-1) + # log directory: ultranest_outputs/ + # derived parameter names: flux_s_1 flux_b_1 flux_s_2 flux_b_2 + # show_status: True + # dlogz: 2. + # frac_remain: 0.5 + # max_num_improvement_loops: 0 + ## Both for MultiNest and UltraNest (number of live points) + # n_live_points: 20 plots: best model: # You can skip the line below - the light curve will be plotted on screen. diff --git a/examples/example_16/ob08092-o4_MN.yaml b/examples/example_16/ob08092-o4_MN.yaml index e15f5890..103d526c 100644 --- a/examples/example_16/ob08092-o4_MN.yaml +++ b/examples/example_16/ob08092-o4_MN.yaml @@ -1,5 +1,6 @@ photometry_files: data/OB08092/phot_ob08092_O4.dat +fit_method: MultiNest prior_limits: t_0: [2455379.4, 2455379.76] u_0: [0.3, 0.65] @@ -7,11 +8,14 @@ prior_limits: fitting_parameters: basename: out_ob08092_O4_MN- multimodal: True + # Default settings of other parameters: + # evidence tolerance: 0.5 n_live_points: 1000 -# Default settings of other parameters: -# evidence tolerance: 0.5 plots: best model: file: out_ob08092_O4_MN_model.png triangle: file: out_ob08092_O4_MN_triangle.png +other_output: + yaml output: + file name: ob08092_O4_MN_all_results.yaml diff --git a/examples/example_16/ob08092-o4_UN.yaml b/examples/example_16/ob08092-o4_UN.yaml new file mode 100644 index 00000000..2cf18914 --- /dev/null +++ b/examples/example_16/ob08092-o4_UN.yaml @@ -0,0 +1,35 @@ +photometry_files: + data/OB08092/phot_ob08092_O4.dat +fit_method: UltraNest +prior_limits: + t_0: [2455379.4, 2455379.76] + u_0: [0.3, 0.65] + t_E: 16. 19.6 +fit_constraints: + negative_blending_flux_sigma_mag: 20. + prior: + t_E: Mroz et al. 2020 +fitting_parameters: + log directory: outputs_ultranest/ + derived parameter names: flux_s_1 flux_b_1 + show_status: True + n_live_points: 1000 + # `n_live_points` can also be named `min_num_live_points` + # If it is smaller than 40, `cluster_num_live_points` is also reduced. + # UltraNest may increase n_live_points if it is too low to achieve the + # logz accuracy (default=0.5). It can be avoided increasing dlogz: + dlogz: 2. + # The parameters below can reduce runtime (default is 0.01 and -1) + # frac_remain: 0.5 + # max_num_improvement_loops: 0 +plots: + best model: + file: out_ob08092_O4_UN_model.png + second Y scale: + magnifications: optimal + triangle: + file: out_ob08092_O4_UN_triangle.png + shift t_0: False +other_output: + yaml output: + file name: ob08092_O4_UN_all_results.yaml diff --git a/examples/example_16/ob08092-o4_minimal_MN.yaml b/examples/example_16/ob08092-o4_minimal_MN.yaml index 9a998a07..7e95e488 100644 --- a/examples/example_16/ob08092-o4_minimal_MN.yaml +++ b/examples/example_16/ob08092-o4_minimal_MN.yaml @@ -1,5 +1,7 @@ photometry_files: data/OB08092/phot_ob08092_O4.dat +fit_method: MultiNest +# fit_method options if prior_limits are given: MultiNest and UltraNest prior_limits: t_0: [2455379.4, 2455379.76] u_0: [0.46, 0.65] diff --git a/examples/example_16/ulens_model_fit.py b/examples/example_16/ulens_model_fit.py index 3beb2495..5bdfe6d2 100644 --- a/examples/example_16/ulens_model_fit.py +++ b/examples/example_16/ulens_model_fit.py @@ -33,13 +33,17 @@ from pymultinest.analyse import Analyzer except Exception: import_failed.add("pymultinest") +try: + import ultranest +except Exception: + import_failed.add("ultranest") try: import MulensModel as mm except Exception: raise ImportError('\nYou have to install MulensModel first!\n') -__version__ = '0.38.0' +__version__ = '0.39.0' class UlensModelFit(object): @@ -69,6 +73,19 @@ class UlensModelFit(object): ``'scale_errorbars': {'factor': kappa, 'minimum': epsilon}`` to scale uncertainties. + fit_method: *str* + Method of fitting. Currently accepted values are ``EMCEE``, + ``MultiNest``, and ``UltraNest``. If not provided, the script + will guess it based on other parameters: ``EMCEE`` is selected + if ``starting_parameters`` are provided. If ``prior_limits`` are + provided, either ``MultiNest`` or ``UltraNest`` will be selected + depending on the ``fitting_parameters``. + + Webpage of each method: + - EMCEE: https://emcee.readthedocs.io/en/stable/ + - MultiNest: https://johannesbuchner.github.io/PyMultiNest/ + - UltraNest: https://johannesbuchner.github.io/UltraNest/ + starting_parameters: *dict* Starting values of the parameters. It also indicates the EMCEE fitting mode. @@ -108,12 +125,12 @@ class UlensModelFit(object): prior_limits: *dict* Upper and lower limits of parameters. - It also indicates the pyMultiNest fitting mode. + It only applies to pyMultiNest and UltraNest fitting. Keys are MulensModel parameters and values are lists of two floats each (alternatively a string giving 2 floats can be provided - see example below). Currently, no informative priors are allowed for - pyMultiNest fitting. Example input: + pyMultiNest and UltraNest fitting. Example input: .. code-block:: python @@ -157,7 +174,7 @@ class UlensModelFit(object): min_values: *dict* Minimum values of parameters that define the prior, e.g., ``{'t_E': 0.}``. Note that the these are only limits of a prior. - Functional form of priors can be defines in ``fit_constraints``. + Functional form of priors can be defined in ``fit_constraints``. It works only for EMCEE fitting. max_values: *dict* @@ -167,7 +184,7 @@ class UlensModelFit(object): fitting_parameters: *dict* Parameters of the fit function. They depend on the method used - - we discuss EMCEE and pyMultiNest below. + we discuss EMCEE, pyMultiNest and UltraNest below. First - EMCEE. The required parameter is ``n_steps``. You can also specify ``n_burn`` and ``n_walkers``. The ``n_burn`` @@ -201,7 +218,7 @@ class UlensModelFit(object): the posterior to be detected and reported separately? ``n_live_points`` (*int*) - number of live points, default value - is 400. + is 400. Also valid for UltraNest. ``sampling efficiency`` (*float*) - requested sampling efficiency. MultiNest documentation suggests 0.8 (default value) for parameter @@ -210,6 +227,37 @@ class UlensModelFit(object): ``evidence tolerance`` (*float*) - requested tolerance of ln(Z) calculation; default is 0.5 and should work well in most cases. + Third - UltraNest. There are no required parameters, but a few + can be provided. Currently accepted ones are: + + ``log directory`` (*str*) - where to store output files. If given, + there is a check if directory exists. If not given, no outputs + are saved. + + ``derived parameter names`` (*str*) - names of additional derived + parameters created by transform. In microlensing, they are usually + the source(s) and blending fluxes. If not given, they are ignored + in the transform function. + + ``show_status`` (*bool*) - whether to show integration progress + as a status line or not. Default is *True*. + + ``min_num_live_points`` (*int*) - minimum number of live points + throughout the run. Default value is 400. + + ``dlogz`` (*float*) - Target evidence uncertainty, in order to + obtain a logz error below a threshold. Default value is 0.5. + It can be increased to allow `min_num_live_points` values below: + sqrt(iterations) / dlogz = sqrt(1000) / 0.5 ~ 64. + + ``frac_remain`` (*float*) - Integrate until this fraction of the + integral is left in the remainder. Numbers smaller than 0.01 + ensure that peaks are discovered, higher numbers can be set if + the posterior is simple. Default value is 0.01. + + ``max_num_improvement_loops`` (*int*) - Limit the number of times + the algorithm is repeated to improve. Default value is -1. + fit_constraints: *dict* Constraints on model other than minimal and maximal values. @@ -337,7 +385,8 @@ def __init__( 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 + fit_constraints=None, plots=None, other_output=None, + fit_method=None ): self._check_MM_version() self._photometry_files = photometry_files @@ -351,11 +400,16 @@ def __init__( self._fit_constraints = fit_constraints self._plots = plots self._other_output = other_output + self._fit_method = fit_method self._which_task() self._set_default_parameters() if self._task == 'fit': - self._guess_fitting_method() + if self._fit_method is None: + self._guess_fitting_method() + else: + self._fit_method = self._fit_method.lower() + self._check_fitting_method() self._check_starting_parameters_type() self._set_fit_parameters_unsorted() self._check_imports() @@ -419,7 +473,7 @@ def _which_task(self): def _check_unnecessary_settings_plot(self): """ - Make sure that there arent' too many parameters specified for: + Make sure that there aren't too many parameters specified for: self._task == 'plot' """ keys = ['_starting_parameters_input', '_min_values', '_max_values', @@ -506,7 +560,7 @@ def _guess_fitting_method(self): "which makes impossible to choose the fitting method. " "These settings indicate EMCEE and pyMultiNest " "respectively, and cannot be both set.") - method = "MultiNest" + method = self._guess_MultiNest_or_UltraNest() if method is None: raise ValueError( "No fitting method chosen. You can chose either 'EMCEE' or " @@ -514,12 +568,50 @@ def _guess_fitting_method(self): "starting_parameters or prior_limits, respectively.") self._fit_method = method + def _guess_MultiNest_or_UltraNest(self): + """ + Guess fit_method between MultiNest or UltraNest, based on the + provided fitting_parameters. + """ + args_MultiNest = ['basename', 'multimodal', 'evidence tolerance', + 'n_live_points'] + if all([key in args_MultiNest for key in self._fitting_parameters]): + return "MultiNest" + + args_UltraNest = ['log directory', 'derived parameter names', + 'show_status', 'dlogz', 'frac_remain', + 'max_num_improvement_loops', 'n_live_points'] + if all([key in args_UltraNest for key in self._fitting_parameters]): + return "UltraNest" + + raise ValueError( + "Cannot guess fitting method. Provide more parameters in " + "fitting_parameters.") + + def _check_fitting_method(self): + """ + Check if fitting method is consistent with the settings. + """ + if self._fit_method == "emcee": + self._fit_method = "EMCEE" + if self._starting_parameters_input is None: + raise ValueError( + "EMCEE fitting method requires starting_parameters.") + elif self._fit_method in ["multinest", "ultranest"]: + self._fit_method = self._fit_method.capitalize() + self._fit_method = self._fit_method.replace("nest", "Nest") + if self._prior_limits is None: + msg = "{:} fitting method requires prior_limits." + raise ValueError(msg.format(self._fit_method)) + else: + raise ValueError("Invalid fitting method was inserted.") + def _check_starting_parameters_type(self): """ Check if starting parameters are read from file or will be drawn from distributions specified. """ - if self._fit_method == "MultiNest": + if self._fit_method in ["MultiNest", "UltraNest"]: return if 'file' in self._starting_parameters_input: @@ -547,7 +639,7 @@ def _set_fit_parameters_unsorted(self): else: raise ValueError( 'unexpected: ' + str(self._starting_parameters_type)) - elif self._fit_method == "MultiNest": + elif self._fit_method in ["MultiNest", "UltraNest"]: unsorted_keys = self._prior_limits.keys() else: raise ValueError('unexpected method error') @@ -576,6 +668,8 @@ def _check_imports(self): required_packages.add('emcee') elif self._fit_method == "MultiNest": required_packages.add('pymultinest') + elif self._fit_method == "UltraNest": + required_packages.add('ultranest') if self._plots is not None and 'triangle' in self._plots: required_packages.add('corner') @@ -612,7 +706,6 @@ def run_fit(self): self._check_ulens_model_parameters() self._get_parameters_ordered() self._get_parameters_latex() - self._parse_fitting_parameters() self._set_prior_limits() self._parse_fit_constraints() if self._fit_method == "EMCEE": @@ -620,6 +713,7 @@ def run_fit(self): self._check_fixed_parameters() self._make_model_and_event() + self._parse_fitting_parameters() if self._fit_method == "EMCEE": self._get_starting_parameters() @@ -844,9 +938,9 @@ def _check_plots_parameters_trace(self): raise ValueError( 'Unknown settings for "trace" plot: {:}'.format(unknown)) - if self._fit_method == "MultiNest": + if self._fit_method in ["MultiNest", "UltraNest"]: raise ValueError( - 'Trace plot cannot be requested for MultiNest fit') + f'Trace plot cannot be requested for {self._fit_method}.') self._parse_plots_parameter_shift_t_0(self._plots['trace']) @@ -897,10 +991,10 @@ def _check_other_fit_parameters(self): """ Check if there aren't any other inconsistencies between settings """ - if self._fit_method == "MultiNest": + if self._fit_method in ["MultiNest", "UltraNest"]: if self._min_values is not None or self._max_values is not None: - raise ValueError("In MultiNest fitting you cannot set " - "min_values or max_values") + msg = "In {:} fitting you cannot set min_values or max_values" + raise ValueError(msg.format(self._fit_method)) def _parse_methods(self, methods): """ @@ -1127,11 +1221,12 @@ def _get_parameters_latex(self): if self._fit_constraints is not None: if 'posterior parsing' in self._fit_constraints: - if 'abs' in self._fit_constraints['posterior parsing']: - settings = self._fit_constraints['posterior parsing']['abs'] - if not isinstance(settings, list): - raise ValueError("Error: fit_constraints -> posterior parsing -> abs - list expected") - for key in self._fit_constraints['posterior parsing']['abs']: + settings = self._fit_constraints['posterior parsing'] + if 'abs' in settings: + if not isinstance(settings['abs'], list): + raise ValueError("Error: fit_constraints -> posterior" + " parsing -> abs - list expected") + for key in settings['abs']: conversion[key] = "|" + conversion[key] + "|" self._fit_parameters_latex = [ @@ -1147,6 +1242,8 @@ def _parse_fitting_parameters(self): self._get_n_walkers() elif self._fit_method == 'MultiNest': self._parse_fitting_parameters_MultiNest() + elif self._fit_method == 'UltraNest': + self._parse_fitting_parameters_UltraNest() else: raise ValueError('internal inconsistency') @@ -1220,8 +1317,8 @@ def _check_required_and_allowed_parameters(self, required, allowed): for required_ in required: if required_ not in settings: - raise ValueError('EMCEE method requires fitting parameter: ' + - required_) + msg = '{:} method requires fitting parameter: {:}' + raise ValueError(msg.format(self._fit_method, required_)) if len(set(settings.keys()) - set(full)) > 0: raise ValueError('Unexpected fitting parameters: ' + @@ -1275,11 +1372,9 @@ def _get_n_walkers(self): if 'n_walkers' in self._fitting_parameters: self._n_walkers = self._fitting_parameters['n_walkers'] else: - if self._starting_parameters_type == 'file': - self._n_walkers = None - elif self._starting_parameters_type == 'draw': + if self._starting_parameters_type == 'draw': self._n_walkers = 4 * self._n_fit_parameters - else: + elif self._starting_parameters_type != 'file': raise ValueError( 'Unexpected: ' + self._starting_parameters_type) @@ -1311,7 +1406,7 @@ def _parse_fitting_parameters_MultiNest(self): same_keys = ["multimodal", "n_live_points"] keys = {**keys, **{key: key for key in same_keys}} - self._set_dict_safetly(self._kwargs_MultiNest, settings, keys) + self._set_dict_safely(self._kwargs_MultiNest, settings, keys) self._kwargs_MultiNest['importance_nested_sampling'] = ( not self._kwargs_MultiNest['multimodal']) @@ -1323,7 +1418,7 @@ def _parse_fitting_parameters_MultiNest(self): self._MN_temporary_files = True self._check_output_files_MultiNest() - def _set_dict_safetly(self, target, source, keys_mapping): + def _set_dict_safely(self, target, source, keys_mapping): """ For each key in keys_mapping (*dict*) check if it is in source (*dict*). If it is, then set @@ -1368,13 +1463,66 @@ def _check_output_files_MultiNest(self): message += "(unless you kill this process)!!!\n" warnings.warn(message + str(existing) + "\n") + def _parse_fitting_parameters_UltraNest(self): + """ + Make sure UltraNest fitting parameters are properly defined + """ + self._kwargs_UltraNest = dict() + self._kwargs_UltraNest['viz_callback'] = False + + settings = self._fitting_parameters.copy() + if settings is None: + settings = dict() + + required = [] + bools = ['show_status'] + ints = ['min_num_live_points', 'max_num_improvement_loops'] + if 'n_live_points' in settings: + ints[0] = 'n_live_points' + strings = ['log directory', 'derived parameter names'] + floats = ['dlogz', 'frac_remain'] + allowed = strings + bools + ints + floats + + self._check_required_and_allowed_parameters(required, allowed) + self._check_parameters_types(settings, bools, ints, floats, strings) + self._log_dir_UltraNest = settings.pop("log directory", None) + value = settings.pop("derived parameter names", "") + self._derived_param_names_UltraNest = value.split() + self._check_dir_and_parameter_names_Ultranest() + + keys = {"n_live_points": "min_num_live_points"} + same_keys = ["min_num_live_points", 'max_num_improvement_loops', + "show_status", "dlogz", "frac_remain"] + keys = {**keys, **{key: key for key in same_keys}} + self._set_dict_safely(self._kwargs_UltraNest, settings, keys) + + def _check_dir_and_parameter_names_Ultranest(self): + """ + Checks if the path to `log directory` exists and is a directory, + and also if the number of `derived parameter names` matches the + number of derived fluxes. + """ + if self._log_dir_UltraNest is not None: + if not path.exists(self._log_dir_UltraNest): + raise ValueError("log directory value in fitting_parameters" + "does not exist.") + elif not path.isdir(self._log_dir_UltraNest): + raise ValueError("log directory value in fitting_parameters" + "exists, but it is a file.") + + n_datasets = len(self._datasets) + n_fluxes = self._n_fluxes_per_dataset + if len(self._derived_param_names_UltraNest) != n_datasets * n_fluxes: + raise ValueError("The number of `derived parameter names` must " + "match the number of derived fluxes.") + def _set_prior_limits(self): """ Set minimum and maximum values of the prior space """ if self._fit_method == 'EMCEE': self._set_prior_limits_EMCEE() - elif self._fit_method == 'MultiNest': + elif self._fit_method in ['MultiNest', 'UltraNest']: self._set_prior_limits_MultiNest() else: raise ValueError('internal bug') @@ -1430,7 +1578,7 @@ def _set_prior_limits_MultiNest(self): max_values = [] for parameter in self._fit_parameters: if parameter not in self._prior_limits: - raise ValueError("interal issue") + raise ValueError("internal issue") values = self._prior_limits[parameter] if isinstance(values, str): values = values.split() @@ -1459,13 +1607,6 @@ def _parse_fit_constraints(self): """ Parse the fitting constraints that are not simple limits on parameters """ - if self._fit_method == 'MultiNest': - if self._fit_constraints is not None: - raise NotImplementedError( - "Currently no fit_constraints are implemented for " - "MultiNest fit. Please contact Radek Poleski with " - "a specific request.") - self._prior_t_E = None self._priors = None @@ -1485,15 +1626,27 @@ def _check_fit_constraints(self): """ Run checks on self._fit_constraints """ - if self._fit_constraints is not None and self._fit_method == 'MultiNest': + if self._fit_method == 'MultiNest': raise NotImplementedError( "Currently no fit_constraints are implemented for MultiNest " "fit. Please contact Radek Poleski with a specific request.") + if self._fit_method == 'UltraNest': + allowed_keys = {'negative_blending_flux_sigma_mag', 'prior'} + used_keys = set(self._fit_constraints.keys()) + if not used_keys.issubset(allowed_keys): + raise NotImplementedError( + "The supported fit_constraints options for UltraNest are" + " `negative_blending_flux_sigma_mag` and `prior`.") + if 'prior' in used_keys: + if set(self._fit_constraints['prior'].keys()) != {'t_E'}: + raise ValueError( + "Only `t_E` is allowed in fit_constraints['prior'].") + if isinstance(self._fit_constraints, list): raise TypeError( "In version 0.5.0 we've changed type of 'fit_constraints' " + - "from list to dict. Please correct you input and re-run " + + "from list to dict. Please correct your input and re-run " + "the code. Most probably what you need is:\n" + "fit_constraints = {'no_negative_blending_flux': True}") @@ -1540,7 +1693,7 @@ def _check_color_constraints_conflict(self, allowed_keys_color): 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(): @@ -1626,10 +1779,8 @@ def _parse_fit_constraints_prior(self): priors = dict() for (key, value) in self._fit_constraints['prior'].items(): if key == 't_E': - if value == "Mroz et al. 2017": - self._prior_t_E = 'Mroz+17' - elif value == "Mroz et al. 2020": - self._prior_t_E = 'Mroz+20' + if value in ["Mroz et al. 2017", "Mroz et al. 2020"]: + self._prior_t_E = value.replace(" et al. 20", "+") else: raise ValueError("Unrecognized t_E prior: " + value) self._read_prior_t_E_data() @@ -1661,18 +1812,19 @@ def _parse_fit_constraints_posterior(self): Parse constraints on what is done with posterior. """ if 'posterior parsing' not in self._fit_constraints: + self._parse_posterior_abs = list() return if self._fit_method != "EMCEE": - raise ValueError('Input in "posterior parsing" is allowed only for EMCEE') + raise ValueError('Input in "posterior parsing" is allowed only' + ' for EMCEE') allowed_keys = {"abs"} settings = self._fit_constraints['posterior parsing'] unknown = set(settings.keys()) - allowed_keys if len(unknown) > 0: - raise KeyError( - "Unrecognized key in fit_constraints -> 'posterior parsing': " + - str(unknown)) + msg = "Unrecognized key in fit_constraints -> 'posterior parsing':" + raise KeyError(msg + ' ' + str(unknown)) if 'abs' in settings: self._parse_posterior_abs = settings['abs'] @@ -1947,7 +2099,7 @@ def _get_example_parameters(self): elif self._task == 'fit': if self._fit_method == 'EMCEE': parameters.update(self._get_example_parameters_EMCEE()) - elif self._fit_method == 'MultiNest': + elif self._fit_method in ['MultiNest', 'UltraNest']: means = 0.5 * (self._max_values + self._min_values) parameters.update(dict(zip(self._fit_parameters, means))) if "x_caustic_in" in self._fit_parameters: @@ -2105,9 +2257,12 @@ def _ln_prob(self, theta): NOTE: we're using np.log(), i.e., natural logarithms. """ - ln_prior = self._ln_prior(theta) - if not np.isfinite(ln_prior): - return self._return_ln_prob(-np.inf) + if self._fit_method == "EMCEE": + ln_prior = self._ln_prior(theta) + if not np.isfinite(ln_prior): + return self._return_ln_prob(-np.inf) + elif self._fit_method == "UltraNest": + ln_prior = self._ln_prior_t_E() if self._prior_t_E else 0. ln_like = self._ln_like(theta) if not np.isfinite(ln_like): @@ -2123,9 +2278,11 @@ def _ln_prob(self, theta): ln_prob += ln_prior_flux - self._update_best_model_EMCEE(ln_prob, theta, fluxes) + if self._fit_method == "EMCEE": + self._update_best_model_EMCEE(ln_prob, theta, fluxes) + return self._return_ln_prob(ln_prob, fluxes) - return self._return_ln_prob(ln_prob, fluxes) + return ln_prob def _return_ln_prob(self, value, fluxes=None): """ @@ -2359,7 +2516,8 @@ def _run_flux_checks_ln_prior(self, fluxes): return outside inside += self._apply_negative_blending_flux_sigma_mag_prior(fluxes) - inside += self._apply_color_prior(fluxes) + if self._fit_method == "EMCEE": + inside += self._apply_color_prior(fluxes) return inside @@ -2418,6 +2576,8 @@ def _setup_fit(self): self._setup_fit_EMCEE() elif self._fit_method == 'MultiNest': self._setup_fit_MultiNest() + elif self._fit_method == 'UltraNest': + self._setup_fit_UltraNest() else: raise ValueError('internal bug') @@ -2441,21 +2601,47 @@ def _setup_fit_MultiNest(self): if self._return_fluxes: self._kwargs_MultiNest['n_params'] += self._n_fluxes + def _setup_fit_UltraNest(self): + """ + Prepare UltraNest fit, declaring sampler instance. + If the names of the derived parameters are not given, the source + and blending fluxes are assigned, with indexes depending on the + number of sources (s1, s2) and datasets (_1, _2). + """ + if self._return_fluxes: + if len(self._derived_param_names_UltraNest) == 0: + if self._flux_names is None: + self._flux_names = self._get_fluxes_names_to_print() + self._derived_param_names_UltraNest = self._flux_names + + n_dims = self._n_fit_parameters + n_params = n_dims + self._n_fluxes_per_dataset * self._return_fluxes + t_kwargs = {'n_dims': n_dims, 'n_params': n_params} + self._sampler = ultranest.ReactiveNestedSampler( + self._fit_parameters, self._ln_prob, + transform=lambda cube: self._transform_unit_cube(cube, **t_kwargs), + derived_param_names=self._derived_param_names_UltraNest, + log_dir=self._log_dir_UltraNest + ) + def _transform_unit_cube(self, cube, n_dims, n_params): """ - Transform MulitNest unit cube to microlensing parameters. + Transform MultiNest/UltraNest unit cube to microlensing parameters. - Based on SafePrior() in + MultiNest: based on SafePrior() in https://github.com/JohannesBuchner/PyMultiNest/blob/master/ pymultinest/solve.py + UltraNest: based on the above and UltraNest documentation + https://johannesbuchner.github.io/UltraNest/example-sine-line.html NOTE: We call self._ln_like() here (and remember the result) because in MultiNest you can add fluxes only in "prior" function, not in likelihood function. """ cube_out = self._min_values + cube[:n_dims] * self._range_values - for i in range(n_dims): - cube[i] = cube_out[i] + if self._fit_method == "MultiNest": + for i in range(n_dims): + cube[i] = cube_out[i] if "x_caustic_in" in self._model.parameters.parameters: self._set_model_parameters(cube_out) @@ -2463,19 +2649,25 @@ def _transform_unit_cube(self, cube, n_dims, n_params): self._last_ln_like = -1.e300 self._last_theta = cube_out if self._return_fluxes: - for i in range(n_dims, n_params): - cube[i] = 0. self._last_fluxes = np.zeros(n_params - n_dims) - return + if self._fit_method == "MultiNest": + for i in range(n_dims, n_params): + cube[i] = 0. + return + cube_out = np.append(cube_out, self._last_fluxes) + return cube_out self._last_ln_like = self._ln_like(cube_out) self._last_theta = cube_out if self._return_fluxes: fluxes = self._get_fluxes() + self._last_fluxes = fluxes + if self._fit_method == "UltraNest": + cube_out = np.append(cube_out, fluxes) + return cube_out for i in range(n_dims, n_params): cube[i] = fluxes[i-n_dims] - self._last_fluxes = fluxes def _ln_like_MN(self, theta, n_dim, n_params, lnew): """ @@ -2512,6 +2704,8 @@ def _run_fit(self): self._run_fit_EMCEE() elif self._fit_method == 'MultiNest': self._run_fit_MultiNest() + elif self._fit_method == 'UltraNest': + self._run_fit_UltraNest() else: raise ValueError('internal bug') @@ -2527,6 +2721,16 @@ def _run_fit_MultiNest(self): """ mn_run(**self._kwargs_MultiNest) + def _run_fit_UltraNest(self): + """ + Run Ultranest fit + """ + min_n_live = self._kwargs_UltraNest.get("min_num_live_points", 400) + cluster_n_live = 40 if min_n_live >= 40 else min_n_live + self._kwargs_UltraNest['cluster_num_live_points'] = cluster_n_live + + self._result_UltraNest = self._sampler.run(**self._kwargs_UltraNest) + def _finish_fit(self): """ Make the things that are necessary after the fit is done. @@ -2602,6 +2806,8 @@ def _parse_results(self): self._save_posterior_EMCEE() elif self._fit_method == "MultiNest": self._parse_results_MultiNest() + elif self._fit_method == "UltraNest": + self._parse_results_UltraNest() else: raise ValueError('internal bug') @@ -2698,7 +2904,7 @@ def _print_results(self, data, names="parameters", mode=None): if self._fit_method == "EMCEE": results = self._get_weighted_percentile(data) - elif self._fit_method == "MultiNest": + elif self._fit_method in ["MultiNest", "UltraNest"]: if mode is None: weights = self._samples_flat_weights else: @@ -2747,7 +2953,7 @@ def _print_yaml_results(self, data, names="parameters", mode=None): if self._fit_method == "EMCEE": results = self._get_weighted_percentile(data) - elif self._fit_method == "MultiNest": + elif self._fit_method in ["MultiNest", "UltraNest"]: if mode is None: weights = self._samples_flat_weights else: @@ -3061,12 +3267,45 @@ def _extract_posterior_samples_MultiNest(self): def _get_fluxes_to_print_MultiNest(self): """ - prepare values to be printed for EMCEE fitting + prepare flux values to be printed for MultiNest and Ultranest fitting """ - index = 2 + self._n_fit_parameters - data = self._analyzer_data[:, index:] + if self._fit_method == "MultiNest": + index = 2 + self._n_fit_parameters + data = self._analyzer_data + elif self._fit_method == "UltraNest": + index = self._n_fit_parameters + data = self._result_UltraNest['weighted_samples']['points'] + + return data[:, index:] + + def _parse_results_UltraNest(self): + """ + Parse results of UltraNest fitting. + Functions that print and save EMCEE results are also called here. + """ + # re-weighted posterior samples: + # self._samples_flat = self._result_UltraNest['samples'][:, :-2] + # weighted samples from the posterior: + weighted_samples = self._result_UltraNest['weighted_samples'] + index = self._n_fit_parameters + self._samples_flat = weighted_samples['points'][:, :index] + self._samples_flat_weights = weighted_samples['weights'] + self._sampler.print_results() + + max_like = self._result_UltraNest['maximum_likelihood'] + self._best_model_ln_prob = max_like['logl'] + self._best_model_theta = max_like['point'][:self._n_fit_parameters] + self._best_model_fluxes = max_like['point'][self._n_fit_parameters:] + self._parse_results_MultiNest_singlemode() - return data + self._shift_t_0_in_samples() + self._print_best_model() + if self._yaml_results: + self._print_yaml_best_model() + ln_ev = self._result_UltraNest['logz_single'] + ln_ev_err = self._result_UltraNest['logzerr_single'] + lns = " ln_ev: [{:.5f}, +{:.5f}, -{:.5f}]" + print(lns.format(ln_ev, ln_ev_err, ln_ev_err), **self._yaml_kwargs) def _write_residuals(self): """