MICOM elasticities change current implementation to ensure same results

[Robaina]

Currently, elasticities are computed by optimizing growth (tradeoff) upon perturbing flux bounds of effector reaction (or taxa abundances).

    for r in fluxes:
        flux = import_fluxes[r]
        with com:
            if flux < -1e-6:
                r.lower_bound *= np.exp(STEP)
            else:
                r.upper_bound *= np.exp(STEP)
            sol = optimize_with_fraction(com, fraction, growth_rate, True)
            after = _get_fluxes(sol, reactions)
        deriv, dirs = _derivatives(before, after)

This strategy is not ideal since 1) we may not get the actual elasticity value and 2) we may run into alternative solutions so different elasticity values for the same optimal growth are possible. We should optimize the reaction under evaluation upon perturbation of the flux bound (taxa abundance) of the effector while constraining community growth (and / or individual growth) to a previously found optimum.

NOTES:

• However doing it that way requires solving way more optimization problems, as many problems as pairs of effectors and evaluated reactions...
• Having issues with the mutable com object across spawned processes. Any other solution besides deep copying com inside the helper function to be passed to multiprocessing? A possible workaround: update all reaction bounds just before solving the system, might still be problematic, tho.

[Robaina]

Have found an unexpected Python behavior that was causing discrepancies between the parallel and non-parallels versions of the algorithm. In this block:

    for rxn, flux in active_import_fluxes.items():
        with com:
            if flux < 0:
                rxn.lower_bound *= np.exp(STEP)
            else:
                rxn.upper_bound *= np.exp(STEP)
            sol = optimize_with_fraction(com, fraction, growth_rate, True)
            after_fluxes = _get_fluxes(sol, reactions)
        deriv, dirs = _derivatives(before_fluxes, after_fluxes, flux_tolerance)
        res = pd.DataFrame(
            {
                "reaction": [r.global_id for r in reactions],
                "taxon": [list(rx.compartments)[0] for rx in reactions],
                "effector": rxn.id,
                "direction": dirs,
                "elasticity": deriv,
            }
        )
        results.append(res)

assigning the same name to the variables inside the two list comprehensions than the one in the for loop, i.e. r, changes the value of the variable defined in the for loop, here causing artificial "effectors". This is an unexpected behavior (to me, at least) since I would expect the scope of variables within list comprehensions to remain within the list comprehension.

[Robaina]

Found another unexpected behavior when running the original elasticities_by_abundance and the modified elasticities_by_abundance_parallel. It looks like elasticities are very sensitive to small numerical variations. Also still not convinced by micom's approach due to alternative optima. To get robust elasticity values we definitely should maximize each reaction's flux value upon the perturbation while constraining community growth and individual growth.

UPDATE: Pretty sure the differences between parallel and sequential are due to gurobi's warm start, as these don't happen when using the "hybrid" (HIGHS + OSQP) solver. These differences are basically due to alternative optima. Need to reimplement elasticities and maximize each target reaction individually to avoid this effect.