diff --git a/demos/true_measure_domain_inclusion.ipynb b/demos/true_measure_domain_inclusion.ipynb new file mode 100644 index 000000000..a2f18c171 --- /dev/null +++ b/demos/true_measure_domain_inclusion.ipynb @@ -0,0 +1,1607 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a4c44dd5", + "metadata": {}, + "source": [ + "# Chaining `TrueMeasure` Transformations When Domains Are Compatible\n", + "\n", + "A `TrueMeasure` transformation takes values from one set, transforms them,\n", + "and passes the result to the next transformation in a chain.\n", + "\n", + "The important question is:\n", + "\n", + "> **Does the next transformation need the previous range to be exactly equal\n", + "> to its domain, or does it only need every incoming value to be valid?**\n", + "\n", + "This notebook shows why **containment** is the correct condition.\n", + "\n", + "We first look at a simple real-world example, reproduce what the previous\n", + "QMCPy compatibility rule would do, and then show the same transformation\n", + "working with the updated rule." + ] + }, + { + "cell_type": "markdown", + "id": "d94f50a6", + "metadata": {}, + "source": [ + "## 1. A real-world picture\n", + "\n", + "Suppose we model the **operating load of a battery system** as a normalized\n", + "fraction between `0` and `1`.\n", + "\n", + "A downstream response model is designed to accept any normalized load\n", + "\n", + "$$\n", + "x \\in [0,1].\n", + "$$\n", + "\n", + "However, under normal operating conditions, the battery is deliberately kept\n", + "between 25% and 75% load:\n", + "\n", + "$$\n", + "x \\in [0.25,0.75].\n", + "$$\n", + "\n", + "There is no incompatibility here.\n", + "\n", + "The downstream model accepts every value from `0` to `1`, and the upstream\n", + "stage produces only values between `0.25` and `0.75`.\n", + "\n", + "Mathematically,\n", + "\n", + "$$\n", + "[0.25,0.75] \\subseteq [0,1].\n", + "$$\n", + "\n", + "The fact that\n", + "\n", + "$$\n", + "[0.25,0.75] \\neq [0,1]\n", + "$$\n", + "\n", + "does **not** make the two stages incompatible.\n", + "\n", + "This is exactly the situation addressed by this change in QMCPy." + ] + }, + { + "cell_type": "markdown", + "id": "8a14030c", + "metadata": {}, + "source": [ + "## 2. The mathematical compatibility condition\n", + "\n", + "Consider two consecutive transformations\n", + "\n", + "$$\n", + "T_{j-1}:D_{j-1}\\rightarrow R_{j-1},\n", + "$$\n", + "\n", + "and\n", + "\n", + "$$\n", + "T_j:D_j\\rightarrow R_j.\n", + "$$\n", + "\n", + "The output of the first transformation becomes the input of the second:\n", + "\n", + "$$\n", + "x \\in R_{j-1}\n", + "\\quad\\longrightarrow\\quad\n", + "T_j(x).\n", + "$$\n", + "\n", + "Therefore, the second transformation is valid whenever **every possible\n", + "output of the first transformation lies inside the domain of the second**:\n", + "\n", + "$$\n", + "\\boxed{R_{j-1}\\subseteq D_j}.\n", + "$$\n", + "\n", + "The previous QMCPy check effectively required\n", + "\n", + "$$\n", + "R_{j-1}=D_j.\n", + "$$\n", + "\n", + "Equality is sufficient, but it is stronger than necessary.\n", + "\n", + "For one-dimensional intervals,\n", + "\n", + "$$\n", + "R_{j-1}=[r_L,r_U],\n", + "\\qquad\n", + "D_j=[d_L,d_U],\n", + "$$\n", + "\n", + "containment means\n", + "\n", + "$$\n", + "\\boxed{\n", + "d_L\\le r_L\n", + "\\quad\\text{and}\\quad\n", + "r_U\\le d_U\n", + "}.\n", + "$$\n", + "\n", + "For our battery example,\n", + "\n", + "$$\n", + "0\\le0.25\n", + "\\qquad\\text{and}\\qquad\n", + "0.75\\le1,\n", + "$$\n", + "\n", + "so\n", + "\n", + "$$\n", + "[0.25,0.75]\\subseteq[0,1].\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "6c4c9d29", + "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/QMCSoftware/QMCSoftware/blob/develop/demos/true_measure_domain_inclusion.ipynb)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "15ab6588", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-05T15:28:35.491526Z", + "iopub.status.busy": "2026-09-05T15:28:35.491162Z", + "iopub.status.idle": "2026-09-05T15:28:35.500428Z", + "shell.execute_reply": "2026-09-05T15:28:35.499277Z" + } + }, + "outputs": [], + "source": [ + "# @title Execute this cell to install dependencies\n", + "try:\n", + " import google.colab\n", + " IN_COLAB = True\n", + "except ImportError:\n", + " IN_COLAB = False\n", + "if IN_COLAB:\n", + " !pip install -q qmcpy\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8d512296", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-05T15:28:35.503295Z", + "iopub.status.busy": "2026-09-05T15:28:35.503074Z", + "iopub.status.idle": "2026-09-05T15:28:36.936871Z", + "shell.execute_reply": "2026-09-05T15:28:36.936105Z" + } + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "from qmcpy import DigitalNetB2, Kumaraswamy, Uniform\n", + "from qmcpy.util import ParameterError" + ] + }, + { + "cell_type": "markdown", + "id": "85f8195a", + "metadata": {}, + "source": [ + "## 3. What happened before this change?\n", + "\n", + "The previous compatibility check compared the two sets using equality.\n", + "\n", + "Conceptually, it asked\n", + "\n", + "$$\n", + "R_{j-1}=D_j\\;?\n", + "$$\n", + "\n", + "For the battery example, that becomes\n", + "\n", + "$$\n", + "[0.25,0.75]=[0,1]\\;?\n", + "$$\n", + "\n", + "which is false.\n", + "\n", + "So the previous rule treated the chain as incompatible even though every\n", + "value produced by the first stage was a perfectly valid input for the second.\n", + "\n", + "The code below reproduces that previous equality rule so that we can compare\n", + "it directly with the mathematically correct containment rule." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ed148c64", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-05T15:28:36.938578Z", + "iopub.status.busy": "2026-09-05T15:28:36.938314Z", + "iopub.status.idle": "2026-09-05T15:28:36.951949Z", + "shell.execute_reply": "2026-09-05T15:28:36.951483Z" + } + }, + "outputs": [], + "source": [ + "def previous_equality_rule(transform_range, next_domain):\n", + " \"\"\"Reproduce the previous equality-based compatibility decision.\"\"\"\n", + " transform_range = np.asarray(transform_range, dtype=float)\n", + " next_domain = np.asarray(next_domain, dtype=float)\n", + "\n", + " try:\n", + " transform_range, next_domain = np.broadcast_arrays(\n", + " transform_range,\n", + " next_domain,\n", + " )\n", + " except ValueError:\n", + " return False\n", + "\n", + " return bool(np.all(transform_range == next_domain))\n", + "\n", + "\n", + "def containment_rule(transform_range, next_domain):\n", + " \"\"\"Check whether every transform output lies inside the next domain.\"\"\"\n", + " transform_range = np.asarray(transform_range, dtype=float)\n", + " next_domain = np.asarray(next_domain, dtype=float)\n", + "\n", + " try:\n", + " transform_range, next_domain = np.broadcast_arrays(\n", + " transform_range,\n", + " next_domain,\n", + " )\n", + " except ValueError:\n", + " return False\n", + "\n", + " lower_bounds_valid = np.all(\n", + " next_domain[:, 0] <= transform_range[:, 0]\n", + " )\n", + " upper_bounds_valid = np.all(\n", + " transform_range[:, 1] <= next_domain[:, 1]\n", + " )\n", + "\n", + " return bool(lower_bounds_valid and upper_bounds_valid)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e20bf2af", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Upstream operating range: [[0.25 0.75]]\n", + "Downstream accepted domain: [[0. 1.]]\n", + "\n", + "Previous equality rule accepts the chain: False\n", + "Containment rule accepts the chain: True\n" + ] + } + ], + "source": [ + "battery_operating_range = np.array([[0.25, 0.75]])\n", + "response_model_domain = np.array([[0.0, 1.0]])\n", + "\n", + "old_result = previous_equality_rule(\n", + " battery_operating_range,\n", + " response_model_domain,\n", + ")\n", + "\n", + "new_result = containment_rule(\n", + " battery_operating_range,\n", + " response_model_domain,\n", + ")\n", + "\n", + "print(\"Upstream operating range:\", battery_operating_range)\n", + "print(\"Downstream accepted domain:\", response_model_domain)\n", + "print()\n", + "print(\"Previous equality rule accepts the chain:\", old_result)\n", + "print(\"Containment rule accepts the chain:\", new_result)" + ] + }, + { + "cell_type": "markdown", + "id": "5fad26b1", + "metadata": {}, + "source": [ + "### What does this output mean?\n", + "\n", + "The two checks are answering different questions.\n", + "\n", + "The previous rule asks:\n", + "\n", + "> “Does the upstream stage produce **exactly the entire domain** expected by\n", + "> the downstream stage?”\n", + "\n", + "and returns False because\n", + "\n", + "$$\n", + "[0.25, 0.75] \\neq [0, 1].\n", + "$$\n", + "\n", + "The containment rule instead asks:\n", + "\n", + "> “Can every value produced upstream safely be passed downstream?”\n", + "\n", + "and returns True because\n", + "\n", + "$$\n", + "0 \\leq 0.25\n", + "\\qquad \\text{and} \\qquad\n", + "0.75 \\leq 1.\n", + "$$\n", + "\n", + "Therefore,\n", + "\n", + "$$\n", + "\\boxed{[0.25, 0.75] \\subseteq [0, 1]}.\n", + "$$\n", + "\n", + "This is the behavior the PR introduces." + ] + }, + { + "cell_type": "markdown", + "id": "6c017194", + "metadata": {}, + "source": [ + "## 4. Now use the actual QMCPy transformation chain\n", + "\n", + "We now represent the same situation using two `TrueMeasure`s.\n", + "\n", + "The first transformation produces normalized operating loads in\n", + "\n", + "$$\n", + "[0.25,0.75].\n", + "$$\n", + "\n", + "The second transformation is a `Kumaraswamy` transformation. Its input domain\n", + "is the unit interval\n", + "\n", + "$$\n", + "[0,1].\n", + "$$\n", + "\n", + "A Kumaraswamy distribution is useful for quantities naturally bounded between\n", + "0 and 1, such as proportions, normalized scores, or utilization fractions.\n", + "\n", + "For this example we use\n", + "\n", + "$$\n", + "a=2,\n", + "\\qquad\n", + "b=2.\n", + "$$\n", + "\n", + "The important point for this PR is not the particular choice of distribution.\n", + "It is that the downstream transformation accepts the whole unit interval while\n", + "the upstream transformation produces only a valid subset of it." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "fd51b128", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-05T15:28:36.966639Z", + "iopub.status.busy": "2026-09-05T15:28:36.966489Z", + "iopub.status.idle": "2026-09-05T15:28:36.970611Z", + "shell.execute_reply": "2026-09-05T15:28:36.970195Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Range produced by the inner transformation:\n", + "[[0.25 0.75]]\n", + "\n", + "Domain accepted by the outer transformation:\n", + "[[0 1]]\n", + "\n", + "QMCPy reports a compatibility error: False\n" + ] + } + ], + "source": [ + "inner_measure = Uniform(\n", + " DigitalNetB2(1, randomize=\"FALSE\"),\n", + " lower_bound=0.25,\n", + " upper_bound=0.75,\n", + ")\n", + "\n", + "battery_response = Kumaraswamy(\n", + " inner_measure,\n", + " a=2.0,\n", + " b=2.0,\n", + ")\n", + "\n", + "print(\"Range produced by the inner transformation:\")\n", + "print(inner_measure.range)\n", + "\n", + "print(\"\\nDomain accepted by the outer transformation:\")\n", + "print(battery_response.domain)\n", + "\n", + "print(\n", + " \"\\nQMCPy reports a compatibility error:\",\n", + " battery_response.sub_compatibility_error,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9e3e56de", + "metadata": {}, + "source": [ + "### Why is `Compatibility error: False` the important result?\n", + "\n", + "The inner transformation produces only\n", + "\n", + "$$\n", + "[0.25,0.75].\n", + "$$\n", + "\n", + "The outer transformation accepts\n", + "\n", + "$$\n", + "[0,1].\n", + "$$\n", + "\n", + "These sets are **not equal**, so the previous equality-based rule would have\n", + "marked this composition as incompatible.\n", + "\n", + "With containment,\n", + "\n", + "$$\n", + "[0.25,0.75]\\subseteq[0,1],\n", + "$$\n", + "\n", + "so QMCPy correctly allows the transformation chain.\n", + "\n", + "This is the direct behavioral change introduced by the PR." + ] + }, + { + "cell_type": "markdown", + "id": "02c17b09", + "metadata": {}, + "source": [ + "## 5. Generate values through the complete chain\n", + "\n", + "Accepting the chain is useful only if it can actually be evaluated.\n", + "\n", + "Let the original DigitalNet coordinate be\n", + "\n", + "$$\n", + "u\\in[0,1].\n", + "$$\n", + "\n", + "The inner `Uniform(0.25, 0.75)` transformation is\n", + "\n", + "$$\n", + "r(u)\n", + "=\n", + "0.25+(0.75-0.25)u\n", + "=\n", + "0.25+0.5u.\n", + "$$\n", + "\n", + "Therefore,\n", + "\n", + "$$\n", + "r(u)\\in[0.25,0.75].\n", + "$$\n", + "\n", + "For a Kumaraswamy distribution with parameters $a=b=2$, the CDF is\n", + "\n", + "$$\n", + "F(y)=1-(1-y^2)^2,\n", + "\\qquad 0\\le y\\le1.\n", + "$$\n", + "\n", + "Its inverse CDF is\n", + "\n", + "$$\n", + "F^{-1}(r)\n", + "=\n", + "\\left[\n", + "1-(1-r)^{1/2}\n", + "\\right]^{1/2}.\n", + "$$\n", + "\n", + "The complete chained transformation is therefore\n", + "\n", + "$$\n", + "\\boxed{\n", + "y(u)\n", + "=\n", + "\\sqrt{\n", + "1-\\sqrt{\n", + "1-(0.25+0.5u)\n", + "}\n", + "}\n", + "}.\n", + "$$\n", + "\n", + "The next cell compares this formula with the values produced by QMCPy." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "f6da6902", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-05T15:28:37.287457Z", + "iopub.status.busy": "2026-09-05T15:28:37.287316Z", + "iopub.status.idle": "2026-09-05T15:28:37.400058Z", + "shell.execute_reply": "2026-09-05T15:28:37.399472Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\Owner\\Downloads\\QMCSoftware\\qmcpy\\discrete_distribution\\digital_net_b2\\digital_net_b2.py:675: ParameterWarning: Without randomization, the first digtial net point is the origin\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
DigitalNet uOperating load r(u)Formula resultQMCPy result
00.0000.25000.3660250.366025
10.5000.50000.5411960.541196
20.2500.37500.4576360.457636
30.7500.62500.6225970.622597
40.1250.31250.4133330.413333
50.6250.56250.5818610.581861
60.3750.43750.5000000.500000
70.8750.68750.6640660.664066
\n", + "
" + ], + "text/plain": [ + " DigitalNet u Operating load r(u) Formula result QMCPy result\n", + "0 0.000 0.2500 0.366025 0.366025\n", + "1 0.500 0.5000 0.541196 0.541196\n", + "2 0.250 0.3750 0.457636 0.457636\n", + "3 0.750 0.6250 0.622597 0.622597\n", + "4 0.125 0.3125 0.413333 0.413333\n", + "5 0.625 0.5625 0.581861 0.581861\n", + "6 0.375 0.4375 0.500000 0.500000\n", + "7 0.875 0.6875 0.664066 0.664066" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sample_count = 8\n", + "\n", + "# Use the same deterministic DigitalNet construction for a manual calculation.\n", + "reference_u = DigitalNetB2(\n", + " 1,\n", + " randomize=\"FALSE\",\n", + ").gen_samples(sample_count)\n", + "\n", + "inner_values = 0.25 + 0.5 * reference_u\n", + "\n", + "expected_response = np.sqrt(\n", + " 1.0 - np.sqrt(1.0 - inner_values)\n", + ")\n", + "\n", + "qmcpy_response = battery_response.gen_samples(sample_count)\n", + "\n", + "comparison = pd.DataFrame(\n", + " {\n", + " \"DigitalNet u\": reference_u[:, 0],\n", + " \"Operating load r(u)\": inner_values[:, 0],\n", + " \"Formula result\": expected_response[:, 0],\n", + " \"QMCPy result\": qmcpy_response[:, 0],\n", + " }\n", + ")\n", + "\n", + "comparison" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e9045766", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Maximum difference between the mathematical formula and QMCPy: 0.000e+00\n", + "All generated values are finite: True\n" + ] + } + ], + "source": [ + "maximum_difference = np.max(\n", + " np.abs(qmcpy_response - expected_response)\n", + ")\n", + "\n", + "print(\n", + " \"Maximum difference between the mathematical formula and QMCPy:\",\n", + " f\"{maximum_difference:.3e}\",\n", + ")\n", + "\n", + "print(\n", + " \"All generated values are finite:\",\n", + " np.isfinite(qmcpy_response).all(),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "f8a0abb5", + "metadata": {}, + "source": [ + "### What did we just verify?\n", + "\n", + "Each row follows the complete path\n", + "\n", + "$$\n", + "u\n", + "\\longrightarrow\n", + "r(u)\n", + "\\longrightarrow\n", + "y(u).\n", + "$$\n", + "\n", + "The first transformation converts the DigitalNet coordinate into the restricted\n", + "operating range\n", + "\n", + "$$\n", + "r(u)\\in[0.25,0.75].\n", + "$$\n", + "\n", + "The second transformation then accepts that value because its domain is\n", + "\n", + "$$\n", + "[0,1].\n", + "$$\n", + "\n", + "The final comparison checks two independent calculations:\n", + "\n", + "1. the closed-form mathematical expression\n", + "\n", + " $$\n", + " y(u)\n", + " =\n", + " \\sqrt{\n", + " 1-\\sqrt{\n", + " 1-(0.25+0.5u)\n", + " }\n", + " },\n", + " $$\n", + "\n", + "2. the value returned by the actual QMCPy chained `TrueMeasure`.\n", + "\n", + "A maximum difference near machine precision means the actual QMCPy chain is\n", + "performing the same transformation predicted by the mathematics.\n", + "\n", + "So this is not only a compatibility-table change: the previously rejected\n", + "strict-subset chain can now be evaluated normally." + ] + }, + { + "cell_type": "markdown", + "id": "443bc49a", + "metadata": {}, + "source": [ + "## 6. What output range should we expect?\n", + "\n", + "Because\n", + "\n", + "$$\n", + "r\\in[0.25,0.75]\n", + "$$\n", + "\n", + "and the Kumaraswamy inverse CDF is increasing, the smallest and largest possible\n", + "outputs occur at the two endpoints.\n", + "\n", + "For $a=b=2$,\n", + "\n", + "$$\n", + "y_{\\min}\n", + "=\n", + "\\sqrt{1-\\sqrt{1-0.25}}\n", + "\\approx0.3660,\n", + "$$\n", + "\n", + "and\n", + "\n", + "$$\n", + "y_{\\max}\n", + "=\n", + "\\sqrt{1-\\sqrt{1-0.75}}\n", + "\\approx0.7071.\n", + "$$\n", + "\n", + "Therefore the composite transformation produces values approximately in\n", + "\n", + "$$\n", + "\\boxed{\n", + "[0.3660,\\;0.7071]\n", + "}.\n", + "$$\n", + "\n", + "We can check this directly with a larger point set." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "158b2803", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Theoretical composite output interval: [0.366025, 0.707107]\n", + "Observed sample minimum: 0.366025\n", + "Observed sample maximum: 0.707020\n", + "Every sample lies inside the theoretical interval: True\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\Owner\\Downloads\\QMCSoftware\\qmcpy\\discrete_distribution\\digital_net_b2\\digital_net_b2.py:675: ParameterWarning: Without randomization, the first digtial net point is the origin\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "many_samples = battery_response.gen_samples(2**12)\n", + "\n", + "theoretical_minimum = np.sqrt(\n", + " 1.0 - np.sqrt(1.0 - 0.25)\n", + ")\n", + "theoretical_maximum = np.sqrt(\n", + " 1.0 - np.sqrt(1.0 - 0.75)\n", + ")\n", + "\n", + "print(\n", + " \"Theoretical composite output interval:\",\n", + " f\"[{theoretical_minimum:.6f}, {theoretical_maximum:.6f}]\",\n", + ")\n", + "\n", + "print(\n", + " \"Observed sample minimum:\",\n", + " f\"{many_samples.min():.6f}\",\n", + ")\n", + "\n", + "print(\n", + " \"Observed sample maximum:\",\n", + " f\"{many_samples.max():.6f}\",\n", + ")\n", + "\n", + "print(\n", + " \"Every sample lies inside the theoretical interval:\",\n", + " bool(\n", + " np.all(many_samples >= theoretical_minimum)\n", + " and np.all(many_samples <= theoretical_maximum)\n", + " ),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6da698c3", + "metadata": {}, + "source": [ + "### Interpreting this output\n", + "\n", + "The upstream `Uniform` transformation deliberately uses only part of the\n", + "outer transformation's valid input domain.\n", + "\n", + "That restriction propagates through the second transformation and produces a\n", + "correspondingly restricted output interval.\n", + "\n", + "The key point is:\n", + "\n", + "$$\n", + "\\text{using only part of a valid domain}\n", + "\\neq\n", + "\\text{being incompatible with that domain}.\n", + "$$\n", + "\n", + "The downstream transformation does not require its caller to generate every\n", + "possible value in $[0,1]$.\n", + "\n", + "It requires only that every value it actually receives belongs to $[0,1]$." + ] + }, + { + "cell_type": "markdown", + "id": "cc9a4e48", + "metadata": {}, + "source": [ + "## 7. Visualizing why containment is enough\n", + "\n", + "The outer transformation accepts the entire interval\n", + "\n", + "$$\n", + "D=[0,1].\n", + "$$\n", + "\n", + "The inner transformation produces\n", + "\n", + "$$\n", + "R=[0.25,0.75].\n", + "$$\n", + "\n", + "The figure below shows the relationship directly." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "4ba66b23", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3kAAAEOCAYAAAA0SN6FAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAO75JREFUeJzt3Qd4FNXex/ETCBAQQhHpTVCQjsIrTUWUoqKADS54EcUOFkBFrCCgomLFiiLYkKYgVxRURFFARURFaVJFFLAgRHqZ9/kd7+yd3ewmm2w2yU6+n+dZZSezs2fOnJk9/zllkhzHcQwAAAAAwBcK5XUCAAAAAAA5hyAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgD/D4+OOPTVJSkpk+fbqv8qVWrVrmsssuy+tkANbGjRvteTZx4sRAjgwfPtwuy006J3RuxMvpp59uX9Fed/T/eFN6GjVqFPfviTVPkD/oHFXZ1DmbX+Xm+QMkEoI8+J4u/tG8Ev0HYtGiRbai/Ndff5mC5t1337X7DsCYX375xZ4P33zzDdmRg8hXAIkkOa8TAMTbq6++GvT+lVdeMR988EG65fXr1zcrV65M6CDv3nvvta0TZcqUCfrb6tWrTaFChXwd5D399NMEegnsrrvuMkOHDs3V73zhhRfMkSNH4rb9999/3+RVMKJrgVopmzVrlidp8CPyNX867bTTzN69e03RokXzOilAvkKQB9/797//HfT+888/t0Fe6HJJ5CAvI8WKFct0nd27d5ujjjrK+N2hQ4dsxT6RKwR79uwxJUqUMH6SnJxsX7mpSJEicd1+IpcxxM6P52l+pBuYKSkpeZ0MIN/x7619IAYKAu677z5TrVo1++Nx5plnmrVr16Zb74svvjBnnXWWKV26tP0xb9eunVm4cGFU37F//34zbNgwc9xxx9kgrHr16mbIkCF2uZe6kl5//fVm5syZdiyN1m3YsKGZM2dOYB11zbr11lvtv4899thAF1R3HEXomDx3nMUnn3xi+vfvbypUqGD31fXee++ZU0891QZ9pUqVMl26dDE//PBDpvsUaVxVuHEdStO5555rWzvU2qB8btCggXnrrbeCPnvw4EHbKnH88cfbdY4++mhzyimn2EBdtF9qxXPzyn15x36NGTPGPP7446ZOnTo2/1asWGH/vmrVKnPRRReZcuXK2W23aNHCzJo1K+j7//zzT3PLLbeYxo0bm5IlS5rU1FRz9tlnm2+//TbsuJCpU6fa9FatWtXmnba/c+dOe1wHDhxo81rbufzyy9Md64zGUC1dutTesVY5u+OOO+zf3n77bXtsqlSpYvdL+zdy5Ehz+PDhsNvQfrdv395uQ+l76KGH0n3fpk2bTNeuXe2xV1oHDRpk5s6dG7ZLcyzlP9qy89prr5nmzZub4sWL2+P0r3/9y2zevDlonR9//NFceOGFplKlSvY4qixrPeV7VsbkecvLuHHjAuXl//7v/8ySJUuCPrt161Z7DPVdWqdy5cqmW7duQWU83Pizn3/+2XTv3j0ofyOVg+zkr46R0itKn3s+eMc/SjRlIdprVCRuHurYnXzyyebTTz8Nu9727dvNFVdcYSpWrGiPX9OmTc3LL78ctM5JJ51kLrjggqBlOie1b999911g2ZQpU+wy94adW650/XZ7OSg/lTcKwrx0TdG1RevoHK1Xr17gXMssXzM6T6PNxwkTJpgzzjjDlgutp+vhs88+my6/3Gun0qRrlvJXeeGen7qG6r3yUufOsmXLojpeusbr+7U9letRo0ZFbOl+5pln7O+Q0qnrz4ABA9INFXDzRMdHZVd5ojxwx7zr96dly5b2+5TXH374YbprkX6f9Deto2v/xRdfnG58YLgxeVm55gF+RUseEMbo0aPt3UFV7lVR1A/DJZdcYitdro8++shW9vUjqh9wre/+SKsyo0pNJPrhVEX6s88+M1dffbXtKrp8+XLz2GOPmTVr1tiAzkvr6YdbP3gKHJ588klbqf3pp5/sD58qP/rcG2+8YbdRvnx5+7ljjjkmw+Or7Wmde+65x7bkibqx9u3b13Tu3Nk8+OCDtiKkioYqP6os5OREFaqc9+zZ01x77bX2O5V/+hFXANuxY8dAJe2BBx4wV155pc3TXbt2ma+++sp8/fXXdp1rrrnGdqMK1wXXpe3u27fP5rUqJQoWVKFp27at/eFXN0FVuhWgqQL+5ptvmvPPP99+dv369fZ4KF0KoLdt22aef/55W2lRBUIVHC+lVRUSbVMVy7Fjx9oWI5WPHTt22P1Ra7Iqh9qe8j4zf/zxhy1rClzUAq3KsGgbqowOHjzY/l9lUttTHj388MNB29B3K2BQWenRo4etaN122222Mqhti8qAyu+vv/5qbrrpJhs0TZo0ycyfPz9dmmIp/9HSjZa7777bplfH/7fffrP5qUq0yqIq4wcOHLBlVRXmG264waZ5y5Yt5p133rGVTlXos0r7nJaWZsuWKo86/5VvKgtu65/OP5UhfafOCQUqKoM6JyOdI+pSphtGWufGG2+0ZUdlVnmZU/mra8mIESNsOVB5180aadOmTZbKQlavUaHGjx9v80/fq5sbyjttT+eeghxvnqhCrnNFN7N0TkybNs0GZDp+Koei/dD1zXvzRfmvfFF+NGnSxC7Xv3VNU3q9tJ/ats5PXTtefPFFG0zpGifalgInbUf5p+uE0uQG1dHka7jzNCv5qOusAietr1bt//znP/YarW0oiPJS2nr37m3zWN+lGxPnnXeeee6552xwqc+J9lf7nlmXfd20UDCkng7u9VBBuq5loXQN042sDh06mOuuu85uW2nXjRDll7eFXGVN+ao80TVU6+nfr7/+ui0XuvZrP3S90g0x3cDRb5xoexqGoPUVdCq40+dVXnTtzayVNJpyDviaAxQwAwYMcCIV/fnz59u/1a9f39m/f39g+RNPPGGXL1++3L4/cuSIc/zxxzudO3e2/3bt2bPHOfbYY52OHTtmmIZXX33VKVSokPPpp58GLX/uuefs9yxcuDCwTO+LFi3qrF27NrDs22+/tcvHjh0bWPbwww/bZRs2bEj3fTVr1nT69u0beD9hwgS77imnnOIcOnQosDwtLc0pU6aMc9VVVwV9fuvWrU7p0qXTLQ81bNiwsHnrfp83bUqTlr355puBZTt37nQqV67snHjiiYFlTZs2dbp06ZKtY6rv0/LU1FRn+/btQX8788wzncaNGzv79u0LLNOxbNOmjT22Lv398OHD6bZbrFgxZ8SIEenKTqNGjZwDBw4Elvfq1ctJSkpyzj777KBttG7d2uZBZtq1a2e3q7IRSuUt1DXXXOOUKFEiaL/cbbzyyiuBZSrflSpVci688MLAskceecSuN3PmzMCyvXv3OieccIJdrn3MifLvHheVi0hlZ+PGjU7hwoWd++67L+izOgeTk5MDy5ctW2Y/N23aNCerdE54j4GbrqOPPtr5888/A8vffvttu/w///mPfb9jxw77XudcRpTverkef/xx+7mpU6cGlu3evds57rjjcjR/lyxZki5/s1oWsnKNCqXyX6FCBadZs2ZB19Fx48bZz4bLk9deey3o8zo/SpYs6ezatcsu0/HVeitWrLDvZ82aZc/Brl27Oj179gx8tkmTJs7555+frlz169cvKI1aR8fZ9dhjj9n1fvvtt5jyNfQ8zUo+hjufVQZq164dtMy9di5atCiwbO7cuXZZ8eLFnU2bNgWWP//880FlK5KBAwfa9b744ovAMl0zdd33Xru1TL9HnTp1CrouPvXUU3a9l156KV2eTJo0KbBs1apVdpny5PPPP0+Xfm/ehsuPxYsXpyu/7rXXu4/RlnPAz+iuCYSh7jje8TTuXVvdjRbNWqdWKN2B1N3b33//3b7UEqI79QsWLMhwQgfdqdYd3RNOOCHwWb10l15CW050x1Tdnly626xug256suuqq64yhQsXDrxXS4Tunvfq1SsoXVpH3WrCtejEQi0ZbouZaJ8uvfRS20qjO8ui1hrdZVd+Z5daXbytmmoFUEuJ7u6qxcbdTx1LtQrpu9QaJLqj794BVzdIreN25VKLQCil33snW/mmWL1fv35B62m57lrrznlmlAaVyVDeu+zufqisqvVVXVG9lGbvOFSVb7UGecuQWlDVsqmWBJe6fKmceMVa/qOhlmttQ8fIWxbVUqeuu25ZdFvq1KU0tPtddql1uWzZshHPf+W78k/dw9RakJUJgtStUy0WLrVGqIUnN/M3mrKQ1WuUl1ra1bKpVhrvdVStc6Etq8oTHVNdc1w6f9TS+ffff9sufd5joH13W+zUfVKt+W43UF27vv/++8C6XkqLl9ZR3qrVW9zJqtQFOrt5G+48zUo+es9n9SDReuoxoOMS2vVYXTlbt24ddD0RbbdGjRrplmf2W6Hj0KpVq6AWYl0z1YPFS10q1XquVjhvy6CuEbp+z549O11ZU0ucS9dN5bXyxE1bpHR680Pd9nW81N1Tnw937c1OOQf8jO6aQBjeH0lxK3xuhc4NONTFMBL9KHsril76vMaMROpOqQpSRulx05SVCmY46r4Umi5xKyCh9COek/SDHToOq27duvb/6pqjyp+6SGmsk5ZrjIW63/Tp0yfQPSsaofuprk4KvNQVUK9Ix0ABjyp8TzzxhB2DsmHDhqDxbuoqGyr0WLmVWm8XNXe5tq1yEm47XkpHuEk8FPxqVkoFrG5l1RVaKVR3p9C8VhnyjmfSGBjdTAhdT8fJK9byHw19h46RArpw3EBax1bdVR999FHbBUyVdwWpqtxlp6tmNOe/KvPq5nfzzTfbLnmqHKtLmgJ8ldlIlL/hyrwqvqH7Hs/8jaYsZPUaFbqfEnrsdMxq166dbl2tF9qV0O1u6W5L+az1FNCpi6L+r+6F6rqrLrOquCu9OqfCBXkZHVNd1xTYqwunugWru6KCaXXzU0Ae7czE4c7TrOSjujqqa+7ixYvT3bDQ8faW56xcZ9z9zIjy2Rt0RSqb7vEIXa791rF1/55RWVOaokmnuvKqu6m6Keum2z8dW/6R2XjbaMs54GcEeUAY3tYtL/dHxr3Tq3EEkaYo113ESPR5jQtQxTSc0B/AzNKTXaHjLdz90jihcJXVzGY/jPQw69CJQLJClbh169bZO+yapEUVMY1n0dgTVchi2U+NuVTLXThuYHP//ffbQFAtcZrURGOKVOnTnexwd/wjHatYjmG4cTFqtdBdflVQFQgrOFOrm+5wa9xJaNpysgzFWv6j/Q6VJ00CFC7t3u0/8sgjtpXILSNqBVLlUGMfvRMKRSuavNLx1xgojalSK6LKiL5TAfeJJ55oYhHv/I1m/7J6jcoNGhc8b948W/nXBCcaH6cbP2rZUdCnYEr5Ei7/M9tnnWNqJVTLmlqj1KqtSVx0w0tlKtLnMztPo81HXeMUWKrFT+tquQIntbDpehft+Ryv34rsiiWdCt4V4OlcU6ulAkFdE9QyGE1ra37LCyC3EeQB2eB2nVQFW10ps/N5zc6oH/VIgVFW5cR23P3ShATZ2S/37rgCEO+z+kLv7oa2qHnTrskIxDt5hQIrdYPSS124FPhp8L8b5GV1393WBLUsZLafGqyvFgNNJOGlfXQnuMkL6iqo7kvq1qj8cKm1Mbtq1qxpJzQIPSahM8vGWv6joe9QOtRS57buZkQVab3UsqnJGjSpjm4EaIbAeFEa1Zqnl1psFJAp4NSMoJHyV90JQ/NXE1eEbjeW/M2pa0F2r1HaT1GeeHsFqMudyqdmz/Suq5YVVdq9LWZud2N3W6IWOlX6J0+ebG8cadITfUbBnxvkaVk0AVk42pb2Vy8FWrrBc+edd9rAT8chO/kabT5qkhVNHqTZfb2tdDndRT4S5XO4LvGhZdM9HlrubZVVF04d25y8Hujaq9ZsnVMuTaAVOosngPAYkwdkg2a804+3ZjRT0BFKswBmROOM1P1ED2MOpbvU7kyXWeE+4y6WH0C1aqliqcqNKmRZ3S+3cuqOmxHtS+h06C7NijljxozAe3U51MPqVVl2WxIVyHjpTr1a2bzTj2d13xXEaoY2zZKpmSQz2k9VGEPv/GqcjTtmL6+4FVlv2lTRUrfSWI6/9sv7GAlVqkLLaazlPxrqKqd91Cx+ofmv9265UJkJHdeoYE8V9min+s8qdaVTvngpPzQrYEbfec4559gy704h725LsxjmZP7mxLUglmuUpvVX90QF2SqTLs0GG5om5YnG36rVzKXjqVlUda6rtdrldsNUV1l113a7+Gm5Wvg0FjBcV81oaJxuKLcV1T2m2cnXaPMx3PmsLokKanODjoNavr/88sugcqYu0F4K4tTCqBmevWnVTTClV490ySnhrr0qF7H0DAEKElrygGxQBVLdBjUNs6a8VguTxmPox1x3XhUo6c5sJBpTpun6NRmA1lerg364dPday9X9SxWlrFDFUHTnWd1Z1Eql7mRZecC50q0pqpU+PZdK21FlTVO+qwuT0vnUU09F/HynTp3sXWg980rP7dOP9EsvvRTYRii10GhdTZWtMTdaV48o8FZsNMGAAjLtn1r0VJFTJVnTrYfuu7rpKVDR93oH+4ejZ+upBUABgSYN0F1pfbfGw+hZZu5z8DTWSt0hdYzVSqDpz1XxCR1blNuUFrWc6k639lutBOpmG0tXJI110vHVJBiaul6ThGhf3QcNuy0RsZb/aCjIUSvc7bffbsdn6tEWCqLUWqAbA5qsRN1t1T1SZUHTs6s8KUBQPqgMaMKdeFBrs1pmVIFX+VQ3ZqVJ5SejcqdypvzV2D11N1T+Kq2hU8HHmr/KO7WkK8hSnukaoPFWoWNTMxLLNUrXHh07lSe15Gm8m46bzuvQ80bHUTdb1N1WeaIWfJ3fGp+mZ1u60+mLbu7o5o9akdSVz6WWbHVRluwGeTrHdXNKQYpaqzRWTjdM1N1X14ns5mu0+ahrp4InXbOVbwruFRjqhlS4G1E5Tc/tU1nUmGed++4jFNyWVpeu5TondfNF62r8q46H8koT4XgnOomVrr1Kk4J5nWe6Nmvil8zGMAP4r7ye3hPIj49QCJ2OPdyU7+707RdccIGdilvTeWtq6x49ejjz5s3LNB2aJvzBBx90GjZsaD9btmxZp3nz5s69995rHyXg0vcqzZk9FkFGjhzpVK1a1U5P7Z32OtIjFDQleKR80NTdmj47JSXFqVOnjnPZZZc5X331Vab7tXTpUqdly5Z2mu0aNWo4jz76aMRHKOjRCJo6W9OeKw80VX9o3o8aNco5+eST7aMdND241tH0+d7HFOgxEDfccINzzDHH2McVuMfXPW6Rprpft26dc+mll9pptYsUKWLz7txzz3WmT58eWEePIrj55pvtox30/W3btrXTeIdOjx+p7ETKa3dq94ymbBd9h8pIOJp+vVWrVjZdVapUcYYMGRKYijx0OvFw2wh9hICsX7/eHhdtU/mpfddjLrRN75TnsZT/aB6h4NJ361EfRx11lH3p+Ot8WL16dSC9mh5fZVRltVy5ck779u2dDz/8MMM0hNv/jMqLliuN8vvvv9s0KC1Kk84TlXnvoxEktIyIprfXtP96zEX58uWdm266yZkzZ07Yae5jub7osQ8NGjSwj5vw5nVWykK016hInnnmGfvIB322RYsWzoIFC8LmybZt25zLL7/c5oeuG3q0SbjHFMjFF19s92fKlClB6VR+6rN65Ec051noNUl52q1bN3seaTv6vx5/smbNmpjyNSv5qMdC6FqoclyrVi37GT2SINK1M1S434rMroFe3333nd0Pfb+uhfo9GT9+fNhH8+iRCSr/um5WrFjRue666+yjRbwi5Um06df23HKhx2noN0mPYAj9PYv0CIVoyzngV0n6jxvwAUBu0R17TZqgh1Yjf1OLyqBBg2wLp1qUAABA/saYPABA0DghL409U3c6TV9PgAcAQGJgTB4AIGjCE42r1KQTmkhBM0Vq/FDoBAwAACD/IsgDAARo4hpN+qGgThNEaMIDTVmvyTMAAEBiYEweAAAAAPgIY/IAAAAAwEcI8gAAAACgoI3JO3LkiPnll1/swz/dh+ECAAAAAHKPnn6XlpZmqlSpYgoVKhRbkKcAr3r16jmZPgAAAABANmzevNlUq1YttiBPLXjuxlJTU01+jGg11Xfp0qVpaQTlDr7G9Q6UOxQUXO9AuUtv165dtvHNjc9iCvLcLpoK8PJrkKeX0kZ3UlDu4Gdc70C5Q0HB9Q6Uu8gyi3mYeAUAAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfCTZJBjHcczf+w8F3q/ZlmY+XLHNHNq72yQXP8p0aFDR1K1YKk/TiILBLYuF9h00SUlJeZ0cFBCUO1DuUFBwvUNuW5NJXFGyWHLC1PmSHJ1Bmdi1a5cpXbq02blzp0lNTTV5KW3fQdN4+Pt5mgYAAAAABcvy4Z1MqZQieZqGaOMyumsCAAAAgI8Q5AEAAACAjxRKxL6yAAAAAJCb1iRQHJJwQd6HK7fldRIAAAAAFDDzVm43iSLhgry0vYfzOgkAAAAACphde/83w39+l3BBXqnihfM6CQAAAAAKmNTiifP0uYQL8jrUr5jXSQAAAABQwHRIoDgk4YI8HnQOAAAAILcdX7FkwmR6wgV5AAAAAIDICPIAAAAAwEcSZ/Tgf5UslmyWD+8U9LyKeSu2mYN7d5sixY8yHRpUSqimVCQux3HMrl27TGpqqklKSsrr5KCAoNyBcoeCgusdctuaTOIKxSGJInFS+l+qTJdKKRJ437xmOXNSjbJm586dpnTp0lS2kas/Pkf2J9vySJAHyh38jOsdKHcoCJr7KK6guyYAAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+AhBHgAAAAD4CEEeAAAAAPgIQR4AAAAA+EhyXicAQIwcx5j9af97/9sqY1a/Z8y+XcakpBpT72xjjjmBbPaN/x7v/UnGGL0Ayh2yKbPfi2KljEniOgMkIoI8INGpwj+6euS/f/ZobqYGcabqVhlyGbmMcldAhP5eDN38T/AHIOHQXRMAAAAAfIQgDwAAAAB8hCAP8MOYCgAA+H0B8F8EeUCi06B5AABy/PdlDnkKJCiCPCDRaVY0AABy/PdlJ3kKJCiCPCDRMfMZACAuvy+lyVcgQRHkAYlOzzUCAIDfFwD/RZAHJDoedA4AiMvvSz3yFUhQBHkAAAAA4CMEeQAAAADgI8l5nQAAMSpWypihm4Ofm6dprzUrmgbNa8weXW58wzGO2blzlyldOtUkmaS8Tg4KCMqdT2X2e6HfFwAJiSAPSHRJScEzbFY/+Z8X/MlxjCmmV+o/xx6g3CG7+L0AfIvumgAAAADgIwR5AAAAAOAjBHkAAAAA4CMEeQAAAADgIwR5AAAAAOAjBHkAAAAA4CMEeQAAAADgIwR5AAAAAOAjBHkAAAAA4CMEeQAAAADgI8l5nQDADxzHMYcOHTKHDx/O66SgAJS1AwcOmH379pmkpKS8Tg4KCModclPhwoVNcjJVVCAWnEFAjFTh3rp1q9mzZw95iVxx5MgR88cff5DbyFWUO+SmEiVKmEqVKpHpQDYR5AEx3t3euHGjvetYpUoVU7RoUVpXEPcypxZjlTla8pBbKHfIzbKmm6e//fab/X2tUKECmQ9kA0EeEIODBw/au9vVq1e3dx2BeKOyjbxAuUNuKl68uClSpIjZtGmTHQoBIOuYeAXIAYUKcSoBAJDTv6u6wQAg66iZAgAAAICPEOQBAAAAgI8Q5AEAAACAjxDkATHSeIG/9x0yafsO5suXn8cznH766WbgwIEmP8hPaYk2vZqdU69vvvkmV7/7sssuC3z3zJkz4/dFKvv7duXfl4/PzYIsL8+tRJFr1wCgAGN2TSBGew46psXoT/JtPi4f3smUSikS8e+bN282w4YNM3PmzDG///67qVy5sunevbu55557zNFHH52lik2zZs3M448/bgqit956y84GF6t27dqZBQsW2H/rMQnHHHOMad26tRkyZIhp1aqVyUlXXXWVGTFihClfvnzQ8qeffto8/PDD9vmPTZs2NWPHjjUnn3xyVNtU2vXZpUuXml9//dXMmDHDlievJ554wowePdqWtbjan2bM6Oom3xq62ZiU1IiV4L/++osKcA7znl86X+vUqWPuvvtu07t373x3bj3wwAP2urJq1So722SbNm3Mgw8+aOrVqxdYZ/jw4ebee+8N+pz+rs9EKytpqlWrlp3xMlT//v3tdqJNU65dA4ACjJY8oABbv369adGihfnxxx/NG2+8YdauXWuee+45M2/ePBtY/Pnnn7meJj0fKRGVK1fOlCpVKqZtqNV12bJltvKjAEnHY9KkSXaWuVNPPdV88MEHJh4PG05O/t/9vilTppjBgwfbwP/rr7+2lb7OnTub7du3R7XN3bt328+4Fb5wSpcuzUOOE0S052MinLfu+TVmzBh7fq1evdqcddZZ5tJLLzUbNmzId+fWJ598YgYMGGA+//xze+7rkT2dOnWy55hXw4YN7f64r88++yzqdGY1TUuWLAn6LveadPHFF2cpTVwDgPgjyAMKMFUg9AD3999/397hrlGjhjn77LPNhx9+aLZs2WLuvPPOwN3b0BY6tdrpjq3b6qAKie7Oul1w9BBb0XMEdUf62GOPtXejVYmYPn16UAvg9ddfb7s66q63KhjhqGKjyljJkiXt3d9HHnkk3Tr79+83N954o314bkpKijnllFNspcRL33fDDTfY7ytbtqypWLGieeGFF+z2L7/8chuoHXfccea9994LfEatnNpWmTJlbOvmueeea9atWxexu6b+rXSo9U3Bnyp7bl5lRMF2WlqaOe200+xnlO/t27c306ZNs0G3ezzi6dFHH7WtEMqLBg0a2KBfFdaXXnopqs+r/IwaNcqcf/75cU9rQZJZmYqmzGV2LmblfIy0XmbnSjTp1DlwySWXmKOOOsqe64899li67tDR7Euk80uBnb5Xn73iiivM4cOHbcCX384t5aWurQqYtH8TJ040P/30k20l91Igqf1xX6GthzmZJvUs8H7XO++8Y1tD9fuRU2kCkDMI8oACSq10c+fOtd1sVEny0o+yKlm6yxvNmD4FdwpCVFlw79zqAfGiitgrr7xiKw8//PCDGTRokPn3v/9tg0LXyy+/bIPNhQsX2vXCufXWW+1n3n77bRuUfvzxx/bOs5cqjm+++abdnv6mYE2Vz9AWSf1dlY4vv/zSBnzXXXedvROt7lD6nO6W9+nTx+zZs8eurwBQd7u/+uor28qpljUFMapoRqLvUCX1iy++MA899JDtupVZS5wqb+qiqQqdl4Lmjh07hh3fc//999vAN6OXKobRtsYoDR06dAgs077q/eLFi6PaBuInszKV2d+jORejPR8jrRfNuZJZOvV5bXPWrFl2+aeffpruXI92X7xUtnVjR8GM/Pzzz/bGSbFixUyTJk3y/bm1c+dO+38Fx6HBa5UqVUzt2rXtdTu30qTPv/baa6Zfv372GpUTaQKQcxiTBxRQ+hFWAFe/fv2wf9fyHTt2mN9++y3TbanrjSp7bhclb8uaKkpqGVQQKPrRV9ed559/PnD39/jjj7eVvUj+/vtvM378eFuhOPPMMwMVxWrVqgXWUeXy2WeftXe71ZokaqFTJVGfVZDoUhB111132X/ffvvttnukgj4FqaLxiNrWd999Z8fBXXjhhUHp0V1u3dFesWKFadSoUdg0q9KoLlDu/j311FO20qtgLRJVZDV2RfkYSvkbbszftddea3r06GEyospWNDQmU60aat300vusjPFBfGRWpjL6e7TnYjTnY0brRXOuZJROtbTp3FY3ZfdcnzBhQlAZzsq+hJ5fCpTUWq9yvm/fPnuDS4FiuHMkP51bCpLVktm2bduga07Lli3tNU/XDd1c01g4de3+/vvvM+0+HmuaNGGKxo2qtdErljQByDkEeUABF8/ZNzWmTK1hoYGN7gCfeOKJgffNmzfPcDvq7qXPqPLg0t1s7wQEWkdjVlQJciko0gQCK1euDNqe9669Ws7Uraxx48aBZW6lxx2XooBYgZ9aHlQxclsldHc6oyDPS93OMhvXpkpopLxQd7ITTjgh3XLlQ+idffhTZmUqo79Hey5Gcz5mtF4050pG6dQ4YZ3H3ok/dBPJe65nZV9Czy91UVd3UQUnt9xyi71ehAYp+fHcUroVJIWObXNvaLn5qmtkzZo1zdSpU21X1HjSzTN9f2igm5dpAvA/BHlAAaWujOpiowAo3PgpLVfXJt2FVxee0GBQFbHMqAVOZs+ebapWrRr0N3WRcqnrVm4KbRFTPniXuV2P3ArqeeedZyspahlUhUbLVWHNaLKJcN+RUfdOtxKq7wqlVkp1XVM3tlBq0dArI2pF0XjLzKg1U0Hvtm3bgpbrvbeFFnkjszKV0d+jPRezcj6GWy+acyU754ZXVvYl9PxSa72uffLMM8/YIETLNP41v55bGvuosW+aFdTbeyEcjYWsW7euDYTjmSbNsKmWVM3+mZmspAlAzmFMHlBAqfVKd8JV0dm7d2/Q3zSV9uuvv2569uxpK2AK9NTtxrVr1650s9GpO6G6/nhp7IsqXbqLr4qV9+WO2YuGBvarYqjWAZe6kq5ZsyZoHXd8kDcQ1cQr7hic7Pjjjz9sK5q6d6r7mNuNNaepBUOtCyeddFLQcuWpuo2pW5nu5ofS3zRWL6NXtF3KlH9qnVHXOZcq3+5sq0hcOXUuxvtcUbdLneveCZPUxdJ7rmdnX9zzy9vyru3ouqGuoeHk9bmlG2sK8PQYko8++shOFBNNAKxeDdE8miCW811daDXBVZcuXXI0TQByDi15QAGmsTCabESTk2hGRFUiNImBxq/pDvl9991n1zvjjDPsGAvdpdddWXXH0h1gL90JVxCmWTU1IYH7SAF1idKkCKo8aNY9VdgUiKWmppq+fftGlU5tT918lC4Fp6pcaMIEtTB6WxU0gYrW0Xfr7rrGC6lbVyxdhNSaqe8cN26craSoYjl06FCT09wZ8/QdCrIVSGuZJrXRswx1J1/d1kJbVHO6S5laC3Vc9GgNdZnTrKruzKPRUIXOe8deNwNUGXaPCfJGTp2L8T5XlE6lxT2Pda5r/J7OdbeFPTv7onNJwaNalLwUjCqIuuOOO9J9JrfPLV2PlRY36NJNHQWgmmxK+6zrgug64E6WpXxwW09/+eUXm1e6Nvfq1StH0hQuXcpzBXn6nPcREa5Y0wQgZxDkAQWYJj3QLHj6EdYEA5qFUt109PBqLXMrOJqcRJV1TYeuCsbIkSPTteTph10/+ro7rpZB/V2Bn9ZVS6Bmw9PddAWJaq0KV6nKiB7WqwBClQdVeG6++ebAbHMuTaCiCohmxtQEDqq4aAZRVT6zS5XLyZMn23E8agXQ2KAnn3zSTumek9zZA1UJVYVI+aTv6tq1q21RyK2xQWq91WQ7CuRVqdSjMjSVu3dyBgX8qgSGG8+p8qTHPrjcLqYqG/oc8k5OnYvxPlc0rb/KvK43Cto0a65udOixKNndF51fut6p9cpLM0lq4hXNtJlZV8h4n1saw+h93IQmf5LQ/FOA5Y4jVLoVPKkVVfmhgFfP1dO/ozlfoznfQ9OlbpoK4DWrZjjRpAlA/CU5Ucy6oDvKqtipQqULbn6jXVDalMbQaXyBeJY7/Tj+9Otv5qI3/nkmXH60fHgnUyol/ayMSNxypy6cCgRjud6p4qgKXejzD6OhGwCaql6PsYiF0q8WAt1UiIt9u4wZnTNdEeNi6GZjUvLfb2o8y112qGVJPQv0bMxEmLgjlnMrHnLqfI2HjK4Bmv1UNwvVOqwAkfodcouTz+OKaOMyWvKAGJUokmS+GtrOFEuJPOA/L5UsxmmO8DQe88UXX7TPxPLOLpoZPSheXbiyS600ehxG3BUr9U8glV8pfUhn2bJldgp/dR9UJUbP0ZNu3br5/tyKh1jP13jItWsAUIDRkgfE2JKnLikay+btSgTk9xaVLVu2BCbc0Vi50G5s8aTp8nUnUjR2K7dnV0X+bslTkHfllVfaSVzcyUHUhTOvg6VEOLcSRTTXAFrykFccWvIAAIkqdOr53KTJNPQCwtGz7tyJiBJRXp5biYJrABB/PEIBAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA/IAVE8bhIAAPC7CuQKgjwgBsnJ/zyDbs+ePeQjAAA5xP1ddX9nAWQNZw4Qg0KFCpnSpUvbZ/5IiRIl4v4MKRRsufm8MoByh7y4xinA0++qfl/1Owsg6wjygBhVqlTJVrbdQA+ItyNHjlDxQa6j3CE3lSlTxlSsWDHw0HQAWUOQB8RIAV7lypXtw10PHjxIfiLud7nT0tJMqVKlaMlDrqHcITcVKVLE9lZgvDuQfQR5QA7RD5JeQDyp0rN//36TkpJCkIdcQ7kDgMRCR2cAAAAA8BGCPAAAAADwEYI8AAAAAChoY/Lcga/5dYYjpU9p0wQYTCkOyh38jOsdKHcoKLjegXKXnhuPZTYxUVRBnmZyk+rVq0ezOgAAAAAgThSf6VmSkSQ5UcxPq2fj/PLLL/l2ym5FtApAN2/ebFJTU/M6OSggKHeg3KGg4HoHyh0Kil35PK5wH2lTpUqVDJ+ZG1VLnjZQrVo1k9/pQOTHgwF/o9yBcoeCgusdKHcoKFLzcVyRUQuei4lXAAAAAMBHCPIAAAAAwEd8EeQVK1bMDBs2zP4foNzBz7jegXKHgoLrHSh32RfVxCsAAAAAgMTgi5Y8AAAAAMA/CPIAAAAAwEcI8gAAAADARxIiyHv66adNrVq1TEpKimnZsqX58ssvM1x/2rRp5oQTTrDrN27c2Lz77ru5llb4S1bK3gsvvGBOPfVUU7ZsWfvq0KFDpmUViLXceU2ePNkkJSWZ7t27k7GIe7n766+/zIABA0zlypXtBBl169bl9xZxL3ePP/64qVevnilevLh9YPWgQYPMvn37yHlEbcGCBea8886zDxPXb+bMmTMz/czHH39sTjrpJHutO+6448zEiRPzfY7n+yBvypQpZvDgwXb2zK+//to0bdrUdO7c2Wzfvj3s+osWLTK9evUyV1xxhVm2bJmt7Oj1/fff53rakdiyWvZ0AVDZmz9/vlm8eLH98enUqZPZsmVLrqcdBafcuTZu3GhuueUWe6MBiHe5O3DggOnYsaMtd9OnTzerV6+2N7qqVq1K5iNu5W7SpElm6NChdv2VK1ea8ePH223ccccd5Dqitnv3blvWdIMhGhs2bDBdunQx7du3N998840ZOHCgufLKK83cuXPzd647+dzJJ5/sDBgwIPD+8OHDTpUqVZwHHngg7Po9evRwunTpErSsZcuWzjXXXBP3tMJfslr2Qh06dMgpVaqU8/LLL8cxlfCb7JQ7lbU2bdo4L774otO3b1+nW7duuZRaFNRy9+yzzzq1a9d2Dhw4kIupREEvd1r3jDPOCFo2ePBgp23btnFPK/zJGOPMmDEjw3WGDBniNGzYMGhZz549nc6dOzv5Wb5uydOdwqVLl9pub65ChQrZ92opCUfLveuL7gpFWh/IqbIXas+ePebgwYOmXLlyZDLiWu5GjBhhKlSoYHswALlR7mbNmmVat25tu2tWrFjRNGrUyNx///3m8OHDHADErdy1adPGfsbt0rl+/XrbRficc84h1xE3ixM0tkg2+djvv/9ufzD0A+Kl96tWrQr7ma1bt4ZdX8uBeJa9ULfddpvt7x16YQBystx99tlntsuSupAAuVXuVLn+6KOPzCWXXGIr2WvXrjX9+/e3N7bUlQ6IR7nr3bu3/dwpp5yinmjm0KFD5tprr6W7JuJqa4TYYteuXWbv3r12fGh+lK9b8oBENXr0aDsJxowZM+xgciAe0tLSTJ8+fexYqPLly5PJyDVHjhyxrcfjxo0zzZs3Nz179jR33nmnee655zgKiBuNfVeL8TPPPGPH8L311ltm9uzZZuTIkeQ6kEgteaq0FC5c2Gzbti1oud5XqlQp7Ge0PCvrAzlV9lxjxoyxQd6HH35omjRpQgYjbuVu3bp1duILzRLmrXxLcnKynQyjTp06HAHkaLkTzahZpEgR+zlX/fr17R1vdcMrWrQouY4cL3d33323vbGlSS9EM6hrEo2rr77a3mRQd08gp1WKEFukpqbm21Y8yddng34kdIdw3rx5QRUYvddYgHC03Lu+fPDBBxHXB3Kq7MlDDz1k7yjOmTPHtGjRgsxFXMudHhWzfPly21XTfXXt2jUwA5hmeAVyutxJ27ZtbRdN96aCrFmzxgZ/BHiIx/XOHeseGsi5Nxr+mUMDyHmtEzW2cPK5yZMnO8WKFXMmTpzorFixwrn66qudMmXKOFu3brV/79OnjzN06NDA+gsXLnSSk5OdMWPGOCtXrnSGDRvmFClSxFm+fHke7gUSUVbL3ujRo52iRYs606dPd3799dfAKy0tLQ/3An4vd6GYXRO5Ue5++uknO3vw9ddf76xevdp55513nAoVKjijRo3iACBu5U51OpW7N954w1m/fr3z/vvvO3Xq1LEzqwPRSktLc5YtW2ZfCoUeffRR++9NmzbZv6vMqey5VNZKlCjh3HrrrTa2ePrpp53ChQs7c+bMydeZnu+DPBk7dqxTo0YNW4HWdLuff/554G/t2rWzlRqvqVOnOnXr1rXra8rT2bNn50Gq4QdZKXs1a9a0F4vQl36UgHiVu1AEecitcrdo0SL7iCJV0vU4hfvuu88+zgOIV7k7ePCgM3z4cBvYpaSkONWrV3f69+/v7Nixg0xH1ObPnx+2vuaWNf1fZS/0M82aNbPlVNe7CRMm5PscT9J/8ro1EQAAAABQAMbkAQAAAACyhiAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AAAAAfIQgDwAAAAB8hCAPAAAAAHyEIA8AkHA+/vhjk5SUZP766y/7fuLEiaZMmTJx/c7LLrvMdO/ePa7fMXz4cNOsWbO4fgcAwP8I8gCgAFPgomBp9OjRQctnzpxplyeKnj17mjVr1uR1MgAAyBcI8gCggEtJSTEPPvig2bFjR45u98CBAya3FC9e3FSoUCHXvg8AgPyMIA8ACrgOHTqYSpUqmQceeCDD9d58803TsGFDU6xYMVOrVi3zyCOPBP1dy0aOHGkuvfRSk5qaaq6++upAN8p33nnH1KtXz5QoUcJcdNFFZs+ePebll1+2nylbtqy58cYbzeHDhwPbevXVV02LFi1MqVKlbNp69+5ttm/fHjFtod01tV21RIa+XJs3bzY9evSwnylXrpzp1q2b2bhxY+DvSsvgwYPt348++mgzZMgQ4zhOxO/ftWuXDTTfe++9oOUzZsyw+6D9ldtuu83UrVvX5kPt2rXN3XffbQ4ePBhxu6effroZOHBg0DJ1GVULrGv//v3mlltuMVWrVjVHHXWUadmype3OCgAouAjyAKCAK1y4sLn//vvN2LFjzc8//xx2naVLl9qg6F//+pdZvny5HTumAEXBldeYMWNM06ZNzbJly+zfRQHOk08+aSZPnmzmzJljA5Dzzz/fvPvuu/algO75558306dPD2xHgY8Cxm+//dZ2HVUA5g1sMrNkyRLz66+/2pf2qVWrVubUU08NbLtz5842+Pr000/NwoULTcmSJc1ZZ50VaH1UAKt9e+mll8xnn31m/vzzTxuwRaKg9txzzzWTJk0KWv7666/boExBneg7td0VK1aYJ554wrzwwgvmscceM7G4/vrrzeLFi23+fvfdd+biiy+2+/Ljjz/GtF0AQAJzAAAFVt++fZ1u3brZf7dq1crp16+f/feMGTPUbBVYr3fv3k7Hjh2DPnvrrbc6DRo0CLyvWbOm071796B1JkyYYLezdu3awLJrrrnGKVGihJOWlhZY1rlzZ7s8kiVLltjtuJ+ZP3++fb9jx47A95QuXTrsZ2+88Uabtu3bt9v3r776qlOvXj3nyJEjgXX279/vFC9e3Jk7d659X7lyZeehhx4K/P3gwYNOtWrVAnkVjvKsZMmSzu7du+37nTt3OikpKc57770X8TMPP/yw07x588D7YcOGOU2bNg28b9eunXPTTTcFfUZp0HGTTZs2OYULF3a2bNkStM6ZZ57p3H777RG/FwDgb7TkAQAsjctTF8qVK1emyxEta9u2bdAyvVdrkbebpbpYhlIrVp06dQLvK1asaLtTqvXMu8zbHVMth+edd56pUaOGbf1q166dXf7TTz9l6WiNGzfOjB8/3syaNcscc8wxdplaB9euXWu3qzTopS6b+/btM+vWrTM7d+60LYDq9uhKTk4Ou29e55xzjilSpIj9Lrd7q1r41B3WNWXKFJtv6oKq773rrruyvE9ealVV/qsLqLsven3yySd2XwAABVNyXicAAJA/nHbaabYb4+23356lrpFeGhMWSoGPl8bGhVt25MgR++/du3fbdOil7o4KzhQI6X1WJnOZP3++ueGGG8wbb7xhmjRpElj+999/m+bNm9tth3IDwewoWrSoHW+oLpvq1qr/a9ZPBYiiLpWXXHKJuffee+2+lC5d2naxDB3b6FWoUKF0YwG9Y/i0L+puq6BY//fyBtEAgIKFIA8AEKBHKeg5bZokxat+/fp27JqX3qsFKTS4iNWqVavMH3/8YdNSvXp1u+yrr77K0jbUUqeA64477jAXXHBB0N9OOukk26Km2TjV0hZO5cqVzRdffGEDXzl06JANpPTZjCiI69ixo/nhhx/MRx99ZEaNGhX426JFi0zNmjXNnXfeGVi2adOmDLenoFOtii612n3//femffv29v2JJ55ol6kV1B1zCAAA3TUBAAGNGze2gYomSvG6+eabzbx58+xkKHoenbp1PvXUU3ZWx5ymLppqFdNEMOvXr7fdH/W90dq7d6/t6qkASDN8bt26NfAS7V/58uXtjJqaeGXDhg12MhjN8OlOPHPTTTfZIFOTvijo7N+/f+DB6xlRUKiumPqOY489NqjL5/HHH29bJNV6p66UyuOMJnORM844w8yePdu+lI7rrrsuKB0KsvVdmtH0rbfesvvy5Zdf2plS9RkAQMFEkAcACDJixIhA10mXWrCmTp1qA5RGjRqZe+65x66X3W6dmbVeaQbKadOmmQYNGthgS7N2Rmvbtm02IFJQWqVKFdsq577cMYILFiywwaRa+dRKecUVV9gxeW7LnoLaPn36mL59+5rWrVvb8XuaETQz6nbaq1cvO+5PwZdX165dzaBBg+xsmGotVcueOwNpJP369bNpUBCncYl67ILbiueaMGGC/bvSrBZYzeap2UW1fwCAgilJs6/kdSIAAAAAADmDljwAAAAA8BGCPAAAAADwEYI8AAAAAPARgjwAAAAA8BGCPAAAAADwEYI8AAAAAPARgjwAAAAA8BGCPAAAAADwEYI8AAAAAPARgjwAAAAA8BGCPAAAAADwEYI8AAAAADD+8f8MEm8WGQW7+wAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots(figsize=(9, 2.8))\n", + "\n", + "ax.plot(\n", + " [0, 1],\n", + " [1, 1],\n", + " linewidth=8,\n", + " label=r\"Outer domain $D=[0,1]$\",\n", + ")\n", + "\n", + "ax.plot(\n", + " [0.25, 0.75],\n", + " [0.65, 0.65],\n", + " linewidth=8,\n", + " label=r\"Inner range $R=[0.25,0.75]$\",\n", + ")\n", + "\n", + "ax.scatter(\n", + " [0, 1],\n", + " [1, 1],\n", + " s=60,\n", + ")\n", + "\n", + "ax.scatter(\n", + " [0.25, 0.75],\n", + " [0.65, 0.65],\n", + " s=60,\n", + ")\n", + "\n", + "ax.set_xlim(-0.05, 1.05)\n", + "ax.set_ylim(0.4, 1.25)\n", + "ax.set_yticks([])\n", + "ax.set_xlabel(\"Normalized value\")\n", + "ax.set_title(\"The entire upstream range lies inside the downstream domain\")\n", + "ax.legend(loc=\"lower center\", ncol=2)\n", + "ax.grid(axis=\"x\", alpha=0.2)\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e53d3628", + "metadata": {}, + "source": [ + "The shorter interval does not need to cover the longer interval.\n", + "\n", + "It only needs to fit inside it:\n", + "\n", + "$$\n", + "\\boxed{\n", + "[0.25,0.75]\\subseteq[0,1]\n", + "}.\n", + "$$\n", + "\n", + "That simple geometric relationship is the core of the compatibility change." + ] + }, + { + "cell_type": "markdown", + "id": "27194fb4", + "metadata": {}, + "source": [ + "## 8. The same idea in multiple dimensions\n", + "\n", + "The rule is applied coordinate by coordinate.\n", + "\n", + "Suppose a system has two normalized operating quantities:\n", + "\n", + "- coordinate 1 is restricted to\n", + "\n", + " $$\n", + " [0.10,0.80],\n", + " $$\n", + "\n", + "- coordinate 2 is restricted to\n", + "\n", + " $$\n", + " [0.20,0.90].\n", + " $$\n", + "\n", + "The next transformation accepts the unit interval in each coordinate:\n", + "\n", + "$$\n", + "D=[0,1]^2.\n", + "$$\n", + "\n", + "The upstream range is therefore\n", + "\n", + "$$\n", + "R=\n", + "[0.10,0.80]\n", + "\\times\n", + "[0.20,0.90].\n", + "$$\n", + "\n", + "Compatibility requires\n", + "\n", + "$$\n", + "[0.10,0.80]\\subseteq[0,1]\n", + "$$\n", + "\n", + "and\n", + "\n", + "$$\n", + "[0.20,0.90]\\subseteq[0,1].\n", + "$$\n", + "\n", + "Since both conditions hold,\n", + "\n", + "$$\n", + "\\boxed{R\\subseteq D}.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "ad1a30f8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Inner range:\n", + "[[0.1 0.8]\n", + " [0.2 0.9]]\n", + "\n", + "Outer domain:\n", + "[[0 1]]\n", + "\n", + "QMCPy reports a compatibility error: False\n", + "Generated sample shape: (8, 2)\n", + "All generated values are finite: True\n" + ] + } + ], + "source": [ + "inner_2d = Uniform(\n", + " DigitalNetB2(2, randomize=\"FALSE\"),\n", + " lower_bound=[0.10, 0.20],\n", + " upper_bound=[0.80, 0.90],\n", + ")\n", + "\n", + "outer_2d = Kumaraswamy(\n", + " inner_2d,\n", + " a=[2.0, 2.0],\n", + " b=[2.0, 2.0],\n", + ")\n", + "\n", + "samples_2d = outer_2d.gen_samples(8)\n", + "\n", + "print(\"Inner range:\")\n", + "print(inner_2d.range)\n", + "\n", + "print(\"\\nOuter domain:\")\n", + "print(outer_2d.domain)\n", + "\n", + "print(\n", + " \"\\nQMCPy reports a compatibility error:\",\n", + " outer_2d.sub_compatibility_error,\n", + ")\n", + "\n", + "print(\n", + " \"Generated sample shape:\",\n", + " samples_2d.shape,\n", + ")\n", + "\n", + "print(\n", + " \"All generated values are finite:\",\n", + " np.isfinite(samples_2d).all(),\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "448ea229", + "metadata": {}, + "source": [ + "### What does the two-dimensional result show?\n", + "\n", + "The inner range has one interval for each coordinate:\n", + "\n", + "$$\n", + "R=\n", + "\\begin{bmatrix}\n", + "0.10 & 0.80\\\\\n", + "0.20 & 0.90\n", + "\\end{bmatrix}.\n", + "$$\n", + "\n", + "The outer `Kumaraswamy` transformation uses the common unit-domain interval\n", + "\n", + "$$\n", + "D=[0,1].\n", + "$$\n", + "\n", + "QMCPy broadcasts that common domain across both coordinates and checks\n", + "\n", + "$$\n", + "0\\le0.10,\\qquad0.80\\le1,\n", + "$$\n", + "\n", + "and\n", + "\n", + "$$\n", + "0\\le0.20,\\qquad0.90\\le1.\n", + "$$\n", + "\n", + "Every coordinate passes, so the chain is compatible and samples can be\n", + "generated normally.\n", + "\n", + "This demonstrates that the change applies not only to a single interval but\n", + "also to multidimensional axis-aligned boxes." + ] + }, + { + "cell_type": "markdown", + "id": "c9a49861", + "metadata": {}, + "source": [ + "## 9. Containment does not mean “accept everything”\n", + "\n", + "The new rule is less restrictive than equality, but it is not weaker in the\n", + "sense of allowing invalid inputs.\n", + "\n", + "Consider an upstream transformation whose range is\n", + "\n", + "$$\n", + "[-0.10,0.75].\n", + "$$\n", + "\n", + "The next transformation still accepts only\n", + "\n", + "$$\n", + "[0,1].\n", + "$$\n", + "\n", + "Now\n", + "\n", + "$$\n", + "-0.10 < 0,\n", + "$$\n", + "\n", + "so part of the upstream range lies outside the downstream domain.\n", + "\n", + "Therefore,\n", + "\n", + "$$\n", + "[-0.10,0.75]\\nsubseteq[0,1].\n", + "$$\n", + "\n", + "This chain must still be rejected." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "be14463b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Invalid inner range:\n", + "[[-0.1 0.75]]\n", + "\n", + "Outer domain:\n", + "[[0 1]]\n", + "\n", + "QMCPy reports a compatibility error: True\n", + "\n", + "Sampling is correctly rejected:\n", + "ParameterError - The sub-transform range must be contained within the transform domain.\n" + ] + } + ], + "source": [ + "invalid_inner = Uniform(\n", + " DigitalNetB2(1, randomize=\"FALSE\"),\n", + " lower_bound=-0.10,\n", + " upper_bound=0.75,\n", + ")\n", + "\n", + "invalid_outer = Kumaraswamy(\n", + " invalid_inner,\n", + " a=2.0,\n", + " b=2.0,\n", + ")\n", + "\n", + "print(\"Invalid inner range:\")\n", + "print(invalid_inner.range)\n", + "\n", + "print(\"\\nOuter domain:\")\n", + "print(invalid_outer.domain)\n", + "\n", + "print(\n", + " \"\\nQMCPy reports a compatibility error:\",\n", + " invalid_outer.sub_compatibility_error,\n", + ")\n", + "\n", + "try:\n", + " invalid_outer.gen_samples(8)\n", + "except ParameterError as error:\n", + " print(\n", + " \"\\nSampling is correctly rejected:\"\n", + " )\n", + " print(type(error).__name__, \"-\", error)" + ] + }, + { + "cell_type": "markdown", + "id": "be62911d", + "metadata": {}, + "source": [ + "### Why is this rejection correct?\n", + "\n", + "This time the problem is not merely that the two intervals are different.\n", + "\n", + "The upstream transformation can actually produce values that the next\n", + "transformation does not accept.\n", + "\n", + "Specifically,\n", + "\n", + "$$\n", + "-0.10\\notin[0,1].\n", + "$$\n", + "\n", + "So the required containment condition fails:\n", + "\n", + "$$\n", + "R_{j-1}\\nsubseteq D_j.\n", + "$$\n", + "\n", + "The updated rule therefore distinguishes the two cases correctly:\n", + "\n", + "| Situation | Relationship | Result |\n", + "|---|---|---|\n", + "| Exact match | $R=D$ | Accepted |\n", + "| Smaller valid range | $R\\subset D$ | Accepted |\n", + "| Range leaves domain | $R\\nsubseteq D$ | Rejected |\n", + "\n", + "This is the safety property we want." + ] + }, + { + "cell_type": "markdown", + "id": "05374cca", + "metadata": {}, + "source": [ + "## 10. General form of the rule\n", + "\n", + "For a $d$-dimensional axis-aligned range,\n", + "\n", + "$$\n", + "R=\n", + "\\prod_{k=1}^{d}\n", + "[r_{k,L},r_{k,U}],\n", + "$$\n", + "\n", + "and domain\n", + "\n", + "$$\n", + "D=\n", + "\\prod_{k=1}^{d}\n", + "[d_{k,L},d_{k,U}],\n", + "$$\n", + "\n", + "the transformation chain is compatible exactly when\n", + "\n", + "$$\n", + "\\boxed{\n", + "d_{k,L}\\le r_{k,L}\n", + "\\quad\\text{and}\\quad\n", + "r_{k,U}\\le d_{k,U}\n", + "\\qquad\n", + "\\text{for every }k=1,\\ldots,d.\n", + "}\n", + "$$\n", + "\n", + "This includes ordinary finite intervals as well as numerical envelopes using\n", + "\n", + "$$\n", + "-\\infty\n", + "\\qquad\\text{or}\\qquad\n", + "+\\infty.\n", + "$$\n", + "\n", + "For example,\n", + "\n", + "$$\n", + "[-5,5]\\subseteq(-\\infty,\\infty).\n", + "$$\n", + "\n", + "QMCPy represents these using numeric bounds such as `-np.inf` and `np.inf`.\n", + "\n", + "This PR does not introduce representations for disconnected sets,\n", + "nonrectangular regions, or explicit open-versus-closed endpoints." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "ad8bd6c3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
CasePrevious equality ruleContainment rule
0Exact equalityTrueTrue
1Strict subsetFalseTrue
2Finite interval inside real lineFalseTrue
3Lower-bound violationFalseFalse
4Two-dimensional subsetFalseTrue
\n", + "
" + ], + "text/plain": [ + " Case Previous equality rule Containment rule\n", + "0 Exact equality True True\n", + "1 Strict subset False True\n", + "2 Finite interval inside real line False True\n", + "3 Lower-bound violation False False\n", + "4 Two-dimensional subset False True" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "examples = [\n", + " (\n", + " \"Exact equality\",\n", + " [[0.0, 1.0]],\n", + " [[0.0, 1.0]],\n", + " ),\n", + " (\n", + " \"Strict subset\",\n", + " [[0.25, 0.75]],\n", + " [[0.0, 1.0]],\n", + " ),\n", + " (\n", + " \"Finite interval inside real line\",\n", + " [[-5.0, 5.0]],\n", + " [[-np.inf, np.inf]],\n", + " ),\n", + " (\n", + " \"Lower-bound violation\",\n", + " [[-0.10, 0.75]],\n", + " [[0.0, 1.0]],\n", + " ),\n", + " (\n", + " \"Two-dimensional subset\",\n", + " [[0.10, 0.80], [0.20, 0.90]],\n", + " [[0.0, 1.0]],\n", + " ),\n", + "]\n", + "\n", + "summary = pd.DataFrame(\n", + " [\n", + " {\n", + " \"Case\": name,\n", + " \"Previous equality rule\": previous_equality_rule(\n", + " transform_range,\n", + " domain,\n", + " ),\n", + " \"Containment rule\": containment_rule(\n", + " transform_range,\n", + " domain,\n", + " ),\n", + " }\n", + " for name, transform_range, domain in examples\n", + " ]\n", + ")\n", + "\n", + "summary" + ] + }, + { + "cell_type": "markdown", + "id": "f591deeb", + "metadata": {}, + "source": [ + "## 11. Reading the final comparison\n", + "\n", + "The table summarizes exactly what changed.\n", + "\n", + "### Exact equality\n", + "\n", + "$$\n", + "[0,1]=[0,1].\n", + "$$\n", + "\n", + "Both rules accept it.\n", + "\n", + "### Strict containment\n", + "\n", + "$$\n", + "[0.25,0.75]\\subset[0,1].\n", + "$$\n", + "\n", + "The old equality rule rejects it, while the containment rule correctly accepts\n", + "it.\n", + "\n", + "### Unbounded containing domain\n", + "\n", + "$$\n", + "[-5,5]\\subset(-\\infty,\\infty).\n", + "$$\n", + "\n", + "Again, equality is unnecessary; containment is sufficient.\n", + "\n", + "### Genuine violation\n", + "\n", + "$$\n", + "[-0.10,0.75]\\nsubseteq[0,1].\n", + "$$\n", + "\n", + "Both the mathematics and QMCPy's updated behavior reject it.\n", + "\n", + "### Multiple dimensions\n", + "\n", + "Each coordinate is checked independently, so a rectangular range can be\n", + "contained in a larger rectangular domain even when their endpoint arrays are\n", + "not identical." + ] + }, + { + "cell_type": "markdown", + "id": "7f80cf0e", + "metadata": {}, + "source": [ + "## 12. Takeaway\n", + "\n", + "The change can be summarized in one line:\n", + "\n", + "$$\n", + "\\boxed{\n", + "R_{j-1}=D_j\n", + "\\quad\\longrightarrow\\quad\n", + "R_{j-1}\\subseteq D_j\n", + "}\n", + "$$\n", + "\n", + "but its practical meaning is important.\n", + "\n", + "A downstream transformation should not reject an upstream transformation\n", + "simply because the upstream stage uses only part of its valid input domain.\n", + "\n", + "In the battery example,\n", + "\n", + "$$\n", + "[0.25,0.75]\\subseteq[0,1],\n", + "$$\n", + "\n", + "so every upstream value is safe to pass to the next transformation.\n", + "\n", + "The previous equality-based rule would reject that composition.\n", + "\n", + "The updated containment rule:\n", + "\n", + "- accepts exact equality;\n", + "- accepts valid strict subsets;\n", + "- works coordinate-wise in multiple dimensions;\n", + "- supports finite and unbounded numeric envelopes;\n", + "- still rejects any range that actually leaves the next domain.\n", + "\n", + "The result is a more mathematically accurate compatibility check while keeping\n", + "the safety condition intact." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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 +} diff --git a/mkdocs.yml b/mkdocs.yml index 792731c47..35edf3f00 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -75,6 +75,7 @@ nav: - Importance Sampling with True Measures: - Statistics for True Measures: demos/statistics_for_TrueMeasure.ipynb - Some True Measures: demos/some_true_measures.ipynb + - TrueMeasure Domain Inclusion: demos/true_measure_domain_inclusion.ipynb - SciPyWrapper dependence and Custom distributions: demos/scipywrapper_dependence_custom/scipywrapper_demo.ipynb - ProductMeasure: demos/product_measure.ipynb - Acceptance-Rejection Sampling: demos/acceptance_rejection.ipynb diff --git a/qmcpy/discrete_distribution/dummy_sampler.py b/qmcpy/discrete_distribution/dummy_sampler.py index cf33720e5..cc8363cc9 100644 --- a/qmcpy/discrete_distribution/dummy_sampler.py +++ b/qmcpy/discrete_distribution/dummy_sampler.py @@ -30,6 +30,16 @@ class DummySampler(AbstractLDDiscreteDistribution): """ def __init__(self, dimension=1, replications=None, seed=None, warn=True): + r""" + Args: + dimension (Union[int, list, tuple, np.ndarray]): Dimension of the placeholder sampler. A list, tuple, or array specifies unique coordinate indices. Defaults to `1`. + replications (Union[None, int]): Replication metadata preserved when spawning placeholders. `None` records no explicit replication axis. Defaults to `None`. + seed (Union[None, int, np.random.SeedSequence]): Seed used to initialize the sampler state and spawn child samplers. Defaults to `None`. + warn (bool): Compatibility argument matching other discrete-distribution constructors. It is ignored because `DummySampler` cannot generate samples. Defaults to `True`. + + Raises: + ParameterError: If an array-like `dimension` is not one-dimensional with unique entries, if it exceeds the dimension limit, or if `replications` is negative. + """ # Keep the same constructor as other discrete distributions. del warn diff --git a/qmcpy/true_measure/abstract_true_measure.py b/qmcpy/true_measure/abstract_true_measure.py index b7611b1f3..abb0be649 100644 --- a/qmcpy/true_measure/abstract_true_measure.py +++ b/qmcpy/true_measure/abstract_true_measure.py @@ -6,6 +6,12 @@ from scipy import sparse +def _clip_unit_interval(u): + """Clip unit-interval values away from endpoints for stable quantiles.""" + eps = np.finfo(float).eps + return np.clip(u, eps, 1.0 - eps) + + class AbstractTrueMeasure(object): def __init__(self): @@ -13,12 +19,12 @@ def __init__(self): if not hasattr(self, "domain"): raise ParameterError( prefix - + "self.domain, 2xd ndarray of domain lower bounds (first col) and upper bounds (second col)" + + "self.domain, (d, 2) ndarray of domain lower bounds (first col) and upper bounds (second col)" ) if not hasattr(self, "range"): raise ParameterError( prefix - + "self.range, 2xd ndarray of range lower bounds (first col) and upper bounds (second col)" + + "self.range, (d, 2) ndarray of range lower bounds (first col) and upper bounds (second col)" ) if not hasattr(self, "parameters"): self.parameters = [] @@ -30,6 +36,52 @@ def _read_only_array(value): array.setflags(write=False) return array + @staticmethod + def _range_in_domain(transform_range, domain): + """Return whether a transform range is contained within a domain.""" + try: + transform_range = np.asarray(transform_range) + domain = np.asarray(domain) + except (TypeError, ValueError): + return False + + if ( + transform_range.ndim != 2 + or domain.ndim != 2 + or transform_range.shape[1] != 2 + or domain.shape[1] != 2 + or transform_range.shape[0] == 0 + or domain.shape[0] == 0 + ): + return False + + if not ( + np.issubdtype(transform_range.dtype, np.number) + and np.issubdtype(domain.dtype, np.number) + and np.isrealobj(transform_range) + and np.isrealobj(domain) + and transform_range.dtype != np.bool_ + and domain.dtype != np.bool_ + ): + return False + + if np.isnan(transform_range).any() or np.isnan(domain).any(): + return False + + if np.any(transform_range[:, 0] > transform_range[:, 1]) or np.any( + domain[:, 0] > domain[:, 1] + ): + return False + + try: + transform_range, domain = np.broadcast_arrays(transform_range, domain) + except ValueError: + return False + + lower_bounds_valid = np.all(domain[:, 0] <= transform_range[:, 0]) + upper_bounds_valid = np.all(transform_range[:, 1] <= domain[:, 1]) + return bool(lower_bounds_valid and upper_bounds_valid) + def _set_moments(self, mean, variance, standard_deviation, covariance): self._mean = self._read_only_array(mean) self._variance = self._read_only_array(variance) @@ -100,11 +152,11 @@ def _parse_sampler(self, sampler): sampler.d ) # take the dimension from the sub-sampler (composed transform) self.discrete_distrib = self.transform.discrete_distrib - if (self.domain != self.transform.range).any(): + if not self._range_in_domain(self.transform.range, self.domain): self.sub_compatibility_error = True if self.transform.sub_compatibility_error: raise ParameterError( - "The sub-transform domain must match the sub-sub-transform range." + "The nested sub-transform range must be contained within its transform domain." ) else: raise ParameterError( @@ -147,7 +199,7 @@ def _jacobian_transform_r(self, x, return_weights): jac = None if self.sub_compatibility_error: raise ParameterError( - "The transform domain must match the sub-transform range." + "The sub-transform range must be contained within the transform domain." ) if self.transform == self: # is \Psi_0 if return_weights: diff --git a/qmcpy/true_measure/bernoulli_cont.py b/qmcpy/true_measure/bernoulli_cont.py index cadf9f727..78dfd0751 100644 --- a/qmcpy/true_measure/bernoulli_cont.py +++ b/qmcpy/true_measure/bernoulli_cont.py @@ -75,6 +75,7 @@ def _transform(self, x): return tf def _weight(self, x): + in_support = np.all((0 <= x) & (x <= 1), axis=-1) w = np.zeros(x.shape, dtype=float) for j in range(self.d): C = ( @@ -83,7 +84,7 @@ def _weight(self, x): else 2 * np.arctanh(1 - 2 * self.l[j]) / (1 - 2 * self.l[j]) ) w[..., j] = C * self.l[j] ** x[..., j] * (1 - self.l[j]) ** (1 - x[..., j]) - return np.prod(w, -1) + return np.where(in_support, np.prod(w, -1), 0.0) def _spawn(self, sampler, dimension): if dimension == self.d: # don't do anything if the dimension doesn't change diff --git a/qmcpy/true_measure/brownian_motion.py b/qmcpy/true_measure/brownian_motion.py index 564223a00..88c85fe1b 100644 --- a/qmcpy/true_measure/brownian_motion.py +++ b/qmcpy/true_measure/brownian_motion.py @@ -1,4 +1,5 @@ from .gaussian import Gaussian +from .abstract_true_measure import _clip_unit_interval from ..discrete_distribution import DigitalNetB2 from ..util import ParameterError, ParameterWarning import warnings @@ -252,7 +253,7 @@ def _spawn(self, sampler, dimension): def _transform(self, x): if self.decomp_type == "BROWNIANBRIDGE": - z = norm.ppf(x) + z = norm.ppf(_clip_unit_interval(x)) w = self._bridge_transform(z) paths = self.drift_time_vec_plus_init + np.sqrt(self.diffusion) * w return paths[..., self._output_order] diff --git a/qmcpy/true_measure/copula.py b/qmcpy/true_measure/copula.py index 67ade5612..3851ebd94 100644 --- a/qmcpy/true_measure/copula.py +++ b/qmcpy/true_measure/copula.py @@ -2,7 +2,7 @@ import numpy as np -from .abstract_true_measure import AbstractTrueMeasure +from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval from ..util import DimensionError, MethodImplementationError, ParameterError @@ -97,11 +97,6 @@ def _unit_weight_with_warning(self, x): return np.ones(x.shape[:-1], dtype=float) -def _clip_unit_interval(u): - eps = np.finfo(float).eps - return np.clip(u, eps, 1.0 - eps) - - def _validate_marginals(marginals): try: parsed = list(marginals) diff --git a/qmcpy/true_measure/gaussian.py b/qmcpy/true_measure/gaussian.py index 37d58982b..7db292440 100644 --- a/qmcpy/true_measure/gaussian.py +++ b/qmcpy/true_measure/gaussian.py @@ -1,4 +1,4 @@ -from .abstract_true_measure import AbstractTrueMeasure +from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval from ..util import DimensionError, ParameterError from ..discrete_distribution import DigitalNetB2 import numpy as np @@ -158,6 +158,7 @@ def mvn_scipy(self, value): self._mvn_scipy_cache = value def _transform(self, x): + x = _clip_unit_interval(x) return self.mu + np.einsum("...ij,kj->...ik", norm.ppf(x), self.a) def _weight(self, t): diff --git a/qmcpy/true_measure/johnsons_su.py b/qmcpy/true_measure/johnsons_su.py index 1d0cc54e2..d04dfbe2e 100644 --- a/qmcpy/true_measure/johnsons_su.py +++ b/qmcpy/true_measure/johnsons_su.py @@ -1,4 +1,4 @@ -from .abstract_true_measure import AbstractTrueMeasure +from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval from ..util import DimensionError, ParameterError from ..discrete_distribution import DigitalNetB2 import numpy as np @@ -92,6 +92,7 @@ def __init__(self, sampler, gamma=1, xi=1, delta=2, lam=2): ) def _transform(self, x): + x = _clip_unit_interval(x) return self._lam * np.sinh((norm.ppf(x) - self._gamma) / self._delta) + self._xi def _weight(self, x): diff --git a/qmcpy/true_measure/kumaraswamy.py b/qmcpy/true_measure/kumaraswamy.py index 47d4bdd9d..3419acec0 100644 --- a/qmcpy/true_measure/kumaraswamy.py +++ b/qmcpy/true_measure/kumaraswamy.py @@ -158,13 +158,17 @@ def _transform(self, x): return (1 - (1 - x) ** (1 / self.beta)) ** (1 / self.alpha) def _weight(self, x): - return np.prod( - self.alpha - * self.beta - * x ** (self.alpha - 1) - * (1 - x**self.alpha) ** (self.beta - 1), - -1, - ) + in_support = np.all((0 <= x) & (x <= 1), axis=-1) + x_in_support = np.clip(x, 0, 1) + with np.errstate(divide="ignore", invalid="ignore"): + weight = np.prod( + self.alpha + * self.beta + * x_in_support ** (self.alpha - 1) + * (1 - x_in_support**self.alpha) ** (self.beta - 1), + -1, + ) + return np.where(in_support, weight, 0.0) def _spawn(self, sampler, dimension): if dimension == self.d: # don't do anything if the dimension doesn't change diff --git a/qmcpy/true_measure/product_measure.py b/qmcpy/true_measure/product_measure.py index 3f860e84b..258bff693 100644 --- a/qmcpy/true_measure/product_measure.py +++ b/qmcpy/true_measure/product_measure.py @@ -105,25 +105,22 @@ def __init__(self, sampler, marginals): """ Initialize a product measure from one sampler and several marginals. - Parameters - ---------- - sampler : AbstractDiscreteDistribution - The sampler for the whole product measure. Its dimension must - equal the sum of the marginal dimensions. + Args: + sampler (AbstractDiscreteDistribution): Sampler for the whole product measure. Its dimension must equal the sum of the marginal dimensions. + marginals (Union[list, tuple]): Nonempty sequence of independent `AbstractTrueMeasure` instances to place side by side. A marginal may itself be multidimensional. - marginals : list or tuple of AbstractTrueMeasure - Independent true measures to place side by side. A marginal may - itself be multidimensional. + Raises: + ParameterError: If `sampler` is not an `AbstractDiscreteDistribution`, or if `marginals` is empty or contains a non-`AbstractTrueMeasure` value. + DimensionError: If a marginal is not dimension-preserving or the sampler dimension differs from the sum of the marginal dimensions. - Why one sampler? - ---------------- - The product measure should be driven by one total-dimensional QMC - point set. We do not generate separate QMC samples from each marginal. - Instead, one sample u in [0,1]^d is split into blocks: + Note: + The product measure is driven by one total-dimensional QMC point + set. It does not generate separate QMC samples from each marginal. + Instead, one sample u in [0,1]^d is split into blocks: - u = (u_marginal_1, u_marginal_2, ..., u_marginal_k). + u = (u_marginal_1, u_marginal_2, ..., u_marginal_k). - This preserves the intended total-dimensional QMC construction. + This preserves the intended total-dimensional QMC construction. """ if not isinstance(marginals, (list, tuple)) or len(marginals) == 0: raise ParameterError("ProductMeasure requires a nonempty list of marginals.") diff --git a/qmcpy/true_measure/scipy_wrapper.py b/qmcpy/true_measure/scipy_wrapper.py index 9f528b175..6fbb24c6c 100644 --- a/qmcpy/true_measure/scipy_wrapper.py +++ b/qmcpy/true_measure/scipy_wrapper.py @@ -1,4 +1,4 @@ -from .abstract_true_measure import AbstractTrueMeasure +from .abstract_true_measure import AbstractTrueMeasure, _clip_unit_interval from ..util import DimensionError, ParameterError from ..discrete_distribution.abstract_discrete_distribution import ( AbstractDiscreteDistribution, @@ -107,8 +107,7 @@ def transform(self, u): ) # Clip so we never hit exactly 0 or 1 inside norm.ppf. - eps = np.finfo(float).eps - u_clip = np.clip(u, eps, 1.0 - eps) + u_clip = _clip_unit_interval(u) # Map to i.i.d. standard normals. z = scipy.stats.norm.ppf(u_clip) @@ -455,6 +454,7 @@ def _transform(self, x): if self._is_joint: return self._joint.transform(x) + x = _clip_unit_interval(x) t = np.empty_like(x, dtype=float) for j in range(self.d): t[..., j] = self.sds[j].ppf(x[..., j]) diff --git a/qmcpy/true_measure/student_t.py b/qmcpy/true_measure/student_t.py index 510b11388..e0e44ae51 100644 --- a/qmcpy/true_measure/student_t.py +++ b/qmcpy/true_measure/student_t.py @@ -2,6 +2,7 @@ import scipy.stats as stats from ..util import ParameterError, DimensionError +from .abstract_true_measure import _clip_unit_interval from .scipy_wrapper import SciPyWrapper @@ -39,8 +40,7 @@ def __init__(self, loc, shape, df): @staticmethod def _clip_u(u): - eps = np.finfo(float).eps - return np.clip(u, eps, 1.0 - eps) + return _clip_unit_interval(u) def transform(self, u): u = np.asarray(u, dtype=float) diff --git a/qmcpy/true_measure/uniform.py b/qmcpy/true_measure/uniform.py index 0dc8c496e..97f25ab70 100644 --- a/qmcpy/true_measure/uniform.py +++ b/qmcpy/true_measure/uniform.py @@ -100,7 +100,8 @@ def _transform(self, x): return x * self.delta + self.a def _weight(self, x): - return np.tile(self.inv_delta_prod, x.shape[:-1]) + in_support = np.all((self.a <= x) & (x <= self.b), axis=-1) + return np.where(in_support, self.inv_delta_prod, 0.0) def _spawn(self, sampler, dimension): if dimension == self.d: # don't do anything if the dimension doesn't change diff --git a/qmcpy/true_measure/zero_inflated_exp_uniform.py b/qmcpy/true_measure/zero_inflated_exp_uniform.py index 11526556f..2fbb6bfd6 100644 --- a/qmcpy/true_measure/zero_inflated_exp_uniform.py +++ b/qmcpy/true_measure/zero_inflated_exp_uniform.py @@ -187,6 +187,17 @@ class ZeroInflatedExpUniform(SciPyWrapper): """ def __init__(self, sampler, p_zero=0.4, lam=1.5, y_split=None): + r""" + Args: + sampler (Union[AbstractDiscreteDistribution, AbstractTrueMeasure]): One-dimensional sampler for the current construction. The deprecated `y_split` construction also accepts a two-dimensional sampler. + p_zero (float): Probability mass at zero, strictly between `0` and `1`. Defaults to `0.4`. + lam (float): Rate of the exponential component. Must be positive. Defaults to `1.5`. + y_split (Union[None, float]): Deprecated split point for the legacy two-dimensional construction. With a two-dimensional sampler, it must lie strictly between `0` and `1`. With a one-dimensional sampler, it is accepted for backward compatibility, emits a `DeprecationWarning`, and is otherwise ignored. Defaults to `None`. + + Raises: + DimensionError: If the sampler dimension is incompatible with the selected construction. + ParameterError: If `p_zero`, `lam`, or a two-dimensional `y_split` is outside its valid range. + """ if y_split is not None: warnings.warn( "`y_split` is deprecated. The 2D zero-inflated " diff --git a/scripts/colab_notebooks_manifest.json b/scripts/colab_notebooks_manifest.json index ed62ba8be..82b4a4be4 100644 --- a/scripts/colab_notebooks_manifest.json +++ b/scripts/colab_notebooks_manifest.json @@ -44,6 +44,7 @@ "demos/talk_paper_demos/Sorokin_random_LD_seq_QMC_fast_kernel_methods_2026/Sorokin_random_LD_seq_QMC_fast_kernel_methods_2026.ipynb", "demos/talk_paper_demos/pydata_chi_2023.ipynb", "demos/talk_paper_demos/why_add_q_to_mc_blog/why_add_q_to_mc_blog.ipynb", + "demos/true_measure_domain_inclusion.ipynb", "demos/vectorized_qmc.ipynb", "demos/vectorized_qmc_bayes.ipynb" ], diff --git a/test/booktests/tb_true_measure_domain_inclusion.py b/test/booktests/tb_true_measure_domain_inclusion.py new file mode 100644 index 000000000..b34abcb60 --- /dev/null +++ b/test/booktests/tb_true_measure_domain_inclusion.py @@ -0,0 +1,18 @@ +import unittest +from testbook import testbook +from __init__ import TB_TIMEOUT, BaseNotebookTest + + +class NotebookTests(BaseNotebookTest): + + @testbook( + "../../demos/true_measure_domain_inclusion.ipynb", + execute=True, + timeout=TB_TIMEOUT, + ) + def test_true_measure_domain_inclusion_notebook(self, tb): + pass + + +if __name__ == "__main__": + unittest.main() diff --git a/test/test_true_measures.py b/test/test_true_measures.py index ae58d9604..9dd78a5b3 100644 --- a/test/test_true_measures.py +++ b/test/test_true_measures.py @@ -1,4 +1,5 @@ from qmcpy import ( + AbstractTrueMeasure, BernoulliCont, BrownianMotion, DigitalNetB2, @@ -10,17 +11,19 @@ Lattice, Lebesgue, MaternGP, + SciPyWrapper, + StudentT, Uniform, ZeroInflatedExpUniform, ) -from qmcpy.util import DimensionError, ParameterError +from qmcpy.util import DimensionError, ParameterError, ParameterWarning import numpy as np +import re import scipy.stats from scipy.sparse import issparse import unittest import warnings from qmcpy.true_measure.uniform_triangle import UniformTriangle, _UniformTriangleAdapter -from qmcpy import SciPyWrapper def dense_covariance(covariance): @@ -42,6 +45,307 @@ def assert_sample_mean_and_covariance(measure): class TestTrueMeasure(unittest.TestCase): """General tests for TrueMeasures""" + def test_range_in_domain(self): + cases = [ + ("exact equality", [[0, 1]], [[0, 1]], True), + ("strict finite inclusion", [[0.25, 0.75]], [[0, 1]], True), + ( + "compatible unbounded interval", + [[-np.inf, np.inf]], + [[-np.inf, np.inf]], + True, + ), + ("infinite domain", [[-5, 5]], [[-np.inf, np.inf]], True), + ("positive half-line", [[1, 3]], [[0, np.inf]], True), + ("lower-bound failure", [[-0.1, 0.75]], [[0, 1]], False), + ("upper-bound failure", [[0.25, 1.1]], [[0, 1]], False), + ( + "multidimensional box", + [[0.1, 0.8], [0.2, 0.9]], + [[0, 1], [0, 1]], + True, + ), + ( + "failure in one coordinate", + [[0.1, 0.8], [0.2, 1.1]], + [[0, 1], [0, 1]], + False, + ), + ( + "broadcast domain", + [[0.1, 0.8], [0.2, 0.9]], + [[0, 1]], + True, + ), + ( + "broadcast transform range", + [[0.25, 0.75]], + [[0, 1], [-1, 2]], + True, + ), + ( + "non-broadcastable rows", + [[0.1, 0.8], [0.2, 0.9]], + [[0, 1], [0, 1], [0, 1]], + False, + ), + ] + + for name, transform_range, domain, expected in cases: + with self.subTest(name=name): + self.assertIs( + AbstractTrueMeasure._range_in_domain(transform_range, domain), + expected, + ) + + def test_range_in_domain_rejects_invalid_bounds(self): + invalid_cases = [ + ("one-dimensional range", [0, 1], [[0, 1]]), + ("three-column range", [[0, 0.5, 1]], [[0, 1]]), + ("ragged range", [[0, 1], [0, 0.5, 1]], [[0, 1]]), + ("one-dimensional domain", [[0, 1]], [0, 1]), + ("three-column domain", [[0, 1]], [[0, 0.5, 1]]), + ( + "reversed transform range", + np.array([[0.8, 0.2]]), + np.array([[0.0, 1.0]]), + ), + ("reversed domain", [[0, 1]], [[0.8, 0.2]]), + ( + "reversed multidimensional transform range", + [[0.1, 0.8], [0.9, 0.2]], + [[0, 1], [0, 1]], + ), + ( + "reversed multidimensional domain", + [[0.1, 0.8], [0.2, 0.9]], + [[0, 1], [0.9, 0.2]], + ), + ("empty transform range", np.empty((0, 2)), [[0, 1]]), + ("empty domain", [[0, 1]], np.empty((0, 2))), + ( + "string bounds", + np.array([["a", "z"]]), + np.array([["a", "z"]]), + ), + ( + "object bounds", + np.array([[0, 1]], dtype=object), + np.array([[0, 1]], dtype=object), + ), + ("complex bounds", [[0 + 0j, 1 + 0j]], [[0 + 0j, 1 + 0j]]), + ("boolean bounds", [[False, True]], [[False, True]]), + ("NaN transform range", [[np.nan, 1]], [[0, 1]]), + ("NaN domain", [[0, 1]], [[np.nan, 1]]), + ] + + for name, transform_range, domain in invalid_cases: + with self.subTest(name=name): + self.assertIs( + AbstractTrueMeasure._range_in_domain(transform_range, domain), + False, + ) + + def test_strict_range_in_domain_chain(self): + inner = Uniform( + DigitalNetB2(1, seed=7), lower_bound=0.25, upper_bound=0.75 + ) + outer = Kumaraswamy(inner) + + self.assertFalse(outer.sub_compatibility_error) + samples = outer.gen_samples(8) + self.assertEqual(samples.shape, (8, 1)) + self.assertTrue(np.isfinite(samples).all()) + + def test_multidimensional_range_in_domain_chain_broadcasts(self): + inner = Uniform( + DigitalNetB2(2, seed=7), + lower_bound=[0.1, 0.2], + upper_bound=[0.8, 0.9], + ) + outer = Kumaraswamy(inner) + + self.assertFalse(outer.sub_compatibility_error) + samples = outer.gen_samples(8) + self.assertEqual(samples.shape, (8, 2)) + self.assertTrue(np.isfinite(samples).all()) + + def test_recursive_transform_applies_each_layer(self): + # This records current recursive execution only; it does not assert + # correctness of nominal range or moment metadata. + points = np.array([[0.1], [0.25], [0.5], [0.75], [0.9]]) + + uniform_inner = Uniform( + DigitalNetB2(1, seed=7), lower_bound=0.25, upper_bound=0.75 + ) + uniform_outer = Uniform( + uniform_inner, lower_bound=0.25, upper_bound=0.75 + ) + + kumaraswamy_inner = Uniform( + DigitalNetB2(1, seed=7), lower_bound=0.25, upper_bound=0.75 + ) + kumaraswamy_outer = Kumaraswamy(kumaraswamy_inner) + + gaussian_inner = Uniform( + DigitalNetB2(1, seed=7), lower_bound=0.25, upper_bound=0.75 + ) + gaussian_outer = Gaussian(gaussian_inner) + + bernoulli_inner = BernoulliCont(DigitalNetB2(1, seed=7), lam=0.9) + bernoulli_outer = Kumaraswamy(bernoulli_inner, a=2.0, b=2.0) + + cases = [ + ("1(a) Uniform -> Uniform", uniform_inner, uniform_outer), + ( + "1(b) Uniform -> Kumaraswamy", + kumaraswamy_inner, + kumaraswamy_outer, + ), + ("1(c) Uniform -> Gaussian", gaussian_inner, gaussian_outer), + ( + "1(d) BernoulliCont -> Kumaraswamy", + bernoulli_inner, + bernoulli_outer, + ), + ] + + for name, inner, outer in cases: + with self.subTest(name=name): + transformed = outer._jacobian_transform_r( + points, return_weights=False + ) + expected = outer._transform(inner._transform(points)) + np.testing.assert_allclose(transformed, expected) + + if name.startswith("1(a)"): + np.testing.assert_allclose( + transformed, 0.375 + 0.25 * points + ) + if name.startswith("1(d)"): + np.testing.assert_array_equal(inner.range, outer.domain) + + def test_unrandomized_inverse_cdf_paths_are_finite(self): + cases = [ + ( + "JohnsonsSU", + JohnsonsSU(DigitalNetB2(1, randomize="FALSE")), + ), + ( + "BrownianMotion BrownianBridge", + BrownianMotion( + DigitalNetB2(2, randomize="FALSE"), + decomp_type="BROWNIANBRIDGE", + ), + ), + ( + "SciPyWrapper marginal", + SciPyWrapper( + DigitalNetB2(1, randomize="FALSE"), + scipy.stats.norm(), + ), + ), + ( + "SciPyWrapper MVN", + SciPyWrapper( + DigitalNetB2(2, randomize="FALSE"), + scipy.stats.multivariate_normal( + mean=[0.0, 0.0], + cov=[[1.0, 0.5], [0.5, 1.0]], + ), + ), + ), + ( + "StudentT", + StudentT( + DigitalNetB2(2, randomize="FALSE"), + loc=[0.0, 0.0], + shape=[[1.0, 0.25], [0.25, 1.0]], + df=5, + ), + ), + ] + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", ParameterWarning) + for name, measure in cases: + with self.subTest(name=name): + samples = measure.gen_samples(2) + self.assertTrue(np.isfinite(samples).all()) + + def test_inverse_cdf_clipping_preserves_interior_values(self): + endpoints_1d = np.array([[0.0], [1.0]]) + points_1d = np.array([[0.1], [0.25], [0.5], [0.75], [0.9]]) + + johnsons_su = JohnsonsSU(DigitalNetB2(1, seed=7)) + self.assertTrue(np.isfinite(johnsons_su._transform(endpoints_1d)).all()) + johnsons_expected = johnsons_su._lam * np.sinh( + (scipy.stats.norm.ppf(points_1d) - johnsons_su._gamma) + / johnsons_su._delta + ) + johnsons_su._xi + np.testing.assert_allclose( + johnsons_su._transform(points_1d), + johnsons_expected, + rtol=0, + atol=0, + ) + + scipy_wrapper = SciPyWrapper(DigitalNetB2(1, seed=7), scipy.stats.norm()) + self.assertTrue(np.isfinite(scipy_wrapper._transform(endpoints_1d)).all()) + np.testing.assert_allclose( + scipy_wrapper._transform(points_1d), + scipy.stats.norm.ppf(points_1d), + rtol=0, + atol=0, + ) + + points_2d = np.array( + [[0.1, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 0.9]] + ) + brownian_bridge = BrownianMotion( + DigitalNetB2(2, seed=7), + decomp_type="BROWNIANBRIDGE", + ) + endpoints_2d = np.array([[0.0, 0.0], [1.0, 1.0]]) + self.assertTrue( + np.isfinite(brownian_bridge._transform(endpoints_2d)).all() + ) + bridge_normals = scipy.stats.norm.ppf(points_2d) + bridge_expected = ( + brownian_bridge.drift_time_vec_plus_init + + np.sqrt(brownian_bridge.diffusion) + * brownian_bridge._bridge_transform(bridge_normals) + )[..., brownian_bridge._output_order] + np.testing.assert_allclose( + brownian_bridge._transform(points_2d), + bridge_expected, + rtol=0, + atol=0, + ) + + def test_out_of_domain_chain_preserves_deferred_errors(self): + inner = Uniform( + DigitalNetB2(1, seed=7), lower_bound=-0.1, upper_bound=0.75 + ) + incompatible = Kumaraswamy(inner) + + self.assertTrue(incompatible.sub_compatibility_error) + with self.assertRaisesRegex( + ParameterError, + re.escape( + "The sub-transform range must be contained within the transform domain." + ), + ): + incompatible.gen_samples(8) + + with self.assertRaisesRegex( + ParameterError, + re.escape( + "The nested sub-transform range must be contained within its transform domain." + ), + ): + Kumaraswamy(incompatible) + def test_abstract_methods(self): d = 2 tms = [ @@ -381,6 +685,34 @@ def test_spawn_recomputes_moment_attributes(self): ) np.testing.assert_allclose(dense_covariance(spawn.covariance), 3.0 * np.eye(4)) + def test_weight_is_zero_outside_support(self): + uniform = Uniform(DigitalNetB2(1, seed=7), 0.25, 0.75) + points = np.array([[0.1], [0.25], [0.3], [0.5], [0.75], [0.9]]) + + np.testing.assert_allclose( + uniform._weight(points), + [0.0, 2.0, 2.0, 2.0, 2.0, 0.0], + ) + + def test_weight_support_mask_handles_dimensions_and_batches(self): + uniform = Uniform( + DigitalNetB2(2, seed=7), + lower_bound=[0.0, -1.0], + upper_bound=[1.0, 1.0], + ) + points = np.array( + [ + [[0.0, -1.0], [0.5, 0.0], [1.0, 1.0]], + [[-0.1, 0.0], [0.5, 1.1], [0.5, 0.0]], + ] + ) + + weights = uniform._weight(points) + + self.assertEqual(weights.shape, (2, 3)) + np.testing.assert_allclose(weights, [[0.5, 0.5, 0.5], [0.0, 0.0, 0.5]]) + self.assertTrue(np.all(weights >= 0)) + class TestKumaraswamy(unittest.TestCase): def test_sample_mean_and_covariance(self): @@ -437,6 +769,36 @@ def test_uniform_special_case(self): dense_covariance(kumaraswamy.covariance), np.diag(expected_variance) ) + def test_weight_is_zero_outside_support(self): + kumaraswamy = Kumaraswamy(DigitalNetB2(1, seed=7), a=2, b=2) + points = np.array([[-0.5], [0.0], [0.5], [1.0], [1.5]]) + + weights = kumaraswamy._weight(points) + + np.testing.assert_allclose(weights, [0.0, 0.0, 1.5, 0.0, 0.0]) + self.assertTrue(np.all(weights >= 0)) + + singular_boundaries = Kumaraswamy( + DigitalNetB2(1, seed=7), a=0.5, b=0.5 + )._weight(np.array([[0.0], [1.0]])) + self.assertTrue(np.isposinf(singular_boundaries).all()) + + def test_weight_support_mask_handles_dimensions_and_batches(self): + kumaraswamy = Kumaraswamy(DigitalNetB2(2, seed=7), a=2, b=2) + points = np.array( + [ + [[0.5, 0.5], [0.25, 0.75]], + [[-0.1, 0.5], [0.5, 1.1]], + ] + ) + + weights = kumaraswamy._weight(points) + + self.assertEqual(weights.shape, (2, 2)) + np.testing.assert_allclose(weights[0], [2.25, 1.23046875]) + np.testing.assert_allclose(weights[1], [0.0, 0.0]) + self.assertTrue(np.all(weights >= 0)) + def test_spawn_recomputes_moment_attributes(self): kumaraswamy = Kumaraswamy( DigitalNetB2(2, seed=7), a=1, b=3 @@ -636,6 +998,45 @@ def test_deprecated_2d_construction_has_no_moment_parameters(self): self.assertNotIn(parameter, tm.parameters) +class TestBernoulliCont(unittest.TestCase): + def test_weight_is_zero_outside_support(self): + lam = 0.9 + bernoulli = BernoulliCont(DigitalNetB2(1, seed=7), lam=lam) + points = np.array([[-0.1], [0.0], [0.5], [1.0], [1.1]]) + + weights = bernoulli._weight(points) + normalizer = 2 * np.arctanh(1 - 2 * lam) / (1 - 2 * lam) + expected_in_support = ( + normalizer + * lam ** points[1:4, 0] + * (1 - lam) ** (1 - points[1:4, 0]) + ) + + np.testing.assert_allclose(weights[[0, 4]], 0.0) + np.testing.assert_allclose(weights[1:4], expected_in_support) + self.assertTrue(np.all(weights >= 0)) + + def test_weight_support_mask_handles_dimensions_and_batches(self): + bernoulli = BernoulliCont( + DigitalNetB2(2, seed=7), + lam=[0.9, 0.8], + ) + points = np.array( + [ + [[0.0, 0.0], [0.5, 0.5], [1.0, 1.0]], + [[-0.1, 0.5], [0.5, 1.1], [0.25, 0.75]], + ] + ) + + weights = bernoulli._weight(points) + + self.assertEqual(weights.shape, (2, 3)) + self.assertTrue(np.all(weights[0] > 0)) + np.testing.assert_allclose(weights[1, :2], 0.0) + self.assertGreater(weights[1, 2], 0.0) + self.assertTrue(np.all(weights >= 0)) + + class TestUniformTriangle(unittest.TestCase): """Tests for UniformTriangle and _UniformTriangleAdapter.""" @@ -702,6 +1103,36 @@ def test_sample_mean_and_covariance(self): assert_sample_mean_and_covariance(gaussian) + def test_transform_clips_unit_interval_endpoints(self): + gaussian = Gaussian(DigitalNetB2(1, seed=7)) + endpoints = gaussian._transform(np.array([[0.0], [1.0]])) + interior = np.array([[0.25], [0.5], [0.75]]) + + self.assertTrue(np.isfinite(endpoints).all()) + np.testing.assert_allclose( + gaussian._transform(interior), + scipy.stats.norm.ppf(interior), + rtol=0, + atol=0, + ) + + def test_unrandomized_direct_and_composed_samples_are_finite(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", ParameterWarning) + direct = Gaussian( + DigitalNetB2(1, randomize="FALSE") + ).gen_samples(2) + composed = Gaussian( + Uniform( + DigitalNetB2(1, randomize="FALSE"), + 0.0, + 0.75, + ) + ).gen_samples(2) + + self.assertTrue(np.isfinite(direct).all()) + self.assertTrue(np.isfinite(composed).all()) + def test_gaussian_basic_output_reproducibility(self): """Test that basic Gaussian sample generation produces expected values with fixed seed.""" gaussian = Gaussian(Lattice(4, seed=self.seed), mean=0, covariance=1) @@ -1154,7 +1585,6 @@ def test_brownian_bridge_output_order(self): def test_brownian_bridge_warning_for_non_power_of_2(self): """BrownianBridge issues ParameterWarning for suboptimal d but still produces valid output.""" - from qmcpy.util import ParameterWarning with self.assertWarns(ParameterWarning): bm = BrownianMotion(DigitalNetB2(6, seed=self.seed), decomp_type='BrownianBridge') samples = bm.gen_samples(4)