simulation.py

import numpy as np
import pandas as pd
from typing import Dict, List, Tuple, Any

class D2CSimulator:
    def __init__(self, config: Dict[str, Any]):
        """
        config contains:
          - materials: Dict[str, Dict] containing 'CurrentStock', 'CostPerUnit', 'LeadTimeDays', 'ShelfLifeDays', 'UnitsPerProduct', 'Category'
          - financial: Dict[str, float] containing 'StartingCash', 'MRP', 'DistributorMarginPercent', 'LogisticsCostPerUnit', 'FixedMonthlyOverheads', 'DailyMarketingSpend', 'ROAS', 'AverageOrderValue'
          - edge_cases: Dict[str, bool] containing 'viral_spike', 'monsoon_block', 'inflation_shock'
        """
        self.config = config
        self.materials = config['materials']
        self.financial = config['financial']
        self.edge_cases = config.get('edge_cases', {
            'viral_spike': False,
            'monsoon_block': False,
            'inflation_shock': False
        })
        
    def run_simulation(self, random_seed: int = None) -> Dict[str, Any]:
        """
        Runs a 90-day deterministic simulation with small daily demand variance.
        """
        if random_seed is not None:
            np.random.seed(random_seed)
            
        days = 90
        
        # Financial variables
        cash = self.financial['StartingCash']
        starting_cash = cash
        mrp = self.financial['MRP']
        margin = self.financial['DistributorMarginPercent'] / 100.0
        asp = mrp * (1 - margin)  # Average Selling Price received by brand
        logistics_cost_per_unit = self.financial['LogisticsCostPerUnit']
        daily_fixed_overhead = self.financial['FixedMonthlyOverheads'] / 30.0
        daily_mkt_spend = self.financial['DailyMarketingSpend']
        roas = self.financial['ROAS']
        aov = self.financial['AverageOrderValue']
        
        # Marketing-driven daily base demand:
        # Base units sold online = (Daily Marketing Spend * ROAS) / Average Order Value
        # Plus some organic base demand
        mkt_demand = (daily_mkt_spend * roas) / aov if aov > 0 else 0
        organic_demand = 50.0  # constant organic demand
        base_daily_demand = mkt_demand + organic_demand
        
        # Materials initialization
        # Inventory is tracked as batches: list of dicts: {'qty': float, 'expiry_day': int, 'unit_cost': float}
        # Packaging is non-perishable, but tracked similarly for consistency (expiry_day = 9999)
        inventory_batches: Dict[str, List[Dict[str, Any]]] = {}
        for name, m in self.materials.items():
            cost = m['CostPerUnit']
            qty = m['CurrentStock']
            shelf_life = m['ShelfLifeDays']
            expiry = 9999 if shelf_life >= 9999 else shelf_life
            inventory_batches[name] = [{'qty': qty, 'expiry_day': expiry, 'unit_cost': cost}]
            
        # Track orders in transit: list of dicts: {'arrival_day': int, 'qty': float, 'cost_per_unit': float, 'item': str}
        orders_in_transit: List[Dict[str, Any]] = []
        
        # Metrics logs
        daily_logs = []
        
        cumulative_revenue = 0.0
        cumulative_cogs = 0.0
        cumulative_mkt_spend = 0.0
        cumulative_overhead = 0.0
        cumulative_logistics = 0.0
        cumulative_purchases = 0.0
        cumulative_spoilage_loss = 0.0
        cumulative_stockouts = 0
        cumulative_stockout_penalty = 0.0
        
        for t in range(1, days + 1):
            # 1. Handle incoming orders
            arrived_orders = [o for o in orders_in_transit if o['arrival_day'] == t]
            orders_in_transit = [o for o in orders_in_transit if o['arrival_day'] != t]
            
            for o in arrived_orders:
                name = o['item']
                shelf_life = self.materials[name]['ShelfLifeDays']
                expiry = t + shelf_life if shelf_life < 9999 else 9999
                inventory_batches[name].append({
                    'qty': o['qty'],
                    'expiry_day': expiry,
                    'unit_cost': o['cost_per_unit']
                })
                
            # 2. Spoilage/Expiration check (at start of day)
            spoilage_today = {}
            for name, batches in inventory_batches.items():
                expired_qty = 0.0
                expired_value = 0.0
                active_batches = []
                for b in batches:
                    if b['expiry_day'] < t:
                        expired_qty += b['qty']
                        expired_value += b['qty'] * b['unit_cost']
                    else:
                        active_batches.append(b)
                inventory_batches[name] = active_batches
                if expired_qty > 0:
                    spoilage_today[name] = {
                        'qty': expired_qty,
                        'value': expired_value
                    }
                    cumulative_spoilage_loss += expired_value

            # 3. Demand calculation
            # Add some randomness (+/- 15%)
            variation = np.random.uniform(-0.15, 0.15)
            day_demand = base_daily_demand * (1 + variation)
            
            # Edge Case A: Viral Demand Spike (Day 15 to 21)
            is_spike = False
            if self.edge_cases['viral_spike'] and (15 <= t <= 21):
                day_demand *= 4.0
                is_spike = True
                
            day_demand = max(0.0, float(day_demand))
            
            # 4. Production & Inventory Consumption
            # How many units of product can we actually fulfill given current inventory?
            # We look at the bottleneck across all ingredients
            possible_production = day_demand
            for name, m in self.materials.items():
                units_needed_per_product = m['UnitsPerProduct']
                if units_needed_per_product > 0:
                    total_stock = sum(b['qty'] for b in inventory_batches[name])
                    possible_for_this = total_stock / units_needed_per_product
                    if possible_for_this < possible_production:
                        possible_production = possible_for_this
            
            actual_sales = possible_production
            stockout_units = max(0.0, day_demand - actual_sales)
            cumulative_stockouts += stockout_units
            
            # Penalty for Q-Commerce stockouts (10% of ASP for every stocked out unit)
            stockout_penalty = stockout_units * asp * 0.10
            cumulative_stockout_penalty += stockout_penalty
            
            # Consume inventory batches (FIFO)
            cogs_today = 0.0
            for name, m in self.materials.items():
                units_needed = actual_sales * m['UnitsPerProduct']
                consumed_qty = 0.0
                while consumed_qty < units_needed and inventory_batches[name]:
                    b = inventory_batches[name][0]
                    needed_from_batch = units_needed - consumed_qty
                    if b['qty'] <= needed_from_batch:
                        # Consume entire batch
                        consumed_qty += b['qty']
                        cogs_today += b['qty'] * b['unit_cost']
                        inventory_batches[name].pop(0)
                    else:
                        # Consume part of batch
                        cogs_today += needed_from_batch * b['unit_cost']
                        b['qty'] -= needed_from_batch
                        consumed_qty += needed_from_batch
            
            cumulative_cogs += cogs_today
            
            # 5. Financial Calculations
            revenue_today = actual_sales * asp
            cumulative_revenue += revenue_today
            
            logistics_today = actual_sales * logistics_cost_per_unit
            cumulative_logistics += logistics_today
            
            cumulative_mkt_spend += daily_mkt_spend
            cumulative_overhead += daily_fixed_overhead
            
            # Daily net non-procurement cash movement
            # (Procurement payments happen immediately upon order placement below)
            cash_outflows = daily_mkt_spend + daily_fixed_overhead + logistics_today + stockout_penalty
            cash_inflows = revenue_today
            
            cash = cash + cash_inflows - cash_outflows
            
            # 6. Reordering Decisions (End of day)
            # ROP = LeadTimeDays * DailyUsageRate * SafetyFactor
            # DailyUsageRate is based on base_daily_demand (average demand)
            reorder_payments_today = 0.0
            for name, m in self.materials.items():
                units_needed_per_product = m['UnitsPerProduct']
                if units_needed_per_product == 0:
                    continue
                    
                daily_usage_base = base_daily_demand * units_needed_per_product
                lead_time = m['LeadTimeDays']
                
                # Check for Edge Case B: Monsoon logistics delay
                # If monsoon_block is active, packaging materials have a 15-day delay for any order arriving between day 30 and 45.
                # Since lead time is variable under monsoon, we account for it.
                current_lead_time = lead_time
                if self.edge_cases['monsoon_block'] and m['Category'].lower() == 'packaging':
                    # If this order would arrive in the monsoon window [30, 45]
                    est_arrival = t + lead_time
                    if 30 <= est_arrival <= 45:
                        current_lead_time += 15
                
                rop = daily_usage_base * lead_time * 1.5  # Safety stock factor of 1.5
                
                # Current stock on hand + in transit
                on_hand = sum(b['qty'] for b in inventory_batches[name])
                in_transit = sum(o['qty'] for o in orders_in_transit if o['item'] == name)
                total_position = on_hand + in_transit
                
                if total_position < rop:
                    # Place order for 15-day supply
                    order_qty = daily_usage_base * 15
                    
                    # Edge Case C: Inflationary Shock (Cost of whey/protein jumps on Day 45)
                    unit_cost = m['CostPerUnit']
                    if self.edge_cases['inflation_shock'] and t >= 45:
                        # Whey/proteins or primary ingredients cost 25% more
                        if "whey" in name.lower() or "protein" in name.lower() or name.lower() == "organic whey protein":
                            unit_cost *= 1.25
                            
                    order_cost = order_qty * unit_cost
                    
                    # Cash outflow happens immediately
                    cash -= order_cost
                    reorder_payments_today += order_cost
                    cumulative_purchases += order_cost
                    
                    orders_in_transit.append({
                        'arrival_day': t + current_lead_time,
                        'qty': order_qty,
                        'cost_per_unit': unit_cost,
                        'item': name
                    })
            
            # Save daily status
            stocks = {name: sum(b['qty'] for b in inventory_batches[name]) for name in self.materials.keys()}
            daily_logs.append({
                'day': t,
                'cash': float(cash),
                'revenue': float(revenue_today),
                'cogs': float(cogs_today),
                'spoilage_loss': sum(v['value'] for v in spoilage_today.values()),
                'stockouts': float(stockout_units),
                'stockout_penalty': float(stockout_penalty),
                'reorder_payments': float(reorder_payments_today),
                'is_spike': is_spike,
                'inventory': stocks
            })
            
        df_logs = pd.DataFrame(daily_logs)
        
        # Calculate summary metrics
        net_profit = cumulative_revenue - cumulative_cogs - cumulative_mkt_spend - cumulative_overhead - cumulative_logistics - cumulative_stockout_penalty - cumulative_spoilage_loss
        gross_profit = cumulative_revenue - cumulative_cogs
        gross_margin = (gross_profit / cumulative_revenue) * 100.0 if cumulative_revenue > 0 else 0.0
        
        # Calculate when cash runs out (runway)
        runway_days = days
        runout_idx = df_logs[df_logs['cash'] < 0].index
        if len(runout_idx) > 0:
            runway_days = int(df_logs.loc[runout_idx[0], 'day'])
            
        summary = {
            'starting_cash': starting_cash,
            'ending_cash': float(cash),
            'cumulative_revenue': float(cumulative_revenue),
            'cumulative_cogs': float(cumulative_cogs),
            'cumulative_mkt_spend': float(cumulative_mkt_spend),
            'cumulative_overhead': float(cumulative_overhead),
            'cumulative_logistics': float(cumulative_logistics),
            'cumulative_spoilage_loss': float(cumulative_spoilage_loss),
            'cumulative_stockouts': float(cumulative_stockouts),
            'cumulative_stockout_penalty': float(cumulative_stockout_penalty),
            'net_profit': float(net_profit),
            'gross_margin': float(gross_margin),
            'runway_days': runway_days
        }
        
        return {
            'logs': df_logs,
            'summary': summary
        }

    def run_monte_carlo(self, n_simulations: int = 100) -> Dict[str, Any]:
        """
        Runs n_simulations to capture probabilistic cash runway and inventory positions.
        """
        runways = []
        ending_cashes = []
        daily_cash_paths = np.zeros((90, n_simulations))
        stockout_rates = []
        spoilage_losses = []
        
        for i in range(n_simulations):
            res = self.run_simulation(random_seed=i)
            runways.append(res['summary']['runway_days'])
            ending_cashes.append(res['summary']['ending_cash'])
            daily_cash_paths[:, i] = res['logs']['cash'].values
            
            # Stockout rate (percentage of demand missed)
            # Total demand = cumulative_revenue / average ASP + cumulative_stockouts
            sales_units = res['summary']['cumulative_revenue'] / (self.financial['MRP'] * (1 - self.financial['DistributorMarginPercent']/100.0))
            tot_demand = sales_units + res['summary']['cumulative_stockouts']
            stockout_rate = (res['summary']['cumulative_stockouts'] / tot_demand * 100) if tot_demand > 0 else 0
            stockout_rates.append(stockout_rate)
            
            spoilage_losses.append(res['summary']['cumulative_spoilage_loss'])
            
        # Compute percentiles for Cash Runway over time
        # 10th percentile (Worst Case), 50th percentile (Expected Case), 90th percentile (Best Case)
        worst_case = np.percentile(daily_cash_paths, 10, axis=1)
        expected_case = np.percentile(daily_cash_paths, 50, axis=1)
        best_case = np.percentile(daily_cash_paths, 90, axis=1)
        
        # Calculate runway under Worst, Expected, Best
        def get_runway_from_path(path):
            days_under_zero = np.where(path < 0)[0]
            if len(days_under_zero) > 0:
                return int(days_under_zero[0] + 1)
            return 90
            
        return {
            'worst_cash_path': worst_case.tolist(),
            'expected_cash_path': expected_case.tolist(),
            'best_cash_path': best_case.tolist(),
            'worst_runway': get_runway_from_path(worst_case),
            'expected_runway': get_runway_from_path(expected_case),
            'best_runway': get_runway_from_path(best_case),
            'avg_stockout_rate': float(np.mean(stockout_rates)),
            'avg_spoilage_loss': float(np.mean(spoilage_losses)),
            'cash_runs': daily_cash_paths.tolist()
        }