Python3 skglm - 'Poisson' object has no attribute 'get_lipschitz'

Question

Im working with count data and want to fit a poisson regression with a L1 norm. I have the following code which throws the error and is reproducable:

import numpy as np
import skglm
import sklearn

X = np.random.rand(100, 10)  
y = np.random.poisson(lam=2.0, size=100)  + 1 #If poisson == 0 at any point it throws an error for count data, so add 1.

# Split the data into training and test sets
X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(X, y, test_size=0.2, random_state=42)

# Create the Poisson regression model with L1 regularization
poisson_model = skglm.GeneralizedLinearEstimator(datafit=skglm.datafits.Poisson(), penalty=skglm.penalties.L1(alpha=1.0))

# Fit the model to the training data
#The error is thrown on this line.
poisson_model.fit(X_train, y_train)

# Predict on the test set
y_pred = poisson_model.predict(X_test)

# Evaluate the model
mse = sklearn.metrics.mean_squared_error(y_test, y_pred)
print(f'Mean Squared Error: {mse}')

The error I get running on skglm 0.3.1 (most recent version I think):

Cell In[58], line 15
     12 poisson_model = skglm.GeneralizedLinearEstimator(datafit=skglm.datafits.Poisson(), penalty=skglm.penalties.L1(alpha=1.0))
     14 # Fit the model to the training data
---> 15 poisson_model.fit(X_train, y_train)
     17 # Step 5: Make predictions
     18 # Predict on the test set
     19 y_pred = poisson_model.predict(X_test)

File c:...\Python\Python311\Lib\site-packages\skglm\estimators.py:252, in GeneralizedLinearEstimator.fit(self, X, y)
    249 self.datafit = self.datafit if self.datafit else Quadratic()
    250 self.solver = self.solver if self.solver else AndersonCD()
--> 252 return _glm_fit(X, y, self, self.datafit, self.penalty, self.solver)

File c:...\Python\Python311\Lib\site-packages\skglm\estimators.py:137, in _glm_fit(X, y, model, datafit, penalty, solver)
    132     if len(penalty.weights) != n_features:
    133         raise ValueError(
    134             "The size of the WeightedL1 penalty weights should be n_features, "
    135             "expected %i, got %i." % (X_.shape[1], len(penalty.weights)))
--> 137 coefs, p_obj, kkt = solver.solve(X_, y, datafit_jit, penalty_jit, w, Xw)
    138 model.coef_, model.stop_crit_ = coefs[:n_features], kkt
    139 if y.ndim == 1:

File c:...\Python\Python311\Lib\site-packages\skglm\solvers\anderson_cd.py:83, in AndersonCD.solve(self, X, y, datafit, penalty, w_init, Xw_init)
     81 else:
     82     datafit.initialize(X, y)
---> 83     lipschitz = datafit.get_lipschitz(X, y)
     85 if len(w) != n_features + self.fit_intercept:
     86     if self.fit_intercept:

AttributeError: 'Poisson' object has no attribute 'get_lipschitz'

I changed the file location names since it has my full name but otherwise it is copy pasted.

This thread on sklearn says that it should be possible to perform a poisson regression with an L1 norm, but I can't figure out what I'm doing wrong.

Any help is greatly appreciated, thank you!