Skip to content

Data Pipeline Code Walkthrough

This tutorial walks through each Python module in the Mirabelle data pipeline. Read it alongside the source files in files/.

Pipeline Overview

The pipeline processes data in sequential stages, each reading from and writing to the scenario PostgreSQL schema:

Source DB (Neon remote)
    ↓ ETL (files/etl/etl.py)
scenario.demand_actuals
    ↓ Characterization (files/characterization/characterization.py)
scenario.series_characteristics
    ↓ Forecast step (run_pipeline.py --only forecast)
        ① Outlier Detection (pre-step, embedded)
           → scenario.demand_corrected  [corrections per pipeline_id]
        ② Forecasting (files/forecasting/statistical_models.py + neural/foundation/ml)
scenario.forecast_results
    ↓ Backtesting (files/evaluation/metrics.py)
scenario.backtest_results
    ↓ Best Method Selection (files/selection/best_method.py)
scenario.series_best_methods
    ↓ Distribution Fitting (files/distribution/fitting.py)
scenario.fitted_distributions
    ↓ MEIO Optimization (files/meio_runner.py + Rust .pyd)
scenario.meio_results

Parallelization for the forecasting stage is handled by the orchestrator (files/utils/orchestrator.py) which distributes series across Dask workers.

ETL — files/etl/etl.py

Purpose

Extracts demand actuals from a remote source database (Neon PostgreSQL), transforms them to a standardised weekly format, and loads them into the local scenario.demand_actuals table.

Architecture

class ETLPipeline:
    """
    Two-database ETL.
    Source DB: dp_plan.calc_dp_actual (Neon remote)
    Target DB: scenario.demand_actuals (local)
    """

If no source DB is configured, the pipeline falls back to self-reload mode (re-aggregates from the local DB itself).

Key Method: run()

The main entry point orchestrates the three phases:

def run(self):
    # 1. Extract: read from source DB with configurable SQL
    raw_df = self._extract()

    # 2. Transform: aggregate to target frequency, create unique_id
    transformed_df = self._transform(raw_df)

    # 3. Load: upsert into scenario.demand_actuals
    n_rows = self._load(transformed_df)
    return n_rows

Unique ID Construction

The unique_id field (e.g. "12_301") is constructed from item_id and site_id:

df['unique_id'] = df['item_id'].astype(str) + '_' + df['site_id'].astype(str)

This convention is still constructed at code edges for legacy paths, but no table stores a unique_id column — every ClickHouse PIPE_* table (and the migrated PG tables) is keyed by item_id + site_id (unique_id was dropped in the v6 migration).

Frequency Aggregation

FREQ_MAP = {
    "D": "D",    # daily
    "W": "W",    # weekly (default)
    "M": "MS",   # month start
    "Q": "QS",   # quarter start
    "Y": "YS",   # year start
}
AGG_MAP = {"sum": "sum", "mean": "mean", "median": "median"}

For weekly data, the ETL groups by (unique_id, week_start_date) and sums quantities. The aggregation.function config key selects the aggregation method.

Series Coverage Safety Gate

Before loading, the ETL compares the incoming series list against the existing series in the DB:

class ETLSeriesLossError(RuntimeError):
    """
    Raised when the incoming dataset would permanently remove series
    that have user-generated data (locked methods, adjustments, overrides).
    """
    def __init__(self, message, report):
        super().__init__(message)
        self.report = report  # {existing, new, dropped, loss_pct, at_risk}

If more than a configurable threshold (default 10%) of user-annotated series would be dropped, the ETL raises ETLSeriesLossError and aborts. This prevents accidental data loss from upstream schema changes. Pass --force to bypass this gate.

Outlier Detection — embedded in the Forecast step

Purpose

Detects statistical outliers in each demand series and writes corrections to scenario.demand_corrected. Outlier detection is not a standalone pipeline step — it runs automatically as the first action inside run_pipeline.py --only forecast, before the statistical/neural models are called. The corrected values are then used as the input to forecasting for that pipeline.

The implementation lives in files/outlier/detection.py but is invoked by run_pipeline.py:step_outlier_detection().

Three Detection Methods

The OutlierDetector class supports three methods, applied in combination:

1. IQR (Interquartile Range)

def _detect_iqr(self, values, multiplier=1.5):
    """
    Classic Tukey fence method.

    Fences: Q1 - k*IQR, Q3 + k*IQR
    Points outside the fences are flagged.

    k = iqr_multiplier (default 1.5, configurable per parameter set)
    Lower k → stricter, catches more outliers.
    """
    Q1 = np.percentile(values, 25)
    Q3 = np.percentile(values, 75)
    IQR = Q3 - Q1
    lower = Q1 - multiplier * IQR
    upper = Q3 + multiplier * IQR
    return (values < lower) | (values > upper)

2. Z-Score

def _detect_zscore(self, values, threshold=3.0):
    """
    Standardise values and flag those beyond the threshold.

    z = (x - mean) / std
    Outlier if |z| > threshold.

    Robust to the choice of location/scale when the series has no trend
    but can miss outliers in trended or seasonal series.
    """
    mean = np.nanmean(values)
    std  = np.nanstd(values)
    if std == 0:
        return np.zeros(len(values), dtype=bool)
    z = np.abs((values - mean) / std)
    return z > threshold

3. STL (Seasonal-Trend Decomposition)

def _detect_stl(self, series, period=52, threshold=3.0):
    """
    Decompose the series with STL, then z-score the residuals.

    STL separates trend + seasonal components, so the residuals represent
    only the irregular component. Outliers in residuals = genuine spikes
    not explained by trend or seasonality.

    Requires statsmodels. Falls back to IQR if unavailable.
    """
    from statsmodels.tsa.seasonal import STL
    stl = STL(series, period=period)
    result = stl.fit()
    residuals = result.resid
    return self._detect_zscore(residuals, threshold)

STL is the most accurate method for seasonal series but requires at least min_obs_for_stl observations (default 24). The final outlier mask is the union of all three methods.

Output

Corrections are written to scenario.demand_corrected (one row per corrected date per series per pipeline):

Column Description
item_id, site_id Series scope
date The corrected week
pipeline_id Isolation key — re-running the same pipeline overwrites only its own rows
original_value Raw demand value before correction
corrected_value Replacement value after applying the correction strategy
detection_method iqr, zscore, or stl
z_score Standardised distance from the series mean
lower_bound, upper_bound Detection fence values used for this data point

The forecasting step subsequently reads demand_corrected for the active pipeline and merges corrections back into the demand series before running models.

Characterization — files/characterization/characterization.py

Purpose

Analyses each time series and classifies it along several dimensions:

  • Seasonality (presence, periods, strength)
  • Trend (direction, strength)
  • Intermittency (ADI, CV², zero ratio)
  • Stationarity (ADF test)
  • Complexity (composite score)
  • Data sufficiency (for ML / deep learning)
  • Recommended forecasting methods

The SeriesCharacteristics Dataclass

@dataclass
class SeriesCharacteristics:
    unique_id: str
    # Basic stats
    n_observations: int = 0
    mean: float = np.nan
    std: float = np.nan
    # Seasonality
    has_seasonality: bool = False
    seasonal_periods: List[int] = field(default_factory=list)
    seasonal_strength: float = 0.0
    # Trend
    has_trend: bool = False
    trend_direction: str = 'none'  # 'up', 'down', 'none'
    trend_strength: float = 0.0
    # Intermittency
    is_intermittent: bool = False
    zero_ratio: float = 0.0
    adi: float = 0.0   # Average Demand Interval
    cov: float = 0.0   # Coefficient of Variation
    # Stationarity
    is_stationary: bool = False
    adf_pvalue: float = np.nan
    # Complexity
    complexity_score: float = 0.0
    complexity_level: str = 'low'  # 'low', 'medium', 'high'
    # Sufficiency
    sufficient_for_ml: bool = False
    sufficient_for_deep_learning: bool = False
    # Output
    recommended_methods: List[str] = field(default_factory=list)

analyze_all() — The Entry Point

def analyze_all(self, df, id_col='unique_id', date_col='ds', value_col='y', save=True):
    """
    Iterate over every unique_id and call analyze_single() for each.
    Results are accumulated into a DataFrame and saved to series_characteristics.
    """
    unique_ids = df[id_col].unique()
    results = []

    for uid in unique_ids:
        series_df = df[df[id_col] == uid].sort_values(date_col)
        try:
            chars = self.analyze_single(series_df, uid, date_col, value_col)
        except Exception:
            chars = SeriesCharacteristics(unique_id=uid)  # default fallback
        results.append(chars.to_dict())

    characteristics_df = pd.DataFrame(results)

    if save:
        bulk_insert(schema + '.series_characteristics', c_cols, c_rows)

    return characteristics_df

analyze_single() — Per-Series Analysis Pipeline

def analyze_single(self, series_df, unique_id, date_col, value_col):
    chars = SeriesCharacteristics(unique_id=unique_id)
    values = series_df[value_col].values.astype(float)
    dates  = pd.to_datetime(series_df[date_col])

    # Basic statistics (always computed)
    chars.n_observations  = len(values)
    chars.date_range_start = str(dates.min())
    chars.date_range_end   = str(dates.max())
    chars.mean = float(np.nanmean(values))
    chars.std  = float(np.nanstd(values))

    # Guard: skip deep analysis for very short series
    if chars.n_observations < 4:
        chars.recommended_methods = self._recommend_methods(chars)
        return chars

    self._detect_seasonality(values, chars, dates)
    self._detect_trend(values, chars)
    self._detect_intermittency(values, chars)
    self._detect_stationarity(values, chars)
    self._compute_complexity(values, chars)
    self._assess_sufficiency(chars)
    chars.recommended_methods = self._recommend_methods(chars)

    return chars

_detect_seasonality() — ACF-Based Seasonal Detection

The method uses the Autocorrelation Function (ACF) to detect seasonality. The core idea: if demand 52 weeks ago is strongly correlated with demand today, the series is seasonal with period 52.

def _detect_seasonality(self, values, chars, dates):
    # Gate: require at least 2 years of history to detect annual cycles
    date_span_days = (dates.max() - dates.min()).days
    if date_span_days < 2 * 365:
        chars.has_seasonality = False
        return

    # Reindex to uniform grid so ACF lags correspond to actual time periods.
    # Irregular gaps would shift lag indices and produce false negatives.
    # (Reindexing introduces NaN for missing dates; statsmodels handles
    # these with missing='conservative')
    acf_source = self._uniform_reindex(values, dates)

    test_periods = self.seasonality_cfg['test_periods']  # e.g. [4, 13, 26, 52]
    min_strength  = self.seasonality_cfg['min_strength']  # e.g. 0.3

    max_lag = max(test_periods) + 1
    acf_values = acf(acf_source, nlags=max_lag, missing='conservative')

    detected = []
    for period in test_periods:
        if period < len(acf_values) and acf_values[period] >= min_strength:
            detected.append(period)

    chars.has_seasonality   = len(detected) > 0
    chars.seasonal_periods  = detected
    chars.seasonal_strength = max((acf_values[p] for p in detected), default=0.0)

_detect_trend() — Mann-Kendall Test

def _mann_kendall_test(values, significance_level=0.05):
    """
    Non-parametric test for monotonic trend.

    S = sum of sign(x_j - x_i) for all j > i
    Var(S) is derived under the null hypothesis (no trend).
    z = S / sqrt(Var(S)) follows approximately N(0,1).

    Returns: (has_trend, direction, strength)
    """
    n = len(values)
    S = sum(
        np.sign(values[j] - values[i])
        for i in range(n - 1)
        for j in range(i + 1, n)
    )
    # Variance under null (accounting for ties)
    var_S = n * (n - 1) * (2 * n + 5) / 18

    if var_S == 0:
        return False, 'none', 0.0

    z = S / np.sqrt(var_S)
    p_value = 2 * (1 - scipy_stats.norm.cdf(abs(z)))
    has_trend = p_value < significance_level
    direction = 'up' if S > 0 else ('down' if S < 0 else 'none')
    strength  = abs(z)  # larger z = stronger trend
    return has_trend, direction, strength

_detect_intermittency() — ADI and CV²

def _detect_intermittency(self, values, chars):
    """
    Intermittency classification based on ADI and CV².

    ADI (Average Demand Interval): mean gap between non-zero periods.
    ADI > threshold → demand is too infrequent for standard methods.

    CV² (squared Coefficient of Variation): variance of non-zero quantities.
    High CV² with high ADI → lumpy demand (Syntetos-Boylan quadrant).
    """
    non_zero = values[values > 0]
    zeros     = values[values == 0]

    chars.zero_ratio = len(zeros) / len(values) if len(values) > 0 else 0.0

    if len(non_zero) == 0:
        chars.is_intermittent = True
        chars.adi = float('inf')
        return

    # ADI: number of periods per non-zero demand event
    chars.adi = len(values) / len(non_zero)

    # CV²: (std / mean)² of non-zero quantities
    mean_nz = np.mean(non_zero)
    std_nz  = np.std(non_zero)
    chars.cov = (std_nz / mean_nz) ** 2 if mean_nz > 0 else 0.0

    # Classification threshold from config
    adi_threshold = self.intermittency_cfg['adi_threshold']  # default 1.3
    chars.is_intermittent = chars.adi > adi_threshold

_recommend_methods() — Method Routing

def _recommend_methods(self, chars):
    """
    Route a series to appropriate forecasting methods based on its characteristics.

    Priority rules:
    1. Intermittent → always CrostonOptimized (optionally IMAPA)
    2. Multi-seasonal → MSTL
    3. Strong trend → AutoARIMA, AutoTheta
    4. Seasonal → AutoETS, MSTL
    5. All series → AutoARIMA, AutoETS as baseline
    6. Sufficient history → optionally add ML / neural models
    """
    methods = []

    if chars.is_intermittent:
        methods.append('CrostonOptimized')
        # Do not include MSTL for intermittent demand — it handles zeros poorly
    else:
        # Always include the two universal baselines
        methods.extend(['AutoARIMA', 'AutoETS'])

        if chars.has_seasonality:
            if len(chars.seasonal_periods) > 1:
                methods.append('MSTL')  # multiple seasonal periods
            else:
                methods.append('AutoETS')  # single seasonal period

        if chars.has_trend:
            if 'AutoARIMA' not in methods:
                methods.append('AutoARIMA')
            methods.append('AutoTheta')  # good for trend + seasonal

    # Neural models: only if the series has enough history
    if chars.sufficient_for_deep_learning:
        for m in ['NHITS', 'NBEATS', 'TFT', 'PatchTST', 'DeepAR']:
            if self._model_enabled(m):
                methods.append(m)
    elif chars.sufficient_for_ml:
        for m in ['LightGBM', 'XGBoost']:
            if self._model_enabled(m):
                methods.append(m)

    return list(dict.fromkeys(methods))  # deduplicate preserving order

Forecasting — files/forecasting/statistical_models.py

Purpose

Runs all enabled statistical forecasting models via the statsforecast library and produces point forecasts plus prediction intervals for each model and series.

StatisticalForecaster Class

class StatisticalForecaster:
    def __init__(self, config_path=None, config_override=None):
        # Load config (file → DB fallback → deep merge override)
        self.config = self._load_config(config_path, config_override)
        self.fc_config = self.config['forecasting']

        # Build the list of statsforecast model objects from config
        self.models = self._build_models()

_build_models() — Config-Driven Model Construction

def _build_models(self):
    """
    Inspect config['forecasting']['models'] and instantiate only the
    enabled models. Returns a list of statsforecast model objects.
    """
    from statsforecast.models import (
        AutoARIMA, AutoETS, AutoTheta, MSTL, CrostonOptimized
    )

    model_defs = self.fc_config.get('models', {})
    models = []

    for name, cfg in model_defs.items():
        if not cfg.get('enabled', False):
            continue
        season_length = cfg.get('season_length', 52)

        if name == 'AutoARIMA':
            models.append(AutoARIMA(season_length=season_length))
        elif name == 'AutoETS':
            models.append(AutoETS(season_length=season_length))
        elif name == 'MSTL':
            # MSTL accepts a list of seasonal periods
            models.append(MSTL(season_length=season_length if isinstance(season_length, list)
                               else [season_length]))
        elif name == 'CrostonOptimized':
            models.append(CrostonOptimized())
        # ... other models

    return models

forecast_multiple_series() — Batch Forecasting

def forecast_multiple_series(self, df, characteristics_df, overrides_map=None, show_progress=True):
    """
    Run all enabled models on all series in df.

    Uses StatsForecast's built-in parallelism (n_jobs=-1) for the
    statistical models. The characteristics_df is used to filter
    which models are appropriate for each series.
    """
    from statsforecast import StatsForecast

    # statsforecast expects columns: unique_id, ds, y
    sf = StatsForecast(
        models=self.models,
        freq=self.fc_config.get('frequency', 'W'),
        n_jobs=-1,    # use all available cores
        fallback_model=AutoETS(),  # if a model fails, use AutoETS instead
    )

    horizon = self.fc_config.get('horizon', 24)

    # Produce point forecasts + prediction intervals (80% and 90%)
    forecasts_df = sf.forecast(
        df=df,
        h=horizon,
        level=[80, 90],  # prediction interval levels
    )

    return forecasts_df

Per-SKU Overrides

When overrides_map is provided (from series_hyperparameters_overrides), the forecaster rebuilds the model configuration for those specific series with the override parameters applied:

# overrides_map: {unique_id: {method: {param: value}}}
for uid, method_overrides in overrides_map.items():
    for method_name, params in method_overrides.items():
        # Rebuild the specific model with overridden params
        # and re-forecast just that series

Backtesting / Evaluation — files/evaluation/metrics.py

Purpose

Implements rolling-window backtesting: for each series, the model is fit on a truncated history and evaluated on the held-out window. This produces unbiased estimates of forecast accuracy.

Rolling Window Logic

History: |-------|-------|-------|-------|-------|
Window 1:       |fit  |eval  |
Window 2:              |fit  |eval  |
Window 3:                     |fit  |eval  |

For each window (origin date): 1. Truncate the series at the origin. 2. Run the forecaster on the truncated series. 3. Compare the forecast against the actual held-out values. 4. Record metrics for this origin.

Final metrics per series × method = average across all origins.

ForecastEvaluator.evaluate()

def evaluate(self, df, forecasts_df, n_windows=3, window_size=12):
    """
    Rolling-window backtesting.

    Returns a DataFrame with columns:
        unique_id, method, forecast_origin, mae, rmse, mase, bias, coverage_90
    """
    origins = self._compute_origins(df, n_windows, window_size)
    all_metrics = []

    for origin_date in origins:
        # Truncate actuals to this origin
        train = df[df['ds'] <= origin_date]
        test  = df[(df['ds'] > origin_date) & 
                   (df['ds'] <= origin_date + window_size_periods)]

        # Forecast from this origin
        fc = self.forecaster.forecast_from_origin(train, horizon=window_size)

        # Compute metrics for each method
        for method in fc.columns:
            metrics = self._compute_metrics(
                actual=test['y'].values,
                predicted=fc[method].values
            )
            metrics.update({'unique_id': uid, 'method': method,
                           'forecast_origin': origin_date})
            all_metrics.append(metrics)

    return pd.DataFrame(all_metrics)

Metric Formulas

def _compute_metrics(self, actual, predicted):
    n = len(actual)
    errors = actual - predicted

    mae  = np.mean(np.abs(errors))
    rmse = np.sqrt(np.mean(errors ** 2))
    bias = np.mean(errors) / (np.mean(actual) + 1e-8)  # relative bias

    # MASE: normalise by naive forecast error (lag-1)
    naive_errors = np.abs(np.diff(actual))
    mase = mae / (np.mean(naive_errors) + 1e-8) if len(naive_errors) > 0 else np.nan

    # Coverage: what fraction of actuals fall within the 90% prediction interval?
    lower = predicted_lower_90
    upper = predicted_upper_90
    coverage_90 = np.mean((actual >= lower) & (actual <= upper))

    return {'mae': mae, 'rmse': rmse, 'mase': mase, 'bias': bias, 'coverage_90': coverage_90}

Best Method Selection — files/selection/best_method.py

Purpose

After backtesting, ranks all forecasting methods for each series using a weighted composite score. The method with the lowest composite score becomes the "best method" for that series.

_rank_methods_for_series()

def _rank_methods_for_series(self, unique_id, group):
    """
    Rank methods for one series. Steps:

    1. Average metrics across all forecast origins (windows).
    2. Min-max normalise each metric to [0, 1] within this series.
    3. Transform metrics so that lower = better:
       - bias → |bias| (both positive and negative are penalised)
       - coverage_90 → |coverage_90 - 0.9| (deviation from target)
    4. Compute weighted composite score.
    5. Rank methods ascending (lower = better).
    """
    # Step 1: average across origins
    avg = group.groupby('method')[list(self.weights.keys())].mean()

    # Step 3: transform bias and coverage
    if 'bias' in avg.columns:
        avg['bias'] = avg['bias'].abs()
    if 'coverage_90' in avg.columns:
        avg['coverage_90'] = (avg['coverage_90'] - 0.9).abs()

    # Step 2: min-max normalise each metric within this series
    for col in avg.columns:
        col_min = avg[col].min()
        col_max = avg[col].max()
        if col_max > col_min:
            avg[col] = (avg[col] - col_min) / (col_max - col_min)
        else:
            avg[col] = 0.0  # all methods tied on this metric

    # Step 4: weighted composite
    avg['composite_score'] = sum(
        avg[metric] * weight
        for metric, weight in self.weights.items()
        if metric in avg.columns
    )

    # Step 5: rank
    ranked = avg.sort_values('composite_score')
    best   = ranked.index[0]
    runner = ranked.index[1] if len(ranked) > 1 else None

    return {
        'unique_id':      unique_id,
        'best_method':    best,
        'best_score':     ranked.iloc[0]['composite_score'],
        'runner_up_method': runner,
        'runner_up_score':  ranked.iloc[1]['composite_score'] if runner else None,
        'all_rankings':   ranked['composite_score'].to_dict(),
    }

Default Weights

_DEFAULT_WEIGHTS = {
    "mae":         0.40,   # primary accuracy metric
    "rmse":        0.20,   # penalises large errors more
    "bias":        0.15,   # systematic over/under-forecasting
    "coverage_90": 0.15,   # prediction interval calibration
    "mase":        0.10,   # scale-independent accuracy
}

All weights sum to 1.0. These are overridable via the best_method parameter set in the DB.

Distribution Fitting — files/distribution/fitting.py

Purpose

Fits parametric probability distributions to each series' demand data. The fitted distribution is used by the MEIO optimizer to compute safety stock at the requested service level.

FittedDistribution Dataclass

@dataclass
class FittedDistribution:
    unique_id: str
    method: str             # which forecast method's residuals were used
    forecast_horizon: int
    distribution_type: str  # 'normal', 'gamma', 'negative_binomial',
                            # 'lognormal', 'weibull', 'poisson'
    mean: float
    std: float
    params: Dict            # distribution-specific MLE parameters
    ks_statistic: float     # Kolmogorov-Smirnov test statistic
    ks_pvalue: float        # KS p-value (higher = better fit)
    service_level_quantiles: Dict[float, float]  # {0.9: ..., 0.95: ..., 0.99: ...}

The Fitting Loop

def fit_distribution(self, values, unique_id, preferred=None):
    """
    Try each candidate distribution in order of preference.
    Return the first one that passes the KS goodness-of-fit test.

    Candidates: lognormal → gamma → normal → weibull → negative_binomial
    """
    preferred = preferred or self.config.get('preferred_distributions',
        ['lognormal', 'gamma', 'normal', 'weibull'])

    for dist_name in preferred:
        try:
            params, ks_stat, ks_pvalue = self._fit_and_test(values, dist_name)

            if ks_pvalue >= self.ks_pvalue_threshold:  # good fit
                return FittedDistribution(
                    unique_id=unique_id,
                    distribution_type=dist_name,
                    params=params,
                    ks_statistic=ks_stat,
                    ks_pvalue=ks_pvalue,
                    mean=np.mean(values),
                    std=np.std(values),
                    service_level_quantiles=self._compute_quantiles(dist_name, params),
                )
        except Exception:
            continue  # distribution failed to fit, try next

    # Fallback: use normal distribution regardless of fit quality
    return self._fit_normal_fallback(values, unique_id)

MLE Fitting and KS Test

def _fit_and_test(self, values, dist_name):
    """
    Fit by Maximum Likelihood Estimation using scipy.stats,
    then test goodness of fit with Kolmogorov-Smirnov.
    """
    dist_map = {
        'normal':   scipy.stats.norm,
        'gamma':    scipy.stats.gamma,
        'lognormal': scipy.stats.lognorm,
        'weibull':  scipy.stats.weibull_min,
    }
    dist = dist_map[dist_name]

    # MLE: scipy fits all parameters simultaneously
    params = dist.fit(values)  # returns (shape..., loc, scale) tuple

    # KS test: compare empirical CDF with fitted CDF
    ks_stat, ks_pvalue = scipy.stats.kstest(values, dist.cdf, args=params)

    return params, ks_stat, ks_pvalue

get_quantile() — Service Level Quantiles

def get_quantile(self, q):
    """Return the demand level at quantile q from the fitted distribution."""
    if self.distribution_type == 'lognormal':
        s, loc, scale = self.params['s'], self.params.get('loc', 0), self.params['scale']
        return scipy.stats.lognorm.ppf(q, s=s, loc=loc, scale=scale)
    elif self.distribution_type == 'gamma':
        a, loc, scale = self.params['shape'], self.params.get('loc', 0), self.params['scale']
        return scipy.stats.gamma.ppf(q, a=a, loc=loc, scale=scale)
    # ... other distributions

The service_level_quantiles dict (e.g. {0.90: 62.1, 0.95: 71.4, 0.99: 89.3}) is passed directly to the MEIO Rust optimizer as the demand uncertainty input.

MEIO Runner — files/meio_runner.py

Purpose

Orchestrates the MEIO (Multi-Echelon Inventory Optimisation) run. Loads all SKU data from PostgreSQL once, applies per-scenario parameter overrides, then calls the Rust optimizer extension.

Architecture Decision: Python → Rust

The heavy optimisation work is done by a compiled Rust extension (meio_optimizer.pyd on Windows). This approach gives:

  • Python: data loading, preprocessing, result persistence — benefits from Python's ecosystem.
  • Rust: the optimisation algorithm — inner loop runs at native speed with Rayon parallelism.

The boundary is a simple JSON interface:

# Python serialises the batch to JSON
rust_batches = [
    {
        "skus_json":          json.dumps(skus),
        "config_json":        json.dumps(config),
        "group_targets_json": json.dumps(targets),
    }
]

# Call Rust — returns JSON string of results
raw_results = json.loads(
    meio_optimizer.run_optimization_batch(json.dumps(rust_batches))
)

run_scenarios() — Main Entry Point

def run_scenarios(scenario_ids=None, base_only=False, config_path=None, pipeline_id=None):
    """
    1. Load scenarios from DB
    2. Split into availability (Erlang-B) vs standard (Rust) scenarios
    3. Load base SKU data once
    4. For each scenario: apply overrides → build batch
    5. Submit all batches to the Rust optimizer in one call
    6. Parse results → save to meio_results + meio_group_results
    """

Key Load Functions

def _load_sku_records(schema, pipeline_id=None):
    """
    Load all (item_id, site_id) combinations with:
    - demand statistics (mean, std, CV from fitted distributions)
    - unit cost (item_site.unit_cost → item.unit_cost fallback)
    - lead time (route.lead_time + transit_time + pick_pack_time)
    - EOQ: from pipeline-linked inventory_scenario K-Curve results (preferred)
            or item.eoq (static fallback) — see EOQ sourcing note below
    - j_target_groups (segment membership for fill-rate targeting)
    Returns (skus: list[dict], unit_cost_map: dict)
    """

EOQ Sourcing — K-Curve Results Take Priority

When a pipeline_id is supplied, _load_sku_records joins against the K-Curve results stored in inventory_scenario.results (JSONB) via the scenario_pipeline link table:

-- Simplified CTE added when pipeline_id is known
eoq_scenario AS (
    SELECT
        (elem->>'unique_id')             AS unique_id,
        (elem->>'eoq')::double precision AS eoq_qty
    FROM scenario.inventory_scenario inv
    JOIN scenario.scenario_pipeline sp
        ON sp.inventory_scenario_id = inv.scenario_id
       AND sp.pipeline_id = %(pipeline_id)s
    CROSS JOIN LATERAL jsonb_array_elements(inv.results->'items') AS elem
    WHERE inv.results ? 'items'
      AND (elem->>'eoq') IS NOT NULL
      AND (elem->>'eoq')::double precision > 0
)
-- Final EOQ expression:
COALESCE(eoq_scenario.eoq_qty, NULLIF(item.eoq, 0.0), 0.0)

Priority order: 1. K-Curve portfolio-optimised EOQ from inventory_scenario.results 2. item.eoq — static field seeded from the simple \(\sqrt{2DS/H}\) formula 3. 0.0 — when no ordering cost is configured

This ensures supply planning uses the jointly-optimised order quantities rather than per-item approximations. See K-Curve — Per-item EOQ output for the storage format.

def _load_sku_records(schema, pipeline_id=None):

def _load_meio_config(schema):
    """
    Load global MEIO parameters from scenario.parameters.
    Also loads meio_group_results targets.
    Returns (config: dict, group_targets: list[dict])
    """

def _load_scenario_segment_params(schema, scenario_id):
    """
    Load per-segment parameter overrides from meio_scenario_segment.
    These override fill_rate_target, max_budget, return_rate, etc.
    for specific item groups within this scenario.
    """

Per-Scenario Override Pipeline

for scenario in rust_scenarios:
    sid       = scenario['scenario_id']
    overrides = scenario.get('param_overrides') or {}
    seg_params = _load_scenario_segment_params(schema, sid)

    # Apply overrides to the base SKU data
    patched_skus    = _apply_sku_overrides(base_skus, overrides, seg_params=seg_params)
    patched_config  = _apply_config_overrides(base_config, overrides, seg_params=seg_params)
    patched_targets = _apply_target_overrides(base_group_targets, overrides, seg_params=seg_params)

    batches.append({
        "skus_json":          json.dumps(patched_skus),
        "config_json":        json.dumps(patched_config),
        "group_targets_json": json.dumps(patched_targets),
    })

Availability Scenarios (Erlang-B)

When a scenario is flagged as an availability scenario (asset-driven demand for MRO):

avail_scenarios = [s for s in scenarios if _is_availability_scenario(s)]
rust_scenarios  = [s for s in scenarios if not _is_availability_scenario(s)]

for s in avail_scenarios:
    # Erlang-B pool sizing (Python path, no Rust)
    sku_results, group_results = _run_causal_asset_scenario(schema, s, pipeline_id)

Erlang-B calculates the minimum pool size (buffer stock) needed to achieve a given availability service level, accounting for the failure rate distribution and repair cycle time.

Dask Orchestrator — files/utils/orchestrator.py

Purpose

Parallelizes the forecasting pipeline across multiple CPU cores using Dask distributed. The orchestrator splits the series into chunks and submits each chunk as a Dask task.

The Top-Level run() Function

def run(config_path, pipeline_id=None, scenario_id=1, segment_id=None):
    """
    Full pipeline run:
    1. Load demand actuals from DB
    2. Load characterization results
    3. Group series by parameter set (ParameterResolver)
    4. For each group: create Dask tasks, submit to cluster
    5. Gather results, persist to DB
    6. Run best method selection
    7. Run distribution fitting
    """

Forecast incremental skip — PIPE_forecast_fingerprints

A series is re-forecast only when its inputs changed. _get_changed_series() keeps a series when ANY of:

  1. No row in PIPE_forecast_point_values for the current archive_date (= planning_monday()) — this week not run yet.
  2. No data_hash in the demand-only / tenant-global PIPE_series_hashes (keyed by (item_id, site_id), no pipeline_id) — new series / ETL not run.
  3. The combined input fingerprint differs from the stored value in the pipeline-scoped PIPE_forecast_fingerprints (keyed by (pipeline_id, scenario_id, item_id, site_id), ReplacingMergeTree(computed_at), read with FINAL).

The combined fingerprint is MD5(data_hash | param_hash). param_hash (_compute_param_hash) folds in every forecast parameter surface: forecast parameter set, hyperparameter set, backtesting/method-selection set, and scenario-level overrides (confidence_levels, demand_horizon_limits, method_selection, frequency, horizon). It is memoised per-uid and cross-series (by assigned set IDs when there are no per-series overrides).

After a successful forecast run, _mark_series_forecast() writes the combined fingerprint to PIPE_forecast_fingerprints (it no longer touches PIPE_series_hashes.forecast_hash). Because the skip state is keyed by pipeline_id, multiple pipelines on one tenant cannot clobber each other's state; PIPE_series_hashes stays demand-only / tenant-global.

Last-resort Naive fallback

Inside StatisticalForecaster.forecast_multiple_series, after eligibility filtering, if a series has no eligible statistical method for its observation count (e.g. a 2-observation series routed to AutoETS/AutoARIMA/ AutoTheta/MSTL, all of which need ≥ 10), it is reassigned to Naive (needs only ≥ 1 observation) rather than aborting the whole pipeline with "0 forecasts produced". This keeps very short series forecastable.

Dask Task Structure

# The module-level function that Dask workers call
def _dask_forecast_batch(config_path, batch_df, batch_chars, config_override=None, overrides_map=None):
    """
    Standalone function (NOT a bound method) so it can be pickled and sent
    to a Dask worker process. Each worker reconstructs the forecasters from
    the config path.
    """
    all_forecasts = []

    # Statistical forecasting
    stat_forecaster = StatisticalForecaster(config_path, config_override=config_override)
    stat_forecasts  = stat_forecaster.forecast_multiple_series(
        df=batch_df,
        characteristics_df=batch_chars,
        overrides_map=overrides_map,
    )
    all_forecasts.append(stat_forecasts)

    # ML forecasting (if sufficient observations)
    if eligible_for_ml:
        ml_forecaster = MLForecaster(config_path, config_override=config_override)
        ml_forecasts  = ml_forecaster.forecast_multiple_series(batch_df, batch_chars)
        all_forecasts.append(ml_forecasts)

    return pd.concat(all_forecasts, ignore_index=True)

Parameter-Aware Batching

The orchestrator uses ParameterResolver.group_series_by_param_set() to group series so that all series using the same parameter set run together — this avoids re-instantiating the forecaster for each series:

resolver = ParameterResolver(config_path)
groups   = resolver.group_series_by_param_set(all_unique_ids, 'forecasting')

tasks = []
for param_id, uids in groups.items():
    override = resolver.build_group_config_override(param_id, 'forecasting')

    # Split this group into chunks of chunk_size
    for chunk_uids in chunked(uids, chunk_size):
        chunk_df    = df[df['unique_id'].isin(chunk_uids)]
        chunk_chars = chars_df[chars_df['unique_id'].isin(chunk_uids)]

        task = client.submit(
            _dask_forecast_batch,
            config_path,
            chunk_df,
            chunk_chars,
            override,
            overrides_map=sku_overrides_for_chunk,
        )
        tasks.append(task)

# Gather all results
results = [task.result() for task in as_completed(tasks)]

Dask Cluster Configuration

The orchestrator starts a dask.distributed.LocalCluster:

cluster = LocalCluster(
    n_workers         = config.get('dask', {}).get('n_workers', 4),
    threads_per_worker = config.get('dask', {}).get('threads_per_worker', 1),
    processes         = True,  # use separate processes for Python GIL avoidance
)
client = Client(cluster)

threads_per_worker=1 is important: statsforecast already uses its own internal parallelism, so giving each Dask worker one thread prevents over-subscription.

Process Logger — files/utils/process_logger.py

Every pipeline step creates a ProcessLogger to track its progress in scenario.process_log:

from utils.process_logger import ProcessLogger, ListHandler

_pl = ProcessLogger(run_id="forecast-2026-04-12", schema="scenario", pipeline_id=42)
step_id = _pl.start_step("forecast")

# ... do the work ...

_pl.end_step(step_id, "success", rows=590, log_tail=_list_handler.get_tail(200))

ListHandler is a Python logging handler that buffers log lines in memory so the API can stream them to the frontend via GET /api/pipeline/jobs/{job_id}/stream (Server-Sent Events).