Changeset 2964

Timestamp:

Aug 7, 2017 8:18:52 AM (6 years ago)

Author:

vondreele

Message:

work on protein validation - not correct yet; gives plot tho.
fix refinement of Dij values in a LeBail? refinement (esp. sequential LeBail?).

Location:

branch/2frame

Files:

• GSASIImath.py (modified) (8 diffs)
• GSASIIphsGUI.py (modified) (2 diffs)
• GSASIIplot.py (modified) (1 diff)
• GSASIIstrIO.py (modified) (1 diff)

branch/2frame/GSASIImath.py

@@ -2584 +2584 @@
         'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL','MSE']
 # data taken from erratv2.ccp
-    lmt = np.array([17.190823041860433,11.526684477428809])
     b1 = np.array([     [0,0,0,0,0,0],
           [0,   5040.279078850848200,   3408.805141583649400,   4152.904423767300600,   4236.200004171890200,   5054.781210204625500],  
@@ -2624 +2623 @@
     chainIntAct = []
     res = []
+    resNames = []
+    resname = []
     newChain = True
     intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
@@ -2632 +2633 @@
                 resIntAct.append(sumintact(intact))
                 chainIntAct.append(resIntAct)
+                resNames += resname
                 res = []
+                resname = []
                 resIntAct = []
                 intact = {'CC':0,'CN':0,'CO':0,'NN':0,'NO':0,'OO':0,'NC':0,'OC':0,'ON':0}
@@ -2638 +2641 @@
         if atom[0] not in res:  #new residue, get residue no.
             res.append(atom[0])
+            resname.append('%s-%s%s'%(atom[2],atom[0],atom[1]))
             if not newChain:
                 resIntAct.append(sumintact(intact))
@@ -2658 +2662 @@
                 mult = 2.*(3.75-dsqt)
             intype = atom[4].strip()+cartAtoms[tgt][4].strip()
-            intact[intype] += mult            
+            intact[intype] += mult
+    resNames += resname
     resIntAct.append(sumintact(intact))
     chainIntAct.append(resIntAct)
@@ -2665 +2670 @@
         IntAct = chainIntAct[ich]
         nRes = len(IntAct)
-        Probs = []
+        Probs = [0.,0.,0.,0.]
         for i in range(4,nRes-4):
             mtrx = np.zeros(6)
@@ -2671 +2676 @@
             for j in range(i-4,i+4):
                 summ += np.sum(np.array(IntAct[j].values()))
+                print 'CC',IntAct[j]['CC'],'CN',IntAct[j]['CN'],'CO',IntAct[j]['CO']
+                print 'NN',IntAct[j]['NN'],'NO',IntAct[j]['NO'],'OO',IntAct[j]['OO']
                 mtrx[1] += IntAct[j]['CC']
                 mtrx[2] += IntAct[j]['CN']
@@ -2678 +2685 @@
             mtrx /= summ
             mtrx -= avg
-            prob = np.inner(np.inner(mtrx,b1),mtrx)
+            prob = np.inner(np.inner(mtrx,b1),mtrx)/36.
+            print i, mtrx
             Probs.append(prob)
-        chainProb.append(Probs)
-            
-        
-    print 'Do VALIDPROTEIN analysis - TBD'
-            
+        Probs += 4*[0.,]
+        chainProb += Probs
+    return resNames,chainProb
     
 ################################################################################

branch/2frame/GSASIIphsGUI.py

@@ -2103 +2103 @@
                 line3Sizer.Add(slope,0,WACV)
             elif 'Basin Hopping' in MCSAdata['Algorithm']:
-                pass        #TODO basinhopping controls here
+                pass        #TODO basinhopping controls here?
             mcsaSizer.Add(line3Sizer)
             mcsaSizer.Add((5,5),)
@@ -3291 +3291 @@
         
     def OnValidProtein(event):
-        G2mth.validProtein(data)
+        resNames,Probs = G2mth.validProtein(data)
+        G2plt.PlotAAProb(G2frame,resNames,Probs,Title='Error score for %s'%(data['General']['Name']),
+            thresh=[17.190823041860433,11.526684477428809])
 
     def OnSetAll(event):

branch/2frame/GSASIIplot.py

@@ -3085 +3085 @@
     else:
         Page.canvas.draw()
+        
+################################################################################
+##### PlotHist
+################################################################################
+def PlotAAProb(G2frame,resNames,Probs,Title='',thresh=None,newPlot=True):
+
+
+    global xylim
+    def OnMotion(event):
+        xpos = event.xdata
+        if xpos:                                        #avoid out of frame mouse position
+            ypos = int(event.ydata)
+            Page.canvas.SetCursor(wx.CROSS_CURSOR)
+            try:
+                if 0 <= ypos < len(resNames):
+                    G2frame.G2plotNB.status.SetStatusText('Residue: %s score: %.2f'%(resNames[ypos],xpos),1)
+            except TypeError:
+                G2frame.G2plotNB.status.SetStatusText('Select AA error plot first',1)
+    
+    def Draw():
+        global xylim
+        Plot.clear()
+        Plot.set_title(Title)
+        Plot.set_xlabel(r'Error score',fontsize=14)
+        Plot.set_ylabel(r'Residue',fontsize=14)
+        colors=['b','r','g','c','m','k']
+        Plot.barh(np.arange(len(resNames)),Probs,tick_label=resNames,)
+        if thresh is not None:
+            for item in thresh:
+                Plot.axvline(item,dashes=(5,5),picker=False)
+        if not newPlot:
+            Page.toolbar.push_current()
+            Plot.set_xlim(xylim[0])
+            Plot.set_ylim(xylim[1])
+            xylim = []
+            Page.toolbar.push_current()
+            Page.toolbar.draw()
+            Page.canvas.draw()
+        else:
+            Page.canvas.draw()
+            
+    
+    new,plotNum,Page,Plot,lim = G2frame.G2plotNB.FindPlotTab(Title,'mpl')
+    Page.Offset = [0,0]
+    if not new:
+        if not newPlot:
+            xylim = lim
+    else:
+        newPlot = True
+#        Page.canvas.mpl_connect('key_press_event', OnKeyPress)
+        Page.canvas.mpl_connect('motion_notify_event', OnMotion)
+        Page.Offset = [0,0]
+        Page.Choice = None
+    
+    Draw()
 
 ################################################################################

branch/2frame/GSASIIstrIO.py

@@ -2226 +2226 @@
                     hapDict[pfx+name] = hapData['HStrain'][0][i]
                     HSvals.append(hapDict[pfx+name])
-                    if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
+                    if hapData['HStrain'][1][i]:
+#                    if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
                         hapVary.append(pfx+name)
                 controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]