all these correlation functions

[kushalkolar]

There are a lot of correlation functions in caiman.summary_images and it's not easy to understand them from the names and docstrings. I think it would be important to make a comparison that considers:

• performance
• some of them seem to work directly on array objects, such as correlation_pnr, which would allow subsampling. Some take in files and have arguments like window and stride and seem to compute them with parallelization, but why? In the CNMFE demo nb it says that subsampling for corr-pnr gives the same results, so what are these patch-parallelized functions for?
• What use cases are each of these suitable for? I think an example gallery with this module would help.

[pgunn]

Usage of all functions in that namespace:

• max_correlation_image(): Unused
• local_correlations_fft(): Used internally in a few places
• local_correlations_multicolor(): Used internally by local_correlations_movie_parallel()
• local_correlations(): used by notebooks/demos as well as internally
• correlation_pnr(): used by notebooks/demos
• iter_chunk_array(): used internally
• correlation_image_ecobost(): Unused
• map_coorr(): map function used internally
• prepare_local_correlations(): Used internally
• update_local_correlations(): Used internally
• local_correlations_movie(): Unused?
• local_correlations_movie_offline(): Used by demos/notebooks
• local_correlations_movie_parallel(): Used internally
• mean_image(): Used by notebooks/demos
• mean_image_parallel(): Used internally

[ryanirl]

To add to the chaos, here's a pretty significant upgrade (speed and memory) to these functions that I wrote some time ago for my research. You might have to look over it a bit to verify it's validity mathematically, but huge gains in terms of speed and memory.

It's also pretty easy to generalize this optimization idea to other local-correlation-like functions.

from tqdm.auto import tqdm
import numpy as np

from scipy.ndimage import convolve

def local_correlation_projection(frames: np.ndarray, kernel_size: int = 5, *, verbose: bool = True) -> np.ndarray:
    """Memory and time efficient local correlation algorithm. It makes use of
    the fact that a correlation between two signals is a normalized dot product,
    and that the dot product is an inherently iterative procedure. From the dot
    product perspective, a correlation between one pixel's temporal trace and 
    anothers is just a running sum of their time points multiplied by each other
    after mean and unit-one normalization. Thus a local correlation is just this
    operation applied with every pixel and it's neighbors within some pre-defined 
    radius, followed by a mean at the end. This can further be optimized when you
    realize that a convolution with a 0 center uniform filter was just described.

    The extra space required for this algorithm is roughly `frames[i] * 5` for 
    some arbitrary `i`. Unlike other algorithms, this implementation does not
    scale with the size of `frames`, but instead with the size of each individual 
    frame given that it's iterative.

    In my testing, for a 20 gigabyte video, running this function didn't raise 
    my RAM over extra gigabyte using a kernel size of 5. Furthermore, it was
    faster than any other implementation I have used. 

    Args:
        frames (np.ndarray): A numpy movie of shape (n_frames, height, width).
        kernel_size (int): The size in diameter (not radius) of the nearest 
            neighbor kernel to use. 
        verbose (bool): Whether to show the TQDM loading bars or not.

    Returns:
        np.ndarray: A two-dimensional array representing the local correlation
            projection of `frames`.

    """
    if kernel_size % 2 == 0:
        raise ValueError("'kernel_size' must be odd.")

    n_frames, height, width = frames.shape

    # Compute the kernel for the per-frame convolution. 
    kernel = np.ones((kernel_size, kernel_size), dtype = np.float32)
    kernel[kernel_size // 2, kernel_size // 2] = 0
    kernel = kernel / kernel.sum()

    mean = frames.mean(axis = 0)

    # Using np.linalg.norm(frames, axis = 0) is not memory efficient so we
    # vectorize it only over one dimension.
    norm = np.zeros((height, width), dtype = np.float32) 
    for i in tqdm(range(frames.shape[1]), disable = not verbose):
        norm[i] = np.linalg.norm(frames[:, i] - mean[i], axis = 0)

    # Now compute the local correlation as a running sum of normalized dot
    # products of each pixel with it's nearest neighbors.
    local_corr = np.zeros((height, width), dtype = np.float32)
    for i in tqdm(range(n_frames), disable = not verbose):
        image = (frames[i].copy() - mean) / norm
        image = image * convolve(image, kernel, mode = "constant")
        local_corr += image

    return local_corr

[pgunn]

Thanks; we'd need to stare at it a bit more and test it out a bit, but if it holds up then giving it as an option could be reasonable, particularly if it's faster and more memory-efficient.