<jupyter_start><jupyter_text>Neuromatch Academy Project - Steinmetz Dataset - Mouse Mind Reading<jupyter_code>session_folder = all_session_folders[0]
session_name =  all_session_names[0]
#session = loader.load_session(session_folder, fast=True)

# load default session = 
session = loader.load_session(all_session_folders[0])
#channels_df = session['channels_df']
spikes_df = session['spikes_df']
clusters_df = session['clusters_df']
trials_df = session['trials_df']
session_name<jupyter_output><empty_output><jupyter_text>Trials:<jupyter_code># select the included trials
trials_df2 = trials_df.loc[trials_df['included']==True]
#trials_select_df = trials_df.loc[trials_df['included']==True, ['start time','end time']] 
##just for the start and end times in the new df

no_contrast_trials_df = trials_df2.loc[(trials_df2['stim contrast left']==0) & (trials_df2['stim contrast right']==0)]


# response choice
'''
response_choice.npy [enumerated type] (nTrials) 
The response registered at the end of the trial, which determines the feedback according to the contrast condition. 
Note that in a small percentage of cases (~4%, see manuscript Methods) 
the initial wheel turn was in the opposite direction. 
-1 for Right choice (i.e. correct when stimuli are on the right); 
+1 for left choice; 
0 for Nogo choice.
'''
# create data frame only with the correct rejection
correctrejection_df = no_contrast_trials_df.loc[(no_contrast_trials_df['response choice']==0)]
''' some of the trials give feedback = -1 ... ??????? we can join the feedback = 1'''
#correctrejection_df = no_contrast_trials_df.loc[(no_contrast_trials_df['response choice']==0) 
#                                                & (no_contrast_trials_df['feedback type']==1)]



# create data frame only with the false alarms
'''same happens here'''
falsealarms_df = no_contrast_trials_df.loc[(no_contrast_trials_df['response choice']==1) 
                                           | (no_contrast_trials_df['response choice']==-1)]
#falsealarms_df = no_contrast_trials_df.loc[((no_contrast_trials_df['response choice']==1) 
#                                           | (no_contrast_trials_df['response choice']==-1))
#                                           & (no_contrast_trials_df['feedback type']==-1)]<jupyter_output><empty_output><jupyter_text>Clusters:<jupyter_code># Drop the rows with clusters labeled as "bad"
clusters_df.drop(clusters_df[clusters_df['label']=='bad'].index, axis=0, inplace=True)
clusters_df
# Filter to visual area
cluster_VISp_df = clusters_df.loc[(clusters_df['recording area']=='VISp')]<jupyter_output><empty_output><jupyter_text>Spikes:<jupyter_code>spikes_df.head()
# Filter to visual areas, clusters = 220, 730, 867, 873, 884 #this for session2? 
spikes_VISp_df = spikes_df.loc[(spikes_df['cluster']==220) | (spikes_df['cluster']==730) 
                               | (spikes_df['cluster']==867) | (spikes_df['cluster']==873)
                               | (spikes_df['cluster']==884)]
spikes_VISp_df.head()<jupyter_output><empty_output><jupyter_text>Plot spike trains:  1) Plotting functions from Max: https://github.com/KceniaB/NMA_project/blob/master/README.md<jupyter_code># 1) Plotting functions from Max: 
# https://github.com/KceniaB/NMA_project/blob/master/README.md 

data = [spikes_df, clusters_df, trials_df]
params = dict()
params['sampling_rate'] = 20000

session_eda = eda(data, all_session_names[0], all_session_folders[0], params)
print((clusters_df['spikes'][3]))
# 1) Plotting functions from Max: 
# https://github.com/KceniaB/NMA_project/blob/master/README.md 
cluster = 4
spikes = clusters_df.loc[cluster, 'spikes']
trials_select_df = trials_df.loc[trials_df['included']==True, ['start time','end time']]
fig, ax = session_eda.plt_spike_train(cluster, spikes, trials_select_df)
plt.show()

'''
This means that Cluster 4 had the following spike times(xx) per trial(yy)
'''
# 1) Plot for many clusters

cluster = [220, 730, 867, 873, 884]

for i in range(len(cluster)): 
    spikes = clusters_df.loc[cluster[i], 'spikes']
    trials_select_df = trials_df.loc[trials_df['included']==True, ['start time','end time']]
    fig, ax = session_eda.plt_spike_train(cluster[i], spikes, trials_select_df)
plt.show()<jupyter_output><empty_output><jupyter_text>Plot spike trains around event with histogram:  2) Plotting functions from Max: https://github.com/KceniaB/NMA_project/blob/master/README.md<jupyter_code># 2) Plotting functions from Max: 
# https://github.com/KceniaB/NMA_project/blob/master/README.md 
title = None
cluster = 8
event = 'stim time'
window = 1
trials_select_df =  session['trials_df'].loc[session['trials_df']['included']==True, 'stim time']
spikes = clusters_df.loc[8, 'spikes']

fig, (ax1, ax2, ax3) = session_eda.plt_spike_train_hist_bar(cluster, event, window, trials_select_df, spikes)
# 2) Plotting functions from Max: 
# https://github.com/KceniaB/NMA_project/blob/master/README.md 
title = None
#cluster = [220, 730, 867, 873, 884]
cluster = [663, 664, 593, 976] #visual primary cortx
#cluster = [264, 271, 273] #MOs

event = 'start time' #select the event to be time = 0 and plot the spikes around it
window = 3 #from -1 to 1 seconds if window = 1

for i in range(len(cluster)): 
    trials_select_df =  session['trials_df'].loc[session['trials_df']['included']==True, 'start time']
    spikes = clusters_df.loc[cluster[i], 'spikes']
    fig, (ax1, ax2, ax3) = session_eda.plt_spike_train_hist_bar(cluster, event, window, trials_select_df, spikes)
plt.show()<jupyter_output><empty_output><jupyter_text>Clusters:<jupyter_code>clusters_df.head()
clusters_df.drop(clusters_df[clusters_df['label']=='bad'].index, axis=0, inplace=True)
clusters_df
# Filter to visual area
cluster_VISp_df = clusters_df.loc[(clusters_df['recording area']=='VISp')]
# cluster_VISp_df = cluster_VISp_df.loc[(cluster_VISp_df['label']=='good')] #already did it above
cluster_VISp_df
# Filter to MOs
cluster_MOs_df = clusters_df.loc[(clusters_df['recording area']=='MOs')]
# cluster_VISp_df = cluster_VISp_df.loc[(cluster_VISp_df['label']=='good')] #already did it above
cluster_MOs_df
# Filter to 
cluster_DG_df = clusters_df.loc[(clusters_df['recording area']=='DG')]
# cluster_VISp_df = cluster_VISp_df.loc[(cluster_VISp_df['label']=='good')] #already did it above
cluster_DG_df<jupyter_output><empty_output><jupyter_text>Trials:<jupyter_code>trials_df
max(trials_df['wheel movement'][8])
# 1. All correct 
# create new variable with true where all the trials are correct = 1
#all_correct = (trials_df['feedback type']==1)
#print(all_correct.head())

# https://stackoverflow.com/questions/17071871/how-to-select-rows-from-a-dataframe-based-on-column-values

# Create data frame only with the correct
correct_df = trials_df.loc[(trials_df['feedback type']==1)]
correct_df.head()

#incorrect_df = trials_df.loc[(trials_df['feedback type']==-1)]
#incorrect_df

no_contrasts2_df = incorrect_df.loc[(incorrect_df['stim contrast left']==0) & (incorrect_df['stim contrast right']==0)]
no_contrasts2_df.head()
# 2. Create one df only for dif contrast trials and another for 0 and 0 

#2.1 diff contract trials
diff_contrasts_df = correct_df.loc[(correct_df['stim contrast left']!=(correct_df['stim contrast right']))]
diff_contrasts_df.head()
# 2. Create one df only for dif contrast trials and another for 0 and 0 

#2.2 same contrast 0
no_contrasts_df = correct_df.loc[(correct_df['stim contrast left']==0) & (correct_df['stim contrast right']==0)]
no_contrasts_df.head()
all_t = len(trials_df)
crrct = len(correct_df)
cont = len(diff_contrasts_df)
nocont = len(no_contrasts_df)

df1 = ''
df1 = {#'session name': [session_name],
      'sessio name': ' CHANGE ',
      'all trials': [all_t], 
      'correct trials': [crrct], 
      'diff contrast t': [cont], 
      'no contras': [nocont], 
      '%Correct': [round(crrct/all_t*100,2)], '%Contrast':[round(cont/crrct*100,2)], '% 0contrast':[round(nocont/crrct*100,2)]}

df2 = pd.DataFrame(data=df1)


print('Correct trials: ', crrct,'/',all_t, '->', round(crrct/all_t*100,2), '%')
print('Dif contrast trials: ', cont,'/',crrct, '->', round(cont/crrct*100,2), '%')
print('0 Contrast trials: ', nocont,'/',crrct, '->', round(nocont/crrct*100,2), '%')
df2
#add the different values per session 
# to initialize
#df3 = df2
#df4 = df2
#to keep adding
df4 = df4.append(df2) 

df4 
test = []

test = trials_df['stim time'] - trials_df['start time']
test
print(max(test))
print(min (test))
#add new column to existing df

# creating a new column TRIAL_NUMBER, from the index
trials_df['trial_number'] = trials_df.index+1
trials_df.head()
df_aligned = trials_df
df_aligned['startA'] = ""
df_aligned['stimA'] = ""
df_aligned['responseA'] = ""
df_aligned['feedbackA'] = ""
df_aligned['endA'] = ""

for i in range(0,len(df_aligned)):
    df_aligned['startA'][i] = df_aligned['start time'][i]-trials_df['start time'][i]
    #print(df_aligned['go cueA'][i])
    df_aligned['stimA'][i] = df_aligned['stim time'][i]-trials_df['start time'][i]
    #print(df_aligned['stim timeA'][i])
    df_aligned['responseA'][i] = df_aligned['response time'][i]-trials_df['start time'][i]
    df_aligned['feedbackA'][i] = df_aligned['feedback time'][i]-trials_df['start time'][i]
    df_aligned['endA'][i] = df_aligned['end time'][i]-trials_df['start time'][i]


    
    
    
    
    
    
#to test 
#df = pd.DataFrame(data) 
#df['C'] = ""
#for i in range(0,len(A)): 
#    df['C'][i] = B[i]-A[i]
#df
min(df_aligned['stimA'])
#scatter plot divided by colors according to visual stim, auditory cue, response (end of rotation), rew/norew, end of the trial
df = df_aligned

ax = plt.gca()
    
for i in range (1,len(df)):
    #df.plot.scatter(x='go cue', y='trial_number', s=0.2, c=np.where((df['response choice']==1), 'r', 'g'), marker='.', grid=False, ax=ax)
    df.plot.scatter(x='startA', y='trial_number', s=0.2, c='b', grid=False, ax=ax)
    df.plot.scatter(x='stimA', y='trial_number', s=0.2, c='c',grid=False, ax=ax)
    df.plot.scatter(x='responseA', y='trial_number', s=0.2, c='y', grid=False, ax=ax)
    df.plot.scatter(x='feedbackA', y='trial_number', s=0.2, c='r', grid=False, ax=ax)
    df.plot.scatter(x='endA', y='trial_number', s=0.2, c='g', grid=False, ax=ax)



#plt.savefig(...)
#fig=plt.figure(figsize=(50, 50))
plt.show()

#how to pick a new variable from dataset "session"
eye_area = session['eye_area']
eye_area
len(eye_area)
wheel_position = session['wheel_position']
wheel_position
plt.plot(wheel_position)
#plt.plot(trials_df['stim contrast left'])
#plt.plot(trials_df['stim contrast right'])
#plt.plot(trials_df['response choice'])

test = np.cumsum(trials_df['wheel movement'])
plt.plot(trials_df['wheel movement'][1])


#scatter plot divided by colors according to 
df = trials_df
df.plot.scatter(x='feedback time', y='response choice', s=20, c=np.where(df['response choice'], 'r', 'g'), grid=False)
df.plot.scatter(x='feedback time', y='response choice', s=20, c=np.where(df['response choice'], 'r', 'g'), grid=False)
plt.show()
go_cues_diff = ((trials_df['stim time'] - trials_df['go cue']))
print(go_cues_diff)
print(min(go_cues_diff)) #0.49
print(max(go_cues_diff)) #1.19 
plt.plot(go_cues_diff)
plt.scatter(trial_df.trials_intervals, trial_df.trials_vStim_cLeft)
# create new variable with true where all the contrast left is 0
is_0 =  trial_df['trials_vStim_cLeft'] ==0


# create new variable with true where all the contrast is 0 in both - NoGo trials
is_0_RL = (trial_df['trials_vStim_cLeft']==0) & (trial_df['trials_vStim_cRight']==0)
print(is_0_RL.head())

len(is_0_RL)
# filter data by creating a new dataframe

# difference between consecutive rows - trials_intervals
# https://pythontic.com/pandas/dataframe-computations/difference

trials_intervals_diff = trial_df.trials_intervals.diff()
plt.plot(trials_intervals_diff)
# https://stackoverflow.com/questions/17071871/how-to-select-rows-from-a-dataframe-based-on-column-values

# Create data frame only with the NoGo
trial2_df = trial_df.loc[(trial_df['trials_vStim_cLeft'] == 0) & (trial_df['trials_vStim_cRight']==0)]
trial2_df.head()
#plt.plot(trial_df.trials_feedbackType)
trial_df.trials_feedbackType

df4 = trial_df['trials_feedbackType']({'a': trial_df['trials_feedbackType'] == 1, 'b': trial_df['trials_feedbackType'] == -1})
plt.figure();
ax = df4.plot.area()

#https://matplotlib.org/3.3.0/tutorials/introductory/pyplot.html

plt.figure(figsize=(9, 3))

plt.subplot(131)
plt.scatter(trial2_df['trials_vStim_times'],trial2_df['trials_feedbackType'])

plt.subplot(132)
plt.bar(trial2_df['trials_feedbackType'], trial2_df['trials_vStim_times'],)

plt.suptitle('NoGo Trials and corresponding decisions')
plt.show()

#plt.scatter(trial2_df['trials_vStim_times'],trial2_df['trials_feedbackType'])

#plotting, within 0 contrast (NoGo trials), how many times they did false positives, i.e., turned the wheel<jupyter_output><empty_output><jupyter_text>Trial history<jupyter_code>#Incomplete
for i in range(3, len(trials_select_df)): 
    if trials_select_df['feedback type'][i] == -1: 
        
        trials_select_df['feedback type'][i]<jupyter_output><empty_output><jupyter_text>Wheel<jupyter_code># add wheelmovement
print(session['wheel_position'].shape)
print(session['wheel_timestamps'].shape)
print(session['wheelMoves_intervals'].shape)
print(session['wheelMoves_type'].shape)

'''
0 for 'flinches' or otherwise unclassified movements, 
1 for left/clockwise turns, 
2 for right/counter-clockwise turns 
(where again "left" means "would be the correct direction for a stimulus presented on the left). 
A detected movement is counted as 'left' or 'right' only if it was sufficient amplitude that it would have registered 
a correct response (and possibly did), within a minimum amount of time from the start of the movement. 
Movements failing those criteria are flinch/unclassified type.
'''
session['wheel_position'][1]
wheel_df = pd.merge( pd.DataFrame(session['wheelMoves_intervals']), pd.DataFrame(session['wheelMoves_type']), how='left', left_index=True, right_index=True)
wheel_df
## Duration of the wheel movement - wheel_mov_dur

#duration of the wheel movement
wheel_mov_dur = wheel_df[1]-wheel_df['0_x'] 
wheel_mov_dur
wheelMoves_type1 = session['wheelMoves_type']


#create empty array of given length; bar plot
##create var wheel_Moves_type2 where 0 is L and 1 is R
wheelMoves_type2 = [0] * (len(wheelMoves_type1))
for i in range(len(wheelMoves_type1)):
    if wheelMoves_type1[i] == 0: 
        wheelMoves_type2[i] = 'unclassified'
        #print('0')
    elif wheelMoves_type1[i] == 1:
        wheelMoves_type2[i] = 'left'
        #print('1')
    else:
        wheelMoves_type2[i] = 'right'
        #print('2')    

        
colors = {'unclassified':'red', 'left':'blue', 'right':'green'}

#plot the figure 
plt.figure(figsize=(9, 3))

plt.subplot(131)
#plt.plot(wheel_mov_time2)

plt.subplot(132)
plt.bar(wheelMoves_type2,wheel_mov_dur)

plt.subplot(133)

#plt.plot(wheel_mov_time2)

plt.suptitle('Title')
plt.show()


#for i in range(len(test1)):
#    if test1[i] == 0: 
#        print('0')
#    else: 
#        print('1')
len(trials_df['wheel movement'][1])

for i in range(len(trials_df['wheel movement'][1])): 

    plt.plot(trials_df['wheel movement'][1][i], 'ro',markersize=1)<jupyter_output><empty_output><jupyter_text>testing other codes sample_dataset_downloadthings I edited/doubts start with  KB<jupyter_code>## Data retrieval
import os, requests

fname = []
for j in range(3):
    fname.append('steinmetz_part%d.npz'%j)
url = ["https://osf.io/agvxh/download"]
url.append("https://osf.io/uv3mw/download")
url.append("https://osf.io/ehmw2/download")

for j in range(len(url)):
    if not os.path.isfile(fname[j]):
        try:
            r = requests.get(url[j])
        except requests.ConnectionError:
            print("!!! Failed to download data !!!")
        else:
            if r.status_code != requests.codes.ok:
                print("!!! Failed to download data !!!")
            else:
                with open(fname[j], "wb") as fid:
                    fid.write(r.content)<jupyter_output><empty_output><jupyter_text>this script is based on the W1D4 tutorial 2 logistic regression decodersimple_decoder<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import os
import requests
import hashlib

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
#Helper functions taken from the tutorial

def plot_weights(models, sharey=True):
    """Draw a stem plot of weights for each model in models dict."""
    n = len(models)
    f = plt.figure(figsize=(10, 2.5 * n))
    axs = f.subplots(n, sharex=True, sharey=sharey)
    axs = np.atleast_1d(axs)
    
    for ax, (title, model) in zip(axs, models.items()):
        ax.margins(x=.02)
        stem = ax.stem(model.coef_.squeeze(), use_line_collection=True)
        stem[0].set_marker(".")
        stem[0].set_color(".2")
        stem[1].set_linewidths(.5)
        stem[1].set_color(".2")
        stem[2].set_visible(False)
        ax.axhline(0, color="C3", lw=3)
        ax.set(ylabel="Weight", title=title)
    ax.set(xlabel="Neuron (a.k.a. feature)")
    f.tight_layout()


def plot_function(f, name, var, points=(-10, 10)):
    """Evaluate f() on linear space between points and plot.

    Args:
      f (callable): function that maps scalar -> scalar
      name (string): Function name for axis labels
      var (string): Variable name for axis labels.
      points (tuple): Args for np.linspace to create eval grid.
    """
    x = np.linspace(*points)
    ax = plt.figure().subplots()
    ax.plot(x, f(x))
    ax.set(
      xlabel=f'${var}$',
      ylabel=f'${name}({var})$'
    )

def plot_model_selection(C_values, accuracies):
    """Plot the accuracy curve over log-spaced C values."""
    ax = plt.figure().subplots()
    ax.set_xscale("log")
    ax.plot(C_values, accuracies, marker="o")
    best_C = C_values[np.argmax(accuracies)]
    ax.set(
      xticks=C_values,
      xlabel="$C$",
      ylabel="Cross-validated accuracy",
      title=f"Best C: {best_C:1g} ({np.max(accuracies):.2%})",
    )

def plot_non_zero_coefs(C_values, non_zero_l1, n_voxels):
    """Plot the accuracy curve over log-spaced C values."""
    ax = plt.figure().subplots()
    ax.set_xscale("log")
    ax.plot(C_values, non_zero_l1, marker="o")
    ax.set(
    xticks=C_values,
    xlabel="$C$",
    ylabel="Number of non-zero coefficients",
    )
    ax.axhline(n_voxels, color=".1", linestyle=":")
    ax.annotate("Total\n# Neurons", (C_values[0], n_voxels * .98), va="top")

def sigmoid(z):
    return 1 / (1+np.exp(-z))<jupyter_output><empty_output><jupyter_text>download the steinmetz dataset used in the tutorial<jupyter_code>## from the tutorial, from the data retrieval and loading cell
url = "https://osf.io/r9gh8/download"
fname = "W1D4_steinmetz_data.npz"
expected_md5 = "d19716354fed0981267456b80db07ea8"

if not os.path.isfile(fname):
    try:
        r = requests.get(url)
    except requests.ConnectionError:
        print("!!! Failed to download data !!!")
    else:
        if r.status_code != requests.codes.ok:
            print("!!! Failed to download data !!!")
        elif hashlib.md5(r.content).hexdigest() != expected_md5:
            print("!!! Data download appears corrupted !!!")
        else:
            with open(fname, "wb") as fid:
                fid.write(r.content)
## from the tutorial, from the data retrieval and loading cell 
def load_steinmetz_data(data_fname=fname):
    with np.load(data_fname) as dobj:
        data = dict(**dobj)
    return data<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code>#this is sample data provided in tutorial 
## KB data with array of spikes (normalized spike rates!, n_trials, n_neurons) and choices (0s and 1s, n_trials)
data = load_steinmetz_data() 
for key, val in data.items():
    print(key, val.shape)
len((data['spikes'][0]))<jupyter_output><empty_output><jupyter_text>* an X matrix with shape (n_samples, n_features) - something about the experiment * a y vector with length n_samples. - neural data<jupyter_code>#DATA #here we need to load our own data!!
y = data["choices"] #these are the binary choices the animal made on each trial (trialnr)
X = data["spikes"] #these are the mean spike rates for each neuron (trialnr x neurons)
plt.figure
plt.plot(X[0,:]) ## KB spikes of the first trial<jupyter_output><empty_output><jupyter_text>fit logistic regression models<jupyter_code># First define the model
log_reg = LogisticRegression(penalty='none',solver='newton-cg')

#Then fit it to data
log_reg.fit(X, y)<jupyter_output><empty_output><jupyter_text>There's two steps here:- We *initialized* the model with a hyperparameter, telling it what penalty to use (we'll focus on this in the second part of the notebook)- We *fit* the model by passing it the `X` and `y` objects.<jupyter_code>#get accuracy of the model on thectraining set (should be high!)

# Fitting the model performs maximum likelihood optimization, learning a set of *feature weights*. 
# We can use those learned weights to *classify* new data, or predict the labels for each sample:
y_pred = log_reg.predict(X) ## classifying the training data
len(y_pred)
def compute_accuracy(X, y, model):
  """Compute accuracy of classifier predictions.
  
  Args:
    X (2D array): Data matrix
    y (1D array): Label vector
    model (sklearn estimator): Classifier with trained weights.

  Returns:
    accuracy (float): Proportion of correct predictions.  
  """

  y_pred = model.predict(X)
  accuracy = np.sum(y_pred == y)/len(y) #or np.mean() in the last part

  return accuracy

# Uncomment and run to test your function:
train_accuracy = compute_accuracy(X, y, log_reg)
print(f"Accuracy on the training data: {train_accuracy:.2%}")
## Sam (same as prev step)
#train_accuracy = np.mean(y_pred == y)
#print(f"Accuracy on the training data: {train_accuracy:.2%}")<jupyter_output><empty_output><jupyter_text>Cross-validatinng using scikit-learn helper functions L2 Rwgularization below<jupyter_code>#get accuracy of the model on a k-fold cross validated set

k=8 #choose k-fold
reg = 'l2' #choose regularization. 'none', 'l1', or 'l2'
lam = 1 #choose lambda for regularizer

#fit model
#log_reg = LogisticRegression(penalty='none').fit(X, y)
log_reg = LogisticRegression(penalty=reg, C=1/lam, solver='lbfgs').fit(X, y)

#get cross-validated accuracies
#accuracies = cross_val_score(LogisticRegression(penalty='none'), X, y, cv=8) ##KB the way it was used in the tutorial 
accuracies = cross_val_score(log_reg, X, y, cv=k)

#plot the resulting accuracies
f, ax = plt.subplots(figsize=(8, 3))
ax.boxplot(accuracies, vert=False, widths=.7)
ax.scatter(accuracies, np.ones_like(accuracies))
ax.set(
  xlabel="Accuracy",
  yticks=[],
  title=f"Average test accuracy: {accuracies.mean():.2%}"
)
ax.spines["left"].set_visible(False)
## KB 
log_reg = LogisticRegression(penalty="none").fit(X, y)
plot_weights({"No regularization": log_reg})
##KB 
log_reg_l2 = LogisticRegression(penalty="l2", C=1).fit(X, y)

# now show the two models
models = {
  "No regularization": log_reg,
  "$L_2$ (C = 1)": log_reg_l2,
}
plot_weights(models)
## KB just to test L1

log_reg_l1 = LogisticRegression(penalty="l1", C=1, solver="saga", max_iter=5000) #too low can we stuck in local min
log_reg_l1.fit(X, y)
models = {
  "$L_2$ (C = 1)": log_reg_l2,
  "$L_1$ (C = 1)": log_reg_l1,
}
plot_weights(models)

#S almost uses L2 
#once fitted we can take out the 0s, it can be faster
def model_selection(X, y, C_values, CV_values):
    """Compute CV accuracy for each C value.

    Args:
    X (2D array): Data matrix
    y (1D array): Label vector
    C_values (1D array): Array of hyperparameter values

    Returns:
    accuracies (1D array): CV accuracy with each value of C

    """
    accuracies_all = []

    for CV in CV_values:

        accuracies = []

        for C in C_values:

            # Initialize and fit the model
            model = LogisticRegression(penalty='l2',C=C, solver='lbfgs',max_iter=5000)
            model.fit(X,y)

            # Get the accuracy for each test split
            accs = cross_val_score(model,X,y,cv=int(CV))

            # Store the average test accuracy for this value of C
            accuracies.append(np.mean(accs))

        accuracies_all.append(accuracies)

    return accuracies_all



# Use log-spaced values for C and CV

C_values = np.logspace(-4, 4, 9)

#CV_values = np.linspace(2,20,10)
CV_values = np.logspace(1,7,7,base=2)
#CV_values = np.logspace(np.arange(2,50,2)

accuracies_all = model_selection(X, y, C_values, CV_values)
ax = plt.figure().subplots()
ax.set_xscale("log")
cmap = plt.cm.get_cmap('coolwarm')
for i,accuracies in enumerate(accuracies_all):
    c = cmap(i/(len(accuracies_all)-1))
    ax.plot(C_values, accuracies, marker="o",color=c)
    best_C = C_values[np.argmax(accuracies)]
    plt.text(1,0.8-0.02*i,f"{CV_values[i]:1g}-fold CV: Cmax={best_C:1g}",color=c)
ax.set(
  xticks=C_values,
  xlabel="$C$",
  ylabel="Cross-validated accuracy",
)
plt.show()<jupyter_output><empty_output><jupyter_text>Old stuff or useful code/links<jupyter_code>##To adjust the size of one figure:

#import matplotlib.pyplot as plt

#fig=plt.figure(figsize=(15, 15))
#To change the default settings, and therefore all your plots:

#import matplotlib.pyplot as plt

#plt.rcParams['figure.figsize'] = [15, 15]<jupyter_output><empty_output><jupyter_text>Useful links: data: https://www.w3resource.com/pandas/dataframe/dataframe-plot-area.phphttps://stackoverflow.com/questions/17071871/how-to-select-rows-from-a-dataframe-based-on-column-valueshttps://www.delftstack.com/howto/python-pandas/how-to-create-an-empty-column-in-pandas-dataframe/loop issues: https://www.datacamp.com/community/tutorials/for-loops-in-python?utm_source=adwords_ppc&utm_campaignid=898687156&utm_adgroupid=48947256715&utm_device=c&utm_keyword=&utm_matchtype=b&utm_network=g&utm_adpostion=&utm_creative=332602034349&utm_targetid=dsa-429603003980&utm_loc_interest_ms=&utm_loc_physical_ms=1011752&gclid=EAIaIQobChMIkfT0pKXj6gIVhIbVCh0ctg9-EAAYASAAEgKCSvD_BwEplot: https://matplotlib.org/3.3.0/tutorials/introductory/pyplot.htmlhttps://pandas.pydata.org/pandas-docs/stable/user_guide/visualization.html<jupyter_code>## JUST TO CHECK THE WHAT THE DIFF VARS MEAN
# pars trials info =================================
# create spike_df tataframe, with each spike time and the cluster it belongs to
trial_df = pd.DataFrame( 
    { 'trials_resp_choice':session['trials_response_choice'][:,0], #-1 for Right choice (i.e. correct when stimuli are on the right); +1 for left choice; 0 for Nogo choice. Note that in a small percentage of cases (~4%, see manuscript Methods) the initial wheel turn was in the opposite direction. 
    'trials_resp_times': session['trials_response_times'][:,0], 
     
    'trials_goCue_times': session['trials_goCue_times'][:,0] ,#auditory cue
     
    'trials_vStim_cLeft': session['trials_visualStim_contrastLeft'][:,0], #A value of 0.5 means 50% contrast. 0 is a blank screen: no change to any pixel values on that side (completely undetectable).
    'trials_vStim_cRight': session['trials_visualStim_contrastRight'][:,0], 
    'trials_vStim_times': session['trials_visualStim_times'][:,0], 

    'trials_intervals': session['trials_intervals'][:,0], 
                                                
    'trials_feedback_times': session['trials_feedback_times'][:,0],
    'trials_feedbackType': session['trials_feedbackType'][:,0], # -1 negative - noise; 1 positive - water
    
    'trials_included': session['trials_included'][:,0], #Importantly, while this variable gives inclusion criteria according to the definition of disengagement (see manuscript Methods), it does not give inclusion criteria based on the time of response, as used for most analyses in the paper.
    'trials_repNum': session['trials_repNum'][:,0], #Trials are repeated if they are "easy" trials (high contrast stimuli with large difference between the two sides, or the blank screen condition) and this keeps track of how many times the current trial's condition has been repeated.
    },
                       )
#????? 
# add spikes to session dictionary
session['trial_df']=trial_df<jupyter_output><empty_output><jupyter_text>Table of Contents1&nbsp;&nbsp;Neuromatch Academy Project - Steinmetz Dataset - Mouse Mind Reading1.0.0.1&nbsp;&nbsp;explore dictionary for session1.1&nbsp;&nbsp;Load Data1.1.1&nbsp;&nbsp;via oneibl by client1.1.2&nbsp;&nbsp;via google drive directly1.1.3&nbsp;&nbsp;cleanup channel infos1.1.4&nbsp;&nbsp;cleaning up spikes1.1.5&nbsp;&nbsp;cleanup trial info<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import os
import pandas as pd
import matplotlib.pyplot as plt

from nma_class_k import NMA_project
#from nma_class_Copy1 import NMA_project

%matplotlib inline
# only if you are on jupyter notebook not if on colab
#import qgrid 
# link to parent folder with all datasets
#main_folder = "/content/gdrive/My Drive/3 Projekte/Neuromatch Summer School/Projekt/steinmetz_dataset"
#main_folder = "C:/Users/User/Google Drive/3 Projekte/Neuromatch Summer School/Projekt/steinmetz_dataset"
#main_folder = "Z:/3 Projekte/Neuromatch Summer School/Projekt/steinmetz_dataset"
main_folder = '/Volumes/GoogleDrive/My Drive/NMA_Dataset/Dataset_1'

# create session object
nma_obj = NMA_project(main_folder) 

# default session is initialized, extract interesting dataframes in dedicated variables
#session = nma_obj.std_session
#channels_df = session['channels_df']
#spikes_df = session['spikes_df']
#clusters_df = session['clusters_df']
#trials_df = session['trials_df']





channels_df = session['channels_df']
spikes_df = session['spikes_df']
clusters_df = session['clusters_df']
trials_df = session['trials_df']
spikes_ar = session['spikes_ar']<jupyter_output><empty_output><jupyter_text>explore dictionary for session- each ession is a dictionary- each element kontains eather a numpy array or a dataframe- the most important dataframes we are going to work with are:    - clusters_df: all the spikes and clusters    - trials_df: all the behaviour data for all trials<jupyter_code>#how to pick a new variable from dataset "session"
some_var = session['spontaneous_intervals']
#print(len(some_var))
#some_var
for key in nma_obj.std_session.keys():
    print(key)
#test
A = [1, 2, 3,4,5,6,7]
for i in range(3,len(A)):
    print(i, A[i])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13969453112169486, 'fraction_numerical': 0.020003284611594677, 'mean_word_length': 4.9639569049951024, 'pattern_counts': {'":': 6, '<': 89, '<?xml version=': 0, '>': 93, 'https://': 20, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '86296', 'license': 'MIT', 'n_tokens_mistral': 11313, 'n_tokens_neox': 10741, 'n_words': 2508, 'path': 'kce/checking the sessions - false alarms and neurodata.ipynb', 'repo_name': 'KceniaB/NMA_project', 'size': '30445'}