{ "cells": [ { "cell_type": "markdown", "id": "intro", "metadata": {}, "source": [ "# Parameter constraints with `LinearConstraint`\n", "\n", "This notebook is a runnable companion to the [Constraints user guide](core/constraints.rst). It shows how to declare a `LinearConstraint`, what `verify_constraints` and `enforce_constraints` actually do, and how the residual parameterization shows up in autograd Jacobians. The final section uses the production GL06INSOUT model to confirm that the constraint system makes every derived M1 parameter invisible to autograd, and discusses how to read the numerical Jacobian comparison on the free parameters." ] }, { "cell_type": "code", "execution_count": 1, "id": "version", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:34.649075Z", "iopub.status.busy": "2026-04-10T07:49:34.648981Z", "iopub.status.idle": "2026-04-10T07:49:35.413578Z", "shell.execute_reply": "2026-04-10T07:49:35.413241Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "macrostat version: 0.5.7.post1.dev7+g35c4d7655.d20260409\n", "torch version: 2.11.0+cu130\n" ] } ], "source": [ "import macrostat\n", "import torch\n", "\n", "print('macrostat version:', macrostat.__version__)\n", "print('torch version: ', torch.__version__)" ] }, { "cell_type": "markdown", "id": "minimal-model", "metadata": {}, "source": [ "## A minimal Parameters subclass\n", "\n", "The smallest demonstration is a three-parameter group whose values must sum to one. We pick the third parameter as the residual." ] }, { "cell_type": "code", "execution_count": 2, "id": "minimal-code", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.414685Z", "iopub.status.busy": "2026-04-10T07:49:35.414534Z", "iopub.status.idle": "2026-04-10T07:49:35.523523Z", "shell.execute_reply": "2026-04-10T07:49:35.522981Z" } }, "outputs": [], "source": [ "from macrostat.core.parameters import Parameters\n", "from macrostat.core.constraints import LinearConstraint\n", "\n", "\n", "class ToyParameters(Parameters):\n", " def get_default_parameters(self):\n", " return {\n", " 'share_a': {'value': 0.4, 'lower bound': 0.0,\n", " 'upper bound': 1.0, 'notation': 'a', 'unit': '.'},\n", " 'share_b': {'value': 0.3, 'lower bound': 0.0,\n", " 'upper bound': 1.0, 'notation': 'b', 'unit': '.'},\n", " 'share_c': {'value': 0.3, 'lower bound': 0.0,\n", " 'upper bound': 1.0, 'notation': 'c', 'unit': '.'},\n", " }\n", "\n", " def get_constraints(self):\n", " return (\n", " LinearConstraint(\n", " param_names=('share_a', 'share_b', 'share_c'),\n", " target=1.0,\n", " ),\n", " )" ] }, { "cell_type": "markdown", "id": "failing-case", "metadata": {}, "source": [ "## What happens when defaults violate the constraint?\n", "\n", "Before the success case, deliberately give the toy parameters a set of defaults that do not sum to the target. The verifier emits a warning at construction time, and `enforce_constraints` then adjusts the residual parameter so the identity holds." ] }, { "cell_type": "code", "execution_count": 3, "id": "failing-code", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.524707Z", "iopub.status.busy": "2026-04-10T07:49:35.524550Z", "iopub.status.idle": "2026-04-10T07:49:35.528041Z", "shell.execute_reply": "2026-04-10T07:49:35.527776Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:macrostat.core.parameters:Constraint violation in defaults: sum(('share_a', 'share_b', 'share_c')) = 0.90000000, target = 1.00000000. Derived parameter 'share_c' will be adjusted.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "after init:\n", " share_a = 0.5\n", " share_b = 0.3\n", " share_c = 0.19999999999999996\n", "sum = 1.0\n" ] } ], "source": [ "import logging\n", "logging.basicConfig(level=logging.WARNING, force=True)\n", "\n", "class ViolatingToy(ToyParameters):\n", " def get_default_parameters(self):\n", " defaults = super().get_default_parameters()\n", " # These sum to 0.9, not 1.0\n", " defaults['share_a']['value'] = 0.5\n", " defaults['share_b']['value'] = 0.3\n", " defaults['share_c']['value'] = 0.1\n", " return defaults\n", "\n", "p_bad = ViolatingToy()\n", "print('after init:')\n", "for k in ('share_a', 'share_b', 'share_c'):\n", " print(f' {k} = {p_bad.values[k][\"value\"]}')\n", "print('sum =', sum(p_bad.values[k]['value'] for k in\n", " ('share_a', 'share_b', 'share_c')))" ] }, { "cell_type": "markdown", "id": "explain-fix", "metadata": {}, "source": [ "The warning fired during construction; the residual `share_c` was then overwritten from `0.1` to `0.2` so the group sums to `1.0`. The free parameters `share_a` and `share_b` were left untouched." ] }, { "cell_type": "markdown", "id": "happy-path", "metadata": {}, "source": [ "## The happy path: free parameters move, derived parameter follows\n", "\n", "Now use the consistent defaults from `ToyParameters`. Mutate one of the free parameters and re-run `enforce_constraints` to see the residual track." ] }, { "cell_type": "code", "execution_count": 4, "id": "happy-code", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.529023Z", "iopub.status.busy": "2026-04-10T07:49:35.528937Z", "iopub.status.idle": "2026-04-10T07:49:35.531126Z", "shell.execute_reply": "2026-04-10T07:49:35.530838Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "initial values:\n", " share_a = 0.4\n", " share_b = 0.3\n", " share_c = 0.30000000000000004\n", "\n", "after share_a := 0.55 and enforce_constraints():\n", " share_a = 0.55\n", " share_b = 0.3\n", " share_c = 0.1499999999999999\n" ] } ], "source": [ "p = ToyParameters()\n", "print('initial values:')\n", "for k in ('share_a', 'share_b', 'share_c'):\n", " print(f' {k} = {p.values[k][\"value\"]}')\n", "\n", "p.values['share_a']['value'] = 0.55\n", "p.enforce_constraints()\n", "print()\n", "print('after share_a := 0.55 and enforce_constraints():')\n", "for k in ('share_a', 'share_b', 'share_c'):\n", " print(f' {k} = {p.values[k][\"value\"]}')" ] }, { "cell_type": "code", "execution_count": 5, "id": "free-names", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.532007Z", "iopub.status.busy": "2026-04-10T07:49:35.531925Z", "iopub.status.idle": "2026-04-10T07:49:35.533641Z", "shell.execute_reply": "2026-04-10T07:49:35.533345Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "all parameters: ['share_a', 'share_b', 'share_c']\n", "free parameters: ['share_a', 'share_b']\n" ] } ], "source": [ "print('all parameters:', list(p.values.keys()))\n", "print('free parameters:', p.get_free_param_names())" ] }, { "cell_type": "markdown", "id": "autograd-section", "metadata": {}, "source": [ "## Autograd consequence: derived parameter has zero gradient\n", "\n", "Because `share_c = 1 - share_a - share_b` is encoded as a tensor operation, autograd never treats `share_c` as a free leaf. Any downstream loss has zero gradient with respect to it; the sensitivity is absorbed into `share_a` and `share_b`." ] }, { "cell_type": "code", "execution_count": null, "id": "autograd-code", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.534475Z", "iopub.status.busy": "2026-04-10T07:49:35.534370Z", "iopub.status.idle": "2026-04-10T07:49:35.538342Z", "shell.execute_reply": "2026-04-10T07:49:35.537929Z" } }, "outputs": [], "source": "p = ToyParameters()\n\n# Build leaf tensors for the FREE parameters; the residual is then\n# derived from them so it is not a leaf and has no gradient of its own.\n# Initialise the residual slot as a zero placeholder of the same dtype;\n# LinearConstraint.apply will overwrite it through the resolver.\nfree_names = ('share_a', 'share_b')\ntensors = {\n name: torch.tensor(p.values[name]['value'], requires_grad=True)\n for name in free_names\n}\ntensors['share_c'] = torch.tensor(0.0)\n\n# The step-time apply() path wants a ConstraintResolver. Build one\n# against the Parameters instance; for this scalar toy model every\n# slot resolves to its own tensor key with no index.\nresolver = p.get_constraint_resolver()\nfor c in p.get_constraints():\n c.apply(tensors, resolver) # writes share_c = 1 - share_a - share_b\n\nloss = (\n 2.0 * tensors['share_a']\n + 3.0 * tensors['share_b']\n + 4.0 * tensors['share_c']\n)\nloss.backward()\n\nfor name in ('share_a', 'share_b'):\n print(f' d(loss)/d({name}) = {tensors[name].grad.item():+.3f}')\nshare_c = tensors['share_c']\nprint(f' share_c is derived: is_leaf={share_c.is_leaf}, '\n f'requires_grad={share_c.requires_grad}')" }, { "cell_type": "markdown", "id": "explain-autograd", "metadata": {}, "source": [ "The free parameters report `2 - 4 = -2` and `3 - 4 = -1` because the residual coefficient `4` flows back through `share_c = 1 - share_a - share_b` with `d(share_c)/d(share_a) = -1`. The derived parameter `share_c` is a function of the leaves and so has no gradient of its own. This is the residual parameterization in action." ] }, { "cell_type": "markdown", "id": "real-section", "metadata": {}, "source": [ "## Real-world example: GL06INSOUT\n", "\n", "The Godley-Lavoie INSOUT model has a 4-asset Tobin portfolio with M1 cash as the buffer-stock residual. Five `LinearConstraint` instances cover the constants and the four rate-sensitivity columns.\n", "\n", "The constraint system's guarantee is narrow and mechanical: every derived M1 parameter is a tensor function of the free parameters, not a leaf, so autograd sees *exactly* zero sensitivity with respect to it. We check that first. Numerical-vs-autograd agreement on the *free* parameters is a separate question that depends on the rest of the model's numerical behaviour, which we discuss after the comparison." ] }, { "cell_type": "code", "execution_count": 7, "id": "real-load", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.539727Z", "iopub.status.busy": "2026-04-10T07:49:35.539571Z", "iopub.status.idle": "2026-04-10T07:49:35.546510Z", "shell.execute_reply": "2026-04-10T07:49:35.546142Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GL06INSOUT declares 5 constraints:\n", " 1. target=1.0\n", " free : ['WealthShareM2_Constant', 'WealthShareBills_Constant', 'WealthShareBonds_Constant']\n", " derived: WealthShareM1_Constant\n", " 2. target=0.0\n", " free : ['WealthShareM2_DepositRate', 'WealthShareBills_DepositRate', 'WealthShareBonds_DepositRate']\n", " derived: WealthShareM1_DepositRate\n", " 3. target=0.0\n", " free : ['WealthShareM2_BillRate', 'WealthShareBills_BillRate', 'WealthShareBonds_BillRate']\n", " derived: WealthShareM1_BillRate\n", " 4. target=0.0\n", " free : ['WealthShareM2_BondYield', 'WealthShareBills_BondYield', 'WealthShareBonds_BondYield']\n", " derived: WealthShareM1_BondYield\n", " 5. target=0.0\n", " free : ['WealthShareM2_Income', 'WealthShareBills_Income', 'WealthShareBonds_Income']\n", " derived: WealthShareM1_Income\n" ] } ], "source": [ "from macrostat.models import get_model\n", "\n", "model = get_model('GL06INSOUT')()\n", "constraints = model.parameters.get_constraints()\n", "print(f'GL06INSOUT declares {len(constraints)} constraints:')\n", "for i, c in enumerate(constraints, 1):\n", " print(f' {i}. target={c.target}')\n", " print(f' free : {list(c.free_params)}')\n", " print(f' derived: {c.derived_param}')" ] }, { "cell_type": "code", "execution_count": 8, "id": "real-jacobians", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:35.547536Z", "iopub.status.busy": "2026-04-10T07:49:35.547443Z", "iopub.status.idle": "2026-04-10T07:49:42.098196Z", "shell.execute_reply": "2026-04-10T07:49:42.097842Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jacobian comparison: autograd vs numerical (log, eps=1e-3)\n", "Parameters compared: 5\n", "Overall max abs diff: 1.899e+00 max rel diff: 2.077e-01\n", "\n", "Parameter max_abs mean_abs max_rel mean_rel close nan/inf\n", "---------------------------------------------------------------------------------------------------------\n", "WealthShareBills_Constant 1.899e+00 1.899e+00 2.077e-01 2.077e-01 0/1 no\n", "WealthShareM2_Constant 5.841e-02 5.841e-02 5.841e-02 5.841e-02 0/1 no\n", "WealthShareBonds_Constant 1.534e-01 1.534e-01 7.729e-03 7.729e-03 0/1 no\n", "InventorySalesRatioBaseline 9.693e-02 9.693e-02 3.941e-03 3.941e-03 0/1 no\n", "TaxRate 3.356e-01 3.356e-01 1.868e-04 1.868e-04 1/1 no\n" ] } ], "source": [ "import sys, os\n", "# Make the doc helper module importable from the notebook\n", "sys.path.insert(0, os.path.abspath('.'))\n", "from _doc_helpers import mean_nominal_output\n", "\n", "from macrostat.diff import (\n", " JacobianAutograd,\n", " JacobianNumerical,\n", " compare_jacobian_dicts,\n", ")\n", "\n", "# Hand-picked subset: a few free parameters from the Tobin matrix\n", "# (where the constraint system matters) plus two unrelated parameters\n", "# as a sanity check. A full sweep over all 42 free parameters takes\n", "# minutes; for a worked example, five is enough to demonstrate the\n", "# pattern.\n", "sample_params = [\n", " 'WealthShareM2_Constant',\n", " 'WealthShareBills_Constant',\n", " 'WealthShareBonds_Constant',\n", " 'TaxRate',\n", " 'InventorySalesRatioBaseline',\n", "]\n", "\n", "jac_auto = JacobianAutograd(model).compute(\n", " loss_fn=mean_nominal_output, mode='rev',\n", ")\n", "jac_num = JacobianNumerical(model).compute(\n", " loss_fn=mean_nominal_output, mode='central', num_workers=1,\n", " param_names=sample_params,\n", ")\n", "\n", "jac_auto_subset = {k: jac_auto[k] for k in sample_params}\n", "\n", "report = compare_jacobian_dicts(\n", " jac_auto_subset, jac_num,\n", " method_a='autograd', method_b='numerical (log, eps=1e-3)',\n", " rtol=1e-3, atol=1e-6,\n", ")\n", "print(report.summary())" ] }, { "cell_type": "markdown", "id": "real-jacobians-interp", "metadata": {}, "source": "`TaxRate` agrees to four significant digits. The two unrelated parameters and `WealthShareBonds_Constant` agree to roughly 1%; the two other Tobin constants show disagreements of 6% and 21%. This residual disagreement is *not* a constraint-system bug. GL06INSOUT's `portfolio_switches` method uses a step-function `torch.where` to select between the M1- and M2-dominated regimes, and any parameter that moves the economy near that boundary (which the Tobin constants do by shifting the baseline portfolio composition) straddles the discontinuity under finite differences. Autograd reports zero gradient at the switch; numerical differences see a finite jump. See the \"Expected disagreements\" section of [Differentiability and Jacobian Tools](diff.rst) and the worked example in [Jacobian Diagnostics](jacobian_diagnostics.ipynb).\n\nThe next cell runs the check that *is* the constraint-system guarantee: every derived M1 parameter should have exactly zero autograd gradient." }, { "cell_type": "code", "execution_count": 9, "id": "derived-zero", "metadata": { "execution": { "iopub.execute_input": "2026-04-10T07:49:42.099320Z", "iopub.status.busy": "2026-04-10T07:49:42.099173Z", "iopub.status.idle": "2026-04-10T07:49:42.101310Z", "shell.execute_reply": "2026-04-10T07:49:42.101030Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " WealthShareM1_Constant max|grad| = 0.000e+00\n", " WealthShareM1_DepositRate max|grad| = 0.000e+00\n", " WealthShareM1_BillRate max|grad| = 0.000e+00\n", " WealthShareM1_BondYield max|grad| = 0.000e+00\n", " WealthShareM1_Income max|grad| = 0.000e+00\n" ] } ], "source": [ "# Confirm that every derived M1 parameter has exactly zero autograd gradient\n", "derived = [c.derived_param for c in constraints]\n", "for name in derived:\n", " g = jac_auto[name]\n", " print(f' {name:<35s} max|grad| = {g.abs().max().item():.3e}')" ] }, { "cell_type": "markdown", "id": "closing", "metadata": {}, "source": "## Closing notes\n\n* The constraint system is opt-in: existing models that do not override `get_constraints` are unaffected.\n* `LinearConstraint.apply` runs inside `Behavior.apply_parameter_shocks` at every step, fetched fresh against the current `Parameters` instance, so scenario shocks to free parameters automatically propagate to the residual. The init-time adjustment happens via `Parameters.enforce_constraints`.\n* For the API and the verification rules, see [Constraints](core/constraints.rst).\n* For diagnosing autograd vs numerical disagreements that are *not* explained by constraints, see [Jacobian Diagnostics](jacobian_diagnostics.ipynb)." } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.5" } }, "nbformat": 4, "nbformat_minor": 5 }