  1# -*- coding: utf-8 -*-
     2"""
     3Created on Thu Aug  1 14:18:42 2024
     4
     5@author: Bart Steemans. Govers Lab.
     6"""
     7
     8import pandas as pd
     9import numpy as np
    10from scipy import stats
    11from skimage import filters
    12from scipy.spatial.distance import pdist, squareform
    13import logging
    14
    15bactoscoop_logger = logging.getLogger("logger")
    16class SignalCorrelation:
    17    def __init__(
    18        self,
    19        df,
    20        channel1,
    21        channel2,
    22        feature,
    23        method_name,
    24        prepared_feature=None,
    25    ):
    26        self.df = df
    27        self.channel1 = channel1
    28        self.channel2 = channel2
    29        self.feature = feature
    30        self.method_name = method_name
    31        self.prepared_feature = prepared_feature
    32        self.method = self.get_method()
    33        self.validate_inputs()
    34
    35    def get_method(self):
    36        methods = {
    37            "manders": self.manders_overlap_coefficient,
    38            "pearson": self.pearson_correlation_coefficient,
    39            "li_icq": self.li_icq,
    40            "ratio": self.ratio,
    41            "spearman": self.spearman_rank_correlation,
    42            "kendall": self.kendall_tau,
    43            "distance_corr": self.distance_correlation,
    44            "covariance": self.covariance,
    45            "n_cross_corr": self.normalized_cross_correlation,
    46            "entropy_diff": self.entropy_difference,
    47            "kurtosis_ratio": self.kurtosis_ratio,
    48            "skewness_product": self.skewness_product,
    49            "zero_crossings_diff": self.zero_crossings_difference,
    50            "fft_peak_ratio": self.fft_peak_ratio,
    51            "fft_energy_ratio": self.fft_energy_ratio,
    52            "histogram_intersection": self.histogram_intersection,
    53            "cosine_similarity": self.cosine_similarity,
    54        }
    55        if self.method_name in methods:
    56            return methods[self.method_name]
    57        else:
    58            raise ValueError(f'Method "{self.method_name}" is not recognized')
    59
    60    def validate_inputs(self):
    61        feature1 = f"{self.channel1}_{self.feature}"
    62        feature2 = f"{self.channel2}_{self.feature}"
    63        if feature1 not in self.df.columns or feature2 not in self.df.columns:
    64            raise ValueError(
    65                f"Features {feature1} and {feature2} must be in the DataFrame columns"
    66            )
    67
    68        if not callable(self.method):
    69            raise ValueError(f'The method "{self.method_name}" is not callable')
    70
    71    def calculate(self):
    72        prepared_feature = self.prepared_feature
    73        if prepared_feature is None:
    74            prepared_feature = self.prepare_feature_pair(
    75                self.df, self.channel1, self.channel2, self.feature
    76            )
    77
    78        signals1 = prepared_feature["signals1"]
    79        signals2 = prepared_feature["signals2"]
    80        valid_mask = prepared_feature["valid_mask"]
    81
    82        results = [np.nan] * len(signals1)
    83        for index, (signal1, signal2, is_valid) in enumerate(
    84            zip(signals1, signals2, valid_mask)
    85        ):
    86            try:
    87                if is_valid:
    88                    results[index] = self.method(signal1, signal2)
    89            except ValueError as e:
    90                bactoscoop_logger.debug(f"Encountered a ValueError: {e}")
    91                continue
    92
    93        method_name = self.method.__name__
    94        new_feature_name = (
    95            f"{self.channel1}_{self.channel2}_{self.feature}_{method_name}"
    96        )
    97        self.df[new_feature_name] = results
    98        return self.df
    99
   100    @classmethod
   101    def prepare_feature_pair(cls, df, channel1, channel2, feature):
   102        feature1 = f"{channel1}_{feature}"
   103        feature2 = f"{channel2}_{feature}"
   104        signals1 = df[feature1].to_numpy(copy=False)
   105        signals2 = df[feature2].to_numpy(copy=False)
   106        valid_mask = np.fromiter(
   107            (
   108                cls.is_valid_signal(signal1) and cls.is_valid_signal(signal2)
   109                for signal1, signal2 in zip(signals1, signals2)
   110            ),
   111            dtype=bool,
   112            count=len(signals1),
   113        )
   114        return {"signals1": signals1, "signals2": signals2, "valid_mask": valid_mask}
   115
   116    @staticmethod
   117    def is_valid_signal(signal):
   118        try:
   119            signal_array = np.asarray(signal, dtype=float)
   120        except (TypeError, ValueError):
   121            return False
   122        return bool(signal_array.size and np.isfinite(signal_array).all())
   123
   124    @staticmethod
   125    def ratio(signal1, signal2):
   126        if (
   127            np.isscalar(signal1)
   128            and np.isscalar(signal2)
   129            and np.isfinite(signal1)
   130            and np.isfinite(signal2)
   131            and signal2 != 0
   132        ):
   133            return signal1 / signal2
   134        return np.nan
   135
   136    @staticmethod
   137    def manders_overlap_coefficient(signal1, signal2):
   138        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   139            signal2, (list, np.ndarray)
   140        ):
   141            if len(signal1) > 1 and len(signal2) > 1:
   142                signal1 = np.array(signal1)
   143                signal2 = np.array(signal2)
   144                return np.sum(signal1 * signal2) / np.sqrt(
   145                    np.sum(signal1**2) * np.sum(signal2**2)
   146                )
   147        return np.nan
   148
   149    @staticmethod
   150    def pearson_correlation_coefficient(signal1, signal2):
   151        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   152            signal2, (list, np.ndarray)
   153        ):
   154            if len(signal1) > 1 and len(signal2) > 1:
   155                return stats.pearsonr(signal1, signal2)[0]
   156        return np.nan
   157
   158    @staticmethod
   159    def li_icq(signal1, signal2):
   160        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   161            signal2, (list, np.ndarray)
   162        ):
   163            if len(signal1) > 1 and len(signal2) > 1:
   164                signal1 = np.array(signal1)
   165                signal2 = np.array(signal2)
   166                mean1, mean2 = np.mean(signal1), np.mean(signal2)
   167                product = ((signal1 - mean1) * (signal2 - mean2)) > 0
   168                return (np.sum(product) / len(signal1) - 0.5) * 2
   169        return np.nan
   170
   171    @staticmethod
   172    def spearman_rank_correlation(signal1, signal2):
   173        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   174            signal2, (list, np.ndarray)
   175        ):
   176            if len(signal1) > 1 and len(signal2) > 1:
   177                return stats.spearmanr(signal1, signal2)[0]
   178        return np.nan
   179
   180    @staticmethod
   181    def kendall_tau(signal1, signal2):
   182        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   183            signal2, (list, np.ndarray)
   184        ):
   185            if len(signal1) > 1 and len(signal2) > 1:
   186                corr, _ = stats.kendalltau(signal1, signal2)
   187                return corr
   188        return np.nan
   189
   190    @staticmethod
   191    def distance_correlation(signal1, signal2):
   192        """Biased sample distance correlation; undefined profiles return NaN."""
   193        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   194            signal2, (list, np.ndarray)
   195        ):
   196            if len(signal1) > 1 and len(signal2) > 1:
   197                x = np.asarray(signal1, dtype=float)
   198                y = np.asarray(signal2, dtype=float)
   199                if (
   200                    x.ndim != 1
   201                    or y.ndim != 1
   202                    or len(x) != len(y)
   203                    or not np.isfinite(x).all()
   204                    or not np.isfinite(y).all()
   205                ):
   206                    return np.nan
   207
   208                # Biased sample distance correlation uses separately centered
   209                # pairwise distances for each signal.
   210                distance_x = squareform(pdist(x[:, None]))
   211                distance_y = squareform(pdist(y[:, None]))
   212                centered_x = (
   213                    distance_x
   214                    - distance_x.mean(axis=0)
   215                    - distance_x.mean(axis=1)[:, None]
   216                    + distance_x.mean()
   217                )
   218                centered_y = (
   219                    distance_y
   220                    - distance_y.mean(axis=0)
   221                    - distance_y.mean(axis=1)[:, None]
   222                    + distance_y.mean()
   223                )
   224                variance_x = np.mean(centered_x**2)
   225                variance_y = np.mean(centered_y**2)
   226                if variance_x == 0 or variance_y == 0:
   227                    return np.nan
   228                correlation_squared = np.mean(centered_x * centered_y) / np.sqrt(
   229                    variance_x * variance_y
   230                )
   231                return float(np.sqrt(np.clip(correlation_squared, 0, 1)))
   232        return np.nan
   233
   234    @staticmethod
   235    def covariance(signal1, signal2):
   236        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   237            signal2, (list, np.ndarray)
   238        ):
   239            if len(signal1) > 1 and len(signal2) > 1:
   240                return np.cov(signal1, signal2)[0, 1]
   241        return np.nan
   242
   243    @staticmethod
   244    def normalized_cross_correlation(signal1, signal2):
   245        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   246            signal2, (list, np.ndarray)
   247        ):
   248            if len(signal1) > 1 and len(signal2) > 1:
   249                return np.correlate(
   250                    signal1 - np.mean(signal1), signal2 - np.mean(signal2)
   251                )[0] / (len(signal1) * np.std(signal1) * np.std(signal2))
   252        return np.nan
   253
   254    @staticmethod
   255    def entropy_difference(signal1, signal2):
   256        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   257            signal2, (list, np.ndarray)
   258        ):
   259            if len(signal1) > 1 and len(signal2) > 1:
   260                entropy1 = stats.entropy(signal1)
   261                entropy2 = stats.entropy(signal2)
   262                return np.abs(entropy1 - entropy2)
   263        return np.nan
   264
   265    @staticmethod
   266    def kurtosis_ratio(signal1, signal2):
   267        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   268            signal2, (list, np.ndarray)
   269        ):
   270            if len(signal1) > 1 and len(signal2) > 1:
   271                kurtosis1 = stats.kurtosis(signal1)
   272                kurtosis2 = stats.kurtosis(signal2)
   273                return np.abs(kurtosis1 / kurtosis2) if kurtosis2 != 0 else np.nan
   274        return np.nan
   275
   276    @staticmethod
   277    def skewness_product(signal1, signal2):
   278        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   279            signal2, (list, np.ndarray)
   280        ):
   281            if len(signal1) > 1 and len(signal2) > 1:
   282                skew1 = stats.skew(signal1)
   283                skew2 = stats.skew(signal2)
   284                return skew1 * skew2
   285        return np.nan
   286
   287    @staticmethod
   288    def zero_crossings_difference(signal1, signal2):
   289        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   290            signal2, (list, np.ndarray)
   291        ):
   292            if len(signal1) > 1 and len(signal2) > 1:
   293                zero_crossings1 = np.sum(np.diff(np.sign(signal1)) != 0)
   294                zero_crossings2 = np.sum(np.diff(np.sign(signal2)) != 0)
   295                return np.abs(zero_crossings1 - zero_crossings2)
   296        return np.nan
   297
   298    @staticmethod
   299    def fft_peak_ratio(signal1, signal2):
   300        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   301            signal2, (list, np.ndarray)
   302        ):
   303            if len(signal1) > 1 and len(signal2) > 1:
   304                fft1 = np.abs(np.fft.fft(signal1))
   305                fft2 = np.abs(np.fft.fft(signal2))
   306                peak1 = np.max(fft1[1:])  # Exclude DC component
   307                peak2 = np.max(fft2[1:])  # Exclude DC component
   308                return peak1 / peak2 if peak2 != 0 else np.nan
   309        return np.nan
   310
   311    @staticmethod
   312    def fft_energy_ratio(signal1, signal2):
   313        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   314            signal2, (list, np.ndarray)
   315        ):
   316            if len(signal1) > 1 and len(signal2) > 1:
   317                # Compute the FFT of the signals
   318                fft1 = np.fft.fft(signal1)
   319                fft2 = np.fft.fft(signal2)
   320
   321                # Compute the energy in each frequency bin
   322                energy1 = np.sum(np.abs(fft1) ** 2)
   323                energy2 = np.sum(np.abs(fft2) ** 2)
   324
   325                # Return the energy ratio
   326                return energy1 / energy2 if energy2 != 0 else np.nan
   327        return np.nan
   328
   329    @staticmethod
   330    def histogram_intersection(signal1, signal2, bins=10):
   331        """Intersect probability histograms using shared bins for both profiles."""
   332        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   333            signal2, (list, np.ndarray)
   334        ):
   335            if len(signal1) > 1 and len(signal2) > 1:
   336                x = np.asarray(signal1, dtype=float)
   337                y = np.asarray(signal2, dtype=float)
   338                if (
   339                    x.ndim != 1
   340                    or y.ndim != 1
   341                    or not np.isfinite(x).all()
   342                    or not np.isfinite(y).all()
   343                ):
   344                    return np.nan
   345
   346                # Both histograms must describe the same intensity intervals.
   347                # Normalize each by its own sample count so different profile
   348                # lengths remain comparable and the score is symmetric.
   349                edges = np.histogram_bin_edges(np.concatenate((x, y)), bins=bins)
   350                hist1, _ = np.histogram(x, bins=edges)
   351                hist2, _ = np.histogram(y, bins=edges)
   352                count1 = hist1.sum()
   353                count2 = hist2.sum()
   354                if count1 == 0 or count2 == 0:
   355                    return np.nan
   356                return float(np.minimum(hist1 / count1, hist2 / count2).sum())
   357        return np.nan
   358
   359    @staticmethod
   360    def cosine_similarity(signal1, signal2):
   361        if isinstance(signal1, (list, np.ndarray)) and isinstance(
   362            signal2, (list, np.ndarray)
   363        ):
   364            if len(signal1) > 1 and len(signal2) > 1:
   365                return np.dot(signal1, signal2) / (
   366                    np.linalg.norm(signal1) * np.linalg.norm(signal2)
   367                )
   368        return np.nan