prophet

  36class Prophet:
  37    """Prophet forecaster.
  38
  39    Parameters
  40    ----------
  41    growth: String 'linear', 'logistic' or 'flat' to specify a linear, logistic or
  42        flat trend.
  43    changepoints: List of dates at which to include potential changepoints. If
  44        not specified, potential changepoints are selected automatically.
  45    n_changepoints: Number of potential changepoints to include. Not used
  46        if input `changepoints` is supplied. If `changepoints` is not supplied,
  47        then n_changepoints potential changepoints are selected uniformly from
  48        the first `changepoint_range` proportion of the history.
  49    changepoint_range: Proportion of history in which trend changepoints will
  50        be estimated. Defaults to 0.8 for the first 80%. Not used if
  51        `changepoints` is specified.
  52    yearly_seasonality: Fit yearly seasonality.
  53        Can be 'auto', True, False, or a number of Fourier terms to generate.
  54    weekly_seasonality: Fit weekly seasonality.
  55        Can be 'auto', True, False, or a number of Fourier terms to generate.
  56    daily_seasonality: Fit daily seasonality.
  57        Can be 'auto', True, False, or a number of Fourier terms to generate.
  58    holidays: pd.DataFrame with columns holiday (string) and ds (date type)
  59        and optionally columns lower_window and upper_window which specify a
  60        range of days around the date to be included as holidays.
  61        lower_window=-2 will include 2 days prior to the date as holidays. Also
  62        optionally can have a column prior_scale specifying the prior scale for
  63        that holiday.
  64    seasonality_mode: 'additive' (default) or 'multiplicative'.
  65    seasonality_prior_scale: Parameter modulating the strength of the
  66        seasonality model. Larger values allow the model to fit larger seasonal
  67        fluctuations, smaller values dampen the seasonality. Can be specified
  68        for individual seasonalities using add_seasonality.
  69    holidays_prior_scale: Parameter modulating the strength of the holiday
  70        components model, unless overridden in the holidays input.
  71    changepoint_prior_scale: Parameter modulating the flexibility of the
  72        automatic changepoint selection. Large values will allow many
  73        changepoints, small values will allow few changepoints.
  74    mcmc_samples: Integer, if greater than 0, will do full Bayesian inference
  75        with the specified number of MCMC samples. If 0, will do MAP
  76        estimation.
  77    interval_width: Float, width of the uncertainty intervals provided
  78        for the forecast. If mcmc_samples=0, this will be only the uncertainty
  79        in the trend using the MAP estimate of the extrapolated generative
  80        model. If mcmc.samples>0, this will be integrated over all model
  81        parameters, which will include uncertainty in seasonality.
  82    uncertainty_samples: Number of simulated draws used to estimate
  83        uncertainty intervals. Settings this value to 0 or False will disable
  84        uncertainty estimation and speed up the calculation.
  85    stan_backend: str as defined in StanBackendEnum default: None - will try to
  86        iterate over all available backends and find the working one.
  87    scaling: 'absmax' (default) or 'minmax'.
  88    holidays_mode: 'additive' or 'multiplicative'. Defaults to seasonality_mode.
  89    """
  90
  91    growth: Literal["linear", "logistic", "flat"]
  92    changepoints: pd.Series[pd.Timestamp] | None
  93    n_changepoints: int
  94    specified_changepoints: bool
  95    changepoint_range: float
  96    yearly_seasonality: Literal["auto"] | int
  97    weekly_seasonality: Literal["auto"] | int
  98    daily_seasonality: Literal["auto"] | int
  99    holidays: pd.DataFrame | None
 100    seasonality_mode: _Mode
 101    holidays_mode: _Mode
 102    seasonality_prior_scale: float
 103    holidays_prior_scale: float
 104    changepoint_prior_scale: float
 105    mcmc_samples: int
 106    interval_width: float
 107    uncertainty_samples: int
 108    scaling: Literal["absmax", "minmax"]
 109
 110    start: pd.Timestamp | None
 111    y_min: float | None
 112    y_scale: float | None
 113    logistic_floor: bool
 114    t_scale: pd.Timedelta | None
 115    changepoints_t: npt.NDArray[np.float64] | None
 116    seasonalities: OrderedDict[str, dict[str, Any]]
 117    extra_regressors: OrderedDict[str, dict[str, Any]]
 118    country_holidays: str | None
 119    stan_fit: Any | None
 120    params: dict[str, Any]
 121    history: pd.DataFrame | None
 122    history_dates: pd.Series[pd.Timestamp] | None
 123    train_component_cols: pd.DataFrame | None
 124    component_modes: dict[_Mode, list[str]] | None
 125    train_holiday_names: pd.Series | None
 126    fit_kwargs: dict[str, Any]
 127    stan_backend: IStanBackend | None
 128    # Optional marker set on nested regressor-predictor models for diagnostics/tests.
 129    _regressor_name: str | None
 130
 131    def __init__(
 132        self,
 133        growth: Literal["linear", "logistic", "flat"] = "linear",
 134        changepoints: pd.Series | list[pd.Timestamp] | None = None,
 135        n_changepoints: int = 25,
 136        changepoint_range: float = 0.8,
 137        yearly_seasonality: Literal["auto"] | int = "auto",
 138        weekly_seasonality: Literal["auto"] | int = "auto",
 139        daily_seasonality: Literal["auto"] | int = "auto",
 140        holidays: pd.DataFrame | None = None,
 141        seasonality_mode: _Mode = "additive",
 142        seasonality_prior_scale: SupportsFloat = 10.0,
 143        holidays_prior_scale: SupportsFloat = 10.0,
 144        changepoint_prior_scale: SupportsFloat = 0.05,
 145        mcmc_samples: int = 0,
 146        interval_width: float = 0.80,
 147        uncertainty_samples: int = 1000,
 148        stan_backend: str | None = None,
 149        scaling: Literal["absmax", "minmax"] = "absmax",
 150        holidays_mode: _Mode | None = None,
 151    ) -> None:
 152        self.growth = growth
 153
 154        if changepoints is not None:
 155            self.changepoints = pd.Series(pd.to_datetime(changepoints), name='ds')
 156            self.n_changepoints = len(self.changepoints)
 157            self.specified_changepoints = True
 158        else:
 159            self.changepoints = changepoints
 160            self.n_changepoints = n_changepoints
 161            self.specified_changepoints = False
 162
 163        self.changepoint_range = changepoint_range
 164        self.yearly_seasonality = yearly_seasonality
 165        self.weekly_seasonality = weekly_seasonality
 166        self.daily_seasonality = daily_seasonality
 167        self.holidays = holidays
 168
 169        self.seasonality_mode = seasonality_mode
 170        self.holidays_mode = holidays_mode or self.seasonality_mode
 171
 172        self.seasonality_prior_scale = float(seasonality_prior_scale)
 173        self.changepoint_prior_scale = float(changepoint_prior_scale)
 174        self.holidays_prior_scale = float(holidays_prior_scale)
 175
 176        self.mcmc_samples = mcmc_samples
 177        self.interval_width = interval_width
 178        self.uncertainty_samples = uncertainty_samples
 179        if scaling not in ("absmax", "minmax"):
 180            raise ValueError("scaling must be one of 'absmax' or 'minmax'")
 181        self.scaling = scaling
 182
 183        # Set during fitting or by other methods
 184        self.start = None
 185        self.y_min = None
 186        self.y_scale = None
 187        self.logistic_floor = False
 188        self.t_scale = None
 189        self.changepoints_t = None
 190        self.seasonalities = OrderedDict({})
 191        self.extra_regressors = OrderedDict({})
 192        self.country_holidays = None
 193        self.stan_fit = None
 194        self.params = {}
 195        self.history = None
 196        self.history_dates = None
 197        self.train_component_cols = None
 198        self.component_modes = None
 199        self.train_holiday_names = None
 200        self.fit_kwargs = {}
 201        self._regressor_name = None
 202        self.validate_inputs()
 203        self._load_stan_backend(stan_backend)
 204
 205    def _load_stan_backend(self, stan_backend: str | None) -> None:
 206        if stan_backend is None:
 207            for i in StanBackendEnum:
 208                try:
 209                    logger.debug("Trying to load backend: %s", i.name)
 210                    return self._load_stan_backend(i.name)
 211                except Exception as e:
 212                    logger.debug("Unable to load backend %s (%s), trying the next one", i.name, e)
 213        else:
 214            # pyrefly:ignore[not-callable]
 215            self.stan_backend = StanBackendEnum.get_backend_class(stan_backend)()
 216
 217        assert self.stan_backend
 218        logger.debug("Loaded stan backend: %s", self.stan_backend.get_type())
 219
 220    def validate_inputs(self) -> None:
 221        """Validates the inputs to Prophet."""
 222        if self.growth not in ('linear', 'logistic', 'flat'):
 223            raise ValueError(
 224                'Parameter "growth" should be "linear", "logistic" or "flat".')
 225        if not isinstance(self.changepoint_range, (int, float)):
 226            raise ValueError("changepoint_range must be a number in [0, 1]'")
 227        if ((self.changepoint_range < 0) or (self.changepoint_range > 1)):
 228            raise ValueError('Parameter "changepoint_range" must be in [0, 1]')
 229        if self.holidays is not None:
 230            if not (
 231                isinstance(self.holidays, pd.DataFrame)
 232                and 'ds' in self.holidays  # noqa W503
 233                and 'holiday' in self.holidays  # noqa W503
 234            ):
 235                raise ValueError('holidays must be a DataFrame with "ds" and '
 236                                 '"holiday" columns.')
 237            self.holidays['ds'] = pd.to_datetime(self.holidays['ds'])
 238            if (
 239                self.holidays['ds'].isnull().any()
 240                or self.holidays['holiday'].isnull().any()
 241            ):
 242                raise ValueError('Found a NaN in holidays dataframe.')
 243            has_lower = 'lower_window' in self.holidays
 244            has_upper = 'upper_window' in self.holidays
 245            if has_lower + has_upper == 1:
 246                raise ValueError('Holidays must have both lower_window and ' +
 247                                 'upper_window, or neither')
 248            if has_lower:
 249                if self.holidays['lower_window'].max() > 0:
 250                    raise ValueError('Holiday lower_window should be <= 0')
 251                if self.holidays['upper_window'].min() < 0:
 252                    raise ValueError('Holiday upper_window should be >= 0')
 253            for h in self.holidays['holiday'].unique():
 254                self.validate_column_name(h, check_holidays=False)
 255        if self.seasonality_mode not in ['additive', 'multiplicative']:
 256            raise ValueError(
 257                'seasonality_mode must be "additive" or "multiplicative"'
 258            )
 259        if self.holidays_mode not in ['additive', 'multiplicative']:
 260            raise ValueError(
 261                'holidays_mode must be "additive" or "multiplicative"'
 262            )
 263
 264    def validate_column_name(
 265        self,
 266        name: str,
 267        check_holidays: bool = True,
 268        check_seasonalities: bool = True,
 269        check_regressors: bool = True,
 270    ) -> None:
 271        """Validates the name of a seasonality, holiday, or regressor.
 272
 273        Parameters
 274        ----------
 275        name: string
 276        check_holidays: bool check if name already used for holiday
 277        check_seasonalities: bool check if name already used for seasonality
 278        check_regressors: bool check if name already used for regressor
 279        """
 280        if '_delim_' in name:
 281            raise ValueError('Name cannot contain "_delim_"')
 282        reserved_names = [
 283            'trend', 'additive_terms', 'daily', 'weekly', 'yearly',
 284            'holidays', 'zeros', 'extra_regressors_additive', 'yhat',
 285            'extra_regressors_multiplicative', 'multiplicative_terms',
 286        ]
 287        rn_l = [n + '_lower' for n in reserved_names]
 288        rn_u = [n + '_upper' for n in reserved_names]
 289        reserved_names.extend(rn_l)
 290        reserved_names.extend(rn_u)
 291        reserved_names.extend([
 292            'ds', 'y', 'cap', 'floor', 'y_scaled', 'cap_scaled'])
 293        if name in reserved_names:
 294            raise ValueError(
 295                'Name {name!r} is reserved.'.format(name=name)
 296            )
 297        if (check_holidays and self.holidays is not None and
 298                name in self.holidays['holiday'].unique()):
 299            raise ValueError(
 300                'Name {name!r} already used for a holiday.'.format(name=name)
 301            )
 302        if (check_holidays and self.country_holidays is not None and
 303                name in get_holiday_names(self.country_holidays)):
 304            raise ValueError(
 305                'Name {name!r} is a holiday name in {country_holidays}.'
 306                .format(name=name, country_holidays=self.country_holidays)
 307            )
 308        if check_seasonalities and name in self.seasonalities:
 309            raise ValueError(
 310                'Name {name!r} already used for a seasonality.'
 311                .format(name=name)
 312            )
 313        if check_regressors and name in self.extra_regressors:
 314            raise ValueError(
 315                'Name {name!r} already used for an added regressor.'
 316                .format(name=name)
 317            )
 318
 319    def setup_dataframe(self, df: pd.DataFrame, initialize_scales: bool = False) -> pd.DataFrame:
 320        """Prepare dataframe for fitting or predicting.
 321
 322        Adds a time index and scales y. Creates auxiliary columns 't', 't_ix',
 323        'y_scaled', and 'cap_scaled'. These columns are used during both
 324        fitting and predicting.
 325
 326        Parameters
 327        ----------
 328        df: pd.DataFrame with columns ds, y, and cap if logistic growth. Any
 329            specified additional regressors must also be present.
 330        initialize_scales: Boolean set scaling factors in self from df.
 331
 332        Returns
 333        -------
 334        pd.DataFrame prepared for fitting or predicting.
 335        """
 336        if 'y' in df:  # 'y' will be in training data
 337            df['y'] = pd.to_numeric(df['y'])
 338            if np.isinf(df['y'].values).any():
 339                raise ValueError('Found infinity in column y.')
 340        if df['ds'].dtype == np.int64:
 341            df['ds'] = df['ds'].astype(str)
 342        df['ds'] = pd.to_datetime(df['ds'])
 343        if df['ds'].dt.tz is not None:
 344            raise ValueError(
 345                'Column ds has timezone specified, which is not supported. '
 346                'Remove timezone.'
 347            )
 348        if df['ds'].isnull().any():
 349            raise ValueError('Found NaN in column ds.')
 350        for name in self.extra_regressors:
 351            if name not in df:
 352                raise ValueError(
 353                    'Regressor {name!r} missing from dataframe'
 354                    .format(name=name)
 355                )
 356            df[name] = pd.to_numeric(df[name])
 357            if df[name].isnull().any():
 358                raise ValueError(
 359                    'Found NaN in column {name!r}'.format(name=name)
 360                )
 361        for props in self.seasonalities.values():
 362            condition_name = props['condition_name']
 363            if condition_name is not None:
 364                if condition_name not in df:
 365                    raise ValueError(
 366                        'Condition {condition_name!r} missing from dataframe'
 367                        .format(condition_name=condition_name)
 368                    )
 369                if not df[condition_name].isin([True, False]).all():
 370                    raise ValueError(
 371                        'Found non-boolean in column {condition_name!r}'
 372                        .format(condition_name=condition_name)
 373                    )
 374                df[condition_name] = df[condition_name].astype('bool')
 375
 376        if df.index.name == 'ds':
 377            df.index.name = None
 378        df = df.sort_values('ds', kind='mergesort')
 379        df = df.reset_index(drop=True)
 380
 381        self.initialize_scales(initialize_scales, df)
 382
 383        if self.logistic_floor:
 384            if 'floor' not in df:
 385                raise ValueError('Expected column "floor".')
 386        else:
 387            if self.scaling == "absmax":
 388                df['floor'] = 0.
 389            elif self.scaling == "minmax":
 390                df['floor'] = self.y_min
 391        if self.growth == 'logistic':
 392            if 'cap' not in df:
 393                raise ValueError(
 394                    'Capacities must be supplied for logistic growth in '
 395                    'column "cap"'
 396                )
 397            if (df['cap'] <= df['floor']).any():
 398                raise ValueError(
 399                    'cap must be greater than floor (which defaults to 0).'
 400                )
 401            df['cap_scaled'] = (df['cap'] - df['floor']) / self.y_scale  # pyrefly:ignore[unsupported-operation]
 402
 403        df['t'] = (df['ds'] - self.start) / self.t_scale  # pyrefly:ignore[unsupported-operation]
 404        if 'y' in df:
 405            df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale  # pyrefly:ignore[unsupported-operation]
 406
 407        for name, props in self.extra_regressors.items():
 408            df[name] = ((df[name] - props['mu']) / props['std'])
 409        return df
 410
 411    def initialize_scales(self, initialize_scales: bool, df: pd.DataFrame) -> None:
 412        """Initialize model scales.
 413
 414        Sets model scaling factors using df.
 415
 416        Parameters
 417        ----------
 418        initialize_scales: Boolean set the scales or not.
 419        df: pd.DataFrame for setting scales.
 420        """
 421        if not initialize_scales:
 422            return
 423
 424        if self.growth == 'logistic' and 'floor' in df:
 425            self.logistic_floor = True
 426            if self.scaling == "absmax":
 427                self.y_min = float((df['y'] - df['floor']).abs().min())
 428                self.y_scale = float((df['y'] - df['floor']).abs().max())
 429            elif self.scaling == "minmax":
 430                self.y_min = df['floor'].min()
 431                self.y_scale = float(df['cap'].max() - self.y_min)
 432        else:
 433            if self.scaling == "absmax":
 434                self.y_min = 0.
 435                self.y_scale = float((df['y']).abs().max())
 436            elif self.scaling == "minmax":
 437                self.y_min = df['y'].min()
 438                self.y_scale =  float(df['y'].max() - self.y_min)
 439        if self.y_scale == 0:
 440            self.y_scale = 1.0
 441
 442        self.start = df['ds'].min()
 443        self.t_scale = df['ds'].max() - self.start
 444        for name, props in self.extra_regressors.items():
 445            standardize = props['standardize']
 446            n_vals = len(df[name].unique())
 447            if n_vals < 2:
 448                standardize = False
 449            if standardize == 'auto':
 450                if set(df[name].unique()) == {1, 0}:
 451                    standardize = False #  Don't standardize binary variables.
 452                else:
 453                    standardize = True
 454            if standardize:
 455                mu = float(df[name].mean())
 456                std = float(df[name].std())
 457                self.extra_regressors[name]['mu'] = mu
 458                self.extra_regressors[name]['std'] = std
 459
 460    def set_changepoints(self) -> None:
 461        """Set changepoints
 462
 463        Sets m$changepoints to the dates of changepoints. Either:
 464        1) The changepoints were passed in explicitly.
 465            A) They are empty.
 466            B) They are not empty, and need validation.
 467        2) We are generating a grid of them.
 468        3) The user prefers no changepoints be used.
 469        """
 470        if self.changepoints is not None:
 471            if len(self.changepoints) == 0:
 472                pass
 473            else:
 474                history = cast(pd.DataFrame, self.history)
 475                too_low = min(self.changepoints) < history['ds'].min()
 476                too_high = max(self.changepoints) > history['ds'].max()
 477                if too_low or too_high:
 478                    raise ValueError('Changepoints must fall within training data.')
 479        else:
 480            # Place potential changepoints evenly through first
 481            # `changepoint_range` proportion of the history
 482            history = cast(pd.DataFrame, self.history)
 483            hist_size = int(np.floor(history.shape[0] * self.changepoint_range))
 484            if self.n_changepoints + 1 > hist_size:
 485                self.n_changepoints = hist_size - 1
 486                logger.info(
 487                    'n_changepoints greater than number of observations. '
 488                    'Using {n_changepoints}.'
 489                    .format(n_changepoints=self.n_changepoints)
 490                )
 491            if self.n_changepoints > 0:
 492                cp_indexes = (
 493                    np.linspace(0, hist_size - 1, self.n_changepoints + 1)
 494                        .round()
 495                        .astype(int)
 496                )
 497                self.changepoints = history.iloc[cp_indexes]['ds'].tail(-1)
 498            else:
 499                # set empty changepoints
 500                self.changepoints = pd.Series(pd.to_datetime([]), name='ds')
 501        if len(self.changepoints) > 0:
 502            self.changepoints_t = np.sort(np.array(
 503                (self.changepoints - self.start) / self.t_scale))  # pyrefly:ignore[unsupported-operation]
 504        else:
 505            self.changepoints_t = np.array([0])  # dummy changepoint
 506
 507    @staticmethod
 508    def fourier_series(
 509        dates: pd.Series,
 510        period: float,
 511        series_order: int,
 512    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
 513        """Provides Fourier series components with the specified frequency
 514        and order.
 515
 516        Parameters
 517        ----------
 518        dates: pd.Series containing timestamps.
 519        period: Number of days of the period.
 520        series_order: Number of components.
 521
 522        Returns
 523        -------
 524        Matrix with seasonality features.
 525        """
 526        if not (series_order >= 1):
 527            raise ValueError("series_order must be >= 1")
 528
 529        epoch = pd.Timestamp("1970-01-01", tz=dates.dt.tz)
 530        t = (dates - epoch).dt.total_seconds() / (24 * 60 * 60)
 531
 532        x_T = np.pi * 2 * t
 533        fourier_components = np.empty((dates.shape[0], 2 * series_order))
 534        for i in range(series_order):
 535            c = (i + 1) / period * x_T
 536            fourier_components[:, 2 * i] = np.sin(c)
 537            fourier_components[:, (2 * i) + 1] = np.cos(c)
 538        return fourier_components
 539
 540    @classmethod
 541    def make_seasonality_features(
 542        cls,
 543        dates: pd.Series[pd.Timestamp],
 544        period: float,
 545        series_order: int,
 546        prefix: str,
 547    ) -> pd.DataFrame:
 548        """Data frame with seasonality features.
 549
 550        Parameters
 551        ----------
 552        cls: Prophet class.
 553        dates: pd.Series containing timestamps.
 554        period: Number of days of the period.
 555        series_order: Number of components.
 556        prefix: Column name prefix.
 557
 558        Returns
 559        -------
 560        pd.DataFrame with seasonality features.
 561        """
 562        features = cls.fourier_series(dates, period, series_order)
 563        columns = [
 564            '{}_delim_{}'.format(prefix, i + 1)
 565            for i in range(features.shape[1])
 566        ]
 567        return pd.DataFrame(features, columns=columns)
 568
 569    def construct_holiday_dataframe(self, dates: pd.Series[pd.Timestamp]) -> pd.DataFrame:
 570        """Construct a dataframe of holiday dates.
 571
 572        Will combine self.holidays with the built-in country holidays
 573        corresponding to input dates, if self.country_holidays is set.
 574
 575        Parameters
 576        ----------
 577        dates: pd.Series containing timestamps used for computing seasonality.
 578
 579        Returns
 580        -------
 581        dataframe of holiday dates, in holiday dataframe format used in
 582        initialization.
 583        """
 584        all_holidays = pd.DataFrame()
 585        if self.holidays is not None:
 586            all_holidays = self.holidays.copy()
 587        if self.country_holidays is not None:
 588            year_list = list({x.year for x in dates})
 589            country_holidays_df = make_holidays_df(
 590                year_list=year_list, country=self.country_holidays
 591            )
 592            all_holidays = pd.concat((all_holidays, country_holidays_df),
 593                                     sort=False)
 594            all_holidays.reset_index(drop=True, inplace=True)
 595        # Drop future holidays not previously seen in training data
 596        if self.train_holiday_names is not None:
 597            # Remove holiday names didn't show up in fit
 598            index_to_drop = all_holidays.index[
 599                np.logical_not(
 600                    all_holidays.holiday.isin(self.train_holiday_names)
 601                )
 602            ]
 603            all_holidays = all_holidays.drop(index_to_drop)
 604            # Add holiday names in fit but not in predict with ds as NA
 605            holidays_to_add = pd.DataFrame({
 606                'holiday': self.train_holiday_names[
 607                    np.logical_not(self.train_holiday_names
 608                                       .isin(all_holidays.holiday))
 609                ]
 610            })
 611            all_holidays = pd.concat((all_holidays, holidays_to_add),
 612                                     sort=False)
 613            all_holidays.reset_index(drop=True, inplace=True)
 614        return all_holidays
 615
 616    def make_holiday_features(
 617        self,
 618        dates: pd.Series[pd.Timestamp],
 619        holidays: pd.DataFrame,
 620    ) -> tuple[pd.DataFrame, list[float], list[str]]:
 621        """Construct a dataframe of holiday features.
 622
 623        Parameters
 624        ----------
 625        dates: pd.Series containing timestamps used for computing seasonality.
 626        holidays: pd.Dataframe containing holidays, as returned by
 627            construct_holiday_dataframe.
 628
 629        Returns
 630        -------
 631        holiday_features: pd.DataFrame with a column for each holiday.
 632        prior_scale_list: List of prior scales for each holiday column.
 633        holiday_names: List of names of holidays
 634        """
 635        # Holds columns of our future matrix.
 636        expanded_holidays = defaultdict(lambda: np.zeros(dates.shape[0]))
 637        prior_scales = {}
 638        # Makes an index so we can perform `get_loc` below.
 639        # Strip to just dates.
 640        row_index = pd.DatetimeIndex(dates.dt.date)
 641
 642        for row in holidays.itertuples():
 643            dt = cast(pd.Timestamp, row.ds).date()
 644            try:
 645                lw = int(getattr(row, 'lower_window', 0))
 646                uw = int(getattr(row, 'upper_window', 0))
 647            except ValueError:
 648                lw = 0
 649                uw = 0
 650            ps = float(getattr(row, 'prior_scale', self.holidays_prior_scale))
 651            if np.isnan(ps):
 652                ps = float(self.holidays_prior_scale)
 653            if row.holiday in prior_scales and prior_scales[row.holiday] != ps:
 654                raise ValueError(
 655                    'Holiday {holiday!r} does not have consistent prior '
 656                    'scale specification.'.format(holiday=row.holiday)
 657                )
 658            if ps <= 0:
 659                raise ValueError('Prior scale must be > 0')
 660            prior_scales[row.holiday] = ps
 661
 662            for offset in range(lw, uw + 1):
 663                occurrence = pd.to_datetime(dt + timedelta(days=offset))
 664                try:
 665                    loc = row_index.get_loc(occurrence)
 666                except KeyError:
 667                    loc = None
 668                key = '{}_delim_{}{}'.format(
 669                    row.holiday,
 670                    '+' if offset >= 0 else '-',
 671                    abs(offset)
 672                )
 673                if loc is not None:
 674                    expanded_holidays[key][loc] = 1.
 675                else:
 676                    expanded_holidays[key]  # Access key to generate value
 677        holiday_features = pd.DataFrame(expanded_holidays)
 678        # Make sure column order is consistent
 679        holiday_features = holiday_features[sorted(holiday_features.columns
 680                                                                   .tolist())]
 681        prior_scale_list = [
 682            prior_scales[h.split('_delim_')[0]]
 683            for h in holiday_features.columns
 684        ]
 685        holiday_names = list(prior_scales.keys())
 686        # Store holiday names used in fit
 687        if self.train_holiday_names is None:
 688            self.train_holiday_names = pd.Series(holiday_names)
 689        return holiday_features, prior_scale_list, holiday_names
 690
 691    def add_regressor(
 692        self,
 693        name: str,
 694        prior_scale: float | None = None,
 695        standardize: Literal['auto'] | bool = 'auto',
 696        mode: _Mode | None = None,
 697        regressor_predictor: bool | dict | None = None,
 698    ) -> Self:
 699        """Add an additional regressor to be used for fitting and predicting.
 700
 701        The dataframe passed to `fit` and `predict` will have a column with the
 702        specified name to be used as a regressor. When standardize='auto', the
 703        regressor will be standardized unless it is binary. The regression
 704        coefficient is given a prior with the specified scale parameter.
 705        Decreasing the prior scale will add additional regularization. If no
 706        prior scale is provided, self.holidays_prior_scale will be used.
 707        Mode can be specified as either 'additive' or 'multiplicative'. If not
 708        specified, self.seasonality_mode will be used. 'additive' means the
 709        effect of the regressor will be added to the trend, 'multiplicative'
 710        means it will multiply the trend.
 711
 712        Parameters
 713        ----------
 714        name: string name of the regressor.
 715        prior_scale: optional float scale for the normal prior. If not
 716            provided, self.holidays_prior_scale will be used.
 717        standardize: optional, specify whether this regressor will be
 718            standardized prior to fitting. Can be 'auto' (standardize if not
 719            binary), True, or False.
 720        mode: optional, 'additive' or 'multiplicative'. Defaults to
 721            self.seasonality_mode.
 722        regressor_predictor: optional. If provided, fits a dedicated Prophet
 723            model to forecast this regressor. Set to True to use default Prophet
 724            parameters, or provide a dict of keyword arguments for the
 725            underlying Prophet constructor. When set, future regressor values
 726            will be generated from this model (with uncertainty if available)
 727            instead of using user-supplied point estimates.
 728
 729        Returns
 730        -------
 731        The prophet object.
 732        """
 733        if self.history is not None:
 734            raise Exception(
 735                "Regressors must be added prior to model fitting.")
 736        self.validate_column_name(name, check_regressors=False)
 737        if prior_scale is None:
 738            prior_scale = float(self.holidays_prior_scale)
 739        if mode is None:
 740            mode = self.seasonality_mode
 741        if prior_scale <= 0:
 742            raise ValueError('Prior scale must be > 0')
 743        if mode not in ['additive', 'multiplicative']:
 744            raise ValueError("mode must be 'additive' or 'multiplicative'")
 745        predictor_spec = None
 746        if regressor_predictor:
 747            if isinstance(regressor_predictor, dict):
 748                predictor_spec = deepcopy(regressor_predictor)
 749            else:
 750                # Truthy non-dict values (typically True) use default Prophet params.
 751                predictor_spec = {}
 752        self.extra_regressors[name] = {
 753            'prior_scale': prior_scale,
 754            'standardize': standardize,
 755            'mu': 0.,
 756            'std': 1.,
 757            'mode': mode,
 758            'predictor_spec': predictor_spec,
 759            'predictor': None,
 760        }
 761        return self
 762
 763    def _fit_regressor_models(self, df: pd.DataFrame):
 764        """Fit regressor-specific Prophet models when configured."""
 765        for name, props in self.extra_regressors.items():
 766            spec = props.get('predictor_spec')
 767            if spec is None:
 768                continue
 769            regressor_df = df[['ds', name]].copy()
 770            regressor_df = regressor_df[regressor_df[name].notnull()]
 771            if regressor_df.shape[0] < 2:
 772                raise ValueError(f"Not enough data to fit regressor model for '{name}'.")
 773            regressor_df = regressor_df.rename(columns={name: 'y'})
 774            logger.info("Fitting regressor model '%s' with %d observations", name, regressor_df.shape[0])
 775            predictor_kwargs = deepcopy(spec)
 776            if 'stan_backend' not in predictor_kwargs and self.stan_backend is not None:
 777                predictor_kwargs['stan_backend'] = self.stan_backend.get_type()
 778            predictor = self.__class__(**predictor_kwargs)
 779            predictor._regressor_name = name  # marker for diagnostics/testing
 780            predictor.fit(regressor_df, **self.fit_kwargs)
 781            props['predictor'] = predictor
 782
 783    def _prepare_regressors_for_predict(
 784        self, df: pd.DataFrame, n_samples: int | None = None
 785    ) -> tuple[pd.DataFrame, dict[str, npt.NDArray[np.float64]] | None]:
 786        """Populate regressor columns using fitted regressor models if present.
 787
 788        Returns updated dataframe and optional regressor sample draws.
 789        """
 790        regressor_samples: dict[str, npt.NDArray[np.float64]] | None = (
 791            {} if n_samples else None
 792        )
 793        last_history_date = self.history['ds'].max() if self.history is not None else None
 794        for name, props in self.extra_regressors.items():
 795            predictor = props.get('predictor')
 796            if props.get('predictor_spec') is None or predictor is None:
 797                continue
 798            if name not in df:
 799                df[name] = np.nan
 800            future_mask = (
 801                df['ds'] > last_history_date
 802                if last_history_date is not None
 803                else df[name].isna()
 804            )
 805            history_mask = ~future_mask
 806            if history_mask.any():
 807                missing_hist = history_mask & df[name].isna()
 808                if missing_hist.any() and self.history is not None and name in self.history:
 809                    hist_lookup = self.history.set_index('ds')[name]
 810                    props = self.extra_regressors[name]
 811                    hist_lookup = hist_lookup * props['std'] + props['mu']
 812                    df.loc[missing_hist, name] = hist_lookup.reindex(
 813                        df.loc[missing_hist, 'ds']
 814                    ).to_numpy()
 815            if future_mask.any():
 816                future_dates = df.loc[future_mask, 'ds']
 817                reg_future = pd.DataFrame({'ds': future_dates})
 818                logger.info(
 819                    "Running regressor model '%s' for %d future dates",
 820                    name,
 821                    len(reg_future),
 822                )
 823                reg_pred = predictor.predict(reg_future)
 824                df.loc[future_mask, name] = reg_pred['yhat'].to_numpy()
 825                if n_samples and regressor_samples is not None:
 826                    sample_dict = predictor.predictive_samples(reg_future)
 827                    reg_samples = sample_dict['yhat']
 828                    samples_full = np.tile(
 829                        df[name].to_numpy()[:, None], (1, reg_samples.shape[1])
 830                    )
 831                    # Use boolean indexing on ndarray, not pandas ExtensionArray 2D slice.
 832                    samples_full[future_mask.to_numpy(), :] = reg_samples
 833                    regressor_samples[name] = samples_full
 834                    logger.info(
 835                        "Collected %d predictive samples for regressor '%s'",
 836                        reg_samples.shape[1],
 837                        name,
 838                    )
 839            elif n_samples and regressor_samples is not None:
 840                n_draws = n_samples if isinstance(n_samples, int) else 1
 841                regressor_samples[name] = np.tile(
 842                    df[name].to_numpy()[:, None], (1, n_draws)
 843                )
 844        if regressor_samples is not None and len(regressor_samples) == 0:
 845            regressor_samples = None
 846        return df, regressor_samples
 847
 848    def add_seasonality(
 849        self,
 850        name: str,
 851        period: float,
 852        fourier_order: int,
 853        prior_scale: float | None = None,
 854        mode: _Mode | None = None,
 855        condition_name: str | None = None,
 856    ) -> Self:
 857        """Add a seasonal component with specified period, number of Fourier
 858        components, and prior scale.
 859
 860        Increasing the number of Fourier components allows the seasonality to
 861        change more quickly (at risk of overfitting). Default values for yearly
 862        and weekly seasonalities are 10 and 3 respectively.
 863
 864        Increasing prior scale will allow this seasonality component more
 865        flexibility, decreasing will dampen it. If not provided, will use the
 866        seasonality_prior_scale provided on Prophet initialization (defaults
 867        to 10).
 868
 869        Mode can be specified as either 'additive' or 'multiplicative'. If not
 870        specified, self.seasonality_mode will be used (defaults to additive).
 871        Additive means the seasonality will be added to the trend,
 872        multiplicative means it will multiply the trend.
 873
 874        If condition_name is provided, the dataframe passed to `fit` and
 875        `predict` should have a column with the specified condition_name
 876        containing booleans which decides when to apply seasonality.
 877
 878        Parameters
 879        ----------
 880        name: string name of the seasonality component.
 881        period: float number of days in one period.
 882        fourier_order: int number of Fourier components to use.
 883        prior_scale: optional float prior scale for this component.
 884        mode: optional 'additive' or 'multiplicative'
 885        condition_name: string name of the seasonality condition.
 886
 887        Returns
 888        -------
 889        The prophet object.
 890        """
 891        if self.history is not None:
 892            raise Exception(
 893                'Seasonality must be added prior to model fitting.')
 894        if name not in ['daily', 'weekly', 'yearly']:
 895            # Allow overwriting built-in seasonalities
 896            self.validate_column_name(name, check_seasonalities=False)
 897        if prior_scale is None:
 898            ps = self.seasonality_prior_scale
 899        else:
 900            ps = float(prior_scale)
 901        if ps <= 0:
 902            raise ValueError('Prior scale must be > 0')
 903        if fourier_order <= 0:
 904            raise ValueError('Fourier Order must be > 0')
 905        if mode is None:
 906            mode = self.seasonality_mode
 907        if mode not in ['additive', 'multiplicative']:
 908            raise ValueError('mode must be "additive" or "multiplicative"')
 909        if condition_name is not None:
 910            self.validate_column_name(condition_name)
 911        self.seasonalities[name] = {
 912            'period': period,
 913            'fourier_order': fourier_order,
 914            'prior_scale': ps,
 915            'mode': mode,
 916            'condition_name': condition_name,
 917        }
 918        return self
 919
 920    def add_country_holidays(self, country_name: str) -> Self:
 921        """Add in built-in holidays for the specified country.
 922
 923        These holidays will be included in addition to any specified on model
 924        initialization.
 925
 926        Holidays will be calculated for arbitrary date ranges in the history
 927        and future. See the online documentation for the list of countries with
 928        built-in holidays.
 929
 930        Built-in country holidays can only be set for a single country.
 931
 932        Parameters
 933        ----------
 934        country_name: Name of the country, like 'UnitedStates' or 'US'
 935
 936        Returns
 937        -------
 938        The prophet object.
 939        """
 940        if self.history is not None:
 941            raise Exception(
 942                "Country holidays must be added prior to model fitting."
 943            )
 944        # Validate names.
 945        for name in get_holiday_names(country_name):
 946            # Allow merging with existing holidays
 947            self.validate_column_name(name, check_holidays=False)
 948        # Set the holidays.
 949        if self.country_holidays is not None:
 950            logger.warning(
 951                'Changing country holidays from {country_holidays!r} to '
 952                '{country_name!r}.'
 953                .format(
 954                    country_holidays=self.country_holidays,
 955                    country_name=country_name,
 956                )
 957            )
 958        self.country_holidays = country_name
 959        return self
 960
 961    def make_all_seasonality_features(self, df: pd.DataFrame) -> tuple[
 962        pd.DataFrame,
 963        list[float],
 964        pd.DataFrame,
 965        dict[_Mode, list[str]],
 966    ]:
 967        """Dataframe with seasonality features.
 968
 969        Includes seasonality features, holiday features, and added regressors.
 970
 971        Parameters
 972        ----------
 973        df: pd.DataFrame with dates for computing seasonality features and any
 974            added regressors.
 975
 976        Returns
 977        -------
 978        pd.DataFrame with regression features.
 979        list of prior scales for each column of the features dataframe.
 980        Dataframe with indicators for which regression components correspond to
 981            which columns.
 982        Dictionary with keys 'additive' and 'multiplicative' listing the
 983            component names for each mode of seasonality.
 984        """
 985        seasonal_features = []
 986        prior_scales = []
 987        modes: dict[_Mode, list[str]] = {'additive': [], 'multiplicative': []}
 988
 989        # Seasonality features
 990        for name, props in self.seasonalities.items():
 991            features = self.make_seasonality_features(
 992                df['ds'],
 993                props['period'],
 994                props['fourier_order'],
 995                name,
 996            )
 997            if props['condition_name'] is not None:
 998                features[~df[props['condition_name']]] = 0
 999            seasonal_features.append(features)
1000            prior_scales.extend(
1001                [props['prior_scale']] * features.shape[1])
1002            modes[props['mode']].append(name)
1003
1004        # Holiday features
1005        holidays = self.construct_holiday_dataframe(df['ds'])
1006        if len(holidays) > 0:
1007            features, holiday_priors, holiday_names = (
1008                self.make_holiday_features(df['ds'], holidays)
1009            )
1010            seasonal_features.append(features)
1011            prior_scales.extend(holiday_priors)
1012            modes[self.holidays_mode].extend(holiday_names)
1013
1014        # Additional regressors
1015        for name, props in self.extra_regressors.items():
1016            seasonal_features.append(pd.DataFrame(df[name]))
1017            prior_scales.append(props['prior_scale'])
1018            modes[props['mode']].append(name)
1019
1020        # Dummy to prevent empty X
1021        if len(seasonal_features) == 0:
1022            seasonal_features.append(
1023                pd.DataFrame({'zeros': np.zeros(df.shape[0])}))
1024            prior_scales.append(1.)
1025
1026        seasonal_features = pd.concat(seasonal_features, axis=1)
1027        component_cols, modes = self.regressor_column_matrix(
1028            seasonal_features, modes
1029        )
1030        return seasonal_features, prior_scales, component_cols, modes
1031
1032    def regressor_column_matrix(
1033        self,
1034        seasonal_features: pd.DataFrame,
1035        modes: dict[_Mode, list[str]],
1036    ) -> tuple[pd.DataFrame, dict[_Mode, list[str]]]:
1037        """Dataframe indicating which columns of the feature matrix correspond
1038        to which seasonality/regressor components.
1039
1040        Includes combination components, like 'additive_terms'. These
1041        combination components will be added to the 'modes' input.
1042
1043        Parameters
1044        ----------
1045        seasonal_features: Constructed seasonal features dataframe
1046        modes: Dictionary with keys 'additive' and 'multiplicative' listing the
1047            component names for each mode of seasonality.
1048
1049        Returns
1050        -------
1051        component_cols: A binary indicator dataframe with columns seasonal
1052            components and rows columns in seasonal_features. Entry is 1 if
1053            that columns is used in that component.
1054        modes: Updated input with combination components.
1055        """
1056        components = pd.DataFrame({
1057            'col': np.arange(seasonal_features.shape[1]),
1058            'component': [
1059                x.split('_delim_')[0] for x in seasonal_features.columns
1060            ],
1061        })
1062        # Add total for holidays
1063        if self.train_holiday_names is not None:
1064            components = self.add_group_component(
1065                components, 'holidays', self.train_holiday_names.unique())
1066        # Add totals additive and multiplicative components, and regressors
1067        for mode in ('additive', 'multiplicative'):
1068            components = self.add_group_component(
1069                components, mode + '_terms', modes[mode]
1070            )
1071            regressors_by_mode = [
1072                r for r, props in self.extra_regressors.items()
1073                if props['mode'] == mode
1074            ]
1075            components = self.add_group_component(
1076                components, 'extra_regressors_' + mode, regressors_by_mode)
1077            # Add combination components to modes
1078            modes[mode].append(mode + '_terms')
1079            modes[mode].append('extra_regressors_' + mode)
1080        # After all of the additive/multiplicative groups have been added,
1081        modes[self.holidays_mode].append('holidays')
1082        # Convert to a binary matrix
1083        component_cols = pd.crosstab(
1084            components['col'], components['component'],
1085        ).sort_index(level='col')
1086        # Add columns for additive and multiplicative terms, if missing
1087        for name in ['additive_terms', 'multiplicative_terms']:
1088            if name not in component_cols:
1089                component_cols[name] = 0
1090        # Remove the placeholder
1091        component_cols.drop('zeros', axis=1, inplace=True, errors='ignore')
1092        # Validation
1093        if (max(component_cols['additive_terms']
1094            + component_cols['multiplicative_terms']) > 1):
1095            raise Exception('A bug occurred in seasonal components.')
1096        # Compare to the training, if set.
1097        if self.train_component_cols is not None:
1098            component_cols = component_cols[self.train_component_cols.columns]
1099            if not component_cols.equals(self.train_component_cols):
1100                raise Exception('A bug occurred in constructing regressors.')
1101        return component_cols, modes
1102
1103    def add_group_component(
1104        self,
1105        components: pd.DataFrame,
1106        name: str,
1107        group: list[str] | np.ndarray,
1108    ) -> pd.DataFrame:
1109        """Adds a component with given name that contains all of the components
1110        in group.
1111
1112        Parameters
1113        ----------
1114        components: Dataframe with components.
1115        name: Name of new group component.
1116        group: List of components that form the group.
1117
1118        Returns
1119        -------
1120        Dataframe with components.
1121        """
1122        new_comp = components[components['component'].isin(set(group))].copy()
1123        group_cols = new_comp['col'].unique()
1124        if len(group_cols) > 0:
1125            new_comp = pd.DataFrame({'col': group_cols, 'component': name})
1126            components = pd.concat([components, new_comp], ignore_index=True)
1127        return components
1128
1129    def parse_seasonality_args(
1130        self,
1131        name: str,
1132        arg: Literal['auto'] | int,
1133        auto_disable: bool,
1134        default_order: int,
1135    ) -> int:
1136        """Get number of fourier components for built-in seasonalities.
1137
1138        Parameters
1139        ----------
1140        name: string name of the seasonality component.
1141        arg: 'auto', True, False, or number of fourier components as provided.
1142        auto_disable: bool if seasonality should be disabled when 'auto'.
1143        default_order: int default fourier order
1144
1145        Returns
1146        -------
1147        Number of fourier components, or 0 for disabled.
1148        """
1149        if arg == 'auto':
1150            fourier_order = 0
1151            if name in self.seasonalities:
1152                logger.info(
1153                    'Found custom seasonality named {name!r}, disabling '
1154                    'built-in {name!r} seasonality.'.format(name=name)
1155                )
1156            elif auto_disable:
1157                logger.info(
1158                    'Disabling {name} seasonality. Run prophet with '
1159                    '{name}_seasonality=True to override this.'
1160                    .format(name=name)
1161                )
1162            else:
1163                fourier_order = default_order
1164        elif arg is True:
1165            fourier_order = default_order
1166        elif arg is False:
1167            fourier_order = 0
1168        else:
1169            fourier_order = int(arg)
1170        return fourier_order
1171
1172    def set_auto_seasonalities(self) -> None:
1173        """Set seasonalities that were left on auto.
1174
1175        Turns on yearly seasonality if there is >=2 years of history.
1176        Turns on weekly seasonality if there is >=2 weeks of history, and the
1177        spacing between dates in the history is <7 days.
1178        Turns on daily seasonality if there is >=2 days of history, and the
1179        spacing between dates in the history is <1 day.
1180        """
1181        history = cast(pd.DataFrame, self.history)
1182        first = history['ds'].min()
1183        last = history['ds'].max()
1184        dt = history['ds'].diff()
1185        min_dt = dt.iloc[cast(np.ndarray, dt.values).nonzero()[0]].min()
1186
1187        # Yearly seasonality
1188        yearly_disable = last - first < pd.Timedelta(days=730)
1189        fourier_order = self.parse_seasonality_args(
1190            'yearly', self.yearly_seasonality, yearly_disable, 10)
1191        if fourier_order > 0:
1192            self.seasonalities['yearly'] = {
1193                'period': 365.25,
1194                'fourier_order': fourier_order,
1195                'prior_scale': self.seasonality_prior_scale,
1196                'mode': self.seasonality_mode,
1197                'condition_name': None
1198            }
1199
1200        # Weekly seasonality
1201        weekly_disable = ((last - first < pd.Timedelta(weeks=2)) or
1202                          (min_dt >= pd.Timedelta(weeks=1)))  # pyrefly:ignore[unsupported-operation]
1203        fourier_order = self.parse_seasonality_args(
1204            'weekly', self.weekly_seasonality, weekly_disable, 3)
1205        if fourier_order > 0:
1206            self.seasonalities['weekly'] = {
1207                'period': 7,
1208                'fourier_order': fourier_order,
1209                'prior_scale': self.seasonality_prior_scale,
1210                'mode': self.seasonality_mode,
1211                'condition_name': None
1212            }
1213
1214        # Daily seasonality
1215        daily_disable = ((last - first < pd.Timedelta(days=2)) or
1216                         (min_dt >= pd.Timedelta(days=1)))  # pyrefly:ignore[unsupported-operation]
1217        fourier_order = self.parse_seasonality_args(
1218            'daily', self.daily_seasonality, daily_disable, 4)
1219        if fourier_order > 0:
1220            self.seasonalities['daily'] = {
1221                'period': 1,
1222                'fourier_order': fourier_order,
1223                'prior_scale': self.seasonality_prior_scale,
1224                'mode': self.seasonality_mode,
1225                'condition_name': None
1226            }
1227
1228    @staticmethod
1229    def linear_growth_init(df: pd.DataFrame) -> tuple[float, float]:
1230        """Initialize linear growth.
1231
1232        Provides a strong initialization for linear growth by calculating the
1233        growth and offset parameters that pass the function through the first
1234        and last points in the time series.
1235
1236        Parameters
1237        ----------
1238        df: pd.DataFrame with columns ds (date), y_scaled (scaled time series),
1239            and t (scaled time).
1240
1241        Returns
1242        -------
1243        A tuple (k, m) with the rate (k) and offset (m) of the linear growth
1244        function.
1245        """
1246        i0, i1 = cast(int, df['ds'].idxmin()), cast(int, df['ds'].idxmax())
1247        T = df['t'].iloc[i1] - df['t'].iloc[i0]
1248        k = (df['y_scaled'].iloc[i1] - df['y_scaled'].iloc[i0]) / T
1249        m = df['y_scaled'].iloc[i0] - k * df['t'].iloc[i0]
1250        return (k, m)
1251
1252    @staticmethod
1253    def logistic_growth_init(df: pd.DataFrame) -> tuple[float, float]:
1254        """Initialize logistic growth.
1255
1256        Provides a strong initialization for logistic growth by calculating the
1257        growth and offset parameters that pass the function through the first
1258        and last points in the time series.
1259
1260        Parameters
1261        ----------
1262        df: pd.DataFrame with columns ds (date), cap_scaled (scaled capacity),
1263            y_scaled (scaled time series), and t (scaled time).
1264
1265        Returns
1266        -------
1267        A tuple (k, m) with the rate (k) and offset (m) of the logistic growth
1268        function.
1269        """
1270        i0, i1 = cast(int, df['ds'].idxmin()), cast(int, df['ds'].idxmax())
1271        T = df['t'].iloc[i1] - df['t'].iloc[i0]
1272
1273        # Force valid values, in case y > cap or y < 0
1274        C0 = df['cap_scaled'].iloc[i0]
1275        C1 = df['cap_scaled'].iloc[i1]
1276        y0 = max(0.01 * C0, min(0.99 * C0, df['y_scaled'].iloc[i0]))
1277        y1 = max(0.01 * C1, min(0.99 * C1, df['y_scaled'].iloc[i1]))
1278
1279        r0 = C0 / y0
1280        r1 = C1 / y1
1281
1282        if abs(r0 - r1) <= 0.01:
1283            r0 = 1.05 * r0
1284
1285        L0 = np.log(r0 - 1)
1286        L1 = np.log(r1 - 1)
1287
1288        # Initialize the offset
1289        m = L0 * T / (L0 - L1)
1290        # And the rate
1291        k = (L0 - L1) / T
1292        return (k, m)
1293
1294    @staticmethod
1295    def flat_growth_init(df: pd.DataFrame) -> tuple[float, float]:
1296        """Initialize flat growth.
1297
1298        Provides a strong initialization for flat growth. Sets the growth to 0
1299        and offset parameter as mean of history y_scaled values.
1300
1301        Parameters
1302        ----------
1303        df: pd.DataFrame with columns ds (date), y_scaled (scaled time series),
1304            and t (scaled time).
1305
1306        Returns
1307        -------
1308        A tuple (k, m) with the rate (k) and offset (m) of the linear growth
1309        function.
1310        """
1311        k = 0
1312        m = df['y_scaled'].mean()
1313        return k, m
1314
1315    def preprocess(self, df: pd.DataFrame, **kwargs: Any) -> ModelInputData:
1316        """
1317        Reformats historical data, standardizes y and extra regressors, sets seasonalities and changepoints.
1318
1319        Saves the preprocessed data to the instantiated object, and also returns the relevant components
1320        as a ModelInputData object.
1321        """
1322        if ('ds' not in df) or ('y' not in df):
1323            raise ValueError(
1324                'Dataframe must have columns "ds" and "y" with the dates and '
1325                'values respectively.'
1326            )
1327        history = df[df['y'].notnull()].copy()
1328        if history.shape[0] < 2:
1329            raise ValueError('Dataframe has less than 2 non-NaN rows.')
1330        self.history_dates = pd.to_datetime(pd.Series(df['ds'].unique(), name='ds')).sort_values()
1331
1332        self.history = self.setup_dataframe(history, initialize_scales=True)
1333        self.set_auto_seasonalities()
1334        seasonal_features, prior_scales, component_cols, modes = (
1335            self.make_all_seasonality_features(self.history))
1336        self.train_component_cols = component_cols
1337        self.component_modes = modes
1338        self.fit_kwargs = deepcopy(kwargs)
1339
1340        self.set_changepoints()
1341
1342        if self.growth in ['linear', 'flat']:
1343            cap = np.zeros(self.history.shape[0])
1344        else:
1345            cap = self.history['cap_scaled']
1346
1347        return ModelInputData(
1348            T=self.history.shape[0],
1349            S=len(cast(npt.NDArray[np.float64], self.changepoints_t)),
1350            K=seasonal_features.shape[1],
1351            tau=self.changepoint_prior_scale,
1352            trend_indicator=TrendIndicator[self.growth.upper()].value,
1353            y=self.history['y_scaled'],
1354            t=self.history['t'],
1355            t_change=cast(npt.NDArray[np.float64], self.changepoints_t),
1356            X=seasonal_features,
1357            sigmas=prior_scales,
1358            s_a=component_cols['additive_terms'],
1359            s_m=component_cols['multiplicative_terms'],
1360            cap=cap,
1361        )
1362
1363    def calculate_initial_params(self, num_total_regressors: int) -> ModelParams:
1364        """
1365        Calculates initial parameters for the model based on the preprocessed history.
1366
1367        Parameters
1368        ----------
1369        num_total_regressors: the count of seasonality fourier components plus holidays plus extra regressors.
1370        """
1371        history = cast(pd.DataFrame, self.history)
1372        if self.growth == 'linear':
1373            k, m = self.linear_growth_init(history)
1374        elif self.growth == 'flat':
1375            k, m = self.flat_growth_init(history)
1376        else:
1377            assert self.growth == "logistic"
1378            k, m = self.logistic_growth_init(history)
1379        return ModelParams(
1380            k=k,
1381            m=m,
1382            delta=np.zeros_like(self.changepoints_t),
1383            beta=np.zeros(num_total_regressors),
1384            sigma_obs=1.0,
1385        )
1386
1387    def fit(self, df: pd.DataFrame, **kwargs: Any) -> Self:
1388        """Fit the Prophet model.
1389
1390        This sets self.params to contain the fitted model parameters. It is a
1391        dictionary parameter names as keys and the following items:
1392            k (Mx1 array): M posterior samples of the initial slope.
1393            m (Mx1 array): The initial intercept.
1394            delta (MxN array): The slope change at each of N changepoints.
1395            beta (MxK matrix): Coefficients for K seasonality features.
1396            sigma_obs (Mx1 array): Noise level.
1397        Note that M=1 if MAP estimation.
1398
1399        Parameters
1400        ----------
1401        df: pd.DataFrame containing the history. Must have columns ds (date
1402            type) and y, the time series. If self.growth is 'logistic', then
1403            df must also have a column cap that specifies the capacity at
1404            each ds.
1405        kwargs: Additional arguments passed to the optimizing or sampling
1406            functions in Stan.
1407
1408        Returns
1409        -------
1410        The fitted Prophet object.
1411        """
1412        if self.history is not None:
1413            raise Exception('Prophet object can only be fit once. '
1414                            'Instantiate a new object.')
1415
1416        model_inputs = self.preprocess(df, **kwargs)
1417        self._fit_regressor_models(df)
1418        initial_params = self.calculate_initial_params(model_inputs.K)
1419
1420        dat = dataclasses.asdict(model_inputs)
1421        stan_init = dataclasses.asdict(initial_params)
1422        stan_backend = cast(IStanBackend, self.stan_backend)
1423
1424        history = cast(pd.DataFrame, self.history)
1425        if history['y'].min() == history['y'].max() and \
1426                (self.growth == 'linear' or self.growth == 'flat'):
1427            self.params = stan_init
1428            self.params['sigma_obs'] = 1e-9
1429            for par in self.params:
1430                self.params[par] = np.array([self.params[par]])
1431        elif self.mcmc_samples > 0:
1432            self.params = stan_backend.sampling(stan_init, dat, self.mcmc_samples, **kwargs)
1433        else:
1434            self.params = stan_backend.fit(stan_init, dat, **kwargs)
1435
1436        self.stan_fit = stan_backend.stan_fit
1437        # If no changepoints were requested, replace delta with 0s
1438        if len(cast(pd.Series, self.changepoints)) == 0:
1439            # Fold delta into the base rate k
1440            self.params['k'] = (
1441                self.params['k'] + self.params['delta'].reshape(-1)
1442            )
1443            self.params['delta'] = (np.zeros(self.params['delta'].shape)
1444                                      .reshape((-1, 1)))
1445
1446        return self
1447
1448    def predict(self, df: pd.DataFrame | None = None, vectorized: bool = True) -> pd.DataFrame:
1449        """Predict using the prophet model.
1450
1451        Parameters
1452        ----------
1453        df: pd.DataFrame with dates for predictions (column ds), and capacity
1454            (column cap) if logistic growth. If not provided, predictions are
1455            made on the history.
1456        vectorized: Whether to use a vectorized method to compute uncertainty intervals. Suggest using
1457            True (the default) for much faster runtimes in most cases,
1458            except when (growth = 'logistic' and mcmc_samples > 0).
1459
1460        Returns
1461        -------
1462        A pd.DataFrame with the forecast components.
1463        """
1464        if self.history is None:
1465            raise Exception('Model has not been fit.')
1466
1467        regressor_samples = None
1468        if df is None:
1469            df = self.history.copy()
1470        else:
1471            if df.shape[0] == 0:
1472                raise ValueError('Dataframe has no rows.')
1473            df, regressor_samples = self._prepare_regressors_for_predict(
1474                df.copy(),
1475                self.uncertainty_samples if self.uncertainty_samples else None,
1476            )
1477            df = self.setup_dataframe(df)
1478
1479        if regressor_samples:
1480            standardized_samples = {}
1481            for name, samples in regressor_samples.items():
1482                props = self.extra_regressors[name]
1483                standardized_samples[name] = (samples - props['mu']) / props['std']
1484            regressor_samples = standardized_samples
1485
1486        df['trend'] = self.predict_trend(df)
1487        seasonal_components = self.predict_seasonal_components(df)
1488        if self.uncertainty_samples:
1489            intervals = self.predict_uncertainty(df, vectorized, regressor_samples)
1490        else:
1491            intervals = None
1492
1493        # Drop columns except ds, cap, floor, and trend
1494        cols = ['ds', 'trend']
1495        if 'cap' in df:
1496            cols.append('cap')
1497        if self.logistic_floor:
1498            cols.append('floor')
1499        # Add in forecast components
1500        df2 = pd.concat((df[cols], intervals, seasonal_components), axis=1)
1501        df2['yhat'] = (
1502                df2['trend'] * (1 + df2['multiplicative_terms'])
1503                + df2['additive_terms']
1504        )
1505        return df2
1506
1507    @staticmethod
1508    def piecewise_linear(
1509        t: np.ndarray,
1510        deltas: np.ndarray,
1511        k: float,
1512        m: float,
1513        changepoint_ts: np.ndarray,
1514    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1515        """Evaluate the piecewise linear function.
1516
1517        Parameters
1518        ----------
1519        t: np.array of times on which the function is evaluated.
1520        deltas: np.array of rate changes at each changepoint.
1521        k: Float initial rate.
1522        m: Float initial offset.
1523        changepoint_ts: np.array of changepoint times.
1524
1525        Returns
1526        -------
1527        Vector y(t).
1528        """
1529        deltas_t = (changepoint_ts[None, :] <= t[..., None]) * deltas
1530        k_t = deltas_t.sum(axis=1) + k
1531        m_t = (deltas_t * -changepoint_ts).sum(axis=1) + m
1532        return k_t * t + m_t
1533
1534    @staticmethod
1535    def piecewise_logistic(
1536        t: np.ndarray,
1537        cap: np.ndarray | pd.Series,
1538        deltas: np.ndarray,
1539        k: float,
1540        m: float,
1541        changepoint_ts: np.ndarray,
1542    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1543        """Evaluate the piecewise logistic function.
1544
1545        Parameters
1546        ----------
1547        t: np.array of times on which the function is evaluated.
1548        cap: np.array of capacities at each t.
1549        deltas: np.array of rate changes at each changepoint.
1550        k: Float initial rate.
1551        m: Float initial offset.
1552        changepoint_ts: np.array of changepoint times.
1553
1554        Returns
1555        -------
1556        Vector y(t).
1557        """
1558        # Compute offset changes
1559        # Ensure k and m are scalars for numpy 2.x compatibility
1560        k_scalar: float = np.asarray(k).item() if np.asarray(k).size == 1 else k
1561        m_scalar: float = np.asarray(m).item() if np.asarray(m).size == 1 else m
1562        k_cum = np.concatenate((np.atleast_1d(k_scalar), np.cumsum(deltas) + k_scalar))
1563        gammas = np.zeros(len(changepoint_ts))
1564        for i, t_s in enumerate(changepoint_ts):
1565            gammas[i] = (
1566                    (t_s - m_scalar - np.sum(gammas))
1567                    * (1 - k_cum[i] / k_cum[i + 1])  # noqa W503
1568            )
1569        # Get cumulative rate and offset at each t
1570        k_t = k_scalar * np.ones_like(t)
1571        m_t = m_scalar * np.ones_like(t)
1572        for s, t_s in enumerate(changepoint_ts):
1573            indx = t >= t_s
1574            k_t[indx] += deltas[s]
1575            m_t[indx] += gammas[s]
1576        return cap / (1 + np.exp(-k_t * (t - m_t)))
1577
1578    @staticmethod
1579    def flat_trend(
1580        t: np.ndarray,
1581        m: float,
1582    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1583        """Evaluate the flat trend function.
1584
1585        Parameters
1586        ----------
1587        t: np.array of times on which the function is evaluated.
1588        m: Float initial offset.
1589
1590        Returns
1591        -------
1592        Vector y(t).
1593        """
1594        m_t = m * np.ones_like(t)
1595        return m_t
1596
1597    def predict_trend(
1598        self,
1599        df: pd.DataFrame,
1600    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1601        """Predict trend using the prophet model.
1602
1603        Parameters
1604        ----------
1605        df: Prediction dataframe.
1606
1607        Returns
1608        -------
1609        Vector with trend on prediction dates.
1610        """
1611        k = np.nanmean(self.params['k'])
1612        m = np.nanmean(self.params['m'])
1613        deltas = np.nanmean(self.params['delta'], axis=0)
1614
1615        t = np.array(df['t'])
1616        changepoints_t = cast(npt.NDArray[np.float64], self.changepoints_t)
1617        if self.growth == 'linear':
1618            trend = self.piecewise_linear(t, deltas, k, m, changepoints_t)
1619        elif self.growth == 'logistic':
1620            cap = df['cap_scaled']
1621            trend = self.piecewise_logistic(t, cap, deltas, k, m, changepoints_t)
1622        else:
1623            # constant trend
1624            assert self.growth == "flat"
1625            trend = self.flat_trend(t, m)
1626
1627        return trend * cast(float, self.y_scale) + df['floor']
1628
1629    def predict_seasonal_components(self, df: pd.DataFrame) -> pd.DataFrame:
1630        """Predict seasonality components, holidays, and added regressors.
1631
1632        Parameters
1633        ----------
1634        df: Prediction dataframe.
1635
1636        Returns
1637        -------
1638        Dataframe with seasonal components.
1639        """
1640        seasonal_features, _, component_cols, _ = (
1641            self.make_all_seasonality_features(df)
1642        )
1643        if self.uncertainty_samples:
1644            lower_p = 100 * (1.0 - self.interval_width) / 2
1645            upper_p = 100 * (1.0 + self.interval_width) / 2
1646
1647        X = seasonal_features.values
1648        data = {}
1649        for component in component_cols.columns:
1650            beta_c = self.params['beta'] * component_cols[component].values
1651
1652            comp = np.matmul(X, beta_c.transpose())
1653            assert self.component_modes is not None
1654            if component in self.component_modes['additive']:
1655                comp *= self.y_scale
1656            data[component] = np.nanmean(comp, axis=1)
1657            if self.uncertainty_samples:
1658                # pyrefly:ignore[unbound-name]
1659                data[component + '_lower'] = self.percentile(comp, lower_p, axis=1)
1660                # pyrefly:ignore[unbound-name]
1661                data[component + '_upper'] = self.percentile(comp, upper_p, axis=1)
1662        return pd.DataFrame(data)
1663
1664    def predict_uncertainty(
1665        self,
1666        df: pd.DataFrame,
1667        vectorized: bool,
1668        regressor_samples: dict[str, npt.NDArray[np.float64]] | None = None,
1669    ) -> pd.DataFrame:
1670        """Prediction intervals for yhat and trend.
1671
1672        Parameters
1673        ----------
1674        df: Prediction dataframe.
1675        vectorized: Whether to use a vectorized method for generating future draws.
1676        regressor_samples: Optional draws for regressor values (already standardized)
1677            aligned with df.
1678
1679        Returns
1680        -------
1681        Dataframe with uncertainty intervals.
1682        """
1683        sim_values = self.sample_posterior_predictive(df, vectorized, regressor_samples)
1684
1685        lower_p = 100 * (1.0 - self.interval_width) / 2
1686        upper_p = 100 * (1.0 + self.interval_width) / 2
1687
1688        series = {}
1689        for key in ['yhat', 'trend']:
1690            series['{}_lower'.format(key)] = self.percentile(
1691                sim_values[key], lower_p, axis=1)
1692            series['{}_upper'.format(key)] = self.percentile(
1693                sim_values[key], upper_p, axis=1)
1694
1695        return pd.DataFrame(series)
1696
1697    def sample_posterior_predictive(
1698        self,
1699        df: pd.DataFrame,
1700        vectorized: bool,
1701        regressor_samples: dict[str, np.ndarray] | None = None,
1702    ) -> dict[str, npt.NDArray[np.float64]]:
1703        """Prophet posterior predictive samples.
1704
1705        Parameters
1706        ----------
1707        df: Prediction dataframe.
1708        vectorized: Whether to use a vectorized method to generate future draws.
1709        regressor_samples: Optional draws for regressors (standardized) keyed by
1710            regressor name. If provided, regressor uncertainty will be propagated
1711            through forecast draws.
1712
1713        Returns
1714        -------
1715        Dictionary with posterior predictive samples for the forecast yhat and
1716        for the trend component.
1717        """
1718        n_iterations = self.params['k'].shape[0]
1719        samp_per_iter = max(1, int(np.ceil(
1720            self.uncertainty_samples / float(n_iterations)
1721        )))
1722        # Generate seasonality features once so we can re-use them.
1723        seasonal_features, _, component_cols, _ = (
1724            self.make_all_seasonality_features(df)
1725        )
1726        sim_values = {'yhat': [], 'trend': []}
1727        regressor_positions: dict[str, int] = {}
1728        regressor_draw_counts: dict[str, int] = {}
1729        if regressor_samples:
1730            vectorized = False
1731            for name in regressor_samples:
1732                if name not in seasonal_features.columns:
1733                    continue
1734                pos = seasonal_features.columns.get_loc(name)
1735                if not isinstance(pos, (int, np.integer)):
1736                    raise ValueError(
1737                        f"Expected a unique column position for regressor '{name}'"
1738                    )
1739                regressor_positions[name] = int(pos)
1740                regressor_draw_counts[name] = regressor_samples[name].shape[1]
1741        sample_counter = 0
1742        for i in range(n_iterations):
1743            if vectorized:
1744                sims = self.sample_model_vectorized(
1745                    df=df,
1746                    seasonal_features=seasonal_features,
1747                    iteration=i,
1748                    s_a=component_cols['additive_terms'],
1749                    s_m=component_cols['multiplicative_terms'],
1750                    n_samples=samp_per_iter
1751                )
1752            else:
1753                sims = []
1754                for _ in range(samp_per_iter):
1755                    seasonal_features_sample = seasonal_features
1756                    if regressor_samples:
1757                        seasonal_features_sample = seasonal_features.copy()
1758                        for name, pos in regressor_positions.items():
1759                            draw_count = regressor_draw_counts[name]
1760                            use_idx = sample_counter % draw_count
1761                            seasonal_features_sample.iloc[:, int(pos)] = (
1762                                regressor_samples[name][:, use_idx]
1763                            )
1764                    sims.append(
1765                        self.sample_model(
1766                            df=df,
1767                            seasonal_features=seasonal_features_sample,
1768                            iteration=i,
1769                            s_a=component_cols['additive_terms'],
1770                            s_m=component_cols['multiplicative_terms'],
1771                        )
1772                    )
1773                    sample_counter += 1
1774            for key in sim_values:
1775                for sim in sims:
1776                    sim_values[key].append(sim[key])
1777
1778        return {k: np.column_stack(v) for k, v in sim_values.items()}
1779
1780    def sample_model(
1781        self,
1782        df: pd.DataFrame,
1783        seasonal_features: pd.DataFrame,
1784        iteration: int,
1785        s_a: pd.Series,
1786        s_m: pd.Series,
1787    ) -> dict[str, npt.NDArray[np.float64]]:
1788        """Simulate observations from the extrapolated generative model.
1789
1790        Parameters
1791        ----------
1792        df: Prediction dataframe.
1793        seasonal_features: pd.DataFrame of seasonal features.
1794        iteration: Int sampling iteration to use parameters from.
1795        s_a: Indicator vector for additive components
1796        s_m: Indicator vector for multiplicative components
1797
1798        Returns
1799        -------
1800        Dictionary with `yhat` and `trend`, each like df['t'].
1801        """
1802        trend = self.sample_predictive_trend(df, iteration)
1803
1804        y_scale = cast(float, self.y_scale)
1805        beta = self.params['beta'][iteration]
1806        Xb_a = np.matmul(seasonal_features.values,
1807                         beta * s_a.values) * y_scale
1808        Xb_m = np.matmul(seasonal_features.values, beta * s_m.values)
1809
1810        sigma = self.params['sigma_obs'][iteration]
1811        noise = np.random.normal(0, sigma, df.shape[0]) * y_scale
1812
1813        return {
1814            'yhat': trend * (1 + Xb_m) + Xb_a + noise,
1815            'trend': trend
1816        }
1817
1818    def sample_model_vectorized(
1819        self,
1820        df: pd.DataFrame,
1821        seasonal_features: pd.DataFrame,
1822        iteration: int,
1823        s_a: pd.Series,
1824        s_m: pd.Series,
1825        n_samples: int,
1826    ) -> list[dict[str, np.ndarray]]:
1827        """Simulate observations from the extrapolated generative model. Vectorized version of sample_model().
1828
1829        Parameters
1830        ----------
1831        df: Prediction dataframe.
1832        seasonal_features: pd.DataFrame of seasonal features.
1833        iteration: Int sampling iteration to use parameters from.
1834        s_a: Indicator vector for additive components.
1835        s_m: Indicator vector for multiplicative components.
1836        n_samples: Number of future paths of the trend to simulate.
1837
1838        Returns
1839        -------
1840        List (length n_samples) of dictionaries with arrays for trend and yhat, each ordered like df['t'].
1841        """
1842        # Get the seasonality and regressor components, which are deterministic per iteration
1843        beta = self.params['beta'][iteration]
1844        Xb_a = np.matmul(seasonal_features.values,
1845                        beta * s_a.values) * self.y_scale
1846        Xb_m = np.matmul(seasonal_features.values, beta * s_m.values)
1847        # Get the future trend, which is stochastic per iteration
1848        trends = self.sample_predictive_trend_vectorized(df, n_samples, iteration)  # already on the same scale as the actual data
1849        sigma = self.params['sigma_obs'][iteration]
1850        noise_terms = np.random.normal(0, sigma, trends.shape) * cast(float, self.y_scale)
1851
1852        simulations = []
1853        for trend, noise in zip(trends, noise_terms):
1854            simulations.append({
1855                'yhat': trend * (1 + Xb_m) + Xb_a + noise,
1856                'trend': trend
1857            })
1858        return simulations
1859
1860    def sample_predictive_trend(self, df: pd.DataFrame, iteration: int) -> np.ndarray:
1861        """Simulate the trend using the extrapolated generative model.
1862
1863        Parameters
1864        ----------
1865        df: Prediction dataframe.
1866        iteration: Int sampling iteration to use parameters from.
1867
1868        Returns
1869        -------
1870        np.array of simulated trend over df['t'].
1871        """
1872        k = self.params['k'][iteration]
1873        m = self.params['m'][iteration]
1874        deltas = self.params['delta'][iteration]
1875
1876        t = np.array(df['t'])
1877        T = t.max()
1878
1879        # New changepoints from a Poisson process with rate S on [1, T]
1880        changepoints_t = cast(np.ndarray, self.changepoints_t)
1881        if T > 1:
1882            S = len(changepoints_t)
1883            n_changes = np.random.poisson(S * (T - 1))
1884        else:
1885            n_changes = 0
1886        if n_changes > 0:
1887            changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
1888            changepoint_ts_new.sort()
1889        else:
1890            changepoint_ts_new = []
1891
1892        # Get the empirical scale of the deltas, plus epsilon to avoid NaNs.
1893        lambda_ = np.mean(np.abs(deltas)) + 1e-8
1894
1895        # Sample deltas
1896        deltas_new = np.random.laplace(0, lambda_, n_changes)
1897
1898        # Prepend the times and deltas from the history
1899        changepoint_ts = np.concatenate((changepoints_t,
1900                                         changepoint_ts_new))
1901        deltas = np.concatenate((deltas, deltas_new))
1902
1903        if self.growth == 'linear':
1904            trend = self.piecewise_linear(t, deltas, k, m, changepoint_ts)
1905        elif self.growth == 'logistic':
1906            cap = df['cap_scaled']
1907            trend = self.piecewise_logistic(t, cap, deltas, k, m,
1908                                            changepoint_ts)
1909        else:
1910            assert self.growth == "flat"
1911            trend = self.flat_trend(t, m)
1912
1913        return trend * cast(float, self.y_scale) + df['floor']
1914
1915    def sample_predictive_trend_vectorized(
1916        self,
1917        df: pd.DataFrame,
1918        n_samples: int,
1919        iteration: int = 0,
1920    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
1921        """Sample draws of the future trend values. Vectorized version of sample_predictive_trend().
1922
1923        Parameters
1924        ----------
1925        df: Prediction dataframe.
1926        iteration: Int sampling iteration to use parameters from.
1927        n_samples: Number of future paths of the trend to simulate.
1928
1929        Returns
1930        -------
1931        Draws of the trend values with shape (n_samples, len(df)). Values are on the scale of the original data.
1932        """
1933        deltas = self.params["delta"][iteration]
1934        m = self.params["m"][iteration]
1935        k = self.params["k"][iteration]
1936        changepoints_t = cast(np.ndarray, self.changepoints_t)
1937        if self.growth == "linear":
1938            expected = self.piecewise_linear(
1939                cast(np.ndarray, df["t"].values), deltas, k, m, changepoints_t
1940            )
1941        elif self.growth == "logistic":
1942            expected = self.piecewise_logistic(
1943                cast(np.ndarray, df["t"].values),
1944                cast(np.ndarray, df["cap_scaled"].values),
1945                deltas,
1946                k,
1947                m,
1948                changepoints_t,
1949            )
1950        elif self.growth == "flat":
1951            expected = self.flat_trend(cast(np.ndarray, df["t"].values), m)
1952        else:
1953            raise NotImplementedError
1954        uncertainty = self._sample_uncertainty(df, n_samples, iteration)
1955        return (
1956            (np.tile(expected, (n_samples, 1)) + uncertainty) * cast(float, self.y_scale) +
1957            np.tile(cast(np.ndarray, df["floor"].values), (n_samples, 1))
1958        )
1959
1960    def _sample_uncertainty(
1961        self,
1962        df: pd.DataFrame,
1963        n_samples: int,
1964        iteration: int = 0,
1965    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
1966        """Sample draws of future trend changes, vectorizing as much as possible.
1967
1968        Parameters
1969        ----------
1970        df: pd.DataFrame with columns `t` (time scaled to the model context), trend, and cap.
1971        n_samples: Number of future paths of the trend to simulate
1972        iteration: The iteration of the parameter set to use. Default 0, the first iteration.
1973
1974        Returns
1975        -------
1976        Draws of the trend changes with shape (n_samples, len(df)). Values are standardized.
1977        """
1978        # handle only historic data
1979        if df["t"].max() <= 1:
1980            # there is no trend uncertainty in historic trends
1981            uncertainties = np.zeros((n_samples, len(df)))
1982        else:
1983            future_df = df.loc[df["t"] > 1]
1984            n_length = len(future_df)
1985            # handle 1 length futures by using history
1986            if n_length > 1:
1987                single_diff = np.diff(future_df["t"]).mean()
1988            else:
1989                single_diff = np.diff(cast(pd.DataFrame, self.history)["t"]).mean()
1990            change_likelihood = len(cast(np.ndarray, self.changepoints_t)) * single_diff
1991            deltas = self.params["delta"][iteration]
1992            m = self.params["m"][iteration]
1993            k = self.params["k"][iteration]
1994            mean_delta = np.mean(np.abs(deltas)) + 1e-8
1995            if self.growth == "linear":
1996                mat = self._make_trend_shift_matrix(mean_delta, change_likelihood, n_length, n_samples=n_samples)
1997                uncertainties = mat.cumsum(axis=1).cumsum(axis=1)  # from slope changes to actual values
1998                uncertainties *= single_diff  # scaled by the actual meaning of the slope
1999            elif self.growth == "logistic":
2000                mat = self._make_trend_shift_matrix(mean_delta, change_likelihood, n_length, n_samples=n_samples)
2001                uncertainties = self._logistic_uncertainty(
2002                    mat=mat,
2003                    deltas=deltas,
2004                    k=k,
2005                    m=m,
2006                    cap=cast(np.ndarray, future_df["cap_scaled"].values),
2007                    t_time=cast(np.ndarray, future_df["t"].values),
2008                    n_length=n_length,
2009                    single_diff=single_diff,
2010                )
2011            elif self.growth == "flat":
2012                # no trend uncertainty when there is no growth
2013                uncertainties = np.zeros((n_samples, n_length))
2014            else:
2015                raise NotImplementedError
2016            # handle past included in dataframe
2017            if df["t"].min() <= 1:
2018                past_uncertainty = np.zeros((n_samples, np.sum(df["t"] <= 1)))
2019                uncertainties = np.concatenate([past_uncertainty, uncertainties], axis=1)
2020        return uncertainties
2021
2022    @staticmethod
2023    def _make_trend_shift_matrix(
2024        mean_delta: float, likelihood: float, future_length: int, n_samples: int
2025    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
2026        """
2027        Creates a matrix of random trend shifts based on historical likelihood and size of shifts.
2028        Can be used for either linear or logistic trend shifts.
2029        Each row represents a different sample of a possible future, and each column is a time step into the future.
2030        """
2031        # create a bool matrix of where these trend shifts should go
2032        bool_slope_change = np.random.uniform(size=(n_samples, future_length)) < likelihood
2033        shift_values = np.random.laplace(0, mean_delta, size=bool_slope_change.shape)
2034        mat = shift_values * bool_slope_change
2035        n_mat = np.hstack([np.zeros((len(mat), 1)), mat])[:, :-1]
2036        mat = (n_mat + mat) / 2
2037        return mat
2038
2039    @staticmethod
2040    def _make_historical_mat_time(
2041        deltas: npt.ArrayLike,
2042        changepoints_t: np.ndarray,
2043        t_time: npt.ArrayLike,
2044        n_row: int = 1,
2045        single_diff: float | None = None,
2046    ) -> tuple[
2047        np.ndarray[tuple[int], np.dtype[np.float64]],
2048        np.ndarray[tuple[int], np.dtype[np.float64]],
2049    ]:
2050        """
2051        Creates a matrix of slope-deltas where these changes occured in training data according to the trained prophet obj
2052        """
2053        if single_diff is None:
2054            single_diff = np.diff(t_time).mean()
2055        prev_time = np.arange(0, 1 + single_diff, single_diff)
2056        idxs = []
2057        for changepoint in changepoints_t:
2058            idxs.append(np.where(prev_time > changepoint)[0][0])
2059        prev_deltas = np.zeros(len(prev_time))
2060        prev_deltas[idxs] = deltas
2061        prev_deltas = np.repeat(prev_deltas.reshape(1, -1), n_row, axis=0)
2062        return prev_deltas, prev_time
2063
2064    def _logistic_uncertainty(
2065        self,
2066        mat: np.ndarray,
2067        deltas: np.ndarray,
2068        k: float,
2069        m: float,
2070        cap: np.ndarray,
2071        t_time: np.ndarray,
2072        n_length: int,
2073        single_diff: float | None = None,
2074    ) -> np.ndarray:
2075        """
2076        Vectorizes prophet's logistic uncertainty by creating a matrix of future possible trends.
2077
2078        Parameters
2079        ----------
2080        mat: A trend shift matrix returned by _make_trend_shift_matrix()
2081        deltas: The size of the trend changes at each changepoint, estimated by the model
2082        k: Float initial rate.
2083        m: Float initial offset.
2084        cap: np.array of capacities at each t.
2085        t_time: The values of t in the model context (i.e. scaled so that anything > 1 represents the future)
2086        n_length: For each path, the number of future steps to simulate
2087        single_diff: The difference between each t step in the model context. Default None, inferred
2088            from t_time.
2089
2090        Returns
2091        -------
2092        A numpy array with shape (n_samples, n_length), representing the width of the uncertainty interval
2093        (standardized, not on the same scale as the actual data values) around 0.
2094        """
2095
2096        def ffill(arr):
2097            mask = arr == 0
2098            idx = np.where(~mask, np.arange(mask.shape[1]), 0)
2099            np.maximum.accumulate(idx, axis=1, out=idx)
2100            return arr[np.arange(idx.shape[0])[:, None], idx]
2101
2102        # for logistic growth we need to evaluate the trend all the way from the start of the train item
2103        historical_mat, historical_time = self._make_historical_mat_time(
2104            deltas, cast(np.ndarray, self.changepoints_t), t_time, len(mat), single_diff
2105        )
2106        mat = np.concatenate([historical_mat, mat], axis=1)
2107        full_t_time = np.concatenate([historical_time, t_time])
2108
2109        # apply logistic growth logic on the slope changes
2110        k_cum = np.concatenate((np.ones((mat.shape[0], 1)) * k, np.where(mat, np.cumsum(mat, axis=1) + k, 0)), axis=1)
2111        k_cum_b = ffill(k_cum)
2112        gammas = np.zeros_like(mat)
2113        for i in range(mat.shape[1]):
2114            x = full_t_time[i] - m - np.sum(gammas[:, :i], axis=1)
2115            ks = 1 - k_cum_b[:, i] / k_cum_b[:, i + 1]
2116            gammas[:, i] = x * ks
2117        # the data before the -n_length is the historical values, which are not needed, so cut the last n_length
2118        k_t = (mat.cumsum(axis=1) + k)[:, -n_length:]
2119        m_t = (gammas.cumsum(axis=1) + m)[:, -n_length:]
2120        sample_trends = cap / (1 + np.exp(-k_t * (t_time - m_t)))
2121        # remove the mean because we only need width of the uncertainty centered around 0
2122        # we will add the width to the main forecast - yhat (which is the mean) - later
2123        return sample_trends - sample_trends.mean(axis=0)
2124
2125    def predictive_samples(self, df: pd.DataFrame, vectorized: bool = True) -> dict[str, npt.NDArray[np.float64]]:
2126        """Sample from the posterior predictive distribution. Returns samples
2127        for the main estimate yhat, and for the trend component. The shape of
2128        each output will be (nforecast x nsamples), where nforecast is the
2129        number of points being forecasted (the number of rows in the input
2130        dataframe) and nsamples is the number of posterior samples drawn.
2131        This is the argument `uncertainty_samples` in the Prophet constructor,
2132        which defaults to 1000.
2133
2134        Parameters
2135        ----------
2136        df: Dataframe with dates for predictions (column ds), and capacity
2137            (column cap) if logistic growth.
2138        vectorized: Whether to use a vectorized method to compute possible draws. Suggest using
2139            True (the default) for much faster runtimes in most cases,
2140            except when (growth = 'logistic' and mcmc_samples > 0).
2141
2142        Returns
2143        -------
2144        Dictionary with keys "trend" and "yhat" containing
2145        posterior predictive samples for that component.
2146        """
2147        regressor_samples = None
2148        df, regressor_samples = self._prepare_regressors_for_predict(
2149            df.copy(),
2150            self.uncertainty_samples if self.uncertainty_samples else None,
2151        )
2152        df = self.setup_dataframe(df)
2153        if regressor_samples:
2154            standardized_samples = {}
2155            for name, samples in regressor_samples.items():
2156                props = self.extra_regressors[name]
2157                standardized_samples[name] = (samples - props['mu']) / props['std']
2158            regressor_samples = standardized_samples
2159        return self.sample_posterior_predictive(df, vectorized, regressor_samples)
2160
2161    def percentile(self, a: npt.ArrayLike, *args: Any, **kwargs: Any) -> np.ndarray:
2162        """
2163        We rely on np.nanpercentile in the rare instances where there
2164        are a small number of bad samples with MCMC that contain NaNs.
2165        However, since np.nanpercentile is far slower than np.percentile,
2166        we only fall back to it if the array contains NaNs. See
2167        https://github.com/facebook/prophet/issues/1310 for more details.
2168        """
2169        fn = np.nanpercentile if np.isnan(a).any() else np.percentile
2170        return fn(a, *args, **kwargs)
2171
2172    def make_future_dataframe(
2173        self, periods: int, freq: str | None = 'D', include_history: bool = True
2174    ) -> pd.DataFrame:
2175        """Simulate the trend using the extrapolated generative model.
2176
2177        Parameters
2178        ----------
2179        periods: Int number of periods to forecast forward.
2180        freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.
2181        include_history: Boolean to include the historical dates in the data
2182            frame for predictions.
2183
2184        Returns
2185        -------
2186        pd.Dataframe that extends forward from the end of self.history for the
2187        requested number of periods.
2188        """
2189        if self.history_dates is None:
2190            raise Exception('Model has not been fit.')
2191        if freq is None:
2192            # taking the tail makes freq inference more reliable
2193            freq = pd.infer_freq(self.history_dates.tail(5))
2194            # returns None if inference failed
2195            if freq is None:
2196                raise Exception('Unable to infer `freq`')
2197        last_date = self.history_dates.max()
2198        dates = pd.date_range(
2199            start=last_date,
2200            periods=periods + 1,  # An extra in case we include start
2201            freq=freq)
2202        dates = dates[dates > last_date]  # Drop start if equals last_date
2203        dates = dates[:periods]  # Return correct number of periods
2204
2205        if include_history:
2206            dates = np.concatenate((np.array(self.history_dates), dates))
2207
2208        return pd.DataFrame({'ds': dates})
2209
2210    def plot(
2211        self,
2212        fcst: pd.DataFrame,
2213        ax: plt.Axes | None = None,
2214        uncertainty: bool = True,
2215        plot_cap: bool = True,
2216        xlabel: str = 'ds',
2217        ylabel: str = 'y',
2218        figsize: tuple[int, int] = (10, 6),
2219        include_legend: bool = False,
2220    ) -> plt.Figure:
2221        """Plot the Prophet forecast.
2222
2223        Parameters
2224        ----------
2225        fcst: pd.DataFrame output of self.predict.
2226        ax: Optional matplotlib axes on which to plot.
2227        uncertainty: Optional boolean to plot uncertainty intervals.
2228        plot_cap: Optional boolean indicating if the capacity should be shown
2229            in the figure, if available.
2230        xlabel: Optional label name on X-axis
2231        ylabel: Optional label name on Y-axis
2232        figsize: Optional tuple width, height in inches.
2233        include_legend: Optional boolean to add legend to the plot.
2234
2235        Returns
2236        -------
2237        A matplotlib figure.
2238        """
2239        return plot(
2240            m=self, fcst=fcst, ax=ax, uncertainty=uncertainty,
2241            plot_cap=plot_cap, xlabel=xlabel, ylabel=ylabel,
2242            figsize=figsize, include_legend=include_legend
2243        )
2244
2245    def plot_components(
2246        self,
2247        fcst: pd.DataFrame,
2248        uncertainty: bool = True,
2249        plot_cap: bool = True,
2250        weekly_start: int = 0,
2251        yearly_start: int = 0,
2252        figsize: tuple[int, int] | None = None,
2253    ) -> plt.Figure:
2254        """Plot the Prophet forecast components.
2255
2256        Will plot whichever are available of: trend, holidays, weekly
2257        seasonality, and yearly seasonality.
2258
2259        Parameters
2260        ----------
2261        fcst: pd.DataFrame output of self.predict.
2262        uncertainty: Optional boolean to plot uncertainty intervals.
2263        plot_cap: Optional boolean indicating if the capacity should be shown
2264            in the figure, if available.
2265        weekly_start: Optional int specifying the start day of the weekly
2266            seasonality plot. 0 (default) starts the week on Sunday. 1 shifts
2267            by 1 day to Monday, and so on.
2268        yearly_start: Optional int specifying the start day of the yearly
2269            seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts
2270            by 1 day to Jan 2, and so on.
2271        figsize: Optional tuple width, height in inches.
2272
2273        Returns
2274        -------
2275        A matplotlib figure.
2276        """
2277        return plot_components(
2278            m=self, fcst=fcst, uncertainty=uncertainty, plot_cap=plot_cap,
2279            weekly_start=weekly_start, yearly_start=yearly_start,
2280            figsize=figsize
2281        )