Requesting some help to setup a complicated constraint #998
Replies: 1 comment 8 replies
-
|
@JulianoGianlupi 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. |
Beta Was this translation helpful? Give feedback.
-
|
@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.
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 |
Beta Was this translation helpful? Give feedback.
-
|
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
8000
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) |
Beta Was this translation helpful? Give feedback.
-
|
@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. |
Beta Was this translation helpful? Give feedback.
-
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
Beta Was this translation helpful? Give feedback.
-
|
As you can read from https://lmfit.github.io/lmfit-py/constraints.html#using-inequality-constraints you could encode as which might translate to Maybe that would be useful to you. Otherwise, I don't know. |
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
-
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.,
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.,
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:
Which is not ideal, as the minimization algorithm can just go
>max_ccindefinitely 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
won't work because there is the multiplication of the betas by the concentration.
I know there is
user_datawith theParameterclass, but I couldn't find any examples on how to use it. Would that be a possible avenue to reach a solution?Many thanks!
Beta Was this translation helpful? Give feedback.
All reactions