ICB-DCM / pyPESTO

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

Hierarchical optimization estimates wrong scaling factors #471

Open stobiblum opened 4 years ago

stobiblum commented 4 years ago

Good morning,

I am using the pypesto 'feature_hierarchical' branch on a WSL. I tried the hierarchical estimation with a model including five observables with three scaling parameters, where everything went fine. However, when I used a larger dataset including 13 observables with 11 scaling parameters the estimation of some of the scaling parameters failed (The SBML and PEtab files are attached). Noise was not estimated in any of the cases.

AhR_hiera_more_data.zip

For estimation I used the analytical inner solver with following code:

import pypesto.optimize as optimize
import pypesto.hierarchical.solver as solver

importer = pypesto.petab.PetabImporter.from_yaml(yaml_config)

obj = importer.create_objective()
obj.amici_solver.setRelativeTolerance(1e-6)
obj.amici_solver.setAbsoluteTolerance(1e-8)

problem = importer.create_problem(obj)
problem.objective.amici_solver.setSensitivityMethod(amici.SensitivityMethod_adjoint)
problem.objective.calculator.inner_solver = solver.AnalyticalInnerSolver()

optimizer = optimize.ScipyOptimizer()
engine = pypesto.engine.SingleCoreEngine()

result = optimize.minimize(problem=problem,
                          optimizer=optimizer,
                          n_starts=n_starts,
                          engine=engine)

As you can see in the plots below, some scaling parameters (namely s_mRNA_AhRR_1 and s_mRNA_Cyp1a1_1) are not estimated correctly. If I manually increase s_mRNA_Cyp1a1_1 from 3.496 to 6.233 the fval decreases from 6526.812775823983 to 6314.67373801991. If I manually increase s_mRNA_AhRR_1 from 0.353 to 5.729 the fval decreases from 6526.812775823983 to 6303.173798458569.

simulation_MC2_5_log_1000_amici_subplots simulation_MC10_log_long_1000_amici_subplots

Do you have any idea what is causing the problem?

Thank you in advance!

LeonardSchmiester commented 4 years ago

Hi @stobiblum, I don't see anything obvious in the PEtab files. Can you do a gradient check (obj.check_gradient(p)), to see if the gradients for the scalings are 0?

stobiblum commented 4 years ago

@LeonardSchmiester The gradient for the scalings is not 0. obj.check_grad(p) yields for s_mRNA_Cyp1a1_1 -1.792550e+00 and for s_mRNA_AhRR_1 -4.339656e-03.