Skip to content

Supply Chain

deepbullwhip.chain.config.EchelonConfig dataclass

Configuration for a single supply chain echelon.

Source code in deepbullwhip/chain/config.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
@dataclass
class EchelonConfig:
    """Configuration for a single supply chain echelon."""

    name: str
    lead_time: int
    holding_cost: float
    backorder_cost: float
    depreciation_rate: float = 0.0
    service_level: float = 0.95
    initial_inventory: float = 50.0

deepbullwhip.chain.config.default_semiconductor_config()

Return the 4-echelon semiconductor config from the paper.

Weekly depreciation derived from 15% quarterly value loss: weekly_dep = 1 - (1 - 0.15) ** (1/13)

Source code in deepbullwhip/chain/config.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def default_semiconductor_config() -> list[EchelonConfig]:
    """Return the 4-echelon semiconductor config from the paper.

    Weekly depreciation derived from 15% quarterly value loss:
        weekly_dep = 1 - (1 - 0.15) ** (1/13)
    """
    weekly_dep = 1 - (1 - 0.15) ** (1 / 13)
    return [
        EchelonConfig(
            "Distributor", lead_time=2, holding_cost=0.15,
            backorder_cost=0.60, depreciation_rate=weekly_dep,
        ),
        EchelonConfig(
            "OSAT", lead_time=4, holding_cost=0.12,
            backorder_cost=0.50, depreciation_rate=weekly_dep * 0.8,
        ),
        EchelonConfig(
            "Foundry", lead_time=12, holding_cost=0.08,
            backorder_cost=0.40, depreciation_rate=weekly_dep * 0.5,
        ),
        EchelonConfig(
            "Supplier", lead_time=8, holding_cost=0.05,
            backorder_cost=0.30, depreciation_rate=weekly_dep * 0.3,
        ),
    ]

deepbullwhip.chain.echelon.SupplyChainEchelon

Single echelon in a serial supply chain.

Parameters:

Name Type Description Default
name str

Human-readable name (e.g. "Distributor").

 required
lead_time int

Replenishment lead time in periods.

 required
policy OrderingPolicy

The ordering policy (must already know its own lead_time).

 required
cost_fn CostFunction

The per-period cost function.

 required
initial_inventory float

Starting on-hand inventory.

 50.0
Source code in deepbullwhip/chain/echelon.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
class SupplyChainEchelon:
    """Single echelon in a serial supply chain.

    Parameters
    ----------
    name : str
        Human-readable name (e.g. "Distributor").
    lead_time : int
        Replenishment lead time in periods.
    policy : OrderingPolicy
        The ordering policy (must already know its own lead_time).
    cost_fn : CostFunction
        The per-period cost function.
    initial_inventory : float
        Starting on-hand inventory.
    """

    def __init__(
        self,
        name: str,
        lead_time: int,
        policy: OrderingPolicy,
        cost_fn: CostFunction,
        initial_inventory: float = 50.0,
    ) -> None:
        self.name = name
        self.lead_time = lead_time
        self.policy = policy
        self.cost_fn = cost_fn
        self.initial_inventory = initial_inventory

        self.inventory: float = 0.0
        self.pipeline: list[float] = []
        self.orders: list[float] = []
        self.inventory_levels: list[float] = []
        self.costs: list[float] = []

    def reset(self) -> None:
        """Reset state to initial conditions."""
        self.inventory = self.initial_inventory
        self.pipeline = [0.0] * self.lead_time
        self.orders = []
        self.inventory_levels = []
        self.costs = []

    def step(self, demand: float, forecast_mean: float, forecast_std: float) -> float:
        """Execute one period: receive, order, satisfy demand, compute cost.

        Returns the order quantity placed this period.
        """
        # Receive oldest pipeline order
        if self.pipeline:
            self.inventory += self.pipeline.pop(0)

        # Inventory position = on-hand + pipeline
        ip = self.inventory + sum(self.pipeline)

        # Ordering decision (delegated to policy)
        order = self.policy.compute_order(ip, forecast_mean, forecast_std)
        self.orders.append(order)
        self.pipeline.append(order)

        # Satisfy demand
        self.inventory -= demand

        # Cost (delegated to cost function)
        cost = self.cost_fn.compute(self.inventory)
        self.inventory_levels.append(self.inventory)
        self.costs.append(cost)

        return order

    @property
    def orders_array(self) -> TimeSeries:
        return np.asarray(self.orders, dtype=np.float64)

    @property
    def inventory_array(self) -> TimeSeries:
        return np.asarray(self.inventory_levels, dtype=np.float64)

    @property
    def costs_array(self) -> TimeSeries:
        return np.asarray(self.costs, dtype=np.float64)

reset()

Reset state to initial conditions.

Source code in deepbullwhip/chain/echelon.py 
47
48
49
50
51
52
53
def reset(self) -> None:
    """Reset state to initial conditions."""
    self.inventory = self.initial_inventory
    self.pipeline = [0.0] * self.lead_time
    self.orders = []
    self.inventory_levels = []
    self.costs = []

step(demand, forecast_mean, forecast_std)

Execute one period: receive, order, satisfy demand, compute cost.

Returns the order quantity placed this period.

Source code in deepbullwhip/chain/echelon.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def step(self, demand: float, forecast_mean: float, forecast_std: float) -> float:
    """Execute one period: receive, order, satisfy demand, compute cost.

    Returns the order quantity placed this period.
    """
    # Receive oldest pipeline order
    if self.pipeline:
        self.inventory += self.pipeline.pop(0)

    # Inventory position = on-hand + pipeline
    ip = self.inventory + sum(self.pipeline)

    # Ordering decision (delegated to policy)
    order = self.policy.compute_order(ip, forecast_mean, forecast_std)
    self.orders.append(order)
    self.pipeline.append(order)

    # Satisfy demand
    self.inventory -= demand

    # Cost (delegated to cost function)
    cost = self.cost_fn.compute(self.inventory)
    self.inventory_levels.append(self.inventory)
    self.costs.append(cost)

    return order

deepbullwhip.chain.serial.SerialSupplyChain

K-echelon serial supply chain.

Parameters:

Name Type Description Default
echelons list[SupplyChainEchelon] or None

Pre-built echelons. If None, uses default_semiconductor_config().

 None
Source code in deepbullwhip/chain/serial.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
class SerialSupplyChain:
    """K-echelon serial supply chain.

    Parameters
    ----------
    echelons : list[SupplyChainEchelon] or None
        Pre-built echelons. If None, uses default_semiconductor_config().
    """

    def __init__(self, echelons: list[SupplyChainEchelon] | None = None) -> None:
        if echelons is None:
            echelons = self._from_config(default_semiconductor_config())
        self.echelons = echelons
        self.K = len(echelons)

    @staticmethod
    def _from_config(configs: list[EchelonConfig]) -> list[SupplyChainEchelon]:
        result = []
        for cfg in configs:
            total_h = cfg.holding_cost + cfg.depreciation_rate
            policy = OrderUpToPolicy(
                lead_time=cfg.lead_time, service_level=cfg.service_level
            )
            cost_fn = NewsvendorCost(
                holding_cost=total_h, backorder_cost=cfg.backorder_cost
            )
            result.append(
                SupplyChainEchelon(
                    name=cfg.name,
                    lead_time=cfg.lead_time,
                    policy=policy,
                    cost_fn=cost_fn,
                    initial_inventory=cfg.initial_inventory,
                )
            )
        return result

    @classmethod
    def from_config(cls, configs: list[EchelonConfig]) -> SerialSupplyChain:
        """Build a chain from a list of EchelonConfig objects."""
        return cls(echelons=cls._from_config(configs))

    def reset(self) -> None:
        for e in self.echelons:
            e.reset()

    def simulate(
        self,
        demand: TimeSeries,
        forecasts_mean: TimeSeries,
        forecasts_std: TimeSeries,
    ) -> SimulationResult:
        """Run the full simulation over T periods.

        Parameters
        ----------
        demand : array, shape (T,)
        forecasts_mean : array, shape (T,)
            Point forecasts for echelon 1 (end-customer demand).
        forecasts_std : array, shape (T,)
            Forecast error std estimates for echelon 1.

        Returns
        -------
        SimulationResult
        """
        T = len(demand)
        self.reset()

        for t in range(T):
            d = demand[t]
            fm = forecasts_mean[t]
            fs = forecasts_std[t]

            # Echelon 1 faces end-customer demand
            order = self.echelons[0].step(d, fm, fs)

            # Each subsequent echelon faces orders from the previous
            for k in range(1, self.K):
                upstream_demand = order
                if t > 0:
                    recent = self.echelons[k - 1].orders[
                        -min(8, len(self.echelons[k - 1].orders)) :
                    ]
                    upstream_fm = float(np.mean(recent))
                    upstream_fs = (
                        float(np.std(recent)) if len(recent) > 1 else fs
                    )
                else:
                    upstream_fm = fm
                    upstream_fs = fs
                order = self.echelons[k].step(upstream_demand, upstream_fm, upstream_fs)

        return self._compute_results(demand)

    def _compute_results(self, demand: TimeSeries) -> SimulationResult:
        var_demand = float(np.var(demand))
        echelon_results: list[EchelonResult] = []

        for k, e in enumerate(self.echelons):
            orders = e.orders_array
            inv = e.inventory_array
            costs = e.costs_array

            if k == 0:
                bw = float(np.var(orders) / var_demand) if var_demand > 0 else 1.0
            else:
                var_prev = float(np.var(self.echelons[k - 1].orders_array))
                bw = float(np.var(orders) / var_prev) if var_prev > 0 else 1.0

            echelon_results.append(
                EchelonResult(
                    name=e.name,
                    orders=orders,
                    inventory_levels=inv,
                    costs=costs,
                    bullwhip_ratio=bw,
                    fill_rate=float(np.mean(inv >= 0)),
                    total_cost=float(costs.sum()),
                )
            )

        all_orders_K = self.echelons[-1].orders_array
        cum_bw = float(np.var(all_orders_K) / var_demand) if var_demand > 0 else 1.0
        total_cost = sum(er.total_cost for er in echelon_results)

        return SimulationResult(
            echelon_results=echelon_results,
            cumulative_bullwhip=cum_bw,
            total_cost=total_cost,
        )

from_config(configs) classmethod

Build a chain from a list of EchelonConfig objects.

Source code in deepbullwhip/chain/serial.py 
49
50
51
52
@classmethod
def from_config(cls, configs: list[EchelonConfig]) -> SerialSupplyChain:
    """Build a chain from a list of EchelonConfig objects."""
    return cls(echelons=cls._from_config(configs))

simulate(demand, forecasts_mean, forecasts_std)

Run the full simulation over T periods.

Parameters:

Name Type Description Default
demand (array, shape(T))
required
forecasts_mean (array, shape(T))

Point forecasts for echelon 1 (end-customer demand).

 required
forecasts_std (array, shape(T))

Forecast error std estimates for echelon 1.

 required

Returns:

Type Description
SimulationResult
Source code in deepbullwhip/chain/serial.py 
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def simulate(
    self,
    demand: TimeSeries,
    forecasts_mean: TimeSeries,
    forecasts_std: TimeSeries,
) -> SimulationResult:
    """Run the full simulation over T periods.

    Parameters
    ----------
    demand : array, shape (T,)
    forecasts_mean : array, shape (T,)
        Point forecasts for echelon 1 (end-customer demand).
    forecasts_std : array, shape (T,)
        Forecast error std estimates for echelon 1.

    Returns
    -------
    SimulationResult
    """
    T = len(demand)
    self.reset()

    for t in range(T):
        d = demand[t]
        fm = forecasts_mean[t]
        fs = forecasts_std[t]

        # Echelon 1 faces end-customer demand
        order = self.echelons[0].step(d, fm, fs)

        # Each subsequent echelon faces orders from the previous
        for k in range(1, self.K):
            upstream_demand = order
            if t > 0:
                recent = self.echelons[k - 1].orders[
                    -min(8, len(self.echelons[k - 1].orders)) :
                ]
                upstream_fm = float(np.mean(recent))
                upstream_fs = (
                    float(np.std(recent)) if len(recent) > 1 else fs
                )
            else:
                upstream_fm = fm
                upstream_fs = fs
            order = self.echelons[k].step(upstream_demand, upstream_fm, upstream_fs)

    return self._compute_results(demand)

deepbullwhip.chain.vectorized.VectorizedSupplyChain

Matrix-based supply chain simulation for Monte Carlo batching.

Instead of Python lists and per-element operations, this engine pre-allocates (N, K, T) arrays and uses NumPy broadcasting for all N paths simultaneously. The pipeline is implemented as a circular buffer with O(1) index arithmetic instead of O(L) list.pop(0).

Parameters:

Name Type Description Default
configs list[EchelonConfig] or None

Echelon configurations. If None, uses default semiconductor config.

 None
Source code in deepbullwhip/chain/vectorized.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
class VectorizedSupplyChain:
    """Matrix-based supply chain simulation for Monte Carlo batching.

    Instead of Python lists and per-element operations, this engine
    pre-allocates (N, K, T) arrays and uses NumPy broadcasting for all
    N paths simultaneously. The pipeline is implemented as a circular
    buffer with O(1) index arithmetic instead of O(L) list.pop(0).

    Parameters
    ----------
    configs : list[EchelonConfig] or None
        Echelon configurations. If None, uses default semiconductor config.
    """

    def __init__(self, configs: list[EchelonConfig] | None = None) -> None:
        if configs is None:
            configs = default_semiconductor_config()
        self.configs = configs
        self.K = len(configs)

        # Pre-compute per-echelon parameters as arrays for vectorized access
        self.lead_times = np.array([c.lead_time for c in configs])
        self.L_max = int(self.lead_times.max())

        self.h = np.array([c.holding_cost + c.depreciation_rate for c in configs])
        self.b = np.array([c.backorder_cost for c in configs])
        self.z_alpha = np.array([stats.norm.ppf(c.service_level) for c in configs])
        self.initial_inv = np.array([c.initial_inventory for c in configs])

    def simulate(
        self,
        demand: np.ndarray,
        forecasts_mean: np.ndarray,
        forecasts_std: np.ndarray,
    ) -> BatchSimulationResult:
        """Run vectorized simulation.

        Parameters
        ----------
        demand : array, shape (N, T) or (T,)
            If 1-D, broadcast to all N paths (N=1).
        forecasts_mean : array, shape (N, T) or (T,)
        forecasts_std : array, shape (N, T) or (T,)

        Returns
        -------
        BatchSimulationResult
        """
        # Normalize to 2-D
        if demand.ndim == 1:
            demand = demand[np.newaxis, :]
            forecasts_mean = forecasts_mean[np.newaxis, :]
            forecasts_std = forecasts_std[np.newaxis, :]

        N, T = demand.shape
        K = self.K

        # ── Pre-allocate all matrices ────────────────────────────────
        orders = np.zeros((N, K, T))
        inventory = np.zeros((N, K, T))
        costs = np.zeros((N, K, T))

        # Pipeline: circular buffer (N, K, L_max)
        pipeline = np.zeros((N, K, self.L_max))
        # On-hand inventory: (N, K)
        inv = np.tile(self.initial_inv, (N, 1))  # (N, K)

        # Pipeline read pointer per echelon (wraps with %)
        ptr = np.zeros(K, dtype=int)

        # Lead time masks for circular buffer indexing
        # lead_time_mask[k, l] = True if l < lead_times[k]
        lead_time_mask = np.arange(self.L_max)[None, :] < self.lead_times[:, None]  # (K, L_max)

        # ── Time loop (sequential — state dependency) ────────────────
        for t in range(T):
            # --- Receive from pipeline: oldest order arrives ---
            # For each echelon k, read pipeline[n, k, ptr[k]]
            incoming = pipeline[:, np.arange(K), ptr]  # (N, K)
            inv += incoming

            # --- Compute inventory position: on-hand + pipeline sum ---
            # Mask out unused pipeline slots (beyond lead_time)
            pipeline_masked = pipeline * lead_time_mask[None, :, :]  # (N, K, L_max)
            pipeline_sum = pipeline_masked.sum(axis=2)  # (N, K)
            ip = inv + pipeline_sum  # (N, K)

            # --- Compute forecasts for each echelon ---
            # Echelon 0: use provided forecasts
            # Echelons 1..K-1: rolling mean/std of downstream orders (last 8)
            fm = np.zeros((N, K))
            fs = np.zeros((N, K))
            fm[:, 0] = forecasts_mean[:, t]
            fs[:, 0] = forecasts_std[:, t]

            if t > 0:
                window = min(8, t)
                for k in range(1, K):
                    recent = orders[:, k - 1, t - window : t]  # (N, window)
                    fm[:, k] = recent.mean(axis=1)
                    fs[:, k] = recent.std(axis=1) if window > 1 else forecasts_std[:, t]
            else:
                fm[:, 1:] = forecasts_mean[:, t : t + 1]
                fs[:, 1:] = forecasts_std[:, t : t + 1]

            # --- Order-Up-To policy (vectorized across N and K) ---
            L_plus_1 = (self.lead_times + 1).astype(float)  # (K,)
            S = L_plus_1[None, :] * fm + (
                self.z_alpha[None, :] * fs * np.sqrt(L_plus_1)[None, :]
            )  # (N, K)
            order_qty = np.maximum(0.0, S - ip)  # (N, K)

            # --- Write order into pipeline at current pointer ---
            pipeline[:, np.arange(K), ptr] = order_qty
            orders[:, :, t] = order_qty

            # --- Satisfy demand ---
            # Echelon 0 faces customer demand; echelon k faces echelon k-1 orders
            echelon_demand = np.zeros((N, K))
            echelon_demand[:, 0] = demand[:, t]
            if K > 1:
                echelon_demand[:, 1:] = order_qty[:, :-1]

            inv -= echelon_demand

            # --- Compute costs (newsvendor: h * inv+ + b * inv-) ---
            costs_t = np.where(inv >= 0,
                               self.h[None, :] * inv,
                               self.b[None, :] * np.abs(inv))
            costs[:, :, t] = costs_t
            inventory[:, :, t] = inv

            # --- Advance pipeline pointer (circular buffer) ---
            ptr = (ptr + 1) % self.L_max

        # ── Compute metrics (fully vectorized) ───────────────────────
        var_demand = np.var(demand, axis=1)  # (N,)

        # Variance of orders: (N, K)
        var_orders = np.var(orders, axis=2)

        # Bullwhip ratios
        bullwhip_ratios = np.ones((N, K))
        safe_var_demand = np.where(var_demand > 0, var_demand, 1.0)
        bullwhip_ratios[:, 0] = var_orders[:, 0] / safe_var_demand

        for k in range(1, K):
            safe_var_prev = np.where(var_orders[:, k - 1] > 0, var_orders[:, k - 1], 1.0)
            bullwhip_ratios[:, k] = var_orders[:, k] / safe_var_prev

        # Fill rates: (N, K)
        fill_rates = np.mean(inventory >= 0, axis=2)

        # Total costs: (N, K)
        total_costs = costs.sum(axis=2)

        # Cumulative bullwhip: (N,)
        var_last = var_orders[:, -1]
        cumulative_bullwhip = var_last / safe_var_demand

        return BatchSimulationResult(
            orders=orders,
            inventory=inventory,
            costs=costs,
            bullwhip_ratios=bullwhip_ratios,
            fill_rates=fill_rates,
            total_costs=total_costs,
            cumulative_bullwhip=cumulative_bullwhip,
        )

simulate(demand, forecasts_mean, forecasts_std)

Run vectorized simulation.

Parameters:

Name Type Description Default
demand (array, shape(N, T) or (T,))

If 1-D, broadcast to all N paths (N=1).

 required
forecasts_mean (array, shape(N, T) or (T,))
required
forecasts_std (array, shape(N, T) or (T,))
required

Returns:

Type Description
BatchSimulationResult
Source code in deepbullwhip/chain/vectorized.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
def simulate(
    self,
    demand: np.ndarray,
    forecasts_mean: np.ndarray,
    forecasts_std: np.ndarray,
) -> BatchSimulationResult:
    """Run vectorized simulation.

    Parameters
    ----------
    demand : array, shape (N, T) or (T,)
        If 1-D, broadcast to all N paths (N=1).
    forecasts_mean : array, shape (N, T) or (T,)
    forecasts_std : array, shape (N, T) or (T,)

    Returns
    -------
    BatchSimulationResult
    """
    # Normalize to 2-D
    if demand.ndim == 1:
        demand = demand[np.newaxis, :]
        forecasts_mean = forecasts_mean[np.newaxis, :]
        forecasts_std = forecasts_std[np.newaxis, :]

    N, T = demand.shape
    K = self.K

    # ── Pre-allocate all matrices ────────────────────────────────
    orders = np.zeros((N, K, T))
    inventory = np.zeros((N, K, T))
    costs = np.zeros((N, K, T))

    # Pipeline: circular buffer (N, K, L_max)
    pipeline = np.zeros((N, K, self.L_max))
    # On-hand inventory: (N, K)
    inv = np.tile(self.initial_inv, (N, 1))  # (N, K)

    # Pipeline read pointer per echelon (wraps with %)
    ptr = np.zeros(K, dtype=int)

    # Lead time masks for circular buffer indexing
    # lead_time_mask[k, l] = True if l < lead_times[k]
    lead_time_mask = np.arange(self.L_max)[None, :] < self.lead_times[:, None]  # (K, L_max)

    # ── Time loop (sequential — state dependency) ────────────────
    for t in range(T):
        # --- Receive from pipeline: oldest order arrives ---
        # For each echelon k, read pipeline[n, k, ptr[k]]
        incoming = pipeline[:, np.arange(K), ptr]  # (N, K)
        inv += incoming

        # --- Compute inventory position: on-hand + pipeline sum ---
        # Mask out unused pipeline slots (beyond lead_time)
        pipeline_masked = pipeline * lead_time_mask[None, :, :]  # (N, K, L_max)
        pipeline_sum = pipeline_masked.sum(axis=2)  # (N, K)
        ip = inv + pipeline_sum  # (N, K)

        # --- Compute forecasts for each echelon ---
        # Echelon 0: use provided forecasts
        # Echelons 1..K-1: rolling mean/std of downstream orders (last 8)
        fm = np.zeros((N, K))
        fs = np.zeros((N, K))
        fm[:, 0] = forecasts_mean[:, t]
        fs[:, 0] = forecasts_std[:, t]

        if t > 0:
            window = min(8, t)
            for k in range(1, K):
                recent = orders[:, k - 1, t - window : t]  # (N, window)
                fm[:, k] = recent.mean(axis=1)
                fs[:, k] = recent.std(axis=1) if window > 1 else forecasts_std[:, t]
        else:
            fm[:, 1:] = forecasts_mean[:, t : t + 1]
            fs[:, 1:] = forecasts_std[:, t : t + 1]

        # --- Order-Up-To policy (vectorized across N and K) ---
        L_plus_1 = (self.lead_times + 1).astype(float)  # (K,)
        S = L_plus_1[None, :] * fm + (
            self.z_alpha[None, :] * fs * np.sqrt(L_plus_1)[None, :]
        )  # (N, K)
        order_qty = np.maximum(0.0, S - ip)  # (N, K)

        # --- Write order into pipeline at current pointer ---
        pipeline[:, np.arange(K), ptr] = order_qty
        orders[:, :, t] = order_qty

        # --- Satisfy demand ---
        # Echelon 0 faces customer demand; echelon k faces echelon k-1 orders
        echelon_demand = np.zeros((N, K))
        echelon_demand[:, 0] = demand[:, t]
        if K > 1:
            echelon_demand[:, 1:] = order_qty[:, :-1]

        inv -= echelon_demand

        # --- Compute costs (newsvendor: h * inv+ + b * inv-) ---
        costs_t = np.where(inv >= 0,
                           self.h[None, :] * inv,
                           self.b[None, :] * np.abs(inv))
        costs[:, :, t] = costs_t
        inventory[:, :, t] = inv

        # --- Advance pipeline pointer (circular buffer) ---
        ptr = (ptr + 1) % self.L_max

    # ── Compute metrics (fully vectorized) ───────────────────────
    var_demand = np.var(demand, axis=1)  # (N,)

    # Variance of orders: (N, K)
    var_orders = np.var(orders, axis=2)

    # Bullwhip ratios
    bullwhip_ratios = np.ones((N, K))
    safe_var_demand = np.where(var_demand > 0, var_demand, 1.0)
    bullwhip_ratios[:, 0] = var_orders[:, 0] / safe_var_demand

    for k in range(1, K):
        safe_var_prev = np.where(var_orders[:, k - 1] > 0, var_orders[:, k - 1], 1.0)
        bullwhip_ratios[:, k] = var_orders[:, k] / safe_var_prev

    # Fill rates: (N, K)
    fill_rates = np.mean(inventory >= 0, axis=2)

    # Total costs: (N, K)
    total_costs = costs.sum(axis=2)

    # Cumulative bullwhip: (N,)
    var_last = var_orders[:, -1]
    cumulative_bullwhip = var_last / safe_var_demand

    return BatchSimulationResult(
        orders=orders,
        inventory=inventory,
        costs=costs,
        bullwhip_ratios=bullwhip_ratios,
        fill_rates=fill_rates,
        total_costs=total_costs,
        cumulative_bullwhip=cumulative_bullwhip,
    )

deepbullwhip.chain.vectorized.BatchSimulationResult dataclass

Results from a vectorized batch simulation.

All arrays have a leading N (number of paths) dimension.

Source code in deepbullwhip/chain/vectorized.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
@dataclass
class BatchSimulationResult:
    """Results from a vectorized batch simulation.

    All arrays have a leading N (number of paths) dimension.
    """

    orders: np.ndarray          # (N, K, T)
    inventory: np.ndarray       # (N, K, T)
    costs: np.ndarray           # (N, K, T)
    bullwhip_ratios: np.ndarray # (N, K)
    fill_rates: np.ndarray      # (N, K)
    total_costs: np.ndarray     # (N, K)
    cumulative_bullwhip: np.ndarray  # (N,)

    @property
    def n_paths(self) -> int:
        return self.orders.shape[0]

    @property
    def n_echelons(self) -> int:
        return self.orders.shape[1]

    @property
    def n_periods(self) -> int:
        return self.orders.shape[2]

    def mean_metrics(self) -> dict[str, float]:
        """Average metrics across all N paths."""
        d: dict[str, float] = {}
        K = self.n_echelons
        for k in range(K):
            d[f"BW_{k + 1}"] = float(self.bullwhip_ratios[:, k].mean())
            d[f"cost_{k + 1}"] = float(self.total_costs[:, k].mean())
            d[f"fill_rate_{k + 1}"] = float(self.fill_rates[:, k].mean())
        d["BW_cumulative"] = float(self.cumulative_bullwhip.mean())
        d["total_cost"] = float(self.total_costs.sum(axis=1).mean())
        return d

    def to_simulation_result(self, path_index: int = 0) -> SimulationResult:
        """Extract a single path as a standard SimulationResult."""
        echelon_results = []
        for k in range(self.n_echelons):
            echelon_results.append(
                EchelonResult(
                    name=f"E{k + 1}",
                    orders=self.orders[path_index, k],
                    inventory_levels=self.inventory[path_index, k],
                    costs=self.costs[path_index, k],
                    bullwhip_ratio=float(self.bullwhip_ratios[path_index, k]),
                    fill_rate=float(self.fill_rates[path_index, k]),
                    total_cost=float(self.total_costs[path_index, k]),
                )
            )
        return SimulationResult(
            echelon_results=echelon_results,
            cumulative_bullwhip=float(self.cumulative_bullwhip[path_index]),
            total_cost=float(self.total_costs[path_index].sum()),
        )

mean_metrics()

Average metrics across all N paths.

Source code in deepbullwhip/chain/vectorized.py 
49
50
51
52
53
54
55
56
57
58
59
def mean_metrics(self) -> dict[str, float]:
    """Average metrics across all N paths."""
    d: dict[str, float] = {}
    K = self.n_echelons
    for k in range(K):
        d[f"BW_{k + 1}"] = float(self.bullwhip_ratios[:, k].mean())
        d[f"cost_{k + 1}"] = float(self.total_costs[:, k].mean())
        d[f"fill_rate_{k + 1}"] = float(self.fill_rates[:, k].mean())
    d["BW_cumulative"] = float(self.cumulative_bullwhip.mean())
    d["total_cost"] = float(self.total_costs.sum(axis=1).mean())
    return d

to_simulation_result(path_index=0)

Extract a single path as a standard SimulationResult.

Source code in deepbullwhip/chain/vectorized.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def to_simulation_result(self, path_index: int = 0) -> SimulationResult:
    """Extract a single path as a standard SimulationResult."""
    echelon_results = []
    for k in range(self.n_echelons):
        echelon_results.append(
            EchelonResult(
                name=f"E{k + 1}",
                orders=self.orders[path_index, k],
                inventory_levels=self.inventory[path_index, k],
                costs=self.costs[path_index, k],
                bullwhip_ratio=float(self.bullwhip_ratios[path_index, k]),
                fill_rate=float(self.fill_rates[path_index, k]),
                total_cost=float(self.total_costs[path_index, k]),
            )
        )
    return SimulationResult(
        echelon_results=echelon_results,
        cumulative_bullwhip=float(self.cumulative_bullwhip[path_index]),
        total_cost=float(self.total_costs[path_index].sum()),
    )