""" Lifecycle-Aware Forecast Engine Generates 90-day per-product daily sales forecasts using analogous lifecycle curves learned from historical brand/category launch patterns. Usage: python forecast_engine.py python forecast_engine.py --backfill 30 Environment variables (from .env): DB_HOST, DB_USER, DB_PASSWORD, DB_NAME, DB_PORT (default 5432) """ import os import sys import json import time import logging from datetime import datetime, date, timedelta import numpy as np import pandas as pd import psycopg2 import psycopg2.extras import psycopg2.extensions from scipy.optimize import curve_fit from statsmodels.tsa.holtwinters import SimpleExpSmoothing, Holt # Register numpy type adapters so psycopg2 can serialize them to SQL psycopg2.extensions.register_adapter(np.float64, lambda x: psycopg2.extensions.AsIs(float(x))) psycopg2.extensions.register_adapter(np.float32, lambda x: psycopg2.extensions.AsIs(float(x))) psycopg2.extensions.register_adapter(np.int64, lambda x: psycopg2.extensions.AsIs(int(x))) psycopg2.extensions.register_adapter(np.int32, lambda x: psycopg2.extensions.AsIs(int(x))) # --------------------------------------------------------------------------- # Config # --------------------------------------------------------------------------- FORECAST_HORIZON_DAYS = 90 CURVE_HISTORY_DAYS = 730 # 2 years of launches to build reference curves CURVE_WINDOW_WEEKS = 13 # Track decay for 13 weeks (91 days) MIN_PRODUCTS_FOR_CURVE = 5 # Minimum launches to fit a brand curve MIN_PRODUCTS_FOR_BRAND_CAT = 10 # Minimum for brand x category curve MATURE_VELOCITY_THRESHOLD = 0.1 # units/day to qualify as "mature" vs "dormant" MATURE_AGE_DAYS = 60 # days since first_received to be considered mature LAUNCH_AGE_DAYS = 14 # days in "launch" phase DECAY_AGE_DAYS = 60 # days in "active decay" phase EXP_SMOOTHING_WINDOW = 60 # days of history for mature product smoothing BATCH_SIZE = 1000 # rows per INSERT batch DOW_LOOKBACK_DAYS = 90 # days of order history for day-of-week indices MIN_R_SQUARED = 0.1 # curves below this are unreliable (fall back to velocity) SEASONAL_LOOKBACK_DAYS = 365 # 12 months of order history for monthly seasonal indices MIN_PREORDER_DAYS = 3 # minimum pre-order accumulation days for reliable scaling MAX_SMOOTHING_MULTIPLIER = 10 # cap exp smoothing forecast at Nx observed velocity logging.basicConfig( level=logging.INFO, format='%(asctime)s [%(levelname)s] %(message)s', datefmt='%H:%M:%S' ) log = logging.getLogger('forecast') # --------------------------------------------------------------------------- # Database helpers # --------------------------------------------------------------------------- def get_connection(): """Create a PostgreSQL connection from environment variables.""" return psycopg2.connect( host=os.environ.get('DB_HOST', 'localhost'), user=os.environ.get('DB_USER', 'inventory_user'), password=os.environ.get('DB_PASSWORD', ''), dbname=os.environ.get('DB_NAME', 'inventory_db'), port=int(os.environ.get('DB_PORT', 5432)), ) def execute_query(conn, sql, params=None): """Execute a query and return a DataFrame.""" import warnings with warnings.catch_warnings(): warnings.filterwarnings("ignore", message=".*pandas only supports SQLAlchemy.*") return pd.read_sql_query(sql, conn, params=params) def cleanup_stale_runs(conn): """Mark any runs stuck in 'running' status as failed (e.g. from prior crashes).""" with conn.cursor() as cur: cur.execute(""" UPDATE forecast_runs SET status = 'failed', finished_at = NOW(), error_message = 'Stale run cleaned up on engine restart' WHERE status = 'running' AND started_at < NOW() - INTERVAL '1 hour' """) cleaned = cur.rowcount conn.commit() if cleaned > 0: log.info(f"Cleaned up {cleaned} stale 'running' forecast run(s)") # --------------------------------------------------------------------------- # Decay curve model: sales(t) = A * exp(-λt) + C # --------------------------------------------------------------------------- def decay_model(t, amplitude, decay_rate, baseline): """Parametric exponential decay with baseline.""" return amplitude * np.exp(-decay_rate * t) + baseline def fit_decay_curve(weekly_medians): """ Fit the decay model to median weekly sales data. Args: weekly_medians: array of median sales per week (index = week number) Returns: (amplitude, decay_rate, baseline, r_squared) or None if fit fails """ weeks = np.arange(len(weekly_medians), dtype=float) y = np.array(weekly_medians, dtype=float) # Skip if all zeros or too few points if len(y) < 3 or np.max(y) == 0: return None # Initial guesses a0 = float(np.max(y)) c0 = float(np.min(y[len(y)//2:])) # baseline from second half lam0 = 0.3 # moderate decay try: popt, _ = curve_fit( decay_model, weeks, y, p0=[a0, lam0, c0], bounds=([0, 0.01, 0], [a0 * 5, 5.0, a0]), maxfev=5000, ) amplitude, decay_rate, baseline = popt # R-squared y_pred = decay_model(weeks, *popt) ss_res = np.sum((y - y_pred) ** 2) ss_tot = np.sum((y - np.mean(y)) ** 2) r_sq = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0.0 return float(amplitude), float(decay_rate), float(baseline), float(r_sq) except (RuntimeError, ValueError) as e: log.debug(f"Curve fit failed: {e}") return None # --------------------------------------------------------------------------- # Day-of-week indices # --------------------------------------------------------------------------- def compute_dow_indices(conn): """ Compute day-of-week revenue indices from recent order history. Returns a dict mapping ISO weekday (1=Mon ... 7=Sun) to a multiplier normalized so they sum to 7.0 (average = 1.0). This means applying them preserves the weekly total while reshaping the daily distribution. """ sql = """ SELECT EXTRACT(ISODOW FROM o.date)::int AS dow, SUM(o.price * o.quantity) AS revenue FROM orders o WHERE o.canceled IS DISTINCT FROM TRUE AND o.date >= CURRENT_DATE - INTERVAL '1 day' * %s GROUP BY 1 ORDER BY 1 """ df = execute_query(conn, sql, [DOW_LOOKBACK_DAYS]) if df.empty or len(df) < 7: log.warning("Insufficient order data for DOW indices, using flat distribution") return {d: 1.0 for d in range(1, 8)} total = df['revenue'].sum() avg = total / 7.0 indices = {} for _, row in df.iterrows(): dow = int(row['dow']) idx = float(row['revenue']) / avg if avg > 0 else 1.0 indices[dow] = round(idx, 4) # Fill any missing days for d in range(1, 8): if d not in indices: indices[d] = 1.0 log.info(f"DOW indices: Mon={indices[1]:.3f} Tue={indices[2]:.3f} Wed={indices[3]:.3f} " f"Thu={indices[4]:.3f} Fri={indices[5]:.3f} Sat={indices[6]:.3f} Sun={indices[7]:.3f}") return indices # --------------------------------------------------------------------------- # Monthly seasonal indices # --------------------------------------------------------------------------- def compute_monthly_seasonal_indices(conn): """ Compute monthly seasonal indices from recent order revenue. Returns a dict mapping month number (1-12) to a multiplier normalized so they average 1.0. Months with above-average revenue get >1, below get <1. """ sql = """ SELECT EXTRACT(MONTH FROM o.date)::int AS month, SUM(o.price * o.quantity) AS revenue FROM orders o WHERE o.canceled IS DISTINCT FROM TRUE AND o.date >= CURRENT_DATE - INTERVAL '1 day' * %s GROUP BY 1 ORDER BY 1 """ df = execute_query(conn, sql, [SEASONAL_LOOKBACK_DAYS]) if df.empty or len(df) < 6: log.warning("Insufficient data for seasonal indices, using flat distribution") return {m: 1.0 for m in range(1, 13)} total = df['revenue'].sum() n_months = len(df) avg = total / n_months indices = {} for _, row in df.iterrows(): month = int(row['month']) idx = float(row['revenue']) / avg if avg > 0 else 1.0 indices[month] = round(idx, 4) # Fill any missing months with 1.0 for m in range(1, 13): if m not in indices: indices[m] = 1.0 present = [f"{m}={indices[m]:.3f}" for m in range(1, 13)] log.info(f"Monthly seasonal indices: {', '.join(present)}") return indices # --------------------------------------------------------------------------- # Phase 1: Build brand-category reference curves # --------------------------------------------------------------------------- DEAL_CATEGORIES = frozenset([ 'Deals', 'Black Friday', 'Week 1', 'Week 2', 'Week 3', '28 Off', '5 Dollar Deals', '10 Dollar Deals', 'Fall Sale', ]) def build_reference_curves(conn): """ Build decay curves for each brand (and brand x category at every hierarchy level). For category curves, we load each product's full set of category assignments (across all hierarchy levels), then fit brand×cat_id curves wherever we have enough products. This gives granular curves like "49 and Market × 12x12 Paper Pads" alongside coarser fallbacks like "49 and Market × Paper". Returns DataFrame of curves written to brand_lifecycle_curves. """ log.info("Building reference curves from historical launches...") # Get daily sales aligned by days-since-first-received for recent launches # (no category join here — we attach categories separately) sales_sql = """ WITH recent_launches AS ( SELECT pm.pid, p.brand FROM product_metrics pm JOIN products p ON p.pid = pm.pid WHERE p.visible = true AND p.brand IS NOT NULL AND pm.date_first_received >= NOW() - INTERVAL '1 day' * %s AND pm.date_first_received < NOW() - INTERVAL '14 days' ), daily_sales AS ( SELECT rl.pid, rl.brand, dps.snapshot_date, COALESCE(dps.units_sold, 0) AS units_sold, (dps.snapshot_date - pm.date_first_received::date) AS day_offset FROM recent_launches rl JOIN product_metrics pm ON pm.pid = rl.pid JOIN daily_product_snapshots dps ON dps.pid = rl.pid WHERE dps.snapshot_date >= pm.date_first_received::date AND dps.snapshot_date < pm.date_first_received::date + INTERVAL '1 week' * %s ) SELECT brand, pid, FLOOR(day_offset / 7)::int AS week_num, SUM(units_sold) AS weekly_sales FROM daily_sales WHERE day_offset >= 0 GROUP BY brand, pid, week_num ORDER BY brand, pid, week_num """ df = execute_query(conn, sales_sql, [CURVE_HISTORY_DAYS, CURVE_WINDOW_WEEKS]) if df.empty: log.warning("No launch data found for reference curves") return pd.DataFrame() log.info(f"Loaded {len(df)} weekly sales records from {df['pid'].nunique()} products across {df['brand'].nunique()} brands") # Load all category assignments for these products (every hierarchy level) launch_pids = df['pid'].unique().tolist() cat_sql = """ SELECT pc.pid, ch.cat_id, ch.name AS cat_name, ch.level AS cat_level FROM product_categories pc JOIN category_hierarchy ch ON ch.cat_id = pc.cat_id WHERE pc.pid = ANY(%s) AND ch.name NOT IN %s ORDER BY pc.pid, ch.level DESC """ cat_df = execute_query(conn, cat_sql, [launch_pids, tuple(DEAL_CATEGORIES)]) # Build pid -> list of (cat_id, cat_name, cat_level) pid_cats = {} for _, row in cat_df.iterrows(): pid = int(row['pid']) if pid not in pid_cats: pid_cats[pid] = [] pid_cats[pid].append((int(row['cat_id']), row['cat_name'], int(row['cat_level']))) # Also get pre-order stats per brand (median pre-order sales AND accumulation window). # Uses de-facto preorders: any product that had orders before date_first_received, # regardless of the preorder_count flag. This gives us 6000+ completed cycles vs ~19 # from the explicit flag alone. preorder_sql = """ WITH preorder_stats AS ( SELECT p.pid, p.brand, COALESCE((SELECT SUM(o.quantity) FROM orders o WHERE o.pid = p.pid AND o.canceled IS DISTINCT FROM TRUE AND o.date < pm.date_first_received), 0) AS preorder_units, GREATEST(EXTRACT(DAY FROM pm.date_first_received - MIN(o.date)), 1) AS preorder_days FROM products p JOIN product_metrics pm ON pm.pid = p.pid LEFT JOIN orders o ON o.pid = p.pid AND o.canceled IS DISTINCT FROM TRUE AND o.date < pm.date_first_received WHERE p.visible = true AND p.brand IS NOT NULL AND pm.date_first_received IS NOT NULL AND pm.date_first_received >= NOW() - INTERVAL '1 day' * %s GROUP BY p.pid, p.brand, pm.date_first_received ) SELECT brand, PERCENTILE_CONT(0.5) WITHIN GROUP (ORDER BY preorder_units) AS median_preorder_sales, PERCENTILE_CONT(0.5) WITHIN GROUP (ORDER BY preorder_days) AS median_preorder_days FROM preorder_stats WHERE preorder_units > 0 GROUP BY brand HAVING COUNT(*) >= 3 """ preorder_df = execute_query(conn, preorder_sql, [CURVE_HISTORY_DAYS]) preorder_map = dict(zip(preorder_df['brand'], preorder_df['median_preorder_sales'])) if not preorder_df.empty else {} preorder_days_map = dict(zip(preorder_df['brand'], preorder_df['median_preorder_days'])) if not preorder_df.empty else {} curves = [] def _fit_and_append(group_df, brand, cat_id=None, cat_name=None, cat_level=None): """Helper: fit a decay curve for a group and append to curves list.""" product_count = group_df['pid'].nunique() min_products = MIN_PRODUCTS_FOR_CURVE if cat_id is None else MIN_PRODUCTS_FOR_BRAND_CAT if product_count < min_products: return False weekly = group_df.groupby('week_num')['weekly_sales'].median() if len(weekly) < 4: return False full_weeks = weekly.reindex(range(CURVE_WINDOW_WEEKS), fill_value=0.0) weekly_arr = full_weeks.values[:CURVE_WINDOW_WEEKS] result = fit_decay_curve(weekly_arr) if result is None: return False amplitude, decay_rate, baseline, r_sq = result # Quality gate: only store curves above the reliability threshold if r_sq < MIN_R_SQUARED: return False first_week = group_df[group_df['week_num'] == 0].groupby('pid')['weekly_sales'].sum() median_fw = float(first_week.median()) if len(first_week) > 0 else 0.0 curves.append({ 'brand': brand, 'root_category': cat_name, # kept for readability; cat_id is the real key 'cat_id': cat_id, 'category_level': cat_level, 'amplitude': amplitude, 'decay_rate': decay_rate, 'baseline': baseline, 'r_squared': r_sq, 'sample_size': product_count, 'median_first_week_sales': median_fw, 'median_preorder_sales': preorder_map.get(brand), 'median_preorder_days': preorder_days_map.get(brand), }) return True # 1. Fit brand-level curves (aggregate across all categories) for brand, brand_df in df.groupby('brand'): _fit_and_append(brand_df, brand) # 2. Fit brand × category curves at every hierarchy level # Build a mapping of (brand, cat_id) -> list of pids brand_cat_pids = {} for pid, cats in pid_cats.items(): brand_rows = df[df['pid'] == pid] if brand_rows.empty: continue brand = brand_rows.iloc[0]['brand'] for cat_id, cat_name, cat_level in cats: key = (brand, cat_id) if key not in brand_cat_pids: brand_cat_pids[key] = {'cat_name': cat_name, 'cat_level': cat_level, 'pids': set()} brand_cat_pids[key]['pids'].add(pid) cat_curves_fitted = 0 for (brand, cat_id), info in brand_cat_pids.items(): group_df = df[(df['brand'] == brand) & (df['pid'].isin(info['pids']))] if _fit_and_append(group_df, brand, cat_id=cat_id, cat_name=info['cat_name'], cat_level=info['cat_level']): cat_curves_fitted += 1 curves_df = pd.DataFrame(curves) if curves_df.empty: log.warning("No curves could be fitted") return curves_df # Write to database with conn.cursor() as cur: cur.execute("TRUNCATE brand_lifecycle_curves") for _, row in curves_df.iterrows(): cur.execute(""" INSERT INTO brand_lifecycle_curves (brand, root_category, cat_id, category_level, amplitude, decay_rate, baseline, r_squared, sample_size, median_first_week_sales, median_preorder_sales, median_preorder_days, computed_at) VALUES (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, NOW()) """, ( row['brand'], None if pd.isna(row.get('root_category')) else row.get('root_category'), None if pd.isna(row.get('cat_id')) else int(row['cat_id']), None if pd.isna(row.get('category_level')) else int(row['category_level']), row['amplitude'], row['decay_rate'], row['baseline'], row['r_squared'], row['sample_size'], row['median_first_week_sales'], None if pd.isna(row.get('median_preorder_sales')) else row.get('median_preorder_sales'), None if pd.isna(row.get('median_preorder_days')) else row.get('median_preorder_days'), )) conn.commit() brand_only = curves_df[curves_df['cat_id'].isna()].shape[0] cat_total = curves_df[curves_df['cat_id'].notna()].shape[0] log.info(f"Wrote {len(curves_df)} reference curves ({cat_total} brand+category across all levels, {brand_only} brand-only)") return curves_df # --------------------------------------------------------------------------- # Phase 2: Classify products and generate forecasts # --------------------------------------------------------------------------- def load_products(conn): """ Load all visible products with their metrics for classification. Also loads each product's full category ancestry (all hierarchy levels), stored as a list of cat_ids ordered deepest-first for hierarchical curve lookup. """ sql = """ SELECT pm.pid, p.brand, pm.current_price, pm.current_stock, pm.sales_velocity_daily, pm.sales_30d, pm.date_first_received, pm.date_last_sold, p.preorder_count, COALESCE(p.baskets, 0) AS baskets, EXTRACT(DAY FROM NOW() - pm.date_first_received) AS age_days FROM product_metrics pm JOIN products p ON p.pid = pm.pid WHERE p.visible = true """ products = execute_query(conn, sql) # Load category assignments for all products (every hierarchy level, deepest first) cat_sql = """ SELECT pc.pid, ch.cat_id, ch.level AS cat_level FROM product_categories pc JOIN category_hierarchy ch ON ch.cat_id = pc.cat_id WHERE ch.name NOT IN %s ORDER BY pc.pid, ch.level DESC """ cat_df = execute_query(conn, cat_sql, [tuple(DEAL_CATEGORIES)]) # Build pid -> [cat_id, ...] ordered deepest-first pid_cat_ids = {} for _, row in cat_df.iterrows(): pid = int(row['pid']) if pid not in pid_cat_ids: pid_cat_ids[pid] = [] pid_cat_ids[pid].append(int(row['cat_id'])) # Attach category list to each product row as a Python object column products['cat_ids'] = products['pid'].apply(lambda p: pid_cat_ids.get(int(p), [])) log.info(f"Loaded {len(products)} products, " f"{sum(1 for c in products['cat_ids'] if len(c) > 0)}/{len(products)} with category data") return products def classify_phase(row): """Classify a product into its lifecycle phase.""" preorder = (row.get('preorder_count') or 0) > 0 age = row.get('age_days') velocity = row.get('sales_velocity_daily') or 0 first_received = row.get('date_first_received') # Pre-order: has preorder_count and either not received or very recently received if preorder and (first_received is None or (age is not None and age <= LAUNCH_AGE_DAYS)): return 'preorder' # No first_received date — can't determine lifecycle if first_received is None or age is None: if velocity > MATURE_VELOCITY_THRESHOLD: return 'mature' if velocity > 0: return 'slow_mover' return 'dormant' if age <= LAUNCH_AGE_DAYS: return 'launch' elif age <= DECAY_AGE_DAYS: if velocity > 0: return 'decay' return 'dormant' else: if velocity > MATURE_VELOCITY_THRESHOLD: return 'mature' if velocity > 0: return 'slow_mover' return 'dormant' def get_curve_for_product(product, curves_df): """ Look up the best matching reference curve for a product. Uses hierarchical category fallback: tries the product's deepest category first, walks up the hierarchy to coarser categories, then falls back to brand-only. This ensures e.g. "49 and Market × 12x12 Paper Pads" is preferred over "49 and Market × Paper Crafts" when available. Skips curves with R² below MIN_R_SQUARED (unreliable fits). Returns (amplitude, decay_rate, baseline, median_first_week, median_preorder, median_preorder_days) or None. """ brand = product.get('brand') if brand is None or curves_df.empty: return None # Filter to this brand's reliable curves once brand_curves = curves_df[ (curves_df['brand'] == brand) & (curves_df['r_squared'] >= MIN_R_SQUARED) ] if brand_curves.empty: return None def _extract(row): return ( float(row['amplitude']), float(row['decay_rate']), float(row['baseline']), float(row['median_first_week_sales'] or 1), float(row['median_preorder_sales']) if pd.notna(row.get('median_preorder_sales')) else None, float(row['median_preorder_days']) if pd.notna(row.get('median_preorder_days')) else None, ) # Try each category from deepest to shallowest cat_ids = product.get('cat_ids') or [] for cat_id in cat_ids: match = brand_curves[brand_curves['cat_id'] == cat_id] if not match.empty: return _extract(match.iloc[0]) # Fall back to brand-only curve (cat_id is NaN/None) brand_only = brand_curves[brand_curves['cat_id'].isna()] if brand_only.empty: return None return _extract(brand_only.iloc[0]) def forecast_from_curve(curve_params, scale_factor, age_days, horizon_days): """ Generate daily forecast from a scaled decay curve. The scale factor is applied only to the decay envelope, NOT the baseline. This prevents hot products from getting inflated long-tail forecasts. Formula: daily_value = (A/7) * exp(-λ * t_weeks) * scale + (C/7) Args: curve_params: (amplitude, decay_rate, baseline, ...) scale_factor: multiplier for the decay envelope age_days: current product age in days horizon_days: how many days to forecast Returns: array of daily forecast values """ amplitude, decay_rate, baseline = curve_params[:3] # The curve is in weekly units; convert to daily daily_amp = amplitude / 7.0 daily_baseline = baseline / 7.0 forecasts = [] for d in range(horizon_days): t_weeks = (age_days + d) / 7.0 daily_value = daily_amp * np.exp(-decay_rate * t_weeks) * scale_factor + daily_baseline forecasts.append(max(0.0, daily_value)) return np.array(forecasts) def forecast_preorder(curve_params, scale_factor, days_until_arrival, preorder_daily_rate, horizon_days): """ Piecewise pre-order forecast: a flat observed pre-order trickle until the product is expected to arrive, then the scaled launch curve from age 0. The launch curve was fit on POST-receipt order history, so running it from today (while the product is still weeks from arriving) front-loads full first-week launch volume that hasn't happened yet — the main driver of the ~2.15x preorder over-forecast. Instead we forecast the slow pre-order rate up to the arrival date, then start the curve's day 0 on that date. See FORECAST_FIX_PLAN F4. Args: curve_params: (amplitude, decay_rate, baseline, ...) weekly curve scale_factor: per-product multiplier for the post-arrival curve envelope days_until_arrival: calendar days from today until expected arrival preorder_daily_rate: observed pre-order units/day (trickle) horizon_days: forecast horizon length Returns: array of daily forecast values of length horizon_days """ amplitude, decay_rate, baseline = curve_params[:3] forecasts = np.zeros(horizon_days) # Clamp the arrival offset into the horizon dua = int(max(0, min(days_until_arrival, horizon_days))) # Pre-arrival segment: flat pre-order trickle, capped at the curve's scaled # week-0 daily value (a pre-order day shouldn't out-sell the launch peak). if dua > 0: week0_daily = (amplitude / 7.0) * scale_factor + (baseline / 7.0) pre_rate = preorder_daily_rate if week0_daily > 0: pre_rate = min(pre_rate, week0_daily) forecasts[:dua] = max(0.0, pre_rate) # Post-arrival segment: scaled launch curve, curve day 0 = arrival date. if dua < horizon_days: curve_part = forecast_from_curve(curve_params, scale_factor, 0, horizon_days - dua) forecasts[dua:] = curve_part return forecasts # --------------------------------------------------------------------------- # Batch data loading (eliminates N+1 per-product queries) # --------------------------------------------------------------------------- def batch_load_product_data(conn, products): """ Batch-load all per-product data needed for forecasting in a few queries instead of one query per product. Returns dict with keys: preorder_sales: {pid: units} — pre-order units (before first received) launch_sales: {pid: units} — first 14 days of sales decay_velocity: {pid: avg} — recent 30-day daily average mature_history: {pid: DataFrame} — daily sales history for SES """ data = { 'preorder_sales': {}, 'preorder_days': {}, 'preorder_arrival_days': {}, 'launch_sales': {}, 'decay_velocity': {}, 'mature_history': {}, 'dormant_rate': {}, } # Pre-order sales: orders placed BEFORE first received date # Also compute the number of days pre-orders accumulated over (for daily-rate normalization) preorder_pids = products[products['phase'] == 'preorder']['pid'].tolist() if preorder_pids: sql = """ SELECT o.pid, COALESCE(SUM(o.quantity), 0) AS preorder_units, GREATEST(EXTRACT(DAY FROM NOW() - MIN(o.date)), 1) AS preorder_days FROM orders o LEFT JOIN product_metrics pm ON pm.pid = o.pid WHERE o.pid = ANY(%s) AND o.canceled IS DISTINCT FROM TRUE AND (pm.date_first_received IS NULL OR o.date < pm.date_first_received) GROUP BY o.pid """ df = execute_query(conn, sql, [preorder_pids]) for _, row in df.iterrows(): data['preorder_sales'][int(row['pid'])] = float(row['preorder_units']) data['preorder_days'][int(row['pid'])] = float(row['preorder_days']) log.info(f"Batch loaded pre-order sales for {len(data['preorder_sales'])}/{len(preorder_pids)} preorder products") # Expected arrival per pre-order product, to time the launch curve. # Prefer the soonest FUTURE expected_date on an open PO; if the only open # PO has a past expected_date assume 7 days; if there's no open PO at all # assume 14 days. See FORECAST_FIX_PLAN F4. arrival_sql = """ SELECT pid, MIN(expected_date) FILTER ( WHERE expected_date IS NOT NULL AND expected_date >= CURRENT_DATE ) AS future_arrival FROM purchase_orders WHERE pid = ANY(%s) AND status IN ('created', 'ordered', 'electronically_sent', 'receiving_started') GROUP BY pid """ adf = execute_query(conn, arrival_sql, [preorder_pids]) today = date.today() for _, row in adf.iterrows(): pid = int(row['pid']) fa = row['future_arrival'] if pd.notna(fa): fa_date = pd.Timestamp(fa).date() data['preorder_arrival_days'][pid] = max(0, (fa_date - today).days) else: data['preorder_arrival_days'][pid] = 7 # open PO, expected_date already past no_po = 0 for pid in preorder_pids: if int(pid) not in data['preorder_arrival_days']: data['preorder_arrival_days'][int(pid)] = 14 # no open PO at all no_po += 1 log.info(f"Batch loaded preorder arrival for " f"{len(data['preorder_arrival_days']) - no_po}/{len(preorder_pids)} via open POs, " f"{no_po} defaulted to 14d") # Launch sales: first 14 days after first received launch_pids = products[products['phase'] == 'launch']['pid'].tolist() if launch_pids: sql = """ SELECT dps.pid, COALESCE(SUM(dps.units_sold), 0) AS total_sold FROM daily_product_snapshots dps JOIN product_metrics pm ON pm.pid = dps.pid WHERE dps.pid = ANY(%s) AND dps.snapshot_date >= pm.date_first_received::date AND dps.snapshot_date < pm.date_first_received::date + INTERVAL '14 days' GROUP BY dps.pid """ df = execute_query(conn, sql, [launch_pids]) for _, row in df.iterrows(): data['launch_sales'][int(row['pid'])] = float(row['total_sold']) log.info(f"Batch loaded launch sales for {len(data['launch_sales'])}/{len(launch_pids)} launch products") # Decay recent velocity: TRUE calendar-daily average over the last 30 days. # We divide the summed units by calendar days (clipped to the product's age), # NOT by the number of snapshot rows. Snapshots are sparse and mostly land on # sold-days, so AVG(units_sold) averages over sold-days only and inflated the # decay rate ~4x (measured 1.353 vs true 0.332 units/day). See FORECAST_FIX_PLAN F1. decay_pids = products[products['phase'] == 'decay']['pid'].tolist() if decay_pids: sql = """ SELECT dps.pid, SUM(COALESCE(dps.units_sold, 0))::float / GREATEST(LEAST(30, (CURRENT_DATE - pm.date_first_received::date)), 1) AS avg_daily FROM daily_product_snapshots dps JOIN product_metrics pm ON pm.pid = dps.pid WHERE dps.pid = ANY(%s) AND dps.snapshot_date >= CURRENT_DATE - INTERVAL '30 days' AND dps.snapshot_date >= pm.date_first_received::date GROUP BY dps.pid, pm.date_first_received """ df = execute_query(conn, sql, [decay_pids]) for _, row in df.iterrows(): data['decay_velocity'][int(row['pid'])] = float(row['avg_daily']) log.info(f"Batch loaded decay velocity for {len(data['decay_velocity'])}/{len(decay_pids)} decay products") # Mature daily history: full time series for exponential smoothing mature_pids = products[products['phase'] == 'mature']['pid'].tolist() if mature_pids: sql = """ SELECT dps.pid, dps.snapshot_date, COALESCE(dps.units_sold, 0) AS units_sold FROM daily_product_snapshots dps WHERE dps.pid = ANY(%s) AND dps.snapshot_date >= CURRENT_DATE - INTERVAL '1 day' * %s ORDER BY dps.pid, dps.snapshot_date """ df = execute_query(conn, sql, [mature_pids, EXP_SMOOTHING_WINDOW]) for pid, group in df.groupby('pid'): data['mature_history'][int(pid)] = group.copy() log.info(f"Batch loaded history for {len(data['mature_history'])}/{len(mature_pids)} mature products") # Dormant trailing order rate: dormant products forecast 0 by default, but # ~11K of them still sell (restocks, promos, long-tail) — ~11% of all demand # currently forecast as a hard zero. Load a trailing-180-day daily order rate # so the dormant branch can carry a small positive rate. See FORECAST_FIX_PLAN F5. dormant_pids = products[products['phase'] == 'dormant']['pid'].tolist() if dormant_pids: sql = """ SELECT o.pid, SUM(o.quantity) / 180.0 AS rate FROM orders o WHERE o.pid = ANY(%s) AND o.canceled IS DISTINCT FROM TRUE AND o.date >= CURRENT_DATE - INTERVAL '180 days' GROUP BY o.pid """ df = execute_query(conn, sql, [dormant_pids]) for _, row in df.iterrows(): data['dormant_rate'][int(row['pid'])] = float(row['rate']) log.info(f"Batch loaded dormant order rate for {len(data['dormant_rate'])}/{len(dormant_pids)} dormant products") return data # --------------------------------------------------------------------------- # Per-product scale factor computation # --------------------------------------------------------------------------- def compute_scale_factor(phase, product, curve_info, batch_data): """ Compute the per-product scale factor for the brand curve. The scale factor captures how much more/less this product sells compared to the brand average. It's applied to the decay envelope only (not baseline). """ if curve_info is None: return 1.0 pid = int(product['pid']) amplitude, decay_rate, baseline, median_fw, median_preorder, med_preorder_days = curve_info if phase == 'preorder': preorder_units = batch_data['preorder_sales'].get(pid, 0) preorder_days = batch_data['preorder_days'].get(pid, 1) baskets = product.get('baskets') or 0 # Too few days of accumulation → noisy signal, use brand average if preorder_days < MIN_PREORDER_DAYS and preorder_units > 0: scale = 1.0 return max(0.1, min(scale, 5.0)) # Use order units as primary signal; fall back to baskets if no orders demand_signal = preorder_units if preorder_units > 0 else baskets signal_days = preorder_days if preorder_units > 0 else max(preorder_days, 14) # Normalize to daily rate before comparing to brand median daily rate. # Use the brand's stored median pre-order window for the denominator # (not the current product's signal_days) to avoid systematic bias. demand_daily = demand_signal / max(signal_days, 1) if median_preorder and median_preorder > 0: brand_preorder_window = max(med_preorder_days or signal_days, 1) median_preorder_daily = median_preorder / brand_preorder_window scale = demand_daily / median_preorder_daily elif median_fw > 0 and demand_daily > 0: median_fw_daily = median_fw / 7.0 scale = demand_daily / median_fw_daily else: scale = 1.0 elif phase == 'launch': actual_sold = batch_data['launch_sales'].get(pid, 0) age = max(0, product.get('age_days') or 0) if median_fw > 0 and actual_sold > 0: days_observed = min(age, 14) if days_observed > 0: projected_first_week = (actual_sold / days_observed) * 7 scale = projected_first_week / median_fw else: scale = 1.0 else: scale = 1.0 elif phase == 'decay': actual_velocity = batch_data['decay_velocity'].get(pid, 0) age = max(0, product.get('age_days') or 0) t_weeks = age / 7.0 # With baseline fix: value = (A/7)*exp(-λt)*scale + C/7 # Solve for scale: scale = (actual - C/7) / ((A/7)*exp(-λt)) decay_part = (amplitude / 7.0) * np.exp(-decay_rate * t_weeks) # Use a higher floor for the denominator at high ages to prevent # extreme scale factors when the decay envelope is nearly zero min_decay = max(0.01, amplitude / 70.0) # at least 10% of week-1 daily value if decay_part > min_decay and actual_velocity > 0: scale = (actual_velocity - baseline / 7.0) / decay_part elif actual_velocity > 0: scale = 1.0 else: scale = 1.0 else: scale = 1.0 # Clamp to avoid extreme values — tighter for preorder since the signal # is noisier (pre-orders accumulate differently than post-launch sales) max_scale = 5.0 if phase == 'preorder' else 8.0 return max(0.1, min(scale, max_scale)) # --------------------------------------------------------------------------- # Mature product forecast (Holt's double exponential smoothing) # --------------------------------------------------------------------------- def forecast_mature(product, history_df): """ Forecast for a mature/evergreen product using Holt's linear trend method on recent daily sales history. Holt's adds a trend component over SES, so it naturally pulls the forecast back down after a sales spike instead of persisting the inflated level. Falls back to SES then flat velocity on failure. """ pid = int(product['pid']) velocity = product.get('sales_velocity_daily') or 0 if history_df is None or history_df.empty or len(history_df) < 7: # Not enough data — flat velocity return np.full(FORECAST_HORIZON_DAYS, velocity) # Reindex over the FULL calendar window ending yesterday, not just the span # between the first and last snapshot. resample() only covers first→last # snapshot, so leading/trailing quiet periods are absent and the Holt level # is fitted only on the product's busy span (can run ~4x too high). An # explicit reindex fills every quiet calendar day with 0. (pid, snapshot_date) # is unique so there is no duplicate-index risk; do NOT use combine_first # (it keeps zeros over real data). See FORECAST_FIX_PLAN F2. hist = history_df.copy() hist['snapshot_date'] = pd.to_datetime(hist['snapshot_date']) hist = hist.set_index('snapshot_date')['units_sold'] full_index = pd.date_range( end=pd.Timestamp(date.today() - timedelta(days=1)), periods=EXP_SMOOTHING_WINDOW, freq='D') series = hist.reindex(full_index, fill_value=0.0).values.astype(float) # Need at least 2 non-zero values for smoothing if np.count_nonzero(series) < 2: return np.full(FORECAST_HORIZON_DAYS, velocity) # Cap: prevent runaway forecasts from one-time spikes. # Use the higher of 30d velocity or the observed mean as the baseline, # so sustained increases are respected. observed_mean = float(np.mean(series)) cap = max(velocity, observed_mean) * MAX_SMOOTHING_MULTIPLIER try: # Holt's with damped trend: the phi parameter dampens the trend over # the horizon so forecasts converge to a level instead of extrapolating # a linear trend indefinitely. model = Holt(series, initialization_method='estimated', damped_trend=True) fit = model.fit(optimized=True) forecast = fit.forecast(FORECAST_HORIZON_DAYS) forecast = np.clip(forecast, 0, cap) return forecast except Exception: # Fall back to SES if Holt's fails (e.g. insufficient data points) try: model = SimpleExpSmoothing(series, initialization_method='estimated') fit = model.fit(optimized=True) forecast = fit.forecast(FORECAST_HORIZON_DAYS) forecast = np.clip(forecast, 0, cap) return forecast except Exception as e: log.debug(f"ExpSmoothing failed for pid {pid}: {e}") return np.full(FORECAST_HORIZON_DAYS, velocity) def forecast_dormant(): """Dormant products get near-zero forecast.""" return np.zeros(FORECAST_HORIZON_DAYS) # --------------------------------------------------------------------------- # Accuracy-driven confidence margins # --------------------------------------------------------------------------- DEFAULT_MARGINS = { 'preorder': 0.4, 'launch': 0.35, 'decay': 0.3, 'mature': 0.35, 'slow_mover': 0.5, 'dormant': 0.5, } MIN_MARGIN = 0.15 # intervals shouldn't be tighter than ±15% MAX_MARGIN = 1.0 # intervals shouldn't exceed ±100% def load_accuracy_margins(conn): """ Load per-phase WMAPE from the most recent forecast accuracy run. Returns a dict of phase -> base_margin, falling back to DEFAULT_MARGINS. WMAPE is already a ratio (e.g. 1.7 = 170%), which we use directly as margin. """ margins = dict(DEFAULT_MARGINS) try: df = execute_query(conn, """ SELECT fa.dimension_value AS phase, fa.wmape FROM forecast_accuracy fa JOIN forecast_runs fr ON fr.id = fa.run_id WHERE fa.metric_type = 'by_phase' AND fr.status IN ('completed', 'backfill') AND fa.wmape IS NOT NULL ORDER BY fr.finished_at DESC """) if df.empty: log.info("No accuracy data available, using default confidence margins") return margins # Take the most recent run's values (they appear first due to ORDER BY) seen = set() for _, row in df.iterrows(): phase = row['phase'] if phase not in seen: wmape = float(row['wmape']) margins[phase] = max(MIN_MARGIN, min(wmape, MAX_MARGIN)) seen.add(phase) log.info(f"Loaded accuracy-based margins: {', '.join(f'{k}={v:.2f}' for k, v in margins.items())}") except Exception as e: log.warning(f"Could not load accuracy margins, using defaults: {e}") return margins # --------------------------------------------------------------------------- # Main orchestration # --------------------------------------------------------------------------- FLUSH_EVERY_PRODUCTS = 5000 # Flush forecast rows to DB every N products def generate_all_forecasts(conn, curves_df, dow_indices, monthly_indices=None, accuracy_margins=None): """Classify all products, batch-load data, generate and stream-write forecasts. Writes forecast rows to product_forecasts in chunks to avoid accumulating millions of rows in memory (37K products × 90 days = 3.3M rows). """ if monthly_indices is None: monthly_indices = {m: 1.0 for m in range(1, 13)} if accuracy_margins is None: accuracy_margins = dict(DEFAULT_MARGINS) log.info("Loading products for classification...") products = load_products(conn) log.info(f"Loaded {len(products)} visible products") # Classify each product products['phase'] = products.apply(classify_phase, axis=1) phase_counts = products['phase'].value_counts().to_dict() log.info(f"Phase distribution: {phase_counts}") # Batch-load per-product data (replaces per-product queries) log.info("Batch loading product data...") batch_data = batch_load_product_data(conn, products) today = date.today() forecast_dates = [today + timedelta(days=i) for i in range(FORECAST_HORIZON_DAYS)] # Pre-compute DOW and seasonal multipliers for each forecast date. # DOW multipliers stay ABSOLUTE — every calibration is a multi-week average # and therefore DOW-neutral, so reshaping by absolute DOW indices is correct. # Seasonal indices must be applied RELATIVE to the calibration period: # each per-product calibration (decay velocity, mature Holt level, launch / # preorder scale) is fitted on raw recent actuals that already embed the # current month's seasonal level. Multiplying by the absolute target-month # index double-counts seasonality (~25% over-forecast at the May→June sale # transition, worse near November). Divide by the trailing-30-day average # index so only the seasonal *change* from calibration to target applies. # See FORECAST_FIX_PLAN F3. dow_multipliers = [dow_indices.get(d.isoweekday(), 1.0) for d in forecast_dates] trailing = [today - timedelta(days=i) for i in range(1, 31)] calibration_index = float(np.mean([monthly_indices.get(d.month, 1.0) for d in trailing])) seasonal_multipliers = [ monthly_indices.get(d.month, 1.0) / max(calibration_index, 0.1) for d in forecast_dates ] # TRUNCATE before streaming writes with conn.cursor() as cur: cur.execute("TRUNCATE product_forecasts") conn.commit() buffer = [] methods = {} processed = 0 errors = 0 total_rows = 0 insert_sql = """ INSERT INTO product_forecasts (pid, forecast_date, forecast_units, forecast_revenue, lifecycle_phase, forecast_method, confidence_lower, confidence_upper) VALUES %s """ def flush_buffer(): nonlocal buffer, total_rows if not buffer: return with conn.cursor() as cur: psycopg2.extras.execute_values( cur, insert_sql, buffer, template="(%s, %s, %s, %s, %s, %s, %s, %s)", page_size=BATCH_SIZE, ) conn.commit() total_rows += len(buffer) buffer = [] for _, product in products.iterrows(): pid = int(product['pid']) phase = product['phase'] price = float(product['current_price'] or 0) age = max(0, product.get('age_days') or 0) try: curve_info = get_curve_for_product(product, curves_df) if phase == 'preorder': if curve_info: scale = compute_scale_factor('preorder', product, curve_info, batch_data) # Time the launch curve to expected arrival instead of # running it from today (F4). Pre-arrival days carry the # observed pre-order trickle rate. days_until_arrival = batch_data['preorder_arrival_days'].get(pid, 14) preorder_units = batch_data['preorder_sales'].get(pid, 0) preorder_days = batch_data['preorder_days'].get(pid, 1) preorder_daily_rate = preorder_units / max(preorder_days, 1) forecasts = forecast_preorder( curve_info, scale, days_until_arrival, preorder_daily_rate, FORECAST_HORIZON_DAYS) method = 'lifecycle_curve' else: # No reliable curve — fall back to velocity if available velocity = product.get('sales_velocity_daily') or 0 if velocity > 0: forecasts = np.full(FORECAST_HORIZON_DAYS, velocity) method = 'velocity' else: forecasts = forecast_dormant() method = 'zero' elif phase == 'launch': if curve_info: scale = compute_scale_factor('launch', product, curve_info, batch_data) forecasts = forecast_from_curve(curve_info, scale, age, FORECAST_HORIZON_DAYS) method = 'lifecycle_curve' else: # No reliable curve — fall back to velocity if available velocity = product.get('sales_velocity_daily') or 0 if velocity > 0: forecasts = np.full(FORECAST_HORIZON_DAYS, velocity) method = 'velocity' else: forecasts = forecast_dormant() method = 'zero' elif phase == 'decay': if curve_info: scale = compute_scale_factor(phase, product, curve_info, batch_data) forecasts = forecast_from_curve(curve_info, scale, age, FORECAST_HORIZON_DAYS) method = 'lifecycle_curve' else: velocity = product.get('sales_velocity_daily') or 0 forecasts = np.full(FORECAST_HORIZON_DAYS, velocity) method = 'velocity' elif phase == 'mature': history = batch_data['mature_history'].get(pid) forecasts = forecast_mature(product, history) method = 'exp_smoothing' elif phase == 'slow_mover': velocity = product.get('sales_velocity_daily') or 0 forecasts = np.full(FORECAST_HORIZON_DAYS, velocity) method = 'velocity' else: # dormant # Carry a small positive rate for dormant products that still # trickle sales (restocks/promos/long-tail); only truly dead # products stay at zero. See FORECAST_FIX_PLAN F5. rate = batch_data['dormant_rate'].get(pid, 0) if rate > 0: forecasts = np.full(FORECAST_HORIZON_DAYS, rate) method = 'velocity' else: forecasts = forecast_dormant() method = 'zero' # Confidence interval: use accuracy-calibrated margins per phase base_margin = accuracy_margins.get(phase, 0.5) for i, d in enumerate(forecast_dates): base_units = float(forecasts[i]) if i < len(forecasts) else 0.0 # Apply day-of-week and seasonal adjustments units = base_units * dow_multipliers[i] * seasonal_multipliers[i] # Widen confidence interval as horizon grows: day 0 = base, day 89 ≈ +50% wider horizon_factor = 1.0 + 0.5 * (i / max(FORECAST_HORIZON_DAYS - 1, 1)) margin = base_margin * horizon_factor buffer.append(( pid, d, round(units, 2), round(units * price, 4), phase, method, round(units * max(1 - margin, 0), 2), round(units * (1 + margin), 2), )) methods[method] = methods.get(method, 0) + 1 except Exception as e: log.warning(f"Error forecasting pid {pid}: {e}") errors += 1 # Write zero forecast so we have complete coverage for d in forecast_dates: buffer.append((pid, d, 0, 0, phase, 'zero', 0, 0)) processed += 1 if processed % FLUSH_EVERY_PRODUCTS == 0: flush_buffer() log.info(f" Processed {processed}/{len(products)} products ({total_rows} rows written)...") # Final flush flush_buffer() log.info(f"Forecast generation complete. {processed} products, {errors} errors, {total_rows} rows") log.info(f"Method distribution: {methods}") return total_rows, processed, phase_counts def archive_forecasts(conn, run_id): """ Copy current product_forecasts into history before they get replaced. Only archives forecast rows for dates that have already passed, so we can later compare them against actuals. """ with conn.cursor() as cur: # Ensure history table exists cur.execute(""" CREATE TABLE IF NOT EXISTS product_forecasts_history ( run_id INT NOT NULL, pid BIGINT NOT NULL, forecast_date DATE NOT NULL, forecast_units NUMERIC(10,2), forecast_revenue NUMERIC(14,4), lifecycle_phase TEXT, forecast_method TEXT, confidence_lower NUMERIC(10,2), confidence_upper NUMERIC(10,2), generated_at TIMESTAMP, PRIMARY KEY (run_id, pid, forecast_date) ) """) cur.execute("CREATE INDEX IF NOT EXISTS idx_pfh_date ON product_forecasts_history(forecast_date)") cur.execute("CREATE INDEX IF NOT EXISTS idx_pfh_pid_date ON product_forecasts_history(pid, forecast_date)") # Naive-baseline column for forecast value-added (FVA). See FORECAST_FIX_PLAN F8. cur.execute("ALTER TABLE product_forecasts_history ADD COLUMN IF NOT EXISTS naive_units NUMERIC(10,2)") # Find the previous completed run (whose forecasts are still in product_forecasts) cur.execute(""" SELECT id FROM forecast_runs WHERE status = 'completed' ORDER BY finished_at DESC LIMIT 1 """) prev_run = cur.fetchone() if prev_run is None: log.info("No previous completed run found, skipping archive") conn.commit() return 0 prev_run_id = prev_run[0] # Archive only past-date forecasts (where actuals now exist). Attach the # naive baseline (flat trailing-28-day daily average) at the same time so # forecast value-added can be measured. See FORECAST_FIX_PLAN F8. cur.execute(""" INSERT INTO product_forecasts_history (run_id, pid, forecast_date, forecast_units, forecast_revenue, lifecycle_phase, forecast_method, confidence_lower, confidence_upper, generated_at, naive_units) SELECT %s, pf.pid, pf.forecast_date, pf.forecast_units, pf.forecast_revenue, pf.lifecycle_phase, pf.forecast_method, pf.confidence_lower, pf.confidence_upper, pf.generated_at, COALESCE(nv.naive_daily, 0) FROM product_forecasts pf LEFT JOIN ( SELECT o.pid, SUM(o.quantity) / 28.0 AS naive_daily FROM orders o WHERE o.canceled IS DISTINCT FROM TRUE AND o.date >= CURRENT_DATE - INTERVAL '28 days' AND o.date < CURRENT_DATE GROUP BY o.pid ) nv ON nv.pid = pf.pid WHERE pf.forecast_date < CURRENT_DATE ON CONFLICT (run_id, pid, forecast_date) DO NOTHING """, (prev_run_id,)) archived = cur.rowcount conn.commit() if archived > 0: log.info(f"Archived {archived} historical forecast rows from run {prev_run_id}") else: log.info("No past-date forecasts to archive") # Prune old history (keep 90 days for accuracy analysis) cur.execute("DELETE FROM product_forecasts_history WHERE forecast_date < CURRENT_DATE - INTERVAL '90 days'") pruned = cur.rowcount if pruned > 0: log.info(f"Pruned {pruned} old history rows (>90 days)") conn.commit() return archived def archive_future_leads(conn, run_id): """ Archive a sampled set of FUTURE-lead forecasts from the just-generated product_forecasts, attributed to the current run. The past-date archive in archive_forecasts() only ever captures the 1-day slice that just elapsed, so every accuracy sample lands in the '1-7d' lead bucket and the 15/30/60/90-day forecasts that purchasing actually rides on are never validated. Here we snapshot the 7/14/30/60/89-day-ahead leads (non-dormant) so that, once each date passes, compute_accuracy() can score them in their lead bucket. The naive baseline is attached the same way as in the past-date path. Future-dated rows survive the 90-day prune until their own date passes. See FORECAST_FIX_PLAN F7. """ with conn.cursor() as cur: cur.execute(""" INSERT INTO product_forecasts_history (run_id, pid, forecast_date, forecast_units, forecast_revenue, lifecycle_phase, forecast_method, confidence_lower, confidence_upper, generated_at, naive_units) SELECT %s, pf.pid, pf.forecast_date, pf.forecast_units, pf.forecast_revenue, pf.lifecycle_phase, pf.forecast_method, pf.confidence_lower, pf.confidence_upper, pf.generated_at, COALESCE(nv.naive_daily, 0) FROM product_forecasts pf LEFT JOIN ( SELECT o.pid, SUM(o.quantity) / 28.0 AS naive_daily FROM orders o WHERE o.canceled IS DISTINCT FROM TRUE AND o.date >= CURRENT_DATE - INTERVAL '28 days' AND o.date < CURRENT_DATE GROUP BY o.pid ) nv ON nv.pid = pf.pid WHERE pf.lifecycle_phase != 'dormant' AND pf.forecast_date - CURRENT_DATE IN (7, 14, 30, 60, 89) ON CONFLICT (run_id, pid, forecast_date) DO NOTHING """, (run_id,)) archived = cur.rowcount conn.commit() log.info(f"Archived {archived} future-lead forecast rows (7/14/30/60/89d) for run {run_id}") return archived def compute_accuracy(conn, run_id): """ Compute forecast accuracy metrics from archived history vs. actual sales. Joins product_forecasts_history with daily_product_snapshots on (pid, forecast_date = snapshot_date) to compare forecasted vs. actual units. Stores results in forecast_accuracy table, broken down by: - overall: two rows — 'all' (non-dormant) and 'all_incl_dormant' (F5) - overall_weekly: per-product weekly-grain WMAPE — the informative headline for intermittent demand (daily grain has a ~190% floor) (F9) - by_phase: per lifecycle phase - by_lead_time: bucketed by how far ahead the forecast was — long-lead buckets populate as the future-lead archives mature (F7) - by_method: per forecast method - daily: per forecast_date (for trend charts) Every dimension also stores naive_wmape (flat trailing-28d baseline) and fva = 1 - wmape/naive_wmape, so the engine can be judged as value-over-naive (F8). Only realized dates (forecast_date < CURRENT_DATE) are scored. """ with conn.cursor() as cur: # Ensure accuracy table exists cur.execute(""" CREATE TABLE IF NOT EXISTS forecast_accuracy ( run_id INT NOT NULL, metric_type TEXT NOT NULL, dimension_value TEXT NOT NULL, sample_size INT, total_actual_units NUMERIC(12,2), total_forecast_units NUMERIC(12,2), mae NUMERIC(10,4), wmape NUMERIC(10,4), bias NUMERIC(10,4), rmse NUMERIC(10,4), computed_at TIMESTAMP NOT NULL DEFAULT NOW(), PRIMARY KEY (run_id, metric_type, dimension_value) ) """) # Naive-baseline WMAPE and forecast value-added (FVA = 1 - wmape/naive_wmape). # See FORECAST_FIX_PLAN F8. cur.execute("ALTER TABLE forecast_accuracy ADD COLUMN IF NOT EXISTS naive_wmape NUMERIC(10,4)") cur.execute("ALTER TABLE forecast_accuracy ADD COLUMN IF NOT EXISTS fva NUMERIC(10,4)") conn.commit() # Check if we have any history to analyze cur.execute("SELECT COUNT(*) FROM product_forecasts_history") history_count = cur.fetchone()[0] if history_count == 0: log.info("No forecast history available for accuracy computation") return # Base CTEs (FORECAST_FIX_PLAN F7): # - Only score realized dates (forecast_date < CURRENT_DATE); future-lead # archives are excluded until their date passes. # - short_lead*: lead 0-6 deduped per (pid, forecast_date) — preserves the # meaning of the existing headline metrics. short_lead_eval keeps the # raw snapshot grid (incl. zero-zero days) for complete-week detection; # `accuracy` drops zero-zero days for daily-grain metrics. # - lead_dedup/lead_accuracy: deduped per (pid, forecast_date, lead_bucket) # so each long-lead bucket gets its own sample (the by_lead_time table). base_cte = """ WITH ranked_all AS ( SELECT pfh.pid, pfh.forecast_date, pfh.forecast_units, pfh.naive_units, pfh.lifecycle_phase, pfh.forecast_method, fr.started_at, (pfh.forecast_date - fr.started_at::date) AS lead_days, CASE WHEN (pfh.forecast_date - fr.started_at::date) BETWEEN 0 AND 6 THEN '1-7d' WHEN (pfh.forecast_date - fr.started_at::date) BETWEEN 7 AND 13 THEN '8-14d' WHEN (pfh.forecast_date - fr.started_at::date) BETWEEN 14 AND 29 THEN '15-30d' WHEN (pfh.forecast_date - fr.started_at::date) BETWEEN 30 AND 59 THEN '31-60d' ELSE '61-90d' END AS lead_bucket FROM product_forecasts_history pfh JOIN forecast_runs fr ON fr.id = pfh.run_id WHERE pfh.forecast_date < CURRENT_DATE ), short_lead AS ( SELECT *, ROW_NUMBER() OVER ( PARTITION BY pid, forecast_date ORDER BY started_at DESC ) AS rn FROM ranked_all WHERE lead_days BETWEEN 0 AND 6 ), short_lead_eval AS ( SELECT sl.pid, sl.lifecycle_phase, sl.forecast_method, sl.forecast_date, sl.forecast_units, sl.naive_units, COALESCE(dps.units_sold, 0) AS actual_units, (sl.forecast_units - COALESCE(dps.units_sold, 0)) AS error, ABS(sl.forecast_units - COALESCE(dps.units_sold, 0)) AS abs_error FROM short_lead sl LEFT JOIN daily_product_snapshots dps ON dps.pid = sl.pid AND dps.snapshot_date = sl.forecast_date WHERE sl.rn = 1 ), accuracy AS ( SELECT * FROM short_lead_eval WHERE NOT (forecast_units = 0 AND actual_units = 0) ), lead_dedup AS ( SELECT *, ROW_NUMBER() OVER ( PARTITION BY pid, forecast_date, lead_bucket ORDER BY started_at DESC ) AS rn FROM ranked_all ), lead_accuracy AS ( SELECT ld.lead_bucket, ld.forecast_units, ld.naive_units, COALESCE(dps.units_sold, 0) AS actual_units, (ld.forecast_units - COALESCE(dps.units_sold, 0)) AS error, ABS(ld.forecast_units - COALESCE(dps.units_sold, 0)) AS abs_error FROM lead_dedup ld LEFT JOIN daily_product_snapshots dps ON dps.pid = ld.pid AND dps.snapshot_date = ld.forecast_date WHERE ld.rn = 1 AND ld.lifecycle_phase != 'dormant' AND NOT (ld.forecast_units = 0 AND COALESCE(dps.units_sold, 0) = 0) ) """ # Daily-grain aggregate over a source CTE aliased `a`, computing the # engine WMAPE plus the naive-baseline WMAPE (NULL-safe: rows archived # before F8 have naive_units NULL and are excluded from the naive sums). def daily_agg(dim_expr, source, where=None, group_by=None): where_sql = f"WHERE {where}" if where else "" group_sql = f"GROUP BY {group_by}" if group_by else "" return f""" SELECT {dim_expr} AS dim, COUNT(*) AS sample_size, COALESCE(SUM(a.actual_units), 0) AS total_actual, COALESCE(SUM(a.forecast_units), 0) AS total_forecast, AVG(a.abs_error) AS mae, CASE WHEN SUM(a.actual_units) > 0 THEN SUM(a.abs_error) / SUM(a.actual_units) ELSE NULL END AS wmape, AVG(a.error) AS bias, SQRT(AVG(POWER(a.error, 2))) AS rmse, CASE WHEN SUM(a.actual_units) FILTER (WHERE a.naive_units IS NOT NULL) > 0 THEN SUM(ABS(a.naive_units - a.actual_units)) FILTER (WHERE a.naive_units IS NOT NULL) / SUM(a.actual_units) FILTER (WHERE a.naive_units IS NOT NULL) ELSE NULL END AS naive_wmape FROM {source} a {where_sql} {group_sql} """ insert_sql = """ INSERT INTO forecast_accuracy (run_id, metric_type, dimension_value, sample_size, total_actual_units, total_forecast_units, mae, wmape, bias, rmse, naive_wmape, fva) VALUES (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s) ON CONFLICT (run_id, metric_type, dimension_value) DO UPDATE SET sample_size = EXCLUDED.sample_size, total_actual_units = EXCLUDED.total_actual_units, total_forecast_units = EXCLUDED.total_forecast_units, mae = EXCLUDED.mae, wmape = EXCLUDED.wmape, bias = EXCLUDED.bias, rmse = EXCLUDED.rmse, naive_wmape = EXCLUDED.naive_wmape, fva = EXCLUDED.fva, computed_at = NOW() """ def _f(x): return float(x) if x is not None else None def run_and_insert(metric_type, sql): cur.execute(base_cte + sql) n = 0 for row in cur.fetchall(): (dim_val, sample_size, total_actual, total_forecast, mae, wmape, bias, rmse, naive_wmape) = row fva = None if wmape is not None and naive_wmape is not None and float(naive_wmape) > 0: fva = 1.0 - float(wmape) / float(naive_wmape) cur.execute(insert_sql, ( run_id, metric_type, dim_val, sample_size, _f(total_actual), _f(total_forecast), _f(mae), _f(wmape), _f(bias), _f(rmse), _f(naive_wmape), _f(fva))) n += 1 return n total_inserted = 0 # overall: two rows — 'all' (non-dormant, the headline) and # 'all_incl_dormant' (everything, so the ~11% dormant demand stops being # invisible). Both are short-lead (lead 0-6). F5. overall_source = """( SELECT a.*, 'all'::text AS dim FROM accuracy a WHERE a.lifecycle_phase != 'dormant' UNION ALL SELECT a.*, 'all_incl_dormant'::text AS dim FROM accuracy a )""" total_inserted += run_and_insert('overall', daily_agg('a.dim', overall_source, group_by='a.dim')) # by_phase / by_method / daily — short-lead daily-grain over `accuracy`. total_inserted += run_and_insert('by_phase', daily_agg('a.lifecycle_phase', 'accuracy', group_by='a.lifecycle_phase')) total_inserted += run_and_insert('by_method', daily_agg('a.forecast_method', 'accuracy', group_by='a.forecast_method')) total_inserted += run_and_insert('daily', daily_agg('a.forecast_date::text', 'accuracy', where="a.lifecycle_phase != 'dormant'", group_by='a.forecast_date')) # by_lead_time — one sample per (pid, date, lead bucket) over `lead_accuracy`. # Buckets beyond '1-7d' populate as the future-lead archives (F7) mature. total_inserted += run_and_insert('by_lead_time', daily_agg('a.lead_bucket', 'lead_accuracy', group_by='a.lead_bucket')) # overall_weekly — the informative headline for intermittent retail demand. # Aggregate the short-lead rows to (pid, complete week), then WMAPE over # pid-weeks. Daily-grain WMAPE has a ~190% floor on this catalog; weekly # grain is ~109% and responds to real improvement. F9. weekly_sql = """, weekly AS ( SELECT pid, date_trunc('week', forecast_date) AS wk, SUM(forecast_units) AS fc_week, SUM(actual_units) AS act_week, SUM(naive_units) AS naive_week, bool_and(naive_units IS NOT NULL) AS naive_complete FROM short_lead_eval WHERE lifecycle_phase != 'dormant' GROUP BY pid, date_trunc('week', forecast_date) HAVING COUNT(*) = 7 ) SELECT 'all'::text AS dim, COUNT(*) AS sample_size, COALESCE(SUM(act_week), 0) AS total_actual, COALESCE(SUM(fc_week), 0) AS total_forecast, AVG(ABS(fc_week - act_week)) AS mae, CASE WHEN SUM(act_week) > 0 THEN SUM(ABS(fc_week - act_week)) / SUM(act_week) ELSE NULL END AS wmape, AVG(fc_week - act_week) AS bias, SQRT(AVG(POWER(fc_week - act_week, 2))) AS rmse, CASE WHEN SUM(act_week) FILTER (WHERE naive_complete) > 0 THEN SUM(ABS(naive_week - act_week)) FILTER (WHERE naive_complete) / SUM(act_week) FILTER (WHERE naive_complete) ELSE NULL END AS naive_wmape FROM weekly WHERE NOT (fc_week = 0 AND act_week = 0) """ total_inserted += run_and_insert('overall_weekly', weekly_sql) conn.commit() # Prune old accuracy data (keep 90 days of runs + any in-progress run) cur.execute(""" DELETE FROM forecast_accuracy WHERE run_id NOT IN ( SELECT id FROM forecast_runs WHERE finished_at >= NOW() - INTERVAL '90 days' OR finished_at IS NULL ) """) pruned = cur.rowcount conn.commit() log.info(f"Accuracy metrics: {total_inserted} rows computed" + (f", {pruned} old rows pruned" if pruned > 0 else "")) def backfill_accuracy_data(conn, backfill_days=30): """ Generate retroactive forecast data for the past N days to bootstrap accuracy metrics. Uses the current brand curves with per-product scaling to approximate what the model would have predicted for each past day, then stores results in product_forecasts_history for comparison against actual snapshots. This is a model backtest (in-sample), not true out-of-sample accuracy, but provides much better initial estimates than unscaled brand curves. """ backfill_start_time = time.time() log.info(f"Backfilling {backfill_days} days of accuracy data with per-product scaling...") # Load DOW indices dow_indices = compute_dow_indices(conn) # Load brand curves (already fitted) curves_df = execute_query(conn, """ SELECT brand, root_category, cat_id, category_level, amplitude, decay_rate, baseline, r_squared, median_first_week_sales, median_preorder_sales, median_preorder_days FROM brand_lifecycle_curves """) # Load products products = load_products(conn) products['phase'] = products.apply(classify_phase, axis=1) # Skip dormant — they forecast 0 and are filtered from accuracy anyway active = products[products['phase'] != 'dormant'].copy() log.info(f"Backfilling for {len(active)} non-dormant products " f"(skipping {len(products) - len(active)} dormant)") # Batch load product data for per-product scaling batch_data = batch_load_product_data(conn, active) today = date.today() backfill_start = today - timedelta(days=backfill_days) # Create a synthetic run entry with conn.cursor() as cur: cur.execute(""" INSERT INTO forecast_runs (started_at, finished_at, status, products_forecast, phase_counts, error_message) VALUES (%s, NOW(), 'backfill', %s, %s, %s) RETURNING id """, ( backfill_start, len(active), json.dumps({'backfill_days': backfill_days}), f'Model backtest: {backfill_days} days with per-product scaling', )) backfill_run_id = cur.fetchone()[0] conn.commit() log.info(f"Created backfill run {backfill_run_id} " f"(simulated start: {backfill_start})") # Generate retroactive forecasts all_rows = [] backfill_dates = [backfill_start + timedelta(days=i) for i in range(backfill_days)] for _, product in active.iterrows(): pid = int(product['pid']) price = float(product['current_price'] or 0) current_age = product.get('age_days') velocity = float(product.get('sales_velocity_daily') or 0) phase = product['phase'] curve_info = get_curve_for_product(product, curves_df) # Compute per-product scale factor (same logic as main forecast) scale = compute_scale_factor(phase, product, curve_info, batch_data) for forecast_date in backfill_dates: # How many days ago was this date? days_ago = (today - forecast_date).days # Product's age on that date past_age = (current_age - days_ago) if current_age is not None else None if past_age is not None and past_age < 0: # Product didn't exist yet on this date continue # Determine what phase the product was likely in if past_age is not None: if past_age <= LAUNCH_AGE_DAYS: past_phase = 'launch' elif past_age <= DECAY_AGE_DAYS: past_phase = 'decay' else: past_phase = phase # use current classification else: past_phase = phase # Compute forecast value for this date if past_phase in ('launch', 'decay', 'preorder') and curve_info: amplitude, decay_rate, baseline = curve_info[:3] age_for_calc = max(0, past_age or 0) t_weeks = age_for_calc / 7.0 # Use corrected formula: scale only the decay envelope, not the baseline daily_value = (amplitude / 7.0) * np.exp(-decay_rate * t_weeks) * scale + (baseline / 7.0) units = max(0.0, float(daily_value)) method = 'lifecycle_curve' elif past_phase == 'mature' and velocity > 0: units = velocity method = 'exp_smoothing' else: units = velocity if velocity > 0 else 0.0 method = 'velocity' if velocity > 0 else 'zero' # Apply DOW multiplier dow_mult = dow_indices.get(forecast_date.isoweekday(), 1.0) units *= dow_mult if units == 0 and method == 'zero': continue # skip zero-zero rows revenue = units * price margin = 0.3 if method == 'lifecycle_curve' else 0.4 all_rows.append(( backfill_run_id, pid, forecast_date, round(float(units), 2), round(float(revenue), 4), past_phase, method, round(float(units * (1 - margin)), 2), round(float(units * (1 + margin)), 2), backfill_start, # generated_at )) log.info(f"Generated {len(all_rows)} backfill forecast rows") # Write to history table if all_rows: with conn.cursor() as cur: sql = """ INSERT INTO product_forecasts_history (run_id, pid, forecast_date, forecast_units, forecast_revenue, lifecycle_phase, forecast_method, confidence_lower, confidence_upper, generated_at) VALUES %s ON CONFLICT (run_id, pid, forecast_date) DO NOTHING """ psycopg2.extras.execute_values( cur, sql, all_rows, template="(%s, %s, %s, %s, %s, %s, %s, %s, %s, %s)", page_size=BATCH_SIZE, ) conn.commit() log.info(f"Wrote {len(all_rows)} rows to product_forecasts_history") # Now compute accuracy on the backfilled data compute_accuracy(conn, backfill_run_id) # Mark the backfill run as completed backfill_duration = time.time() - backfill_start_time with conn.cursor() as cur: cur.execute(""" UPDATE forecast_runs SET finished_at = NOW(), status = 'backfill', duration_seconds = %s WHERE id = %s """, (round(backfill_duration, 2), backfill_run_id)) conn.commit() log.info(f"Backfill complete in {backfill_duration:.1f}s") return backfill_run_id def main(): start_time = time.time() conn = get_connection() # Clean up any stale "running" entries from prior crashes cleanup_stale_runs(conn) # Check for --backfill flag if '--backfill' in sys.argv: idx = sys.argv.index('--backfill') days = int(sys.argv[idx + 1]) if idx + 1 < len(sys.argv) else 30 log.info("=" * 60) log.info(f"Backfill mode: {days} days") log.info("=" * 60) try: backfill_accuracy_data(conn, days) finally: conn.close() return log.info("=" * 60) log.info("Forecast Engine starting") log.info("=" * 60) run_id = None try: # Record run start with conn.cursor() as cur: cur.execute( "INSERT INTO forecast_runs (started_at, status) VALUES (NOW(), 'running') RETURNING id" ) run_id = cur.fetchone()[0] conn.commit() # Phase 0: Compute day-of-week and monthly seasonal indices dow_indices = compute_dow_indices(conn) monthly_indices = compute_monthly_seasonal_indices(conn) # Phase 1: Build reference curves curves_df = build_reference_curves(conn) # Phase 2: Archive historical forecasts (before TRUNCATE in generation) archive_forecasts(conn, run_id) # Phase 3: Compute accuracy from archived history vs actuals compute_accuracy(conn, run_id) # Phase 3b: Load accuracy-calibrated confidence margins accuracy_margins = load_accuracy_margins(conn) # Phase 4: Generate and stream-write forecasts (TRUNCATE + chunked INSERT) total_rows, products_forecast, phase_counts = generate_all_forecasts( conn, curves_df, dow_indices, monthly_indices, accuracy_margins ) # Phase 4b: Snapshot sampled future-lead forecasts (7/14/30/60/89d) from # the fresh run so long-lead accuracy populates once those dates pass (F7). archive_future_leads(conn, run_id) duration = time.time() - start_time # Record run completion (include DOW indices in metadata) with conn.cursor() as cur: cur.execute(""" UPDATE forecast_runs SET finished_at = NOW(), status = 'completed', products_forecast = %s, phase_counts = %s, curve_count = %s, duration_seconds = %s WHERE id = %s """, ( products_forecast, json.dumps({ **phase_counts, '_dow_indices': {str(k): v for k, v in dow_indices.items()}, '_seasonal_indices': {str(k): v for k, v in monthly_indices.items()}, }), len(curves_df) if not curves_df.empty else 0, round(duration, 2), run_id, )) conn.commit() log.info("=" * 60) log.info(f"Forecast complete in {duration:.1f}s") log.info(f" Products: {products_forecast}") log.info(f" Curves: {len(curves_df) if not curves_df.empty else 0}") log.info(f" Phases: {phase_counts}") log.info(f" Rows written: {total_rows}") log.info("=" * 60) except Exception as e: duration = time.time() - start_time log.error(f"Forecast engine failed: {e}", exc_info=True) if run_id: try: with conn.cursor() as cur: cur.execute(""" UPDATE forecast_runs SET finished_at = NOW(), status = 'failed', error_message = %s, duration_seconds = %s WHERE id = %s """, (str(e), round(duration, 2), run_id)) conn.commit() except Exception: pass sys.exit(1) finally: conn.close() if __name__ == '__main__': main()