Changeset 4079

Timestamp:

Aug 5, 2019 3:33:15 PM (3 years ago)

Author:

vondreele

Message:

modify SHAPES for better output to GSAS-II
begin plotting capability for SHAPES output

Location:

trunk

Files:

• G2shapes.py (modified) (8 diffs)
• GSASIIplot.py (modified) (7 diffs)
• GSASIIpwdGUI.py (modified) (5 diffs)

trunk/G2shapes.py

@@ -1083 +1083 @@
     print (aString)
     
-    aList_summary = []
-    aList_summary.append(version_aString)
-    aList_summary.append(str(localtime))
+#    aList_summary = []
+#    aList_summary.append(version_aString)
+#    aList_summary.append(str(localtime))
     
     ######################
@@ -1377 +1377 @@
                          
         aString  = 'Trial:' + file_no
-        aList_summary.append(aString)
+#        aList_summary.append(aString)
     
         file_beads = prefix + 'beads_' + file_no + '.pdb'
-        file_pr = prefix + 'pr_calc_' + file_no + '.dat'
+#        file_pr = prefix + 'pr_calc_' + file_no + '.dat'
         file_psv = prefix + 'psv_shape_' + file_no + '.pdb'
 #        file_intensity = prefix + 'intensity_' + file_no + '.dat'
@@ -1780 +1780 @@
         print (aString)
     
+        Phases.append([file_beads.split('.')[0],aList_beads_x,aList_beads_y,aList_beads_z])
         # Write out beads as pseudo a PDB file
-        Phases.append({'pname':file_beads.split('.')[0],
-            'Atoms':[aList_beads_x,aList_beads_y,aList_beads_z,file_beads]})
         if pdbOut:
             pdb_writer(aList_beads_x,aList_beads_y,aList_beads_z,file_beads,aString)
@@ -1830 +1829 @@
         # Write input and model P(r)
 #        pr_writer(aList_pr,aList_r,aList_pr_model,file_pr)
-        PRcalc.append([aList_r,aList_pr,aList_pr_model,])
+        PRcalc.append([aList_r,aList_pr,aList_pr_model,delta_hist_sum])
     
         # Calculate comparison versus intensities
@@ -1848 +1847 @@
     
             # Write output intensity file
-            Patterns.append([aList_q,aList_i,aList_i_calc,aString])
+            Patterns.append([aList_q,aList_i,aList_i_calc,rvalue])
 #            write_all_data(file_intensity,aList_q,aList_i,aList_i_calc,aString)
     
@@ -1869 +1868 @@
                   ' PSV:' + str(fraction_psv)
     
-        Phases.append({'pname':file_psv.split('.')[0],
-            'Atoms':[aList_box_x,aList_box_y,aList_box_z]})
+        Phases.append([file_psv.split('.')[0],aList_box_x,aList_box_y,aList_box_z])
         if pdbOut:
             pdb_writer(aList_box_x,aList_box_y,aList_box_z,file_psv,aString)
@@ -1878 +1876 @@
         aString = 'P(r) dif:' + str(delta_hist_sum) + ' E:' \
                          + str(vdw_all) + ' CHISQ:' + str(chi_sq) + ' PSV:' + str(fraction_psv)
-        aList_summary.append(aString)               
+#        aList_summary.append(aString)               
     
         i_soln = i_soln + 1
@@ -1901 +1899 @@
     print (aString)
                          
-    aList_summary.append(str(localtime))
-    
-    # Create summary
-    
-    aFile_log = prefix + 'shapes_summary.log'
-    num_lines = len(aList_summary)
-    
-    file = open(aFile_log,'w')
-    i = 0
-    while i < num_lines:
-        aString = str(aList_summary[i])
-        file.write(aString)
-        file.write('\n')
-        i = i + 1
-    file.close()
+#    aList_summary.append(str(localtime))
+#    
+#    # Create summary
+#    
+#    aFile_log = prefix + 'shapes_summary.log'
+#    num_lines = len(aList_summary)
+#    
+#    file = open(aFile_log,'w')
+#    i = 0
+#    while i < num_lines:
+#        aString = str(aList_summary[i])
+#        file.write(aString)
+#        file.write('\n')
+#        i = i + 1
+#    file.close()
 
     return Phases,Patterns,PRcalc

trunk/GSASIIplot.py

@@ -4679 +4679 @@
         Ymax = np.amax(XYlist.T[1])
         dy = 0.02*(Ymax-Ymin)
-        Plot.set_xlim(Xmin-dx,Xmax+dx)
-        Plot.set_ylim(Ymin-dy,Ymax+dy)
+        try:
+            Plot.set_xlim(Xmin-dx,Xmax+dx)
+            Plot.set_ylim(Ymin-dy,Ymax+dy)
+        except:
+            pass
         if Peaks == None:
             normcl = mpcls.Normalize(Ymin,Ymax)
@@ -5225 +5228 @@
         elif 'Pair' in PlotText:
             PlotSASDPairDist(G2frame)
+            
     
     def OnMotion(event):
@@ -9134 +9138 @@
 
 ################################################################################
+#### Plot Bead Model
+################################################################################
+
+def PlotBeadModel(G2frame,Atoms,defaults):
+    '''Bead modelplotting package. For bead models from SHAPES
+    '''
+
+    Mydir = G2frame.dirname
+    Rd = np.array([255,0,0])
+    Gr = np.array([0,255,0])
+    Bl = np.array([0,0,255])
+    uBox = np.array([[0,0,0],[1,0,0],[0,1,0],[0,0,1]])
+    uEdges = np.array([[uBox[0],uBox[1]],[uBox[0],uBox[2]],[uBox[0],uBox[3]]])
+    uColors = [Rd,Gr,Bl]
+    XYZ = np.array(Atoms[1:]).T      #don't mess with original!
+
+#    def SetRBOrigin():
+#        page = getSelection()
+#        if page:
+#            if G2frame.GetPageText(page) == 'Rigid bodies':
+#                G2frame.MapPeaksTable.SetData(mapPeaks)
+#                panel = G2frame.GetPage(page).GetChildren()
+#                names = [child.GetName() for child in panel]
+#                panel[names.index('grid window')].Refresh()
+            
+    def OnMouseDown(event):
+        xy = event.GetPosition()
+        defaults['oldxy'] = list(xy)
+
+    def OnMouseMove(event):
+        newxy = event.GetPosition()
+                                
+        if event.Dragging():
+            if event.LeftIsDown():
+                SetRotation(newxy)
+                Q = defaults['Quaternion']
+                G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
+#            elif event.RightIsDown():
+#                SetRBOrigin(newxy)
+            elif event.MiddleIsDown():
+                SetRotationZ(newxy)
+                Q = defaults['Quaternion']
+                G2frame.G2plotNB.status.SetStatusText('New quaternion: %.2f+, %.2fi+ ,%.2fj+, %.2fk'%(Q[0],Q[1],Q[2],Q[3]),1)
+            Draw('move')
+        
+    def OnMouseWheel(event):
+        defaults['cameraPos'] += event.GetWheelRotation()/24
+        defaults['cameraPos'] = max(10,min(500,defaults['cameraPos']))
+        G2frame.G2plotNB.status.SetStatusText('New camera distance: %.2f'%(defaults['cameraPos']),1)
+        Draw('wheel')
+        
+    def SetBackground():
+        R,G,B,A = Page.camera['backColor']
+        GL.glClearColor(R,G,B,A)
+        GL.glClear(GL.GL_COLOR_BUFFER_BIT | GL.GL_DEPTH_BUFFER_BIT)
+        
+    def SetLights():
+        GL.glEnable(GL.GL_DEPTH_TEST)
+        GL.glShadeModel(GL.GL_FLAT)
+        GL.glEnable(GL.GL_LIGHTING)
+        GL.glEnable(GL.GL_LIGHT0)
+        GL.glLightModeli(GL.GL_LIGHT_MODEL_TWO_SIDE,0)
+        GL.glLightfv(GL.GL_LIGHT0,GL.GL_AMBIENT,[1,1,1,.8])
+        GL.glLightfv(GL.GL_LIGHT0,GL.GL_DIFFUSE,[1,1,1,1])
+        
+    def SetRotation(newxy):
+#first get rotation vector in screen coords. & angle increment        
+        oldxy = defaults['oldxy']
+        if not len(oldxy): oldxy = list(newxy)
+        dxy = newxy-oldxy
+        defaults['oldxy'] = list(newxy)
+        V = np.array([dxy[1],dxy[0],0.])
+        A = 0.25*np.sqrt(dxy[0]**2+dxy[1]**2)
+# next transform vector back to xtal coordinates via inverse quaternion
+# & make new quaternion
+        Q = defaults['Quaternion']
+        V = G2mth.prodQVQ(G2mth.invQ(Q),V)
+        DQ = G2mth.AVdeg2Q(A,V)
+        Q = G2mth.prodQQ(Q,DQ)
+        defaults['Quaternion'] = Q
+# finally get new view vector - last row of rotation matrix
+        VD = G2mth.Q2Mat(Q)[2]
+        VD /= np.sqrt(np.sum(VD**2))
+        defaults['viewDir'] = VD
+        
+    def SetRotationZ(newxy):                        
+#first get rotation vector (= view vector) in screen coords. & angle increment        
+        View = GL.glGetIntegerv(GL.GL_VIEWPORT)
+        cent = [View[2]/2,View[3]/2]
+        oldxy = defaults['oldxy']
+        if not len(oldxy): oldxy = list(newxy)
+        dxy = newxy-oldxy
+        defaults['oldxy'] = list(newxy)
+        V = defaults['viewDir']
+        A = [0,0]
+        A[0] = dxy[1]*.25
+        A[1] = dxy[0]*.25
+        if newxy[0] > cent[0]:
+            A[0] *= -1
+        if newxy[1] < cent[1]:
+            A[1] *= -1        
+# next transform vector back to xtal coordinates & make new quaternion
+        Q = defaults['Quaternion']
+        Qx = G2mth.AVdeg2Q(A[0],V)
+        Qy = G2mth.AVdeg2Q(A[1],V)
+        Q = G2mth.prodQQ(Q,Qx)
+        Q = G2mth.prodQQ(Q,Qy)
+        defaults['Quaternion'] = Q
+
+    def RenderUnitVectors(x,y,z):
+        GL.glEnable(GL.GL_COLOR_MATERIAL)
+        GL.glLineWidth(1)
+        GL.glPushMatrix()
+        GL.glTranslate(x,y,z)
+        GL.glBegin(GL.GL_LINES)
+        for line,color in zip(uEdges,uColors):
+            GL.glColor3ubv(color)
+            GL.glVertex3fv(-line[1])
+            GL.glVertex3fv(line[1])
+        GL.glEnd()
+        GL.glPopMatrix()
+        GL.glColor4ubv([0,0,0,0])
+        GL.glDisable(GL.GL_COLOR_MATERIAL)
+                
+    def RenderSphere(x,y,z,radius,color):
+        GL.glMaterialfv(GL.GL_FRONT_AND_BACK,GL.GL_DIFFUSE,color)
+        GL.glPushMatrix()
+        GL.glTranslate(x,y,z)
+        q = GLU.gluNewQuadric()
+        GLU.gluSphere(q,radius,20,10)
+        GL.glPopMatrix()
+        
+    def Draw(caller=''):
+#useful debug?        
+#        if caller:
+#            print caller
+# end of useful debug
+        cPos = defaults['cameraPos']
+        VS = np.array(Page.canvas.GetSize())
+        aspect = float(VS[0])/float(VS[1])
+        Q = defaults['Quaternion']
+        SetBackground()
+        GL.glInitNames()
+        GL.glPushName(0)
+        
+        GL.glMatrixMode(GL.GL_PROJECTION)
+        GL.glLoadIdentity()
+        GL.glViewport(0,0,VS[0],VS[1])
+        GLU.gluPerspective(20.,aspect,1.,500.)
+        GLU.gluLookAt(0,0,cPos,0,0,0,0,1,0)
+        SetLights()            
+            
+        GL.glMatrixMode(GL.GL_MODELVIEW)
+        GL.glLoadIdentity()
+        matRot = G2mth.Q2Mat(Q)
+        matRot = np.concatenate((np.concatenate((matRot,[[0],[0],[0]]),axis=1),[[0,0,0,1],]),axis=0)
+        GL.glMultMatrixf(matRot.T)
+        RenderUnitVectors(0.,0.,0.)
+        radius = 2.0
+        for iat,atom in enumerate(XYZ):
+            x,y,z = atom
+            color = (144,144,144)
+            RenderSphere(x,y,z,radius,color)
+        try:
+            if Page.context: Page.canvas.SetCurrent(Page.context)
+        except:
+            pass
+        Page.canvas.SwapBuffers()
+
+    def OnSize(event):
+        Draw('size')
+        
+    def OnFocus(event):
+        Draw('focus')
+        
+    def OnKeyBox(event):
+        mode = cb.GetValue()
+        if mode in ['jpeg','bmp','tiff',]:
+            try:
+                import Image as Im
+            except ImportError:
+                try:
+                    from PIL import Image as Im
+                except ImportError:
+                    print ("PIL/pillow Image module not present. Cannot save images without this")
+                    raise Exception("PIL/pillow Image module not found")
+            
+            Fname = os.path.join(Mydir,Page.name+'.'+mode)
+            print (Fname+' saved')
+            size = Page.canvas.GetSize()
+            GL.glPixelStorei(GL.GL_UNPACK_ALIGNMENT, 1)
+            Pix = GL.glReadPixels(0,0,size[0],size[1],GL.GL_RGB, GL.GL_UNSIGNED_BYTE)
+            im = Im.new("RGB", (size[0],size[1]))
+            try:
+                im.frombytes(Pix)
+            except AttributeError:
+                im.fromstring(Pix)
+            im = im.transpose(Im.FLIP_TOP_BOTTOM)
+            im.save(Fname,mode)
+            cb.SetValue(' save as/key:')
+            G2frame.G2plotNB.status.SetStatusText('Drawing saved to: '+Fname,1)
+
+    # PlotRigidBody execution starts here (N.B. initialization above)
+    new,plotNum,Page,Plot,lim = G2frame.G2plotNB.FindPlotTab('Bead model','ogl')
+    if new:
+        Page.views = False
+    Page.name = Atoms[0]
+    Page.Choice = None
+    choice = [' save as:','jpeg','tiff','bmp',]
+    cb = wx.ComboBox(G2frame.G2plotNB.status,style=wx.CB_DROPDOWN|wx.CB_READONLY,choices=choice)
+    cb.Bind(wx.EVT_COMBOBOX, OnKeyBox)
+    cb.SetValue(' save as/key:')
+    
+    Page.canvas.Bind(wx.EVT_MOUSEWHEEL, OnMouseWheel)
+    Page.canvas.Bind(wx.EVT_LEFT_DOWN, OnMouseDown)
+    Page.canvas.Bind(wx.EVT_RIGHT_DOWN, OnMouseDown)
+    Page.canvas.Bind(wx.EVT_MIDDLE_DOWN, OnMouseDown)
+    Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
+    Page.canvas.Bind(wx.EVT_SIZE, OnSize)
+    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
+    Page.camera['position'] = defaults['cameraPos']
+    Page.camera['backColor'] = np.array([0,0,0,0])
+    try:
+        if Page.context: Page.canvas.SetCurrent(Page.context)
+    except:
+        pass
+    Draw('main')
+    Draw('main')    #to fill both buffers so save works
+
+################################################################################
 #### Plot Rigid Body
 ################################################################################
@@ -9377 +9611 @@
         Draw('size')
         
+    def OnFocus(event):
+        Draw('focus')
+        
     def OnKeyBox(event):
         mode = cb.GetValue()
@@ -9422 +9659 @@
     Page.canvas.Bind(wx.EVT_MOTION, OnMouseMove)
     Page.canvas.Bind(wx.EVT_SIZE, OnSize)
+    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
     Page.camera['position'] = defaults['cameraPos']
     Page.camera['backColor'] = np.array([0,0,0,0])
@@ -9831 +10069 @@
         Draw('size')
         
+    def OnFocus(event):
+        Draw('focus')
+        
     # PlotLayers execution starts here
     cell = Layers['Cell'][1:7]
@@ -9876 +10117 @@
     Page.canvas.Bind(wx.EVT_KEY_UP, OnPlotKeyPress)
     Page.canvas.Bind(wx.EVT_SIZE, OnSize)
+    Page.canvas.Bind(wx.EVT_SET_FOCUS, OnFocus)
     Page.camera['position'] = defaults['cameraPos']
     Page.camera['backColor'] = np.array([0,0,0,0])

trunk/GSASIIpwdGUI.py

@@ -358 +358 @@
         'Pair':{'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User',
                 'Percent error':2.5,'Background':[0,False],'Distribution':[],
-                'Moore':20,'Dist G':100.,},            
+                'Moore':20,'Dist G':100.,'Result':[],},            
         'Particle':{'Matrix':{'Name':'vacuum','VolFrac':[0.0,False]},'Levels':[],},
         'Shapes':{'outName':'run','NumAA':100,'Niter':1,'AAscale':1.0,'Symm':1,'bias-z':0.0,
@@ -5224 +5224 @@
         data['BackFile'] = ''
     if 'Pair' not in data:
-        data['Pair'] = {'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User',
+        data['Pair'] = {'Method':'Regularization','MaxRadius':100.,'NBins':100,'Errors':'User','Result':[],
             'Percent error':2.5,'Background':[0,False],'Distribution':[],'Moore':[20,False],'Dist G':100.,}  
     if 'Shapes' not in data:
         data['Shapes'] = {'outName':'run','NumAA':100,'Niter':1,'AAscale':1.0,'Symm':1,'bias-z':0.0,
             'inflateV':1.0,'AAglue':0.0,'pdbOut':False}
+    plotDefaults = {'oldxy':[0.,0.],'Quaternion':[0.,0.,0.,1.],'cameraPos':150.,'viewDir':[0,0,1],}
 
     #end patches
@@ -5428 +5429 @@
       doi: xxx''',
       caption='Program Shapes',style=wx.ICON_INFORMATION)
-            Phases,Patterns,PRcalc = G2shapes.G2shapes(Profile,ProfDict,Limits,data)
-            
-            
+            data['Pair']['Result'] = []       #clear old results (if any) for now
+            data['Pair']['Result'] = G2shapes.G2shapes(Profile,ProfDict,Limits,data)
+            wx.CallAfter(UpdateModelsGrid,G2frame,data)
             
     def OnUnDo(event):
@@ -5654 +5655 @@
         def OnPDBout(event):
             data['Shapes']['pdbOut'] = not data['Shapes']['pdbOut']
+            
+        def OnShapeSelect(event):
+            r,c =  event.GetRow(),event.GetCol()
+            shapeTable.SetValue(r,c,False)
+            selAtoms = Atoms[2*r+(c-1)]
+            pattern = Patterns[r]
+            PRvals = PRcalc[r]
+            print('%s %d'%('num. beads',len(selAtoms[1])))
+            print('%s %.3f'%('selected r value',pattern[-1]))
+            print('%s %.3f'%('selected Delta P(r)',PRvals[-1]))
+            G2plt.PlotSASDPairDist(G2frame)
+            RefreshPlots(True)
+            
+            G2plt.PlotBeadModel(G2frame,selAtoms,plotDefaults)
+            
+            
         
         shapeSizer = wx.BoxSizer(wx.VERTICAL)
@@ -5693 +5710 @@
         
         shapeSizer.Add(parmSizer)
+        
+        if len(data['Pair'].get('Result',[])):
+            shapeSizer.Add(wx.StaticText(G2frame.dataWindow,label=' SHAPES run results:'),0,WACV)
+            Atoms,Patterns,PRcalc = data['Pair']['Result']
+            colLabels = ['name','show beads','show shape','Rvalue','P(r) dif','Nbeads','Nshape']
+            Types = [wg.GRID_VALUE_STRING,]+2*[wg.GRID_VALUE_BOOL,]+2*[wg.GRID_VALUE_FLOAT+':10,3',]+2*[wg.GRID_VALUE_LONG,]
+            rowLabels = [str(i) for i in range(len(Patterns))]
+            tableVals = []
+            for i in range(len(Patterns)):
+                tableVals.append([Atoms[2*i][0],False,False,Patterns[i][-1],PRcalc[i][-1],len(Atoms[2*1][1]),len(Atoms[2*i+1][1])])
+            shapeTable = G2G.Table(tableVals,rowLabels=rowLabels,colLabels=colLabels,types=Types)
+            ShapesResult = G2G.GSGrid(G2frame.dataWindow)
+            ShapesResult.SetTable(shapeTable,True)
+            ShapesResult.AutoSizeColumns(False)
+            ShapesResult.Bind(wg.EVT_GRID_CELL_CHANGED, OnShapeSelect)
+            for r in range(len(Patterns)):
+                for c in range(7):
+                    if c in [1,2]:
+                        ShapesResult.SetReadOnly(r,c,isReadOnly=False)
+                    else:
+                        ShapesResult.SetReadOnly(r,c,isReadOnly=True)
+    
+            shapeSizer.Add(ShapesResult,0,WACV)
         return shapeSizer