Skip to content
Snippets Groups Projects
dualmicro_functions.py 88.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • glia's avatar
    glia committed
    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
    # -*- coding: utf-8 -*-
    """
    Functions for the analysis of two-brain microstates
    
    @author: Qianliang Li (glia@dtu.dk)
    """
    import os
    import numpy as np
    import pandas as pd
    import mne
    import mat73
    import matplotlib.pyplot as plt
    import nolds
    
    from eeg_microstates3 import (locmax, T_empirical, H_1, H_2)
    
    def load_epoch_from_fieldtrip(i, files, event_id):
        """Loads the preprocessed EEG data """
        # Get environmental variables
        Subject_id = os.environ.get('Subject_id')
        # Load the file
        tmp_dict = mat73.loadmat(files[i])
        s_id = Subject_id[i]
        eeg = tmp_dict["eeg_ppn"+str(s_id)[-1]+"_seg"]
        # Get the epoched time series data (epoch, ch, time)
        data = np.stack(eeg["trial"])
        # Re-scale data
        data *= 1e-6 # micro Volt to Volt
        # Create info instance
        ch_names = [ele[0] for ele in eeg["hdr"]["label"]]
        sfreq = int(eeg["fsample"])
        eeg_info = mne.create_info(ch_names, sfreq,ch_types="eeg")
        # Update filter info
        mne.filter._filt_update_info(eeg_info, True, 1.0, 40)
        # Make events array (epoch number, _, event id)
        events = np.zeros((len(data),3)).astype(int)
        events[:,0] = np.arange(0,len(data))
        events[:,2] = eeg["trialinfo"][:,0]
        # Save trialinfo object for synchronization of pairs based on timings
        trialinfo = eeg["trialinfo"][:,0:3]
        # Create the Epochs
        epoch = mne.EpochsArray(data, eeg_info, events, 0, event_id)
        # Set montage
        epoch.set_montage("biosemi64")
        return epoch, trialinfo
    
    # The function takes the data as a numpy array (n_t, n_ch)
    def prepare_1P_micro_arr(i, ppn2_correction, sfreq, freq_range=None, standardize:bool=True):
        """ Prepare single-brain EEG data for microstate analysis """
        # Get environmental variables
        Subject_id = os.environ.get('Subject_id')
        # Load epochs
        epoch, trialinfo = load_epoch_from_fieldtrip(i)
        # Ensure it is averaged referenced
        epoch = epoch.set_eeg_reference("average")
        # Get numpy arrays
        micro_data = epoch.get_data() # get data
        if standardize == True:
            # Standardize data - will make EEG amplitudes comparable between
            micro_data = micro_data - micro_data.mean(axis=1, keepdims=True)
            micro_data /= micro_data.std(axis=1, keepdims=True)
        # Correct the event_id for ppn2 in each pair, to swap all 6 with 4
        # and 7 with 5 and vice versa
        if str(Subject_id[i])[-1]=="2":
            trialinfo_df = pd.DataFrame(trialinfo.copy())
            trialinfo_df.iloc[:,0].replace(ppn2_correction, inplace=True)
            trialinfo = trialinfo_df.to_numpy()
    
        # Transform data to correct shape
        arr_shape = micro_data.shape # get shape
        micro_data = micro_data.swapaxes(1,2) # swap ch and time axis
        micro_data = micro_data.reshape(arr_shape[0]*arr_shape[2],arr_shape[1]) # reshape by combining epochs and times
        
        # Filter the data
        if freq_range == None:
            micro_data_filtered = micro_data # do not perform filtering
        elif not freq_range == None:
            # Notice it is done only after combining epochs and time, as filter length
            # would be too long for 1s epochs. The filter function wants time on the last axis
            micro_data = micro_data.transpose()
            micro_data_filtered = mne.filter.filter_data(micro_data, sfreq, freq_range[0], freq_range[1])
            # Reverse the shape
            micro_data_filtered = micro_data_filtered.transpose()
        
        return micro_data_filtered, trialinfo
    
    def plot_microstates(n_maps, maps, gev, epoch_info):
        """ Plot the microstates """
        fig, axarr = plt.subplots(1, n_maps, figsize=(5*n_maps,5))
        fig.patch.set_facecolor('white')
        for imap in range(n_maps):
            mne.viz.plot_topomap(maps[imap,:], pos = epoch_info, axes = axarr[imap]) # plot
            axarr[imap].set_title("GEV: {:.2f}".format(gev[imap]), fontsize=16, fontweight="bold") # title
        fig.suptitle("Microstates: {}".format(n_maps), fontsize=20, fontweight="bold")
        plt.tight_layout()
        return fig
    
    def reorder_microstate_results(new_order, maps, gev, m_labels):
        """ Function to re-order microstates based on manual order """
        reordered_maps = maps[new_order,:]
        reordered_gev = gev[new_order]
        # Make directory to find and replace map labels
        dic0 = {value:key for key, value in enumerate(new_order)}
        reordered_m_labels = np.array([dic0.get(n, n) for n in m_labels]) # re-order labels
        return reordered_maps, reordered_gev, reordered_m_labels
    
    def microstate_run_length_encoding(m_labels):
        """ Take a 1D stream of microstates and returns the label and duration """
        # Adapted from RLE Matlab code by Abdulrahman Ikram Siddiq, 2011
        # Initialize
        counter = 0
        label = [9999] * len(m_labels)
        duration = [9999] * len(m_labels)
        # Starting the lists
        label[counter] = m_labels[counter] # first element
        duration[counter] = 1 # first element
        for i in range(1,len(m_labels)):
            # Add 1 to duration, if the next element is the same as current microstate
            if m_labels[i-1] == m_labels[i]:
                duration[counter] += 1
            # If it is not the same, get the new microstate label and restart duration
            else:
                counter += 1
                label[counter] = m_labels[i]
                duration[counter] = 1
        # Trim the unused parts of the list
        label = label[:counter+1]
        duration = duration[:counter+1]
        assert not any(np.array(label) == 9999) # all fillers removed
        assert not any(np.array(duration) == 9999) # all fillers removed
        
        return np.array(label), np.array(duration)
    
    def single_micro_fit_all_feature_computation(i, n_maps, microstate_results, trialinfo_list, sfreq, event_id,
                                                 single_brain_event_id):
        """
        Estimates the common (intrabrain) microstate features:
            1. Average duration a given microstate remains stable (Dur)
            2. Frequency occurrence, independent of individual duration (Occ)
                Average number of times a microstate becomes dominant per second
            3. Ratio of total Time Covered (TCo)
            4. Transition probabilities (TMx)
            5. Ratio of shannon entropy relative to theoretical max chaos (Ent)
        
        Parameters
        ----------
        i : int
            The index.
        n_maps : int
            The number of maps (clusters) used.
        microstate_results : list
            The estimated microstates.
        trialinfo_list : list
            List with trial informations.
        sfreq : int
            The sampling frequency.
        event_id : dict
            The mappings between event id and condition.
        single_brain_event_id : dict
            The renamed mappings between event id and condition for single-brain
            analysis.
    
        Returns
        -------
        m_labels : np.array
            The microstate sequence (labels) time series in the format (epoch, time).
        events : pd.DataFrame
            The events corresponding to each epoch.
        MFeatures : list
            List of arrays of each microstate feature.
    
        """
        # Get environmental variables
        Subject_id = os.environ.get('Subject_id')
        # Get the microstate labels
        sub_idx = microstate_results[5]
        subject_indices = sub_idx[i], sub_idx[i+1]
        m_labels = microstate_results[1][subject_indices[0]:subject_indices[1]]
        
        # Get the trialinfo with conditions
        # Notice that ppn2 events have been corrected, so we are using separate
        # observer and actor conditions
        # and also separate follower and leader conditions
        Subject0, trialinfo = trialinfo_list[i]
        assert Subject_id[i] == Subject0
        assert len(trialinfo) == len(m_labels)/sfreq
        # Convert to dataframe
        event_id_inv = {v: k for k, v in event_id.items()} # Inverse the event id
        events = pd.DataFrame(trialinfo,columns=["Event_id","Trial_number","Trial_start_time"])
        events["Event_label"] = events["Event_id"].replace(event_id_inv)
        events = events.reset_index().rename(columns={"index":"Epoch_idx"})
        
        # Reshape m_labels to (epoch, time)
        m_labels = m_labels.reshape(len(trialinfo),sfreq)
        
        # Remove pre-trial epochs
        pre_trial_epochs = events["Trial_start_time"] < 0
        m_labels = m_labels[np.invert(pre_trial_epochs)]
        events = events.loc[np.invert(pre_trial_epochs)].reset_index(drop=True)
        events["Epoch_idx"] = events.index
        
        # Initialize arrays
        Dur_arr = np.zeros((len(single_brain_event_id),n_maps)); Dur_arr.fill(np.nan)
        Occ_arr = np.zeros((len(single_brain_event_id),n_maps)); Occ_arr.fill(np.nan)
        TCo_arr = np.zeros((len(single_brain_event_id),n_maps)); TCo_arr.fill(np.nan)
        TMx_arr = np.zeros((len(single_brain_event_id),n_maps,n_maps))
        Ent_arr = np.zeros(len(single_brain_event_id))
        for e in range(len(single_brain_event_id)):
            ev_idx = list(single_brain_event_id.values())[e]
            ep_idx = events["Epoch_idx"][events["Event_id"] == ev_idx]
            trial_numbers0 = np.unique(events["Trial_number"][ep_idx])
            if len(trial_numbers0) > 8:
                # Compute the first 8 trials and the last 8 separately, before
                # taking the average to be consistent with asymmetric trials
                Dur_arr0 = np.zeros((2,n_maps))
                Occ_arr0 = np.zeros((2,n_maps))
                TCo_arr0 = np.zeros((2,n_maps))
                TMx_arr0 = np.zeros((2,n_maps,n_maps))
                Ent_arr0 = np.zeros(2)
                trial_numbers_split = np.array_split(trial_numbers0, 2)
                # Array_split is used in cases where a trial might have been
                # dropped, hence the split is only 8/7
                for s in range(len(trial_numbers_split)):
                    ep_idx0 = events.loc[(events["Event_id"] == ev_idx)&
                                        (events["Trial_number"].isin(trial_numbers_split[s])),
                                        "Epoch_idx"]
                    m_labels0 = m_labels[ep_idx0,:]
                    m_labels0_flat = m_labels0.reshape(m_labels0.shape[0]*m_labels0.shape[1])
                    # Compute duration, occurrence and time covered
                    l_, d_ = microstate_run_length_encoding(m_labels0_flat)
                    # For each microstate
                    for ii in range(n_maps):
                        if np.isnan(np.nanmean(d_[l_==ii])):
                            # The specific microstate did not occur at all for this event
                            Dur_arr0[s,ii] = 0
                            Occ_arr0[s,ii] = 0
                            TCo_arr0[s,ii] = 0
                        else:
                            Dur_arr0[s,ii] = np.mean(d_[l_==ii]) * 1000/sfreq # convert to ms
                            Occ_arr0[s,ii] = len(d_[l_==ii])/len(d_) * sfreq
                            TCo_arr0[s,ii] = np.sum(d_[l_==ii])/np.sum(d_)
                    # Compute transition matrix
                    TMx_arr0[s] = T_empirical(m_labels0_flat, n_maps)
                    # Compute Shannon Entropy relative to max
                    Ent_arr0[s] = H_1(m_labels0_flat, n_maps)/np.log(float(n_maps))
                # Average over the splits
                Dur_arr[e] = np.mean(Dur_arr0, axis=0)
                Occ_arr[e] = np.mean(Occ_arr0, axis=0)
                TCo_arr[e] = np.mean(TCo_arr0, axis=0)
                TMx_arr[e] = np.mean(TMx_arr0, axis=0)
                Ent_arr[e] = np.mean(Ent_arr0, axis=0)
                
            else:
                m_labels0 = m_labels[ep_idx,:]
                m_labels0_flat = m_labels0.reshape(m_labels0.shape[0]*m_labels0.shape[1])
                # Compute duration, occurrence and time covered
                l_, d_ = microstate_run_length_encoding(m_labels0_flat)
                # For each microstate
                for ii in range(n_maps):
                    if np.isnan(np.nanmean(d_[l_==ii])):
                        # The specific microstate did not occur at all for this event
                        Dur_arr[e,ii] = 0
                        Occ_arr[e,ii] = 0
                        TCo_arr[e,ii] = 0
                    else:
                        Dur_arr[e,ii] = np.mean(d_[l_==ii]) * 1000/sfreq # convert to ms
                        Occ_arr[e,ii] = len(d_[l_==ii])/len(d_) * sfreq
                        TCo_arr[e,ii] = np.sum(d_[l_==ii])/np.sum(d_)
                # Compute transition matrix
                TMx_arr[e] = T_empirical(m_labels0_flat, n_maps)
                # Compute Shannon Entropy relative to max
                Ent_arr[e] = H_1(m_labels0_flat, n_maps)/np.log(float(n_maps))
        # Save all microstate features together
        MFeatures = [Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr]
        return m_labels, events, MFeatures
    
    def interbrain_microstate_run_length_encoding(m_labels1, m_labels2):
        """
        Takes two 1D streams of microstate labels and returns the label
        and duration for the common microstate. If there are no common
        microstate, the label -1 is returned with a corresponding duration.
        
        Adapted from RLE Matlab code by Abdulrahman Ikram Siddiq, 2011
        """
        assert len(m_labels1) == len(m_labels2), "The provided microstate labels do not have equal length"
        # Initialize
        counter = 0
        label = [9999] * len(m_labels1)
        duration = [9999] * len(m_labels1)
        # Starting the lists
        if m_labels1[counter] == m_labels2[counter]:
            label[counter] = m_labels1[counter]
        else:
            label[counter] = -1
        duration[counter] = 1 # first element
        for i in range(1,len(m_labels1)):
            # If the previous timepoint was not common microstate, then add
            # 1 to the duration if the next timepoint is also not common
            if label[counter] == -1:
                if m_labels1[i] != m_labels2[i]:
                    duration[counter] += 1
                # If it becomes a common microstate, then get new microstate
                # label and restart duration
                elif m_labels1[i] == m_labels2[i]:
                    counter += 1
                    label[counter] = m_labels1[i]
                    duration[counter] = 1
            # If the previous time point was a common microstate, then add
            # 1 to duration if the next time point is also common and the same label
            elif label[counter] != -1:
                if m_labels1[i] == m_labels2[i]:
                    if label[counter] == m_labels1[i]:
                        duration[counter] += 1
                    # If still common, but a new label for both timeseries
                    # then restart with new label and duration
                    else:
                        counter += 1
                        label[counter] = m_labels1[i]
                        duration[counter] = 1
                # If the two timeseries have different microstates, then
                # set lable to -1 and restart duration
                elif m_labels1[i] != m_labels2[i]:
                    counter += 1
                    label[counter] = -1
                    duration[counter] = 1
        
        # Trim the unused parts of the list
        label = label[:counter+1]
        duration = duration[:counter+1]
        assert not any(np.array(label) == 9999) # all fillers removed
        assert not any(np.array(duration) == 9999) # all fillers removed
        
        return np.array(label), np.array(duration)
    
    def interbrain_T_matrix(m_labels1, m_labels2, n_maps):
        """
        Takes two 1D microstate labels and returns the transition matrix
        The first row and column correspond to label -1, which indicates
        no common microstate
        
        Adapted from von Wegner F and Laufs H, 2018
        """
        assert len(m_labels1) == len(m_labels2), "The provided microstate labels do not have equal length"
        T = np.zeros((n_maps+1, n_maps+1))
        n = len(m_labels1)
        for i in range(n-1):
            # If common microstate, then use that label
            if m_labels1[i] == m_labels2[i]:
                row_idx = m_labels1[i]+1
            # Or set row as 0, the not common microstate
            else:
                row_idx = 0
            # If next state is also common, set it to that state
            if m_labels1[i+1] == m_labels2[i+1]:
                col_idx = m_labels1[i+1]+1
            # Or set col as 0, if it is not common microstate
            else:
                col_idx = 0
            # Add 1 count to the transition matrix
            T[row_idx,col_idx] += 1.0
        assert n - np.sum(T) == 1 # check whether all transitions have been counted
        # Normalize row sums to 1
        T /= T.sum(axis=1, keepdims=True)
        
        return T
    
    def interbrain_microstate_feature_computation(i, n_maps, microstate_results, trialinfo_list, sfreq, event_id,
                                                  collapsed_event_id):
        """
        Interbrain features:
            1. Average duration of common interbrain microstates (IBDur)
            2. Frequency occurrence of common interbrain microstates in the pair (IBOcc)
            3. Ratio of total time covered by interbrain common microstates in the pair (IBCov)
            4. Transition probability towards common interbrain microstates in the pair (IBTMx)
            5. Ratio of joint shannon entropy relative to theoretical max chaos (IBEnt)
        
        Parameters
        ----------
        i : int
            The index.
        n_maps : int
            The number of maps (clusters) used.
        microstate_results : list
            The estimated microstates.
        trialinfo_list : list
            List with trial informations.
        sfreq : int
            The sampling frequency.
        event_id : dict
            The mappings between event id and condition.
        collapsed_brain_event_id : dict
            The renamed mappings between event id and collapsed conditions
    
        Returns
        -------
        m_labels : list of two np.array
            The microstate sequence (labels) time series in the format (epoch, time).
        events : list of two pd.DataFrame
            The events corresponding to each epoch.
        MFeatures : list
            List of arrays of each microstate feature.
            
        """
        # Here i refers to the pair in range(n_subjects//2)
        sub_idx = microstate_results[5]
        
        # Get the microstate labels and events for participant 1
        subject_indices1 = sub_idx[2*i], sub_idx[2*i+1]
        m_labels1 = microstate_results[1][subject_indices1[0]:subject_indices1[1]]
        # Get the trialinfo with conditions
        Subject1, trialinfo1 = trialinfo_list[2*i]
        # Convert to dataframe
        event_id_inv = {v: k for k, v in event_id.items()} # Inverse the event id
        
        events1 = pd.DataFrame(trialinfo1,columns=["Event_id","Trial_number","Trial_start_time"])
        events1["Event_label"] = events1["Event_id"].replace(event_id_inv)
        events1 = events1.reset_index().rename(columns={"index":"Epoch_idx"})
        # Reshape m_labels to (epoch, time)
        m_labels1 = m_labels1.reshape(len(trialinfo1),sfreq)
        
        # Get the microstate labels and events for participant 2
        subject_indices2 = sub_idx[2*i+1], sub_idx[2*i+2]
        m_labels2 = microstate_results[1][subject_indices2[0]:subject_indices2[1]]
        # Get the trialinfo with conditions
        Subject2, trialinfo2 = trialinfo_list[2*i+1]
        # Convert to dataframe
        events2 = pd.DataFrame(trialinfo2,columns=["Event_id","Trial_number","Trial_start_time"])
        events2["Event_label"] = events2["Event_id"].replace(event_id_inv)
        events2 = events2.reset_index().rename(columns={"index":"Epoch_idx"})
        # Reshape m_labels to (epoch, time)
        m_labels2 = m_labels2.reshape(len(trialinfo2),sfreq)
        
        # check that the participants from a pair was loaded
        assert Subject2-1 == Subject1
        # Check that the same amount of event types are present
        assert trialinfo1[-1,1] == trialinfo2[-1,1]
        # Synchronize the events from the pair based on timing info
        # By trimming the epochs to only include the epochs that are present
        # in both participants of the pair
        # Initialize with an array filled with a unique number not in use
        sync_m_labels1 = np.zeros_like(m_labels1); sync_m_labels1.fill(9999)
        sync_m_labels2 = np.zeros_like(m_labels2); sync_m_labels2.fill(9999)
        for t in np.unique(trialinfo1[:,1]):
            t_idx1 = np.where(trialinfo1[:,1] == t)[0]
            t_idx2 = np.where(trialinfo2[:,1] == t)[0]
            # Get the timings for the epochs for the specific trial t
            t_timings1 = trialinfo1[t_idx1,2]
            t_timings2 = trialinfo2[t_idx2,2]
            timings_intersect = np.intersect1d(t_timings1,t_timings2)
            # Get the indices where the timings matches
            t_idx_match1 = t_idx1[pd.Series(t_timings1).isin(timings_intersect)]
            t_idx_match2 = t_idx2[pd.Series(t_timings2).isin(timings_intersect)]
            # Get the actual values from the synchronized epochs
            sync_m_labels1[t_idx_match1] = m_labels1[t_idx_match1]
            sync_m_labels2[t_idx_match2] = m_labels2[t_idx_match2]
        # Find the epochs that were asynchronous, which have to be trimmed
        asynch_epochs1 = np.unique(np.where(sync_m_labels1==9999)[0])
        asynch_epochs2 = np.unique(np.where(sync_m_labels2==9999)[0])
        # Trim/delete the asynchronous epochs
        sync_m_labels1 = np.delete(sync_m_labels1,asynch_epochs1,axis=0)
        sync_m_labels2 = np.delete(sync_m_labels2,asynch_epochs2,axis=0)
        assert len(sync_m_labels1) == len(sync_m_labels2) # check the amount of synchronized epochs are equal
        # Fix events
        sync_events1 = events1.drop(asynch_epochs1,axis=0).reset_index(drop=True)
        sync_events2 = events2.drop(asynch_epochs2,axis=0).reset_index(drop=True)
        
        # Since they are synchronized already, just take the first and fix
        # epoch idx column with the new idx for the synchronized labels
        events = sync_events1
        events["Epoch_idx"] = events.index
        # Notice that for the intrabrain fit all alpha v6 I only corrected
        # ppn2. So by taking ppn1 I get the original event labels.
        
        # Update in v7
        # # Collapse event_id 6 and 7 to 4 and 5
        # events.loc[events["Event_id"] == 6,["Event_id","Event_label"]] = [4, "observe, actor"]
        # events.loc[events["Event_id"] == 7,["Event_id","Event_label"]] = [5, "imitate, leader"]
    
        # Due to the hard constraint for both participants being in the same
        # microstate. It does not matter who will be treated as ppn1 and ppn2
        # since the prototypical microstate topographies are exactly the same!
        
        # Remove pre-trial epochs
        pre_trial_epochs = events["Trial_start_time"] < 0
        sync_m_labels1 = sync_m_labels1[np.invert(pre_trial_epochs)]
        sync_m_labels2 = sync_m_labels2[np.invert(pre_trial_epochs)]
        events = events.loc[np.invert(pre_trial_epochs)].reset_index(drop=True)
        events["Epoch_idx"] = events.index
        
        # Pre-allocate memory
        Dur_arr = np.zeros((len(event_id),n_maps+1)); Dur_arr.fill(np.nan)
        Occ_arr = np.zeros((len(event_id),n_maps+1)); Occ_arr.fill(np.nan)
        TCo_arr = np.zeros((len(event_id),n_maps+1)); TCo_arr.fill(np.nan)
        TMx_arr = np.zeros((len(event_id),n_maps+1,n_maps+1))
        Ent_arr = np.zeros(len(event_id))
        for e in range(len(event_id)):
            ev_idx = list(event_id.values())[e]
            ep_idx = events["Epoch_idx"][events["Event_id"] == ev_idx]
            trial_numbers0 = np.unique(events["Trial_number"][ep_idx])
            if len(trial_numbers0) > 8:
                # Compute the first 8 trials and the last 8 separately, before
                # taking the average to be consistent with asymmetric trials
                Dur_arr0 = np.zeros((2,n_maps+1))
                Occ_arr0 = np.zeros((2,n_maps+1))
                TCo_arr0 = np.zeros((2,n_maps+1))
                TMx_arr0 = np.zeros((2,n_maps+1,n_maps+1))
                Ent_arr0 = np.zeros(2)
                trial_numbers_split = np.array_split(trial_numbers0, 2)
                # Array_split is used in cases where a trial might have been
                # dropped, hence the split is only 8/7
                for s in range(len(trial_numbers_split)):
                    ep_idx0 = events.loc[(events["Event_id"] == ev_idx)&
                                        (events["Trial_number"].isin(trial_numbers_split[s])),
                                        "Epoch_idx"]
                    # Get the microstate labels
                    m_labels10 = sync_m_labels1[ep_idx0,:]
                    m_labels20 = sync_m_labels2[ep_idx0,:]
                    # Flatten the labels
                    m_labels10_flat = m_labels10.reshape(m_labels10.shape[0]*m_labels10.shape[1])
                    m_labels20_flat = m_labels20.reshape(m_labels20.shape[0]*m_labels20.shape[1])
                    # Compute average duration of common microstate
                    # Output: label and duration of common microstate. Label -1 is used
                    # for not common microstate
                    l_, d_ = interbrain_microstate_run_length_encoding(m_labels10_flat,m_labels20_flat)
                    # For each microstate
                    for ii in range(n_maps+1):
                        if np.isnan(np.nanmean(d_[l_==ii-1])):
                            # The specific microstate did not occur at all for this event
                            Dur_arr0[s,ii] = 0
                            Occ_arr0[s,ii] = 0
                            TCo_arr0[s,ii] = 0
                        else:
                            Dur_arr0[s,ii] = np.mean(d_[l_==ii-1]) * 1000/sfreq # convert to ms
                            Occ_arr0[s,ii] = len(d_[l_==ii-1])/len(d_) * sfreq
                            TCo_arr0[s,ii] = np.sum(d_[l_==ii-1])/np.sum(d_)
                    # Compute transition matrix
                    TMx_arr0[s] = interbrain_T_matrix(m_labels10_flat,m_labels20_flat,n_maps)
                    # Compute Joint Shannon Entropy relative to max
                    # Max is the sum of max individual entropies
                    Ent_arr0[s] = H_2(m_labels10_flat,m_labels20_flat,n_maps)/(2*np.log(float(n_maps)))
                # Average over the splits
                Dur_arr[e] = np.mean(Dur_arr0, axis=0)
                Occ_arr[e] = np.mean(Occ_arr0, axis=0)
                TCo_arr[e] = np.mean(TCo_arr0, axis=0)
                TMx_arr[e] = np.mean(TMx_arr0, axis=0)
                Ent_arr[e] = np.mean(Ent_arr0, axis=0)
                
            else:
                # Get the microstate labels
                m_labels10 = sync_m_labels1[ep_idx,:]
                m_labels20 = sync_m_labels2[ep_idx,:]
                # Flatten the labels
                m_labels10_flat = m_labels10.reshape(m_labels10.shape[0]*m_labels10.shape[1])
                m_labels20_flat = m_labels20.reshape(m_labels20.shape[0]*m_labels20.shape[1])
                # Compute average duration of common microstate
                # Output: label and duration of common microstate. Label -1 is used
                # for not common microstate
                l_, d_ = interbrain_microstate_run_length_encoding(m_labels10_flat,m_labels20_flat)
                # For each microstate
                for ii in range(n_maps+1):
                    if np.isnan(np.nanmean(d_[l_==ii-1])):
                        # The specific microstate did not occur at all for this event
                        Dur_arr[e,ii] = 0
                        Occ_arr[e,ii] = 0
                        TCo_arr[e,ii] = 0
                    else:
                        Dur_arr[e,ii] = np.mean(d_[l_==ii-1]) * 1000/sfreq # convert to ms
                        Occ_arr[e,ii] = len(d_[l_==ii-1])/len(d_) * sfreq
                        TCo_arr[e,ii] = np.sum(d_[l_==ii-1])/np.sum(d_)
                # Compute transition matrix
                TMx_arr[e] = interbrain_T_matrix(m_labels10_flat,m_labels20_flat,n_maps)
                # Compute Joint Shannon Entropy relative to max
                # Max is the sum of max individual entropies
                Ent_arr[e] = H_2(m_labels10_flat,m_labels20_flat,n_maps)/(2*np.log(float(n_maps)))
                
        # Combine all features in a list
        MFeatures = [Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr]
        
        # Update in v7 - Collapse after computations in single events
        Dur_arr2 = np.zeros((len(collapsed_event_id),n_maps+1))
        Occ_arr2 = np.zeros((len(collapsed_event_id),n_maps+1))
        TCo_arr2 = np.zeros((len(collapsed_event_id),n_maps+1))
        TMx_arr2 = np.zeros((len(collapsed_event_id),n_maps+1,n_maps+1))
        Ent_arr2 = np.zeros(len(collapsed_event_id))
        MFeatures2 = [Dur_arr2,Occ_arr2,TCo_arr2,TMx_arr2,Ent_arr2]
        
        for f in range(len(MFeatures2)):
            tmp_feat = MFeatures[f]
            tmp_feat2 = MFeatures2[f]
            for e in range(len(collapsed_event_id)):
                ee = list(collapsed_event_id.keys())[e]
                if (ee == 'observer_actor'):
                    old_ev_idx1 = list(event_id.keys()).index("observe, actor")
                    old_ev_idx2 = list(event_id.keys()).index("observe, observer")
                    tmp_feat2[e] = np.mean([tmp_feat[old_ev_idx1],tmp_feat[old_ev_idx2]], axis=0)
                elif ee == 'follower_leader':
                    old_ev_idx1 = list(event_id.keys()).index("imitate, leader")
                    old_ev_idx2 = list(event_id.keys()).index("imitate, follower")
                    tmp_feat2[e] = np.mean([tmp_feat[old_ev_idx1],tmp_feat[old_ev_idx2]], axis=0)
                else:
                    new_ev_idx = list(collapsed_event_id.values())[e]
                    old_ev_idx = list(event_id.values()).index(new_ev_idx)
                    tmp_feat2[e] = tmp_feat[old_ev_idx]
        
        return [sync_m_labels1, sync_m_labels2], [sync_events1, sync_events2], MFeatures2
    
    
    def get_synch_events_from_pairs(i, trialinfo1, trialinfo2, event_id):
        """ Get the synchronized events for each pair """
        # Get environmental variables
        Subject_id = os.environ.get('Subject_id')
        # Check that the participants from a pair was loaded
        assert Subject_id[2*i+1]-1 == Subject_id[2*i]
        # Synchronize the events from the pair based on timing info
        # By trimming the epochs to only include the epochs that are present
        # in both participants of the pair
        event_id_inv = {v: k for k, v in event_id.items()} # Inverse the event id
        events1 = pd.DataFrame(trialinfo1,columns=["Event_id","Trial_number","Trial_start_time"])
        events1["Event_label"] = events1["Event_id"].replace(event_id_inv)
        events1 = events1.reset_index().rename(columns={"index":"Epoch_idx"})
        
        events2 = pd.DataFrame(trialinfo2,columns=["Event_id","Trial_number","Trial_start_time"])
        events2["Event_label"] = events2["Event_id"].replace(event_id_inv)
        events2 = events2.reset_index().rename(columns={"index":"Epoch_idx"})
        # Concatenate all the synched events
        sync_events1 = np.zeros((len(events1)))
        sync_events2 = np.zeros((len(events2)))
        for t in np.unique(trialinfo1[:,1]):
            t_idx1 = np.where(trialinfo1[:,1] == t)[0]
            t_idx2 = np.where(trialinfo2[:,1] == t)[0]
            # Get the timings for the epochs for the specific trial t
            t_timings1 = trialinfo1[t_idx1,2]
            t_timings2 = trialinfo2[t_idx2,2]
            timings_intersect = np.intersect1d(t_timings1,t_timings2)
            # Get the indices where the timings matches
            t_idx_match1 = t_idx1[pd.Series(t_timings1).isin(timings_intersect)]
            t_idx_match2 = t_idx2[pd.Series(t_timings2).isin(timings_intersect)]
            # Save the trials that should be kept
            sync_events1[t_idx_match1] = 1
            sync_events2[t_idx_match2] = 1
        # Find the epochs that were asynchronous, which have to be trimmed
        asynch_epochs1 = np.where(sync_events1==0)[0]
        asynch_epochs2 = np.where(sync_events2==0)[0]
        # Trim/delete the asynchronous epochs
        events1 = events1.drop(asynch_epochs1,axis=0).reset_index(drop=True)
        events2 = events2.drop(asynch_epochs2,axis=0).reset_index(drop=True)
        assert len(events1) == len(events2) # check the amount of synchronized epochs are equal
        return events1, events2
    
    def prepare_2P_micro_arr_collapsed_events(i, sfreq, event_id, freq_range=None, standardize:bool=True):
        # Here i refers to the pair in range(n_subjects//2)
        # Load epochs
        epoch1, trialinfo1 = load_epoch_from_fieldtrip(2*i)
        epoch2, trialinfo2 = load_epoch_from_fieldtrip(2*i+1)
        # Ensure it is averaged referenced
        epoch1 = epoch1.set_eeg_reference("average")
        epoch2 = epoch2.set_eeg_reference("average")
        # Get numpy arrays
        micro_data1 = epoch1.get_data() # get data
        micro_data2 = epoch2.get_data() # get data
        # Get only the synchronized epochs
        events = get_synch_events_from_pairs(i, trialinfo1, trialinfo2, event_id)
        micro_data1 = micro_data1[events[0]["Epoch_idx"]]
        micro_data2 = micro_data2[events[1]["Epoch_idx"]]
        assert micro_data1.shape == micro_data2.shape
        # Get event labels for each epoch
        # Since they are already synchronized I just take the first
        event_labels0 = events[0]["Event_label"]
        if standardize == True:
            # Standardize data - will make EEG amplitudes comparable between
            # recordings/trials. Standardized along the ch axis, meaning the mean
            # and std of the topomap at each timepoint is normalized
            # It will be normalized separately for each participant
            # in order to ensure both participants are treated equally by the
            # kmeans
            micro_data10 = micro_data1 - micro_data1.mean(axis=1, keepdims=True)
            micro_data10 /= micro_data10.std(axis=1, keepdims=True)
            micro_data20 = micro_data2 - micro_data2.mean(axis=1, keepdims=True)
            micro_data20 /= micro_data20.std(axis=1, keepdims=True)
        elif standardize == False:
            micro_data10 = micro_data1
            micro_data20 = micro_data2
        # Concatenate along the channel axis, to treat the paired EEG as one entity
        micro_data = np.concatenate([micro_data10,micro_data20],axis=1)
    
        # Flip the micro_data between ppn1 and ppn2 to fix first microstate
        # always being observer/follower, and second microstate being actor/leader
        cond6_idx = event_labels0 == "observe, observer"
        cond7_idx = event_labels0 == "imitate, follower"
        cond67_idx = cond6_idx+cond7_idx
        micro_data[cond67_idx] = np.concatenate([micro_data20[cond67_idx],
                                                micro_data10[cond67_idx]],axis=1)
        # Transform data to correct shape
        arr_shape = micro_data.shape # get shape
        micro_data = micro_data.swapaxes(1,2) # swap ch and time axis
        micro_data = micro_data.reshape(arr_shape[0]*arr_shape[2],arr_shape[1]) # reshape by combining epochs and times
        
        # Filter the data
        if freq_range == None:
            micro_data_filtered = micro_data # do not perform filtering
        elif not freq_range == None:
            # Notice it is done only after combining epochs and time, as filter length
            # would be too long for 1s epochs. The filter function wants time on the last axis
            micro_data = micro_data.transpose()
            micro_data_filtered = mne.filter.filter_data(micro_data, sfreq, freq_range[0], freq_range[1])
            # Reverse the shape
            micro_data_filtered = micro_data_filtered.transpose()
        
        return micro_data_filtered, [trialinfo1, trialinfo2], events
        
    def plot_dualmicro(n_maps, maps, gev, epoch_info, vlims=(None, None)):
        # Split the map from the two participants
        n_channels = epoch_info["nchan"]
        maps1, maps2, _ = np.split(maps, [n_channels,2*n_channels], axis=1)
        fig, axarr = plt.subplots(2, n_maps, figsize=(20,10))
        fig.patch.set_facecolor('white')
        for imap in range(n_maps):
            mne.viz.plot_topomap(maps1[imap,:], pos = epoch_info, vlim = vlims, axes = axarr[0,imap]) # plot
            axarr[0,imap].set_title("GEV: {:.2f}".format(gev[imap]), fontsize=16, fontweight="bold") # title
            mne.viz.plot_topomap(maps2[imap,:], pos = epoch_info, vlim = vlims, axes = axarr[1,imap])
        fig.suptitle("Microstates: {}".format(n_maps), fontsize=20, fontweight="bold")
        plt.tight_layout()
        return fig
    
    def sign_swap_microstates(sign_swap, maps, n_maps, n_channels):
        # Split maps
        maps1, maps2, _ = np.split(maps, [n_channels,2*n_channels], axis=1)
        for m in range(n_maps):
            # Sign swap each
            maps1[m] *= sign_swap[0][m]
            maps2[m] *= sign_swap[1][m]
        # Concatenate maps
        maps = np.concatenate([maps1,maps2],axis=1)
        return maps
    
    def dualmicro_fit_all_feature_computation(i, n_maps, microstate_results, trialinfo_list, sfreq, event_id,
                                                  collapsed_event_id):
        """
        Overview of common microstate features:
            1. Average duration a given microstate remains stable (Dur)
            2. Frequency occurrence, independent of individual duration (Occ)
                Average number of times a microstate becomes dominant per second
            3. Ratio of total Time Covered (TCo)
            4. Transition probabilities (TMx)
            5. Ratio of shannon entropy relative to theoretical max chaos (Ent)
        
        Parameters
        ----------
        i : int
            The index.
        n_maps : int
            The number of maps (clusters) used.
        microstate_results : list
            The estimated microstates.
        trialinfo_list : list
            List with trial informations.
        sfreq : int
            The sampling frequency.
        event_id : dict
            The mappings between event id and condition.
        collapsed_brain_event_id : dict
            The renamed mappings between event id and collapsed conditions
    
        Returns
        -------
        m_labels : list of two np.array
            The microstate sequence (labels) time series in the format (epoch, time).
        events : list of two pd.DataFrame
            The events corresponding to each epoch.
        MFeatures : list
            List of arrays of each microstate feature.
    
        """
        
        pair_idx = microstate_results[5]
        pair_indices = pair_idx[i], pair_idx[i+1]
        m_labels = microstate_results[1][pair_indices[0]:pair_indices[1]]
        
        events = trialinfo_list[2][i]
        # Since they are synchronized already, just take the first and fix
        # epoch idx column with the new idx for the synchronized labels
        events = events[0]
        
        # Update in v7
        # # Collapse event_id 6 and 7 to 4 and 5
        # events.loc[events["Event_id"] == 6,["Event_id","Event_label"]] = [4, "observe, actor"]
        # events.loc[events["Event_id"] == 7,["Event_id","Event_label"]] = [5, "imitate, leader"]
        # The microstate clustering was performed on flipped (collapsed) events,
        # but I will compute the features on the 8 trials to avoid the flipping
        # effect and then collapse by averaging afterwards
        
        # Reshape m_labels to (epoch, time)
        m_labels = m_labels.reshape(len(events),sfreq)
        
        # Remove pre-trial epochs
        pre_trial_epochs = events["Trial_start_time"] < 0
        m_labels = m_labels[np.invert(pre_trial_epochs)]
        events = events.loc[np.invert(pre_trial_epochs)].reset_index(drop=True)
        events["Epoch_idx"] = events.index
        
        # Pre-allocate memory
        Dur_arr = np.zeros((len(event_id),n_maps)); Dur_arr.fill(np.nan)
        Occ_arr = np.zeros((len(event_id),n_maps)); Occ_arr.fill(np.nan)
        TCo_arr = np.zeros((len(event_id),n_maps)); TCo_arr.fill(np.nan)
        TMx_arr = np.zeros((len(event_id),n_maps,n_maps))
        Ent_arr = np.zeros(len(event_id))
        for e in range(len(event_id)):
            ev_idx = list(event_id.values())[e]
            ep_idx = events["Epoch_idx"][events["Event_id"] == ev_idx]
            trial_numbers0 = np.unique(events["Trial_number"][ep_idx])
            if len(trial_numbers0) > 8:
                # Compute the first 8 trials and the last 8 separately, before
                # taking the average to be consistent with asymmetric trials
                Dur_arr0 = np.zeros((2,n_maps))
                Occ_arr0 = np.zeros((2,n_maps))
                TCo_arr0 = np.zeros((2,n_maps))
                TMx_arr0 = np.zeros((2,n_maps,n_maps))
                Ent_arr0 = np.zeros(2)
                trial_numbers_split = np.array_split(trial_numbers0, 2)
                # Array_split is used in cases where a trial might have been
                # dropped, hence the split is only 8/7
                for s in range(len(trial_numbers_split)):
                    ep_idx0 = events.loc[(events["Event_id"] == ev_idx)&
                                        (events["Trial_number"].isin(trial_numbers_split[s])),
                                        "Epoch_idx"]
                    m_labels0 = m_labels[ep_idx0,:]
                    m_labels0_flat = m_labels0.reshape(m_labels0.shape[0]*m_labels0.shape[1])
                    # Compute duration, occurrence and time covered
                    l_, d_ = microstate_run_length_encoding(m_labels0_flat)
                    # For each microstate
                    for ii in range(n_maps):
                        if np.isnan(np.nanmean(d_[l_==ii])):
                            # The specific microstate did not occur at all for this event
                            Dur_arr0[s,ii] = 0
                            Occ_arr0[s,ii] = 0
                            TCo_arr0[s,ii] = 0
                        else:
                            Dur_arr0[s,ii] = np.mean(d_[l_==ii]) * 1000/sfreq # convert to ms
                            Occ_arr0[s,ii] = len(d_[l_==ii])/len(d_) * sfreq
                            TCo_arr0[s,ii] = np.sum(d_[l_==ii])/np.sum(d_)
                    # Compute transition matrix
                    TMx_arr0[s] = T_empirical(m_labels0_flat, n_maps)
                    # Compute Shannon Entropy relative to max
                    Ent_arr0[s] = H_1(m_labels0_flat, n_maps)/np.log(float(n_maps))
                # Average over the splits
                Dur_arr[e] = np.mean(Dur_arr0, axis=0)
                Occ_arr[e] = np.mean(Occ_arr0, axis=0)
                TCo_arr[e] = np.mean(TCo_arr0, axis=0)
                TMx_arr[e] = np.mean(TMx_arr0, axis=0)
                Ent_arr[e] = np.mean(Ent_arr0, axis=0)
                
            else:
                m_labels0 = m_labels[ep_idx,:]
                m_labels0_flat = m_labels0.reshape(m_labels0.shape[0]*m_labels0.shape[1])
                # Compute duration, occurrence and time covered
                l_, d_ = microstate_run_length_encoding(m_labels0_flat)
                # For each microstate
                for ii in range(n_maps):
                    if np.isnan(np.nanmean(d_[l_==ii])):
                        # The specific microstate did not occur at all for this event
                        Dur_arr[e,ii] = 0
                        Occ_arr[e,ii] = 0
                        TCo_arr[e,ii] = 0
                    else:
                        Dur_arr[e,ii] = np.mean(d_[l_==ii]) * 1000/sfreq # convert to ms
                        Occ_arr[e,ii] = len(d_[l_==ii])/len(d_) * sfreq
                        TCo_arr[e,ii] = np.sum(d_[l_==ii])/np.sum(d_)
                # Compute transition matrix
                TMx_arr[e] = T_empirical(m_labels0_flat, n_maps)
                # Compute Shannon Entropy relative to max
                Ent_arr[e] = H_1(m_labels0_flat, n_maps)/np.log(float(n_maps))
        # Combine all features in a list
        MFeatures = [Dur_arr,Occ_arr,TCo_arr,TMx_arr,Ent_arr]
        
        # Update in v7 - Collapse after computations in single events
        Dur_arr2 = np.zeros((len(collapsed_event_id),n_maps))
        Occ_arr2 = np.zeros((len(collapsed_event_id),n_maps))
        TCo_arr2 = np.zeros((len(collapsed_event_id),n_maps))
        TMx_arr2 = np.zeros((len(collapsed_event_id),n_maps,n_maps))
        Ent_arr2 = np.zeros(len(collapsed_event_id))
        MFeatures2 = [Dur_arr2,Occ_arr2,TCo_arr2,TMx_arr2,Ent_arr2]
        
        for f in range(len(MFeatures2)):
            tmp_feat = MFeatures[f]
            tmp_feat2 = MFeatures2[f]
            for e in range(len(collapsed_event_id)):
                ee = list(collapsed_event_id.keys())[e]
                if (ee == 'observer_actor'):
                    old_ev_idx1 = list(event_id.keys()).index("observe, actor")
                    old_ev_idx2 = list(event_id.keys()).index("observe, observer")
                    tmp_feat2[e] = np.mean([tmp_feat[old_ev_idx1],tmp_feat[old_ev_idx2]], axis=0)
                elif ee == 'follower_leader':
                    old_ev_idx1 = list(event_id.keys()).index("imitate, leader")
                    old_ev_idx2 = list(event_id.keys()).index("imitate, follower")
                    tmp_feat2[e] = np.mean([tmp_feat[old_ev_idx1],tmp_feat[old_ev_idx2]], axis=0)
                else:
                    new_ev_idx = list(collapsed_event_id.values())[e]
                    old_ev_idx = list(event_id.values()).index(new_ev_idx)
                    tmp_feat2[e] = tmp_feat[old_ev_idx]
        
        return m_labels, events, MFeatures2
    
    def load_microstate_arrays(i, standardize:bool=True):
        """ Loads the EEG data and convert to np.array """
        epoch1, trialinfo1 = load_epoch_from_fieldtrip(i)
        # Ensure it is averaged referenced
        epoch1 = epoch1.set_eeg_reference("average")
        # Get numpy arrays
        micro_data1 = epoch1.get_data() # get data
        if standardize == True:
            # Standardize data - will make EEG amplitudes comparable between
            # recordings/trials. Standardized along the ch axis, meaning the mean
            # and std of the topomap at each timepoint is normalized
            # It will be normalized separately for each participant
            micro_data10 = micro_data1 - micro_data1.mean(axis=1, keepdims=True)
            micro_data10 /= micro_data10.std(axis=1, keepdims=True)
        elif standardize == False:
            micro_data10 = micro_data1
        return micro_data10, trialinfo1
    
    def get_synch_events_from_pseudo_pairs(trialinfo1, trialinfo2, event_id):
        """
        Takes two subject indices from a pseudo pair and align their events (based
        on when they occur, e.g. first coupled in ppn1 with first coupled in ppn2)
        After their microstates have been synchronized, the microstate features
        are computed to serve as baseline as a non-interacting pair
        """
        # Synchronize the events from the pair based on timing info
        # By trimming the epochs to only include the epochs that are present
        # in both participants of the pair
        event_id_inv = {v: k for k, v in event_id.items()} # Inverse the event id
        events1 = pd.DataFrame(trialinfo1,columns=["Event_id","Trial_number","Trial_start_time"])
        events1["Event_label"] = events1["Event_id"].replace(event_id_inv)
        events1 = events1.reset_index().rename(columns={"index":"Epoch_idx"})
        
        events2 = pd.DataFrame(trialinfo2,columns=["Event_id","Trial_number","Trial_start_time"])
        events2["Event_label"] = events2["Event_id"].replace(event_id_inv)
        events2 = events2.reset_index().rename(columns={"index":"Epoch_idx"})
        
        # I should just make sure to align the timings of the trials for each event type
        counter_idx = 0
        # Initialize events as dataframe filled with unique number
        sync_events1 = events1.applymap(lambda x: 9999)
        sync_events2 = events2.applymap(lambda x: 9999)
    
        for e in range(len(event_id)):
            # Get event number for the specific event
            ev_idx = list(event_id.values())[e]
            # Get the indices corresponding to the event for each participant
            ep_idx1 = events1["Epoch_idx"][events1["Event_id"] == ev_idx]
            ep_idx2 = events2["Epoch_idx"][events2["Event_id"] == ev_idx]
            # Get the first trial from ppn1 and ppn2 that correspond to the event
            event_trial_idx1 = np.unique(trialinfo1[ep_idx1][:,1])
            event_trial_idx2 = np.unique(trialinfo2[ep_idx2][:,1])
            if len(event_trial_idx1) != len(event_trial_idx2):
                print(f"Not equal number of {list(event_id.keys())[e]} trials")
                min_trials = np.min([len(event_trial_idx1),len(event_trial_idx2)])
                print(f"Using the first {min_trials} trials")
                event_trial_idx1 = event_trial_idx1[:min_trials]
                event_trial_idx2 = event_trial_idx2[:min_trials]
    
            for t1, t2 in zip(event_trial_idx1,event_trial_idx2):
                t_idx1 = np.where(trialinfo1[:,1] == t1)[0]
                t_idx2 = np.where(trialinfo2[:,1] == t2)[0]
                # Get the timings for the epochs for the specific trial t
                t_timings1 = trialinfo1[t_idx1,2]
                t_timings2 = trialinfo2[t_idx2,2]
                timings_intersect = np.intersect1d(t_timings1,t_timings2)
                # Get the indices where the timings matches
                t_idx_match1 = t_idx1[pd.Series(t_timings1).isin(timings_intersect)]
                t_idx_match2 = t_idx2[pd.Series(t_timings2).isin(timings_intersect)]
                # Save the event info
                sync_events1.iloc[counter_idx:(counter_idx+len(timings_intersect))] = events1.iloc[t_idx_match1]
                sync_events2.iloc[counter_idx:(counter_idx+len(timings_intersect))] = events2.iloc[t_idx_match2]
                # Move counter
                counter_idx += len(timings_intersect)
        # Find the epochs that were asynchronous, which have to be trimmed
        asynch_epochs1 = np.unique(np.where(sync_events1==9999)[0])
        asynch_epochs2 = np.unique(np.where(sync_events2==9999)[0])
        # Fix events
        sync_events1 = sync_events1.drop(asynch_epochs1,axis=0).reset_index(drop=True)
        sync_events2 = sync_events2.drop(asynch_epochs2,axis=0).reset_index(drop=True)
        assert len(sync_events1) == len(sync_events2) # check the amount of synchronized epochs are equal
        return [sync_events1, sync_events2]
    
    def combine_two_person_microstate_arrays(micro_data1, micro_data2, events, sfreq, freq_range=None):
        """ Combine two EEG array data into one array for microstate estimation """
        # Get the synchronized epochs
        micro_data1 = micro_data1[events[0]["Epoch_idx"]]
        micro_data2 = micro_data2[events[1]["Epoch_idx"]]
        # Get event labels for each epoch
        # Since they are already synchronized I just take the first
        event_labels0 = events[0]["Event_label"]
        # Concatenate along the channel axis, to treat the paired EEG as one entity
        micro_data = np.concatenate([micro_data1,micro_data2],axis=1)
        
        ### Update in v4.1
        # Flip the micro_data between ppn1 and ppn2 to fix first microstate
        # always being observer/follower, and second microstate being actor/leader
        cond6_idx = event_labels0 == "observe, observer"