KIAPBCI/kiap_bci/analytics/lfp_analysis3.py

'''
description: offline analysis of kiap data
author: Ioannis Vlachos
date:   16.04.2019

Copyright (c) 2019 Ioannis Vlachos.
All rights reserved.'''

import re
import os
import glob
import importlib
import urllib
import time

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy import signal
from scipy import stats
# from cachetools import cached, LRUCache, TTLCache
from neo.io.blackrockio import BlackrockIO
# from mne.connectivity import spectral_connectivity

import aux
from analytics import analytics1
importlib.reload(analytics1)



# @cached(cache={})
def read_data():
    '''returns first array1 and then array2 channels'''

    data = reader.get_analogsignal_chunk(i_stop=eval(str(params.lfp.i_stop)),
                                channel_indexes=params.lfp.array1 + params.lfp.array2)

    return data



params = aux.load_config()
ch_nr1 = len(params.lfp.array1)
ch_nr2 = len(params.lfp.array2)
psth_win = params.lfp.psth_win

# file_name = '/home/vlachos/data_kiap/20190321-135311/20190321-135311-001'
# path = '/home/kiap/Documents/kiap/k01/20190414-123429/'
path = '/media/vlachos/kiap_backup/Recordings/K01/Recordings/20190414-123429/'
file_name0 = '20190414-123429-001.ns2'
file_name = path + file_name0

plt.ioff()

# if not 'data0' in locals():

Sxx1_tot = np.zeros((999, ch_nr1, 479, 75))
for fid in range(1, 76):
# for fid in range(1,76):
# for fid in range(62,63):

    file_name0 = f'20190414-123429-{fid:03d}.ns2'
    file_name = path + file_name0
    print(file_name)

    reader = BlackrockIO(filename=file_name)
    data0 = read_data()  # (samples, channels)

    data0 = data0.copy().astype(np.float64)
    # trials_msk = get_valid_trials(triggers_no_ns2, triggers_yes_ns2)

    array1 = params.lfp.array1
    array2 = params.lfp.array2
    # data1[data1 > params.lfp.artifact_thr] = 0
    # data1[data1 < -params.lfp.artifact_thr] = 0

    # data0 = data0 - np.mean(data0, 1)[:, np.newaxis]

    # reader.nev_data['Comments'][0]

    data0 = stats.zscore(data0)
    data01 = data0[:, :len(array1)]
    data02 = data0[:, len(array1):len(array1) + len(array2)]


    data1 = np.zeros((data01.shape[0], data01.shape[1], 3))     # (samples, channels, 3)
    data2 = np.zeros((data02.shape[0], data02.shape[1], 3))     # (samples, channels, 3)
    # Y1 = np.zeros((data01.shape[0], data01.shape[1], 3))

    # data1[:, :, 0], _, _ = analytics1.bw_filter(data01.T, params.lfp.filter_fc_lb, params.lfp.fs, params.lfp.filter_order_lb, plot=params.lfp.plot.filters)
    # data1[:, :, 1], _, _ = analytics1.bw_filter(data01.T, params.lfp.filter_fc_mb, params.lfp.fs, params.lfp.filter_order_mb, plot=params.lfp.plot.filters)
    # data1[:, :, 2], _, _ = analytics1.bw_filter(data01.T, params.lfp.filter_fc_hb, params.lfp.fs, params.lfp.filter_order_hb, plot=params.lfp.plot.filters)

    # data2[:, :, 0], _, _ = analytics1.bw_filter(data02.T, params.lfp.filter_fc_lb, params.lfp.fs, params.lfp.filter_order_lb, plot=params.lfp.plot.filters)
    # data2[:, :, 1], _, _ = analytics1.bw_filter(data02.T, params.lfp.filter_fc_mb, params.lfp.fs, params.lfp.filter_order_mb, plot=params.lfp.plot.filters)
    # data2[:, :, 2], _, _ = analytics1.bw_filter(data02.T, params.lfp.filter_fc_hb, params.lfp.fs, params.lfp.filter_order_hb, plot=params.lfp.plot.filters)

    # Y01, ff01 = analytics1.calc_fft(data01, params.lfp.fs, i_stop=1e5, axis=0)
    # Y02, ff02 = analytics1.calc_fft(data02, params.lfp.fs, i_stop=1e5, axis=0)
    # Y1, ff1 = analytics1.calc_fft(data1, params.lfp.fs, i_stop=1e5, axis=0)
    # Y2, ff2 = analytics1.calc_fft(data2, params.lfp.fs, i_stop=1e5, axis=0)

    #  SPECTROGRAMS
    # ----------------------------------

    nperseg = 5000

    ff, tt, Sxx1 = signal.spectrogram(data01, params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.5))
    # ff, tt, Sxx2 = signal.spectrogram(data02, params.lfp.fs, axis=0, nperseg=nperseg, noverlap=int(nperseg * 0.5))

    print(f'{fid}, Sxx1 shape: {Sxx1.shape}, {Sxx1_tot.shape}')
    f_idx = (ff>0) & (ff<200)
    Sxx1_tot[:, :, :Sxx1.shape[2], fid - 1] = Sxx1[f_idx, :]
    # fig1 = plt.figure(1, figsrize=(14,10))

    # plt.clf()


    # for ii in range(1):
    #     # plt.plot()
    #     ax = plt.subplot(8,8,ii+1)
    #     # plt.plot(ff01[fidx], Y01[fidx])
    #     plt.pcolormesh(tt, ff[f_idx], np.log10(Sxx1[f_idx, ii, :]))
    #     # plt.pcolormesh(tt, ff[f_idx], np.log10(Sxx2[f_idx, ii, :]))
    #     # plt.pcolormesh(tt, ff[f_idx], np.log10(Sxx1[f_idx, ii, :]/Sxx1[f_idx,:,:].mean(1)))
    #     print(ii)
    #     if ii<56:
    #         ax.set_xticks([])
    #     if np.mod(ii,8)!=0:
    #         ax.set_yticks([])

    # fname = path + os.path.splitext(file_name0)[0]+'-arr3.png'
    # fig1.savefig(fname)
# plt.show()


# con, freqs, times, epochs, tapers = mne.connectivity.spectral_connectivity(
        # sig.values[None, :, :], sfreq=freq, indices=indices)

# save individual plots
fig1 = plt.figure(1)
plt.ioff()
for ii in range(4):
    aa = Sxx1_tot[:, ii, :, :]
    # bb = np.reshape(aa,(1999,239*75), order='F')
    bb = np.reshape(aa, (999, 479 * 75), order='F')
    plt.clf()
    # plt.pcolormesh(np.arange(239*75),ff[f_idx],np.log10(bb))
    plt.pcolormesh(np.arange(479 * 75), ff[f_idx], np.log10(bb))
    fname = path + f'total/ch_{ii}_2.png'
    fig1.savefig(fname)
    print(ii)