def harmonic_aggregation(window_scores: torch.Tensor, output_size: tuple, masks: torch.Tensor) -> torch.Tensor:
    """Perform harmonic aggregation on window scores.

    Computes a single score for each patch location by aggregating the scores of all
    windows that cover the patch. Scores are aggregated using the harmonic mean.

    Args:
        window_scores (torch.Tensor): Scores for each sliding window location.
            Shape: ``(batch_size, n_masks)``.
        output_size (tuple): Output dimensions ``(H, W)``.
        masks (torch.Tensor): Binary masks indicating which patches are covered by each
            window. Shape: ``(n_patches_per_mask, n_masks)``.

    Returns:
        torch.Tensor: Aggregated scores. Shape: ``(batch_size, H, W)``.

    Example:
        Example for a 3x3 patch grid with 4 sliding windows of size 2x2:

        >>> window_scores = torch.tensor([[1.0, 0.75, 0.5, 0.25]])
        >>> output_size = (3, 3)
        >>> masks = torch.Tensor([[0, 1, 3, 4],
        ...                      [1, 2, 4, 5],
        ...                      [3, 4, 6, 7],
        ...                      [4, 5, 7, 8]])
        >>> harmonic_aggregation(window_scores, output_size, masks)
        tensor([[[1.0000, 0.8571, 0.7500],
                 [0.6667, 0.4800, 0.3750],
                 [0.5000, 0.3333, 0.2500]]])

    Note:
        The harmonic mean is used instead of arithmetic mean as it is more sensitive to
        low scores, making it better suited for anomaly detection where we want to
        emphasize potential defects.
    """
    batch_size = window_scores.shape[0]
    height, width = output_size

    scores = []
    for idx in range(height * width):
        patch_mask = torch.any(masks == idx, dim=0)  # boolean tensor indicating which masks contain the patch
        scores.append(sum(patch_mask) / (1 / window_scores.T[patch_mask]).sum(dim=0))

    return torch.stack(scores).T.reshape(batch_size, height, width).nan_to_num(posinf=0.0)
Perform harmonic aggregation on window scores.

    Computes a single score for each patch location by aggregating the scores of all
    windows that cover the patch. Scores are aggregated using the harmonic mean.

    Args:
        window_scores (torch.Tensor): Scores for each sliding window location.
            Shape: ``(batch_size, n_masks)``.
        output_size (tuple): Output dimensions ``(H, W)``.
        masks (torch.Tensor): Binary masks indicating which patches are covered by each
            window. Shape: ``(n_patches_per_mask, n_masks)``.

    Returns:
        torch.Tensor: Aggregated scores. Shape: ``(batch_size, H, W)``.

    Example:
        Example for a 3x3 patch grid with 4 sliding windows of size 2x2:

        >>> window_scores = torch.tensor([[1.0, 0.75, 0.5, 0.25]])
        >>> output_size = (3, 3)
        >>> masks = torch.Tensor([[0, 1, 3, 4],
        ...                      [1, 2, 4, 5],
        ...                      [3, 4, 6, 7],
        ...                      [4, 5, 7, 8]])
        >>> harmonic_aggregation(window_scores, output_size, masks)
        tensor([[[1.0000, 0.8571, 0.7500],
                 [0.6667, 0.4800, 0.3750],
                 [0.5000, 0.3333, 0.2500]]])

    Note:
        The harmonic mean is used instead of arithmetic mean as it is more sensitive to
        low scores, making it better suited for anomaly detection where we want to
        emphasize potential defects.
    
src/anomalib/models/image/winclip/utils.py