import numpy as np
from scipy import optimize
import pdfo
import x3cflux
import hopsy

Interoperation with Optimization Packages

Via Python, a huge variety of optimization packages are available. Still, some restrictions have to be met as these software packages are usually implemented from a methodological point of view. 13CFLUX, however, was designed to match this perspective and therefore seamlessly operates with many packages. This notebook demonstrates, how 13CFLUX integrates into two optimization packages.

simulator = x3cflux.create_simulator_from_fml("spiralus.fml", "ms+uptake")
mle = x3cflux.get_parameters(simulator.parameter_space, simulator.configurations[0].parameter_entries)

ineq_sys = simulator.parameter_space.inequality_system

Generate starting point and constraints using hopsy:

problem = hopsy.Problem(ineq_sys.matrix, ineq_sys.bound)
problem = hopsy.add_box_constraints(problem, 2*[0], 2*[5])
starting_point = hopsy.compute_chebyshev_center(problem)

Using solvers from scipy

scipy.minimize is one of the most frequently used framework for general purpose optimization. With the interface functions and definition of linear inequality constraints, 13CFLUX integrates seamlessly. Here, trust-constr was used, but others are available to handle linear inequalities too (e.g. SLSQP and COBYQA).

simulator.parameter_space.constraint_violation_tolerance = 1e-3
x3cflux.omp.num_max_threads = 10
result = optimize.minimize(simulator.compute_loss, starting_point.flatten(), jac=simulator.compute_loss_gradient,
                           constraints=[optimize.LinearConstraint(ineq_sys.matrix, ub=ineq_sys.bound, keep_feasible=True)],
                           method="trust-constr", tol=1e-5, options={"maxiter": 50, "verbose": 2})

| niter |f evals|CG iter|  obj func   |tr radius |   opt    |  c viol  |
|-------|-------|-------|-------------|----------|----------|----------|
|   1   |   1   |   0   | +3.2424e+03 | 1.00e+00 | 8.25e+02 | 0.00e+00 |
|   2   |   2   |   1   | +1.8018e+03 | 2.00e+00 | 6.82e+02 | 0.00e+00 |
|   3   |   3   |   3   | +9.4316e+02 | 9.23e+00 | 5.03e+02 | 0.00e+00 |
|   4   |   4   |   5   | +6.2678e+02 | 9.23e+00 | 3.52e+02 | 0.00e+00 |
|   5   |   5   |   7   | +4.4453e+02 | 9.23e+00 | 2.58e+02 | 0.00e+00 |
|   6   |   6   |   9   | +3.1440e+01 | 1.01e+01 | 1.05e+02 | 0.00e+00 |
|   7   |   7   |  11   | +3.1440e+01 | 1.01e+00 | 1.05e+02 | 0.00e+00 |
|   8   |   8   |  13   | +1.8182e+01 | 1.01e+00 | 7.21e+01 | 0.00e+00 |
|   9   |   9   |  15   | +1.8182e+01 | 2.85e-01 | 7.21e+01 | 0.00e+00 |
|  10   |  10   |  17   | +3.1015e+00 | 1.23e+00 | 2.85e+01 | 0.00e+00 |
|  11   |  11   |  19   | +5.2174e-01 | 1.23e+00 | 1.05e+01 | 0.00e+00 |
|  12   |  12   |  21   | +1.2135e-02 | 1.23e+00 | 1.47e+00 | 0.00e+00 |
|  13   |  13   |  23   | +1.3566e-03 | 1.23e+00 | 9.39e-02 | 0.00e+00 |
|  14   |  14   |  25   | +4.2115e-04 | 1.23e+00 | 6.06e-03 | 0.00e+00 |
|  15   |  14   |  25   | +4.2115e-04 | 6.13e+00 | 5.58e-02 | 0.00e+00 |
|  16   |  15   |  27   | +1.6714e-05 | 6.13e+00 | 4.17e-03 | 0.00e+00 |
|  17   |  15   |  27   | +1.6714e-05 | 3.06e+01 | 1.34e-02 | 0.00e+00 |
|  18   |  16   |  29   | +6.3684e-07 | 3.06e+01 | 1.34e-04 | 0.00e+00 |
|  19   |  16   |  29   | +6.3684e-07 | 1.53e+02 | 2.10e-03 | 0.00e+00 |
|  20   |  17   |  31   | +2.5505e-08 | 1.53e+02 | 1.66e-05 | 0.00e+00 |
|  21   |  17   |  31   | +2.5505e-08 | 7.66e+02 | 3.85e-04 | 0.00e+00 |
|  22   |  18   |  33   | +1.0195e-09 | 7.66e+02 | 1.27e-06 | 0.00e+00 |

`gtol` termination condition is satisfied.
Number of iterations: 22, function evaluations: 18, CG iterations: 33, optimality: 1.27e-06, constraint violation: 0.00e+00, execution time: 0.13 s.

result.fun, np.linalg.norm(result.x - mle)

(1.019462704647814e-09, 1.7322214344899365e-06)

Using gradient-free solver LINCOA

Gradient-based optimization is not always viable if gradient costs are high. This is especially the case for high numbers of parameters in large network models. LINCOA is one of Powel’s derivative-free optimizers that is specifically designed for linear inequality constraints and is available via the Python package pdfo. Similar to the scipy example, 13CFLUX integrates without notable overhead. Warnings are issued when the solver hits numerical edge-cases (here: evaluating at the boundary of the polytope), which are handeled automatically.

simulator.parameter_space.constraint_violation_tolerance = 1
result = pdfo.pdfo(lambda x: simulator.compute_loss(x) if simulator.parameter_space.contains(x) else 1e20,
                   starting_point.flatten(),
                   constraints=optimize.LinearConstraint(ineq_sys.matrix, ub=ineq_sys.bound, lb=-np.infty),
                   method="lincoa", options={'maxfev': 20_000, 'radius_final': 1e-15})

Operation "SparseLES Solver" failed: Singular system
Operation "SparseLES Solver" failed: Singular system
Operation "SparseLES Solver" failed: Singular system
Operation "SparseLES Solver" failed: Singular system
Operation "SparseLES Solver" failed: Singular system
Operation "SparseLES Solver" failed: Singular system

result.fun, np.linalg.norm(result.x - mle)

(2.0830858278492343e-28, 7.447602459741819e-16)