diff --git a/pypesto/objective/amici/amici.py b/pypesto/objective/amici/amici.py index 42ce9269c..1693e948e 100644 --- a/pypesto/objective/amici/amici.py +++ b/pypesto/objective/amici/amici.py @@ -10,6 +10,7 @@ from typing import TYPE_CHECKING, Union import numpy as np +import pandas as pd from ...C import ( FVAL, @@ -36,6 +37,8 @@ ) if TYPE_CHECKING: + import amici.importers.petab + from ...hierarchical import InnerCalculatorCollector try: @@ -718,3 +721,67 @@ def update_from_problem( ) in condition_mapping.map_preeq_fix.items(): if (val := id_to_val.get(mapped_to_par)) is not None: condition_mapping.map_preeq_fix[model_par] = val + + +class AmiciPetabV2Objective(AmiciObjective): + """An AMICI objective constructed from a PEtab v2 problem.""" + + def __init__( + self, + petab_importer: amici.importers.petab.PetabImporter, + **kwargs, + ) -> None: + from .amici_calculator import AmiciCalculatorPetabV2 + + self._petab_simulator: amici.petab.petab_importer.PetabSimulator = ( + petab_importer.create_simulator() + ) + self.petab_problem = petab_importer.petab_problem + amici_model = self._petab_simulator.model + amici_solver = self._petab_simulator.solver + edatas = self._petab_simulator.exp_man.create_edatas() + + super().__init__( + amici_model=amici_model, + amici_solver=amici_solver, + edatas=edatas, + calculator=AmiciCalculatorPetabV2(self._petab_simulator), + **kwargs, + ) + + def __deepcopy__(self, memo=None): + """Override AmiciObjective.__deepcopy__.""" + if memo is None: + memo = {} + cls = self.__class__ + result = cls.__new__(cls) + memo[id(self)] = result + for k, v in self.__dict__.items(): + setattr(result, k, copy.deepcopy(v, memo)) + return result + + def __getstate__(self) -> dict: + """Use Python's default pickling semantics (shallow copy of instance dict).""" + return dict(self.__dict__) + + def __setstate__(self, state: dict) -> None: + """Restore state using the instance dict (default unpickling behaviour).""" + self.__dict__.update(state) + + def rdatas_to_simulation_df( + self, + rdatas: Sequence[amici.ReturnData], + ) -> pd.DataFrame: + """ + See :meth:`rdatas_to_measurement_df`. + + Except a petab simulation dataframe is created, i.e. the measurement + column label is adjusted. + """ + from amici.importers.petab import rdatas_to_simulation_df + + return rdatas_to_simulation_df( + rdatas, + self._petab_simulator._model, + self._petab_simulator._petab_problem, + ) diff --git a/pypesto/objective/amici/amici_calculator.py b/pypesto/objective/amici/amici_calculator.py index f1ab22c3f..03ecf1126 100644 --- a/pypesto/objective/amici/amici_calculator.py +++ b/pypesto/objective/amici/amici_calculator.py @@ -154,6 +154,163 @@ def __call__( ) +class AmiciCalculatorPetabV2(AmiciCalculator): + """Class to perform the AMICI call and obtain objective function values.""" + + def __init__( + self, + petab_simulator: amici.petab.petab_importer.PetabSimulator, + **kwargs, + ): + super().__init__(**kwargs) + self.petab_simulator = petab_simulator + + def __call__( + self, + x_dct: dict, + sensi_orders: tuple[int], + mode: ModeType, + amici_model: AmiciModel, + amici_solver: AmiciSolver, + edatas: list[amici.ExpData], + n_threads: int, + x_ids: Sequence[str], + parameter_mapping: ParameterMapping, + fim_for_hess: bool, + ): + """Perform the actual AMICI call. + + Called within the :func:`AmiciObjective.__call__` method. + + Parameters + ---------- + x_dct: + Parameters for which to compute function value and derivatives. + sensi_orders: + Tuple of requested sensitivity orders. + mode: + Call mode (function value or residual based). + amici_model: + The AMICI model. + amici_solver: + The AMICI solver. + edatas: + The experimental data. + n_threads: + Number of threads for AMICI call. + x_ids: + Ids of optimization parameters. + parameter_mapping: + Mapping of optimization to simulation parameters. + fim_for_hess: + Whether to use the FIM (if available) instead of the Hessian (if + requested). + """ + amici_solver = self.petab_simulator._solver + + if mode != MODE_FUN: + raise NotImplementedError( + "Only function value mode is currently supported for " + f"PEtab v2. Got mode {mode}." + ) + + # TODO: -> method + # set order in solver + sensi_order = 0 + if sensi_orders: + sensi_order = max(sensi_orders) + + if sensi_order == 2 and fim_for_hess: + # we use the FIM + amici_solver.set_sensitivity_order(sensi_order - 1) + else: + amici_solver.set_sensitivity_order(sensi_order) + + dim = len(x_ids) + + # run amici simulation + result = self.petab_simulator.simulate(x_dct) + rdatas = result.rdatas + + # check if the simulation failed + if any(rdata["status"] < 0.0 for rdata in rdatas): + return get_error_output( + amici_model, edatas, rdatas, sensi_orders, mode, dim + ) + + nllh, snllh, s2nllh, chi2, res, sres = init_return_values( + sensi_orders, mode, dim + ) + nllh = -result.llh + + if ( + not self._known_least_squares_safe + and mode == MODE_RES + and 1 in sensi_orders + ): + if not amici_model.get_add_sigma_residuals() and any( + ( + (r["ssigmay"] is not None and np.any(r["ssigmay"])) + or (r["ssigmaz"] is not None and np.any(r["ssigmaz"])) + ) + for r in rdatas + ): + raise RuntimeError( + "Cannot use least squares solver with" + "parameter dependent sigma! Support can be " + "enabled via " + "amici_model.setAddSigmaResiduals()." + ) + self._known_least_squares_safe = True # don't check this again + + # TODO: compute res, sres + + if 1 in sensi_orders: + if result.sllh is None and np.isnan(result["llh"]): + # TODO: to amici -- set sllh even if llh is nan? + snllh = np.full(len(x_ids), np.nan) + else: + try: + # llh to nllh, dict to array + snllh = -np.array( + [ + result.sllh[ + x_id + ] # if x_id in res["sllh"] else 0.0 + for x_id in x_ids + if x_id in x_dct.keys() + ] + ) + except KeyError as e: + # A requested sensitivity is missing. + # Probably the affected parameter is a fixed parameter + # in amici instead of a sensitivity parameter + # (non_estimated_parameters_as_constants=True ?). + # In this case, only max(sensi_orders) == 0 is supported + # unless this parameter is fixed in the pypesto problem. + raise ValueError( + f"Cannot compute gradient, missing entry for " + f"{set(result.sllh) - set(x_dct.keys())}." + ) from e + if 2 in sensi_orders: + if result.s2llh is None and np.isnan(result.llh): + # TODO: to amici -- set s2llh even if llh is nan? + s2nllh = np.full((len(x_ids), len(x_ids)), np.nan) + else: + s2nllh = -result.s2llh + + ret = { + FVAL: nllh, + GRAD: snllh, + HESS: s2nllh, + RES: res, + SRES: sres, + RDATAS: rdatas, + } + + return filter_return_dict(ret) + + def calculate_function_values( rdatas, sensi_orders: tuple[int, ...], diff --git a/pypesto/petab/importer.py b/pypesto/petab/importer.py index eaaca9338..57b121233 100644 --- a/pypesto/petab/importer.py +++ b/pypesto/petab/importer.py @@ -20,6 +20,8 @@ except ImportError: roadrunner = None +from petab import v2 + from ..C import ( AMICI, CENSORED, @@ -39,6 +41,7 @@ from ..startpoint import StartpointMethod from .objective_creator import ( AmiciObjectiveCreator, + AmiciPetabV2ObjectiveCreator, ObjectiveCreator, PetabSimulatorObjectiveCreator, RoadRunnerObjectiveCreator, @@ -69,7 +72,7 @@ class PetabImporter: def __init__( self, - petab_problem: petab.Problem, + petab_problem: petab.Problem | v2.Problem, output_folder: str | None = None, model_name: str | None = None, validate_petab: bool = True, @@ -116,6 +119,12 @@ def __init__( self.petab_problem = petab_problem self._hierarchical = hierarchical + if self._hierarchical and isinstance(self.petab_problem, v2.Problem): + raise NotImplementedError( + "Hierarchical optimization is not yet supported " + "for PEtab v2 problems." + ) + self._non_quantitative_data_types = ( get_petab_non_quantitative_data_types(petab_problem) ) @@ -148,8 +157,19 @@ def __init__( self.validate_petab = validate_petab if self.validate_petab: - if petab.lint_problem(petab_problem): + if isinstance(petab_problem, petab.Problem) and petab.lint_problem( + petab_problem + ): raise ValueError("Invalid PEtab problem.") + if ( + isinstance(petab_problem, v2.Problem) + and ( + validation_result := petab_problem.validate() + ).has_errors() + ): + validation_result.log(logger=logger) + raise ValueError("Invalid PEtab v2 problem.") + if self._hierarchical and validate_petab_hierarchical: from ..hierarchical.petab import ( validate_hierarchical_petab_problem, @@ -185,7 +205,15 @@ def from_yaml( simulator_type: str = AMICI, ) -> PetabImporter: """Simplified constructor using a petab yaml file.""" - petab_problem = petab.Problem.from_yaml(yaml_config) + from petab.versions import get_major_version + + match get_major_version(yaml_config): + case 1: + petab_problem = petab.Problem.from_yaml(yaml_config) + case 2: + petab_problem = v2.Problem.from_yaml(yaml_config) + case _: + raise ValueError("Only PEtab v1 and v2 are supported.") return PetabImporter( petab_problem=petab_problem, @@ -277,6 +305,42 @@ def create_prior(self) -> NegLogParameterPriors | None: else: return None + def _create_prior_v2(self) -> NegLogParameterPriors | None: + """Create prior for PEtab v2 problem.""" + import petab.v2.C as petab_c + + import pypesto.C as pypesto_c + + petab_to_pypesto = { + petab_c.LAPLACE: pypesto_c.LAPLACE, + petab_c.LOG_LAPLACE: pypesto_c.LOG_LAPLACE, + petab_c.LOG_NORMAL: pypesto_c.LOG_NORMAL, + petab_c.LOG_UNIFORM: pypesto_c.LOG_UNIFORM, + petab_c.NORMAL: pypesto_c.NORMAL, + petab_c.UNIFORM: pypesto_c.UNIFORM, + } + + prior_list = [] + for parameter in self.petab_problem.parameters: + if not parameter.estimate or parameter.prior_distribution is None: + continue + + prior_list.append( + get_parameter_prior_dict( + index=len(prior_list), + # TODO map names + prior_type=petab_to_pypesto.get( + str(parameter.prior_distribution), + str(parameter.prior_distribution), + ), + prior_parameters=parameter.prior_parameters, + parameter_scale="lin", + ) + ) + if prior_list: + return NegLogParameterPriors(prior_list) + return None + def create_startpoint_method(self, **kwargs) -> StartpointMethod: """Create a startpoint method. @@ -306,6 +370,22 @@ def create_objective_creator( has to be provided. Otherwise the argument is not used. """ + if isinstance(self.petab_problem, v2.Problem): + if simulator_type != AMICI: + raise ValueError( + "Only 'amici' simulator type is supported for PEtab v2 " + "problems." + ) + return AmiciPetabV2ObjectiveCreator( + petab_problem=self.petab_problem, + output_folder=self.output_folder, + model_name=self.model_name, + hierarchical=self._hierarchical, + inner_options=self.inner_options, + non_quantitative_data_types=self._non_quantitative_data_types, + validate_petab=self.validate_petab, + ) + if simulator_type == AMICI: return AmiciObjectiveCreator( petab_problem=self.petab_problem, @@ -364,10 +444,19 @@ def create_problem( objective = self.objective_constructor.create_objective(**kwargs) x_fixed_indices = self.petab_problem.x_fixed_indices - x_fixed_vals = self.petab_problem.x_nominal_fixed_scaled x_ids = self.petab_problem.x_ids - lb = self.petab_problem.lb_scaled - ub = self.petab_problem.ub_scaled + + if isinstance(self.petab_problem, petab.Problem): + # PEtab v1 + x_fixed_vals = self.petab_problem.x_nominal_fixed_scaled + lb = self.petab_problem.lb_scaled + ub = self.petab_problem.ub_scaled + prior = self.create_prior() + else: + x_fixed_vals = self.petab_problem.x_nominal_fixed + lb = self.petab_problem.lb + ub = self.petab_problem.ub + prior = self._create_prior_v2() # Raise error if the correct calculator is not used. if self._hierarchical: @@ -395,10 +484,18 @@ def create_problem( map(x_ids.index, self.petab_problem.x_fixed_ids) ) - x_scales = [ - self.petab_problem.parameter_df.loc[x_id, petab.PARAMETER_SCALE] - for x_id in x_ids - ] + # parameter scales -- only for outer parameters + if isinstance(self.petab_problem, petab.Problem): + # PEtab v1 + x_scales = [ + self.petab_problem.parameter_df.loc[ + x_id, petab.PARAMETER_SCALE + ] + for x_id in x_ids + ] + else: + # PEtab v2 -- no parameter scaling + x_scales = [petab.LIN for x_id in x_ids] if problem_kwargs is None: problem_kwargs = {} @@ -406,8 +503,6 @@ def create_problem( if startpoint_kwargs is None: startpoint_kwargs = {} - prior = self.create_prior() - if prior is not None: if self._hierarchical: raise NotImplementedError( diff --git a/pypesto/petab/objective_creator.py b/pypesto/petab/objective_creator.py index 261a010fb..d406894f2 100644 --- a/pypesto/petab/objective_creator.py +++ b/pypesto/petab/objective_creator.py @@ -18,6 +18,7 @@ import numpy as np import pandas as pd import petab.v1 as petab +from petab import v2 from petab.v1.C import ( OBSERVABLE_FORMULA, PREEQUILIBRATION_CONDITION_ID, @@ -606,6 +607,453 @@ def prediction_to_petab_simulation_df( ).rename(columns={petab.MEASUREMENT: petab.SIMULATION}) +class AmiciPetabV2ObjectiveCreator(ObjectiveCreator, AmiciObjectBuilder): + """ObjectiveCreator for creating an amici objective function.""" + + def __init__( + self, + petab_problem: v2.Problem, + hierarchical: bool = False, + non_quantitative_data_types: Iterable[str] | None = None, + inner_options: dict[str, Any] | None = None, + output_folder: str | None = None, + model_name: str | None = None, + validate_petab: bool = True, + ): + """ + Initialize the creator. + + Parameters + ---------- + petab_problem: + The PEtab problem. + hierarchical: + Whether to use hierarchical optimization. + non_quantitative_data_types: + The non-quantitative data types to consider. + inner_options: + Options for the inner optimization. + output_folder: + The output folder for the compiled model. + model_name: + The name of the AMICI model module. + validate_petab: + Whether to check the PEtab problem for errors. + """ + if hierarchical or non_quantitative_data_types or inner_options: + raise NotImplementedError() + + self.petab_problem = petab_problem + self._hierarchical = hierarchical + self._non_quantitative_data_types = non_quantitative_data_types + self.inner_options = inner_options + self.output_folder = output_folder + self.model_name = model_name + self.validate_petab = validate_petab + + def create_model( + self, + force_compile: bool = False, + verbose: bool = True, + **kwargs, + ) -> amici.Model: + """ + Import amici model. + + Parameters + ---------- + force_compile: + If False, the model is compiled only if the output folder does not + exist yet. If True, the output folder is deleted and the model + (re-)compiled in either case. + + .. warning:: + If `force_compile`, then an existing folder of that name will + be deleted. + verbose: + Passed to AMICI's model compilation. If True, the compilation + progress is printed. + kwargs: + Extra arguments passed to amici.SbmlImporter.sbml2amici + """ + # courtesy check whether target is folder + if os.path.exists(self.output_folder) and not os.path.isdir( + self.output_folder + ): + raise AssertionError( + f"Refusing to remove {self.output_folder} for model " + f"compilation: Not a folder." + ) + + # compile + if self._must_compile(force_compile): + logger.info( + f"Compiling amici model to folder {self.output_folder}." + ) + if self.petab_problem.model.type_id == MODEL_TYPE_SBML: + self.compile_model( + validate=self.validate_petab, + verbose=verbose, + **kwargs, + ) + else: + self.compile_model(verbose=verbose, **kwargs) + else: + logger.debug( + f"Using existing amici model in folder {self.output_folder}." + ) + + return self._create_model() + + def _create_model(self) -> amici.Model: + """Load model module and return the model, no checks/compilation.""" + # load moduĺe + module = amici.import_model_module( + module_name=self.model_name, module_path=self.output_folder + ) + model = module.get_model() + + from amici.importers.petab.v1._import_helpers import check_model + + check_model( + amici_model=model, + petab_problem=self.petab_problem, + ) + + return model + + def _must_compile(self, force_compile: bool): + """Check whether the model needs to be compiled first.""" + # asked by user + if force_compile: + return True + + # folder does not exist + if not os.path.exists(self.output_folder) or not os.listdir( + self.output_folder + ): + return True + + # try to import (in particular checks version) + try: + # importing will already raise an exception if version wrong + amici.import_model_module(self.model_name, self.output_folder) + except ModuleNotFoundError: + return True + except amici.AmiciVersionError as e: + logger.info( + f"amici model will be re-imported due to version mismatch: {e}" + ) + return True + + # no need to (re-)compile + return False + + def compile_model(self, **kwargs): + """ + Compile the model. + + If the output folder exists already, it is first deleted. + + Parameters + ---------- + kwargs: + Extra arguments passed to :meth:`amici.sbml_import.SbmlImporter.sbml2amici` + or :func:`amici.pysb_import.pysb2amici`. + """ + # delete output directory + if os.path.exists(self.output_folder): + shutil.rmtree(self.output_folder) + + petab_importer = self._create_amici_importer() + petab_importer.import_module(force_import=True) + + def _create_amici_importer( + self, + ) -> amici.petab.petab_importer.PetabImporter: + """Create an AMICI PEtab importer.""" + from amici.importers.petab import PetabImporter + + petab_importer = PetabImporter( + self.petab_problem, + output_dir=self.output_folder, + module_name=self.model_name, + ) + return petab_importer + + def create_solver( + self, + model: amici.Model = None, + verbose: bool = True, + ) -> amici.Solver: + """Return model solver.""" + # create model + if model is None: + model = self.create_model(verbose=verbose) + + solver = model.create_solver() + return solver + + def create_edatas( + self, + model: amici.Model = None, + simulation_conditions=None, + verbose: bool = True, + ) -> list[amici.ExpData]: + """Create list of :class:`amici.amici.ExpData` objects.""" + # create model + # if model is None: + # model = self.create_model(verbose=verbose) + # + # return amici.petab.conditions.create_edatas( + # amici_model=model, + # petab_problem=self.petab_problem, + # simulation_conditions=simulation_conditions, + # ) + raise NotImplementedError() + + def create_objective( + self, + model: amici.Model = None, + solver: amici.Solver = None, + edatas: Sequence[amici.ExpData] = None, + force_compile: bool = False, + verbose: bool = True, + **kwargs, + ) -> AmiciObjective: + """Create a :class:`pypesto.objective.AmiciObjective`. + + Parameters + ---------- + model: + The AMICI model. + solver: + The AMICI solver. + edatas: + The experimental data in AMICI format. + force_compile: + Whether to force-compile the model if not passed. + verbose: + Passed to AMICI's model compilation. If True, the compilation + progress is printed. + **kwargs: + Additional arguments passed on to the objective. In case of ordinal + or semiquantitative measurements, ``inner_options`` can optionally + be passed here. If none are given, ``inner_options`` given to the + importer constructor (or inner defaults) will be chosen. + + Returns + ------- + A :class:`pypesto.objective.AmiciObjective` for the model and the data. + """ + from ..objective.amici.amici import AmiciPetabV2Objective + + petab_importer = self._create_amici_importer() + + return AmiciPetabV2Objective( + petab_importer=petab_importer, + amici_object_builder=self, + x_ids=self.petab_problem.x_ids, + ) + + # TODO + + simulation_conditions = petab.get_simulation_conditions( + self.petab_problem.measurement_df + ) + if model is None: + model = self.create_model( + force_compile=force_compile, verbose=verbose + ) + if solver is None: + solver = self.create_solver(model) + # create conditions and edatas from measurement data + if edatas is None: + edatas = self.create_edatas( + model=model, simulation_conditions=simulation_conditions + ) + parameter_mapping = ( + amici.petab.parameter_mapping.create_parameter_mapping( + petab_problem=self.petab_problem, + simulation_conditions=simulation_conditions, + scaled_parameters=True, + amici_model=model, + fill_fixed_parameters=False, + ) + ) + par_ids = self.petab_problem.x_ids + + # fill in dummy parameters (this is needed since some objective + # initialization e.g. checks for preeq parameters) + problem_parameters = dict( + zip( + self.petab_problem.x_ids, + self.petab_problem.x_nominal_scaled, + strict=True, + ) + ) + amici.petab.conditions.fill_in_parameters( + edatas=edatas, + problem_parameters=problem_parameters, + scaled_parameters=True, + parameter_mapping=parameter_mapping, + amici_model=model, + ) + + calculator = None + amici_reporting = None + + if ( + self._non_quantitative_data_types is not None + and self._hierarchical + ): + inner_options = kwargs.pop("inner_options", None) + inner_options = ( + inner_options + if inner_options is not None + else self.inner_options + ) + calculator = InnerCalculatorCollector( + self._non_quantitative_data_types, + self.petab_problem, + model, + edatas, + inner_options, + ) + amici_reporting = amici.RDataReporting.full + + # FIXME: currently not supported with hierarchical + if "guess_steadystate" in kwargs and kwargs["guess_steadystate"]: + warnings.warn( + "`guess_steadystate` not supported with hierarchical " + "optimization. Disabling `guess_steadystate`.", + stacklevel=1, + ) + kwargs["guess_steadystate"] = False + inner_parameter_ids = calculator.get_inner_par_ids() + par_ids = [x for x in par_ids if x not in inner_parameter_ids] + + max_sensi_order = kwargs.get("max_sensi_order", None) + + if ( + self._non_quantitative_data_types is not None + and any( + data_type in self._non_quantitative_data_types + for data_type in [ORDINAL, CENSORED, SEMIQUANTITATIVE] + ) + and max_sensi_order is not None + and max_sensi_order > 1 + ): + raise ValueError( + "Ordinal, censored and semiquantitative data cannot be " + "used with second order sensitivities. Use a up to first order " + "method or disable ordinal, censored and semiquantitative " + ) + + # create objective + obj = AmiciObjective( + amici_model=model, + amici_solver=solver, + edatas=edatas, + x_ids=par_ids, + x_names=par_ids, + parameter_mapping=parameter_mapping, + amici_object_builder=self, + calculator=calculator, + amici_reporting=amici_reporting, + **kwargs, + ) + + return obj + + def create_predictor( + self, + objective: AmiciObjective = None, + amici_output_fields: Sequence[str] = None, + post_processor: Callable | None = None, + post_processor_sensi: Callable | None = None, + post_processor_time: Callable | None = None, + max_chunk_size: int | None = None, + output_ids: Sequence[str] = None, + condition_ids: Sequence[str] = None, + ) -> AmiciPredictor: + """Create a :class:`pypesto.predict.AmiciPredictor`.""" + raise NotImplementedError() + + def rdatas_to_measurement_df( + self, + rdatas: Sequence[amici.ReturnData], + model: amici.Model = None, + verbose: bool = True, + ) -> pd.DataFrame: + """ + Create a measurement dataframe in the petab format. + + Parameters + ---------- + rdatas: + A list of rdatas as produced by + ``pypesto.AmiciObjective.__call__(x, return_dict=True)['rdatas']``. + model: + The amici model. + verbose: + Passed to AMICI's model compilation. If True, the compilation + progress is printed. + + Returns + ------- + A dataframe built from the rdatas in the format as in + ``self.petab_problem.measurement_df``. + """ + # create model + if model is None: + model = self.create_model(verbose=verbose) + + measurement_df = self.petab_problem.measurement_df + + return amici.petab.simulations.rdatas_to_measurement_df( + rdatas, model, measurement_df + ) + + def rdatas_to_simulation_df( + self, + rdatas: Sequence[amici.ReturnData], + model: amici.Model = None, + ) -> pd.DataFrame: + """ + See :meth:`rdatas_to_measurement_df`. + + Except a petab simulation dataframe is created, i.e. the measurement + column label is adjusted. + """ + return self.rdatas_to_measurement_df(rdatas, model).rename( + columns={petab.MEASUREMENT: petab.SIMULATION} + ) + + def prediction_to_petab_measurement_df( + self, + prediction: PredictionResult, + predictor: AmiciPredictor = None, + ) -> pd.DataFrame: + """Not implemented yet.""" + raise NotImplementedError() + + def prediction_to_petab_simulation_df( + self, + prediction: PredictionResult, + predictor: AmiciPredictor = None, + ) -> pd.DataFrame: + """ + See :meth:`prediction_to_petab_measurement_df`. + + Except a PEtab simulation dataframe is created, i.e. the measurement + column label is adjusted. + """ + return self.prediction_to_petab_measurement_df( + prediction, predictor + ).rename(columns={petab.MEASUREMENT: petab.SIMULATION}) + + class PetabSimulatorObjectiveCreator(ObjectiveCreator): """ObjectiveCreator for creating an objective based on a PEtabSimulator.""" diff --git a/pypesto/petab/util.py b/pypesto/petab/util.py index 92e2a0a0a..ef280619e 100644 --- a/pypesto/petab/util.py +++ b/pypesto/petab/util.py @@ -4,10 +4,13 @@ try: import petab.v1 as petab + from petab import v2 from petab.v1.C import ( ESTIMATE, + LIN, NOISE_PARAMETERS, OBSERVABLE_ID, + PARAMETER_SCALE_UNIFORM, ) except ImportError: petab = None @@ -103,9 +106,9 @@ class PetabStartpoints(CheckedStartpoints): provided PEtab problem. The PEtab-problem is copied. """ - def __init__(self, petab_problem: petab.Problem, **kwargs): + def __init__(self, petab_problem: petab.Problem | v2.Problem, **kwargs): super().__init__(**kwargs) - self._parameter_df = petab_problem.parameter_df.copy() + self._petab_problem = petab_problem self._priors: list[tuple] | None = None self._free_ids: list[str] | None = None @@ -132,16 +135,33 @@ def _setup( # update priors self._free_ids = current_free_ids - id_to_prior = dict( - zip( - self._parameter_df.index[self._parameter_df[ESTIMATE] == 1], - petab.parameters.get_priors_from_df( - self._parameter_df, mode=petab.INITIALIZATION - ), - strict=True, + if isinstance(self._petab_problem, petab.Problem): + parameter_df = self._petab_problem.parameter_df + id_to_prior = dict( + zip( + parameter_df.index[parameter_df[ESTIMATE] == 1], + petab.parameters.get_priors_from_df( + parameter_df, mode=petab.INITIALIZATION + ), + strict=True, + ) ) - ) - + else: + id_to_prior = { + parameter.id: ( + # prior_type, prior_pars, par_scale, par_bounds + PARAMETER_SCALE_UNIFORM + if parameter.prior_distribution is None + else parameter.prior_distribution, + (parameter.lb, parameter.ub) + if parameter.prior_distribution is None + else parameter.prior_parameters, + LIN, + (parameter.lb, parameter.ub), + ) + for parameter in self._petab_problem.parameters + if parameter.estimate + } self._priors = list(map(id_to_prior.__getitem__, current_free_ids)) def __call__( diff --git a/test/petab/test_petab_import.py b/test/petab/test_petab_import.py index 77e6a2ec1..76ed05681 100644 --- a/test/petab/test_petab_import.py +++ b/test/petab/test_petab_import.py @@ -10,7 +10,7 @@ import amici.sim.sundials as asd import benchmark_models_petab as models import numpy as np -import petab.v1 as petab +import petab import petabtests import pytest @@ -266,6 +266,7 @@ def test_max_sensi_order(): def test_petab_pysb_optimization(): + """Test optimization for a PySB-based PEtab 2.0 problem.""" test_case = "0001" test_case_dir = petabtests.get_case_dir( test_case, version="v2.0.0", format_="pysb" @@ -276,7 +277,7 @@ def test_petab_pysb_optimization(): test_case, format="pysb", version="v2.0.0" ) - petab_problem = petab.Problem.from_yaml(petab_yaml) + petab_problem = petab.v2.Problem.from_yaml(petab_yaml) importer = PetabImporter(petab_problem) problem = importer.create_problem() @@ -298,6 +299,120 @@ def test_petab_pysb_optimization(): assert np.all(fvals <= -solution[petabtests.LLH]) +def test_petab_v2_boehm(): + import copy + import pickle + + from pypesto.optimize.optimizer import ScipyOptimizer + + # load test problem + problem_id = "Boehm_JProteomeRes2014" + petab_problem = petab.v2.Problem.from_yaml( + models.get_problem_yaml_path(problem_id) + ) + expected_fval_nominal = 138.22199693517703 + + # create model + importer = PetabImporter(petab_problem) + problem = importer.create_problem() + assert problem.x_names == petab_problem.x_ids + assert problem.dim == petab_problem.n_estimated == 9 + assert problem.dim_full == len(petab_problem.parameters) == 11 + # petab-non-estimated parameters are fixed in the pypesto problem + assert problem.x_fixed_indices == [ + i for i, p in enumerate(petab_problem.parameters) if not p.estimate + ] + assert isinstance( + problem.objective, pypesto.objective.amici.amici.AmiciPetabV2Objective + ) + + # evaluate objective at nominal parameters + fval = problem.objective(np.asarray(petab_problem.x_nominal_free)) + assert np.isclose(fval, expected_fval_nominal) + + # deepcopy works? + problem_deepcopy = copy.deepcopy(problem) + fval = problem_deepcopy.objective(np.asarray(petab_problem.x_nominal_free)) + assert np.isclose(fval, expected_fval_nominal) + + # pickling works? + problem_pickled = pickle.loads(pickle.dumps(problem)) # noqa: S301 + fval = problem_pickled.objective(np.asarray(petab_problem.x_nominal_free)) + assert np.isclose(fval, expected_fval_nominal) + + # gradient works? + fval, grad = problem.objective( + np.asarray(petab_problem.x_nominal_free), sensi_orders=(0, 1) + ) + assert np.isclose(fval, expected_fval_nominal) + assert len(grad) == petab_problem.n_estimated + + # fixing parameters, ... + problem.unfix_parameters(petab_problem.x_fixed_indices) + assert problem.dim == len(petab_problem.parameters) + with pytest.raises(ValueError, match="Cannot compute gradient"): + # cannot compute sensitivities for fixed parameters + problem.objective( + np.asarray(petab_problem.x_nominal), sensi_orders=(0, 1) + ) + fval = problem.objective(np.asarray(petab_problem.x_nominal)) + assert np.isclose(fval, expected_fval_nominal) + # re-fixing parameters + problem.fix_parameters( + petab_problem.x_fixed_indices, petab_problem.x_nominal_fixed + ) + assert problem.dim == petab_problem.n_estimated + fval, grad = problem.objective( + np.asarray(petab_problem.x_nominal_free), sensi_orders=(0, 1) + ) + assert np.isclose(fval, expected_fval_nominal) + assert len(grad) == petab_problem.n_estimated + + # TODO mode=res/fun, + # TODO hess/fim... + + # single optimization works? + optimizer = ScipyOptimizer() + result = optimizer.minimize( + id="1", + problem=problem, + x0=np.asarray(petab_problem.x_nominal_free) + 0.1, + ) + print(result) + assert result.fval0 is not None, result.message + assert result.fval < result.fval0 + + # multi-processing optimization works? + result = pypesto.optimize.minimize( + problem=problem, + optimizer=optimizer, + n_starts=4, + engine=pypesto.engine.MultiProcessEngine(), + progress_bar=False, + ) + assert len(result.optimize_result.list) == 4 + for local_result in result.optimize_result.list: + assert local_result.fval0 is not None, local_result.message + assert local_result.fval < local_result.fval0 + + +@pytest.mark.filterwarnings( + "ignore:.*distribution.*is not supported in PEtab v2.*:UserWarning" +) +def test_petab_v2_schwen(): + problem_id = "Schwen_PONE2014" + petab_problem = petab.v2.Problem.from_yaml( + models.get_problem_yaml_path(problem_id) + ) + assert petab_problem.n_priors + importer = PetabImporter(petab_problem) + problem = importer.create_problem() + assert problem.x_names == petab_problem.x_ids + assert isinstance(problem.objective, pypesto.objective.AggregatedObjective) + fval = problem.objective(np.asarray(petab_problem.x_nominal_free)) + assert np.isfinite(fval) + + if __name__ == "__main__": suite = unittest.TestSuite() suite.addTest(PetabImportTest()) diff --git a/test/petab/test_petab_suite.py b/test/petab/test_petab_suite.py index a0cc19dec..aa0ecc16d 100644 --- a/test/petab/test_petab_suite.py +++ b/test/petab/test_petab_suite.py @@ -5,7 +5,6 @@ import petab.v1 as petab import petabtests import pytest -from amici.sim.sundials.petab.v1 import rdatas_to_measurement_df import pypesto import pypesto.petab @@ -20,8 +19,8 @@ (case, model_type, version) for (model_type, version) in ( ("sbml", "v1.0.0"), - # FIXME: disabled until there is proper PEtab v2 support https://github.com/ICB-DCM/pyPESTO/pull/1634 - # ("pysb", "v2.0.0") + ("sbml", "v2.0.0"), + ("pysb", "v2.0.0"), ) for case in petabtests.get_cases(format_=model_type, version=version) ], @@ -69,29 +68,33 @@ def _execute_case(case, model_type, version): model_name = ( f"petab_test_case_{case}_{model_type}_{version.replace('.', '_')}" ) - output_folder = f"amici_models/{model_name}" # import and create objective function if case.startswith("0006"): + if version == "v2.0.0": + pytest.skip("TODO") petab_problem = petab.Problem.from_yaml(yaml_file) petab.flatten_timepoint_specific_output_overrides(petab_problem) importer = pypesto.petab.PetabImporter( petab_problem=petab_problem, - output_folder=output_folder, model_name=model_name, ) petab_problem = petab.Problem.from_yaml(yaml_file) else: importer = pypesto.petab.PetabImporter.from_yaml( - yaml_file, output_folder=output_folder + yaml_file, model_name=model_name ) petab_problem = importer.petab_problem + factory = importer.create_objective_creator() model = factory.create_model(generate_sensitivity_code=False) obj = factory.create_objective(model=model) - # the scaled parameters - problem_parameters = factory.petab_problem.x_nominal_scaled + if version == "v1.0.0": + # the scaled parameters + problem_parameters = factory.petab_problem.x_nominal_scaled + else: + problem_parameters = importer.petab_problem.x_nominal # simulate ret = obj(problem_parameters, sensi_orders=(0,), return_dict=True) @@ -100,20 +103,24 @@ def _execute_case(case, model_type, version): rdatas = ret["rdatas"] chi2 = sum(rdata["chi2"] for rdata in rdatas) llh = -ret["fval"] - simulation_df = rdatas_to_measurement_df( - rdatas, model, importer.petab_problem.measurement_df - ) + if version == "v1.0.0": + from amici.sim.sundials.petab.v1 import rdatas_to_measurement_df - if case.startswith("0006"): - simulation_df = petab.unflatten_simulation_df( - simulation_df, petab_problem + simulation_df = rdatas_to_measurement_df( + rdatas, model, importer.petab_problem.measurement_df ) + if case.startswith("0006"): + simulation_df = petab.unflatten_simulation_df( + simulation_df, petab_problem + ) - petab.check_measurement_df(simulation_df, petab_problem.observable_df) - simulation_df = simulation_df.rename( - columns={petab.MEASUREMENT: petab.SIMULATION} - ) - simulation_df[petab.TIME] = simulation_df[petab.TIME].astype(int) + petab.check_measurement_df(simulation_df, petab_problem.observable_df) + simulation_df = simulation_df.rename( + columns={petab.MEASUREMENT: petab.SIMULATION} + ) + simulation_df[petab.TIME] = simulation_df[petab.TIME].astype(int) + else: + simulation_df = obj.rdatas_to_simulation_df(rdatas) # check if matches chi2s_match = petabtests.evaluate_chi2(chi2, gt_chi2, tol_chi2) diff --git a/tox.ini b/tox.ini index 8942f8cf4..bb24b03fb 100644 --- a/tox.ini +++ b/tox.ini @@ -54,11 +54,7 @@ description = extras = test,amici,petab,emcee,dynesty,mltools,pymc,jax,fides,roadrunner deps = git+https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab.git@master\#subdirectory=src/python -# FIXME: Until we have proper PEtab v2 import -# we need `10ffb050380933b60ac192962bf9550a9df65a9c` -# with the pseudo-v2 format for pysb tests instead of `main` -# git+https://github.com/PEtab-dev/petab_test_suite@main - git+https://github.com/PEtab-dev/petab_test_suite@10ffb050380933b60ac192962bf9550a9df65a9c + git+https://github.com/PEtab-dev/petab_test_suite@a80dd936766d347a86b991d140d3e4a2eeda26bf commands = pytest --cov=pypesto --cov-report=xml --cov-append \ test/base --durations=0 \ @@ -85,17 +81,15 @@ deps = # always install amici from main branch, avoid caching # to skip re-installation, run `tox -e petab --override testenv:petab.commands_pre=` commands_pre = - python3 -m pip uninstall -y amici + python3 -m pip uninstall -y amici petabtests python3 -m pip install git+https://github.com/AMICI-dev/amici.git@main\#egg=amici&subdirectory=python/sdist + python3 -m pip install git+https://github.com/PEtab-dev/petab_test_suite@a80dd936766d347a86b991d140d3e4a2eeda26bf python -m pip list commands = -# FIXME: Until we have proper PEtab v2 import -# we need `10ffb050380933b60ac192962bf9550a9df65a9c` -# with the pseudo-v2 format for pysb tests instead of `main` -# python3 -m pip install git+https://github.com/PEtab-dev/petab_test_suite@main - python3 -m pip install git+https://github.com/PEtab-dev/petab_test_suite@10ffb050380933b60ac192962bf9550a9df65a9c python3 -m pip install git+https://github.com/pysb/pysb@master python3 -m pip install -U copasi-basico[petab] + python3 -m pip uninstall -y petab + python3 -m pip install git+https://github.com/petab-dev/libpetab-python@main\#egg=petab # upgrade after installing pysb which requires an older sympy python3 -m pip install -U sympy pytest --cov=pypesto --cov-report=xml --cov-append \