From ff088cf7e136b9db196828fbab5b242a2f60f12e Mon Sep 17 00:00:00 2001 From: nepslor Date: Fri, 26 Apr 2024 20:49:18 +0200 Subject: [PATCH] refactoring fast adaptive models --- .../fast_adaptive_models.py | 506 +++++++++++++++++ .../forecasting_models/holtwinters.py | 537 +----------------- .../random_fourier_features.py | 227 ++++++++ tests/test_models.py | 6 +- 4 files changed, 760 insertions(+), 516 deletions(-) create mode 100644 pyforecaster/forecasting_models/fast_adaptive_models.py create mode 100644 pyforecaster/forecasting_models/random_fourier_features.py diff --git a/pyforecaster/forecasting_models/fast_adaptive_models.py b/pyforecaster/forecasting_models/fast_adaptive_models.py new file mode 100644 index 0000000..7d80122 --- /dev/null +++ b/pyforecaster/forecasting_models/fast_adaptive_models.py @@ -0,0 +1,506 @@ +import numpy as np +import pandas as pd +from pyforecaster.forecaster import ScenarioGenerator +from pyforecaster.utilities import kalman +from pyforecaster.forecasting_models.holtwinters import tune_hyperpars, hankel + +def get_basis(t, l, n_h, frequencies=None): + """ + Get the first n_h sine and cosine basis functions and the projection + matrix P + """ + frequencies = np.arange(n_h)+1 if frequencies is None else frequencies + trigonometric = np.vstack([np.vstack([np.sin(2*np.pi*k*t/l), np.cos(2*np.pi*k*t/l)]) for k in frequencies]) + P = np.vstack([np.ones(len(t))/np.sqrt(2), trigonometric]).T * np.sqrt(2 / l) + return P + +class Fourier_es(ScenarioGenerator): + + def __init__(self, target_name='target', targets_names=None, n_sa=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, + val_ratio=0.8, nodes_at_step=None, q_vect=None, periodicity=None, optimize_hyperpars=True, + optimization_budget=100, diagnostic_plots=True, verbose=True, **scengen_kwgs): + """ + :param y: + :param h: + :param alpha: + :param m: + :return: + """ + assert m>0, 'm must be positive' + assert 00, 'n_harmonics must be positive' + + self.init_pars = {'target_name': target_name, 'n_sa': n_sa, 'alpha': alpha, 'm': m, 'omega': omega, + 'n_harmonics': n_harmonics, 'val_ratio': val_ratio, 'nodes_at_step': nodes_at_step, + 'q_vect': q_vect, 'periodicity': periodicity, 'optimize_hyperpars': optimize_hyperpars, + 'optimization_budget': optimization_budget, 'diagnostic_plots': diagnostic_plots, + 'verbose': verbose} + self.init_pars.update(scengen_kwgs) + self.targets_names = [target_name] if targets_names is None else targets_names + self.optimize_hyperpars = optimize_hyperpars + self.optimization_budget = optimization_budget + self.periodicity = periodicity if periodicity is not None else n_sa + self.target_name = target_name + self.verbose = verbose + # precompute basis over all possible periods + self.alpha = alpha + self.omega = omega + self.n_sa=n_sa + self.m = m + self.n_harmonics = n_harmonics + self.n_harmonics_int = int(np.minimum(n_harmonics, m // 2)) + self.coeffs = None + self.eps = None + self.last_1sa_preds = 0 + self.w = np.zeros(m) + self.P_past = None + self.P_future = None + self.reinint_pars() + self.diagnostic_plots = diagnostic_plots + super().__init__(q_vect, nodes_at_step=nodes_at_step, val_ratio=val_ratio, **scengen_kwgs) + + def reinint_pars(self): + self.store_basis() + + def store_basis(self): + self.n_harmonics_int = int(np.minimum(self.n_harmonics, self.m // 2)) + t_f = np.arange(2 * self.m + np.maximum(self.n_sa, self.periodicity)) + self.P_future = get_basis(t_f, self.m, self.n_harmonics_int) + + def fit(self, x_pd, y_pd=None, **kwargs): + if self.optimize_hyperpars: + pars_opt = tune_hyperpars(x_pd, Fourier_es, hyperpars={'alpha':[0, 1], 'omega':[0, 1], 'n_harmonics':[1, self.m//2]}, + n_trials=self.optimization_budget, targets_names=self.targets_names, return_model=False, + **self.init_pars) + self.set_params(**pars_opt) + self.reinint_pars() + + y_present = x_pd[self.target_name].values + x = x_pd.values + + + # exclude the last n_sa, we need them to create the target + preds = self.run(x, y_present, start_from=0, fit=True)[:-self.n_sa] + #preds = self.predict(x_pd)[:-self.n_sa] + + # hankelize the target + hw_target = hankel(y_present[1:], self.n_sa) + resid = hw_target - preds + self.err_distr = np.quantile(resid, self.q_vect, axis=0).T + return self + + def predict(self, x_pd, **kwargs): + x = x_pd.values + y = x_pd[self.target_name].values + return self.run(x, y, start_from=0, fit=False) + + + def run(self, x, y, return_coeffs=False, start_from=0, fit=True): + if self.coeffs is None: + coeffs = np.zeros(2 * self.n_harmonics_int + 1) + coeffs[0] = y[0]*np.sqrt(self.m) + eps = 0 + preds = [[0]] + else: + eps = self.eps + coeffs = self.coeffs + preds = [[y[start_from]+eps]] + + coeffs_t_history = [] + last_1sa_preds = np.copy(self.last_1sa_preds) + w = np.copy(self.w) + for i in range(start_from, len(y)): + P_past = self.P_future[i % self.m:(i % self.m + self.m), :] + # this is copy-pasting the Fourier smoothing in the last periodicity + P_f = self.P_future[i % self.m + self.m - self.periodicity:i % self.m + self.m -self.periodicity + self.n_sa, :] + start = np.maximum(0, i-self.m+1) + + # deal with the prediction case: no y, we use stored window values + y_w = y[start:i+1] + if len(y_w)==self.m: + w = y_w + else: + w = np.roll(np.copy(self.w), -len(y_w)) + w[-len(y_w):] = y_w + + #eps = self.omega * (w[-1] - preds[-1][0]) + (1-self.omega) * eps + eps_obs = w[-1] - last_1sa_preds + eps = self.omega * eps_obs + (1-self.omega) * eps + coeffs_t = P_past[-len(w):,:].T@w + coeffs = self.alpha*coeffs_t + (1-self.alpha)*coeffs + last_preds = (P_f@coeffs).ravel() + last_1sa_preds = last_preds[0] + preds.append((last_preds + eps*(self.n_sa-np.arange(self.n_sa))**2/self.n_sa**2).ravel() ) + if return_coeffs: + coeffs_t_history.append(coeffs) + + # only store states if we are fitting + if fit: + self.coeffs = coeffs + self.eps = eps + self.last_1sa_preds = last_1sa_preds + self.w = w + + if return_coeffs: + return np.vstack(preds[1:]), np.vstack(coeffs_t_history) + return np.vstack(preds[1:]) + + def predict_quantiles(self, x, **kwargs): + preds = self.predict(x) + return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) + + def __getstate__(self): + return (self.coeffs, self.eps, self.last_1sa_preds, self.w) + def __setstate__(self, state): + self.coeffs, self.eps, self.last_1sa_preds, self.w = state + + + +#@njit +def update_predictions(coeffs_t_history, start_from, y, Ps_future, period, h, m, omega, last_1sa_preds, eps, x, P, F, Q, H, R, n_harmonics, w_init): + preds = [] + for i in range(start_from, len(y)): + P_past = Ps_future[i % m:(i % m + m), :] + # this is copy-pasting the Fourier smoothing in the last periodicity + P_f = Ps_future[i % m + m -period: i % m + m -period+h, :] + + start = max(0, i - m + 1) + + # deal with the prediction case: no y, we use stored window values + y_w = y[start:i + 1].copy() + if len(y_w) == m: + w = y_w + else: + w = np.roll(w_init, -len(y_w)) + w[-len(y_w):] = y_w + + + coeffs_t = np.dot(P_past[-len(w):, :].T, w) + x, P = kalman(x, P, F, Q, coeffs_t, H, R) # Assuming kalman is a Numba-compatible function + coeffs = x + + last_preds = np.dot(P_f, coeffs).ravel() + last_1sa_preds = last_preds[0].copy() + + eps = omega * (w[-1] - last_1sa_preds) + (1 - omega) * eps + preds.append((last_preds + eps * (h - np.arange(h)) ** 2 / h ** 2).ravel()) + coeffs_t_history[:, i] = coeffs + + if i % m == 0 and i >= m: + Q = np.corrcoef(coeffs_t_history[:, :i]) + R = Q * 0.01 + #P = np.eye(n_harmonics * 2 + 1) * 1000 + + return preds, last_1sa_preds, eps, x, P, Q, R, coeffs_t_history, w + + + +class FK(ScenarioGenerator): + + + def __init__(self, target_name='target', targets_names=None, n_sa=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, val_ratio=0.8, nodes_at_step=None, q_vect=None, periodicity=None, + optimize_hyperpars=True, optimization_budget=100, r=0.1, q=0.1, verbose=True, **scengen_kwgs): + """ + :param y: + :param h: + :param alpha: + :param m: + :return: + """ + assert m>0, 'm must be positive' + assert 00, 'n_harmonics must be positive' + + + self.targets_names = [target_name] if targets_names is None else targets_names + self.init_pars = {'target_name': target_name, 'n_sa': n_sa, 'alpha': alpha, 'm': m, 'omega': omega, + 'n_harmonics': n_harmonics, 'val_ratio': val_ratio, 'nodes_at_step': nodes_at_step, + 'q_vect': q_vect, 'periodicity': periodicity, 'optimize_hyperpars': optimize_hyperpars, + 'optimization_budget': optimization_budget, 'r': r, 'q': q, 'verbose':verbose} + self.init_pars.update(scengen_kwgs) + self.verbose=verbose + self.optimize_hyperpars = optimize_hyperpars + self.optimization_budget = optimization_budget + self.periodicity = periodicity if periodicity is not None else n_sa + if self.periodicity < n_sa: + print('WARNING: periodicity is smaller than n_sa, this may lead to suboptimal results.') + + self.target_name = target_name + # precompute basis over all possible periods + self.alpha = alpha + self.omega = omega + self.n_sa=n_sa + self.m = m + + self.eps = None + self.last_1sa_preds = 0 + self.w = np.zeros(m) + + self.n_harmonics = n_harmonics + self.r = r + self.q = q + + self.reinit_pars() + self.coeffs_t_history = [] + super().__init__(q_vect, nodes_at_step=nodes_at_step, val_ratio=val_ratio, **scengen_kwgs) + + def reinit_pars(self): + self.n_harmonics_int = int(np.minimum(self.n_harmonics, self.m // 2)) + # precompute basis over all possible periods + self.x = np.zeros(2 * self.n_harmonics_int + 1) + + self.F = np.eye(self.n_harmonics_int * 2 + 1) + + self.H = np.eye(self.n_harmonics_int * 2 + 1) + self.P = np.eye(self.n_harmonics_int * 2 + 1) * 1000 + + self.R = np.eye(self.n_harmonics_int * 2 + 1) * self.r + self.Q = np.eye(self.n_harmonics_int * 2 + 1) * self.q + self.P_future = None + self.store_basis() + + def store_basis(self): + t_f = np.arange(2 * self.m + np.maximum(self.n_sa, self.periodicity)) + self.P_future = get_basis(t_f, self.m, self.n_harmonics_int) + + def fit(self, x_pd, y_pd=None, **kwargs): + if self.optimize_hyperpars: + pars_opt = tune_hyperpars(x_pd, FK, hyperpars={'alpha':[0, 1], 'omega':[0, 1], + 'n_harmonics':[1, self.m//2], + 'r':[0.001, 0.8], 'q':[0.001, 0.8]}, + n_trials=self.optimization_budget, targets_names=self.targets_names, return_model=False, + **self.init_pars) + self.set_params(**pars_opt) + self.reinit_pars() + + + y_present = x_pd[self.target_name].values + x = x_pd.values + + # exclude the last n_sa, we need them to create the target + preds = self.run(x, y_present, start_from=0, fit=True)[:-self.n_sa] + #preds = self.predict(x_pd)[:-self.n_sa] + + # hankelize the target + hw_target = hankel(y_present[1:], self.n_sa) + resid = hw_target - preds + self.err_distr = np.quantile(resid, self.q_vect, axis=0).T + return self + + def predict(self, x_pd, **kwargs): + x = x_pd.values + y = x_pd[self.target_name].values + return self.run(x, y, start_from=0, fit=False) + + + def run(self, x, y, return_coeffs=False, start_from=0, fit=True): + if self.eps is None: + eps = 0 + preds = [[0]] + else: + eps = self.eps + preds = [[y[start_from]+eps]] + + if len(self.coeffs_t_history)>0: + coeffs_t_history = np.vstack([self.coeffs_t_history, np.zeros((len(y) - start_from, self.n_harmonics_int * 2 + 1))]) + else: + coeffs_t_history = np.zeros((len(y) - start_from, self.n_harmonics_int * 2 + 1)) + w_init = np.copy(self.w) + preds_updated, last_1sa_preds, eps, x, P, Q, R, coeffs_t_history, w = update_predictions(coeffs_t_history.T, start_from, y, self.P_future, self.periodicity, self.n_sa, self.m, self.omega, self.last_1sa_preds, eps, self.x, self.P, self.F, self.Q, self.H, self.R, + self.n_harmonics_int, w_init) + if fit: + self.last_1sa_preds = last_1sa_preds + self.eps = eps + self.x = x + self.P = P + self.Q = Q + self.R = R + self.w = w + + preds = preds + preds_updated + self.coeffs_t_history = coeffs_t_history.T + + if return_coeffs: + return np.vstack(preds[1:]), coeffs_t_history.T + return np.vstack(preds[1:]) + + def predict_quantiles(self, x, **kwargs): + preds = self.predict(x) + return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) + + def __getstate__(self): + return (self.eps,self.last_1sa_preds, self.x, self.P, self.R, self.w) + def __setstate__(self, state): + self.eps, self.last_1sa_preds, self.x, self.P, self.R, self.w = state + + +class FK_multi(ScenarioGenerator): + """ + Multistep ahead forecasting with multiple Fourier-Kalman regressors + """ + + def __init__(self, target_name='target', targets_names=None, n_sa=1, n_predictors=4, alphas=None, m=24, omegas=None, + ns_harmonics=None, val_ratio=0.8, nodes_at_step=None, q_vect=None, periodicity=None, + base_predictor=Fourier_es, optimize_hyperpars=False, optimize_submodels_hyperpars=True, + optimization_budget=100, r=0.1, q=0.1, verbose=True, **scengen_kwgs): + """ + :param y: + :param h: this is the numebr of steps ahead to be predicted. + :param alpha: + :param m: + :return: + """ + n_predictors = int(n_predictors) + self.base_predictor = base_predictor + self.targets_names = targets_names + self.init_pars = {'target_name': target_name, 'n_sa': n_sa, 'n_predictors': n_predictors, 'alpha': alphas, 'm': m, 'omegas': omegas, + 'ns_harmonics': ns_harmonics, 'val_ratio': val_ratio, 'nodes_at_step': nodes_at_step, + 'q_vect': q_vect, 'periodicity': periodicity, 'optimize_hyperpars': optimize_hyperpars, + 'optimization_budget': optimization_budget, 'r': r, 'q': q,'optimize_submodels_hyperpars':optimize_submodels_hyperpars, + 'verbose':verbose, 'base_predictor': base_predictor,'targets_names':targets_names} + self.init_pars.update(scengen_kwgs) + self.verbose=verbose + self.optimize_hyperpars = optimize_hyperpars + self.optimization_budget = optimization_budget + self.optimize_submodels_hyperpars = optimize_submodels_hyperpars + + self.periodicity = periodicity if periodicity is not None else n_sa + assert self.periodicity < m, 'Periodicity must be smaller than history m' + if self.periodicity < n_sa: + print('WARNING: periodicity is smaller than n_sa, this may lead to suboptimal results.') + self.alphas = np.ones(n_predictors) * 0.01 if alphas is None else alphas + self.omegas = np.ones(n_predictors) * 0.01 if omegas is None else omegas + self.n_sa = n_sa + self.m = m + self.target_name = target_name + self.ns_harmonics = (np.ones(n_predictors) * 10).astype(int) if ns_harmonics is None else ns_harmonics + self.n_predictors = n_predictors + self.r = r + self.q = q + self.reinit_pars() + self.coeffs_t_history = [] + + super().__init__(q_vect, nodes_at_step=nodes_at_step, val_ratio=val_ratio, **scengen_kwgs) + + def reinit_pars(self): + ms = np.linspace(1, self.m, self.n_predictors + 1).astype(int)[1:] + ms = np.maximum(ms, self.n_sa) + + if np.any([m= self.n_sa: + # average last point error over different prediction times for all the models + last_obs = x_pd[self.target_name].iloc[i] + prev_obs = x_pd[self.target_name].iloc[i-1] + # sometimes persistence error is too small, sometimes signal could be small. We normalize by the average + norm_avg_err = [np.mean([np.abs(preds_experts[i-sa, sa, predictor]-last_obs) for sa in range(self.n_sa)]) for predictor in range(self.n_predictors)] + norm_avg_err = np.array(norm_avg_err) / np.mean(norm_avg_err) + coeffs_t = np.exp(-np.array(norm_avg_err)) / np.exp(-np.array(norm_avg_err)).sum() + else: + coeffs_t = x + try: + x, P = kalman(x, P, self.F, Q, coeffs_t, self.H, R) + except Exception as e: + print(e) + if not np.any(np.isnan(x)): + coeffs = x + else: + coeffs = coeffs_t + + coeffs = np.abs(coeffs / np.abs(coeffs).sum()) + preds.append(preds_experts[i] @ coeffs) + + if return_coeffs: + coeffs_t_history.append(coeffs_t) + if i % self.m == 0 and i >= self.m: + Q = np.corrcoef(np.vstack(coeffs_t_history).T) + R = np.corrcoef(np.vstack(coeffs_t_history).T) * 10 + #P = np.eye(self.n_predictors) * 1000 + + if fit: + # if we were fitting, store states. Do nothing if we're predicting + self.states = [self.models[j].__getstate__() for j in range(self.n_predictors)] + self.Q = Q + self.R = R + self.P = P + self.x = x + + if return_coeffs: + return np.vstack(preds[1:]), np.vstack(coeffs_t_history) + return np.vstack(preds[1:]) + + def predict_quantiles(self, x, **kwargs): + preds = self.predict(x) + return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) + + def __getstate__(self): + return (self.coeffs, self.eps,self.last_1sa_preds) + def __setstate__(self, state): + self.coeffs, self.eps, self.last_1sa_preds = state \ No newline at end of file diff --git a/pyforecaster/forecasting_models/holtwinters.py b/pyforecaster/forecasting_models/holtwinters.py index 3669221..b8219ea 100644 --- a/pyforecaster/forecasting_models/holtwinters.py +++ b/pyforecaster/forecasting_models/holtwinters.py @@ -23,14 +23,15 @@ def fit_sample(pars_dict,model_class, model_init_kwargs, x, return_model=False): targets_names = model_init_kwargs['targets_names'] if 'targets_names' in model_init_kwargs.keys() else None model_init_kwargs.update(pars_dict) model = model_class(**model_init_kwargs) - if targets_names is not None and model_class not in [HoltWintersMulti, FK_multi]: + if targets_names is not None and model_class.__name__ not in ["HoltWintersMulti", "FK_multi"]: subscores = [] for c in targets_names: model_init_kwargs.update({'target_name': c, 'targets_names': None}) model = model_class(**model_init_kwargs) # features are all x columns except the targets and current target c features_names = [n for n in x.columns if n not in set(targets_names) - {c}] - subscores.append(score_autoregressive(model, x[features_names], target_name=c)) + model.target_name = c + subscores.append(score_autoregressive(model, x[features_names], target_name=c, n_sa=model.n_sa)) score = np.mean(subscores) else: model.fit(x) @@ -68,25 +69,32 @@ def tune_hyperpars(x, model_class, hyperpars, n_trials=100, verbose=True, return if len(pars_cartridge)>1: - from concurrent.futures import ProcessPoolExecutor - from multiprocessing import cpu_count - with ProcessPoolExecutor(max_workers=np.minimum(cpu_count(), 10)) as executor: - scores = list(tqdm(executor.map(partial(fit_sample, model_class=model_class, - model_init_kwargs=model_init_kwargs, x=x), - pars_cartridge), total=n_trials, desc='Tuning hyperpars for {}'.format(model_class.__name__))) - if verbose: - plt.figure() - t = np.linspace(0.01, 0.99, 30) - plt.plot(np.quantile(scores, t), t) + if parallel: + from concurrent.futures import ProcessPoolExecutor + from multiprocessing import cpu_count + with ProcessPoolExecutor(max_workers=np.minimum(cpu_count(), 10)) as executor: + scores = list(tqdm(executor.map(partial(fit_sample, model_class=model_class, + model_init_kwargs=model_init_kwargs, x=x), + pars_cartridge), total=n_trials, desc='Tuning hyperpars for {}'.format(model_class.__name__))) + if verbose: + plt.figure() + t = np.linspace(0.01, 0.99, 30) + plt.plot(np.quantile(scores, t), t) + else: + scores = [] + for i in tqdm(range(n_trials), desc='Tuning hyperpars for {}'.format(model_class.__name__)): + scores.append(fit_sample(pars_cartridge[i], model_class, model_init_kwargs, x)) + + else: fitted_model = fit_sample(pars_cartridge[0], model_class, model_init_kwargs, x, return_model=True) - if model_class in [HoltWintersMulti]: + if model_class.__name__ in ["HoltWintersMulti"]: alphas = np.array([m.alpha for m in fitted_model.models]) gammas_1 = ([m.gamma_1 for m in fitted_model.models]) gammas_2 = ([m.gamma_2 for m in fitted_model.models]) best_pars = {'alphas': alphas, 'gammas_1': gammas_1, 'gammas_2': gammas_2, 'optimize_submodels_hyperpars':False} - elif model_class in [FK_multi]: + elif model_class.__name__ in ["FK_multi"]: alphas = np.array([m.alpha for m in fitted_model.models]) omegas = ([m.omega for m in fitted_model.models]) ns_harmonics = ([m.n_harmonics for m in fitted_model.models]) @@ -105,11 +113,7 @@ def tune_hyperpars(x, model_class, hyperpars, n_trials=100, verbose=True, return else: return model.get_params() -def score_autoregressive(model, x, tr_ratio=0.7, target_name=None): - target_name = model.target_name if target_name is None else target_name - model.target_name = target_name - n_sa = model.n_sa - +def score_autoregressive(model, x, tr_ratio=0.7, target_name=None, n_sa=1): # training test split n_training = int(tr_ratio*len(x)) x_tr, x_te = x.iloc[:n_training], x.iloc[n_training:] @@ -452,500 +456,5 @@ def predict_quantiles(self, x, **kwargs): return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) -def get_basis(t, l, n_h, frequencies=None): - """ - Get the first n_h sine and cosine basis functions and the projection - matrix P - """ - frequencies = np.arange(n_h)+1 if frequencies is None else frequencies - trigonometric = np.vstack([np.vstack([np.sin(2*np.pi*k*t/l), np.cos(2*np.pi*k*t/l)]) for k in frequencies]) - P = np.vstack([np.ones(len(t))/np.sqrt(2), trigonometric]).T * np.sqrt(2 / l) - return P - - -class Fourier_es(ScenarioGenerator): - - - def __init__(self, target_name='target', targets_names=None, n_sa=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, - val_ratio=0.8, nodes_at_step=None, q_vect=None, periodicity=None, optimize_hyperpars=True, - optimization_budget=100, diagnostic_plots=True, verbose=True, **scengen_kwgs): - """ - :param y: - :param h: - :param alpha: - :param m: - :return: - """ - assert m>0, 'm must be positive' - assert 00, 'n_harmonics must be positive' - - self.init_pars = {'target_name': target_name, 'n_sa': n_sa, 'alpha': alpha, 'm': m, 'omega': omega, - 'n_harmonics': n_harmonics, 'val_ratio': val_ratio, 'nodes_at_step': nodes_at_step, - 'q_vect': q_vect, 'periodicity': periodicity, 'optimize_hyperpars': optimize_hyperpars, - 'optimization_budget': optimization_budget, 'diagnostic_plots': diagnostic_plots, - 'verbose': verbose} - self.init_pars.update(scengen_kwgs) - self.targets_names = [target_name] if targets_names is None else targets_names - self.optimize_hyperpars = optimize_hyperpars - self.optimization_budget = optimization_budget - self.periodicity = periodicity if periodicity is not None else n_sa - self.target_name = target_name - self.verbose = verbose - # precompute basis over all possible periods - self.alpha = alpha - self.omega = omega - self.n_sa=n_sa - self.m = m - self.n_harmonics = n_harmonics - self.n_harmonics_int = int(np.minimum(n_harmonics, m // 2)) - self.coeffs = None - self.eps = None - self.last_1sa_preds = 0 - self.w = np.zeros(m) - self.P_past = None - self.P_future = None - self.store_basis() - self.diagnostic_plots = diagnostic_plots - super().__init__(q_vect, nodes_at_step=nodes_at_step, val_ratio=val_ratio, **scengen_kwgs) - - def store_basis(self): - t_f = np.arange(2 * self.m + np.maximum(self.n_sa, self.periodicity)) - self.P_future = get_basis(t_f, self.m, self.n_harmonics_int) - def fit(self, x_pd, y_pd=None, **kwargs): - if self.optimize_hyperpars: - pars_opt = tune_hyperpars(x_pd, Fourier_es, hyperpars={'alpha':[0, 1], 'omega':[0, 1], 'n_harmonics':[1, self.m//2]}, - n_trials=self.optimization_budget, targets_names=self.targets_names, return_model=False, - **self.init_pars) - self.set_params(**pars_opt) - self.store_basis() - - y_present = x_pd[self.target_name].values - x = x_pd.values - - - # exclude the last n_sa, we need them to create the target - preds = self.run(x, y_present, start_from=0, fit=True)[:-self.n_sa] - #preds = self.predict(x_pd)[:-self.n_sa] - - # hankelize the target - hw_target = hankel(y_present[1:], self.n_sa) - resid = hw_target - preds - self.err_distr = np.quantile(resid, self.q_vect, axis=0).T - return self - - def predict(self, x_pd, **kwargs): - x = x_pd.values - y = x_pd[self.target_name].values - return self.run(x, y, start_from=0, fit=False) - - - def run(self, x, y, return_coeffs=False, start_from=0, fit=True): - if self.coeffs is None: - coeffs = np.zeros(2 * self.n_harmonics_int + 1) - coeffs[0] = y[0]*np.sqrt(self.m) - eps = 0 - preds = [[0]] - else: - eps = self.eps - coeffs = self.coeffs - preds = [[y[start_from]+eps]] - - coeffs_t_history = [] - last_1sa_preds = np.copy(self.last_1sa_preds) - w = np.copy(self.w) - for i in range(start_from, len(y)): - P_past = self.P_future[i % self.m:(i % self.m + self.m), :] - # this is copy-pasting the Fourier smoothing in the last periodicity - P_f = self.P_future[i % self.m + self.m - self.periodicity:i % self.m + self.m -self.periodicity + self.n_sa, :] - start = np.maximum(0, i-self.m+1) - - # deal with the prediction case: no y, we use stored window values - y_w = y[start:i+1] - if len(y_w)==self.m: - w = y_w - else: - w = np.roll(np.copy(self.w), -len(y_w)) - w[-len(y_w):] = y_w - - #eps = self.omega * (w[-1] - preds[-1][0]) + (1-self.omega) * eps - eps_obs = w[-1] - last_1sa_preds - eps = self.omega * eps_obs + (1-self.omega) * eps - coeffs_t = P_past[-len(w):,:].T@w - coeffs = self.alpha*coeffs_t + (1-self.alpha)*coeffs - last_preds = (P_f@coeffs).ravel() - last_1sa_preds = last_preds[0] - preds.append((last_preds + eps*(self.n_sa-np.arange(self.n_sa))**2/self.n_sa**2).ravel() ) - if return_coeffs: - coeffs_t_history.append(coeffs) - - # only store states if we are fitting - if fit: - self.coeffs = coeffs - self.eps = eps - self.last_1sa_preds = last_1sa_preds - self.w = w - - if return_coeffs: - return np.vstack(preds[1:]), np.vstack(coeffs_t_history) - return np.vstack(preds[1:]) - - def predict_quantiles(self, x, **kwargs): - preds = self.predict(x) - return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) - - def __getstate__(self): - return (self.coeffs, self.eps, self.last_1sa_preds, self.w) - def __setstate__(self, state): - self.coeffs, self.eps, self.last_1sa_preds, self.w = state - - - -#@njit -def update_predictions(coeffs_t_history, start_from, y, Ps_future, period, h, m, omega, last_1sa_preds, eps, x, P, F, Q, H, R, n_harmonics, w_init): - preds = [] - for i in range(start_from, len(y)): - P_past = Ps_future[i % m:(i % m + m), :] - # this is copy-pasting the Fourier smoothing in the last periodicity - P_f = Ps_future[i % m + m -period: i % m + m -period+h, :] - - start = max(0, i - m + 1) - - # deal with the prediction case: no y, we use stored window values - y_w = y[start:i + 1].copy() - if len(y_w) == m: - w = y_w - else: - w = np.roll(w_init, -len(y_w)) - w[-len(y_w):] = y_w - - - coeffs_t = np.dot(P_past[-len(w):, :].T, w) - x, P = kalman(x, P, F, Q, coeffs_t, H, R) # Assuming kalman is a Numba-compatible function - coeffs = x - - last_preds = np.dot(P_f, coeffs).ravel() - last_1sa_preds = last_preds[0].copy() - - eps = omega * (w[-1] - last_1sa_preds) + (1 - omega) * eps - preds.append((last_preds + eps * (h - np.arange(h)) ** 2 / h ** 2).ravel()) - coeffs_t_history[:, i] = coeffs - - if i % m == 0 and i >= m: - Q = np.corrcoef(coeffs_t_history[:, :i]) - R = Q * 0.01 - #P = np.eye(n_harmonics * 2 + 1) * 1000 - - return preds, last_1sa_preds, eps, x, P, Q, R, coeffs_t_history, w - - - -class FK(ScenarioGenerator): - - - def __init__(self, target_name='target', targets_names=None, n_sa=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, val_ratio=0.8, nodes_at_step=None, q_vect=None, periodicity=None, - optimize_hyperpars=True, optimization_budget=100, r=0.1, q=0.1, verbose=True, **scengen_kwgs): - """ - :param y: - :param h: - :param alpha: - :param m: - :return: - """ - assert m>0, 'm must be positive' - assert 00, 'n_harmonics must be positive' - - - self.targets_names = [target_name] if targets_names is None else targets_names - self.init_pars = {'target_name': target_name, 'n_sa': n_sa, 'alpha': alpha, 'm': m, 'omega': omega, - 'n_harmonics': n_harmonics, 'val_ratio': val_ratio, 'nodes_at_step': nodes_at_step, - 'q_vect': q_vect, 'periodicity': periodicity, 'optimize_hyperpars': optimize_hyperpars, - 'optimization_budget': optimization_budget, 'r': r, 'q': q, 'verbose':verbose} - self.init_pars.update(scengen_kwgs) - self.verbose=verbose - self.optimize_hyperpars = optimize_hyperpars - self.optimization_budget = optimization_budget - self.periodicity = periodicity if periodicity is not None else n_sa - if self.periodicity < n_sa: - print('WARNING: periodicity is smaller than n_sa, this may lead to suboptimal results.') - - self.target_name = target_name - # precompute basis over all possible periods - self.alpha = alpha - self.omega = omega - self.n_sa=n_sa - self.m = m - - self.eps = None - self.last_1sa_preds = 0 - self.w = np.zeros(m) - - self.n_harmonics = n_harmonics - self.r = r - self.q = q - - self.reinit_pars() - self.coeffs_t_history = [] - super().__init__(q_vect, nodes_at_step=nodes_at_step, val_ratio=val_ratio, **scengen_kwgs) - def reinit_pars(self): - self.n_harmonics_int = int(np.minimum(self.n_harmonics, self.m // 2)) - # precompute basis over all possible periods - self.x = np.zeros(2 * self.n_harmonics_int + 1) - - self.F = np.eye(self.n_harmonics_int * 2 + 1) - - self.H = np.eye(self.n_harmonics_int * 2 + 1) - self.P = np.eye(self.n_harmonics_int * 2 + 1) * 1000 - - self.R = np.eye(self.n_harmonics_int * 2 + 1) * self.r - self.Q = np.eye(self.n_harmonics_int * 2 + 1) * self.q - self.P_future = None - self.store_basis() - - def store_basis(self): - t_f = np.arange(2 * self.m + np.maximum(self.n_sa, self.periodicity)) - self.P_future = get_basis(t_f, self.m, self.n_harmonics_int) - - def fit(self, x_pd, y_pd=None, **kwargs): - if self.optimize_hyperpars: - pars_opt = tune_hyperpars(x_pd, FK, hyperpars={'alpha':[0, 1], 'omega':[0, 1], - 'n_harmonics':[1, self.m//2], - 'r':[0.001, 0.8], 'q':[0.001, 0.8]}, - n_trials=self.optimization_budget, targets_names=self.targets_names, return_model=False, - **self.init_pars) - self.set_params(**pars_opt) - self.reinit_pars() - - - y_present = x_pd[self.target_name].values - x = x_pd.values - - # exclude the last n_sa, we need them to create the target - preds = self.run(x, y_present, start_from=0, fit=True)[:-self.n_sa] - #preds = self.predict(x_pd)[:-self.n_sa] - - # hankelize the target - hw_target = hankel(y_present[1:], self.n_sa) - resid = hw_target - preds - self.err_distr = np.quantile(resid, self.q_vect, axis=0).T - return self - - def predict(self, x_pd, **kwargs): - x = x_pd.values - y = x_pd[self.target_name].values - return self.run(x, y, start_from=0, fit=False) - - - def run(self, x, y, return_coeffs=False, start_from=0, fit=True): - if self.eps is None: - eps = 0 - preds = [[0]] - else: - eps = self.eps - preds = [[y[start_from]+eps]] - - if len(self.coeffs_t_history)>0: - coeffs_t_history = np.vstack([self.coeffs_t_history, np.zeros((len(y) - start_from, self.n_harmonics_int * 2 + 1))]) - else: - coeffs_t_history = np.zeros((len(y) - start_from, self.n_harmonics_int * 2 + 1)) - w_init = np.copy(self.w) - preds_updated, last_1sa_preds, eps, x, P, Q, R, coeffs_t_history, w = update_predictions(coeffs_t_history.T, start_from, y, self.P_future, self.periodicity, self.n_sa, self.m, self.omega, self.last_1sa_preds, eps, self.x, self.P, self.F, self.Q, self.H, self.R, - self.n_harmonics_int, w_init) - if fit: - self.last_1sa_preds = last_1sa_preds - self.eps = eps - self.x = x - self.P = P - self.Q = Q - self.R = R - self.w = w - - preds = preds + preds_updated - self.coeffs_t_history = coeffs_t_history.T - - if return_coeffs: - return np.vstack(preds[1:]), coeffs_t_history.T - return np.vstack(preds[1:]) - - def predict_quantiles(self, x, **kwargs): - preds = self.predict(x) - return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) - - def __getstate__(self): - return (self.eps,self.last_1sa_preds, self.x, self.P, self.R, self.w) - def __setstate__(self, state): - self.eps, self.last_1sa_preds, self.x, self.P, self.R, self.w = state - - -class FK_multi(ScenarioGenerator): - """ - Multistep ahead forecasting with multiple Fourier-Kalman regressors - """ - - def __init__(self, target_name='target', targets_names=None, n_sa=1, n_predictors=4, alphas=None, m=24, omegas=None, - ns_harmonics=None, val_ratio=0.8, nodes_at_step=None, q_vect=None, periodicity=None, - base_predictor=Fourier_es, optimize_hyperpars=False, optimize_submodels_hyperpars=True, - optimization_budget=100, r=0.1, q=0.1, verbose=True, **scengen_kwgs): - """ - :param y: - :param h: this is the numebr of steps ahead to be predicted. - :param alpha: - :param m: - :return: - """ - n_predictors = int(n_predictors) - self.base_predictor = base_predictor - self.targets_names = targets_names - self.init_pars = {'target_name': target_name, 'n_sa': n_sa, 'n_predictors': n_predictors, 'alpha': alphas, 'm': m, 'omegas': omegas, - 'ns_harmonics': ns_harmonics, 'val_ratio': val_ratio, 'nodes_at_step': nodes_at_step, - 'q_vect': q_vect, 'periodicity': periodicity, 'optimize_hyperpars': optimize_hyperpars, - 'optimization_budget': optimization_budget, 'r': r, 'q': q,'optimize_submodels_hyperpars':optimize_submodels_hyperpars, - 'verbose':verbose, 'base_predictor': base_predictor,'targets_names':targets_names} - self.init_pars.update(scengen_kwgs) - self.verbose=verbose - self.optimize_hyperpars = optimize_hyperpars - self.optimization_budget = optimization_budget - self.optimize_submodels_hyperpars = optimize_submodels_hyperpars - - self.periodicity = periodicity if periodicity is not None else n_sa - assert self.periodicity < m, 'Periodicity must be smaller than history m' - if self.periodicity < n_sa: - print('WARNING: periodicity is smaller than n_sa, this may lead to suboptimal results.') - self.alphas = np.ones(n_predictors) * 0.01 if alphas is None else alphas - self.omegas = np.ones(n_predictors) * 0.01 if omegas is None else omegas - self.n_sa = n_sa - self.m = m - self.target_name = target_name - self.ns_harmonics = (np.ones(n_predictors) * 10).astype(int) if ns_harmonics is None else ns_harmonics - self.n_predictors = n_predictors - self.r = r - self.q = q - self.reinit_pars() - self.coeffs_t_history = [] - - super().__init__(q_vect, nodes_at_step=nodes_at_step, val_ratio=val_ratio, **scengen_kwgs) - - def reinit_pars(self): - ms = np.linspace(1, self.m, self.n_predictors + 1).astype(int)[1:] - ms = np.maximum(ms, self.n_sa) - - if np.any([m= self.n_sa: - # average last point error over different prediction times for all the models - last_obs = x_pd[self.target_name].iloc[i] - prev_obs = x_pd[self.target_name].iloc[i-1] - # sometimes persistence error is too small, sometimes signal could be small. We normalize by the average - norm_avg_err = [np.mean([np.abs(preds_experts[i-sa, sa, predictor]-last_obs) for sa in range(self.n_sa)]) for predictor in range(self.n_predictors)] - norm_avg_err = np.array(norm_avg_err) / np.mean(norm_avg_err) - coeffs_t = np.exp(-np.array(norm_avg_err)) / np.exp(-np.array(norm_avg_err)).sum() - else: - coeffs_t = x - try: - x, P = kalman(x, P, self.F, Q, coeffs_t, self.H, R) - except Exception as e: - print(e) - if not np.any(np.isnan(x)): - coeffs = x - else: - coeffs = coeffs_t - - coeffs = np.abs(coeffs / np.abs(coeffs).sum()) - preds.append(preds_experts[i] @ coeffs) - - if return_coeffs: - coeffs_t_history.append(coeffs_t) - if i % self.m == 0 and i >= self.m: - Q = np.corrcoef(np.vstack(coeffs_t_history).T) - R = np.corrcoef(np.vstack(coeffs_t_history).T) * 10 - #P = np.eye(self.n_predictors) * 1000 - - if fit: - # if we were fitting, store states. Do nothing if we're predicting - self.states = [self.models[j].__getstate__() for j in range(self.n_predictors)] - self.Q = Q - self.R = R - self.P = P - self.x = x - - if return_coeffs: - return np.vstack(preds[1:]), np.vstack(coeffs_t_history) - return np.vstack(preds[1:]) - - def predict_quantiles(self, x, **kwargs): - preds = self.predict(x) - return np.expand_dims(preds, -1) + np.expand_dims(self.err_distr, 0) - def __getstate__(self): - return (self.coeffs, self.eps,self.last_1sa_preds) - def __setstate__(self, state): - self.coeffs, self.eps, self.last_1sa_preds = state \ No newline at end of file diff --git a/pyforecaster/forecasting_models/random_fourier_features.py b/pyforecaster/forecasting_models/random_fourier_features.py new file mode 100644 index 0000000..bc67340 --- /dev/null +++ b/pyforecaster/forecasting_models/random_fourier_features.py @@ -0,0 +1,227 @@ +# MOST OF THIS CODE IS TAKEN FROM https://github.com/tiskw/random-fourier-features/blob/main/rfflearn/cpu/rfflearn_cpu_regression.py +# Author: Tetsuya Ishikawa + +import functools +import numpy as np +import scipy.stats +from sklearn.linear_model import LinearRegression + +def seed(seed): + """ + Fix random seed used in this script. + + Args: + seed (int): Random seed. + """ + # Now it is enough to fix the random seed of Numpy. + np.random.seed(seed) + + +def get_rff_matrix(dim_in, dim_out, std): + """ + Generates random matrix of random Fourier features. + + Args: + dim_in (int) : Input dimension of the random matrix. + dim_out (int) : Output dimension of the random matrix. + std (float): Standard deviation of the random matrix. + + Returns: + (np.ndarray): Random matrix with shape (dim_out, dim_in). + """ + return std * np.random.randn(dim_in, dim_out) + + +def get_orf_matrix(dim_in, dim_out, std): + """ + Generates random matrix of orthogonal random features. + + Args: + dim_in (int) : Input dimension of the random matrix. + dim_out (int) : Output dimension of the random matrix. + std (float): Standard deviation of the random matrix. + + Returns: + (np.ndarray): Random matrix with shape (dim_out, dim_in). + """ + # Initialize matrix W. + W = None + + for _ in range(dim_out // dim_in + 1): + s = scipy.stats.chi.rvs(df = dim_in, size = (dim_in, )) + Q = np.linalg.qr(np.random.randn(dim_in, dim_in))[0] + V = std * np.dot(np.diag(s), Q) + W = V if W is None else np.concatenate([W, V], axis = 1) + + # Trim unnecessary part. + return W[:dim_in, :dim_out] + + +def get_matrix_generator(rand_type, std, dim_kernel): + """ + This function returns a function which generate RFF/ORF matrix. + The usage of the returned value of this function are: + f(dim_input:int) -> np.array with shape (dim_input, dim_kernel) + """ + if rand_type == "rff": return functools.partial(get_rff_matrix, std = std, dim_out = dim_kernel) + elif rand_type == "orf": return functools.partial(get_orf_matrix, std = std, dim_out = dim_kernel) + else : raise RuntimeError("matrix_generator: 'rand_type' must be 'rff', 'orf', or 'qrf'.") + + +class Base: + """ + Base class of the following RFF/ORF related classes. + """ + def __init__(self, rand_type, dim_kernel, std_kernel, W, b): + """ + Constractor of the Base class. + Create random matrix generator and random matrix instance. + + Args: + rand_type (str) : Type of random matrix ("rff", "orf", "qrf", etc). + dim_kernel (int) : Dimension of the random matrix. + std_kernel (float) : Standard deviation of the random matrix. + W (np.ndarray): Random matrix for the input `X`. If None then generated automatically. + b (np.ndarray): Random bias for the input `X`. If None then generated automatically. + + Notes: + If `W` is None then the appropriate matrix will be set just before the training. + """ + self.dim = dim_kernel + self.s_k = std_kernel + self.mat = get_matrix_generator(rand_type, std_kernel, dim_kernel) + self.W = W + self.b = b + + def conv(self, X, index=None): + """ + Applies random matrix to the given input vectors `X` and create feature vectors. + + Args: + X (np.ndarray): Input matrix with shape (n_samples, n_features). + index (int) : Index of the random matrix. This value should be specified only + when multiple random matrices are used. + + Notes: + This function can manipulate multiple random matrix. If argument 'index' is given, + then use self.W[index] as a random matrix, otherwise use self.W itself. + Also, computation of `ts` is equivarent with ts = X @ W, however, for reducing memory + consumption, split X to smaller matrices and concatenate after multiplication wit W. + """ + W = self.W if index is None else self.W[index] + b = self.b if index is None else self.b[index] + return np.cos(X @ W + b) + + def set_weight(self, dim_in): + """ + Set the appropriate random matrix to 'self.W' if 'self.W' is None (i.e. empty). + + Args: + dim_in (int): Input dimension of the random matrix. + + Notes: + This function can manipulate multiple random matrix. If argument 'dim_in' is + a list/tuple of integers, then generate multiple random matrixes. + """ + # Generate matrix W. + if self.W is not None : pass + elif hasattr(dim_in, "__iter__"): self.W = tuple([self.mat(d) for d in dim_in]) + else : self.W = self.mat(dim_in) + + # Generate vector b. + if self.b is not None : pass + elif hasattr(dim_in, "__iter__"): self.b = tuple([np.random.uniform(0, 2*np.pi, size=self.W.shape[1]) for _ in dim_in]) + else : self.b = np.random.uniform(0, 2*np.pi, size=self.W.shape[1]) + + + +import sklearn + + +class Regression(Base): + """ + Regression with random matrix (RFF/ORF). + """ + def __init__(self, rand_type, dim_kernel=16, std_kernel=0.1, W=None, b=None, **args): + """ + Constractor. Save hyper parameters as member variables and create LinearRegression instance. + + Args: + rand_type (str) : Type of random matrix ("rff", "orf", "qrf", etc). + dim_kernel (int) : Dimension of the random matrix. + std_kernel (float) : Standard deviation of the random matrix. + W (np.ndarray): Random matrix for the input `X`. If None then generated automatically. + b (np.ndarray): Random bias for the input `X`. If None then generated automatically. + args (dict) : Extra arguments. This arguments will be passed to the constructor of sklearn's LinearRegression model. + """ + super().__init__(rand_type, dim_kernel, std_kernel, W, b) + self.reg = LinearRegression(**args) + + def fit(self, X, y, **args): + """ + Trains the RFF regression model according to the given data. + + Args: + X (np.ndarray): Input matrix with shape (n_samples, n_features_input). + y (np.ndarray): Output vector with shape (n_samples,). + args (dict) : Extra arguments. This arguments will be passed to the sklearn's `fit` function. + + Returns: + (rfflearn.cpu.Regression): Fitted estimator. + """ + self.set_weight(X.shape[1]) + self.reg.fit(self.conv(X), y, **args) + return self + + def predict(self, X, **args): + """ + Performs prediction on the given data. + + Args: + X (np.ndarray): Input matrix with shape (n_samples, n_features_input). + args (dict) : Extra arguments. This arguments will be passed to the sklearn's `predict` function. + + Returns: + (np.ndarray): Predicted vector. + """ + self.set_weight(X.shape[1]) + return self.reg.predict(self.conv(X), **args) + + def score(self, X, y, **args): + """ + Returns the R2 score (coefficient of determination) of the prediction. + + Args: + X (np.ndarray): Input matrix with shape (n_samples, n_features_input). + y (np.ndarray): Output vector with shape (n_samples,). + args (dict) : Extra arguments. This arguments will be passed to sklearn's `score` function. + + Returns: + (float): R2 score of the prediction. + """ + self.set_weight(X.shape[1]) + return self.reg.score(self.conv(X), y, **args) + + +# The above functions/classes are not visible from users of this library, becasue the usage of +# the function is a bit complicated. The following classes are simplified version of the above +# classes. The following classes are visible from users. + + +class RFFRegression(Regression): + """ + Regression with RFF. + """ + def __init__(self, *pargs, **kwargs): + super().__init__("rff", *pargs, **kwargs) + + +class ORFRegression(Regression): + """ + Regression with ORF. + """ + def __init__(self, *pargs, **kwargs): + super().__init__("orf", *pargs, **kwargs) + + + diff --git a/tests/test_models.py b/tests/test_models.py index 050acc5..afaf2f0 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -4,7 +4,8 @@ import pandas as pd import numpy as np import logging -from pyforecaster.forecasting_models.holtwinters import HoltWinters, HoltWintersMulti, FK, FK_multi, Fourier_es +from pyforecaster.forecasting_models.holtwinters import HoltWinters, HoltWintersMulti +from pyforecaster.forecasting_models.fast_adaptive_models import Fourier_es, FK, FK_multi from pyforecaster.forecasting_models.randomforests import QRF from pyforecaster.forecaster import LinearForecaster, LGBForecaster from pyforecaster.plot_utils import plot_quantiles @@ -100,6 +101,7 @@ def test_hw_multi(self): self.data = self.data.resample('1h').mean() df_tr, df_te = self.data.iloc[:1200], self.data.iloc[1200:1500] steps_day = 24 + fes = Fourier_es(n_sa=steps_day, m=steps_day*7, target_name='all', optimization_budget=20).fit(df_tr, df_tr['all']) fks_multi = FK_multi(n_predictors=4, n_sa=steps_day, m=steps_day*7, target_name='all', periodicity=steps_day*2, optimize_hyperpars=True, optimization_budget=3, targets_names=df_tr.columns[:2]).fit(df_tr) @@ -110,7 +112,7 @@ def test_hw_multi(self): fks = FK(n_sa=steps_day, m=steps_day, target_name='all', optimization_budget=2).fit(df_tr, df_tr['all']) - fes = Fourier_es(n_sa=steps_day, m=steps_day*7, target_name='all', optimization_budget=2).fit(df_tr, df_tr['all']) + hw = HoltWinters(periods=[steps_day, steps_day * 7], n_sa=steps_day, optimization_budget=2, q_vect=np.arange(11) / 10, target_name='all').fit(df_tr, df_tr['all'])