How to I build a linear additive model in Scikit?

I'm trying to build a model which takes the predication of a simplified intermediate-model f_I and multiplies it by some coefficient c_p and adds the result of some general model f_g:

The coefficients C are then optimized to fit the data.

The general model I am choosing is a RBF model, so far I have fitted the RBF to the data:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import  RBF

# The true formula
beta = np.pi
omega = 10*np.pi
alpha = 1
phi = .5*np.pi

def f(t):
    return alpha*np.tanh(beta*t)*np.sin(omega*t + phi)

# Generate some test data
t = np.linspace(0, 1, 1000)

X_train = np.random.choice(t.reshape(-1), size=12).reshape(-1, 1)
y_train = f(X_train)

kernel = RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0))

gpr = GaussianProcessRegressor(kernel=kernel, random_state=1)
gpr.fit(X_train, y_train)

fig, ax = plt.subplots(dpi=200)
ax.plot(t, f(t), 'grey', label="True")
ax.plot(X_train, y_train, 'ko', label="Training")
ax.plot(t, gpr.predict(t.reshape(-1, 1)), 'b--', label="GLM")
fig.legend()

However, I do not know how to implement this linear additive model. I have attempted to make a function for the intermediate model:

def IM(t):
    return t*np.sin(omega*t)

and then change the kernel in the GaussianProcessRegressor to something like:

kernel = IM()*ConstantKernel() + RBF(length_scale=0.5, length_scale_bounds=(0.01, 100.0))

However, this doesn't work because the IM() needs to be a class of some sort that works with the Scikit library. However I cant find much information online on how to do this.

Is making a custom kernel for the intermediate model, f_I the right approach? If so, How do I build a custom kernel to work with Scikit?

The paper that I'm following to implement this is available here, in section 2.2

from How to I build a linear additive model in Scikit?