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

# Beginner Example 1: Load the retail demand data (UCI Online Retail)
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))

# Beginner Example 2: Filter for a single item in a single country (SKU-level analysis)
item_code = df['StockCode'].value_counts().index[0]
country = 'United Kingdom'
df_item = df[(df['StockCode'] == item_code) & (df['Country'] == country)]
print(f'Selected item code: {item_code} from {country}')
print(df_item.head(3))

# Beginner Example 3: Aggregate item demand per day
daily_demand = df_item.groupby(df_item['InvoiceDate'].dt.date)['Quantity'].sum()
print(daily_demand.head(7))
print(f'Number of days in series: {daily_demand.shape[0]}')

# Intermediate Example 1: Handle missing demand days (fill zeros)
idx = pd.date_range(start=daily_demand.index.min(), end=daily_demand.index.max(), freq='D')
daily_demand_filled = daily_demand.reindex(idx, fill_value=0)
print(daily_demand_filled.head(10))
print(f'Number of expected days (with filled gaps): {daily_demand_filled.shape[0]}')

# Intermediate Example 2: Plot the demand profile
plt.figure(figsize=(10,4))
plt.plot(daily_demand_filled.index, daily_demand_filled.values, marker='o', linestyle='-', alpha=0.8)
plt.title('Daily demand for SKU ' + str(item_code))
plt.xlabel('Date')
plt.ylabel('Demand (units)')
plt.grid(True)
plt.tight_layout()
plt.show()

# Intermediate Example 3: Calculate average demand and standard deviation per day
mean_demand = daily_demand_filled.mean()
std_demand = daily_demand_filled.std()
print(f'Average daily demand: {mean_demand:.2f}')
print(f'Standard deviation of daily demand: {std_demand:.2f}')

# Intermediate Example 4: Estimate lead time (assume 7 days for demo purposes)
lead_time_days = 7
print(f'Assumed lead time (order to receipt) is {lead_time_days} days.')

# Intermediate Example 5: Calculate safety stock (basic model, fixed lead time)
service_level = 0.99  # Target 99% in-stock performance
from scipy.stats import norm
z = norm.ppf(service_level)
safety_stock = z * std_demand * np.sqrt(lead_time_days)
print(f'Calculated safety stock: {safety_stock:.2f} units for a {service_level*100:.0f}% service level.')

# Intermediate Example 6: Calculate average lead time demand (reorder point base)
avg_lead_time_demand = mean_demand * lead_time_days
print(f'Average lead time demand: {avg_lead_time_demand:.2f} units')

# Intermediate Example 7: Calculate reorder point
reorder_point = avg_lead_time_demand + safety_stock
print(f'Reorder point (ROP): {reorder_point:.2f} units')

# Advanced Example 1: Model variable lead times using a random sample from real-world data
np.random.seed(42)
# For demonstration, create a simulated lead time history
lead_time_history = np.random.choice([5,6,7,8,9], size=30, p=[0.1,0.2,0.4,0.2,0.1])
mean_lt, std_lt = np.mean(lead_time_history), np.std(lead_time_history)
print(f'Empirical mean lead time: {mean_lt:.2f} days, stdev: {std_lt:.2f} days')

# Advanced Example 2: Safety stock formula including both demand and lead time variability
combined_std = np.sqrt(lead_time_history.mean() * std_demand ** 2 + mean_demand ** 2 * np.var(lead_time_history))
safety_stock_var_lt = z * combined_std
print(f'Advanced safety stock (variable lead time): {safety_stock_var_lt:.2f} units')

# Advanced Example 3: Visualize actual demand histogram during lead time
window = int(round(mean_lt))
rolling_sums = daily_demand_filled.rolling(window).sum().dropna()
plt.figure(figsize=(8,4))
plt.hist(rolling_sums, bins=20, alpha=0.7, color='steelblue', edgecolor='black')
plt.title(f'Empirical {window}-day demand distribution')
plt.xlabel('Demand during lead time')
plt.ylabel('Frequency')
plt.show()

# Advanced Example 4: Export a summary report of recommended inventory parameters
summary = pd.DataFrame({
    'SKU': [item_code],
    'ServiceLevel': [service_level],
    'MeanDailyDemand': [mean_demand],
    'StdDailyDemand': [std_demand],
    'MeanLeadTime': [mean_lt],
    'StdLeadTime': [std_lt],
    'SafetyStock': [safety_stock_var_lt],
    'ReorderPoint': [mean_demand * mean_lt + safety_stock_var_lt]
})
summary_filename = 'inventory_parameters.csv'
summary.to_csv(summary_filename, index=False)
print('Inventory recommendation report saved as', summary_filename)

# Error handling: Detect missing or negative demand entries
missing_qty = df_item['Quantity'].isna().sum()
negative_qty = (df_item['Quantity'] < 0).sum()
print(f'Missing daily Quantity values: {missing_qty}')
print(f'Negative Quantity values: {negative_qty}')

# Error handling: Ensure all dates are expected type and continuous
if not np.issubdtype(df_item['InvoiceDate'].dtype, np.datetime64):
    print('InvoiceDate is not datetime! Convert it.')
else:
    print('InvoiceDate is already datetime type.')
expected_start = pd.to_datetime(daily_demand.index.min())
expected_end = pd.to_datetime(daily_demand.index.max())
all_days = pd.date_range(expected_start, expected_end, freq='D')
if daily_demand.index.shape[0] != all_days.shape[0]:
    print('Gap in daily dates detected! Fill missing days for correct analytics.')

# Error handling: Simulate a broken join on wrong keys
orders_df = df[['Invoice', 'InvoiceDate', 'StockCode', 'Quantity']]
inventories_df = pd.DataFrame({'StockCode': [item_code], 'InventoryOnHand': [50]})
# Simulate wrong join: joining on InvoiceCode instead of StockCode
joined_bad = pd.merge(orders_df, inventories_df, left_on='Invoice', right_on='StockCode', how='left')
print('Joined shape (bad key):', joined_bad.shape)
print('Number of matched rows:', joined_bad['InventoryOnHand'].notna().sum())

# Best practice: time-based demand aggregation (weekly)
weekly_demand = daily_demand_filled.resample('W-MON').sum()
print('First 4 weeks of demand:')
print(weekly_demand.head(4))

# Best practice: calculate in-stock (cycle service) and fill rate KPIs for a simulated policy
begin_inventory = reorder_point + 10
days_with_stockout = (daily_demand_filled.cumsum() > begin_inventory).sum()
service_level_emp = 1 - days_with_stockout / len(daily_demand_filled)
total_demand = daily_demand_filled.sum()
unfilled_demand = max(0, daily_demand_filled.cumsum().max() - begin_inventory)
fill_rate_emp = 1 - unfilled_demand / total_demand
print(f'Empirical cycle service level: {service_level_emp:.3f}')
print(f'Empirical fill rate: {fill_rate_emp:.3f}')

# Step 1: Get daily demand for a new SKU and fill missing days
new_item = df['StockCode'].value_counts().index[1]
df_new = df[(df['StockCode'] == new_item) & (df['Country'] == country)]
dd_new = df_new.groupby(df_new['InvoiceDate'].dt.date)['Quantity'].sum()
idx_new = pd.date_range(start=dd_new.index.min(), end=dd_new.index.max(), freq='D')
dd_new_filled = dd_new.reindex(idx_new, fill_value=0)

# Step 2: Repeat all calculations for safety stock and reorder point
mean_d = dd_new_filled.mean()
std_d = dd_new_filled.std()
lead_hist = np.random.choice([5,6,7,8,9], size=30, p=[0.1,0.2,0.4,0.2,0.1])
mean_lt2, std_lt2 = np.mean(lead_hist), np.std(lead_hist)
serv_level = 0.97
z2 = norm.ppf(serv_level)
comb_std2 = np.sqrt(lead_hist.mean() * std_d ** 2 + mean_d ** 2 * np.var(lead_hist))
safe_stock2 = z2 * comb_std2
rop2 = mean_d * mean_lt2 + safe_stock2
print(f'Item: {new_item}, Safety stock: {safe_stock2:.1f}, Reorder point: {rop2:.1f}')