Replaced pystan with cmdstanpy (#1083)

* changes

* added actual tests for fit method

* precision

* syntax

* sampling not working

* sampling seems to work

* sampling not working again

* sampling works, tests to be removed

* replaced data with rmse

* replace pystan with cmdstanpy

* cleanup

* cleanup

* test for newton

* added support for multiple backends

* minor fixes

* fixed comment

* added support for --test-slow flag

* fixed import

* reverted style change

* specify backend based on env variable

* fixes

* PR fixes
This commit is contained in:
Christopher Suchanek 2020-02-07 14:34:08 -08:00 committed by GitHub
parent 496facb152
commit 1d18adc0ba
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 391 additions and 92 deletions

View file

@ -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

View file

@ -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)

View file

@ -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:

View file

@ -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))

View file

@ -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',
},
},
)

View file

@ -1,4 +1,5 @@
Cython>=0.22
cmdstanpy==0.4
pystan>=2.14
numpy>=1.10.0
pandas>=0.23.4

View file

@ -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.
"""