From 582f1a653d140e7b196a14ac40032401153cfbc6 Mon Sep 17 00:00:00 2001
From: Armin Etemad <130641768+ArminGEtemad@users.noreply.github.com>
Date: Thu, 27 Jun 2024 16:04:26 +0200
Subject: [PATCH 1/4] added cross-correlation calculation up to 4 data inputs

---
 src/signalsnap/spectrum_calculator.py | 309 +++++++++++++++++++++-----
 1 file changed, 258 insertions(+), 51 deletions(-)

diff --git a/src/signalsnap/spectrum_calculator.py b/src/signalsnap/spectrum_calculator.py
index 1d5029a..87a476a 100644
--- a/src/signalsnap/spectrum_calculator.py
+++ b/src/signalsnap/spectrum_calculator.py
@@ -759,8 +759,76 @@ def c3(self, a_w1, a_w2, a_w3):
                                                 2 * d_1_mean * d_2_mean * d_3_mean)
         return s3
 
-# ==================== new compact algorithm for c4 =================================
-    def c4(self, a_w, a_w_corr):
+# ==================== c4 for combinations =================================
+    def c4(self, a_w1, a_w2, a_w3, a_w4):
+        """
+            Calculation of c4 for trispectrum based on equation 60 in arXiv:1904.12154.
+
+            Parameters
+            ----------
+            a_w1 : array
+                Fourier coefficients of the signal.
+            a_w2 : array
+                Fourier coefficients of the signal or a second signal.
+            a_w3 : array
+                Fourier coefficients of the signal or a third signal.
+            a_w4 : array
+                Fourier coefficients of the signal or a fourth signal.
+
+            Returns
+            -------
+            s4 : array
+                The c4 estimator as a matrix.
+
+            Notes
+            -----
+            The value of `m`, the number of windows used for the calculation of one spectrum,
+            is obtained from the `config` object associated with this instance.
+            """
+
+        m = self.config.m
+
+        a_w2_conj = conj(a_w2)
+        a_w4_conj = conj(a_w4)
+
+        if self.config.coherent:
+            sum_11c22c = af.matmulNT(a_w1 * a_w2_conj, a_w3 * a_w4_conj)
+            sum_11c22c_m = mean(sum_11c22c, dim=2)
+            s4 = sum_11c22c_m
+        else:
+            a_w1_mean = a_w1 - af.tile(mean(a_w1, dim=2), d0=1, d1=1, d2=a_w1.shape[2])
+            a_w2_conj_mean = a_w2_conj - af.tile(mean(a_w2_conj, dim=2), d0=1, d1=1, d2=a_w2_conj.shape[2])
+            a_w3_mean = a_w3 - af.tile(mean(a_w3, dim=2), d0=1, d1=1, d2=a_w3.shape[2])
+            a_w4_conj_mean = a_w4_conj - af.tile(mean(a_w4_conj, dim=2), d0=1, d1=1, d2=a_w4_conj.shape[2])
+
+            xyzw = af.matmulNT(a_w1_mean * a_w2_conj_mean, a_w3_mean * a_w4_conj_mean)
+            xyzw_mean = mean(xyzw, dim=2)
+            
+            xy_mean = mean(a_w1_mean * a_w2_conj_mean, dim=2)
+            zw_mean = mean(a_w3_mean * a_w4_conj_mean, dim=2)
+            xy_zw_mean = af.matmulNT(xy_mean, zw_mean)
+
+            xz_mean = mean(af.matmulNT(a_w1_mean, a_w3_mean), dim=2)
+            yw_mean = mean(af.matmulNT(a_w2_conj_mean, a_w4_conj_mean), dim=2)
+            xz_yw_mean = xz_mean * yw_mean
+
+            xw_mean = mean(af.matmulNT(a_w1_mean, a_w4_conj_mean), dim=2)
+            yz_mean = mean(af.matmulNT(a_w2_conj_mean, a_w3_mean), dim=2)
+            xw_yz_mean = xw_mean * yz_mean
+
+            s4 = m ** 2 / ((m - 1) * (m - 2) * (m - 3)) * (
+                    (m + 1) * xyzw_mean -
+                    (m - 1) * (
+                            xy_zw_mean + xz_yw_mean + xw_yz_mean
+                    )
+            )
+        return s4
+
+# ======================================================================================
+
+
+# ==================== compact algorithm for c4 =================================
+    def c4_old2(self, a_w, a_w_corr):
         """
             Calculation of c4 for trispectrum based on equation 60 in arXiv:1904.12154.
 
@@ -974,12 +1042,19 @@ def import_data(self):
         else:
             return main_data
 
-    def import_corr_data(self):
+    def import_corr_data(self, which_data):
         """
         Helper function to load data from h5 file into a numpy array.
         Imports data in the .h5 format with structure group_key -> data + attrs[dt].
         Parameters such as 'full_import', 'path', 'group_key', and 'dataset' are obtained from the config attribute.
 
+        Parameters
+        ----------
+        which data:
+            2 for corr_data
+            3 for data3
+            4 for data4
+
         Returns
         -------
         numpy.ndarray
@@ -991,19 +1066,47 @@ def import_corr_data(self):
         -----
         Ensure that the config attribute is properly set before calling this method, and that the paths and keys correspond to existing elements in the h5 file.
         """
-
-        main = h5py.File(self.config.corr_path, 'r')
-        if self.config.corr_group_key == '':
-            main_data = main[self.config.corr_dataset]
-        else:
-            main_group = main[self.config.corr_group_key]
-            main_data = main_group[self.config.corr_dataset]
-        if self.config.delta_t is None:
-            self.config.delta_t = main_data.attrs['dt']
-        if self.config.full_import:
-            return main_data[()]
-        else:
-            return main_data
+        if which_data == 2:
+            main = h5py.File(self.config.corr_path, 'r')
+            if self.config.corr_group_key == '':
+                main_data = main[self.config.corr_dataset]
+            else:
+                main_group = main[self.config.corr_group_key]
+                main_data = main_group[self.config.corr_dataset]
+            if self.config.delta_t is None:
+                self.config.delta_t = main_data.attrs['dt']
+            if self.config.full_import:
+                return main_data[()]
+            else:
+                return main_data
+            
+        elif which_data == 3:
+            main = h5py.File(self.config.path3, 'r')
+            if self.config.group_key3 == '':
+                main_data = main[self.config.dataset3]
+            else:
+                main_group = main[self.config.group_key3]
+                main_data = main_group[self.config.dataset3]
+            if self.config.delta_t is None:
+                self.config.delta_t = main_data.attrs['dt']
+            if self.config.full_import:
+                return main_data[()]
+            else:
+                return main_data
+            
+        elif which_data == 4:
+            main = h5py.File(self.config.path4, 'r')
+            if self.config.group_key4 == '':
+                main_data = main[self.config.dataset4]
+            else:
+                main_group = main[self.config.group_key4]
+                main_data = main_group[self.config.dataset4]
+            if self.config.delta_t is None:
+                self.config.delta_t = main_data.attrs['dt']
+            if self.config.full_import:
+                return main_data[()]
+            else:
+                return main_data
 
     def save_spec(self, save_path, remove_S_stationarity=False):
         """
@@ -1210,8 +1313,8 @@ def calc_overlap(self, t_unit=None, imag=False, scale_t=1):
         return t, t_main, overlap_s2, overlap_s3, overlap_s4
 
     def fourier_coeffs_to_spectra(self, orders, a_w_all_gpu, f_max_ind, f_min_ind,
-                                  single_window, window=None, chunk_corr_gpu=None,
-                                  window_points=None):
+                                  single_window, window=None, a_w_all_corr=None,
+                                  a_w3_all_gpu=None, a_w4_all_gpu=None, window_points=None):
         """
         Helper function to calculate the (1,2,3,4)-order cumulant from the Fourier coefficients of the windows in
         one frame.
@@ -1223,7 +1326,7 @@ def fourier_coeffs_to_spectra(self, orders, a_w_all_gpu, f_max_ind, f_min_ind,
         ----------
         orders : list of {1, 2, 3, 4}
             Orders of the spectra to be calculated.
-        a_w_all_gpu : array_like
+        a_w_all_corr/a_w3_all_gpu/a_w4_all_gpu : array_like
             A matrix containing the Fourier coefficients of the windows.
         f_max_ind : int
             Index of the maximum frequency in the frequency array to calculate the spectral values at.
@@ -1259,45 +1362,90 @@ def fourier_coeffs_to_spectra(self, orders, a_w_all_gpu, f_max_ind, f_min_ind,
                     a_w = a_w_all_gpu
 
                 if self.config.corr_data is not None:
-                    a_w_all_corr = fft_r2c(window * chunk_corr_gpu, dim0=0, scale=1)
+                    #a_w_all_corr = fft_r2c(window * chunk_corr_gpu, dim0=0, scale=1)
                     a_w_corr = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
                     single_spectrum = self.c2(a_w, a_w_corr) / (
                             self.config.delta_t * (single_window ** order).sum())
-
                 else:
                     single_spectrum = self.c2(a_w, a_w) / (
                             self.config.delta_t * (single_window ** order).sum())
-
+                    
             elif order == 3:
-                a_w1 = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind // 2))), dim=0)
-                a_w2 = a_w1
-
-                # a_w3 = to_gpu(calc_a_w3(a_w_all_gpu.to_ndarray(), f_max_ind, self.config.m))
-
-                # ======== New algorithm divides the steps ==========
-                t0 = a_w_all_gpu.to_ndarray()
-                t1 = calc_a_w3(t0, f_max_ind, self.config.m, self.a_w3_init, self.indi, self.config.backend)
-                a_w3 = to_gpu(t1)
-                # ===================================================
-
-                single_spectrum = self.c3(a_w1, a_w2, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                if self.config.corr_data is not None and self.config.data3 is not None:
+                    a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind//2))), dim=0)
+                    a_w2 = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind//2))), dim=0)
+                    t0 = a_w3_all_gpu.to_ndarray()
+                    t1 = calc_a_w3(t0, f_max_ind, self.config.m, self.a_w3_init, self.indi, self.config.backend)
+                    a_w3 = to_gpu(t1)
+                    single_spectrum = self.c3(a_w, a_w2, a_w3) / (self.config.delta_t * (single_window ** order).sum())
 
+                elif self.config.corr_data is not None:
+                    if self.config.combination3 == '2_122':
+                        a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind//2))), dim=0)
+                        a_w2 = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind//2))), dim=0)
+                        t0 = a_w_all_corr.to_ndarray()
+                        t1 = calc_a_w3(t0, f_max_ind, self.config.m, self.a_w3_init, self.indi, self.config.backend)
+                        a_w3 = to_gpu(t1)
+                        single_spectrum = self.c3(a_w, a_w2, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                
+                    elif self.config.combination3 == '2_112':
 
-            else:  # order 4
+                        a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind//2))), dim=0)
+                        t0 = a_w_all_corr.to_ndarray()
+                        t1 = calc_a_w3(t0, f_max_ind, self.config.m, self.a_w3_init, self.indi, self.config.backend)
+                        a_w3 = to_gpu(t1)
+                        single_spectrum = self.c3(a_w, a_w, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                else:
+                    a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind//2))), dim=0)
+                    t0 = a_w_all_gpu.to_ndarray()
+                    t1 = calc_a_w3(t0, f_max_ind, self.config.m, self.a_w3_init, self.indi, self.config.backend)
+                    a_w3 = to_gpu(t1)
+                    single_spectrum = self.c3(a_w, a_w, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+
+            else: # order 4
+                if self.config.corr_data is not None and self.config.data3 is not None and self.config.data4 is not None:
+                    a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    a_w2 = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    a_w3 = af.lookup(a_w3_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    a_w4 = af.lookup(a_w4_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    if self.config.combination4 == '4_1234':
+                        single_spectrum = self.c4(a_w, a_w2, a_w3, a_w4) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '4_1324':
+                        single_spectrum = self.c4(a_w, a_w3, a_w2, a_w4) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '4_1423':
+                        single_spectrum = self.c4(a_w, a_w4, a_w2, a_w3) / (self.config.delta_t * (single_window ** order).sum())
                 
+                elif self.config.corr_data is not None and self.config.data3 is not None:
+                    a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    a_w2 = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    a_w3 = af.lookup(a_w3_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    if self.config.combination4 == '3_1233':
+                        single_spectrum = self.c4(a_w, a_w2, a_w3, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '3_1323':
+                        single_spectrum = self.c4(a_w, a_w3, a_w2, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '3_1322':
+                        single_spectrum = self.c4(a_w, a_w3, a_w2, a_w2) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '3_1223':
+                        single_spectrum = self.c4(a_w, a_w2, a_w2, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '3_2311':
+                        single_spectrum = self.c4(a_w2, a_w3, a_w, a_w) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '3_1213':
+                        single_spectrum = self.c4(a_w, a_w2, a_w, a_w3) / (self.config.delta_t * (single_window ** order).sum())
+                    else:
+                        print('Wrong combination!') # due to reoccuring Heisenbugs if you delete it... the program would crash.. most probably
                 
-                a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
-
-                if self.config.corr_data is not None:
-                    a_w_all_corr = fft_r2c(window * chunk_corr_gpu, dim0=0, scale=1)
-                    if self.config.random_phase:
-                        a_w_all_corr = self.add_random_phase(a_w_all_corr, window_points)
-
-                    a_w_corr = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                elif self.config.corr_data is not None:
+                    a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    a_w2 = af.lookup(a_w_all_corr, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    if self.config.combination4 == '2_1212':
+                        single_spectrum = self.c4(a_w, a_w2, a_w, a_w2) / (self.config.delta_t * (single_window ** order).sum())
+                    elif self.config.combination4 == '2_1122':
+                        single_spectrum = self.c4(a_w, a_w, a_w2, a_w2) / (self.config.delta_t * (single_window ** order).sum())
+                    else:
+                        print('Wrong combination!') # due to reoccuring Heisenbugs if you delete it... the program would crash.. most probably
                 else:
-                    a_w_corr = a_w
-
-                single_spectrum = self.c4(a_w, a_w_corr) / (self.config.delta_t * (single_window ** order).sum())
+                    a_w = af.lookup(a_w_all_gpu, af.Array(list(range(f_min_ind, f_max_ind))), dim=0)
+                    single_spectrum = self.c4(a_w, a_w, a_w, a_w) / (self.config.delta_t * (single_window ** order).sum())
 
             self.store_single_spectrum(single_spectrum, order)
             
@@ -1754,7 +1902,11 @@ def calc_spec(self):
         n_chunks = 0
 
         if self.config.corr_path is not None:
-            self.config.corr_data = self.import_corr_data()
+            self.config.corr_data = self.import_corr_data(which_data=2)
+        if self.config.path3 is not None:
+            self.config.data3 = self.import_corr_data(which_data=3)
+        if self.config.path4 is not None:
+            self.config.data4 = self.import_corr_data(which_data=4)
 
         m, window_points, freq_all_freq, f_max_ind, f_min_ind, n_windows = self.setup_data_calc_spec(orders)
         self.window_points = window_points
@@ -1798,12 +1950,44 @@ def calc_spec(self):
                 if self.config.corr_default == 'white_noise':  # use white noise to check for false correlations
                     chunk_corr = np.random.randn(window_points, 1, m)
                     chunk_corr_gpu = to_gpu(chunk_corr)
+
+                    chunk_a_w3 = np.random.randn(window_points, 1, m)
+                    chunk_a_w3_gpu = to_gpu(chunk_a_w3)
+                    chunk_a_w4 = np.random.randn(window_points, 1, m)
+                    chunk_a_w4_gpu = to_gpu(chunk_a_w4)
+
+                elif self.config.corr_data is not None and self.config.data3 is not None and self.config.data4 is not None:
+                    chunk_corr = self.config.corr_data[int(i * (window_points * m) + self.config.corr_shift): int(
+                        (i + 1) * (window_points * m) + self.config.corr_shift)]
+                    chunk_corr_gpu = to_gpu(chunk_corr.reshape((window_points, 1, m), order='F'))
+
+                    chunk_a_w3 = self.config.data3[int(i * (window_points * m) + self.config.corr_shift): int(
+                        (i + 1) * (window_points * m) + self.config.corr_shift)]
+                    chunk_a_w3_gpu = to_gpu(chunk_a_w3.reshape((window_points, 1, m), order='F'))
+                    chunk_a_w4 = self.config.data4[int(i * (window_points * m) + self.config.corr_shift): int(
+                        (i + 1) * (window_points * m) + self.config.corr_shift)]
+                    chunk_a_w4_gpu = to_gpu(chunk_a_w4.reshape((window_points, 1, m), order='F'))
+
+                elif self.config.corr_data is not None and self.config.data3 is not None:
+                    chunk_corr = self.config.corr_data[int(i * (window_points * m) + self.config.corr_shift): int(
+                        (i + 1) * (window_points * m) + self.config.corr_shift)]
+                    chunk_corr_gpu = to_gpu(chunk_corr.reshape((window_points, 1, m), order='F'))
+
+                    chunk_a_w3 = self.config.data3[int(i * (window_points * m) + self.config.corr_shift): int(
+                        (i + 1) * (window_points * m) + self.config.corr_shift)]
+                    chunk_a_w3_gpu = to_gpu(chunk_a_w3.reshape((window_points, 1, m), order='F'))
+                    chunk_a_w4_gpu = None
+                
                 elif self.config.corr_data is not None:
                     chunk_corr = self.config.corr_data[int(i * (window_points * m) + self.config.corr_shift): int(
                         (i + 1) * (window_points * m) + self.config.corr_shift)]
                     chunk_corr_gpu = to_gpu(chunk_corr.reshape((window_points, 1, m), order='F'))
+                    chunk_a_w3_gpu = None
+                    chunk_a_w4_gpu = None
                 else:
                     chunk_corr_gpu = None
+                    chunk_a_w3_gpu = None
+                    chunk_a_w4_gpu = None
 
                 # ---------count windows-----------
                 n_chunks += 1
@@ -1814,7 +1998,28 @@ def calc_spec(self):
                         np.array(m * [np.ones_like(single_window)]).flatten().reshape((window_points, 1, m), order='F'))
                     a_w_all_gpu = fft_r2c(ones * chunk_gpu, dim0=0, scale=1)
                 else:
-                    a_w_all_gpu = fft_r2c(window * chunk_gpu, dim0=0, scale=1)
+                    if self.config.corr_data is not None and self.config.data3 is not None and self.config.data4 is not None:
+                        a_w_all_gpu = fft_r2c(window * chunk_gpu, dim0=0, scale=1)
+                        a_w_all_corr = fft_r2c(window * chunk_corr_gpu, dim0=0, scale=1)
+                        a_w3_all_gpu = fft_r2c(window * chunk_a_w3_gpu, dim0=0, scale=1)
+                        a_w4_all_gpu = fft_r2c(window * chunk_a_w4_gpu, dim0=0, scale=1)
+
+                    elif self.config.corr_data is not None and self.config.data3 is not None:
+                        a_w_all_gpu = fft_r2c(window * chunk_gpu, dim0=0, scale=1)
+                        a_w_all_corr = fft_r2c(window * chunk_corr_gpu, dim0=0, scale=1)
+                        a_w3_all_gpu = fft_r2c(window * chunk_a_w3_gpu, dim0=0, scale=1)
+                        a_w4_all_gpu = None
+                        
+                    elif self.config.corr_data is not None:
+                        a_w_all_gpu = fft_r2c(window * chunk_gpu, dim0=0, scale=1)
+                        a_w_all_corr = fft_r2c(window * chunk_corr_gpu, dim0=0, scale=1)
+                        a_w3_all_gpu = None
+                        a_w4_all_gpu = None
+                    else:
+                        a_w_all_gpu = fft_r2c(window * chunk_gpu, dim0=0, scale=1)
+                        a_w_all_corr = None
+                        a_w3_all_gpu = None
+                        a_w4_all_gpu = None
 
                 # --------- modify data ---------
                 if self.config.filter_func:
@@ -1828,7 +2033,9 @@ def calc_spec(self):
 
                 # --------- calculate spectra ----------
                 self.fourier_coeffs_to_spectra(orders, a_w_all_gpu, f_max_ind, f_min_ind, single_window,
-                                               window, chunk_corr_gpu=chunk_corr_gpu, window_points=window_points)
+                                               window=window, a_w_all_corr=a_w_all_corr,
+                                               a_w3_all_gpu=a_w3_all_gpu, a_w4_all_gpu=a_w4_all_gpu,
+                                               window_points=window_points)
                 if n_chunks == self.config.break_after:
                     break
 
@@ -2076,7 +2283,7 @@ def plot(self, plot_config: PlotConfig = None):
             self.plot_config = PlotConfig()  # Use default plot configuration
 
         from .spectrum_plotter import SpectrumPlotter
-        plotter = SpectrumPlotter(self, self.plot_config)
+        plotter = SpectrumPlotter(self, self.config, self.plot_config)
         return plotter.plot()
 
     def stationarity_plot(self, plot_config: PlotConfig = None):
@@ -2086,5 +2293,5 @@ def stationarity_plot(self, plot_config: PlotConfig = None):
             self.plot_config = PlotConfig()  # Use default plot configuration
 
         from .spectrum_plotter import SpectrumPlotter
-        plotter = SpectrumPlotter(self, self.plot_config)
+        plotter = SpectrumPlotter(self,  self.plot_config)
         return plotter.stationarity_plot()

From a43c64b477423ccb360833f7b47400b7915a31ea Mon Sep 17 00:00:00 2001
From: Armin Etemad <130641768+ArminGEtemad@users.noreply.github.com>
Date: Thu, 27 Jun 2024 16:05:16 +0200
Subject: [PATCH 2/4] added correct combinations and data inputs

---
 src/signalsnap/spectrum_config.py | 57 ++++++++++++++++++++++++++++---
 1 file changed, 52 insertions(+), 5 deletions(-)

diff --git a/src/signalsnap/spectrum_config.py b/src/signalsnap/spectrum_config.py
index f73618d..d5b2ed3 100644
--- a/src/signalsnap/spectrum_config.py
+++ b/src/signalsnap/spectrum_config.py
@@ -29,14 +29,25 @@ class SpectrumConfig:
 
     Correlation Parameters
     ----------------------
-    corr_path : str, optional
+    (corr data, data3 and data4)
+    corr_path/path3/path4 : str, optional
         Path to h5 file with second signal for correlation. Default is None.
-    corr_group_key : str, optional
+    corr_group_key/group_key3/group_key4 : str, optional
         Group key for h5 file with correlation signal. Default is None.
-    corr_dataset : str, optional
+    corr_dataset/dataset3/dataset4 : str, optional
         Name of the dataset in h5 file with correlation signal. Default is None.
     corr_shift : int, optional
-        Non-negative integer or None. Default is 0.
+        Non-negative integer or None. Default is 0. This shift will be applied for
+        data3 and data4 as well
+    combination3 : str, needed for cross correlations only
+        The combination of data to caluclate bispectrum.
+        examples:
+            2_122 means 2 data sets for correlation using c_3(X, Y, Y) for computations
+            or
+            2_112 means 2 data sets for correlation using c_3(X, X, Y) for computations
+            for three data sets there is only one combination allowed c_3(X, Y, Z)
+    combination4: str, needed for cross correlations only
+        The combination of data to caluclate trispectrum.
 
     Frequency Parameters
     --------------------
@@ -77,6 +88,21 @@ class SpectrumConfig:
     -----
     Ensure that the provided paths, group keys, and dataset names are valid, as this class does not handle file reading.
     """
+    
+    ALLOWED_COMBINATION3 = {'2_122': '1, 2, 2',
+                            '2_112': '1, 1, 2'}
+
+    ALLOWED_COMBINATION4 = {'4_1234': '1, 2, 3, 4',
+                            '4_1324': '1, 3, 2, 4',
+                            '4_1423': '1, 4, 2, 3',
+                            '3_1233': '1, 2, 3, 3',
+                            '3_1323': '1, 3, 2, 3',
+                            '3_1322': '1, 3, 2, 2',
+                            '3_1223': '1, 2, 2, 3',
+                            '3_2311': '2, 3, 1, 1', 
+                            '3_1213': '1, 2, 1, 3',
+                            '2_1212': '1, 2, 1, 2',
+                            '2_1122': '1, 2, 1, 2'}
 
     def __init__(self, path=None, group_key=None, dataset=None, delta_t=None, data=None,
                  corr_data=None, corr_path=None, corr_group_key=None, corr_dataset=None,
@@ -84,7 +110,9 @@ def __init__(self, path=None, group_key=None, dataset=None, delta_t=None, data=N
                  corr_shift=0, filter_func=False, verbose=True, coherent=False, corr_default=None,
                  break_after=1e6, m=10, m_var=10, m_stationarity=None, interlaced_calculation=True,
                  random_phase=False,
-                 rect_win=False, full_import=True, show_first_frame=True, turbo_mode=False):
+                 rect_win=False, full_import=True, show_first_frame=True, turbo_mode=False,
+                 combination3=None, combination4=None, path3=None, group_key3=None, dataset3=None, data3=None,
+                 data4=None, path4=None, group_key4=None, dataset4=None):
 
         if path is not None and not isinstance(path, str):
             raise ValueError("path must be a string or None.")
@@ -173,6 +201,12 @@ def __init__(self, path=None, group_key=None, dataset=None, delta_t=None, data=N
                 if order_in == 'all':
                     order_in = [1, 2, 4]
                     print("Order 3 has been removed from order_in as f_min must be 0 to calculate the bispectrum.")
+        
+        if combination3 is not None and combination3 not in self.ALLOWED_COMBINATION3:
+            raise ValueError(f'Invalid input {combination3}. Must be None or one of {list(self.ALLOWED_COMBINATION3.keys())}.')
+
+        if combination4 is not None and combination4 not in self.ALLOWED_COMBINATION4:
+            raise ValueError(f'Invalid input {combination4}. Must be None or one of {list(self.ALLOWED_COMBINATION4.keys())}.')
 
         self.path = path
         self.group_key = group_key
@@ -206,3 +240,16 @@ def __init__(self, path=None, group_key=None, dataset=None, delta_t=None, data=N
         self.full_import = full_import
         self.show_first_frame = show_first_frame
         self.turbo_mode = turbo_mode
+
+        self.combination3 = combination3
+        self.combination4 = combination4
+
+        self.path3 = path3
+        self.group_key3 = group_key3
+        self.dataset3 = dataset3
+        self.data3 = data3
+
+        self.path4 = path4
+        self.group_key4 = group_key4
+        self.dataset4 = dataset4
+        self.data4 = data4

From 3bf380c07cfd0078ff83fff01e3385654e7aa161 Mon Sep 17 00:00:00 2001
From: Armin Etemad <130641768+ArminGEtemad@users.noreply.github.com>
Date: Thu, 27 Jun 2024 16:06:02 +0200
Subject: [PATCH 3/4] added correct title depending on the chosen cross
 correlation

---
 src/signalsnap/spectrum_plotter.py | 77 +++++++++++++++++++++++++++---
 1 file changed, 70 insertions(+), 7 deletions(-)

diff --git a/src/signalsnap/spectrum_plotter.py b/src/signalsnap/spectrum_plotter.py
index 92e282f..071d3ef 100644
--- a/src/signalsnap/spectrum_plotter.py
+++ b/src/signalsnap/spectrum_plotter.py
@@ -17,8 +17,9 @@
 
 
 class SpectrumPlotter:
-    def __init__(self, spectrum_calculator: SpectrumCalculator, plot_config: PlotConfig):
+    def __init__(self, spectrum_calculator: SpectrumCalculator, config: SpectrumConfig, plot_config: PlotConfig):
         self.spectrum_calculator = spectrum_calculator
+        self.config = config
         self.plot_config = plot_config
 
     def __import_spec_data_for_plotting(self, s_data, s_err, order):
@@ -172,13 +173,20 @@ def plot_s2(self, ax, order, s_f_plot, s_data_plot, s_err_plot):
                            linewidth=2, label=r"$%i\sigma$" % (i + 1), alpha=0.5)
 
         ax[0].plot(s_f_plot[order], s_data_plot[order], color=[0, 0.5, 0.9], linewidth=3)
-
         ax[0].tick_params(axis='both', direction='in')
-        ax[0].set_ylabel(r"$S^{(2)}_z$ (" + self.spectrum_calculator.config.f_unit + r"$^{-1}$)", labelpad=13,
-                         fontdict={'fontsize': 14})
-        ax[0].set_xlabel(r"$\omega / 2\pi$ (" + self.spectrum_calculator.config.f_unit + r")", labelpad=13,
-                         fontdict={'fontsize': 14})
-        ax[0].set_title(r"$S^{(2)}_z$ (" + self.spectrum_calculator.config.f_unit + r"$^{-1}$)", fontdict={'fontsize': 16})
+        if self.config.corr_data is not None:
+            ax[0].set_ylabel(r"$S^{(2)}_{1, 2}$ (" + self.spectrum_calculator.config.f_unit + r"$^{-1}$)", labelpad=13,
+                            fontdict={'fontsize': 14})
+            ax[0].set_xlabel(r"$\omega / 2\pi$ (" + self.spectrum_calculator.config.f_unit + r")", labelpad=13,
+                            fontdict={'fontsize': 14})
+            ax[0].set_title(r"$S^{(2)}_{1, 2}$ (" + self.spectrum_calculator.config.f_unit + r"$^{-1}$)", fontdict={'fontsize': 16})
+
+        else:
+            ax[0].set_ylabel(r"$S^{(2)}_z$ (" + self.spectrum_calculator.config.f_unit + r"$^{-1}$)", labelpad=13,
+                            fontdict={'fontsize': 14})
+            ax[0].set_xlabel(r"$\omega / 2\pi$ (" + self.spectrum_calculator.config.f_unit + r")", labelpad=13,
+                            fontdict={'fontsize': 14})
+            ax[0].set_title(r"$S^{(2)}_z$ (" + self.spectrum_calculator.config.f_unit + r"$^{-1}$)", fontdict={'fontsize': 16})
 
         if self.plot_config.broken_lims is not None:
             ylims = ax[0].get_ylim()
@@ -262,7 +270,61 @@ def plot_s3_s4(self, fig, ax, axis, order, s_data_plot, s_err_plot, s_f_plot):
         ax[axis].set_ylabel(r"$\omega_2 / 2 \pi$ (" + self.spectrum_calculator.config.f_unit + r")",
                             fontdict={'fontsize': 14})
         ax[axis].tick_params(axis='both', direction='in')
+        
+        if self.plot_config.green_alpha == 0 and order == 3:
+            if self.config.data3 is not None and self.config.corr_data is not None:
+                ax[axis].set_title(
+                    fr'$S^{{(3)}}_{{1, 2, 3}} $ (' + self.spectrum_calculator.config.f_unit + r'$^{-2}$)',
+                    fontdict={'fontsize': 16})
+            elif self.config.corr_data is not None and self.config.combination3 is not None:
+                ax[axis].set_title(
+                    fr'$S^{{(3)}}_{{{self.config.ALLOWED_COMBINATION3[self.config.combination3]}}} $ (' + self.spectrum_calculator.config.f_unit + r'$^{-2}$)',
+                    fontdict={'fontsize': 16})
+            else:
+                ax[axis].set_title(
+                    r'$S^{(3)}_z $ (' + self.spectrum_calculator.config.f_unit + r'$^{-2}$)',
+                    fontdict={'fontsize': 16})
+        elif self.plot_config.green_alpha != 0 and order == 3:
+            if self.config.data3 is not None and self.config.corr_data is not None:
+                ax[axis].set_title(
+                    fr'$S^{{(3)}}_{{1, 2, 3}} $ (' + self.spectrum_calculator.config.f_unit + r'$^{-2}$) (%i$\sigma$ confidence)' % (
+                    self.plot_config.sigma),
+                    fontdict={'fontsize': 16})
+            elif self.config.corr_data is not None:
+                ax[axis].set_title(
+                    fr'$S^{{(3)}}_{{{self.config.ALLOWED_COMBINATION3[self.config.combination3]}}} $ (' + self.spectrum_calculator.config.f_unit + r'$^{-2}$) (%i$\sigma$ confidence)' % (
+                    self.plot_config.sigma),
+                    fontdict={'fontsize': 16})
+            else:
+                ax[axis].set_title(
+                    r'$S^{(3)}_z $ (' + self.spectrum_calculator.config.f_unit + r'$^{-2}$) (%i$\sigma$ confidence)' % (
+                    self.plot_config.sigma),
+                    fontdict={'fontsize': 16})
+
 
+        if self.plot_config.green_alpha == 0 and order == 4:
+            if self.config.data3 is not None or self.config.corr_data is not None:
+                ax[axis].set_title(
+                    fr'$S^{{(4)}}_{{{self.config.ALLOWED_COMBINATION4[self.config.combination4]}}} $ (' + self.spectrum_calculator.config.f_unit + r'$^{-3}$)',
+                    fontdict={'fontsize': 16})
+            else:
+                ax[axis].set_title(
+                    r'$S^{(4)}_z $ (' + self.spectrum_calculator.config.f_unit + r'$^{-4}$)',
+                    fontdict={'fontsize': 16})
+        
+        
+        if self.plot_config.green_alpha != 0 and order == 4:
+            if self.config.data3 is not None or self.config.corr_data is not None:
+                ax[axis].set_title(
+                    fr'$S^{{(4)}}_{{{self.config.ALLOWED_COMBINATION4[self.config.combination4]}}} $ (' + self.spectrum_calculator.config.f_unit + r'$^{-3}$)(%i$\sigma$ confidence)' % (
+                    self.plot_config.sigma),
+                    fontdict={'fontsize': 16})
+            else:
+                ax[axis].set_title(
+                    r'$S^{(4)}_z $ (' + self.spectrum_calculator.config.f_unit + r'$^{-4}$) (%i$\sigma$ confidence)' % (
+                    self.plot_config.sigma),
+                    fontdict={'fontsize': 16})
+        """
         if self.plot_config.green_alpha == 0:
             ax[axis].set_title(
                 r'$S^{(' + f'{order}' + r')}_z $ (' + self.spectrum_calculator.config.f_unit + r'$^{-' + f'{order - 1}' + r'}$)',
@@ -272,6 +334,7 @@ def plot_s3_s4(self, fig, ax, axis, order, s_data_plot, s_err_plot, s_f_plot):
                 r'$S^{(' + f'{order}' + r')}_z $ (' + self.spectrum_calculator.config.f_unit + r'$^{-' + f'{order - 1}' + r'}$) (%i$\sigma$ confidence)' % (
                     self.plot_config.sigma),
                 fontdict={'fontsize': 16})
+        """
         fig.colorbar(c, ax=(ax[axis]))
 
         if self.plot_config.broken_lims is not None:

From ec4c3c83e3f16bf2e9aea62674d337bb223a59c5 Mon Sep 17 00:00:00 2001
From: Armin Etemad <130641768+ArminGEtemad@users.noreply.github.com>
Date: Mon, 1 Jul 2024 09:03:25 +0200
Subject: [PATCH 4/4] corrected a minor issue in ALLOWED_COMBINATION4 dict

---
 src/signalsnap/spectrum_config.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/signalsnap/spectrum_config.py b/src/signalsnap/spectrum_config.py
index d5b2ed3..c1b5e3d 100644
--- a/src/signalsnap/spectrum_config.py
+++ b/src/signalsnap/spectrum_config.py
@@ -102,7 +102,7 @@ class SpectrumConfig:
                             '3_2311': '2, 3, 1, 1', 
                             '3_1213': '1, 2, 1, 3',
                             '2_1212': '1, 2, 1, 2',
-                            '2_1122': '1, 2, 1, 2'}
+                            '2_1122': '1, 1, 2, 2'}
 
     def __init__(self, path=None, group_key=None, dataset=None, delta_t=None, data=None,
                  corr_data=None, corr_path=None, corr_group_key=None, corr_dataset=None,