From 0517f9f3959bcdcb90101d24863364e9f0ebcff5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Moran?= Date: Fri, 6 Mar 2026 19:35:34 +0100 Subject: [PATCH 1/3] Fix BundledLeontief input usage to distribute within bundles by stock BundledLeontief.compute_intermediate_inputs_used() and compute_capital_inputs_used() ignored the substitution bundle matrix, using fixed Leontief coefficients for each good independently. This caused stock to deplete for expensive goods and accumulate for cheap goods, since buying was adjusted by substitution weights but usage was not. Now distributes total bundle usage proportionally to available stock within each multi-member bundle, so firms use more of whichever substitute they have on hand. Co-Authored-By: Claude Opus 4.6 --- macromodel/agents/firms/func/production.py | 56 ++++++++-- .../test_firms/func/test_production.py | 101 ++++++++++++++++++ 2 files changed, 147 insertions(+), 10 deletions(-) diff --git a/macromodel/agents/firms/func/production.py b/macromodel/agents/firms/func/production.py index 10b7390a..ae23aea5 100644 --- a/macromodel/agents/firms/func/production.py +++ b/macromodel/agents/firms/func/production.py @@ -945,8 +945,10 @@ def compute_intermediate_inputs_used( ) -> np.ndarray: """Calculate intermediate inputs used under bundled Leontief technology. - Input usage is proportional to production with fixed input-output coefficients, - adjusted by the substitution bundle matrix. + Input usage is proportional to production with fixed input-output coefficients. + Within substitution bundles, the total bundle usage is distributed + proportionally to available stock, reflecting that firms use more of + whichever substitute they have on hand. Args: realised_production (np.ndarray): Actual production achieved @@ -960,7 +962,7 @@ def compute_intermediate_inputs_used( Returns: np.ndarray: Intermediate inputs used in production """ - # Calculate base input usage + # Calculate base input usage per good (Leontief coefficients) used_inputs = np.divide( realised_production[:, None], intermediate_inputs_productivity_matrix, @@ -968,9 +970,27 @@ def compute_intermediate_inputs_used( where=intermediate_inputs_productivity_matrix != 0.0, ) - # Apply substitution bundles to determine actual input usage - # This is a simplified approach - in a full implementation, - # you might want to optimize input usage across bundles + # For multi-member bundles, distribute total bundle usage + # proportionally to stock within the bundle + bundle_members_per_col = (substitution_bundle_matrix > 0).sum(axis=0) + multi_member_bundles = np.where(bundle_members_per_col > 1)[0] + + for b in multi_member_bundles: + member_indices = np.where(substitution_bundle_matrix[:, b] > 0)[0] + # Total usage for this bundle = sum of Leontief usage across member goods + total_bundle_usage = used_inputs[:, member_indices].sum(axis=1, keepdims=True) + # Stock of member goods + member_stock = intermediate_inputs_stock[:, member_indices] + member_stock_sum = member_stock.sum(axis=1, keepdims=True) + # Distribute proportionally to stock; if no stock, distribute equally + has_stock = member_stock_sum > 0 + stock_share = np.where( + has_stock, + member_stock / np.maximum(member_stock_sum, 1e-30), + 1.0 / len(member_indices), + ) + used_inputs[:, member_indices] = total_bundle_usage * stock_share + return used_inputs def compute_capital_inputs_used( @@ -984,7 +1004,8 @@ def compute_capital_inputs_used( """Calculate capital depreciation under bundled Leontief technology. Capital depreciation is proportional to production with fixed - depreciation rates, adjusted by the substitution bundle matrix. + depreciation rates. Within substitution bundles, the total bundle + depreciation is distributed proportionally to available stock. Args: realised_production (np.ndarray): Actual production achieved @@ -1003,7 +1024,22 @@ def compute_capital_inputs_used( used_capital_inputs[used_capital_inputs == np.inf] = 0.0 used_capital_inputs[used_capital_inputs == -np.inf] = 0.0 - # Apply substitution bundles to determine actual capital usage - # This is a simplified approach - in a full implementation, - # you might want to optimize capital usage across bundles + # For multi-member bundles, distribute total bundle depreciation + # proportionally to stock within the bundle + bundle_members_per_col = (substitution_bundle_matrix > 0).sum(axis=0) + multi_member_bundles = np.where(bundle_members_per_col > 1)[0] + + for b in multi_member_bundles: + member_indices = np.where(substitution_bundle_matrix[:, b] > 0)[0] + total_bundle_usage = used_capital_inputs[:, member_indices].sum(axis=1, keepdims=True) + member_stock = capital_inputs_stock[:, member_indices] + member_stock_sum = member_stock.sum(axis=1, keepdims=True) + has_stock = member_stock_sum > 0 + stock_share = np.where( + has_stock, + member_stock / np.maximum(member_stock_sum, 1e-30), + 1.0 / len(member_indices), + ) + used_capital_inputs[:, member_indices] = total_bundle_usage * stock_share + return used_capital_inputs diff --git a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py index 6eb8be9c..93240373 100644 --- a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py +++ b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py @@ -131,3 +131,104 @@ def test_target_intermediate_inputs_bundle_empty(self): expected = np.full((n_industries, n_industries), 1) assert np.array_equal(result, expected) + + +class TestBundledLeontiefUsage: + """Test that BundledLeontief distributes input usage within bundles + proportionally to stock, not at fixed Leontief rates.""" + + def test_usage_distributed_by_stock_within_bundle(self): + """If two goods are in a bundle and one has 3x the stock, + it should bear 3x the usage. Before the fix, both goods + were used at the same fixed Leontief rate regardless of stock.""" + n_firms = 2 + n_goods = 3 # goods 0 and 1 in a bundle, good 2 is singleton + + # Leontief coefficients: each firm needs 1 unit of each good per unit of output + productivity = np.full((n_firms, n_goods), 1.0) + + # Firm 0 has substituted toward good 0 (3x stock vs good 1) + # Firm 1 has equal stock + stock = np.array([ + [3.0, 1.0, 2.0], # firm 0: 3x more of good 0 than good 1 + [2.0, 2.0, 2.0], # firm 1: equal stock + ]) + + production = np.array([2.0, 2.0]) + + # Bundle matrix: goods 0,1 in bundle 0; good 2 in bundle 1 + # Shape (n_goods, n_bundles) = (3, 2) + bundle_matrix = np.array([ + [0.5, 0.0], # good 0 in bundle 0 + [0.5, 0.0], # good 1 in bundle 0 + [0.0, 1.0], # good 2 in bundle 1 (singleton) + ]) + + criticality = np.ones((n_firms, n_goods)) + + bundled = BundledLeontief() + used = bundled.compute_intermediate_inputs_used( + realised_production=production, + intermediate_inputs_productivity_matrix=productivity, + intermediate_inputs_stock=stock, + goods_criticality_matrix=criticality, + substitution_bundle_matrix=bundle_matrix, + ) + + # Total bundle usage for bundle 0: production/coeff[0] + production/coeff[1] = 2+2 = 4 + # Firm 0: stock shares are 3/4 and 1/4 → used = [3.0, 1.0] + # Firm 1: stock shares are 2/4 and 2/4 → used = [2.0, 2.0] + # Good 2 (singleton): unchanged at production/coeff = 2.0 + expected = np.array([ + [3.0, 1.0, 2.0], + [2.0, 2.0, 2.0], + ]) + assert np.allclose(used, expected) + + def test_usage_total_preserved_within_bundle(self): + """Total usage across bundle members should equal the sum of + individual Leontief requirements, just redistributed.""" + n_firms = 3 + n_goods = 4 # goods 0,1,2 in a bundle, good 3 singleton + + productivity = np.array([ + [2.0, 4.0, 5.0, 3.0], + [2.0, 4.0, 5.0, 3.0], + [2.0, 4.0, 5.0, 3.0], + ]) + + stock = np.array([ + [10.0, 5.0, 5.0, 8.0], + [1.0, 1.0, 18.0, 8.0], + [0.0, 0.0, 0.0, 8.0], # no stock at all + ]) + + production = np.array([10.0, 10.0, 10.0]) + + bundle_matrix = np.array([ + [1 / 3, 0.0], + [1 / 3, 0.0], + [1 / 3, 0.0], + [0.0, 1.0], + ]) + + criticality = np.ones((n_firms, n_goods)) + + bundled = BundledLeontief() + used = bundled.compute_intermediate_inputs_used( + realised_production=production, + intermediate_inputs_productivity_matrix=productivity, + intermediate_inputs_stock=stock, + goods_criticality_matrix=criticality, + substitution_bundle_matrix=bundle_matrix, + ) + + # Total Leontief usage per firm for bundle goods: 10/2 + 10/4 + 10/5 = 5+2.5+2 = 9.5 + bundle_total = used[:, :3].sum(axis=1) + assert np.allclose(bundle_total, 9.5) + + # Singleton good 3 unchanged: 10/3 + assert np.allclose(used[:, 3], 10.0 / 3.0) + + # Firm 2 has zero stock → equal distribution: 9.5/3 each + assert np.allclose(used[2, :3], 9.5 / 3.0) From b31a1fac8b321fbca2019f0ea67c8eb362cda4b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Moran?= Date: Fri, 6 Mar 2026 19:44:47 +0100 Subject: [PATCH 2/3] Style --- .../test_firms/func/test_production.py | 70 +++++++++++-------- 1 file changed, 41 insertions(+), 29 deletions(-) diff --git a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py index 2ae6ee54..42a8e468 100644 --- a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py +++ b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py @@ -162,20 +162,24 @@ def test_usage_distributed_by_stock_within_bundle(self): # Firm 0 has substituted toward good 0 (3x stock vs good 1) # Firm 1 has equal stock - stock = np.array([ - [3.0, 1.0, 2.0], # firm 0: 3x more of good 0 than good 1 - [2.0, 2.0, 2.0], # firm 1: equal stock - ]) + stock = np.array( + [ + [3.0, 1.0, 2.0], # firm 0: 3x more of good 0 than good 1 + [2.0, 2.0, 2.0], # firm 1: equal stock + ] + ) production = np.array([2.0, 2.0]) # Bundle matrix: goods 0,1 in bundle 0; good 2 in bundle 1 # Shape (n_goods, n_bundles) = (3, 2) - bundle_matrix = np.array([ - [0.5, 0.0], # good 0 in bundle 0 - [0.5, 0.0], # good 1 in bundle 0 - [0.0, 1.0], # good 2 in bundle 1 (singleton) - ]) + bundle_matrix = np.array( + [ + [0.5, 0.0], # good 0 in bundle 0 + [0.5, 0.0], # good 1 in bundle 0 + [0.0, 1.0], # good 2 in bundle 1 (singleton) + ] + ) criticality = np.ones((n_firms, n_goods)) @@ -192,10 +196,12 @@ def test_usage_distributed_by_stock_within_bundle(self): # Firm 0: stock shares are 3/4 and 1/4 → used = [3.0, 1.0] # Firm 1: stock shares are 2/4 and 2/4 → used = [2.0, 2.0] # Good 2 (singleton): unchanged at production/coeff = 2.0 - expected = np.array([ - [3.0, 1.0, 2.0], - [2.0, 2.0, 2.0], - ]) + expected = np.array( + [ + [3.0, 1.0, 2.0], + [2.0, 2.0, 2.0], + ] + ) assert np.allclose(used, expected) def test_usage_total_preserved_within_bundle(self): @@ -204,26 +210,32 @@ def test_usage_total_preserved_within_bundle(self): n_firms = 3 n_goods = 4 # goods 0,1,2 in a bundle, good 3 singleton - productivity = np.array([ - [2.0, 4.0, 5.0, 3.0], - [2.0, 4.0, 5.0, 3.0], - [2.0, 4.0, 5.0, 3.0], - ]) + productivity = np.array( + [ + [2.0, 4.0, 5.0, 3.0], + [2.0, 4.0, 5.0, 3.0], + [2.0, 4.0, 5.0, 3.0], + ] + ) - stock = np.array([ - [10.0, 5.0, 5.0, 8.0], - [1.0, 1.0, 18.0, 8.0], - [0.0, 0.0, 0.0, 8.0], # no stock at all - ]) + stock = np.array( + [ + [10.0, 5.0, 5.0, 8.0], + [1.0, 1.0, 18.0, 8.0], + [0.0, 0.0, 0.0, 8.0], # no stock at all + ] + ) production = np.array([10.0, 10.0, 10.0]) - bundle_matrix = np.array([ - [1 / 3, 0.0], - [1 / 3, 0.0], - [1 / 3, 0.0], - [0.0, 1.0], - ]) + bundle_matrix = np.array( + [ + [1 / 3, 0.0], + [1 / 3, 0.0], + [1 / 3, 0.0], + [0.0, 1.0], + ] + ) criticality = np.ones((n_firms, n_goods)) From 3b980d86d7734fdc0b32a03791f015a4fe6403fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Moran?= Date: Mon, 9 Mar 2026 22:48:45 +0100 Subject: [PATCH 3/3] Fix overconsumption in bundle --- macromodel/agents/firms/func/production.py | 70 +++++++--------- .../test_firms/func/test_production.py | 82 ++++++++++++++----- 2 files changed, 89 insertions(+), 63 deletions(-) diff --git a/macromodel/agents/firms/func/production.py b/macromodel/agents/firms/func/production.py index ae23aea5..7d5a1f26 100644 --- a/macromodel/agents/firms/func/production.py +++ b/macromodel/agents/firms/func/production.py @@ -945,15 +945,18 @@ def compute_intermediate_inputs_used( ) -> np.ndarray: """Calculate intermediate inputs used under bundled Leontief technology. - Input usage is proportional to production with fixed input-output coefficients. - Within substitution bundles, the total bundle usage is distributed - proportionally to available stock, reflecting that firms use more of - whichever substitute they have on hand. + For singleton goods, usage follows the standard Leontief formula: + used = production / coefficient. + + Within multi-member substitution bundles, total usage is set so that + the effective input consumed (sum of weight * coefficient * used) + equals production, and is distributed proportionally to stock: + used[j] = production * stock[j] / bundle_capacity. Args: realised_production (np.ndarray): Actual production achieved intermediate_inputs_productivity_matrix (np.ndarray): Input-output - coefficients + coefficients (productivity: output per unit of input) intermediate_inputs_stock (np.ndarray): Available input stocks goods_criticality_matrix (np.ndarray): Input criticality levels substitution_bundle_matrix (np.ndarray): Matrix defining substitution @@ -962,7 +965,7 @@ def compute_intermediate_inputs_used( Returns: np.ndarray: Intermediate inputs used in production """ - # Calculate base input usage per good (Leontief coefficients) + # Base Leontief usage for singleton goods: used = production / coefficient used_inputs = np.divide( realised_production[:, None], intermediate_inputs_productivity_matrix, @@ -970,26 +973,30 @@ def compute_intermediate_inputs_used( where=intermediate_inputs_productivity_matrix != 0.0, ) - # For multi-member bundles, distribute total bundle usage - # proportionally to stock within the bundle + # For multi-member bundles, distribute usage proportionally to stock + # such that the effective bundle contribution equals production bundle_members_per_col = (substitution_bundle_matrix > 0).sum(axis=0) multi_member_bundles = np.where(bundle_members_per_col > 1)[0] for b in multi_member_bundles: member_indices = np.where(substitution_bundle_matrix[:, b] > 0)[0] - # Total usage for this bundle = sum of Leontief usage across member goods - total_bundle_usage = used_inputs[:, member_indices].sum(axis=1, keepdims=True) - # Stock of member goods member_stock = intermediate_inputs_stock[:, member_indices] - member_stock_sum = member_stock.sum(axis=1, keepdims=True) - # Distribute proportionally to stock; if no stock, distribute equally - has_stock = member_stock_sum > 0 - stock_share = np.where( - has_stock, - member_stock / np.maximum(member_stock_sum, 1e-30), - 1.0 / len(member_indices), + member_coeff = intermediate_inputs_productivity_matrix[:, member_indices] + member_weights = substitution_bundle_matrix[member_indices, b] + + # Only include goods with finite positive coefficients in bundle capacity. + # inf coefficient means the firm doesn't use that input at all. + finite_mask = np.isfinite(member_coeff) & (member_coeff > 0) + effective = np.where(finite_mask, member_weights * member_coeff * member_stock, 0.0) + bundle_cap = effective.sum(axis=1, keepdims=True) + + # used[j] = production * stock[j] / bundle_capacity + # Only for goods the firm actually uses (finite coefficient) + used_inputs[:, member_indices] = np.where( + finite_mask & (bundle_cap > 0), + realised_production[:, None] * member_stock / np.maximum(bundle_cap, 1e-30), + 0.0, ) - used_inputs[:, member_indices] = total_bundle_usage * stock_share return used_inputs @@ -1004,8 +1011,9 @@ def compute_capital_inputs_used( """Calculate capital depreciation under bundled Leontief technology. Capital depreciation is proportional to production with fixed - depreciation rates. Within substitution bundles, the total bundle - depreciation is distributed proportionally to available stock. + depreciation rates. Bundle adjustment is not applied here because + depreciation (production * rate) does not have the same + overconsumption issue as intermediate inputs (production / coefficient). Args: realised_production (np.ndarray): Actual production achieved @@ -1019,27 +1027,7 @@ def compute_capital_inputs_used( Returns: np.ndarray: Capital inputs depreciated in production """ - # Calculate base capital depreciation used_capital_inputs = realised_production[:, None] * capital_inputs_depreciation_matrix used_capital_inputs[used_capital_inputs == np.inf] = 0.0 used_capital_inputs[used_capital_inputs == -np.inf] = 0.0 - - # For multi-member bundles, distribute total bundle depreciation - # proportionally to stock within the bundle - bundle_members_per_col = (substitution_bundle_matrix > 0).sum(axis=0) - multi_member_bundles = np.where(bundle_members_per_col > 1)[0] - - for b in multi_member_bundles: - member_indices = np.where(substitution_bundle_matrix[:, b] > 0)[0] - total_bundle_usage = used_capital_inputs[:, member_indices].sum(axis=1, keepdims=True) - member_stock = capital_inputs_stock[:, member_indices] - member_stock_sum = member_stock.sum(axis=1, keepdims=True) - has_stock = member_stock_sum > 0 - stock_share = np.where( - has_stock, - member_stock / np.maximum(member_stock_sum, 1e-30), - 1.0 / len(member_indices), - ) - used_capital_inputs[:, member_indices] = total_bundle_usage * stock_share - return used_capital_inputs diff --git a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py index 42a8e468..fe2eaeac 100644 --- a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py +++ b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py @@ -151,17 +151,14 @@ class TestBundledLeontiefUsage: proportionally to stock, not at fixed Leontief rates.""" def test_usage_distributed_by_stock_within_bundle(self): - """If two goods are in a bundle and one has 3x the stock, - it should bear 3x the usage. Before the fix, both goods - were used at the same fixed Leontief rate regardless of stock.""" + """Within a bundle, usage is proportional to stock. The effective + input consumed (sum of weight * coeff * used) must equal production.""" n_firms = 2 n_goods = 3 # goods 0 and 1 in a bundle, good 2 is singleton - # Leontief coefficients: each firm needs 1 unit of each good per unit of output + # Equal coefficients: 1 unit of output per unit of input productivity = np.full((n_firms, n_goods), 1.0) - # Firm 0 has substituted toward good 0 (3x stock vs good 1) - # Firm 1 has equal stock stock = np.array( [ [3.0, 1.0, 2.0], # firm 0: 3x more of good 0 than good 1 @@ -172,12 +169,11 @@ def test_usage_distributed_by_stock_within_bundle(self): production = np.array([2.0, 2.0]) # Bundle matrix: goods 0,1 in bundle 0; good 2 in bundle 1 - # Shape (n_goods, n_bundles) = (3, 2) bundle_matrix = np.array( [ - [0.5, 0.0], # good 0 in bundle 0 - [0.5, 0.0], # good 1 in bundle 0 - [0.0, 1.0], # good 2 in bundle 1 (singleton) + [0.5, 0.0], + [0.5, 0.0], + [0.0, 1.0], ] ) @@ -192,10 +188,10 @@ def test_usage_distributed_by_stock_within_bundle(self): substitution_bundle_matrix=bundle_matrix, ) - # Total bundle usage for bundle 0: production/coeff[0] + production/coeff[1] = 2+2 = 4 - # Firm 0: stock shares are 3/4 and 1/4 → used = [3.0, 1.0] - # Firm 1: stock shares are 2/4 and 2/4 → used = [2.0, 2.0] - # Good 2 (singleton): unchanged at production/coeff = 2.0 + # Bundle cap firm 0: 0.5*(1*3 + 1*1) = 2.0 + # used[j] = production * stock[j] / cap = 2*stock[j]/2 = stock[j] + # Bundle cap firm 1: 0.5*(1*2 + 1*2) = 2.0, same result + # Good 2 (singleton): production / coeff = 2.0 expected = np.array( [ [3.0, 1.0, 2.0], @@ -204,9 +200,46 @@ def test_usage_distributed_by_stock_within_bundle(self): ) assert np.allclose(used, expected) - def test_usage_total_preserved_within_bundle(self): - """Total usage across bundle members should equal the sum of - individual Leontief requirements, just redistributed.""" + # Verify effective contribution = production for the bundle + bundle_contrib = (used[:, :2] * productivity[:, :2] * 0.5).sum(axis=1) + assert np.allclose(bundle_contrib, production) + + def test_usage_never_exceeds_stock(self): + """With varying coefficients, usage must never exceed available stock. + This was the overconsumption bug: summing production/coeff across all + bundle members produced usage exceeding total stock.""" + n_firms = 1 + n_goods = 2 # both in one bundle + + # Good 0 very productive, good 1 not + productivity = np.array([[2.0, 0.5]]) + stock = np.array([[3.0, 1.0]]) + + # Bundle capacity: 0.5*(2*3 + 0.5*1) = 3.25 + production = np.array([3.25]) # bundle is binding + + bundle_matrix = np.array([[0.5], [0.5]]) + criticality = np.ones((n_firms, n_goods)) + + bundled = BundledLeontief() + used = bundled.compute_intermediate_inputs_used( + realised_production=production, + intermediate_inputs_productivity_matrix=productivity, + intermediate_inputs_stock=stock, + goods_criticality_matrix=criticality, + substitution_bundle_matrix=bundle_matrix, + ) + + # Must not exceed stock + assert np.all(used <= stock + 1e-10) + + # Effective contribution must equal production + bundle_contrib = (used * productivity * 0.5).sum(axis=1) + assert np.allclose(bundle_contrib, production) + + def test_usage_effective_contribution_equals_production(self): + """The effective input contribution (sum of weight*coeff*used) + must equal production for each bundle, with varying coefficients.""" n_firms = 3 n_goods = 4 # goods 0,1,2 in a bundle, good 3 singleton @@ -248,12 +281,17 @@ def test_usage_total_preserved_within_bundle(self): substitution_bundle_matrix=bundle_matrix, ) - # Total Leontief usage per firm for bundle goods: 10/2 + 10/4 + 10/5 = 5+2.5+2 = 9.5 - bundle_total = used[:, :3].sum(axis=1) - assert np.allclose(bundle_total, 9.5) + # Effective bundle contribution must equal production for firms with stock + weights = np.array([1 / 3, 1 / 3, 1 / 3]) + for firm in [0, 1]: + bundle_contrib = (used[firm, :3] * productivity[firm, :3] * weights).sum() + assert np.allclose(bundle_contrib, production[firm]) # Singleton good 3 unchanged: 10/3 assert np.allclose(used[:, 3], 10.0 / 3.0) - # Firm 2 has zero stock → equal distribution: 9.5/3 each - assert np.allclose(used[2, :3], 9.5 / 3.0) + # Firm 2 has zero stock -> zero usage + assert np.allclose(used[2, :3], 0.0) + + # Usage must never exceed stock + assert np.all(used <= stock + 1e-10)