From 284c6e95fe2d7938ee6ee35cbf35590f9098e3de Mon Sep 17 00:00:00 2001 From: jcopperm Date: Thu, 5 Oct 2023 16:41:13 -0700 Subject: [PATCH] koopman! --- celltraj/.trajectory_legacy.swp | Bin 0 -> 4096 bytes celltraj/model.py | 475 ++++++++++++++++++++++---------- celltraj/trajectory.py | 73 ++++- celltraj/utilities.py | 4 +- environment.yml | 1 + 5 files changed, 403 insertions(+), 150 deletions(-) create mode 100644 celltraj/.trajectory_legacy.swp diff --git a/celltraj/.trajectory_legacy.swp b/celltraj/.trajectory_legacy.swp new file mode 100644 index 0000000000000000000000000000000000000000..7614192886f634a05ee0683838fc733c8e609511 GIT binary patch literal 4096 zcmYc?2=nw+u+TGN00IFJ0RaxI3=GNn1qG=^xrs0_Ty#))uz`iSkr^&|n7Dp!UWtCP zerR!OQL#SEaDD&0RF|U6vQ&MS{N&Qy)Vva)Txw1Zm?$Yq%tA0gywv29{G!VEoYeHh z0: - Mt[iR,:]=Mt[iR,:]/sM[iR] - if sM[iR]==0.0: - Mt[iR,iR]=1.0 - return Mt + n_clusters=clusters.clustercenters.shape[0] + indc0=clusters.assign(x0) + indc1=clusters.assign(x1) + Cm=np.zeros((n_clusters,n_clusters)) + for itt in range(x0.shape[0]): + Cm[indc0[itt],indc1[itt]]=Cm[indc0[itt],indc1[itt]]+1 + Mt=Cm.copy() + sM=np.sum(Mt,1) + for iR in range(n_clusters): + if sM[iR]>0: + Mt[iR,:]=Mt[iR,:]/sM[iR] + if sM[iR]==0.0: + Mt[iR,iR]=1.0 + return Mt def get_transition_matrix_CG(x0,x1,clusters,states): - n_clusters=clusters.clustercenters.shape[0] - n_states=np.max(states)+1 - indc0=states[clusters.assign(x0)] - indc1=states[clusters.assign(x1)] - Cm=np.zeros((n_states,n_states)) - for itt in range(x0.shape[0]): - Cm[indc0[itt],indc1[itt]]=Cm[indc0[itt],indc1[itt]]+1 - Mt=Cm.copy() - sM=np.sum(Mt,1) - for iR in range(n_states): - if sM[iR]>0: - Mt[iR,:]=Mt[iR,:]/sM[iR] - if sM[iR]==0.0: - Mt[iR,iR]=1.0 - return Mt + n_clusters=clusters.clustercenters.shape[0] + n_states=np.max(states)+1 + indc0=states[clusters.assign(x0)] + indc1=states[clusters.assign(x1)] + Cm=np.zeros((n_states,n_states)) + for itt in range(x0.shape[0]): + Cm[indc0[itt],indc1[itt]]=Cm[indc0[itt],indc1[itt]]+1 + Mt=Cm.copy() + sM=np.sum(Mt,1) + for iR in range(n_states): + if sM[iR]>0: + Mt[iR,:]=Mt[iR,:]/sM[iR] + if sM[iR]==0.0: + Mt[iR,iR]=1.0 + return Mt def clean_clusters(clusters,P): - centers=clusters.clustercenters.copy() - graph = csr_matrix(P>0.) - n_components, labels = connected_components(csgraph=graph, directed=False, return_labels=True) - unique, counts = np.unique(labels, return_counts=True) - icc=unique[np.argmax(counts)] - indcc=np.where(labels==icc)[0] - centers=centers[indcc,:] - clusters_clean=coor.clustering.AssignCenters(centers, metric='euclidean') - return clusters_clean + centers=clusters.clustercenters.copy() + graph = csr_matrix(P>0.) + n_components, labels = connected_components(csgraph=graph, directed=False, return_labels=True) + unique, counts = np.unique(labels, return_counts=True) + icc=unique[np.argmax(counts)] + indcc=np.where(labels==icc)[0] + centers=centers[indcc,:] + clusters_clean=coor.clustering.AssignCenters(centers, metric='euclidean') + return clusters_clean def get_path_entropy_2point(self,x0,x1,clusters,Mt,exclude_stays=False): - indc0=clusters.assign(x0) - indc1=clusters.assign(x1) - entp=0.0 - itt=0 - ntraj=indc0.size - try: - for itraj in range(ntraj): - if exclude_stays: - if Mt[indc0[itraj],indc1[itraj]]>0. and indc1[itraj]!=indc0[itraj]: - itt=itt+1 - pt=Mt[indc0[itraj],indc1[itraj]] - entp=entp-pt*np.log(pt) - else: - if Mt[indc0[itraj],indc1[itraj]]>0.: # and Mt[indc1[itraj],indc0[itraj]]>0.: - itt=itt+1 - pt=Mt[indc0[itraj],indc1[itraj]] - entp=entp-pt*np.log(pt) - entp=entp/(1.*itt) - except: - sys.stdout.write('empty arrays or failed calc\n') - entp=np.nan - return entp + indc0=clusters.assign(x0) + indc1=clusters.assign(x1) + entp=0.0 + itt=0 + ntraj=indc0.size + try: + for itraj in range(ntraj): + if exclude_stays: + if Mt[indc0[itraj],indc1[itraj]]>0. and indc1[itraj]!=indc0[itraj]: + itt=itt+1 + pt=Mt[indc0[itraj],indc1[itraj]] + entp=entp-pt*np.log(pt) + else: + if Mt[indc0[itraj],indc1[itraj]]>0.: # and Mt[indc1[itraj],indc0[itraj]]>0.: + itt=itt+1 + pt=Mt[indc0[itraj],indc1[itraj]] + entp=entp-pt*np.log(pt) + entp=entp/(1.*itt) + except: + sys.stdout.write('empty arrays or failed calc\n') + entp=np.nan + return entp def get_path_ll_2point(self,x0,x1,exclude_stays=False): - indc0=clusters.assign(x0) - indc1=clusters.assign(x1) - ll=0.0 - itt=0 - ntraj=indc0.size - try: - for itraj in range(ntraj): - if exclude_stays: - if Mt[indc0[itraj],indc1[itraj]]>0. and indc1[itraj]!=indc0[itraj]: - itt=itt+1 - pt=Mt[indc0[itraj],indc1[itraj]] - ll=ll+np.log(pt) - else: - if Mt[indc0[itraj],indc1[itraj]]>0.: # and Mt[indc1[itraj],indc0[itraj]]>0.: - itt=itt+1 - pt=Mt[indc0[itraj],indc1[itraj]] - ll=ll+np.log(pt) - ll=ll/(1.*itt) - except: - sys.stdout.write('empty arrays or failed calc\n') - ll=np.nan - return ll + indc0=clusters.assign(x0) + indc1=clusters.assign(x1) + ll=0.0 + itt=0 + ntraj=indc0.size + try: + for itraj in range(ntraj): + if exclude_stays: + if Mt[indc0[itraj],indc1[itraj]]>0. and indc1[itraj]!=indc0[itraj]: + itt=itt+1 + pt=Mt[indc0[itraj],indc1[itraj]] + ll=ll+np.log(pt) + else: + if Mt[indc0[itraj],indc1[itraj]]>0.: # and Mt[indc1[itraj],indc0[itraj]]>0.: + itt=itt+1 + pt=Mt[indc0[itraj],indc1[itraj]] + ll=ll+np.log(pt) + ll=ll/(1.*itt) + except: + sys.stdout.write('empty arrays or failed calc\n') + ll=np.nan + return ll def get_kscore(self,Mt,eps=1.e-3): #,nw=10): - indeye=np.where(np.eye(Mt.shape[0])) - diag=Mt[indeye] - indgood=np.where(diag<1.)[0] - Mt=Mt[indgood,:] - Mt=Mt[:,indgood] - w,v=np.linalg.eig(np.transpose(Mt)) - w=np.real(w) - if np.sum(np.abs(w-1.)0: - indw=np.where(np.logical_and(np.logical_and(np.abs(w-1.)>eps,w>0.),w<1.)) - tw=w[indw] - tw=np.sort(tw) - #tw=tw[-nw:] - tw=1./(1.-tw) - kscore=np.sum(tw) - else: - kscore=np.nan - return kscore + indeye=np.where(np.eye(Mt.shape[0])) + diag=Mt[indeye] + indgood=np.where(diag<1.)[0] + Mt=Mt[indgood,:] + Mt=Mt[:,indgood] + w,v=np.linalg.eig(np.transpose(Mt)) + w=np.real(w) + if np.sum(np.abs(w-1.)0: + indw=np.where(np.logical_and(np.logical_and(np.abs(w-1.)>eps,w>0.),w<1.)) + tw=w[indw] + tw=np.sort(tw) + #tw=tw[-nw:] + tw=1./(1.-tw) + kscore=np.sum(tw) + else: + kscore=np.nan + return kscore def get_traj_ll_gmean(self,xt,exclude_stays=False,states=None): - if states is None: - states=np.arange(Mt.shape[0]).astype(int) - x0=xt[0:-1] - x1=xt[1:] - indc0=states[clusters.assign(x0)] - indc1=states[clusters.assign(x1)] - llSet=np.array([]) - itt=0 - ntraj=indc0.size - try: - for itraj in range(ntraj): - if exclude_stays: - if Mt[indc0[itraj],indc1[itraj]]>0. and indc1[itraj]!=indc0[itraj]: - itt=itt+1 - pt=Mt[indc0[itraj],indc1[itraj]] - llSet=np.append(llSet,pt) - else: - if Mt[indc0[itraj],indc1[itraj]]>0.: # and Mt[indc1[itraj],indc0[itraj]]>0.: - itt=itt+1 - pt=Mt[indc0[itraj],indc1[itraj]] - llSet=np.append(llSet,pt) - ll_mean=scipy.stats.mstats.gmean(llSet) - except: - sys.stdout.write('empty arrays or failed calc\n') - ll_mean=np.nan - return ll_mean + if states is None: + states=np.arange(Mt.shape[0]).astype(int) + x0=xt[0:-1] + x1=xt[1:] + indc0=states[clusters.assign(x0)] + indc1=states[clusters.assign(x1)] + llSet=np.array([]) + itt=0 + ntraj=indc0.size + try: + for itraj in range(ntraj): + if exclude_stays: + if Mt[indc0[itraj],indc1[itraj]]>0. and indc1[itraj]!=indc0[itraj]: + itt=itt+1 + pt=Mt[indc0[itraj],indc1[itraj]] + llSet=np.append(llSet,pt) + else: + if Mt[indc0[itraj],indc1[itraj]]>0.: # and Mt[indc1[itraj],indc0[itraj]]>0.: + itt=itt+1 + pt=Mt[indc0[itraj],indc1[itraj]] + llSet=np.append(llSet,pt) + ll_mean=scipy.stats.mstats.gmean(llSet) + except: + sys.stdout.write('empty arrays or failed calc\n') + ll_mean=np.nan + return ll_mean def get_H_eigs(Mt): @@ -312,19 +312,214 @@ def plot_eig(v,x_clusters,ncomp): plt.pause(1); def get_dmat(x1,x2=None): #adapted to python from Russell Fung matlab implementation (github.com/ki-analysis/manifold-ga dmat.m) - x1=np.transpose(x1) #default from Fung folks is D x N - if x2 is None: - nX1 = x1.shape[1]; - y = np.matlib.repmat(np.sum(np.power(x1,2),0),nX1,1) - y = y - np.matmul(np.transpose(x1),x1) - y = y + np.transpose(y); - y = np.abs( y + np.transpose(y) ) / 2. # Iron-out numerical wrinkles - else: - x2=np.transpose(x2) - nX1 = x1.shape[1] - nX2 = x2.shape[1] - y = np.matlib.repmat( np.expand_dims(np.sum( np.power(x1,2), 0 ),1), 1, nX2 ) - y = y + np.matlib.repmat( np.sum( np.power(x2,2), 0 ), nX1, 1 ) - y = y - 2 * np.matmul(np.transpose(x1),x2) - return np.sqrt(y) + x1=np.transpose(x1) #default from Fung folks is D x N + if x2 is None: + nX1 = x1.shape[1]; + y = np.matlib.repmat(np.sum(np.power(x1,2),0),nX1,1) + y = y - np.matmul(np.transpose(x1),x1) + y = y + np.transpose(y); + y = np.abs( y + np.transpose(y) ) / 2. # Iron-out numerical wrinkles + else: + x2=np.transpose(x2) + nX1 = x1.shape[1] + nX2 = x2.shape[1] + y = np.matlib.repmat( np.expand_dims(np.sum( np.power(x1,2), 0 ),1), 1, nX2 ) + y = y + np.matlib.repmat( np.sum( np.power(x2,2), 0 ), nX1, 1 ) + y = y - 2 * np.matmul(np.transpose(x1),x2) + return np.sqrt(y) + +####################################feature tuned kernel DMD a la aristoff######################## +def get_kernel_sigmas(X,M,s=.05): + """Get sigmas from observation matrix + + Parameters + ---------- + X : ndarray + observation matrix, samples x features + M : ndarray + Mahalanobis scaling matrix, features x features + s : float + bandwidth scaling factor + Returns + ------- + h : ndarray, n_features (float) + vector of sigmas to scale observations in kernel + """ + XM=np.matmul(X,M) + h=s*np.std(scipy.spatial.distance.cdist(XM,XM,metric='euclidean'),axis=0) + return h + +def get_gaussianKernelM(X,Y,M,h): + """Get Malanobis scaled gaussian kernel from observation matrices X,Y + + Parameters + ---------- + X : ndarray + observation matrix, samples x features + Y : ndarray + observation matrix at t+1, samples x features + M : ndarray + Mahalanobis scaling matrix, features x features + h : ndarray + vector of sigma scalings for gaussian from get_kernel_sigmas + Returns + ------- + k : ndarray, samples x samples (float) + Malanobis scaled kernel for X,Y + """ + XM=np.matmul(X,M) + YM=np.matmul(Y,M) + k=np.exp(-np.divide(np.power(scipy.spatial.distance.cdist(YM,XM),2),2*h)) + return k + +def get_koopman_eig(X,Y,M=None,s=.05,bta=1.e-5,h=None,psi_X=None,psi_Y=None): + """Get linear matrix solution for Koopman operator from X,Y paired observation Y=F(X) with F the forward operator + + Parameters + ---------- + X : ndarray + observation matrix, samples x features + Y : ndarray + observation matrix at t+1, samples x features + M : ndarray + Mahalanobis scaling matrix, features x features + s : float + kernel bandwidth scaling parameter + bta : float + regularization parameter for linear solve + Returns + ------- + K : ndarray + Koopman operator matrix, samples x samples (float) + Xi : ndarray + Koopman left eigenvectors, samples x samples + Lam : ndarray + Koopman eigenvalue matrix (diagonal), samples x samples + W : ndarray + Koopman right eigenvectors, samples x samples + """ + nsamples=X.shape[0] + if M is None: + M=np.eye(X.shape[1]) + if h is None: + print('getting kernel sigmas...') + h=get_kernel_sigmas(X,M) + if psi_X is None: + print('applying kernel to X...') + psi_X=get_gaussianKernelM(X,X,M,h) + if psi_Y is None: + print('applying kernel to Y...') + psi_Y=get_gaussianKernelM(X,Y,M,h) + print('solving linear system for approximate Koopman...') + A = (psi_X+bta*np.eye(nsamples)) + K,residuals,rank,s = np.linalg.lstsq(A, psi_Y) + print('getting Koopman eigendecomposition...') + Lam, W, Xi = scipy.linalg.eig(K,left=True,right=True) + indsort=np.argsort(np.abs(Lam))[::-1] + Xi=Xi[:,indsort] + W=W[:,indsort] + Lam=np.diag(Lam[indsort]) + return K,Xi,Lam,W + +def get_koopman_modes(psi_X,Xi,W,X_obs): + """Get Koopman modes of an observable + + Parameters + ---------- + psi_X : ndarray + kernel matrix, samples x samples + Xi : ndarray + right eigenvectors of Koopman, samples x samples + W : ndarray + left eigenvectors of Koopman, samples x samples + X : ndarray + observation matrix, samples x features + X_obs : ndarray, samples x observables + observables of interest in same time order as X, can be X or features of X + Returns + ------- + phi_X : ndarray + Koopman eigenfunctions + V : ndarray + Koopman modes of observables + """ + phi_X=np.matmul(psi_X,Xi) + B = np.matmul(np.linalg.pinv(psi_X),X_obs) + Wprime = np.divide(np.conj(W.T),np.diag(np.matmul(np.conj(W.T),Xi))[:,np.newaxis]) + V=np.matmul(Wprime,B) + return phi_X,V + +def get_koopman_inference(start,steps,phi_X,V,Lam): + """Get Koopman prediction of an observable + + Parameters + ---------- + start : int + sample index for start point + steps : int + number of steps of inference to perform + phi_X : ndarray + Koopman eigenfunctions, samples x samples + V : ndarray + Koopman modes of observables, must be calculated samples x samples + Lam : ndarray + Koopman eigenvalues matrix (diagonal), samples x samples + Returns + ------- + X_pred : ndarray + predicted trajectory steps x observables + """ + d = V.shape[1] + nmodes=V.shape[0] + lam = Lam[0:nmodes,:] + lam = lam[:,0:nmodes] + D = np.eye(nmodes) + X_pred = np.zeros((steps,d)).astype('complex128') + for step in range(steps): + #X_pred[step,:] = np.matmul(np.matmul(phi_X[start,:],D),V) + lambdas=np.diag(D) + X_pred[step,:] = np.matmul(np.multiply(phi_X[start,:],lambdas)[np.newaxis,:],V) + D = np.matmul(D,lam) + return X_pred + +def update_mahalanobis_matrix(Mprev,X,phi_X,nmodes=2,h=None,s=.05): + """Update estimation of mahalanobis matrix for kernel tuning + + Parameters + ---------- + Mprev : ndarray, features x features + Koopman eigenfunctions + X : ndarray + samples by features + phi_X : ndarray + Koopman eigenfunctions, samples x samples + Returns + ------- + M : ndarray + updated mahalanobis matrix using Koopman eigenfunction gradients + """ + #define gradient of Koopman eigenfunctions + #dphi = @(x,efcn) sum((X(1:N-1,:)-x)*M.*(k(x)'.*Phi_x(:,efcn))); + #compute M as the gradient outerproduct + if h is None: + h = get_kernel_sigmas(X,Mprev,s=s) + M = np.zeros_like(Mprev) + N = X.shape[0] + for imode in range(nmodes): + print(f'updating M with gradient of mode {imode}') + for n in range(N): + x=X[n,:] + x=x[np.newaxis,:] + xMX=np.matmul(X-x,Mprev) + kxX=get_gaussianKernelM(X,x,Mprev,h) + xMX_kxX=np.multiply(xMX,np.conj(kxX).T) + xMX_kxX_phi=np.multiply(xMX_kxX,phi_X[:,[imode]]) + grad = np.sum(xMX_kxX_phi,axis=0)[np.newaxis,:] + M = M + np.matmul(np.conj(grad.T),grad) + #get square root and regularize M + M = scipy.linalg.sqrtm(M) + M = np.real(M) + svdnorm=np.max(scipy.linalg.svdvals(M)) + M = M/svdnorm + return M diff --git a/celltraj/trajectory.py b/celltraj/trajectory.py index cb86156..8894921 100644 --- a/celltraj/trajectory.py +++ b/celltraj/trajectory.py @@ -504,12 +504,15 @@ def get_cell_data(self,ic,frametype='boundingbox',boundary_expansion=None,return labelc = mskc[...,relabel_mskchannels] else: labelc = mskc[...,np.newaxis] + label_table=regionprops_table(msk[...,self.mskchannel],intensity_image=None,properties=['label']) #table of labels for cells in frame + labelIDs=label_table['label'][indcells] #integer labels of cells in place for ichannel in range(len(relabel_mskchannels)): - label_table=regionprops_table(msk[...,relabel_mskchannels[ichannel]],intensity_image=None,properties=['label']) #table of labels for cells in frame - labelIDs=label_table['label'][indcells] #integer labels of cells in place + label_table_channel=regionprops_table(msk[...,relabel_mskchannels[ichannel]],intensity_image=None,properties=['label']) #table of labels for cells in frame + labelIDs_channel=label_table_channel['label'] #integer labels of cells in place + labelIDs_common=np.intersect1d(labelIDs,labelIDs_channel) mskc_channel=np.zeros_like(labelc[...,ichannel]).astype(int) - for ic_relabel in range(labelIDs.size): - icell=labelIDs[ic_relabel] + for ic_relabel in range(labelIDs_common.size): + icell=labelIDs_common[ic_relabel] icell_global=indcells_global[ic_relabel] mskc_channel[labelc[...,ichannel]==icell]=icell_global mskc[...,relabel_mskchannels[ichannel]]=mskc_channel @@ -688,18 +691,24 @@ def get_cell_compartment_ratio(self,indcells=None,imgchannel=None,mskchannel1=No msk1=msk[...,mskchannel1] msk2=msk[...,mskchannel2] if combined_and_disjoint: - msk2=imprep.get_cyto_minus_nuc_labels(msk2,msk1) + msk2=imprep.get_cyto_minus_nuc_labels(msk1,msk2) elif make_disjoint: - msk2[msk1>0]=0 + msk1[msk2>0]=0 if not np.all(np.unique(msk1)==np.unique(msk2)): - print(f'frame {self.cells_indimgSet[ic]} mask1 and mask2 yielding different indices. try combined_and_disjoint=True, or indexing may be off') - return 1 + print(f'warning: frame {self.cells_indimgSet[ic]} mask1 and mask2 yielding different indices') if erosion_footprint1 is not None: fmsk1=skimage.morphology.binary_erosion(msk1>0,footprint=erosion_footprint1); msk1[fmsk1]=0 if erosion_footprint2 is not None: fmsk2=skimage.morphology.binary_erosion(msk2>0,footprint=erosion_footprint2); msk2[fmsk2]=0 props1 = regionprops_table(msk1, intensity_image=img,properties=('label','intensity_mean','area')) props2 = regionprops_table(msk2, intensity_image=img,properties=('label','intensity_mean','area')) + commonlabels,indcommon1,indcommon2=np.intersect1d(props1['label'],props2['label'],return_indices=True) + props2_matched=props1.copy() + props2_matched['intensity_mean']=np.ones_like(props1['intensity_mean'])*np.nan + props2_matched['intensity_mean'][indcommon1]=props2['intensity_mean'][indcommon2] + props2_matched['area']=np.ones_like(props1['area'])*np.nan + props2_matched['area'][indcommon1]=props2['area'][indcommon2] + props2=props2_matched if intensity_sum: if intensity_ztransform: cratio_frame=np.divide(self.img_zstds[imgchannel]*np.multiply(props1['intensity_mean']+self.img_zmeans[imgchannel],props1['area']),self.img_zstds[imgchannel]*np.multiply(props2['intensity_mean']+self.img_zmeans[imgchannel],props2['area'])) @@ -1113,6 +1122,54 @@ def get_unique_trajectories(self,cell_inds=None,verbose=False,extra_depth=None): sys.stdout.write('tracked cell '+str(indc)+', '+str(cell_traj.size)+' tracks, '+str(n_untracked)+' left\n') self.trajectories=trajectories + def get_traj_segments(self,seg_length): + ntraj=len(self.trajectories) + traj_segSet=np.zeros((0,seg_length)).astype(int) + for itraj in range(ntraj): + cell_traj=self.trajectories[itraj] + traj_len=cell_traj.size + if traj_len>=seg_length: + for ic in range(traj_len-seg_length+1): #was -1, think that was an error, changed 2june21 because ended up missing data + traj_seg=cell_traj[ic:ic+seg_length] + traj_segSet=np.append(traj_segSet,traj_seg[np.newaxis,:],axis=0) + return traj_segSet + + def get_trajectory_steps(self,inds=None,traj=None,Xtraj=None,get_trajectories=True,nlag=1): #traj and Xtraj should be indexed same + if inds is None: + inds=np.arange(self.cells_indSet.size).astype(int) + if get_trajectories: + self.get_unique_trajectories(cell_inds=inds) + if traj is None: + traj=self.traj + if Xtraj is None: + x=self.Xtraj + else: + x=Xtraj + trajp1=self.get_traj_segments(self.trajl+nlag) + inds_nlag=np.flipud(np.arange(self.trajl+nlag-1,-1,-nlag)).astype(int) #keep indices every nlag + trajp1=trajp1[:,inds_nlag] + ntraj=trajp1.shape[0] + neigen=x.shape[1] + x0=np.zeros((0,neigen)) + x1=np.zeros((0,neigen)) + inds_trajp1=np.zeros((0,2)).astype(int) + for itraj in range(ntraj): + test0=trajp1[itraj,0:-1] + test1=trajp1[itraj,1:] + res0 = (traj[:, None] == test0[np.newaxis,:]).all(-1).any(-1) + res1 = (traj[:, None] == test1[np.newaxis,:]).all(-1).any(-1) + if np.sum(res0)==1 and np.sum(res1)==1: + indt0=np.where(res0)[0][0] + indt1=np.where(res1)[0][0] + x0=np.append(x0,np.array([x[indt0,:]]),axis=0) + x1=np.append(x1,np.array([x[indt1,:]]),axis=0) + inds_trajp1=np.append(inds_trajp1,np.array([[indt0,indt1]]),axis=0) + if itraj%100==0: + sys.stdout.write('matching up trajectory '+str(itraj)+'\n') + self.Xtraj0=x0 + self.Xtraj1=x1 + self.inds_trajp1=inds_trajp1 + def get_pair_rdf(self,cell_indsA=None,cell_indsB=None,rbins=None,nr=50,rmax=500): if cell_indsA is None: cell_indsA=np.arange(self.cells_indSet.shape[0]).astype(int) diff --git a/celltraj/utilities.py b/celltraj/utilities.py index e943e60..12fc36b 100644 --- a/celltraj/utilities.py +++ b/celltraj/utilities.py @@ -122,7 +122,7 @@ def recursively_save_dict_contents_to_group( h5file, path, dic): try: #c h5file[path + element_name] = element #c except: #c - element = np.array(element).astype('|S9') #c + element = np.array(element).astype('|S32') #c h5file[path + element_name] = element #c else: #c raise ValueError('Cannot save %s type within a list.' % type(element)) #c @@ -140,7 +140,7 @@ def recursively_save_dict_contents_to_group( h5file, path, dic): try: h5file[path + key] = item except: - item = np.array(item).astype('|S9') + item = np.array(item).astype('|S32') h5file[path + key] = item if not np.array_equal(h5file[path + key][()], item): raise ValueError('The data representation in the HDF5 file does not match the original dict.') diff --git a/environment.yml b/environment.yml index df86fbc..8a919bd 100644 --- a/environment.yml +++ b/environment.yml @@ -21,6 +21,7 @@ dependencies: - pytorch - pandas - dask + - msmtools - pip: - csaps