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.

To solve this problem, we need to compute the pairwise Euclidean distances between all points in a given 2D array. The result should be a symmetric matrix where each element [i, j] represents the distance between the i-th and j-th points. The solution must use NumPy vectorization to avoid explicit Python loops for efficiency.

### Approach
1. **Problem Analysis**: Given an array of points, each point is represented as a row in a 2D array. The goal is to compute the Euclidean distance between every pair of points and return these distances in a symmetric matrix with zeros on the diagonal.

2. **Intuition**: The Euclidean distance between two points in d-dimensional space can be computed using the norm of the difference vector between the points. By leveraging NumPy's broadcasting and vectorization, we can efficiently compute these distances for all pairs without explicit loops.

3. **Algorithm**:
   - **Difference Matrix**: Create a 3D array where each element [i, j] is the difference vector between the i-th and j-th points.
   - **Norm Calculation**: Compute the Euclidean norm along the third axis of the difference matrix to get the pairwise distances.

4. **Complexity**: The approach efficiently computes the distances using NumPy operations, which are implemented in optimized C code, making the solution both concise and efficient.

### Solution Code
```python
import numpy as np

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    diff = points[:, np.newaxis, :] - points[np.newaxis, :, :]
    distances = np.linalg.norm(diff, axis=2)
    return distances
```

### Explanation
- **Difference Matrix**: The expression `points[:, np.newaxis, :] - points[np.newaxis, :, :]` creates a 3D array where each element [i, j] is the vector difference between the i-th and j-th points. This is done using NumPy's broadcasting feature.
- **Norm Calculation**: The `np.linalg.norm` function is used with `axis=2` to compute the Euclidean norm along the third axis of the difference matrix, resulting in a 2D array of distances.

This approach ensures that the solution is both efficient and concise, leveraging NumPy's capabilities to handle vectorized operations and broadcasting.

Okay, I need to write a Python function called pairwise_distances that takes a numpy array of points and returns a symmetric matrix of their pairwise Euclidean distances. Hmm, let's think about how to approach this without using for loops, since the problem specifies to use NumPy vectorization.

So, the input is a 2D array where each row is a point, and each column is a dimension. For example, if the input is [[0,0], [3,4]], the output should be [[0,5], [5,0]]. 

I remember that in NumPy, there's a way to compute pairwise distances efficiently. Oh right, there's a function called np.linalg.norm that can compute the norm of an array. But how do I apply this to all pairs of points?

Wait, another approach: for each point, subtract every other point and compute the norm of the difference. But how to do this efficiently.

Let me think about the structure. If I have points as a matrix of shape (n, d), I can compute the difference matrix by subtracting each point from every other. But how?

Oh, right! I can use broadcasting. If I have points1 as points[:, np.newaxis, :], and points2 as points[np.newaxis, :, :], then subtracting them would give a 3D array where each element [i,j] is the difference between point i and point j. Then, I can compute the norm along the third axis.

Wait, but that might be a bit memory-intensive for large n, but given the problem constraints, it's acceptable.

So the steps would be:

1. Compute the difference between each pair of points. This can be done by subtracting the matrix from its transpose in a way that creates a 3D array.

Wait, no, more accurately, I can create two views of the points array: one where each point is repeated along the second axis (points[:, np.newaxis, :]) and another where each point is repeated along the first axis (points[np.newaxis, :, :]). Then subtracting these gives a 3D array of shape (n, n, d), where each element is the difference between point i and point j.

2. Compute the Euclidean norm along the third axis (the dimensions) of this 3D array. This will give a 2D array of shape (n, n), where each element [i,j] is the distance between points i and j.

3. Since the distance matrix is symmetric, we can just compute it once, but in this approach, it's automatically symmetric because the subtraction is done for all i and j.

4. The diagonal elements will be zero because the distance from a point to itself is zero.

So putting this together, the code would be something like:

Compute the difference matrix:
diff = points[:, np.newaxis, :] - points[np.newaxis, :, :]

Then compute the norm along axis 2:
distances = np.linalg.norm(diff, axis=2)

Wait, but wait: np.linalg.norm returns a 2D array, which is exactly what we need.

So the function can be written as:

def pairwise_distances(points: np.ndarray) -> np.ndarray:
    n, d = points.shape
    diff = points[:, np.newaxis, :] - points[np.newaxis, :, :]
    distances = np.linalg.norm(diff, axis=2)
    return distances

Let me test this with the sample input.

Sample input: points = np.array([[0,0], [3,4]])
Compute diff:
points[:, np.newaxis, :] is shape (2,1,2)
points[np.newaxis, :, :] is shape (1,2,2)
Subtracting gives a 3D array of shape (2,2,2):
For i=0, j=0: [0,0] - [0,0] = [0,0]
i=0, j=1: [0,0] - [3,4] = [-3,-4]
i=1, j=0: [3,4] - [0,0] = [3,4]
i=1, j=1: [3,4] - [3,4] = [0,0]

Then, taking the norm along axis 2:
For [0,0], norm is 0.
For [-3,-4], norm is 5.
For [3,4], norm is 5.
For [0,0], norm is 0.

So the resulting distances matrix is:
[[0,5],
 [5,0]]

Which matches the sample output. So this approach works.

Another test case: what if there are more points?

Suppose points = [[1,2], [4,5], [7,8]]
The distances should be:
Between (1,2) and (4,5): sqrt( (3)^2 + (3)^2 ) = sqrt(18) ≈4.2426
Between (1,2) and (7,8): sqrt(6^2 +6^2)=sqrt(72)=~8.485
Between (4,5) and (7,8): sqrt(3^2 +3^2)=sqrt(18)=~4.2426

So the matrix should be:
[
 [0, 4.2426, 8.485],
 [4.2426, 0, 4.2426],
 [8.485, 4.2426, 0]
]

Let me compute using the code:

points = np.array([[1,2], [4,5], [7,8]])
diff = points[:, np.newaxis, :] - points[np.newaxis, :, :]
diff will have shape (3,3,2)
Each element [i,j] is points[i] - points[j]

Then, distances = np.linalg.norm(diff, axis=2)

Yes, that should give the correct matrix.

So the code seems to handle all cases correctly.

I think this should work. Now, I'll write the function as described.