diff --git a/doc/conf.py b/doc/conf.py index c96cf87b0..0cca2da1d 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -73,6 +73,8 @@ autodoc_mock_imports = ["amici"] autodoc_class_signature = "separated" +# napoleon options +napoleon_use_rtype = False # links for intersphinx intersphinx_mapping = { diff --git a/pypesto/profile/approximate.py b/pypesto/profile/approximate.py index ccddd99d5..34363691d 100644 --- a/pypesto/profile/approximate.py +++ b/pypesto/profile/approximate.py @@ -50,8 +50,7 @@ def approximate_parameter_profile( Returns ------- - result: - The profile results are filled into `result.profile_result`. + The profile results are filled into `result.profile_result`. """ # Handling defaults # profiling indices diff --git a/pypesto/profile/profile.py b/pypesto/profile/profile.py index 116ef21d8..2d0fcceb6 100644 --- a/pypesto/profile/profile.py +++ b/pypesto/profile/profile.py @@ -30,7 +30,7 @@ def parameter_profile( overwrite: bool = False, ) -> Result: """ - Call to do parameter profiling. + Compute parameter profiles. Parameters ---------- @@ -45,6 +45,7 @@ def parameter_profile( The optimizer to be used along each profile. engine: The engine to be used. + Defaults to :class:`pypesto.engine.SingleCoreEngine`. profile_index: List with the parameter indices to be profiled (by default all free indices). @@ -61,12 +62,13 @@ def parameter_profile( :func:`pypesto.profile.profile_next_guess.next_guess`. profile_options: Various options applied to the profile optimization. + See :class:`pypesto.profile.options.ProfileOptions`. progress_bar: Whether to display a progress bar. filename: Name of the hdf5 file, where the result will be saved. Default is None, which deactivates automatic saving. If set to - "Auto" it will automatically generate a file named + ``Auto`` it will automatically generate a file named ``year_month_day_profiling_result.hdf5``. Optionally a method, see docs for :func:`pypesto.store.auto.autosave`. overwrite: @@ -75,8 +77,7 @@ def parameter_profile( Returns ------- - result: - The profile results are filled into `result.profile_result`. + The profile results are filled into `result.profile_result`. """ # Copy the problem to avoid side effects problem = copy.deepcopy(problem) diff --git a/pypesto/profile/profile_next_guess.py b/pypesto/profile/profile_next_guess.py index 7a2cf5364..1d3d89f1b 100644 --- a/pypesto/profile/profile_next_guess.py +++ b/pypesto/profile/profile_next_guess.py @@ -1,5 +1,4 @@ -import copy -from typing import Callable, List, Literal, Tuple, Union +from typing import Callable, Literal import numpy as np @@ -42,10 +41,12 @@ def next_guess( profile_options: Various options applied to the profile optimization. update_type: - Type of update for next profile point: - ``fixed_step`` (see :func:`fixed_step`), - ``adaptive_step_order_0``, ``adaptive_step_order_1``, or ``adaptive_step_regression`` - (see :func:`adaptive_step`). + Type of update for next profile point. Available options are: + + * ``fixed_step`` (see :func:`fixed_step`) + * ``adaptive_step_order_0`` (see :func:`adaptive_step`). + * ``adaptive_step_order_1`` (see :func:`adaptive_step`). + * ``adaptive_step_regression`` (see :func:`adaptive_step`). current_profile: The profile which should be computed. problem: @@ -55,8 +56,7 @@ def next_guess( Returns ------- - next_guess: - The next initial guess as base for the next profile point. + The next initial guess as base for the next profile point. """ if update_type == 'fixed_step': return fixed_step( @@ -96,16 +96,16 @@ def fixed_step( """Most simple method to create the next guess. Computes the next point based on the fixed step size given by - ``default_step_size`` in :class:`ProfileOptions`. + :attr:`pypesto.profile.ProfileOptions.default_step_size`. Parameters ---------- x: The current position of the profiler, size `dim_full`. par_index: - The index of the parameter of the current profile + The index of the parameter of the current profile. par_direction: - The direction, in which the profiling is done (``1`` or ``-1``) + The direction, in which the profiling is done (``1`` or ``-1``). options: Various options applied to the profile optimization. problem: @@ -113,8 +113,7 @@ def fixed_step( Returns ------- - x_new: - The updated parameter vector, of size `dim_full`. + The updated parameter vector, of size `dim_full`. """ delta_x = np.zeros(len(x)) delta_x[par_index] = par_direction * options.default_step_size @@ -150,41 +149,42 @@ def adaptive_step( x: The current position of the profiler, size `dim_full`. par_index: - The index of the parameter of the current profile + The index of the parameter of the current profile. par_direction: - The direction, in which the profiling is done (1 or -1) + The direction, in which the profiling is done (``1`` or ``-1``). options: Various options applied to the profile optimization. current_profile: - The profile which should be computed + The profile which should be computed. problem: The problem to be solved. global_opt: - log-posterior value of the global optimum + Log-posterior value of the global optimum. order: - Specifies the precise algorithm for extrapolation: can be ``0`` ( - just one parameter is updated), ``1`` (last two points used to - extrapolate all parameters), and ``np.nan`` (indicates that a more - complex regression should be used) + Specifies the precise algorithm for extrapolation. + Available options are: + + * ``0``: just one parameter is updated + * ``1``: the last two points are used to extrapolate all parameters + * ``np.nan``: indicates that a more complex regression should be used + as determined by :attr:`pypesto.profile.ProfileOptions.reg_order`. + Returns ------- - x_new: - The updated parameter vector, of size `dim_full`. + The updated parameter vector, of size `dim_full`. """ # restrict step proposal to minimum and maximum step size def clip_to_minmax(step_size_proposal): - return clip( + return np.clip( step_size_proposal, options.min_step_size, options.max_step_size ) # restrict step proposal to bounds def clip_to_bounds(step_proposal): - return clip(step_proposal, problem.lb_full, problem.ub_full) + return np.clip(step_proposal, problem.lb_full, problem.ub_full) - # check if this is the first step - n_profile_points = len(current_profile.fval_path) problem.fix_parameters(par_index, x[par_index]) # Get update directions and first step size guesses @@ -197,7 +197,6 @@ def clip_to_bounds(step_proposal): x, par_index, par_direction, - n_profile_points, global_opt, order, current_profile, @@ -208,43 +207,47 @@ def clip_to_bounds(step_proposal): # check whether we must make a minimum step anyway, since we're close to # the next bound min_delta_x = x[par_index] + par_direction * options.min_step_size + if par_direction == -1 and (min_delta_x < problem.lb_full[par_index]): step_length = problem.lb_full[par_index] - x[par_index] return x + step_length * delta_x_dir - elif par_direction == 1 and (min_delta_x > problem.ub_full[par_index]): + + if par_direction == 1 and (min_delta_x > problem.ub_full[par_index]): step_length = problem.ub_full[par_index] - x[par_index] return x + step_length * delta_x_dir # parameter extrapolation function - def par_extrapol(step_length): - # Do we have enough points to do a regression? - if np.isnan(order) and n_profile_points > 2: - x_step_tmp = [] + n_profile_points = len(current_profile.fval_path) + + # Do we have enough points to do a regression? + if np.isnan(order) and n_profile_points > 2: + + def par_extrapol(step_length): + x_step = [] # loop over parameters, extrapolate each one for i_par in range(problem.dim_full): if i_par == par_index: # if we meet the profiling parameter, just increase, # don't extrapolate - x_step_tmp.append( - x[par_index] + step_length * par_direction - ) + x_step.append(x[par_index] + step_length * par_direction) elif i_par in problem.x_fixed_indices: # common fixed parameter: will be ignored anyway later - x_step_tmp.append(np.nan) + x_step.append(np.nan) else: # extrapolate cur_par_extrapol = np.poly1d(reg_par[i_par]) - x_step_tmp.append( + x_step.append( cur_par_extrapol( x[par_index] + step_length * par_direction ) ) - x_step = np.array(x_step_tmp) - else: - # if we do simple extrapolation - x_step = x + step_length * delta_x_dir + return clip_to_bounds(x_step) - return clip_to_bounds(x_step) + else: + # if not, we do simple extrapolation + def par_extrapol(step_length): + x_step = x + step_length * delta_x_dir + return clip_to_bounds(x_step) # compute proposal next_x = par_extrapol(step_size_guess) @@ -265,7 +268,6 @@ def par_extrapol(step_length): return do_line_search( next_x, step_size_guess, - "decrease" if next_obj_target < next_obj else "increase", par_extrapol, next_obj, next_obj_target, @@ -280,19 +282,31 @@ def par_extrapol(step_length): def handle_profile_history( x: np.ndarray, par_index: int, - par_direction: int, - n_profile_points: int, + par_direction: Literal[1, -1], global_opt: float, order: int, current_profile: ProfilerResult, problem: Problem, options: ProfileOptions, -) -> Tuple: +) -> tuple[float, np.array, list[float], float]: """Compute the very first step direction update guesses. Check whether enough steps have been taken for applying regression, computes regression or simple extrapolation. + + Returns + ------- + step_size_guess: + Guess for the step size. + delta_x_dir: + Parameter update direction. + reg_par: + The regression polynomial for profile extrapolation. + delta_obj_value: + The difference of the objective function value between the last point and `global_opt`. """ + n_profile_points = len(current_profile.fval_path) + # set the update direction delta_x_dir = np.zeros(len(x)) delta_x_dir[par_index] = par_direction @@ -320,29 +334,28 @@ def handle_profile_history( delta_x_dir = last_delta_x / step_size_guess elif np.isnan(order): # compute the regression polynomial for parameter extrapolation - reg_par = get_reg_polynomial( - n_profile_points, par_index, current_profile, problem, options + par_index, current_profile, problem, options ) return step_size_guess, delta_x_dir, reg_par, delta_obj_value def get_reg_polynomial( - n_profile_points: int, par_index: int, current_profile: ProfilerResult, problem: Problem, options: ProfileOptions, -) -> List[float]: +) -> list[float]: """Compute the regression polynomial. - Used to step proposal extrapolation from the last profile points + Used to step proposal extrapolation from the last profile points. """ # determine interpolation order + n_profile_points = len(current_profile.fval_path) reg_max_order = np.floor(n_profile_points / 2) - reg_order = np.min([reg_max_order, options.reg_order]) - reg_points = np.min([n_profile_points, options.reg_points]) + reg_order = min(reg_max_order, options.reg_order) + reg_points = min(n_profile_points, options.reg_points) # set up matrix of regression parameters reg_par = [] @@ -380,7 +393,6 @@ def get_reg_polynomial( def do_line_search( next_x: np.ndarray, step_size_guess: float, - direction: Literal['increase', 'decrease'], par_extrapol: Callable, next_obj: float, next_obj_target: float, @@ -394,36 +406,64 @@ def do_line_search( Based on the objective function we want to reach, based on the current position in parameter space and on the first guess for the proposal. + + Parameters + ---------- + next_x: + Starting parameters for the line search. + step_size_guess: + First guess for the step size. + par_extrapol: + Parameter extrapolation function. + next_obj: + Objective function value at `next_x`. + next_obj_target: + Objective function value we want to reach. + clip_to_minmax: + Function to clip the step size to minimum and maximum step size. + clip_to_bounds: + Function to clip the parameters to the bounds. + par_index: + Index of the parameter we are profiling. + problem: + The parameter estimation problem. + options: + Profile likelihood options. + + Returns + ------- + Parameter vector that is expected to yield the objective function value + closest to `next_obj_target`. """ # Was the initial step too big or too small? + direction = "decrease" if next_obj_target < next_obj else "increase" if direction == 'increase': adapt_factor = options.step_size_factor else: adapt_factor = 1 / options.step_size_factor # Loop until correct step size was found - stop_search = False - while not stop_search: + while True: # Adapt step size of guess - last_x = copy.copy(next_x) + last_x = next_x step_size_guess = clip_to_minmax(step_size_guess * adapt_factor) next_x = clip_to_bounds(par_extrapol(step_size_guess)) # Check if we hit the bounds - hit_bounds = ( + if ( direction == 'decrease' and step_size_guess == options.min_step_size - ) or ( + ): + return next_x + if ( direction == 'increase' and step_size_guess == options.max_step_size - ) - - if hit_bounds: + ): return next_x # compute new objective value problem.fix_parameters(par_index, next_x[par_index]) - last_obj = copy.copy(next_obj) + last_obj = next_obj next_obj = problem.objective(problem.get_reduced_vector(next_x)) # check for root crossing and compute correct step size in case @@ -448,22 +488,3 @@ def next_x_interpolate( # fix final guess and return return last_x + add_x - - -def clip( - vector_guess: Union[float, np.ndarray], - lower: Union[float, np.ndarray], - upper: Union[float, np.ndarray], -) -> Union[float, np.ndarray]: - """Restrict a scalar or a vector to given bounds. - - ``vector_guess`` is modified in-place if it is an array. - """ - if isinstance(vector_guess, float): - return np.max([np.min([vector_guess, upper]), lower]) - - for i_par, i_guess in enumerate(vector_guess): - vector_guess[i_par] = np.max( - [np.min([i_guess, upper[i_par]]), lower[i_par]] - ) - return vector_guess diff --git a/pypesto/profile/task.py b/pypesto/profile/task.py index 4c4ea08f0..0cdbfaddc 100644 --- a/pypesto/profile/task.py +++ b/pypesto/profile/task.py @@ -1,5 +1,5 @@ import logging -from typing import Callable +from typing import Any, Callable import pypesto.optimize @@ -55,7 +55,7 @@ def __init__( self.i_par = i_par self.options = options - def execute(self) -> 'pypesto.profile.ProfilerResult': + def execute(self) -> dict[str, Any]: """Compute profile in descending and ascending direction.""" logger.debug(f"Executing task {self.i_par}.") diff --git a/pypesto/profile/util.py b/pypesto/profile/util.py index 350d64459..ceac7ffde 100644 --- a/pypesto/profile/util.py +++ b/pypesto/profile/util.py @@ -1,5 +1,5 @@ """Utility function for profile module.""" -from typing import Any, Dict, Iterable, Tuple +from typing import Any, Iterable import numpy as np import scipy.stats @@ -25,8 +25,7 @@ def chi2_quantile_to_ratio(alpha: float = 0.95, df: int = 1): Returns ------- - ratio: - Corresponds to a likelihood ratio. + The computed likelihood ratio threshold. """ quantile = scipy.stats.chi2.ppf(alpha, df=df) ratio = np.exp(-quantile / 2) @@ -35,11 +34,11 @@ def chi2_quantile_to_ratio(alpha: float = 0.95, df: int = 1): def calculate_approximate_ci( xs: np.ndarray, ratios: np.ndarray, confidence_ratio: float -) -> Tuple[float, float]: +) -> tuple[float, float]: """ Calculate approximate confidence interval based on profile. - Interval bounds are linerly interpolated. + Interval bounds are linearly interpolated. Parameters ---------- @@ -50,12 +49,11 @@ def calculate_approximate_ci( The likelihood ratios corresponding to the parameter values. confidence_ratio: Minimum confidence ratio to base the confidence interval upon, as - obtained via `pypesto.profile.chi2_quantile_to_ratio`. + obtained via :func:`pypesto.profile.chi2_quantile_to_ratio`. Returns ------- - lb, ub: - Bounds of the approximate confidence interval. + Bounds of the approximate confidence interval. """ # extract indices where the ratio is larger than the minimum ratio (indices,) = np.where(ratios >= confidence_ratio) @@ -147,7 +145,7 @@ def initialize_profile( def fill_profile_list( profile_result: ProfileResult, - optimizer_result: Dict[str, Any], + optimizer_result: dict[str, Any], profile_index: Iterable[int], profile_list: int, problem_dimension: int, diff --git a/pypesto/profile/validation_intervals.py b/pypesto/profile/validation_intervals.py index b77bcbd3d..aff7f6aa3 100644 --- a/pypesto/profile/validation_intervals.py +++ b/pypesto/profile/validation_intervals.py @@ -33,8 +33,8 @@ def validation_profile_significance( The reasoning behind their approach is, that a validation data set is outside the validation interval, if fitting the full data set - would lead to a fit $\theta_{new}$, that does not contain the old - fit $\theta_{train}$ in their (Profile-Likelihood) based + would lead to a fit :math:`\theta_{new}`, that does not contain the old + fit :math:`\theta_{train}` in their (Profile-Likelihood) based parameter-confidence intervals. (I.e. the old fit would be rejected by the fit of the full data.) @@ -50,34 +50,28 @@ def validation_profile_significance( problem_full_data: pypesto.problem, such that the objective is the negative-log-likelihood of the training and validation data set. - result_training_data: - result object from the fitting of the training data set only. - + Result object from the fitting of the training data set only. result_full_data - pypesto.result object that contains the result of fitting + Result object that contains the result of fitting training and validation data combined. - n_starts number of starts for fitting the full data set - (if result_full_data is not provided). - + (if `result_full_data` is not provided). optimizer: - optimizer used for refitting the data (if result_full_data is not + optimizer used for refitting the data (if `result_full_data` is not provided). - - engine - engine for refitting (if result_full_data is not provided). - + engine: + engine for refitting (if `result_full_data` is not provided). lsq_objective: - indicates if the objective of problem_full_data corresponds to a nllh - (False), or a chi^2 value (True). + indicates if the objective of `problem_full_data` corresponds to a nllh + (``False``), or a :math:`\chi^2` value (``True``). return_significance: - indicates, if the function should return the significance (True) (i.e. + indicates, if the function should return the significance (``True``) (i.e. the probability, that the new data set lies outside the Confidence Interval for the validation experiment, as given by the method), or the largest alpha, such that the validation experiment still lies - within the Confidence Interval (False). I.e. alpha = 1-significance. + within the Confidence Interval (``False``). I.e. :math:`\alpha = 1-significance`. .. [#Kreutz] Kreutz, Clemens, Raue, Andreas and Timmer, Jens. diff --git a/pypesto/profile/walk_along_profile.py b/pypesto/profile/walk_along_profile.py index 14ae47a3d..202019ea7 100644 --- a/pypesto/profile/walk_along_profile.py +++ b/pypesto/profile/walk_along_profile.py @@ -1,5 +1,5 @@ import logging -from typing import Callable +from typing import Callable, Literal import numpy as np @@ -15,7 +15,7 @@ def walk_along_profile( current_profile: ProfilerResult, problem: Problem, - par_direction: int, + par_direction: Literal[1, -1], optimizer: Optimizer, options: ProfileOptions, create_next_guess: Callable, @@ -54,11 +54,10 @@ def walk_along_profile( Returns ------- - current_profile: - The current profile, modified in-place. + The current profile, modified in-place. """ - # create variables which are needed during iteration - stop_profile = False + if par_direction not in (-1, 1): + raise AssertionError("par_direction must be -1 or 1") # while loop for profiling (will be exited by break command) while True: @@ -67,18 +66,16 @@ def walk_along_profile( # check if the next profile point needs to be computed # ... check bounds - if par_direction == -1: - stop_profile = x_now[i_par] <= problem.lb_full[[i_par]] - elif par_direction == 1: - stop_profile = x_now[i_par] >= problem.ub_full[[i_par]] - else: - raise AssertionError("par_direction must be -1 or 1") + if par_direction == -1 and x_now[i_par] <= problem.lb_full[[i_par]]: + break + if par_direction == 1 and x_now[i_par] >= problem.ub_full[[i_par]]: + break # ... check likelihood ratio - if not options.whole_path: - stop_profile |= current_profile.ratio_path[-1] < options.ratio_min - - if stop_profile: + if ( + not options.whole_path + and current_profile.ratio_path[-1] < options.ratio_min + ): break # compute the new start point for optimization @@ -92,10 +89,9 @@ def walk_along_profile( global_opt, ) - # fix current profiling parameter to current value and set - # start point + # fix current profiling parameter to current value and set start point problem.fix_parameters(i_par, x_next[i_par]) - startpoint = np.array([x_next[i] for i in problem.x_free_indices]) + startpoint = x_next[problem.x_free_indices] # run optimization if startpoint.size > 0: @@ -113,12 +109,8 @@ def walk_along_profile( if np.isfinite(optimizer_result.fval): break - profiled_par_id = problem.x_names[i_par] - profiled_par_value = startpoint[ - problem.x_free_indices.index(i_par) - ] logger.warning( - f"Optimization at {profiled_par_id}={profiled_par_value} failed." + f"Optimization at {problem.x_names[i_par]}={x_next[i_par]} failed." ) # sample a new starting point for another attempt # might be preferable to stay close to the previous point, at least initially,