123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155 |
- #!/usr/bin/env python
- # coding: utf-8
- database_path = '/media/andrey/My Passport/GIN_new/Anesthesia_CA1/meta_data/meta_recordings - anesthesia.xlsx'
- # ### Select the range of recordings for the analysis (see "Number" row in the meta data file)
- # In[4]:
- rec = [x for x in range(0,188)]
- # In[1]:
- import numpy as np
- import numpy.ma as ma
- import matplotlib.pyplot as plt
- import matplotlib.ticker as ticker
- import pandas as pd
- import seaborn as sns
- import pickle
- import os
- sns.set()
- sns.set_style("whitegrid")
- from scipy.signal import medfilt
- from scipy.stats import skew, kurtosis, zscore
- from scipy import signal
- from sklearn.linear_model import LinearRegression, TheilSenRegressor
- plt.rcParams['figure.figsize'] = [8, 8]
- # In[2]:
- from capipeline import *
- motion_index = np.zeros((len(rec),100000),dtype='float')
- recording_length = np.zeros((len(rec)),dtype='int')
- number_of_quiet_periods = 6
- start_quiet_period = np.zeros((len(rec),number_of_quiet_periods),dtype='int')
- stop_quiet_period = np.zeros((len(rec),number_of_quiet_periods),dtype='int')
- n_quiet_periods_in_recording = np.zeros((len(rec)),dtype='int')
- for i, r in enumerate(rec):
- print("Reconrding # ", r)
-
- animal = get_animal_from_recording(r, database_path)
-
- #condition = get_condition(r, database_path)
-
- #print("#" + str(r) + " " + str(animal) + " " + str(condition) + " ")
- meta_data = pd.read_excel(database_path)
-
- path_excel_rec = str(meta_data['Folder'][r]) + str(meta_data['Subfolder'][r]) + str(meta_data['Recording idx'][r]) + '/suite2p/'
- stat = np.load(path_excel_rec+ '/plane0/ops.npy', allow_pickle=True)
- #plt.plot(stat.item(0)['yoff'],alpha=0.5)
- #plt.plot(stat.item(0)['xoff'],alpha=0.5)
- motion_index[i,:len(stat.item(0)['yoff'])] = np.sqrt(stat.item(0)['yoff']**2 + stat.item(0)['xoff']**2)
- print("Min motion index:", min(motion_index[i,:]))
- print("Max motion index:", max(motion_index[i,:]))
- recording_length[i] = len(stat.item(0)['yoff'])
- print("Recording length:", recording_length[i])
-
- if (type(meta_data['Quiet periods'][r]) == type('str')):
- print("Quiet periods: ", meta_data['Quiet periods'][r].split(','))
- n_quiet_periods_in_recording[i] = int(len(meta_data['Quiet periods'][r].split(','))/2)
- for k in range(n_quiet_periods_in_recording[i]):
- start_quiet_period[i,k] = int(meta_data['Quiet periods'][r].split(',')[k*2])
- stop_quiet_period[i,k] = int(meta_data['Quiet periods'][r].split(',')[k*2+1])
- print(start_quiet_period[i,k])
- print(stop_quiet_period[i,k] )
- else:
- print("Nan for recording:", r)
- start_quiet_period[i,0] = 0
- stop_quiet_period[i,0] = recording_length[i]
- #np.save("./xy-motion.npy",motion_index)
- #mi = motion_index
- mi = np.load("./xy-motion.npy")
- print(np.max(mi))
- print(np.min(mi))
- mi_av = np.mean(mi[:,:18000].reshape(mi.shape[0],60,300), axis=2)
- plt.rcParams["axes.grid"] = False
- plt.figure(figsize = (10,15))
- pos = plt.imshow(mi_av,cmap='Purples',vmin = 0, vmax = 10) #,aspect='equal'
- plt.colorbar(pos)
- for i in range(len(rec)):
- for k in range(int(n_quiet_periods_in_recording[i])):
- line = plt.scatter(start_quiet_period[i,k]/300,i,marker='>',color='k',s = 10,label='start of quiet period')
- line.set_clip_on(False)
- line = plt.scatter(stop_quiet_period[i,k]/300,i,marker='<',color='k',s = 10, label='end of quiet period')
- line.set_clip_on(False)
- line = plt.scatter(recording_length[0:len(rec)]/300, np.arange(len(rec)) ,marker='|',color='k', s = 10, label='end of the recording')
- line.set_clip_on(False)
- #plt.legend(loc='upper right')
- plt.xlabel('x 300 frames')
- plt.ylabel('recording')
- plt.title('Anesthesia dataset motion validation')
- #plt.gcf().set_facecolor("white")
- plt.xlim([0,60])
- line.set_clip_on(False)
- plt.savefig("Validation_motion.png")
- plt.savefig("Validation_motion.svg")
- plt.show()
|