<jupyter_start><jupyter_text>trying tetrode 10 with KS output.<jupyter_code>animal = 'Li'
task = 'T3g'
date = '062718'
fn = animal+'_'+task+'_'+date+'_KSClusters'
orig_time =np.arange(t[0],t[-1]+1,1/32000)

good_tts = [3,8,10]
datPath = Path('/Volumes/groups/giocomo/alexg/Clustered/',animal,fn)

nTotalCells=0
cell_sp = []
cell_ids = {}
nTimePoints=len(t)
for tt in good_tts:
    ttPath = Path(datPath,'tt_'+str(tt))
    spikes = np.load(str(ttPath/'spike_times.npy'))
    clusters = np.load(str(ttPath/'spike_clusters.npy'))
    cl_rating = pd.read_csv(str(ttPath/'cluster_group.tsv'),sep='\t')
    good_cl_ids = np.where(cl_rating['group']=='good')[0]
    nGoodCells = len(good_cl_ids)
    bin_spikes = np.zeros((nGoodCells,nTimePoints))
    cnt=0
    for cl in good_cl_ids:
        cell_ids[cnt] = [tt,cl]
        sp_per_bin, _=np.histogram(orig_time[spikes[clusters==cl]],bins=nTimePoints)    
        bin_spikes[cnt,:],_=np.histogram(orig_time[spikes[clusters==cl]],bins=nTimePoints)    
        cnt+=1
    nTotalCells += nGoodCells
    cell_sp.append(bin_spikes)
cell_sp=np.concatenate(cell_sp)
nCells = nTotalCells
minMark=43
window=np.where(np.logical_and(t>=minMark*60,t<(minMark+1)*60))[0]
f,(a1,a2)=plt.subplots(2,1, figsize=(12,12))
t = PosDat['t']
a1.set_position([0.125, 0.4, 0.775, 0.4])
a1=plotBehWindow(t[window],PosDat['PosMat'].iloc[window].T,TMF.ZonesNames,a1)
a1.set_xticks([])

a2=plotBehWindow(t[window],cell_sp[:,window],np.arange(nCells).astype(str),a2)
a2.set_xlabel('Time[s]')
a2.set_ylabel('Cell Number')
a2.set_position([0.125, 0.15, 0.775, 0.23])
yt=np.linspace(0,nCells,5, endpoint=True).astype(int)
a2.set_yticks(yt)
a2.set_yticklabels(yt.astype(str))
for tick in a2.get_yticklabels():
    tick.set_rotation(0)

f.savefig('/Users/alexgonzalez/Google Drive/PostDoc/FloorMeeting/LocSpikesWindow.png',dpi=500,
          bbox_inches='tight',pad_inches=0.2)    
def PDeviance(y,mu):
    y = y+1e-10
    mu = mu+1e-10
    return 2*np.sum(y*np.log(y/mu)-(y-mu))
def plotFitComp(orig,fit,spks,time,a):    
    a.plot(time,orig,color='k',alpha=0.7)
    a.plot(time,fit,color='r',alpha=1)
    a.plot(time,spks,'k',alpha=0.5)
    a.set_ylabel('FR[sp/s]')
    a.set_xlabel('Time[s]')
    a.set_ylim([0,a.get_ylim()[1]])
    return a
#selected_params = ['Home','SegA','Center','SegB','I1','SegC','G1','SegD',
#                   'G2','SegE','I2','SegF','G3','SegG','G4','RH','RC','R1','R2','R3','R4','CL','CR']
selected_params = ['SegA','Center','SegB','I1','SegC','G1','SegD',
                   'G2','SegE','I2','SegF','G3','SegG','G4','RH','RC','R1','R2','R3','R4','CL','CR']

pos_params = ['SegA','Center','SegB','I1','SegC','G1','SegD',
              'G2','SegE','I2','SegF','G3','SegG','G4']

reward_params = ['RH','RC','R1','R2','R3','R4']
cue_params = ['CL','CR']
dir_params = TMF.SegDirNames
goal_params = ['Home','Center','G1','G2','G3','G4']

FullMat=pd.concat([PosDat['PosMat'],PosDat['SegDirMat'],EventMat],axis=1)

#selected_params =goal_params+dir_params+reward_params+cue_params
selected_params =pos_params+reward_params+cue_params
#selected_params = pos_params+reward_params+cue_params+dir_params
M2 = sm.add_constant(FullMat[selected_params])

y=cell_sp[0,:]
xTrain, xTest, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)
time_test = t[xTest.index]
window_test=np.where(np.logical_and(time_test>=minMark,time_test<(minMark+60)))[0]
    
glm_train_fits = {}
PR2_Train = np.zeros(nCells)
PR2_Test = np.zeros(nCells)
R2_Test = np.zeros(nCells)
# Turn interactive plotting off
pltshow = 1
if not pltshow:
    plt.ioff()
    
for cell in [9]:#np.arange(nCells):
    
    y=cell_sp[cell,:]
    if cell>0:
        _, _, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)

    glm_train = sm.GLM(yTrain,xTrain,family = sm.families.Poisson())
    glm_null = sm.GLM(yTrain,xTrain['const'],family = sm.families.Poisson())
    glm_null_fit = glm_null.fit()
    glm_train_fits[cell] = glm_train.fit()

    PR2_Train[cell] = 1-glm_train_fits[cell].deviance/glm_null_fit.deviance
    muTest =glm_train_fits[cell].predict(xTest)
    TestDev = PDeviance(yTest,muTest)
    muNullTest = glm_null_fit.predict(xTest['const'])
    TestNullDev = PDeviance(yTest,muNullTest)
    PR2_Test[cell] = 1-TestDev/TestNullDev 
    
    f,a= plot_eFR(glm_train_fits[cell],step,selected_params)
    a.set_title('Cell_ID={0}, Train_PR2={1:0.2f}, Test_PR2={2:0.2f}'.format(cell,PR2_Train[cell],PR2_Test[cell]))
    f.savefig('/Users/alexgonzalez/Documents/Data/Li_T3g_062718_Results/GLM_Plots/GLM_Fits_Cell_{}.pdf'.format(cell),
              bbox_inches='tight',pad_inches=0.2)
    if not pltshow:
        plt.close(f)

    lfwin = np.round(1.0/step).astype(int)
    sm_yTest = signal.filtfilt(np.ones(lfwin)/lfwin,1,yTest)/step
    sm_muTest = signal.filtfilt(np.ones(lfwin)/lfwin,1,muTest.as_matrix())/step

    R2_Test[cell] = np.corrcoef(sm_yTest,sm_muTest)[0,1]**2
    
    
    f,(a1,a2)=plt.subplots(2,1, figsize=(12,10))
    a1.set_position([0.125, 0.38, 0.775, 0.53])
    a1=plotBehWindow(time_test[window_test],xTest.iloc[window_test].T,xTest.columns,a1)
    a1.set_xticks([])    
    a1.set_title('Cell_ID={0}, OoS_R2={1:0.2f}'.format(cell,R2_Test[cell]))

    a2.set_position([0.125, 0.2, 0.775, 0.15])
    a2=plotFitComp(sm_yTest[window_test],sm_muTest[window_test],
                   yTest[window_test]>0,time_test[window_test],a2 )
    a2.set_xlim(time_test[window_test][0],time_test[window_test][-1])
    f.savefig('/Users/alexgonzalez/Documents/Data/Li_T3g_062718_Results/GLM_Plots/GLM_FR_Plot_Cell_{}.png'.format(cell),
              dpi=500,bbox_inches='tight',pad_inches=0.2)  
    if not pltshow:
        plt.close(f)
def plot_eFR2(res,baseRate,step,selected_params):
    ci=(np.exp(res.conf_int())-baseRate)/step
    estimatedFRs = (np.exp(res.params)-baseRate)/step
    nParams=len(res.params)
    x=np.tile(np.arange(nParams),[2,1])
    #sns.barplot(x=np.arange(57),y=np.exp(res.params)-1,ci=np.exp(res.conf_int()-1),capsize=.2)
    
    f,a=plt.subplots(1,1, figsize=(14,5))
    ab=sns.barplot(x=np.arange(nParams),y=estimatedFRs,ax=a, 
                   facecolor=(0.5, 0.5, 0.5, 1))
    for p in np.arange(nParams):
        a.plot(x[:,p],ci.iloc[p],color='k',linewidth=3)
    a.set_ylabel('FR[sp/s]')
    a.set_xticklabels(selected_params)
    #a.grid()
    sns.set_style("whitegrid")
    sns.despine(left=True)
    for tick in a.get_xticklabels():
        tick.set_rotation(45)
    return f,a
def plot_params(params,ci,selected_params):

    nParams=len(params)
    x=np.tile(np.arange(nParams),[2,1])
    
    f,a=plt.subplots(1,1, figsize=(14,5))
    ab=sns.barplot(x=np.arange(nParams),y=params,ax=a, 
                   facecolor=(0.5, 0.5, 0.5, 1))
    for p in np.arange(nParams):
        a.plot(x[:,p],ci.iloc[p],color='k',linewidth=3)
        
    a.set_ylabel('PE [au]')
    a.set_xticklabels(selected_params)
    #a.grid()
    sns.set_style("whitegrid")
    sns.despine(left=True)
    for tick in a.get_xticklabels():
        tick.set_rotation(45)
    return f,a
model = 'Complete'

pos_params = ['Home','SegA','Center','SegB','I1','SegC','G1','SegD',
              'G2','SegE','I2','SegF','G3','SegG','G4']

reward_params = ['RH','RC','R1','R2','R3','R4']
cue_params = ['CL','CR']
dir_params = TMF.SegDirNames
goal_params = ['Home','Center','G1','G2','G3','G4']

sp = '/Users/alexgonzalez/Documents/Data/Li_T3g_062718_Results/GLM_Plots'

PlotDir = Path(sp,model)
PlotDir.mkdir(parents=True, exist_ok=True)

FullMat=pd.concat([PosDat['PosMat'],PosDat['SegDirMat'],EventMat],axis=1)

if model=='Complete':
    selected_params = pos_params + reward_params + cue_params
    M2 = sm.add_constant(FullMat[selected_params])
elif model=='Rewards':
    selected_params = reward_params
    M2 = sm.add_constant(FullMat[reward_params])
elif model=='Cues':
    selected_params = cue_params
    M2 = sm.add_constant(FullMat[reward_params])

y=cell_sp[0,:]
xTrain, xTest, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)
time_test = t[xTest.index]
window_test=np.where(np.logical_and(time_test>=minMark,time_test<(minMark+60)))[0]
    
glm_train_fits = {}
PR2_Train = np.zeros(nCells)
PR2_Test = np.zeros(nCells)
R2_Test = np.zeros(nCells)
TrainBaseRate = np.zeros(nCells)
nParams = len(selected_params)
ZVals = pd.DataFrame(np.zeros((nCells,nParams)),columns=selected_params)
ParamsEst = pd.DataFrame(np.zeros((nCells,nParams)),columns=selected_params)
PVals = pd.DataFrame(np.zeros((nCells,nParams)),columns=selected_params)
plotFlag = 1
pltshow = 0
if not pltshow:
    plt.ioff()
    
for cell in np.arange(nCells):
    
    y=cell_sp[cell,:]
    if cell>0:
        _, _, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)

    TrainBaseRate[cell] = np.mean(yTrain)
    
    # train model and firt to test set
    if model=='Complete':
        glm_train = sm.GLM(yTrain,xTrain.drop(['const'],axis=1),family = sm.families.Poisson())
        glm_train_fits[cell] = glm_train.fit()
        muTest =glm_train_fits[cell].predict(xTest.drop(['const'],axis=1))
    else:
        glm_train = sm.GLM(yTrain,xTrain,family = sm.families.Poisson())
        glm_train_fits[cell] = glm_train.fit()
        muTest =glm_train_fits[cell].predict(xTest)
        
    
    # get null / constant model fits
    glm_null = sm.GLM(yTrain,xTrain['const'],family = sm.families.Poisson())
    glm_null_fit = glm_null.fit()
    muNullTest = glm_null_fit.predict(xTest['const'])
   
    # get performance metrics
    PR2_Train[cell] = 1-glm_train_fits[cell].deviance/glm_null_fit.deviance
    TestDev = PDeviance(yTest,muTest)
    TestNullDev = PDeviance(yTest,muNullTest)
    PR2_Test[cell] = 1-TestDev/TestNullDev 
    
    # Store parameters
    ZVals.iloc[cell] = glm_train_fits[cell].tvalues
    ParamsEst.iloc[cell] = np.exp(glm_train_fits[cell].params)
    PVals.iloc[cell]=glm_train_fits[cell].pvalues

    # smooth binned spikes to obtain FRs
    lfwin = np.round(0.5/step).astype(int)
    sm_yTest = signal.filtfilt(np.ones(lfwin)/lfwin,1,yTest)/step
    sm_muTest = signal.filtfilt(np.ones(lfwin)/lfwin,1,muTest.as_matrix())/step

    R2_Test[cell] = np.corrcoef(sm_yTest,sm_muTest)[0,1]**2
    
    # plot results
    if plotFlag:
        f,a= plot_params(glm_train_fits[cell].params,glm_train_fits[cell].conf_int(),selected_params)
        a.set_title('Cell_ID={0}, Train_PR2={1:0.2f}, Test_PR2={2:0.2f}'.format(cell,PR2_Train[cell],PR2_Test[cell]))
        f.savefig(('{}/{}_GLM_Fits_RawParams_Cell_{}.pdf').format(sp,model,cell),bbox_inches='tight',pad_inches=0.2)
        if not pltshow:
            plt.close(f)
        
        f,a= plot_eFR2(glm_train_fits[cell],TrainBaseRate[cell],step,selected_params)
        a.set_title('Cell_ID={0}, Train_PR2={1:0.2f}, Test_PR2={2:0.2f}'.format(cell,PR2_Train[cell],PR2_Test[cell]))
        f.savefig(('{}/{}_GLM_Fits_eFRs_Cell_{}.pdf').format(sp,model,cell),bbox_inches='tight',pad_inches=0.2)
        if not pltshow:
            plt.close(f)
        
        f,(a1,a2)=plt.subplots(2,1, figsize=(12,10))
        a1.set_position([0.125, 0.38, 0.775, 0.53])
        a1=plotBehWindow(time_test[window_test],xTest.iloc[window_test].T,xTest.columns,a1)
        a1.set_xticks([])    
        a1.set_title('Cell_ID={0}, OoS_R2={1:0.2f}'.format(cell,R2_Test[cell]))

        a2.set_position([0.125, 0.2, 0.775, 0.15])
        a2=plotFitComp(sm_yTest[window_test],sm_muTest[window_test],
                       yTest[window_test]>0,time_test[window_test],a2 )
        a2.set_xlim(time_test[window_test][0],time_test[window_test][-1])
        f.savefig(('{}/{}_GLM_TestWindow_Cell_{}.png').format(sp,model,cell),dpi=500,bbox_inches='tight',pad_inches=0.2)

        if not pltshow:
            plt.close(f)
cell=14
y=cell_sp[cell,:]
_, _, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)
glm_train = sm.GLM(yTrain,xTrain.drop(['const'],axis=1),family = sm.families.Poisson())
glm_null = sm.GLM(yTrain,xTrain['const'],family = sm.families.Poisson())
glm_null_fit = glm_null.fit()
glm_train_fits[cell] = glm_train.fit()

pa='Center'
np.mean(yTrain[xTrain[pa]])/step,ParamsEst[pa].loc[cell]/step

np.exp(glm_train_fits[cell].params['Home'])
#X=(np.ParamsEst[reward_params]-np.tile(TrainBaseRate,[len(reward_params),1]).T>0)*(PVals[reward_params]*nParams*nCells<0.001)

X=np.log(ParamsEst[reward_params])
X
ZVals[reward_params].loc[14]
(ParamsEst[['CL','CR']].as_matrix()>0)*(PVals[['CL','CR']]*nParams*nCells<0.001).as_matrix()
pa = 'Center'
((ParamsEst[pa]-TrainBaseRate)>0)*(PVals[pa]<0.001)*R2_Test>0.3
np.sqrt(R2_Test)[R2_Test>0.5]
200
#selected_params = ['Home','SegA','Center','SegB','I1','SegC','G1','SegD',
#                   'G2','SegE','I2','SegF','G3','SegG','G4','RH','RC','R1','R2','R3','R4','CL','CR']
selected_params = ['SegA','Center','SegB','I1','SegC','G1','SegD',
                   'G2','SegE','I2','SegF','G3','SegG','G4','RH','RC','R1','R2','R3','R4','CL','CR']

pos_params = ['SegA','Center','SegB','I1','SegC','G1','SegD',
              'G2','SegE','I2','SegF','G3','SegG','G4']

reward_params = ['RH','RC','R1','R2','R3','R4']
cue_params = ['CL','CR']
dir_params = TMF.SegDirNames
goal_params = ['Home','Center','G1','G2','G3','G4']

FullMat=pd.concat([PosDat['PosMat'],PosDat['SegDirMat'],EventMat],axis=1)

#selected_params =goal_params+dir_params+reward_params+cue_params
selected_params =pos_params+reward_params+cue_params
#selected_params = pos_params+reward_params+cue_params+dir_params
M2 = sm.add_constant(FullMat[selected_params])

y=cell_sp[0,:]
xTrain, xTest, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)
RewardMat = EventMat[reward_params]

detect_params = ['DH','DC','D1','D2','D3','D4']
DetecMat = EventMat[detect_params] 
reLength = 1/step # one second
deLength = 1.5*1/step # two secs

for r in reward_params:
    RewardMat[r] = (signal.lfilter(np.arange(reLength),1,RewardMat[r])>0).astype(int)
    
for d in detect_params:
    DetecMat[d] = (signal.lfilter(np.arange(deLength),1,DetecMat[d])>0).astype(int)
    
cnt = 0

ReDeMat=pd.concat([DetecMat,RewardMat],axis=1)
selected_params = detect_params+reward_params
M2 = sm.add_constant(ReDeMat)

y=cell_sp[0,:]
xTrain, xTest, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)
time_test = t[xTest.index]
window_test=np.where(np.logical_and(time_test>=minMark,time_test<(minMark+60)))[0]
    
glm_train_fits = {}
PR2_Train = np.zeros(nCells)
PR2_Test = np.zeros(nCells)
R2_Test = np.zeros(nCells)
# Turn interactive plotting off
plt.ioff()
for cell in np.arange(nCells):
    
    if cell>0:
        y=cell_sp[cell,:]
        _, _, yTrain, yTest = train_test_split(M2,y, test_size = 0.5, random_state = 0,shuffle=False)

    glm_train = sm.GLM(yTrain,xTrain,family = sm.families.Poisson())
    glm_null = sm.GLM(yTrain,xTrain['const'],family = sm.families.Poisson())
    glm_null_fit = glm_null.fit()
    glm_train_fits[cell] = glm_train.fit()

    PR2_Train[cell] = 1-glm_train_fits[cell].deviance/glm_null_fit.deviance
    muTest =glm_train_fits[cell].predict(xTest)
    TestDev = PDeviance(yTest,muTest)
    muNullTest = glm_null_fit.predict(xTest['const'])
    TestNullDev = PDeviance(yTest,muNullTest)
    PR2_Test[cell] = 1-TestDev/TestNullDev 
    
    f,a= plot_eFR(glm_train_fits[cell],step,selected_params)
    a.set_title('Cell_ID={0}, Train_PR2={1:0.2f}, Test_PR2={2:0.2f}'.format(cell,PR2_Train[cell],PR2_Test[cell]))
    f.savefig('/Users/alexgonzalez/Documents/Data/Li_T3g_062718_Results/GLM_Plots/GLM_FitsReDe_Cell_{}.pdf'.format(cell),
              bbox_inches='tight',pad_inches=0.2)    
    plt.close(f)

    lfwin = np.round(1.0/step).astype(int)
    sm_yTest = signal.filtfilt(np.ones(lfwin)/lfwin,1,yTest)/step
    sm_muTest = signal.filtfilt(np.ones(lfwin)/lfwin,1,muTest.as_matrix())/step

    R2_Test[cell] = np.corrcoef(sm_yTest,sm_muTest)[0,1]**2
     
    f,(a1,a2)=plt.subplots(2,1, figsize=(12,10))
    a1.set_position([0.125, 0.38, 0.775, 0.53])
    a1=plotBehWindow(time_test[window_test],xTest.iloc[window_test].T,xTest.columns,a1)
    for tick in a1.get_yticklabels():
        tick.set_rotation(0)
    
    a1.set_xticks([])    
    a1.set_title('Cell_ID={0}, OoS_R2={1:0.2f}'.format(cell,R2_Test[cell]))

    a2.set_position([0.125, 0.2, 0.775, 0.15])
    a2=plotFitComp(sm_yTest[window_test],sm_muTest[window_test],
                   yTest[window_test]>0,time_test[window_test],a2 )
    a2.set_xlim(time_test[window_test][0],time_test[window_test][-1])
    f.savefig('/Users/alexgonzalez/Documents/Data/Li_T3g_062718_Results/GLM_Plots/GLM_FR_ReDe_Plot_Cell_{}.png'.format(cell),
              dpi=500,bbox_inches='tight',pad_inches=0.2)  
    plt.close(f)
plt.show()
M2<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18719591788299514, 'fraction_numerical': 0.03809184763260947, 'mean_word_length': 4.918188202247191, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '389128', 'license': 'MIT', 'n_tokens_mistral': 7557, 'n_tokens_neox': 7121, 'n_words': 716, 'path': 'Notebooks/GLM_Tests.ipynb', 'repo_name': 'alexgonzl/TreeMazeAnalyses', 'size': '16854'}