MIsc. updates

Entropy/Entropy.py

@@ -359,6 +359,7 @@ def v_avg(self):
             out=[]
             for v,chi,mult in zip(self.vt,self.chi,self.mult):
                 v=copy(v)
+                vout=np.zeros(v.shape[:-1])
 
                 for q in [2,3]:
                     for angle in [m*360/q for m in range(1,q)]:
@@ -367,9 +368,23 @@ def v_avg(self):
                         v1=-v[0,i]*np.sin(2*np.pi/q)+v[1,i]*np.cos(2*np.pi/q)
                         v[0,i]=v0
                         v[1,i]=v1
-                v=v.mean(-1)
-                v/=np.sqrt(v[0]**2+v[1]**2) #Renormalization
-                out.append(v)
+
+                for q in [2,3]:
+                    i=mult==q
+                    angle=np.arctan2(v[1,i],v[0,i])*q
+                    v0=np.median(np.cos(angle),axis=-1)
+                    v1=np.median(np.sin(angle),axis=-1)
+
+                    angle=np.arctan2(v1,v0)/q
+
+                    vout[0,i]=-np.cos(angle)
+                    vout[1,i]=-np.sin(angle)
+
+
+                # v=v.mean(-1)
+                # v/=np.sqrt(v[0]**2+v[1]**2) #Renormalization
+                # out.append(v)
+                out.append(vout)
             self._v_avg=out
 
         return self._v_avg
@@ -399,6 +414,7 @@ def state0(self):
                         overlap.append(v[0,i].T*vref[0]+v[1,i].T*vref[1])
 
                     state[-1][i]=np.argmax(overlap,axis=0).T
+
             self._state0=state
         return self._state0
 
@@ -831,6 +847,52 @@ def format_func(value,tick_number):
 
         return plt.FuncFormatter(format_func)
 
+
+    def plotChi(self,index:int,ax:list=None,step:int=1):
+        """
+        Creates one or more Ramachandran histogram plots, depending on the number
+        of chi angles, i.e. 1D histogram for Valine, 1 Ramachandran plot for
+        Isoleucine, 2 plots for Glutamine, etc.
+
+        Parameters
+        ----------
+        index : int
+            DESCRIPTION.
+
+        Returns
+        -------
+        fig
+
+        """
+        index=np.array(index)
+        if index.dtype==bool:index=np.argmax(index)
+
+        N=self.index[3:,index].sum()
+        nplots=max([1,N-1])
+        if ax is None:
+            fig,ax=plt.subplots(1,nplots)
+        ax=np.atleast_1d(ax).flatten()
+        assert len(ax)==nplots,f"Residue {self.resid[index].resname}{self.resid[index].resid} has {N} chi angles, and therefore requires {nplots} plots"
+
+        chi=[self.chi[k][self.index[k+3,:index].sum()][::step] for k in range(N-1,-1,-1)]
+
+        nstates=self.total_mult[index]
+        cmap=plt.get_cmap('turbo').resampled(nstates)
+
+        for a,k in zip(ax,range(N-1)):
+            for m in range(nstates):
+                i=self.state[index][::step]==m
+                a.scatter(chi[k][i],chi[k+1][i],color=cmap(m),s=1)
+            a.set_xlim([0,360])
+            a.set_ylim([0,360])
+            a.set_xlabel(rf'$\chi_{k+1}$ / $^\circ$')
+            a.set_ylabel(rf'$\chi_{k+2}$ / $^\circ$')
+
+        fig.tight_layout()        
+        return ax
+
+
+
 #%% Chimera functions
     def chimera(self,index=None,scaling:float=None,norm:bool=True,color=[1,0,0,1]):
         """

PCA/PCAclean.py

@@ -419,7 +419,7 @@ def align(ref0,ref,atoms):
         Ut=U.T
 
         R=V@Ut
-        return (R@pos.T).T
+        return (R.T@pos.T).T
 
 
     # def align(self):

PCA/PCAsubs.py

@@ -1209,11 +1209,13 @@ def project(self):
 
     @property
     def n_clusters(self):
+        #Changing order here because some algorithms do not take n_clusters as argument
+        if self._output is not None and hasattr(self._output,'n_clusters'):
+            return self._output.n_clusters
+
         if 'n_clusters' in self.cluster_kwargs:
             return self.cluster_kwargs['n_clusters']
 
-        if self._output is not None and hasattr(self._output,'n_clusters'):
-            return self._output.n_clusters
         return None
 
     @n_clusters.setter
@@ -1308,11 +1310,12 @@ def plot(self,ax=None,skip:int=10,maxbin:float=None,nbins:int=None,cmap='binary'
             for q in range(self.n_clusters):
                 a.scatter(self.pca.PCamp[self.index[k]][self.state==q][::skip],
                           self.pca.PCamp[self.index[k+1]][self.state==q][::skip],s=.1,color=cmap0(q))
-                a.scatter(self.pca.PCamp[self.index[k]][self.state==q].mean(),
-                       self.pca.PCamp[self.index[k+1]][self.state==q].mean(),s=25,
-                       color='black',marker='^')
-
                 if percent:
+                    a.scatter(self.pca.PCamp[self.index[k]][self.state==q].mean(),
+                           self.pca.PCamp[self.index[k+1]][self.state==q].mean(),s=25,
+                           color='black',marker='^')
+
+
                     a.text(self.pca.PCamp[self.index[k]][self.state==q].mean(),
                            self.pca.PCamp[self.index[k+1]][self.state==q].mean(),
                            f'{self.populations[q]*100:.1f}%',

Selection/MolSys.py

@@ -427,7 +427,7 @@ def lengths(self):
         l=(self.mda_traj.total_times/self.mda_traj.dt).astype(int)
         l[0]-=self.t0
         l[-1]-=self.__tf-self.tf
-        return l
+        return (l/self.step).astype(int)
 
 
     @property
@@ -1083,6 +1083,8 @@ def chimera(self,color:tuple=(1.,0.,0.,1.),x=None,index=None,norm:bool=False):
         if index is None:index=np.arange(len(self))
 
         if np.max(color)>1:color=[float(c/255) for c in color]
+        if len(color)==3:color=[*color,1.]
+
 
         if self.project is not None and self.project.chimera.saved_commands is not None:
             for cmd in self.project.chimera.saved_commands: