Requesting some help to setup a complicated constraint

[JulianoGianlupi]

Hello! First of all, I am sorry if I am a bit vague with my descriptions, proprietary information and all that.

So, I have a logistic growth model for cell populations that I am trying to parametrize through fitting. I.e.,

$$\frac{dN}{dt} = r N \left(1-N/C\right)$$

I have data for growth under different conditions of "chemical mixtures," and r (growth rate), C (carry capacity), and the initial population (N0) are a function of the concentrations of the different chemicals in these mixtures. E.g.,

$$C(\vec{\kappa}) = C_{base} + \sum_{i} \beta_i \kappa_i\,\,,$$

with $\vec{\kappa}$ the vector of chemical concentrations, and $\beta_i$ the contribution of the $i$-th chemical to the carry capacity. Now, cells occupy space, so the carry capacity has a limit, and here's my problem, how can I define bounds for $C_{base}$ and the $\beta$-s?

For now what I am doing is calculating C during the cost difference calculation like this:

def calc_carry_cap(params, concs, 
                    max_cc:float|None=None) -> float:
    cc = (params["base_cc"] +
          params["b1_cc"] * concs["1"] +
          params["b2_cc"] * concs["2"] +
          ...)
    if max_cc is not None:
        return min(cc, max_cc)
    return cc

Which is not ideal, as the minimization algorithm can just go >max_cc indefinitely and not explore actually reasonable values for the betas. My issue is that I can't simply constrain the sum of my parameters to an inequality as in the examples because I have influence from the input (data) itself.

In other words doing

from lmfit import Parameters
params = Parameters()
params.add("max_cc", value=max_cc, vary=False)
params.add("base_cc", expr="max_cc-base_cc-b1_cc-...")
...

won't work because there is the multiplication of the betas by the concentration.

I know there is user_data with the Parameter class, but I couldn't find any examples on how to use it. Would that be a possible avenue to reach a solution?

Many thanks!

[newville]

@JulianoGianlupi
Are you looking to minimize the returned value of the objective function calc_cc_cap()? Sorry to be dense, but are the equations you provide related to the Python code? You have two equations, both with letters C in them - are those related? Is there data involved? Should there be some growth rate? If any of that is meant to be obvious, well, I am not following you.

We always recommend a complete, minimal example. We want to see how you define parameters, what your objective function or model function is, and how these are used to do the fit. With incomplete snippets of code, we do not know what you are doing.

[newville]

@JulianoGianlupi Please read the lmfit documentation on constraint expressions: they are designed to set a value for a Parameter from a mathematical expression of other Parameters. They cannot be self-referential, and they do not by themselves impose limits on Parameter values.

If there is a way to define the constraints as in scipy (ie, as Callables) I believe this would be trivial to solve

No.

As said earlier, if you want help, you really have to provide a working, minimal example. If the code does not run, it is not complete. If you don't have import statements, it cannot run. The example you posted appears to not even use lmfit. I have no clue what you are doing. Why is that?

[JulianoGianlupi]

Aight, here's a mwe followed by the working scipy alternative to the constraint

# %%
import numpy as np
import lmfit as lf
import pandas as pd
from scipy.integrate import solve_ivp
from scipy.optimize import minimize
from typing import Callable
from copy import deepcopy
rng = np.random.default_rng()

# %%
# define the growth function and the linear combinations of 
# parameters
def logistic_growth(t, u, params):
    du = np.zeros_like(u)
    du[0] = params["rate"] * u[0] * (1- u[0]/params["cc"])
    return du


def calc_carry_cap(params, concs, 
                    max_cc:float|None=None) -> float:
    cc = (params["base_cc"] +
          params["A_cc"] * concs["A"] +
          params["B_cc"] * concs["B"] +
          params["C_cc"] * concs["C"] )
    if max_cc is not None:
        return min(cc, max_cc)
    return cc

def calc_growth_rate(params, concs) -> float:
    r = (params["base_r"] +
          params["A_r"] * concs["A"] +
          params["B_r"] * concs["B"] +
          params["C_r"] * concs["C"] )
    
    return r

def calc_n0(params, concs) -> float:
    fraction = (params["base_n0"] +
          params["A_n"] * concs["A"] +
          params["B_n"] * concs["B"] +
          params["C_n"] * concs["C"] )
    
    # params["n0"] is fixed (vary=False) and is the number of cells that 
    # where delivered to the petri dish. So this function is calculating what
    # fraction of them successfully "attached" to the dish
    n0 = params["n0"] * min(1, fraction)
    return n0


# %%
# define (artificial) chemical names, mixtures, cells
chem_names = ["A", "B", "C"]
mixes = {"X": [1,2,3], "Y":[4,5,6], "Z":[.1,.1,.1]}

chem_df = pd.DataFrame(mixes, index=chem_names)

cell_line_names = ["bac1","bac2","bac3",]


# %%
# define artificial data and parameters
data = {}
target_parameters = {}
t = np.array([6, 24, 72]) # h

experiment_id = 0
experimental_results = pd.DataFrame()

for cell_line in cell_line_names:

    # "actual" parameters for this cell
    parameters = lf.Parameters()
    parameters.add(name="base_cc", value=1e3*rng.uniform(.8,1.2), min=0)
    parameters.add(name="n0", value=10*rng.uniform(.8,1.2), vary=False)

    parameters.add(name="base_r", value=0.1*rng.uniform(.8,1.2), min=0)
    n0_frac_tally = rng.uniform() # to keep the final max fraction <=1
    parameters.add(name="base_n0", value=n0_frac_tally)

    for chem_name in chem_names:
        parameters.add(name=f"{chem_name}_cc", value=1e3*rng.uniform(.8,1.2), min=0)
        parameters.add(name=f"{chem_name}_r", value=rng.uniform(.8,1.2), min=0)
        n0_frac = rng.uniform(0,n0_frac_tally)
        n0_frac /= chem_df["Y"][chem_name] # to keep the final max fraction <=1
        n0_frac_tally += n0_frac
        parameters.add(name=f"{chem_name}_n", value=n0_frac, min=0)
    
    max_cc = 0.9*calc_carry_cap(parameters, chem_df["Y"]) 
    max_cc = 0.9*max_cc # say the geometry would only allow for 90% of the amount calculated using the highest concentration mix
    parameters.add(name="max_cc", value=max_cc, vary=False) 

    target_parameters[cell_line] = parameters

    # solve the system

    for mix_name in chem_df.columns: 

        n0 = calc_n0(parameters, chem_df[mix_name])
        cc = calc_carry_cap(parameters, chem_df[mix_name], max_cc=max_cc)
        r = calc_growth_rate(parameters, chem_df[mix_name])

        sol = solve_ivp(logistic_growth, [0,73], [n0], 
                        args=({"rate": r, "cc": cc},), t_eval=t)
        


        for experimental_repeat in range(2):
            noisy_y = rng.normal(sol.y[0], 15) # we don't need that much noise 
            experimental_results = pd.concat(
                [experimental_results, pd.DataFrame(
                    {
                        "t": sol.t,
                        "y": noisy_y,
                        "exp": [experiment_id]*len(noisy_y),
                        "mix": [mix_name]*len(noisy_y),
                        "cell_line": [cell_line]*len(noisy_y)
                    }
                )]
            )
            experiment_id +=1
experimental_results.reset_index(inplace=True, drop=True)

# for later
parameter_indexes_bac1 = {
    k: i for i, (k, _) in enumerate(target_parameters["bac1"].valuesdict().items())
}

# %%
# now we can actually start

# As mentioned, I am using a Cost class to minimize iteration-to-iteration
# manipulations. Not super relevant in this reduced example, but this provided
# considerable speedup for me

class Cost:
    def __init__(self, data:pd.DataFrame, chem_df:pd.DataFrame, model:Callable):
        
        self.data = deepcopy(data)
        self.model = model
        self.chem_df = deepcopy(chem_df)
        self.y = self.data["y"].to_numpy()
        self.t = self.data["t"].to_numpy()
        self.base_error_vec = np.zeros(len(self.y))
        return
    
    def __call__(self, parameters):
        
        return self.calc_error_vec(parameters)
    
    def calc_error_vec(self, parameters):
        error = self.base_error_vec.copy()
        for mix_name, mix_result in self.data.groupby("mix"):
            n0 = calc_n0(parameters, chem_df[mix_name])
            cc = calc_carry_cap(parameters, chem_df[mix_name], 
                                max_cc=parameters["max_cc"])
            r = calc_growth_rate(parameters, chem_df[mix_name])

            sol = solve_ivp(logistic_growth, [0,73], [n0], 
                            args=({"rate": r, "cc": cc},), t_eval=t)
            
            if np.any(~np.isfinite(sol.y)):
                    # Return large error so optimizer sees it as invalid
                    return 1e20 * np.ones_like(error)
            
            for repeat_id, repeat_result in mix_result.groupby("exp"):
                error[repeat_result.index] = \
                    repeat_result["y"].to_numpy() - sol.y[0]
        return error



# %%
# now let's solve for each cell type
fit_results = {}
for cell_line in cell_line_names:
    initial_guess = deepcopy(target_parameters[cell_line])
    cell_line_results = experimental_results.loc[experimental_results["cell_line"] == cell_line].copy()
    cell_line_results.reset_index(inplace=True, drop=True)
    cost_fcn = Cost(cell_line_results,
                    chem_df,
                    logistic_growth)
    res = lf.minimize(fcn=cost_fcn, params=initial_guess)
    fit_results[cell_line] = res.params.copy()
    break # you get the idea

# %%
# How to restrict the parameters that are fed into `calc_carry_cap`, `calc_growth_rate`,
# and `calc_n0` using numpy:
# (I'll only do n0 and bac1)

# turn lmfit Parameters into what scipy expects

parameters_sp = [
    v for _, v in fit_results["bac1"].valuesdict().items()
]

bounds_sp = [
    (0, None) if fit_results["bac1"][k].vary else (v,v) 
        for k, v in fit_results["bac1"].valuesdict().items() 
]

# or use the original parameters
# parameters_sp = [
#     v for _, v in target_parameters["bac1"].valuesdict().items()
# ]

# bounds_sp = [
#     (0, None) if target_parameters["bac1"][k].vary else (v,v) 
#         for k, v in target_parameters["bac1"].valuesdict().items() 
# ]

# define constraint

def n0_fraction_ineq(params, concs, idx_dict):
    n0f = (params[idx_dict["base_n0"]] +
          params[idx_dict["A_n"]] * concs["A"] +
          params[idx_dict["B_n"]] * concs["B"] +
          params[idx_dict["C_n"]] * concs["C"] )
    return 1 - n0f

n0_fraction_ineq_bac1 = lambda params: n0_fraction_ineq(params, 
                                                        chem_df["Y"], 
                                                        parameter_indexes_bac1)
cons = ({'type': 'ineq', 'fun': n0_fraction_ineq_bac1},)



# %%
# re-define `calc_carry_cap`, `calc_growth_rate`,
# and `calc_n0` to be usable without Parameters

def logistic_growth(t, u, params):
    du = np.zeros_like(u)
    du[0] = params["rate"] * u[0] * (1- u[0]/params["cc"])
    # du[0] = params[0] * u[0] (1- u[0]/params[1])
    return du


def calc_carry_cap(params, concs, idx_dict, 
                    max_cc:float|None=None) -> float:
    cc = (params[idx_dict["base_cc"]] +
          params[idx_dict["A_cc"]] * concs["A"] +
          params[idx_dict["B_cc"]] * concs["B"] +
          params[idx_dict["C_cc"]] * concs["C"] )
    if max_cc is not None:
        return min(cc, max_cc)
    return cc

def calc_growth_rate(params, concs, idx_dict) -> float:
    r = ( params[idx_dict["base_r"]] +
          params[idx_dict["A_r"]] * concs["A"] +
          params[idx_dict["B_r"]] * concs["B"] +
          params[idx_dict["C_r"]] * concs["C"] )
    
    return r

def calc_n0(params, concs, idx_dict) -> float:
    fraction = (params[idx_dict["base_n0"]] +
                params[idx_dict["A_n"]] * concs["A"] +
                params[idx_dict["B_n"]] * concs["B"] +
                params[idx_dict["C_n"]] * concs["C"] )
    
    # params["n0"] is fixed (vary=False) and is the number of cells that 
    # where delivered to the petri dish. So this function is calculating what
    # fraction of them successfully "attached" to the dish
    n0 = params[idx_dict["n0"]] * fraction# * min(1, fraction); min no longer necessary because of the constraint
    return n0

# %%
# Will use SciPy minimize which expects a scalar returning function, and we are now using the new 
# `calc_carry_cap`, `calc_growth_rate`,and `calc_n0` so reworking the cost class

class CostSp(Cost):
    def __init__(self, data, chem_df, model, idx_dict):
        super().__init__(data, chem_df, model)
        self.idx_dict  = idx_dict
    def __call__(self, parameters):
        return (self.calc_error_vec(parameters)**2).sum()
    def calc_error_vec(self, parameters):
        error = self.base_error_vec.copy()
        for mix_name, mix_result in self.data.groupby("mix"):
            n0 = calc_n0(parameters, chem_df[mix_name], self.idx_dict)
            cc = calc_carry_cap(parameters, chem_df[mix_name], self.idx_dict, 
                                max_cc=parameters[self.idx_dict["max_cc"]])
            r = calc_growth_rate(parameters, chem_df[mix_name], self.idx_dict)

            sol = solve_ivp(logistic_growth, [0,73], [n0], 
                            args=({"rate": r, "cc": cc},), t_eval=t)
            
            if np.any(~np.isfinite(sol.y)):
                    # Return large error so optimizer sees it as invalid
                    return 1e20 * np.ones_like(error)
            
            for repeat_id, repeat_result in mix_result.groupby("exp"):
                error[repeat_result.index] = \
                    repeat_result["y"].to_numpy() - sol.y[0]
        return error

# %%
bac1_res = experimental_results.loc[experimental_results["cell_line"] == "bac1"].copy()
bac1_res.reset_index(inplace=True, drop=True)
cost_fcn = CostSp(bac1_res, chem_df, logistic_growth, idx_dict=parameter_indexes_bac1)

# %%
# this can take much longer than lmfit even when using the already minimized parameters
res_sp = minimize(cost_fcn, parameters_sp, 
                  bounds=bounds_sp, 
                  constraints=cons)

# %%
print(res_sp)

[newville]

@JulianoGianlupi It may be that saying " I am not following you." and then "I have no clue what you are doing" may not have properly expressed my confusion about what you are trying to do or even what your question is.

You can say things like "you get the idea" and I can assure you that I do not understand what you are trying to do or what your actual question is. You started with two mathematical functions -- are those even related to what you are doing? I truly do not know.

At 300 lines, your definition of "Minimal" is different from mine. It looks like you intend to do a bunch of fits and then do something (I do not know what) with those results. If I run your code, it just hangs for me with no output at all - not a clue where or how or what it is doing. You seem to expect that I am following what you are trying to do so that you do not need to explicitly say, "I did this (show code), I expected this (show/describe expected result), and I got this (show actual output or error messages)". It seems like you want to set some constraint on a sum of some parameters or a sum of products of parameters. This is described in the documentation. Please read that.

I do not know what to tell you or how to help you. If the docs and examples don't give you some insight, I doubt I can be of much help. Maybe you need some other tools that can better match your needs. Best of luck.

[JulianoGianlupi]

I understood that you have questions and I have been trying to honestly answer them. Let me try to put it this way:

• Forget the equation in my initial post
• Can this constraint definition defined for SciPy's minimize

x = np.linspace(0,1,100)

y = 0.5 * x
y = rng.normal(y, .01) 

# print(y)

def constraint_fun(params):
    f = (0.1*x[0] + 10 * x[1]+ 5 * x[2])
    return 1 - f # ie, 0.1*x[0] + 10 * x[1] + 5 * x[2] < 1

def error(params, x, yobs):
    y = (params[0] + params[1] + params[2]) * x
    return ((y - yobs)**2).sum()

l_error = lambda p: error(p, x, y)
p0=[9, .9, 1/5]
constraints = {'type': 'ineq', 'fun': constraint_fun}
res = minimize(l_error, p0, 
               constraints=constraints, 
               bounds=((0, None),(0, None),(0, None),))
print(res)
print(res.x.sum())

be translated into lmfit? If yes how? If no why?

I have read your docs and couldn't figure out how. Your example for inequality are for only 2 variables and it is not clear how to expand it into more. As mentioned before I tried

import lmfit as lf
chem_names = ["A", "B", "C"]
max_concs = [4,5,6]
parameters =  lf.Parameters()
parameters.add(name="seed_frac", value=1, min=0, max=1)

complete_seed_frac_expr_str = "seed_frac-base_n"
parameters.add(name="base_n", value=1, min=0) 
for i, chem_name in enumerate(chem_names):
   complete_seed_frac_expr_str += f"-beta_{chem_name}_n * {max_concs[i]}"
   parameters.add(name=f"beta_{chem_name}_n", value=0, min=0) 
print(complete_seed_frac_expr_str)
parameters["base_n"].expr= complete_seed_frac_expr_str.replace("-base_n", "")

for chem_name in chem_names:
   parameters[f"beta_{chem_name}_n"].expr= complete_seed_frac_expr_str.replace(f"-beta_{chem_name}_n", "")
parameters.pretty_print()

I.e.,

base_n + beta_A_n * 4 + beta_B_n * 5 + beta_C_n * 6 < seed_frac

But that runs into a recursion error when trying to minimize

[newville]

As you can read from https://lmfit.github.io/lmfit-py/constraints.html#using-inequality-constraints you could encode

     a + b + c*3 < sum

as

    a + b + c*3 + delta = sum`, with `delta > 0`.

which might translate to

pars = Parameters()
pars.add('b', value=1, vary=True)
pars.add('c', value=2, max=10, vary=True)
pars.add('sum', value = 9, vary = True)
pars.add('delta', value=0.1, min=0, vary=True) # forced positive
pars.add('a', expr = 'sum - (b + c*3  + delta)')

Maybe that would be useful to you. Otherwise, I don't know.