inventory_simulator.py

"""
Inventory Simulator — EOQ / ROP
================================
Calculates optimal order quantity (EOQ) and reorder point (ROP)
for supply chain inventory management.

Author: Emmanuel Beristain Guzmán
GitHub: https://github.com/net421
"""

import math
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import os

# ─────────────────────────────────────────────
#  EOQ CALCULATION
# ─────────────────────────────────────────────

def calculate_eoq(annual_demand: float, ordering_cost: float, holding_cost: float) -> float:
    """
    Economic Order Quantity (EOQ) formula.

    Parameters:
        annual_demand  : Units demanded per year (D)
        ordering_cost  : Cost per order placed    (S)
        holding_cost   : Annual holding cost/unit  (H)

    Returns:
        EOQ in units (float)
    """
    if holding_cost <= 0:
        raise ValueError("Holding cost must be greater than zero.")
    eoq = math.sqrt((2 * annual_demand * ordering_cost) / holding_cost)
    return round(eoq, 2)


# ─────────────────────────────────────────────
#  ROP CALCULATION
# ─────────────────────────────────────────────

def calculate_rop(daily_demand: float, lead_time_days: float, safety_stock: float = 0) -> float:
    """
    Reorder Point (ROP) formula.

    Parameters:
        daily_demand    : Average daily demand (units)
        lead_time_days  : Supplier lead time in days
        safety_stock    : Optional buffer stock (units)

    Returns:
        ROP in units (float)
    """
    rop = (daily_demand * lead_time_days) + safety_stock
    return round(rop, 2)


# ─────────────────────────────────────────────
#  SAFETY STOCK CALCULATION
# ─────────────────────────────────────────────

def calculate_safety_stock(z_score: float, std_demand: float, lead_time_days: float) -> float:
    """
    Safety Stock = Z * σ_demand * sqrt(lead_time)

    Parameters:
        z_score         : Service level z-score (e.g. 1.65 = 95%)
        std_demand      : Standard deviation of daily demand
        lead_time_days  : Lead time in days

    Returns:
        Safety stock in units (float)
    """
    ss = z_score * std_demand * math.sqrt(lead_time_days)
    return round(ss, 2)


# ─────────────────────────────────────────────
#  TOTAL ANNUAL COST
# ─────────────────────────────────────────────

def total_annual_cost(annual_demand: float, eoq: float,
                      ordering_cost: float, holding_cost: float) -> float:
    """
    Total Annual Inventory Cost = Ordering Cost + Holding Cost

    Returns:
        Total cost in currency units (float)
    """
    orders_per_year = annual_demand / eoq
    avg_inventory   = eoq / 2
    cost = (orders_per_year * ordering_cost) + (avg_inventory * holding_cost)
    return round(cost, 2)


# ─────────────────────────────────────────────
#  SIMULATION: INVENTORY LEVEL OVER TIME
# ─────────────────────────────────────────────

def simulate_inventory(annual_demand: float, eoq: float, rop: float,
                       lead_time_days: int, days: int = 365) -> pd.DataFrame:
    """
    Simulates day-by-day inventory levels using EOQ/ROP logic.

    Returns:
        DataFrame with columns: day, inventory, order_placed, order_received
    """
    daily_demand    = annual_demand / 365
    inventory       = eoq        # start with one full order
    pending_order   = False
    order_arrival   = 0
    records         = []

    for day in range(1, days + 1):
        order_placed    = False
        order_received  = False

        # Receive order if it arrives today
        if pending_order and day == order_arrival:
            inventory      += eoq
            pending_order   = False
            order_received  = True

        # Consume daily demand (floor to int for realistic simulation)
        consumed   = min(inventory, daily_demand)
        inventory  = max(0, inventory - consumed)

        # Place new order if inventory hits ROP and no order pending
        if inventory <= rop and not pending_order:
            pending_order  = True
            order_arrival  = day + lead_time_days
            order_placed   = True

        records.append({
            "day":            day,
            "inventory":      round(inventory, 1),
            "order_placed":   order_placed,
            "order_received": order_received,
        })

    return pd.DataFrame(records)


# ─────────────────────────────────────────────
#  CHART: INVENTORY LEVEL OVER TIME
# ─────────────────────────────────────────────

def plot_inventory(df: pd.DataFrame, rop: float, eoq: float,
                   product_name: str = "Product", output_path: str = None):
    """Generates and saves the inventory level chart."""

    fig, ax = plt.subplots(figsize=(14, 6))
    fig.patch.set_facecolor("#0F1117")
    ax.set_facecolor("#0F1117")

    # Inventory line
    ax.plot(df["day"], df["inventory"],
            color="#00D4AA", linewidth=1.8, label="Inventory level", zorder=3)

    # ROP line
    ax.axhline(y=rop, color="#FF6B6B", linewidth=1.2,
               linestyle="--", label=f"ROP = {rop} units", zorder=2)

    # EOQ reference band
    ax.axhline(y=eoq, color="#FFD93D", linewidth=0.8,
               linestyle=":", alpha=0.5, label=f"EOQ = {eoq} units", zorder=2)

    # Order placed markers
    placed = df[df["order_placed"]]
    ax.scatter(placed["day"], placed["inventory"],
               color="#FF6B6B", s=60, zorder=5, label="Order placed")

    # Order received markers
    received = df[df["order_received"]]
    ax.scatter(received["day"], received["inventory"],
               color="#FFD93D", s=60, marker="^", zorder=5, label="Order received")

    # Fill under the line
    ax.fill_between(df["day"], df["inventory"], alpha=0.08, color="#00D4AA")

    # Styling
    ax.set_title(f"Inventory Simulation — {product_name}",
                 color="white", fontsize=14, fontweight="bold", pad=14)
    ax.set_xlabel("Day", color="#AAAAAA", fontsize=11)
    ax.set_ylabel("Units in stock", color="#AAAAAA", fontsize=11)
    ax.tick_params(colors="#AAAAAA", labelsize=9)
    for spine in ax.spines.values():
        spine.set_edgecolor("#333333")

    ax.legend(facecolor="#1A1D27", edgecolor="#333333",
              labelcolor="white", fontsize=9, loc="upper right")
    ax.grid(color="#222222", linewidth=0.5, linestyle="-")

    plt.tight_layout()

    if output_path:
        os.makedirs(os.path.dirname(output_path) if os.path.dirname(output_path) else ".", exist_ok=True)
        plt.savefig(output_path, dpi=150, bbox_inches="tight",
                    facecolor=fig.get_facecolor())
        print(f"  Chart saved → {output_path}")

    plt.show()
    plt.close()


# ─────────────────────────────────────────────
#  CHART: EOQ COST CURVES
# ─────────────────────────────────────────────

def plot_cost_curves(annual_demand: float, ordering_cost: float,
                     holding_cost: float, eoq: float, output_path: str = None):
    """Plots ordering, holding and total cost curves vs. order quantity."""

    quantities = np.linspace(1, eoq * 3, 300)
    order_costs  = (annual_demand / quantities) * ordering_cost
    holding_costs = (quantities / 2) * holding_cost
    total_costs   = order_costs + holding_costs

    fig, ax = plt.subplots(figsize=(10, 5))
    fig.patch.set_facecolor("#0F1117")
    ax.set_facecolor("#0F1117")

    ax.plot(quantities, order_costs,  color="#FF6B6B", linewidth=1.8, label="Ordering cost")
    ax.plot(quantities, holding_costs, color="#4ECDC4", linewidth=1.8, label="Holding cost")
    ax.plot(quantities, total_costs,   color="#FFD93D", linewidth=2.2, label="Total cost")

    # EOQ optimal point
    min_cost = total_costs.min()
    ax.axvline(x=eoq, color="white", linewidth=1, linestyle="--", alpha=0.6)
    ax.scatter([eoq], [min_cost], color="white", s=80, zorder=6)
    ax.annotate(f"EOQ = {eoq:.0f}", xy=(eoq, min_cost),
                xytext=(eoq + eoq * 0.08, min_cost * 1.05),
                color="white", fontsize=9,
                arrowprops=dict(arrowstyle="->", color="white", lw=0.8))

    ax.set_title("EOQ — Cost Curves", color="white", fontsize=13, fontweight="bold", pad=12)
    ax.set_xlabel("Order Quantity (units)", color="#AAAAAA", fontsize=10)
    ax.set_ylabel("Annual Cost ($)", color="#AAAAAA", fontsize=10)
    ax.tick_params(colors="#AAAAAA", labelsize=9)
    for spine in ax.spines.values():
        spine.set_edgecolor("#333333")

    ax.legend(facecolor="#1A1D27", edgecolor="#333333", labelcolor="white", fontsize=9)
    ax.grid(color="#222222", linewidth=0.5)

    plt.tight_layout()

    if output_path:
        os.makedirs(os.path.dirname(output_path) if os.path.dirname(output_path) else ".", exist_ok=True)
        plt.savefig(output_path, dpi=150, bbox_inches="tight",
                    facecolor=fig.get_facecolor())
        print(f"  Chart saved → {output_path}")

    plt.show()
    plt.close()


# ─────────────────────────────────────────────
#  SUMMARY REPORT
# ─────────────────────────────────────────────

def print_report(params: dict, eoq: float, rop: float,
                 safety_stock: float, annual_cost: float):
    """Prints a formatted summary to the console."""

    sep = "─" * 52
    print(f"\n{'═' * 52}")
    print(f"  INVENTORY SIMULATOR — EOQ / ROP REPORT")
    print(f"{'═' * 52}")
    print(f"  Product          : {params['product_name']}")
    print(sep)
    print(f"  INPUT PARAMETERS")
    print(sep)
    print(f"  Annual demand    : {params['annual_demand']:>10,.0f}  units/year")
    print(f"  Daily demand     : {params['annual_demand']/365:>10.1f}  units/day")
    print(f"  Ordering cost    : ${params['ordering_cost']:>9,.2f}  per order")
    print(f"  Holding cost     : ${params['holding_cost']:>9,.2f}  per unit/year")
    print(f"  Lead time        : {params['lead_time_days']:>10}  days")
    print(f"  Service level    : {params['service_level']*100:>9.0f}%")
    print(sep)
    print(f"  RESULTS")
    print(sep)
    print(f"  EOQ              : {eoq:>10,.0f}  units per order")
    print(f"  Safety stock     : {safety_stock:>10,.0f}  units")
    print(f"  ROP              : {rop:>10,.0f}  units (reorder point)")
    print(f"  Orders/year      : {params['annual_demand']/eoq:>10.1f}  times")
    print(f"  Cycle time       : {eoq/params['annual_demand']*365:>10.1f}  days between orders")
    print(f"  Annual cost      : ${annual_cost:>9,.2f}")
    print(f"{'═' * 52}\n")


# ─────────────────────────────────────────────
#  EXPORT TO CSV
# ─────────────────────────────────────────────

def export_results(df: pd.DataFrame, params: dict, eoq: float,
                   rop: float, safety_stock: float, annual_cost: float):
    """Exports simulation data and KPI summary to CSV files."""

    os.makedirs("data", exist_ok=True)

    # Simulation data
    sim_path = "data/simulation_results.csv"
    df.to_csv(sim_path, index=False)
    print(f"  Simulation data  → {sim_path}")

    # KPI summary
    kpi = {
        "product":        [params["product_name"]],
        "annual_demand":  [params["annual_demand"]],
        "eoq":            [eoq],
        "rop":            [rop],
        "safety_stock":   [safety_stock],
        "annual_cost":    [annual_cost],
        "orders_per_year":[round(params["annual_demand"] / eoq, 1)],
        "cycle_days":     [round(eoq / params["annual_demand"] * 365, 1)],
    }
    kpi_path = "data/kpi_summary.csv"
    pd.DataFrame(kpi).to_csv(kpi_path, index=False)
    print(f"  KPI summary      → {kpi_path}")


# ─────────────────────────────────────────────
#  MAIN — EXAMPLE SCENARIO
# ─────────────────────────────────────────────

if __name__ == "__main__":

    # ── Define your product parameters here ──────────────────
    params = {
        "product_name":  "Industrial Filter A-200",
        "annual_demand":  4800,    # units per year
        "ordering_cost":  150.0,   # $ per order
        "holding_cost":   12.0,    # $ per unit per year
        "lead_time_days": 7,       # days
        "service_level":  0.95,    # 95% → z = 1.65
        "std_daily_demand": 3.5,   # standard deviation of daily demand
    }
    # ─────────────────────────────────────────────────────────

    Z_95 = 1.65  # z-score for 95% service level

    # 1. Calculate EOQ
    eoq = calculate_eoq(
        annual_demand  = params["annual_demand"],
        ordering_cost  = params["ordering_cost"],
        holding_cost   = params["holding_cost"],
    )

    # 2. Calculate safety stock
    safety_stock = calculate_safety_stock(
        z_score        = Z_95,
        std_demand     = params["std_daily_demand"],
        lead_time_days = params["lead_time_days"],
    )

    # 3. Calculate ROP
    daily_demand = params["annual_demand"] / 365
    rop = calculate_rop(
        daily_demand   = daily_demand,
        lead_time_days = params["lead_time_days"],
        safety_stock   = safety_stock,
    )

    # 4. Total annual cost
    annual_cost = total_annual_cost(
        annual_demand  = params["annual_demand"],
        eoq            = eoq,
        ordering_cost  = params["ordering_cost"],
        holding_cost   = params["holding_cost"],
    )

    # 5. Print report
    print_report(params, eoq, rop, safety_stock, annual_cost)

    # 6. Run simulation (365 days)
    print("  Running simulation...")
    df_sim = simulate_inventory(
        annual_demand  = params["annual_demand"],
        eoq            = eoq,
        rop            = rop,
        lead_time_days = params["lead_time_days"],
        days           = 365,
    )

    # 7. Export data
    print("\n  Exporting data...")
    export_results(df_sim, params, eoq, rop, safety_stock, annual_cost)

    # 8. Generate charts
    print("\n  Generating charts...")
    plot_inventory(
        df           = df_sim,
        rop          = rop,
        eoq          = eoq,
        product_name = params["product_name"],
        output_path  = "data/inventory_level.png",
    )

    plot_cost_curves(
        annual_demand  = params["annual_demand"],
        ordering_cost  = params["ordering_cost"],
        holding_cost   = params["holding_cost"],
        eoq            = eoq,
        output_path    = "data/cost_curves.png",
    )

    print("\n  Done. Check the /data folder for outputs.\n")