search

ICB-DCM / pyPESTO

python Parameter EStimation TOolbox
https://pypesto.readthedocs.io
BSD 3-Clause "New" or "Revised" License
217 stars 47 forks source link

Multistart optimization with parameter constraints #1334

Closed aidinbii closed 6 months ago

aidinbii commented 6 months ago

Hello, is there a way to specify parameter constraints for parameter estimation. For example, if I have a parameter vector: (p1, p2, p3, p4). And I want a constraint: p2 <= p3.

Similar to this

Thank you

Doresic commented 6 months ago

Hi aidinbii,

In the current implementation, parameter-dependent constraints are not allowed.

However, there is an easy way to go around this problem. In your model, you can define the p3 parameter using the following assignment rule p3 := p2 + p23_diff. Then if you estimate the two parameters p2 and p23_diff, you can satisfy the p2<=p3 bound by setting the lower bound for p23_diff to 0.

Hope this helps.

Best, Domagoj

aidinbii commented 6 months ago

@Doresic, thanks for the reply But, will it work together with conditions?

For instance:

odes:
    - stateId: mu_1
      rightHandSide: -Ka*mu_1
      initialValue: 0

    - stateId: mu_2
      rightHandSide: Ka*mu_1 - Ke*mu_2
      initialValue: 0

parameters:
    - parameterId: Ka
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 0

    - parameterId: Ke
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 0

    - parameterId: VdF
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 0

    - parameterId: Ka_3
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: Ka_30
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: Ka_100
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1 

    - parameterId: Ke_3
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: Ke_30
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: Ke_100
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: VdF_3
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: VdF_30
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: VdF_100
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.0001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

    - parameterId: diff_Ke1
      nominalValue: 0.1
      parameterScale: log
      lowerBound: 0.000001
      upperBound: 10
      initializationPriorType: uniform
      initializationPriorParameters: 0.05; 1
      estimate: 1

assignments:
    - assignmentId: Ke_3
      formula: Ke_30 + diff_Ke1

conditions:
    - conditionId: dose_3
      Ka: Ka_3
      Ke: Ke_3
      VdF: VdF_3
      mu_1: 3e+6 

    - conditionId: dose_30
      Ka: Ka_30
      Ke: Ke_30
      VdF: VdF_30
      mu_1: 30e+6

    - conditionId: dose_100
      Ka: Ka_100
      Ke: Ke_100
      VdF: VdF_100
      mu_1: 100e+6

I get an error:

Checking model...
libSBML Error (SBML identifier consistency): The value of the 'id' field on every instance of the following type of object in a model must be unique: <model>, <functionDefinition>, <compartmentType>, <compartment>, <speciesType>, <species>, <reaction>, <speciesReference>, <modifierSpeciesReference>, <event>, and model-wide <parameter>s. Note that <unitDefinition> and parameters defined inside a reaction are treated separately. 
Reference: L3V1 Section 3.3
   The <parameter> id 'Ke_3' conflicts with the previously defined <parameter> id 'Ke_3' at line 19.

WARNING: Generated invalid SBML model. Check messages above.
Doresic commented 6 months ago

Yes, it should work. The error you're getting now is because you're defining Ke_3 twice: once in the parameter table, and once in the assignment rule. Notice that if you define Ke_3 := Ke_30 + diff_Ke1, then you should not define Ke_3 in the parameter table as well.

So just removing 

- parameterId: Ka_3
  nominalValue: 0.1
  parameterScale: log
  lowerBound: 0.0001
  upperBound: 10
  initializationPriorType: uniform
  initializationPriorParameters: 0.05; 1
  estimate: 1

Should make it work.

aidinbii commented 6 months ago

@Doresic , there is another issue when running

problem = importer.create_problem()

I get

KeyError                                  Traceback (most recent call last)
[~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py) in _apply_parameter_table(par_mapping, scale_mapping, parameter_df, scaled_parameters, fill_fixed_parameters)
    595             # the overridee is a model parameter
--> 596             par_mapping[problem_par] = par_mapping[sim_par]
    597             scale_mapping[problem_par] = scale_mapping[sim_par]

KeyError: 'Ke_3'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
[~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/indexes/base.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/indexes/base.py) in get_loc(self, key)
   3789         try:
-> 3790             return self._engine.get_loc(casted_key)
   3791         except KeyError as err:

index.pyx in pandas._libs.index.IndexEngine.get_loc()

index.pyx in pandas._libs.index.IndexEngine.get_loc()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item()

KeyError: 'Ke_3'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
<ipython-input-44-d2644dde1c84> in <module>
      1 # Create a problem for param.est
----> 2 problem = importer.create_problem()

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/pypesto/petab/importer.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/pypesto/petab/importer.py) in create_problem(self, objective, x_guesses, problem_kwargs, startpoint_kwargs, **kwargs)
    690         """
    691         if objective is None:
--> 692             objective = self.create_objective(**kwargs)
    693 
    694         x_fixed_indices = self.petab_problem.x_fixed_indices

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/pypesto/petab/importer.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/pypesto/petab/importer.py) in create_objective(self, model, solver, edatas, force_compile, **kwargs)
    439             )
    440 
--> 441         parameter_mapping = amici.petab_objective.create_parameter_mapping(
    442             petab_problem=self.petab_problem,
    443             simulation_conditions=simulation_conditions,

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/amici/petab/parameter_mapping.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/amici/petab/parameter_mapping.py) in create_parameter_mapping(petab_problem, simulation_conditions, scaled_parameters, amici_model, **parameter_mapping_kwargs)
    370 
    371     prelim_parameter_mapping = (
--> 372         petab.get_optimization_to_simulation_parameter_mapping(
    373             condition_df=petab_problem.condition_df,
    374             measurement_df=petab_problem.measurement_df,

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py) in get_optimization_to_simulation_parameter_mapping(condition_df, measurement_df, parameter_df, observable_df, mapping_df, sbml_model, simulation_conditions, warn_unmapped, scaled_parameters, fill_fixed_parameters, allow_timepoint_specific_numeric_noise_parameters, model)
    174             ),
    175         )
--> 176         return list(mapping)
    177 
    178     # Run multi-threaded

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py) in _map_condition(packed_args)
    287         )
    288 
--> 289     par_map_sim, scale_map_sim = get_parameter_mapping_for_condition(
    290         condition_id=condition[SIMULATION_CONDITION_ID],
    291         is_preeq=False,

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py) in get_parameter_mapping_for_condition(condition_id, is_preeq, cur_measurement_df, sbml_model, condition_df, parameter_df, mapping_df, simulation_parameters, warn_unmapped, scaled_parameters, fill_fixed_parameters, allow_timepoint_specific_numeric_noise_parameters, model)
    420         mapping_df,
    421     )
--> 422     _apply_parameter_table(
    423         par_mapping,
    424         scale_mapping,

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/petab/parameter_mapping.py) in _apply_parameter_table(par_mapping, scale_mapping, parameter_df, scaled_parameters, fill_fixed_parameters)
    602             # or the overridee is only defined in the parameter table
    603             scale = (
--> 604                 parameter_df.loc[sim_par, PARAMETER_SCALE]
    605                 if PARAMETER_SCALE in parameter_df
    606                 else LIN

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/indexing.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/indexing.py) in __getitem__(self, key)
   1144             key = tuple(com.apply_if_callable(x, self.obj) for x in key)
   1145             if self._is_scalar_access(key):
-> 1146                 return self.obj._get_value(*key, takeable=self._takeable)
   1147             return self._getitem_tuple(key)
   1148         else:

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/frame.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/frame.py) in _get_value(self, index, col, takeable)
   4013             #  results if our categories are integers that dont match our codes
   4014             # IntervalIndex: IntervalTree has no get_loc
-> 4015             row = self.index.get_loc(index)
   4016             return series._values[row]
   4017 

[~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/indexes/base.py](https://file+.vscode-resource.vscode-cdn.net/home/aidin/Documents/tum_courses/hiwi_Atefeh/alivexbiotech/projects/igg4_python/notebooks/case_studies_1/~/miniconda3/envs/modeling/lib/python3.9/site-packages/pandas/core/indexes/base.py) in get_loc(self, key)
   3795             ):
   3796                 raise InvalidIndexError(key)
-> 3797             raise KeyError(key) from err
   3798         except TypeError:
   3799             # If we have a listlike key, _check_indexing_error will raise

KeyError: 'Ke_3'

I guess I also need to change the conditions?

Doresic commented 6 months ago

Did you re-compile the model since changing the petab problem (removing the parameter from the parameter df)? You can do so by adding a force_compile=True into the model creation

problem = importer.create_problem(force_compile=True)

or deleting the compiled model (by default it is in the amici_models folder).

aidinbii commented 6 months ago

Yes, I did. The model is compiled without problem

The issue is with the problem

Doresic commented 6 months ago

Hmm tested it myself now. It should work nicely. What is your current parameter.tsv file, and model SBML file?

aidinbii commented 6 months ago

@Doresic

aidinbii commented 6 months ago

I removed Ka,Ke,Vdf from parameters table because after: !petablint -v -y ./petab_tables/case_studies_1/PK_1C_constr_1_petab/problem_PK_1C_constr_1.yml I got 

Checking model...
Checking measurement table...
Checking condition table...
Checking observable table...
Checking parameter table...
Extraneous parameter(s) in parameter table: {'Ke', 'Ka', 'VdF'}
{'Ka', 'VdF', 'Ke'} are not allowed to occur in the parameters table.
Visualization table not available. Skipping.
Not OK
Doresic commented 6 months ago

My mistake. I found the error.

The definition of Ke_3 in the model should not be with an assignment Ke_3 := Ke_30 + diff_Ke1, but with an initial assignment Ke_3 = Ke_30 + diff_Ke1 (notice different equation sign, this is the example for Antimony). The thing that's not supported is having time-varying parameters, which an assignment can possibly be.

Removing Ka, Ke, Vdf from the parameter table is correct, as those are not the parameters you will estimate.

aidinbii commented 6 months ago

My mistake. I found the error.

The definition of Ke_3 in the model should not be with an assignment Ke_3 := Ke_30 + diff_Ke1, but with an initial assignment Ke_3 = Ke_30 + diff_Ke1 (notice different equation sign, this is the example for Antimony). The thing that's not supported is having time-varying parameters, which an assignment can possibly be.

Removing Ka, Ke, Vdf from the parameter table is correct, as those are not the parameters you will estimate.

But how to define an initial assignment in the yaml file?

assignments:
    - assignmentId: Ke_3
      formula: Ke_30 + diff_Ke1

Also, where is the example for Antimony?

Doresic commented 6 months ago

Not sure. I haven't used yaml2sbml much myself. If it's possible it should be somewhere in the documentation https://yaml2sbml.readthedocs.io/en/latest/

Alternatively, you can use Antimony from Tellurium. A simple and similar way of creating SBML models. https://tellurium.readthedocs.io/en/latest/antimony.html . With Antimony it definitely works, tried it myself.

As this is no longer pypesto's problem, I'll close this issue. But feel free to open up if there are any more problems with pypesto.

Best, Domagoj