import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
import warnings
warnings.filterwarnings('ignore')

# Load real retail demand data from UCI Online Retail (2010-2011)
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))

# Aggregate total demand per day for the entire store
daily_demand = df.groupby(df['InvoiceDate'].dt.date)['Quantity'].sum()
print(daily_demand.head(7))

# Plot store-wide daily demand
plt.figure(figsize=(12,4))
plt.plot(daily_demand.index, daily_demand.values)
plt.title('Total Daily Demand (All Products)')
plt.ylabel('Units Sold')
plt.xlabel('Date')
plt.tight_layout()
plt.show()

# Select a commonly sold stock code (find the most frequent SKU)
sku_counts = df['StockCode'].value_counts()
top_sku = sku_counts.idxmax()
product = df[df['StockCode'] == top_sku]
product_daily = product.groupby(product['InvoiceDate'].dt.date)['Quantity'].sum()
print(f'Top SKU: {top_sku}')
print(product_daily.head(7))

# Plot daily demand for the most common product
plt.figure(figsize=(12,4))
plt.plot(product_daily.index, product_daily.values, color='orange')
plt.title(f'Daily Demand for SKU {top_sku}')
plt.ylabel('Units Sold')
plt.xlabel('Date')
plt.tight_layout()
plt.show()

# Fill missing dates with zero demand for the top SKU
min_date = product_daily.index.min()
max_date = product_daily.index.max()
all_dates = pd.date_range(min_date, max_date)
product_daily_filled = product_daily.reindex(all_dates, fill_value=0)
print(product_daily_filled.head(10))

# Apply seasonal_decompose (daily frequency assumed) to single product demand
decomp = seasonal_decompose(product_daily_filled, model='additive', period=7)
decomp.plot()
plt.suptitle(f'Decomposition of Demand for SKU {top_sku}', y=1.02)
plt.show()

# Highlight unusually high residuals (potential demand spikes)
resid = decomp.resid.fillna(0)
threshold = 2 * resid.std()
spikes = resid[abs(resid) > threshold]
print(spikes)

# Calculate 7-day rolling mean
rolling_avg = product_daily_filled.rolling(window=7).mean()
plt.figure(figsize=(12,4))
plt.plot(product_daily_filled.index, product_daily_filled.values, label='Actual', alpha=0.5)
plt.plot(rolling_avg.index, rolling_avg.values, label='7-Day Rolling Average', color='red')
plt.title(f'Smoothed Demand for SKU {top_sku}')
plt.legend()
plt.tight_layout()
plt.show()

# Fill missing dates for all-store demand
all_min_date = daily_demand.index.min()
all_max_date = daily_demand.index.max()
all_dates_full = pd.date_range(all_min_date, all_max_date)
daily_demand_filled = daily_demand.reindex(all_dates_full, fill_value=0)

# Decompose total store demand (additive, weekly pattern)
total_decomp = seasonal_decompose(daily_demand_filled, model='additive', period=7)
total_decomp.plot()
plt.suptitle('Decomposition of Total Store Demand', y=1.02)
plt.show()

# Combine and export trend and seasonality to CSV for downstream analytics
trend_seasonal = pd.DataFrame({
    'trend': total_decomp.trend.values,
    'seasonal': total_decomp.seasonal.values,
}, index=daily_demand_filled.index)
trend_seasonal.to_csv('demand_decomposition.csv')
print('File saved: demand_decomposition.csv')

# Find the second most common SKU
second_sku = sku_counts.index[1]
product2 = df[df['StockCode'] == second_sku]
product2_daily = product2.groupby(product2['InvoiceDate'].dt.date)['Quantity'].sum()
product2_filled = product2_daily.reindex(all_dates, fill_value=0)

# Decompose both SKUs and compare their seasonality
decomp1 = seasonal_decompose(product_daily_filled, model='additive', period=7)
decomp2 = seasonal_decompose(product2_filled, model='additive', period=7)

plt.figure(figsize=(12,5))
plt.plot(decomp1.seasonal.index, decomp1.seasonal.values, label=f'SKU {top_sku}')
plt.plot(decomp2.seasonal.index, decomp2.seasonal.values, label=f'SKU {second_sku}')
plt.title('Weekly Seasonality Patterns for Two SKUs')
plt.ylabel('Seasonal Effect')
plt.xlabel('Date')
plt.legend()
plt.tight_layout()
plt.show()

# Example: Attempting decomposition with missing dates (will show issues)
try:
    _ = seasonal_decompose(product_daily, model='additive', period=7)
except Exception as e:
    print('Error:', e)

# Identify and quantify negative sales (returns) in the data
returns = df[df['Quantity'] < 0]
print(f'Number of return transactions: {returns.shape[0]}')
print(returns[['InvoiceDate', 'Quantity', 'StockCode']].head())

# Create a demand-only version of the dataset (no returns)
demand_df = df[df['Quantity'] > 0]
daily_demand_demand_only = demand_df.groupby(demand_df['InvoiceDate'].dt.date)['Quantity'].sum()
print(daily_demand_demand_only.head())

# Calculate KPIs for daily demand (no negative sales)
filled = daily_demand_demand_only.reindex(pd.date_range(daily_demand_demand_only.index.min(), daily_demand_demand_only.index.max()), fill_value=0)
average_demand = filled.mean()
peak_demand = filled.max()
cv = filled.std() / average_demand  # coefficient of variation
print(f'Average daily demand: {average_demand:.2f}')
print(f'Peak daily demand: {peak_demand}')
print(f'Coefficient of variation: {cv:.2%}')

# Get last 26 weeks of store-wide daily demand (demand-only, filled)
recent = filled[-26*7:]
mean = recent.mean()
std = recent.std()
z = 1.65  # 95% service level (approximate z-score)
safety_stock = z * std
print(f'Estimated daily mean: {mean:.2f}')
print(f'Estimated daily std: {std:.2f}')
print(f'95% Service Level Safety Stock: {safety_stock:.2f}')