import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# Load real-world retail demand data
url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/00502/online_retail_II.xlsx'
df = pd.read_excel(url, sheet_name='Year 2010-2011')
df['InvoiceDate'] = pd.to_datetime(df['InvoiceDate'])
print(df.shape)
print(df.head(3))

# Quick data cleaning: remove cancellations and errors
df = df[df['Quantity'] > 0]
df = df[df['Invoice'].astype(str).str.startswith('C') == False]
print(df.shape)

# Select a single product for focused analysis
item_code = df['StockCode'].value_counts().idxmax()
item_df = df[df['StockCode'] == item_code]
print(f'Analyzing product {item_code}')

# Aggregate daily demand for the selected product
daily_sales = (item_df.groupby(item_df['InvoiceDate'].dt.date)['Quantity']
              .sum()
              .rename('sales'))
daily_sales = daily_sales.asfreq('D', fill_value=0)
print(daily_sales.head())

# Visualize daily demand
plt.figure(figsize=(12,4))
plt.plot(daily_sales, marker='o', linestyle='-')
plt.title('Daily Demand for Product ' + str(item_code))
plt.xlabel('Date')
plt.ylabel('Units Sold')
plt.show()

# Calculate and plot 7-day moving average
ma7 = daily_sales.rolling(window=7).mean()
plt.figure(figsize=(12,4))
plt.plot(daily_sales.index, daily_sales.values, label='Actual', alpha=0.5)
plt.plot(ma7.index, ma7.values, label='7-Day MA', linewidth=2)
plt.legend()
plt.title('7-Day Simple Moving Average Forecast')
plt.show()

# 14-day simple moving average forecast
ma14 = daily_sales.rolling(window=14).mean()
plt.figure(figsize=(12,4))
plt.plot(daily_sales.index, daily_sales.values, label='Actual', alpha=0.5)
plt.plot(ma14.index, ma14.values, label='14-Day MA', linewidth=2, color='orange')
plt.legend()
plt.title('14-Day Simple Moving Average Forecast')
plt.show()

# Create a true next-day moving average prediction
forecast = ma7.shift(1)
plt.figure(figsize=(12,4))
plt.plot(daily_sales.index, daily_sales.values, label='Actual', alpha=0.5)
plt.plot(forecast.index, forecast.values, label='Next-Day Moving Avg Forecast', linestyle='--')
plt.legend()
plt.title('Operational 7-Day Moving Average Forecast (Shifted)')
plt.show()

# Calculate MAE for 7-day moving average forecast
mae = np.mean(np.abs(daily_sales[7:] - forecast[7:]))
print(f'Mean Absolute Error: {mae:.2f} units per day')

from statsmodels.tsa.holtwinters import SimpleExpSmoothing
# Fit simple exponential smoothing with alpha=0.3
model = SimpleExpSmoothing(daily_sales).fit(smoothing_level=0.3, optimized=False)
exp_forecast = model.fittedvalues.shift(1)
plt.figure(figsize=(12,4))
plt.plot(daily_sales.index, daily_sales.values, label='Actual', alpha=0.5)
plt.plot(exp_forecast.index, exp_forecast.values, label='Exp Smoothing', linestyle='--', color='g')
plt.legend()
plt.title('Simple Exponential Smoothing Forecast (alpha=0.3)')
plt.show()

# Compute MAE for exponential smoothing
exp_mae = np.mean(np.abs(daily_sales[1:] - exp_forecast[1:]))
print(f'Exponential Smoothing MAE: {exp_mae:.2f} units per day')

# Try different alpha values and compare MAE
maes = []
alphas = np.arange(0.1, 1.0, 0.1)
for alpha in alphas:
    model = SimpleExpSmoothing(daily_sales).fit(smoothing_level=alpha, optimized=False)
    forecast = model.fittedvalues.shift(1)
    mae = np.mean(np.abs(daily_sales[1:] - forecast[1:]))
    maes.append(mae)
plt.figure(figsize=(8,4))
plt.plot(alphas, maes, marker='o')
plt.xlabel('Alpha (Smoothing Level)')
plt.ylabel('MAE')
plt.title('Forecast Error vs. Alpha Parameter')
plt.show()

# Aggregate forecast for weekly replenishment
weekly_fcst = forecast.resample('W').sum()
weekly_actual = daily_sales.resample('W').sum()
plt.figure(figsize=(12,4))
plt.plot(weekly_actual.index, weekly_actual.values, label='Actual Weekly Sales')
plt.plot(weekly_fcst.index, weekly_fcst.values, label='Forecasted Weekly Sales')
plt.title('Weekly Demand vs Weekly Forecast (7-Day MA)')
plt.xlabel('Week')
plt.ylabel('Units')
plt.legend()
plt.show()

# Fit exponential smoothing and let library tune alpha
auto_model = SimpleExpSmoothing(daily_sales).fit(optimized=True)
auto_fcst = auto_model.fittedvalues.shift(1)
auto_alpha = auto_model.model.params['smoothing_level']
auto_mae = np.mean(np.abs(daily_sales[1:] - auto_fcst[1:]))
print(f'Best-fit alpha: {auto_alpha:.2f}, MAE: {auto_mae:.2f}')

# Visualize optimal forecast vs. actual demand
plt.figure(figsize=(12,4))
plt.plot(daily_sales.index, daily_sales.values, label='Actual', alpha=0.5)
plt.plot(auto_fcst.index, auto_fcst.values, label='Optimized Exp Smoothing', linestyle='--', color='r')
plt.legend()
plt.title(f'Best Alpha ({auto_alpha:.2f})  Automated Exponential Smoothing Forecast')
plt.show()

# Identify missing dates in the demand series
all_days = pd.date_range(daily_sales.index.min(), daily_sales.index.max(), freq='D')
missing_days = all_days.difference(pd.to_datetime(daily_sales.index))
print(f'Missing days: {missing_days}')

# Sort and aggregate sales by date in case of issues
daily_sales_clean = daily_sales.groupby(daily_sales.index).sum().sort_index()
print(daily_sales_clean.head())

# Simulate an operational lead time of 3 days
lead_time = 3
future_forecast = ma7.shift(lead_time)
print('Lead-time adjusted forecast:')
print(future_forecast.head(10))

# Aggregate sales by month for KPI reporting
monthly_sales = daily_sales.resample('M').sum()
print(monthly_sales.head())

# Calculate a monthly growth rate KPI
monthly_growth = monthly_sales.pct_change().fillna(0)
print('Monthly demand growth rates (%):')
print((monthly_growth * 100).round(2))

# Practice: Insert fake missing holiday and fill
holiday = pd.to_datetime('2010-12-25')
if holiday not in daily_sales.index:
    daily_sales.loc[holiday] = 0
daily_sales = daily_sales.sort_index().asfreq('D', fill_value=0)
print(daily_sales.loc['2010-12-24':'2010-12-26'])

# Generate a forecast and compare to actual weekly sales
fcst = SimpleExpSmoothing(daily_sales).fit(optimized=True).fittedvalues.shift(1)
orders = fcst.resample('W').sum().round().astype(int)
actuals = daily_sales.resample('W').sum().astype(int)
order_df = pd.DataFrame({'OrderQTY': orders, 'ActualQTY': actuals})
order_df['MetDemand'] = np.minimum(order_df['OrderQTY'], order_df['ActualQTY'])
order_df['ServiceLevel'] = order_df['MetDemand'] / order_df['ActualQTY'].replace(0, np.nan)
print(order_df.head())

# Create a KPI report of average service level by month
order_df['Month'] = order_df.index.to_period('M')
monthly_service = order_df.groupby('Month')['ServiceLevel'].mean().round(3)
print('Average monthly service level:')
print(monthly_service)