From eba0575cc15881ab4d6cfd9a834272548e3da0c3 Mon Sep 17 00:00:00 2001
From: alsinmr <51168380+alsinmr@users.noreply.github.com>
Date: Sun, 28 Jan 2024 12:11:27 +0100
Subject: [PATCH] PCA CC

---
 PCA/PCAclean.py | 10 ++++++--
 PCA/PCAsubs.py  | 64 +++++++++++++++++++++++++++++++++++++++++++++++++
 PCA/testPCA.py  |  8 +++++++
 3 files changed, 80 insertions(+), 2 deletions(-)

diff --git a/PCA/PCAclean.py b/PCA/PCAclean.py
index e3231eb..b951c37 100644
--- a/PCA/PCAclean.py
+++ b/PCA/PCAclean.py
@@ -18,7 +18,7 @@
 from copy import copy,deepcopy
 from ..misc.tools import linear_ex
 from .PCAmovies import PCAmovies
-from .PCAsubs import PCA_Ct,PCA_S2,PCAvecs,PCA2Data,Weighting
+from .PCAsubs import PCA_Ct,PCA_S2,PCAvecs,PCA2Data,Weighting,Impulse
 
 
 
@@ -82,6 +82,12 @@ def Weighting(self):
             self._Weighting=Weighting(self)
         return self._Weighting
     
+    @property
+    def Impulse(self):
+        if self._Impulse is None:
+            self._Impulse=Impulse(self)
+        return self._Impulse
+    
     #%% Misc.
     def clear(self):
         """
@@ -93,7 +99,7 @@ def clear(self):
 
         """
         keys=['_pos','_covar','_sel1index','_sel2index','_lambda','_PC','_pcamp','_mean',
-              '_S2','_Ct','_Vecs','_Data','_Movie','_Weighting']
+              '_S2','_Ct','_Vecs','_Data','_Movie','_Weighting','_Impulse']
         for k in keys:setattr(self,k,None)
         return self
     
diff --git a/PCA/PCAsubs.py b/PCA/PCAsubs.py
index 2221f1b..d6764ba 100644
--- a/PCA/PCAsubs.py
+++ b/PCA/PCAsubs.py
@@ -455,6 +455,8 @@ def S2m(self):
             self._S2m=S20m[1:]/S20m[:-1]
         return self._S2m[self.pca.Ct.tc_rev_index]
     
+
+    
     @property
     def Am(self):
         """
@@ -472,6 +474,54 @@ def Am(self):
             self._Am=S20m[:-1]-S20m[1:]
         return self._Am[self.pca.Ct.tc_rev_index]
     
+    def S20mCC(self,q):
+        """
+        Calculates the contribution to the total order parameter arising from
+        the mth principal component, noting that we number according to sorting
+        of the correlation times (so zero is the fastest PC, and -1 is the 
+        slowest PC)
+        
+        Note this is the product of all S2 up to mode m
+
+        Returns
+        -------
+        2D array (n PCs x n bonds)
+
+        """
+        
+
+        pca=self.pca
+        d20m=self.d2_0m
+        
+        PC=pca.PCxyz[:,:,pca.Ct.tc_index]
+        Lambda=pca.Lambda[pca.Ct.tc_index]
+        
+                
+        
+        X=np.zeros([pca.PC.shape[1],pca.nbonds])
+        for k in range(3):
+            for j in range(3):
+                P=(pca.mean[pca.sel1index[q],k]-pca.mean[pca.sel2index[q],k])*\
+                    (pca.mean[pca.sel1index,j]-pca.mean[pca.sel2index,j])
+                    
+                    
+                a=Lambda*(PC[k,pca.sel1index[q]]-PC[k,pca.sel2index[q]])*\
+                    (PC[j,pca.sel1index]-PC[j,pca.sel2index])
+                    
+                b=(P+np.cumsum(a.T,axis=0))/np.sqrt(d20m[:,q]*d20m.T).T
+                
+                # b=P/d20m[0]
+                
+                X+=(b**2)
+        S20m=-1/2+3/2*X
+        #Cleanup: make the S20m sorted
+        while np.any(S20m[:-1]<S20m[1:]):
+            i=S20m[:-1]<S20m[1:]
+            S20m[:-1][i]=S20m[1:][i]+1e-7
+            
+            
+        return S20m
+    
     def plotS2_direct_v_prod(self,ax=None):
         """
         Make a plot of S2 comparing the direct calculation and the product of
@@ -668,6 +718,20 @@ def bond(self,index:int,rho_index:int=None,PCAfit=None,frac:float=0.75):
         out=out[i]
         
         return out
+
+#%% Impulse response: A last attemp
+class Impulse():
+    def __init__(self,pca):
+        self.pca=pca
+        self._PCamp=None
+        
+    @property
+    def PCamp(self):
+        if self._PCamp is None:
+            l2=self.pca.PCamp.shape[1]//2
+            i=np.argmax(self.pca.PCamp[:,:l2],axis=-1)
+            self._PCamp=np.array([a[i0:i0+l2] for a,i0 in zip(self.pca.PCamp,i)])
+        return self._PCamp
     
 
 #%% Port various correlation functions to data objects        
diff --git a/PCA/testPCA.py b/PCA/testPCA.py
index b699426..a102b00 100644
--- a/PCA/testPCA.py
+++ b/PCA/testPCA.py
@@ -18,3 +18,11 @@
 pca.select_bond('15N')
 pca.select_atoms('name N C CA')
 
+
+
+
+ired=pyDR.md2iRED(select).iRED2data()
+ired.detect.r_no_opt(10)
+noopt=ired.fit()
+noopt.detect.r_auto(6)
+fit=noopt.fit().modes2bonds()