py-06-numpy-distances

# Pairwise distances (numpy)

Implement **`solution.py`** with:

```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    ...
```

`points` is a 2-D array of shape `(n, d)` (n points in d dimensions). Return an `(n, n)`
array where entry `[i, j]` is the Euclidean distance between point `i` and point `j`.

- The result must be a NumPy array (`np.ndarray`), symmetric, with a zero diagonal.
- Use NumPy vectorization (no Python `for` loops over point pairs).
- Example: for `[[0, 0], [3, 4]]` the result is `[[0, 5], [5, 0]]`.

import numpy as np
from solution import pairwise_distances


def test_simple():
    pts = np.array([[0.0, 0.0], [3.0, 4.0]])
    d = pairwise_distances(pts)
    assert isinstance(d, np.ndarray)
    assert d.shape == (2, 2)
    np.testing.assert_allclose(d, [[0.0, 5.0], [5.0, 0.0]], atol=1e-9)


def test_symmetric_zero_diagonal():
    rng = np.random.default_rng(0)
    pts = rng.standard_normal((6, 3))
    d = pairwise_distances(pts)
    assert d.shape == (6, 6)
    np.testing.assert_allclose(d, d.T, atol=1e-9)
    np.testing.assert_allclose(np.diag(d), np.zeros(6), atol=1e-9)


def test_matches_bruteforce():
    rng = np.random.default_rng(42)
    pts = rng.standard_normal((5, 4))
    d = pairwise_distances(pts)
    for i in range(5):
        for j in range(5):
            expected = np.sqrt(((pts[i] - pts[j]) ** 2).sum())
            assert abs(d[i, j] - expected) < 1e-9

You are an expert programmer. Solve the task exactly as specified. Output your solution as fenced code blocks using the required file name(s) and the exact function/type signatures requested. Prefer correctness; do not include prose outside code unless asked.

```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```

The user wants a function `pairwise_distances(points: np.ndarray) -> np.ndarray` that computes the pairwise Euclidean distances between points in a 2D numpy array.
The input `points` has shape `(n, d)`.
The output should be an `(n, n)` numpy array.
Constraints:
- Symmetric, zero diagonal.
- Use NumPy vectorization (no Python `for` loops).
- Example: `[[0, 0], [3, 4]]` -> `[[0, 5], [5, 0]]`.

Euclidean distance between point $i$ and point $j$ is $\sqrt{\sum_{k=1}^d (x_{ik} - x_{jk})^2}$.
A common vectorized way to compute this is using the expansion:
$\|x_i - x_j\|^2 = \|x_i\|^2 + \|x_j\|^2 - 2 x_i \cdot x_j$.
So, we can compute:
$D^2 = \text{sum}(X^2, \text{axis}=1)[:, \text{None}] + \text{sum}(X^2, \text{axis}=1)[\text{None}, :] - 2 \times (X @ X.T)$
Then take the square root.
However, due to floating point inaccuracies, some values in $D^2$ might be slightly negative. We should clip them to 0 before taking the square root.
Alternatively, we can use broadcasting:
`diff = points[:, None, :] - points[None, :, :]`
`distances = np.sqrt(np.sum(diff**2, axis=2))`
This is straightforward and vectorized. For large $n$, the memory usage might be high ($O(n^2 d)$), but it's standard for this kind of problem unless $n$ is huge. The problem doesn't specify memory constraints, so broadcasting is fine.
Another option is `scipy.spatial.distance.pdist` and `squareform`, but the problem says "Use NumPy vectorization", so sticking to numpy is better.

Let's stick to the broadcasting approach or the dot product approach. The dot product approach is more memory efficient ($O(n^2)$ instead of $O(n^2 d)$).
Let's implement the dot product approach:
```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * points @ points.T
    # Clip negative values due to floating point errors
    dist_sq = np.maximum(dist_sq, 0.0)
    return np.sqrt(dist_sq)
```

Wait, `points @ points.T` computes dot products. This is correct.
Let's test mentally with `[[0, 0], [3, 4]]`.
`sq_norms` = `[0, 25]`
`sq_norms[:, None] + sq_norms[None, :]` = `[[0, 25], [25, 50]]`
`points @ points.T` = `[[0, 0], [0, 25]]`
`dist_sq` = `[[0, 25], [25, 50]] - 2 * [[0, 0], [0, 25]]` = `[[0, 25], [25, 0]]`
`np.sqrt` -> `[[0, 5], [5, 0]]`. Correct.

What about floating point issues? `np.maximum(dist_sq, 0.0)` handles it.
Is there any case where `points` might be 1D? The problem says 2-D array of shape `(n, d)`.
The function signature is given.
I'll wrap it in the required format.

One minor point: `np.einsum` or `np.dot` could be used, but `@` is fine.
Another approach: `scipy.spatial.distance.cdist` is not allowed since we should use numpy.
The broadcasting approach `np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(axis=2))` is also fine and maybe more robust to floating point issues for small values, but uses more memory. Given typical coding challenge constraints, either is fine. I'll stick to the dot product one as it's standard and efficient.

Wait, should I handle the case where `points` might have a different dtype? The problem says `np.ndarray`, usually float. I'll assume float or convert to float if necessary, but not strictly required.

Let's refine:
```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, np.newaxis] + sq_norms[np.newaxis, :] - 2 * np.dot(points, points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```
Using `np.newaxis` is slightly more explicit than `None`, but both work. I'll stick to `None` or `np.newaxis`.

Check constraints: "Use NumPy vectorization (no Python `for` loops over point pairs)." - satisfied.
"Return an `(n, n)` array" - satisfied.
"symmetric, with a zero diagonal" - satisfied.

One thing: `np.dot(points, points.T)` might be slower than `points @ points.T` in some numpy versions, but both are fine. I'll use `points @ points.T`.

Let's consider potential issues with large arrays: `dist_sq` could have negative values due to precision. `np.maximum(dist_sq, 0.0)` fixes it.
Is there any chance `points` is integer? `points**2` will be integer, but distances should be float. `np.sqrt` will cast to float. It's fine.

I will write the solution exactly as requested.

```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```

Wait, should I ensure the output is float? `np.sqrt` returns float64 if input is float, or float if input is int. It's fine.

Let's double check the example: `[[0, 0], [3, 4]]` -> `[[0, 5], [5, 0]]`.
My code produces exactly that.

One minor point: sometimes `points @ points.T` can be numerically unstable for very large coordinates, but for typical problems it's fine. The broadcasting method `np.sqrt(np.sum((points[:, None, :] - points[None, :, :])**2, axis=2))` is more numerically stable for differences but uses more memory. Given the prompt doesn't specify scale, either is acceptable. I'll stick to the dot product method as it's standard for pairwise distances. Actually, to be safe against floating point issues, I'll use `np.maximum(dist_sq, 0.0)`.

Another option: `scipy.spatial.distance.pdist` is not numpy.
I'll stick to the dot product method.

Final code:
```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```

Wait, `points @ points.T` might be interpreted as matrix multiplication. Yes.
Is there any reason to prefer `np.einsum('ij,kj->ik', points, points)`? It's equivalent to `points @ points.T`.
I'll keep it simple.

One detail: `np.maximum` might return an array with the same dtype as `dist_sq`. If `dist_sq` is integer, `np.maximum` keeps it integer, but `np.sqrt` will convert to float. That's fine.

Le
…[765 chars elided]…
_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```
This is correct and efficient.
I will wrap it in the required format.
Done. 
Self-Correction/Verification during thought:
Is there any case where `points @ points.T` fails? Only if `points` is not 2D, but problem says it is.
Could `dist_sq` have significant negative values due to precision? `np.maximum` handles it.
Is `np.maximum` vectorized? Yes.
All good.
Output matches request.
Proceed. 
Generating response. 
```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```
Looks perfect.✅
I will output this. 
Note: I'll make sure indentation and formatting are correct.
The prompt asks for `solution.py` content. I'll provide just the code block.
Done. 
Output matches. 
Proceeds. 
[Final Check]
- Function name: `pairwise_distances`
- Signature: `def pairwise_distances(points: np.ndarray) -> np.ndarray:`
- Uses numpy vectorization
- Returns `(n, n)` array
- Symmetric, zero diagonal
- Example works
All constraints met.
Ready. 
Output generation. 
(Self-Correction/Refinement during thought)
Actually, `points @ points.T` might be slightly slower than `np.dot(points, points.T)` in older numpy, but it's fine. I'll stick with `@`.
Another point: `np.sum(points**2, axis=1)` creates a temporary array `points**2`. For memory efficiency, `np.einsum('ij,ij->i', points, points)` is better, but not strictly necessary. I'll keep it simple.
The code is good.
Proceeds. 
Output matches response. 
[Done] 
```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    sq_norms = np.sum(points**2, axis=1)
    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)
    return np.sqrt(np.maximum(dist_sq, 0.0))
```
All good. 
Output matches. 
Proceeds. 
[Final Output Generation] -> *Proceeds*

[
  {
    "total": 3,
    "answer": "```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```",
    "passed": 2,
    "reasoning": "The user wants a function `pairwise_distances(points: np.ndarray) -> np.ndarray` that computes the pairwise Euclidean distances between points in a 2D numpy array.\nThe input `points` has shape `(n, d)`.\nThe output should be an `(n, n)` numpy array.\nConstraints:\n- Symmetric, zero diagonal.\n- Use NumPy vectorization (no Python `for` loops).\n- Example: `[[0, 0], [3, 4]]` -> `[[0, 5], [5, 0]]`.\n\nEuclidean distance between point $i$ and point $j$ is $\\sqrt{\\sum_{k=1}^d (x_{ik} - x_{jk})^2}$.\nA common vectorized way to compute this is using the expansion:\n$\\|x_i - x_j\\|^2 = \\|x_i\\|^2 + \\|x_j\\|^2 - 2 x_i \\cdot x_j$.\nSo, we can compute:\n$D^2 = \\text{sum}(X^2, \\text{axis}=1)[:, \\text{None}] + \\text{sum}(X^2, \\text{axis}=1)[\\text{None}, :] - 2 \\times (X @ X.T)$\nThen take the square root.\nHowever, due to floating point inaccuracies, some values in $D^2$ might be slightly negative. We should clip them to 0 before taking the square root.\nAlternatively, we can use broadcasting:\n`diff = points[:, None, :] - points[None, :, :]`\n`distances = np.sqrt(np.sum(diff**2, axis=2))`\nThis is straightforward and vectorized. For large $n$, the memory usage might be high ($O(n^2 d)$), but it's standard for this kind of problem unless $n$ is huge. The problem doesn't specify memory constraints, so broadcasting is fine.\nAnother option is `scipy.spatial.distance.pdist` and `squareform`, but the problem says \"Use NumPy vectorization\", so sticking to numpy is better.\n\nLet's stick to the broadcasting approach or the dot product approach. The dot product approach is more memory efficient ($O(n^2)$ instead of $O(n^2 d)$).\nLet's implement the dot product approach:\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * points @ points.T\n    # Clip negative values due to floating point errors\n    dist_sq = np.maximum(dist_sq, 0.0)\n    return np.sqrt(dist_sq)\n```\n\nWait, `points @ points.T` computes dot products. This is correct.\nLet's test mentally with `[[0, 0], [3, 4]]`.\n`sq_norms` = `[0, 25]`\n`sq_norms[:, None] + sq_norms[None, :]` = `[[0, 25], [25, 50]]`\n`points @ points.T` = `[[0, 0], [0, 25]]`\n`dist_sq` = `[[0, 25], [25, 50]] - 2 * [[0, 0], [0, 25]]` = `[[0, 25], [25, 0]]`\n`np.sqrt` -> `[[0, 5], [5, 0]]`. Correct.\n\nWhat about floating point issues? `np.maximum(dist_sq, 0.0)` handles it.\nIs there any case where `points` might be 1D? The problem says 2-D array of shape `(n, d)`.\nThe function signature is given.\nI'll wrap it in the required format.\n\nOne minor point: `np.einsum` or `np.dot` could be used, but `@` is fine.\nAnother approach: `scipy.spatial.distance.cdist` is not allowed since we should use numpy.\nThe broadcasting approach `np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(axis=2))` is also fine and maybe more robust to floating point issues for small values, but uses more memory. Given typical coding challenge constraints, either is fine. I'll stick to the dot product one as it's standard and efficient.\n\nWait, should I handle the case where `points` might have a different dtype? The problem says `np.ndarray`, usually float. I'll assume float or convert to float if necessary, but not strictly required.\n\nLet's refine:\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, np.newaxis] + sq_norms[np.newaxis, :] - 2 * np.dot(points, points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```\nUsing `np.newaxis` is slightly more explicit than `None`, but both work. I'll stick to `None` or `np.newaxis`.\n\nCheck constraints: \"Use NumPy vectorization (no Python `for` loops over point pairs).\" - satisfied.\n\"Return an `(n, n)` array\" - satisfied.\n\"symmetric, with a zero diagonal\" - satisfied.\n\nOne thing: `np.dot(points, points.T)` might be slower than `points @ points.T` in some numpy versions, but both are fine. I'll use `points @ points.T`.\n\nLet's consider potential issues with large arrays: `dist_sq` could have negative values due to precision. `np.maximum(dist_sq, 0.0)` fixes it.\nIs there any chance `points` is integer? `points**2` will be integer, but distances should be float. `np.sqrt` will cast to float. It's fine.\n\nI will write the solution exactly as requested.\n\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```\n\nWait, should I ensure the output is float? `np.sqrt` returns float64 if input is float, or float if input is int. It's fine.\n\nLet's double check the example: `[[0, 0], [3, 4]]` -> `[[0, 5], [5, 0]]`.\nMy code produces exactly that.\n\nOne minor point: sometimes `points @ points.T` can be numerically unstable for very large coordinates, but for typical problems it's fine. The broadcasting method `np.sqrt(np.sum((points[:, None, :] - points[None, :, :])**2, axis=2))` is more numerically stable for differences but uses more memory. Given the prompt doesn't specify scale, either is acceptable. I'll stick to the dot product method as it's standard for pairwise distances. Actually, to be safe against floating point issues, I'll use `np.maximum(dist_sq, 0.0)`.\n\nAnother option: `scipy.spatial.distance.pdist` is not numpy.\nI'll stick to the dot product method.\n\nFinal code:\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```\n\nWait, `points @ points.T` might be interpreted as matrix multiplication. Yes.\nIs there any reason to prefer `np.einsum('ij,kj->ik', points, points)`? It's equivalent to `points @ points.T`.\nI'll keep it simple.\n\nOne detail: `np.maximum` might return an array with the same dtype as `dist_sq`. If `dist_sq` is integer, `np.maximum` keeps it integer, but `np.sqrt` will convert to float. That's fine.\n\nLe\n…[765 chars elided]…\n_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```\nThis is correct and efficient.\nI will wrap it in the required format.\nDone. \nSelf-Correction/Verification during thought:\nIs there any case where `points @ points.T` fails? Only if `points` is not 2D, but problem says it is.\nCould `dist_sq` have significant negative values due to precision? `np.maximum` handles it.\nIs `np.maximum` vectorized? Yes.\nAll good.\nOutput matches request.\nProceed. \nGenerating response. \n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```\nLooks perfect.✅\nI will output this. \nNote: I'll make sure indentation and formatting are correct.\nThe prompt asks for `solution.py` content. I'll provide just the code block.\nDone. \nOutput matches. \nProceeds. \n[Final Check]\n- Function name: `pairwise_distances`\n- Signature: `def pairwise_distances(points: np.ndarray) -> np.ndarray:`\n- Uses numpy vectorization\n- Returns `(n, n)` array\n- Symmetric, zero diagonal\n- Example works\nAll constraints met.\nReady. \nOutput generation. \n(Self-Correction/Refinement during thought)\nActually, `points @ points.T` might be slightly slower than `np.dot(points, points.T)` in older numpy, but it's fine. I'll stick with `@`.\nAnother point: `np.sum(points**2, axis=1)` creates a temporary array `points**2`. For memory efficiency, `np.einsum('ij,ij->i', points, points)` is better, but not strictly necessary. I'll keep it simple.\nThe code is good.\nProceeds. \nOutput matches response. \n[Done] \n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    return np.sqrt(np.maximum(dist_sq, 0.0))\n```\nAll good. \nOutput matches. \nProceeds. \n[Final Output Generation] -> *Proceeds*\n",
    "test_error": ".F.                                                                      [100%]\n=================================== FAILURES ===================================\n_________________________ test_symmetric_zero_diagonal _________________________\n\n    def test_symmetric_zero_diagonal():\n        rng = np.random.default_rng(0)\n        pts = rng.standard_normal((6, 3))\n        d = pairwise_distances(pts)\n        assert d.shape == (6, 6)\n        np.testing.assert_allclose(d, d.T, atol=1e-9)\n>       np.testing.assert_allclose(np.diag(d), np.zeros(6), atol=1e-9)\nE       AssertionError: \nE       Not equal to tolerance rtol=1e-07, atol=1e-09\nE       \nE       Mismatched elements: 1 / 6 (16.7%)\nE       Mismatch at index:\nE        [5]: 1.4901161193847656e-08 (ACTUAL), 0.0 (DESIRED)\nE       Max absolute difference among violations: 1.49011612e-08\nE       Max relative difference among violations: inf\nE        ACTUAL: array([0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,\nE              0.000000e+00, 1.490116e-08])\nE        DESIRED: array([0., 0., 0., 0., 0., 0.])\n\ntest_distances.py:19: AssertionError\n=========================== short test summary info ============================\nFAILED test_distances.py::test_symmetric_zero_diagonal - AssertionError: \n1 failed, 2 passed in 0.06s"
  },
  {
    "total": 3,
    "answer": "```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    dist = np.sqrt(np.maximum(dist_sq, 0.0))\n    np.fill_diagonal(dist, 0.0)\n    return dist\n```",
    "passed": 3,
    "reasoning": "The test failure shows that the diagonal elements are not exactly zero due to floating-point precision errors.\nThe maximum absolute difference is `1.49e-08`, which is slightly above the `atol=1e-9` tolerance.\nTo fix this, I can explicitly set the diagonal to zero after computing the distances.\nAlternatively, I can use `np.fill_diagonal(d, 0.0)` or `np.diag_indices_from`.\n\nLet's modify the function to ensure the diagonal is exactly zero.\n\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    dist = np.sqrt(np.maximum(dist_sq, 0.0))\n    np.fill_diagonal(dist, 0.0)\n    return dist\n```\n\nWait, `np.fill_diagonal` modifies in-place. That's fine.\nAlternatively, `dist[np.arange(n), np.arange(n)] = 0.0`.\n\nLet's check if there are any other issues. The problem says \"Use NumPy vectorization (no Python for loops over point pairs).\" My solution uses vectorization.\nThe diagonal fix is standard for this kind of numerical instability.\n\nLet's verify the code:\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    dist = np.sqrt(np.maximum(dist_sq, 0.0))\n    np.fill_diagonal(dist, 0.0)\n    return dist\n```\n\nThis should pass the test. The tolerance `atol=1e-9` is strict, so forcing the diagonal to 0 is the correct approach.\nI will output the complete corrected file.\nOne minor point: `np.fill_diagonal` works on the array in place. It's fine.\nAnother option: `dist[np.diag_indices_from(dist)] = 0.0`\nBoth are fine.\n\nI'll stick with `np.fill_diagonal(dist, 0.0)`.\n\nLet's double check the exact function signature requested:\n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    ...\n```\n\nAll good.\nI will produce the final answer.\nOne thing to consider: sometimes `scipy.spatial.distance.pdist` and `squareform` are used, but the prompt says \"Use NumPy vectorization\". So my approach is correct.\nThe code is ready.\nI will wrap it in the required format.\nDone. \nLet's generate the response. \nSelf-Correction/Verification during thought:\nIs there any chance `np.maximum(dist_sq, 0.0)` causes issues? No, it prevents negative values due to floating point errors before sqrt.\nSetting diagonal to 0 explicitly handles the remaining precision error on the diagonal.\nThe test expects `atol=1e-9`, so `1.49e-08` fails. Setting to 0 fixes it.\nCode looks solid.\nOutput matches request.\nProceed. \nOutput generation. \n```python\nimport numpy as np\n\ndef pairwise_distances(points: np.ndarray) -> np.ndarray:\n    sq_norms = np.sum(points**2, axis=1)\n    dist_sq = sq_norms[:, None] + sq_norms[None, :] - 2 * (points @ points.T)\n    dist = np.sqrt(np.maximum(dist_sq, 0.0))\n    np.fill_diagonal(dist, 0.0)\n    return dist\n```\nThis is correct. I'll present it. \nAll steps verified. No extra fluff.\nDone. \nOutput matches. \nProceeds. \n[Final Check]\nThe prompt asks for `solution.py` with the function.\nI will provide exactly that.\nReady. \nOutput Generation. \n(Self-Correction/Note during generation prep)\nWait, `np.fill_diagonal` modifies in-place, which is fine. But to be strictly functional or avoid side effects if `dist` is shared (it isn't here), it's fine.\nAnother way: `dist[np.arange(len(points)), np.arange(len(points))] = 0.0`\nBoth are fine. I'll stick with `np.fill_diagonal`.\nAll good. \nProceeds. \n[Output Generation] -> *Proceeds*\n",
    "test_error": ""
  }
]