diff --git a/.travis.yml b/.travis.yml index 5d4454a..8557bd7 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,7 +1,6 @@ dist: xenial language: python python: - - "2.7" - "3.7" install: @@ -10,3 +9,8 @@ install: script: - cd python && python setup.py develop test + - python setup.py clean + - wget https://github.com/stan-dev/cmdstan/releases/download/v2.22.1/cmdstan-2.22.1.tar.gz -O /tmp/cmdstan.tar.gz > /dev/null + - tar -xvf /tmp/cmdstan.tar.gz -C /tmp > /dev/null + - make -C /tmp/cmdstan-2.22.1/ build > /dev/null + - CMDSTAN=/tmp/cmdstan-2.22.1 STAN_BACKEND=CMDSTANPY python setup.py develop test diff --git a/python/fbprophet/diagnostics.py b/python/fbprophet/diagnostics.py index 3f37198..1461d0f 100644 --- a/python/fbprophet/diagnostics.py +++ b/python/fbprophet/diagnostics.py @@ -185,6 +185,7 @@ def prophet_copy(m, cutoff=None): mcmc_samples=m.mcmc_samples, interval_width=m.interval_width, uncertainty_samples=m.uncertainty_samples, + stan_backend=m.stan_backend.get_type() ) m2.extra_regressors = deepcopy(m.extra_regressors) m2.seasonalities = deepcopy(m.seasonalities) diff --git a/python/fbprophet/forecaster.py b/python/fbprophet/forecaster.py index 31f6648..f38bc8d 100644 --- a/python/fbprophet/forecaster.py +++ b/python/fbprophet/forecaster.py @@ -13,14 +13,10 @@ from datetime import timedelta, datetime import numpy as np import pandas as pd -import pystan # noqa F401 -from fbprophet.diagnostics import prophet_copy from fbprophet.make_holidays import get_holiday_names, make_holidays_df -from fbprophet.models import prophet_stan_model -from fbprophet.plot import (plot, plot_components, plot_forecast_component, - plot_seasonality, plot_weekly, plot_yearly, - seasonality_plot_df) +from fbprophet.models import StanBackendEnum +from fbprophet.plot import (plot, plot_components) logger = logging.getLogger('fbprophet') logger.addHandler(logging.NullHandler()) @@ -77,6 +73,9 @@ class Prophet(object): uncertainty_samples: Number of simulated draws used to estimate uncertainty intervals. Settings this value to 0 or False will disable uncertainty estimation and speed up the calculation. + uncertainty intervals. + stan_backend: str as defined in StanBackendEnum default: None - will try to + iterate over all available backends and find the working one """ def __init__( @@ -96,6 +95,7 @@ class Prophet(object): mcmc_samples=0, interval_width=0.80, uncertainty_samples=1000, + stan_backend=None ): self.growth = growth @@ -140,6 +140,20 @@ class Prophet(object): self.train_holiday_names = None self.fit_kwargs = {} self.validate_inputs() + self._load_stan_backend(stan_backend) + + def _load_stan_backend(self, stan_backend): + if stan_backend is None: + for i in StanBackendEnum: + try: + logger.debug("Trying to load backend: %s", i.name) + return self._load_stan_backend(i.name) + except Exception as e: + logger.info("Unable to load backend %s (%s), trying the next one", i.name, e) + else: + self.stan_backend = StanBackendEnum.get_backend_class(stan_backend)(logger) + + logger.info("Loaded stan backend: %s", self.stan_backend.get_type()) def validate_inputs(self): """Validates the inputs to Prophet.""" @@ -476,7 +490,7 @@ class Prophet(object): all_holidays = pd.concat((all_holidays, country_holidays_df), sort=False) all_holidays.reset_index(drop=True, inplace=True) - # Drop future holidays not previously seen in training data + # Drop future holidays not previously seen in training data if self.train_holiday_names is not None: # Remove holiday names didn't show up in fit index_to_drop = all_holidays.index[ @@ -1108,58 +1122,25 @@ class Prophet(object): dat['cap'] = history['cap_scaled'] kinit = self.logistic_growth_init(history) - model = prophet_stan_model - - def stan_init(): - return { - 'k': kinit[0], - 'm': kinit[1], - 'delta': np.zeros(len(self.changepoints_t)), - 'beta': np.zeros(seasonal_features.shape[1]), - 'sigma_obs': 1, - } + stan_init = { + 'k': kinit[0], + 'm': kinit[1], + 'delta': np.zeros(len(self.changepoints_t)), + 'beta': np.zeros(seasonal_features.shape[1]), + 'sigma_obs': 1, + } if (history['y'].min() == history['y'].max() and self.growth == 'linear'): # Nothing to fit. - self.params = stan_init() + self.params = stan_init self.params['sigma_obs'] = 1e-9 for par in self.params: self.params[par] = np.array([self.params[par]]) elif self.mcmc_samples > 0: - args = dict( - data=dat, - init=stan_init, - iter=self.mcmc_samples, - ) - args.update(kwargs) - self.stan_fit = model.sampling(**args) - for par in self.stan_fit.model_pars: - self.params[par] = self.stan_fit[par] - # Shape vector parameters - if (par in ['delta', 'beta'] - and len(self.params[par].shape) < 2): - self.params[par] = self.params[par].reshape((-1, 1)) + self.params = self.stan_backend.sampling(stan_init, dat, self.mcmc_samples, **kwargs) else: - args = dict( - data=dat, - init=stan_init, - algorithm='Newton' if dat['T'] < 100 else 'LBFGS', - iter=1e4, - ) - args.update(kwargs) - try: - self.stan_fit = model.optimizing(**args) - except RuntimeError: - logger.warning( - 'Optimization terminated abnormally. ' - 'Falling back to Newton.' - ) - args['algorithm'] = 'Newton' - self.stan_fit = model.optimizing(**args) - - for par in self.stan_fit: - self.params[par] = self.stan_fit[par].reshape((1, -1)) + self.params = self.stan_backend.fit(stan_init, dat, **kwargs) # If no changepoints were requested, replace delta with 0s if len(self.changepoints) == 0: diff --git a/python/fbprophet/models.py b/python/fbprophet/models.py index caa79f5..7b7ddd7 100644 --- a/python/fbprophet/models.py +++ b/python/fbprophet/models.py @@ -5,20 +5,269 @@ # LICENSE file in the root directory of this source tree. from __future__ import absolute_import, division, print_function - +from abc import abstractmethod, ABC +from typing import Tuple +from collections import OrderedDict +from enum import Enum import pickle - import pkg_resources +import numpy as np +import os -def get_prophet_stan_model(): - """Load compiled Stan model""" - model_file = pkg_resources.resource_filename( - 'fbprophet', - 'stan_model/prophet_model.pkl', - ) - with open(model_file, 'rb') as f: - return pickle.load(f) +class IStanBackend(ABC): + def __init__(self, logger): + self.model = self.load_model() + self.logger = logger + + @staticmethod + @abstractmethod + def get_type(): + pass + + @abstractmethod + def load_model(self): + pass + + @abstractmethod + def fit(self, stan_init, stan_data, **kwargs) -> dict: + pass + + @abstractmethod + def sampling(self, stan_init, stan_data, samples, **kwargs) -> dict: + pass + + @staticmethod + @abstractmethod + def build_model(target_dir, model_dir): + pass -prophet_stan_model = get_prophet_stan_model() +class CmdStanPyBackend(IStanBackend): + + @staticmethod + def get_type(): + return StanBackendEnum.CMDSTANPY.name + + @staticmethod + def build_model(target_dir, model_dir): + from shutil import copy + import cmdstanpy + model_name = 'prophet.stan' + target_name = 'prophet_model.bin' + + sm = cmdstanpy.Model(stan_file=os.path.join(model_dir, model_name)) + sm.compile() + copy(sm.exe_file, os.path.join(target_dir, target_name)) + + def load_model(self): + import cmdstanpy + model_file = pkg_resources.resource_filename( + 'fbprophet', + 'stan_model/prophet_model.bin', + ) + return cmdstanpy.Model(exe_file=model_file) + + def fit(self, stan_init, stan_data, **kwargs): + (stan_init, stan_data) = self.prepare_data(stan_init, stan_data) + if 'algorithm' not in kwargs: + kwargs['algorithm'] = 'Newton' if stan_data['T'] < 100 else 'LBFGS' + iterations = int(1e4) + try: + stan_fit = self.model.optimize(data=stan_data, + inits=stan_init, + iter=iterations, + **kwargs) + except RuntimeError as e: + # Fall back on Newton + if kwargs['algorithm'] != 'Newton': + self.logger.warning( + 'Optimization terminated abnormally. Falling back to Newton.' + ) + kwargs['algorithm'] = 'Newton' + stan_fit = self.model.optimize(data=stan_data, + inits=stan_init, + iter=iterations, + **kwargs) + else: + raise e + + params = self.stan_to_dict_numpy(stan_fit.column_names, stan_fit.optimized_params_np) + for par in params: + params[par] = params[par].reshape((1, -1)) + return params + + def sampling(self, stan_init, stan_data, samples, **kwargs) -> dict: + (stan_init, stan_data) = self.prepare_data(stan_init, stan_data) + + if 'chains' not in kwargs: + kwargs['chains'] = 4 + if 'warmup_iters' not in kwargs: + kwargs['warmup_iters'] = samples // 2 + + stan_fit = self.model.sample(data=stan_data, + inits=stan_init, + sampling_iters=samples, + **kwargs) + res = stan_fit.sample + (samples, c, columns) = res.shape + res = res.reshape((samples * c, columns)) + params = self.stan_to_dict_numpy(stan_fit.column_names, res) + + for par in params: + s = params[par].shape + if s[1] == 1: + params[par] = params[par].reshape((s[0],)) + + if par in ['delta', 'beta'] and len(s) < 2: + params[par] = params[par].reshape((-1, 1)) + + return params + + @staticmethod + def prepare_data(init, data) -> Tuple[dict, dict]: + cmdstanpy_data = { + 'T': data['T'], + 'S': data['S'], + 'K': data['K'], + 'tau': data['tau'], + 'trend_indicator': data['trend_indicator'], + 'y': data['y'].tolist(), + 't': data['t'].tolist(), + 'cap': data['cap'].tolist(), + 't_change': data['t_change'].tolist(), + 's_a': data['s_a'].tolist(), + 's_m': data['s_m'].tolist(), + 'X': data['X'].to_numpy().tolist(), + 'sigmas': data['sigmas'] + } + + cmdstanpy_init = { + 'k': init['k'], + 'm': init['m'], + 'delta': init['delta'].tolist(), + 'beta': init['beta'].tolist(), + 'sigma_obs': 1 + } + return (cmdstanpy_init, cmdstanpy_data) + + @staticmethod + def stan_to_dict_numpy(column_names: Tuple[str, ...], data: np.array): + output = OrderedDict() + + prev = None + + start = 0 + end = 0 + two_dims = True if len(data.shape) > 1 else False + for cname in column_names: + parsed = cname.split(".") + + curr = parsed[0] + if prev is None: + prev = curr + + if curr != prev: + if prev in output: + raise RuntimeError( + "Found repeated column name" + ) + if two_dims: + output[prev] = np.array(data[:, start:end]) + else: + output[prev] = np.array(data[start:end]) + prev = curr + start = end + end += 1 + else: + end += 1 + + if prev in output: + raise RuntimeError( + "Found repeated column name" + ) + if two_dims: + output[prev] = np.array(data[:, start:end]) + else: + output[prev] = np.array(data[start:end]) + return output + + +class PyStanBackend(IStanBackend): + + @staticmethod + def get_type(): + return StanBackendEnum.PYSTAN.name + + @staticmethod + def build_model(target_dir, model_dir): + import pystan + model_name = 'prophet.stan' + target_name = 'prophet_model.pkl' + with open(os.path.join(model_dir, model_name)) as f: + model_code = f.read() + sm = pystan.StanModel(model_code=model_code) + with open(os.path.join(target_dir, target_name), 'wb') as f: + pickle.dump(sm, f, protocol=pickle.HIGHEST_PROTOCOL) + + def sampling(self, stan_init, stan_data, samples, **kwargs) -> dict: + + args = dict( + data=stan_data, + init=lambda: stan_init, + iter=samples, + ) + args.update(kwargs) + stan_fit = self.model.sampling(**args) + out = dict() + for par in stan_fit.model_pars: + out[par] = stan_fit[par] + # Shape vector parameters + if par in ['delta', 'beta'] and len(out[par].shape) < 2: + out[par] = out[par].reshape((-1, 1)) + return out + + def fit(self, stan_init, stan_data, **kwargs) -> dict: + + args = dict( + data=stan_data, + init=lambda: stan_init, + algorithm='Newton' if stan_data['T'] < 100 else 'LBFGS', + iter=1e4, + ) + args.update(kwargs) + try: + params = self.model.optimizing(**args) + except RuntimeError: + # Fall back on Newton + self.logger.warning( + 'Optimization terminated abnormally. Falling back to Newton.' + ) + args['algorithm'] = 'Newton' + params = self.model.optimizing(**args) + + for par in params: + params[par] = params[par].reshape((1, -1)) + + return params + + def load_model(self): + """Load compiled Stan model""" + model_file = pkg_resources.resource_filename( + 'fbprophet', + 'stan_model/prophet_model.pkl', + ) + with open(model_file, 'rb') as f: + return pickle.load(f) + + +class StanBackendEnum(Enum): + PYSTAN = PyStanBackend + CMDSTANPY = CmdStanPyBackend + + @staticmethod + def get_backend_class(name: str) -> IStanBackend: + try: + return StanBackendEnum[name].value + except KeyError: + raise ValueError("Unknown stan backend: {}".format(name)) diff --git a/python/fbprophet/tests/test_prophet.py b/python/fbprophet/tests/test_prophet.py index df62131..994d755 100644 --- a/python/fbprophet/tests/test_prophet.py +++ b/python/fbprophet/tests/test_prophet.py @@ -9,13 +9,14 @@ from __future__ import print_function from __future__ import unicode_literals import os -from unittest import TestCase +import sys +from unittest import TestCase, skipUnless import numpy as np import pandas as pd - from fbprophet import Prophet + DATA = pd.read_csv( os.path.join(os.path.dirname(__file__), 'data.csv'), parse_dates=['ds'], @@ -28,14 +29,59 @@ DATA2 = pd.read_csv( class TestProphet(TestCase): + @staticmethod + def rmse(predictions, targets): + return np.sqrt(np.mean((predictions - targets) ** 2)) + def test_fit_predict(self): + days = 30 N = DATA.shape[0] - train = DATA.head(N // 2) - future = DATA.tail(N // 2) + train = DATA.head(N - days) + test = DATA.tail(days) + test.reset_index(inplace=True) forecaster = Prophet() - forecaster.fit(train) - forecaster.predict(future) + forecaster.fit(train, seed=1237861298) + np.random.seed(876543987) + future = forecaster.make_future_dataframe(days, include_history=False) + future = forecaster.predict(future) + # this gives ~ 10.64 + res = self.rmse(future['yhat'], test['y']) + self.assertTrue(15 > res > 5, msg="backend: {}".format(forecaster.stan_backend)) + + def test_fit_predict_newton(self): + days = 30 + N = DATA.shape[0] + train = DATA.head(N - days) + test = DATA.tail(days) + test.reset_index(inplace=True) + + forecaster = Prophet() + forecaster.fit(train, seed=1237861298, algorithm='Newton') + np.random.seed(876543987) + future = forecaster.make_future_dataframe(days, include_history=False) + future = forecaster.predict(future) + # this gives ~ 10.64 + res = self.rmse(future['yhat'], test['y']) + self.assertAlmostEqual(res, 23.44, places=2, msg="backend: {}".format(forecaster.stan_backend)) + + @skipUnless("--test-slow" in sys.argv, "Skipped due to the lack of '--test-slow' argument") + def test_fit_sampling_predict(self): + days = 30 + N = DATA.shape[0] + train = DATA.head(N - days) + test = DATA.tail(days) + test.reset_index(inplace=True) + + forecaster = Prophet(mcmc_samples=500) + + forecaster.fit(train, seed=1237861298, chains=4) + np.random.seed(876543987) + future = forecaster.make_future_dataframe(days, include_history=False) + future = forecaster.predict(future) + # this gives ~ 215.77 + res = self.rmse(future['yhat'], test['y']) + self.assertTrue(236 > res > 193, msg="backend: {}".format(forecaster.stan_backend)) def test_fit_predict_no_seasons(self): N = DATA.shape[0] @@ -63,7 +109,8 @@ class TestProphet(TestCase): def test_fit_changepoint_not_in_history(self): train = DATA[(DATA['ds'] < '2013-01-01') | (DATA['ds'] > '2014-01-01')] future = pd.DataFrame({'ds': DATA['ds']}) - forecaster = Prophet(changepoints=['2013-06-06']) + prophet = Prophet(changepoints=['2013-06-06']) + forecaster = prophet forecaster.fit(train) forecaster.predict(future) @@ -121,7 +168,7 @@ class TestProphet(TestCase): self.assertTrue('y_scaled' in history) self.assertEqual(history['y_scaled'].max(), 1.0) - + def test_setup_dataframe_ds_column(self): "Test case where 'ds' exists as an index name and column" df = DATA.copy() @@ -228,7 +275,7 @@ class TestProphet(TestCase): true_values = np.array([ 0.7818315, 0.6234898, 0.9749279, -0.2225209, 0.4338837, -0.9009689 ]) - self.assertAlmostEqual(np.sum((mat[0] - true_values)**2), 0.0) + self.assertAlmostEqual(np.sum((mat[0] - true_values) ** 2), 0.0) def test_fourier_series_yearly(self): mat = Prophet.fourier_series(DATA['ds'], 365.25, 3) @@ -236,7 +283,7 @@ class TestProphet(TestCase): true_values = np.array([ 0.7006152, -0.7135393, -0.9998330, 0.01827656, 0.7262249, 0.6874572 ]) - self.assertAlmostEqual(np.sum((mat[0] - true_values)**2), 0.0) + self.assertAlmostEqual(np.sum((mat[0] - true_values) ** 2), 0.0) def test_growth_init(self): model = Prophet(growth='logistic') @@ -593,7 +640,7 @@ class TestProphet(TestCase): # Test priors m = Prophet( holidays=holidays, yearly_seasonality=False, - seasonality_mode='multiplicative', + seasonality_mode='multiplicative' ) m.add_seasonality(name='monthly', period=30, fourier_order=5, prior_scale=2., mode='additive') @@ -654,7 +701,6 @@ class TestProphet(TestCase): self.assertTrue(np.array_equal((seasonal_features[conditional_weekly_columns] != 0).any(axis=1).values, df['is_conditional_week'].values)) - def test_added_regressors(self): m = Prophet() m.add_regressor('binary_feature', prior_scale=0.2) @@ -740,7 +786,7 @@ class TestProphet(TestCase): self.assertAlmostEqual( fcst['additive_terms'][0], fcst['yearly'][0] + fcst['weekly'][0] - + fcst['extra_regressors_additive'][0], + + fcst['extra_regressors_additive'][0], ) self.assertAlmostEqual( fcst['multiplicative_terms'][0], @@ -749,7 +795,7 @@ class TestProphet(TestCase): self.assertAlmostEqual( fcst['yhat'][0], fcst['trend'][0] * (1 + fcst['multiplicative_terms'][0]) - + fcst['additive_terms'][0], + + fcst['additive_terms'][0], ) # Check works if constant extra regressor at 0 df['constant_feature'] = 0 @@ -800,5 +846,5 @@ class TestProphet(TestCase): {'weekly', 'yearly', 'xmas', 'numeric_feature', 'multiplicative_terms', 'extra_regressors_multiplicative', 'holidays', - }, + }, ) diff --git a/python/requirements.txt b/python/requirements.txt index f474113..baf7893 100644 --- a/python/requirements.txt +++ b/python/requirements.txt @@ -1,4 +1,5 @@ Cython>=0.22 +cmdstanpy==0.4 pystan>=2.14 numpy>=1.10.0 pandas>=0.23.4 diff --git a/python/setup.py b/python/setup.py index 4a70b17..8f3433f 100644 --- a/python/setup.py +++ b/python/setup.py @@ -4,10 +4,9 @@ # LICENSE file in the root directory of this source tree. import os.path -import pickle import platform import sys - +import os from pkg_resources import ( normalize_path, working_set, @@ -18,7 +17,8 @@ from setuptools import setup, find_packages from setuptools.command.build_py import build_py from setuptools.command.develop import develop from setuptools.command.test import test as test_command - +from fbprophet.models import StanBackendEnum +from typing import List PLATFORM = 'unix' if platform.platform().startswith('Win'): @@ -28,15 +28,14 @@ MODEL_DIR = os.path.join('stan', PLATFORM) MODEL_TARGET_DIR = os.path.join('fbprophet', 'stan_model') -def build_stan_model(target_dir, model_dir=MODEL_DIR): - from pystan import StanModel - model_name = 'prophet.stan' - target_name = 'prophet_model.pkl' - with open(os.path.join(model_dir, model_name)) as f: - model_code = f.read() - sm = StanModel(model_code=model_code) - with open(os.path.join(target_dir, target_name), 'wb') as f: - pickle.dump(sm, f, protocol=pickle.HIGHEST_PROTOCOL) +def get_backends_from_env() -> List[str]: + import os + return os.environ.get("STAN_BACKEND", StanBackendEnum.PYSTAN.name).split(",") + + +def build_models(target_dir): + for backend in get_backends_from_env(): + StanBackendEnum.get_backend_class(backend).build_model(target_dir, MODEL_DIR) class BuildPyCommand(build_py): @@ -46,7 +45,7 @@ class BuildPyCommand(build_py): if not self.dry_run: target_dir = os.path.join(self.build_lib, MODEL_TARGET_DIR) self.mkpath(target_dir) - build_stan_model(target_dir) + build_models(target_dir) build_py.run(self) @@ -58,12 +57,29 @@ class DevelopCommand(develop): if not self.dry_run: target_dir = os.path.join(self.setup_path, MODEL_TARGET_DIR) self.mkpath(target_dir) - build_stan_model(target_dir) + build_models(target_dir) develop.run(self) class TestCommand(test_command): + user_options = [ + ('test-module=', 'm', "Run 'test_suite' in specified module"), + ('test-suite=', 's', + "Run single test, case or suite (e.g. 'module.test_suite')"), + ('test-runner=', 'r', "Test runner to use"), + ('test-slow', 'w', "Test slow suites (default off)"), + ] + + def initialize_options(self): + super(TestCommand, self).initialize_options() + self.test_slow = False + + def finalize_options(self): + super(TestCommand, self).finalize_options() + if self.test_slow is None: + self.test_slow = getattr(self.distribution, 'test_slow', False) + """We must run tests on the build directory, not source.""" def with_project_on_sys_path(self, func): @@ -97,6 +113,7 @@ class TestCommand(test_command): sys.modules.update(old_modules) working_set.__init__() + with open('requirements.txt', 'r') as f: install_requires = f.read().splitlines() @@ -125,7 +142,7 @@ setup( 'Programming Language :: Python :: 2.7', 'Programming Language :: Python :: 3', 'Programming Language :: Python :: 3.7', - ], + ], long_description=""" Implements a procedure for forecasting time series data based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality, plus holiday effects. It works best with time series that have strong seasonal effects and several seasons of historical data. Prophet is robust to missing data and shifts in the trend, and typically handles outliers well. """