Bug: `predict()` method fails with new data and pyfixest special syntax - `i()`, `C()`, `var1^var2` operators

[s3alfisc]

Examples:

%load_ext autoreload
%autoreload 2 

import pyfixest as pf 
data = pf.get_data()

#fit = pf.feols("Y ~ i(X1) | f1", data = data)
#fit.predict(newdata = data.iloc[0:100])
#FactorEvaluationError: Unable to evaluate factor `i(X1)`. [NameError: name 'i' is not defined]

fit = pf.feols("Y ~ X1 | f1^f2", data = data)
fit.predict(newdata = data.iloc[0:100])
# KeyError: "None of [Index(['f1^f2'], dtype='object')] are in the [columns]"

[s3alfisc]

Now throws more helpful error message. Needs to be moved to use internal model_matrix_fixest function and not just call formulaic.

[s3alfisc]

This should actually be closed with the recent changes @leostimpfle? At least it should work for the i() and C() operator? interacted fixed effects I am still not sure. Maybe we could add small tests in case we haven't yet?