fit simulation data of SIR model and return with 0

[penguinaugustus]

Hi, I run into a problem: I tried to use sindy to find SIR model. The data looks like this:

I sample with 10000 points in 100 seconds (dt = 0.01) and create the fitting model like this:

feature_names = ['S','I','R']

optimizer = ps.STLSQ(threshold=0)
model = ps.SINDy(feature_names=feature_names, optimizer=optimizer)
model.fit(X,t=0.01)
model.print()

however, the result turns out to be:

(S)' = 0.000
(I)' = 0.000
(R)' = 0.000 

I am sure my data are written in correct format. Could you give me some hint on what i can do to fix the problem? Thank you for your time!

Sincerely, Yu

[penguinaugustus]

I am thinking about whether it is because my odes will reach an equilibrium state at a larger time scale

I also notice that the example here are about dynamic systems with periodic behaviors: https://arxiv.org/abs/2301.02673

[Jacob-Stevens-Haas]

Interesting problem! It may be that most of your data is at equilibrium, and the default optimizer pushes for sparsity. If you can't post a minimal working example here, try an optimizer that uses a sparsity constraint (e.g. MIOSR optimizer with target_sparsity=4).

[akaptano]

Your issue might just be the precision of model.print(). Try model.print(precision=8)

[penguinaugustus]

Your issue might just be the precision of model.print(). Try model.print(precision=8)

Thank you for you suggestions! I tried on the same model like this:

feature_names = ['S','I','R']
poly_order = 5
threshold = 0.0000001

noise_levels = [0]

for eps in noise_levels:
    print("_______________")
    print(eps)

    optimizer = ps.STLSQ(threshold=threshold)
    model = ps.SINDy(feature_names=feature_names, optimizer=optimizer,feature_library=ps.PolynomialLibrary(degree=poly_order))
    model.fit(X+ np.random.normal(scale=eps, size=X.shape),t=t_end/num_points)
    model.print(precision=8)

The result is:

_______________
0
(S)' = 0.00000000
(I)' = 0.00000000
(R)' = 0.00000000
/Users/yufu/opt/anaconda3/lib/python3.9/site-packages/pysindy/optimizers/stlsq.py:191: UserWarning: Sparsity parameter is too big (1e-07) and eliminated all coefficients
  warnings.warn(

it seems that the issue is not a matter of precision.

[penguinaugustus]

Interesting problem! It may be that most of your data is at equilibrium, and the default optimizer pushes for sparsity. If you can't post a minimal working example here, try an optimizer that uses a sparsity constraint (e.g. MIOSR optimizer with target_sparsity=4).

Thank you for your reply! I tried MIOSR optimizer:

feature_names = ['S','I','R']
poly_order = 5
threshold = 0.0000001

noise_levels = [0]

for eps in noise_levels:
    print("_______________")
    print(eps)

    optimizer = ps.MIOSR(target_sparsity=4)
    model = ps.SINDy(feature_names=feature_names, optimizer=optimizer,feature_library=ps.PolynomialLibrary(degree=poly_order))
    model.fit(X+ np.random.normal(scale=eps, size=X.shape),t=t_end/num_points)
    model.print(precision=8)

and it returned with

Restricted license - for non-production use only - expires 2024-10-28
---------------------------------------------------------------------------
GurobiError                               Traceback (most recent call last)
/var/folders/cd/5hbmg50x21dcv5w59_3cwdrm0000gn/T/ipykernel_89288/2120609143.py in <module>
     11     optimizer = ps.MIOSR(target_sparsity=4)
     12     model = ps.SINDy(feature_names=feature_names, optimizer=optimizer,feature_library=ps.PolynomialLibrary(degree=poly_order))
---> 13     model.fit(X+ np.random.normal(scale=eps, size=X.shape),t=t_end/num_points)
     14     model.print(precision=8)

~/opt/anaconda3/lib/python3.9/site-packages/pysindy/pysindy.py in fit(self, x, t, x_dot, u, multiple_trajectories, unbias, quiet, ensemble, library_ensemble, replace, n_candidates_to_drop, n_subset, n_models, ensemble_aggregator)
    412             warnings.filterwarnings(action, category=LinAlgWarning)
    413             warnings.filterwarnings(action, category=UserWarning)
--> 414             self.model.fit(x, x_dot)
    415 
    416         # New version of sklearn changes attribute name

~/opt/anaconda3/lib/python3.9/site-packages/sklearn/pipeline.py in fit(self, X, y, **fit_params)
    392             if self._final_estimator != "passthrough":
    393                 fit_params_last_step = fit_params_steps[self.steps[-1][0]]
--> 394                 self._final_estimator.fit(Xt, y, **fit_params_last_step)
    395 
    396         return self

~/opt/anaconda3/lib/python3.9/site-packages/pysindy/optimizers/sindy_optimizer.py in fit(self, x, y)
     58         )
     59 
---> 60         self.optimizer.fit(x, y)
     61         if not hasattr(self.optimizer, "coef_"):
     62             raise AttributeError("optimizer has no attribute coef_")

~/opt/anaconda3/lib/python3.9/site-packages/pysindy/optimizers/base.py in fit(self, x_, y, sample_weight, **reduce_kws)
    176         x_normed = np.asarray(x_normed)
    177 
--> 178         self._reduce(x_normed, y, **reduce_kws)
    179         self.ind_ = np.abs(self.coef_) > 1e-14
    180 

~/opt/anaconda3/lib/python3.9/site-packages/pysindy/optimizers/miosr.py in _reduce(self, x, y)
    246             self.target_sparsity is not None or self.constraint_lhs is not None
    247         ):  # Regress jointly
--> 248             coefs = self._regress(x, y, self.target_sparsity, self.initial_guess)
    249             # Remove nonzero terms due to numerical error
    250             non_active_ixs = np.argsort(np.abs(coefs))[: -int(self.target_sparsity)]

~/opt/anaconda3/lib/python3.9/site-packages/pysindy/optimizers/miosr.py in _regress(self, X, y, k, warm_start)
    227         """
    228         m, coeff_var = self._make_model(X, y, k, warm_start)
--> 229         m.optimize()
    230         return coeff_var.X
    231 

src/gurobipy/model.pxi in gurobipy.Model.optimize()

GurobiError: Model too large for size-limited license; visit https://www.gurobi.com/free-trial for a full license

I found this: https://pysindy.readthedocs.io/en/latest/examples/1_feature_overview/example.html#MIOSR It seems that MIOSR is supported by a private company and my academic free trial has limited size. Could you recommend another optimizer for me? I have tried to search for optimizer with sparsity constraints in this way: https://pysindy.readthedocs.io/en/latest/search.html?q=+sparsity+constraint+&check_keywords=yes&area=default Thank you for you time!

[akaptano]

To make your life easier (and I think MIOSR should work now), use poly_order = 2 here. If that still doesn't do anything useful, feel free to attach your full code & data and we can take a closer look at what the issue is here.

[penguinaugustus]

To make your life easier (and I think MIOSR should work now), use poly_order = 2 here. If that still doesn't do anything useful, feel free to attach your full code & data and we can take a closer look at what the issue is here.

It works now! The result looks reasonable. Thanks a lot!

_______________
0
(S)' = -0.00040158 S I
(I)' = -0.03874030 I + 0.00040139 S I
(R)' = 0.04024876 I
_______________
0.0001
(S)' = -0.00040158 S I
(I)' = -0.03874031 I + 0.00040139 S I
(R)' = 0.04024876 I
_______________
0.001
(S)' = -0.00040158 S I
(I)' = -0.03874028 I + 0.00040139 S I
(R)' = 0.04024874 I
_______________
0.01
(S)' = -0.00040159 S I
(I)' = -0.03873997 I + 0.00040140 S I
(R)' = 0.04024851 I
_______________
0.1
(S)' = -0.00040161 S I
(I)' = -0.03874822 I + 0.00040147 S I
(R)' = 0.04024934 I
_______________
1.0
(S)' = -0.00040157 S I
(I)' = -0.03863365 I + 0.00040048 S I
(R)' = 0.04023719 I
(