From 038e04c4164c006c06ea8d0846c01e30bb743838 Mon Sep 17 00:00:00 2001 From: jianyoulin Date: Mon, 10 Jun 2019 11:55:36 -0400 Subject: [PATCH 1/3] merge splitted clusters --- blech_post_process.py | 162 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 140 insertions(+), 22 deletions(-) diff --git a/blech_post_process.py b/blech_post_process.py index 68261ad2..5a8759ab 100644 --- a/blech_post_process.py +++ b/blech_post_process.py @@ -106,9 +106,11 @@ class unit_descriptor(tables.IsDescription): # If the user asked to split/re-cluster, ask them for the clustering parameters and perform clustering split_predictions = [] chosen_split = 0 - if re_cluster: + if re_cluster: #I want to split this cluster # Get clustering parameters from user - clustering_params = easygui.multenterbox(msg = 'Fill in the parameters for re-clustering (using a GMM)', fields = ['Number of clusters', 'Maximum number of iterations (1000 is more than enough)', 'Convergence criterion (usually 0.0001)', 'Number of random restarts for GMM (10 is more than enough)']) + clustering_params = easygui.multenterbox(msg = 'Fill in the parameters for re-clustering (using a GMM)', + fields = ['Number of clusters', 'Maximum number of iterations (1000 is more than enough)', + 'Convergence criterion (usually 0.0001)', 'Number of random restarts for GMM (10 is more than enough)']) n_clusters = int(clustering_params[0]) n_iter = int(clustering_params[1]) thresh = float(clustering_params[2]) @@ -152,10 +154,15 @@ class unit_descriptor(tables.IsDescription): # Ask the user for the split clusters they want to choose chosen_split = easygui.multchoicebox(msg = 'Which split cluster do you want to choose? Hit cancel to exit', choices = tuple([str(i) for i in range(n_clusters)])) try: - chosen_split = int(chosen_split[0]) + chosen_split = [int(chosen_split[i]) for i in range(len(chosen_split))] + split_merge = False +# if len(chosen_split) > 1: +# split_merge = easygui.multchoicebox(msg = 'I want to merge these splited-clusters into one unit (True = Yes, False = No)', choices = ('True', 'False')) +# split_merge = ast.literal_eval(split_merge[0]) except: continue + # Get list of existing nodes/groups under /sorted_units node_list = hf5.list_nodes('/sorted_units') @@ -179,25 +186,135 @@ class unit_descriptor(tables.IsDescription): # If the user re-clustered/split clusters, add the chosen clusters in split_clusters if re_cluster: - hf5.create_group('/sorted_units', unit_name) - unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster - unit_waveforms = unit_waveforms[np.where(split_predictions == chosen_split)[0], :] # Subsetting this set of waveforms to include only the chosen split - unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] # Do the same thing for the spike times - unit_times = unit_times[np.where(split_predictions == chosen_split)[0]] - waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) - times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) - unit_description['electrode_number'] = electrode_num - single_unit = easygui.multchoicebox(msg = 'I am almost-SURE that this is a beautiful single unit (True = Yes, False = No)', choices = ('True', 'False')) - unit_description['single_unit'] = int(ast.literal_eval(single_unit[0])) - # If the user says that this is a single unit, ask them whether its regular or fast spiking - unit_description['regular_spiking'] = 0 - unit_description['fast_spiking'] = 0 - if int(ast.literal_eval(single_unit[0])): - unit_type = easygui.multchoicebox(msg = 'What type of unit is this (Regular spiking = Pyramidal cells, Fast spiking = PV+ interneurons)?', choices = ('regular_spiking', 'fast_spiking')) - unit_description[unit_type[0]] = 1 - unit_description.append() - table.flush() - hf5.flush() +# split_merge = False + if len(chosen_split) == 1: + hf5.create_group('/sorted_units', unit_name) + unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster + unit_waveforms = unit_waveforms[np.where(split_predictions == chosen_split[0])[0], :] # Subsetting this set of waveforms to include only the chosen split + unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] # Do the same thing for the spike times + unit_times = unit_times[np.where(split_predictions == chosen_split[0])[0]] + waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) + times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) + unit_description['electrode_number'] = electrode_num + single_unit = easygui.multchoicebox(msg = 'I am almost-SURE that this is a beautiful single unit (True = Yes, False = No)', choices = ('True', 'False')) + unit_description['single_unit'] = int(ast.literal_eval(single_unit[0])) + # If the user says that this is a single unit, ask them whether its regular or fast spiking + unit_description['regular_spiking'] = 0 + unit_description['fast_spiking'] = 0 + if int(ast.literal_eval(single_unit[0])): + unit_type = easygui.multchoicebox(msg = 'What type of unit is this (Regular spiking = Pyramidal cells, Fast spiking = PV+ interneurons)?', choices = ('regular_spiking', 'fast_spiking')) + unit_description[unit_type[0]] = 1 + unit_description.append() + table.flush() + hf5.flush() + + elif len(chosen_split) > 1: + split_merge = easygui.multchoicebox(msg = 'I want to merge these splited-clusters into one unit (True = Yes, False = No)', choices = ('True', 'False')) + split_merge = ast.literal_eval(split_merge[0]) + + if split_merge: #I want to merge these splited-clusters + unit_waveforms = [] + unit_times = [] + cluster_unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] + cluster_unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] + for split_cluster in chosen_split: + if unit_waveforms == []: + unit_waveforms = cluster_unit_waveforms[np.where(split_predictions == split_cluster)[0], :] + unit_times = cluster_unit_times[np.where(split_predictions == split_cluster)[0]] + else: + unit_waveforms = np.concatenate((unit_waveforms, cluster_unit_waveforms[np.where(split_predictions == split_cluster)[0], :])) + unit_times = np.concatenate((unit_times, cluster_unit_times[np.where(split_predictions == split_cluster)[0]])) + + # Show the merged cluster to the user, and ask if they still want to merge + x = np.arange(len(unit_waveforms[0])/10) + 1 + fig, ax = blech_waveforms_datashader.waveforms_datashader(unit_waveforms, x) + # plt.plot(x - 15, unit_waveforms[:, ::10].T, linewidth = 0.01, color = 'red') + ax.set_xlabel('Sample (30 samples per ms)') + ax.set_ylabel('Voltage (microvolts)') + ax.set_title('Merged cluster, No. of waveforms={:d}'.format(unit_waveforms.shape[0])) + plt.show() + + # Warn the user about the frequency of ISI violations in the merged unit + ISIs = np.ediff1d(np.sort(unit_times))/30.0 + violations1 = 100.0*float(np.sum(ISIs < 1.0)/len(unit_times)) + violations2 = 100.0*float(np.sum(ISIs < 2.0)/len(unit_times)) + proceed = easygui.multchoicebox(msg = 'My merged cluster has %.1f percent (<2ms) and %.1f percent (<1ms) ISI violations out of %i total waveforms. I want to still merge these clusters into one unit (True = Yes, False = No)' % (violations2, violations1, len(unit_times)), choices = ('True', 'False')) + proceed = ast.literal_eval(proceed[0]) + + # Create unit if the user agrees to proceed, else include each split_cluster as a separate unit + if proceed: + hf5.create_group('/sorted_units', unit_name) + waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) + times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) + unit_description['electrode_number'] = electrode_num + single_unit = easygui.multchoicebox(msg = 'I am almost-SURE that this is a beautiful single unit (True = Yes, False = No)', choices = ('True', 'False')) + unit_description['single_unit'] = int(ast.literal_eval(single_unit[0])) + # If the user says that this is a single unit, ask them whether its regular or fast spiking + unit_description['regular_spiking'] = 0 + unit_description['fast_spiking'] = 0 + if int(ast.literal_eval(single_unit[0])): + unit_type = easygui.multchoicebox(msg = 'What type of unit is this (Regular spiking = Pyramidal cells, Fast spiking = PV+ interneurons)?', choices = ('regular_spiking', 'fast_spiking')) + unit_description[unit_type[0]] = 1 + unit_description.append() + table.flush() + hf5.flush() +############ + else: # if not merge, then include each split_cluster as a separate unit + for split_cluster in chosen_split: + hf5.create_group('/sorted_units', unit_name) + unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster + unit_waveforms = unit_waveforms[np.where(split_predictions == chosen_split)[0], :] # Subsetting this set of waveforms to include only the chosen split + unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] # Do the same thing for the spike times + unit_times = unit_times[np.where(split_predictions == chosen_split)[0]] + waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) + times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) + unit_description['electrode_number'] = electrode_num + single_unit = easygui.multchoicebox(msg = 'I am almost-SURE that electrode: %i cluster: %i split: %i is a beautiful single unit (True = Yes, False = No)' % (electrode_num, int(cluster), int(split_cluster)), choices = ('True', 'False')) + unit_description['single_unit'] = int(ast.literal_eval(single_unit[0])) + # If the user says that this is a single unit, ask them whether its regular or fast spiking + unit_description['regular_spiking'] = 0 + unit_description['fast_spiking'] = 0 + if int(ast.literal_eval(single_unit[0])): + unit_type = easygui.multchoicebox(msg = 'What type of unit is this (Regular spiking = Pyramidal cells, Fast spiking = PV+ interneurons)?', choices = ('regular_spiking', 'fast_spiking')) + unit_description[unit_type[0]] = 1 + unit_description.append() + table.flush() + hf5.flush() + + # Finally increment max_unit and create a new unit name + max_unit += 1 + unit_name = 'unit%03d' % int(max_unit + 1) + + # Get a new unit_descriptor table row for this new unit + unit_description = table.row + #continue +########### + else: + continue +# ============================================================================= +# +# else:#if not split_merge: +# hf5.create_group('/sorted_units', unit_name) +# unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster +# unit_waveforms = unit_waveforms[np.where(split_predictions == chosen_split[0])[0], :] # Subsetting this set of waveforms to include only the chosen split +# unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] # Do the same thing for the spike times +# unit_times = unit_times[np.where(split_predictions == chosen_split[0])[0]] +# waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) +# times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) +# unit_description['electrode_number'] = electrode_num +# single_unit = easygui.multchoicebox(msg = 'I am almost-SURE that this is a beautiful single unit (True = Yes, False = No)', choices = ('True', 'False')) +# unit_description['single_unit'] = int(ast.literal_eval(single_unit[0])) +# # If the user says that this is a single unit, ask them whether its regular or fast spiking +# unit_description['regular_spiking'] = 0 +# unit_description['fast_spiking'] = 0 +# if int(ast.literal_eval(single_unit[0])): +# unit_type = easygui.multchoicebox(msg = 'What type of unit is this (Regular spiking = Pyramidal cells, Fast spiking = PV+ interneurons)?', choices = ('regular_spiking', 'fast_spiking')) +# unit_description[unit_type[0]] = 1 +# unit_description.append() +# table.flush() +# hf5.flush() +# +# ============================================================================= # If only 1 cluster was chosen (and it wasn't split), add that as a new unit in /sorted_units. Ask if the isolated unit is an almost-SURE single unit @@ -220,6 +337,7 @@ class unit_descriptor(tables.IsDescription): table.flush() hf5.flush() + else: # If the chosen units are going to be merged, merge them if merge: From 6fb1dd82d428125cb7b99f676ba2aaec4b1b7a8a Mon Sep 17 00:00:00 2001 From: jianyoulin Date: Tue, 18 Jun 2019 12:46:42 -0400 Subject: [PATCH 2/3] correct typos --- blech_post_process.py | 78 +++++++++++++------------------------------ 1 file changed, 24 insertions(+), 54 deletions(-) diff --git a/blech_post_process.py b/blech_post_process.py index 5a8759ab..06292795 100644 --- a/blech_post_process.py +++ b/blech_post_process.py @@ -70,6 +70,27 @@ class unit_descriptor(tables.IsDescription): # Run an infinite loop as long as the user wants to pick clusters from the electrodes while True: + # Get list of existing nodes/groups under /sorted_units + node_list = hf5.list_nodes('/sorted_units') + + # If node_list is empty, start naming units from 000 + unit_name = '' + max_unit = 0 + if node_list == []: + unit_name = 'unit%03d' % 0 + # Else name the new unit by incrementing the last unit by 1 + else: + unit_numbers = [] + for node in node_list: + unit_numbers.append(node._v_pathname.split('/')[-1][-3:]) + unit_numbers[-1] = int(unit_numbers[-1]) + unit_numbers = np.array(unit_numbers) + max_unit = np.max(unit_numbers) + unit_name = 'unit%03d' % int(max_unit + 1) + + # Get a new unit_descriptor table row for this new unit + unit_description = table.row + # Get electrode number from user electrode_num = easygui.multenterbox(msg = 'Which electrode do you want to choose? Hit cancel to exit', fields = ['Electrode #']) # Break if wrong input/cancel command was given @@ -152,41 +173,15 @@ class unit_descriptor(tables.IsDescription): plt.show() # Ask the user for the split clusters they want to choose - chosen_split = easygui.multchoicebox(msg = 'Which split cluster do you want to choose? Hit cancel to exit', choices = tuple([str(i) for i in range(n_clusters)])) + chosen_split = easygui.multchoicebox(msg = 'Which split clusters do you want to choose? Hit cancel to exit', choices = tuple([str(i) for i in range(n_clusters)])) try: chosen_split = [int(chosen_split[i]) for i in range(len(chosen_split))] split_merge = False -# if len(chosen_split) > 1: -# split_merge = easygui.multchoicebox(msg = 'I want to merge these splited-clusters into one unit (True = Yes, False = No)', choices = ('True', 'False')) -# split_merge = ast.literal_eval(split_merge[0]) except: continue - # Get list of existing nodes/groups under /sorted_units - node_list = hf5.list_nodes('/sorted_units') - - # If node_list is empty, start naming units from 000 - unit_name = '' - max_unit = 0 - if node_list == []: - unit_name = 'unit%03d' % 0 - # Else name the new unit by incrementing the last unit by 1 - else: - unit_numbers = [] - for node in node_list: - unit_numbers.append(node._v_pathname.split('/')[-1][-3:]) - unit_numbers[-1] = int(unit_numbers[-1]) - unit_numbers = np.array(unit_numbers) - max_unit = np.max(unit_numbers) - unit_name = 'unit%03d' % int(max_unit + 1) - - # Get a new unit_descriptor table row for this new unit - unit_description = table.row - # If the user re-clustered/split clusters, add the chosen clusters in split_clusters - if re_cluster: -# split_merge = False if len(chosen_split) == 1: hf5.create_group('/sorted_units', unit_name) unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster @@ -263,9 +258,9 @@ class unit_descriptor(tables.IsDescription): for split_cluster in chosen_split: hf5.create_group('/sorted_units', unit_name) unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster - unit_waveforms = unit_waveforms[np.where(split_predictions == chosen_split)[0], :] # Subsetting this set of waveforms to include only the chosen split + unit_waveforms = unit_waveforms[np.where(split_predictions == split_cluster)[0], :] # Subsetting this set of waveforms to include only the chosen split unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] # Do the same thing for the spike times - unit_times = unit_times[np.where(split_predictions == chosen_split)[0]] + unit_times = unit_times[np.where(split_predictions == split_cluster)[0]] waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) unit_description['electrode_number'] = electrode_num @@ -291,31 +286,6 @@ class unit_descriptor(tables.IsDescription): ########### else: continue -# ============================================================================= -# -# else:#if not split_merge: -# hf5.create_group('/sorted_units', unit_name) -# unit_waveforms = spike_waveforms[np.where(predictions == int(clusters[0]))[0], :] # Waveforms of originally chosen cluster -# unit_waveforms = unit_waveforms[np.where(split_predictions == chosen_split[0])[0], :] # Subsetting this set of waveforms to include only the chosen split -# unit_times = spike_times[np.where(predictions == int(clusters[0]))[0]] # Do the same thing for the spike times -# unit_times = unit_times[np.where(split_predictions == chosen_split[0])[0]] -# waveforms = hf5.create_array('/sorted_units/%s' % unit_name, 'waveforms', unit_waveforms) -# times = hf5.create_array('/sorted_units/%s' % unit_name, 'times', unit_times) -# unit_description['electrode_number'] = electrode_num -# single_unit = easygui.multchoicebox(msg = 'I am almost-SURE that this is a beautiful single unit (True = Yes, False = No)', choices = ('True', 'False')) -# unit_description['single_unit'] = int(ast.literal_eval(single_unit[0])) -# # If the user says that this is a single unit, ask them whether its regular or fast spiking -# unit_description['regular_spiking'] = 0 -# unit_description['fast_spiking'] = 0 -# if int(ast.literal_eval(single_unit[0])): -# unit_type = easygui.multchoicebox(msg = 'What type of unit is this (Regular spiking = Pyramidal cells, Fast spiking = PV+ interneurons)?', choices = ('regular_spiking', 'fast_spiking')) -# unit_description[unit_type[0]] = 1 -# unit_description.append() -# table.flush() -# hf5.flush() -# -# ============================================================================= - # If only 1 cluster was chosen (and it wasn't split), add that as a new unit in /sorted_units. Ask if the isolated unit is an almost-SURE single unit elif len(clusters) == 1: From c6e394967d4b365cca50e4c921a7e602516df056 Mon Sep 17 00:00:00 2001 From: jianyoulin Date: Fri, 5 Jul 2019 14:33:57 -0400 Subject: [PATCH 3/3] fix a bug that takes tranisition time of palatability epoch within first 100 ms after taste delivery --- variational_HMM_line_up_palatability_plot.py | 154 +++++++++++++++++-- 1 file changed, 140 insertions(+), 14 deletions(-) diff --git a/variational_HMM_line_up_palatability_plot.py b/variational_HMM_line_up_palatability_plot.py index 04caea28..fb8b9852 100644 --- a/variational_HMM_line_up_palatability_plot.py +++ b/variational_HMM_line_up_palatability_plot.py @@ -6,7 +6,7 @@ import easygui import sys import os -import pylab as plt +#import pylab as plt from scipy.stats import rankdata from scipy.stats import spearmanr from scipy.stats import pearsonr @@ -56,6 +56,11 @@ p_spearman_laser_off = [] r_spearman_laser_on = [] p_spearman_laser_on = [] + +# Empty lists for switch times for all tastes +switchtimes_laser_off = [[] for i in range(num_tastes)] +switchtimes_laser_on = [[] for i in range(num_tastes)] + for dir_name in dirs: os.chdir(dir_name) # Locate the hdf5 file @@ -130,18 +135,23 @@ # Run through the laser off trials for trial in range(len(laser_off_trials)): # Find the time that the dominant state first comes on - state_onset = np.where(posterior_proba_laser_on[laser_on_pal_state, trial, :] > 0.5)[0] + state_onset = np.where(posterior_proba_laser_off[laser_off_pal_state, trial, :] > 0.5)[0] later_onset = np.where(np.ediff1d(state_onset) > 1)[0] # If the dominant state does go over 0.8 in probability during the trial, pick its onset if len(state_onset) > 0: if state_onset[0]*bin_size > 100: state_onset = state_onset[0] elif len(later_onset) > 0: - state_onset = later_onset[0] + if state_onset[later_onset[0]+1]*bin_size > 100: + state_onset = state_onset[later_onset[0]+1] + else: + continue else: continue print(state_onset) - # Append spiking data for a total of 2.5s - 1s before to 1.5s after the state onset + # Append switch time of this trial to the array + switchtimes_laser_off[dig_in.index(taste)].append(state_onset*bin_size) + # Append spiking data for a total of 3.5s - 1s before to 2.5s after the state onset this_taste_laser_off_aligned.append(taste.spike_array[laser_off_trials[trial], single_units, pre_stim + state_onset*bin_size + pre_stim_hmm - 1000:pre_stim + state_onset*bin_size + pre_stim_hmm + 2500]) # Append an equal 2.5s worth of unaligned spiking data this_taste_laser_off_unaligned.append(taste.spike_array[laser_off_trials[trial], single_units, pre_stim - 1000:pre_stim + 2500]) @@ -157,10 +167,15 @@ if state_onset[0]*bin_size > 100: state_onset = state_onset[0] elif len(later_onset) > 0: - state_onset = later_onset[0] + if state_onset[later_onset[0]+1]*bin_size > 100: + state_onset = state_onset[later_onset[0]+1] + else: + continue else: continue print(state_onset) + # Append switch time of this trial to the array + switchtimes_laser_on[dig_in.index(taste)].append(state_onset*bin_size) # Append spiking data for a total of 2.5s - 1s before to 1.5s after the state onset this_taste_laser_on_aligned.append(taste.spike_array[laser_on_trials[trial], single_units, pre_stim + state_onset*bin_size + pre_stim_hmm - 1000:pre_stim + state_onset*bin_size + pre_stim_hmm + 2500]) # Append an equal 2.5s worth of unaligned spiking data @@ -206,20 +221,39 @@ for time in range(response_laser_off_aligned.shape[0]): ranks = rankdata(response_laser_off_aligned[time, :, unit]) r, p = spearmanr(ranks, palatability_laser_off) - this_unit_r_laser_off[0].append(r) - this_unit_p_laser_off[0].append(p) + if np.isnan(r): + this_unit_r_laser_off[0].append(0.0) + this_unit_p_laser_off[0].append(1.0) + else: + this_unit_r_laser_off[0].append(r) + this_unit_p_laser_off[0].append(p) + ranks = rankdata(response_laser_off_unaligned[time, :, unit]) r, p = spearmanr(ranks, palatability_laser_off) - this_unit_r_laser_off[1].append(r) - this_unit_p_laser_off[1].append(p) + if np.isnan(r): + this_unit_r_laser_off[1].append(0.0) + this_unit_p_laser_off[1].append(1.0) + else: + this_unit_r_laser_off[1].append(r) + this_unit_p_laser_off[1].append(p) + ranks = rankdata(response_laser_on_aligned[time, :, unit]) r, p = spearmanr(ranks, palatability_laser_on) - this_unit_r_laser_on[0].append(r) - this_unit_p_laser_on[0].append(p) + if np.isnan(r): + this_unit_r_laser_on[0].append(0.0) + this_unit_p_laser_on[0].append(1.0) + else: + this_unit_r_laser_on[0].append(r) + this_unit_p_laser_on[0].append(p) + ranks = rankdata(response_laser_on_unaligned[time, :, unit]) r, p = spearmanr(ranks, palatability_laser_on) - this_unit_r_laser_on[1].append(r) - this_unit_p_laser_on[1].append(p) + if np.isnan(r): + this_unit_r_laser_on[1].append(0.0) + this_unit_p_laser_on[1].append(1.0) + else: + this_unit_r_laser_on[1].append(r) + this_unit_p_laser_on[1].append(p) # Append the unit specific lists to the main lists of r and p r_spearman_laser_off.append(this_unit_r_laser_off) @@ -229,8 +263,100 @@ # Close the hdf5 file hf5.close() +### +#r_spearman_laser_off[np.isnan(r_spearman_laser_off)] = 0 +#r_spearman_laser_on[np.isnan(r_spearman_laser_on)] = 0 r_spearman_laser_off = np.array(r_spearman_laser_off) p_spearman_laser_off = np.array(p_spearman_laser_off) r_spearman_laser_on = np.array(r_spearman_laser_on) -p_spearman_laser_on = np.array(p_spearman_laser_on) +p_spearman_laser_on = np.array(p_spearman_laser_on) + +########## +# Ask the user for the directory to save arrays and plots +print('$$$$Choose the output directory for restuls of HMM alignment coorelation analysis$$$$') +dir_name = easygui.diropenbox(msg = 'Choose the output directory for saving correlation results') +os.chdir(dir_name) +########## + +# Saving arrays +np.save('r_spearman_laser_off.npy', r_spearman_laser_off) +np.save('p_spearman_laser_off.npy', p_spearman_laser_off) +np.save('r_spearman_laser_on.npy', r_spearman_laser_on) +np.save('p_spearman_laser_on.npy', p_spearman_laser_on) + +# Plot results +import matplotlib.pyplot as plt +from scipy.ndimage.filters import gaussian_filter1d +import seaborn as sns +sns.set() + +x_axis = np.arange(0, laser_on_aligned.shape[-1]-params[0]+params[1], params[1]) - 1000 +shift_time=750 # the time shift to match aligned and unaligned palatability curve +aligned_plot_times = np.where((x_axis >= x_axis[0])*(x_axis <= 1500))[0] +unaligned_plot_times = np.where((x_axis >= x_axis[0] + shift_time)*(x_axis <= 1500 + shift_time))[0] +plot_times = [aligned_plot_times, unaligned_plot_times] +#x_s, x_e = int(np.where(x_axis == -1000)[0]), int(np.where(x_axis == 1500)[0]) +labels = ['aligned', 'unaligned'] +colors = ['r', 'k'] + +# Plot laser off trials +fig, ax = plt.subplots() +for i in range(2): + r_squared = r_spearman_laser_off[:, i, :] ** 2 + ax.plot(np.mean(r_squared[:, plot_times[i]], axis = 0), + color = colors[i], label = labels[i]) +ax.legend(loc = 'best', facecolor = 'grey') +ax.set_facecolor('lightgrey') +ax.set( + xlabel = f'Time from transitions or {shift_time}ms after taste delivery', + xticks = [0, 20, 40, 60, 80, 100], + xticklabels = [-1000, -500, 0, 500, 1000, 1500], + ylabel = 'Average Spearman $rho^2$', + title = 'Laser OFF Trials') +plt.savefig('Spearman_laser_OFF.png', bbox_inches="tight", facecolor=fig.get_facecolor(), edgecolor='none') +plt.close('all') + +# Plot laser on trials +fig, ax = plt.subplots() +for i in range(2): + r_squared = r_spearman_laser_on[:, i, :] ** 2 + ax.plot(np.mean(r_squared[:, plot_times[i]], axis = 0), + color = colors[i], label = labels[i]) +ax.legend(loc = 'best', facecolor = 'grey') +ax.set_facecolor('lightgrey') +ax.set( + xlabel = f'Time from transitions or {shift_time}ms after taste delivery', + xticks = [0, 20, 40, 60, 80, 100], + xticklabels = [-1000, -500, 0, 500, 1000, 1500], + ylabel = 'Average Spearman $rho^2$', + title = 'Laser ON Trials') +plt.savefig('Spearman_laser_ON.png', bbox_inches="tight", facecolor=fig.get_facecolor(), edgecolor='none') +plt.close('all') + +# Plot switch times for each taste +# Laser Off trials +fig, ax = plt.subplots() +means = [np.mean(switchtimes_laser_off[taste]) for taste in range(num_tastes)] +errors = [np.std(switchtimes_laser_off[taste])/np.sqrt(len(switchtimes_laser_off[taste]))\ + for taste in range(num_tastes)] +ax.barh(np.arange(4) + 0.5, means, xerr = errors, color = 'grey', height=0.5) +ax.set_yticks(np.arange(4) + 0.5) +ax.set_yticklabels(["NaCl", "Sucrose", "Citric Acid", "QHCl"], rotation=40) +ax.set_xlabel("Time from Taste Delivery (ms)") +ax.set_xticks([0, 250, 500, 750, 1000, 1250])#, [0.0, 0.5, 1.0]) +plt.savefig('Switchtimes_laser_OFF.png', bbox_inches="tight", facecolor=fig.get_facecolor(), edgecolor='none') +plt.close('all') + +# Laser On Trials +fig, ax = plt.subplots() +means = [np.mean(switchtimes_laser_on[taste]) for taste in range(num_tastes)] +errors = [np.std(switchtimes_laser_on[taste])/np.sqrt(len(switchtimes_laser_on[taste]))\ + for taste in range(num_tastes)] +ax.barh(np.arange(4) + 0.5, means, xerr = errors, color = 'grey', height=0.5) +ax.set_yticks(np.arange(4) + 0.5) +ax.set_yticklabels(["NaCl", "Sucrose", "Citric Acid", "QHCl"], rotation=40) +ax.set_xlabel("Time from Taste Delivery (ms)") +ax.set_xticks([0, 250, 500, 750, 1000, 1250])#, [0.0, 0.5, 1.0]) +plt.savefig('Switchtimes_laser_ON.png', bbox_inches="tight", facecolor=fig.get_facecolor(), edgecolor='none') +plt.close('all') \ No newline at end of file