Changeset 961

Timestamp:

Jun 20, 2013 11:57:22 AM (8 years ago)

Author:

vondreele

Message:

fix rb constraints - useCount
fix Uij2Ueqv
mods to anneal algorithms
allow save of MC/SA results from one set of runs to another
use flat shading for polyhedra
a callafter in UpdatePDFGrid - may need more
fixes for when RBs defined but not used

Location:

trunk

Files:

• GSASIIconstrGUI.py (modified) (1 diff)
• GSASIIlattice.py (modified) (1 diff)
• GSASIImath.py (modified) (13 diffs)
• GSASIIphsGUI.py (modified) (7 diffs)
• GSASIIplot.py (modified) (2 diffs)
• GSASIIpwdGUI.py (modified) (1 diff)
• GSASIIstrMath.py (modified) (1 diff)

trunk/GSASIIconstrGUI.py

@@ -987 +987 @@
                 del data['Vector'][rbId]
                 data['RBIds']['Vector'].remove(rbId)
+                rbData['useCount'] -= 1
                 wx.CallAfter(UpdateVectorRB)
                 

trunk/GSASIIlattice.py

@@ -284 +284 @@
     U = np.inner(Amat,np.inner(U,Amat).T)
     E,R = nl.eigh(U)
-    return (E[0]+E[2]+E[5])/3.      #lower triangle?
+    return np.sum(E)/3.
         
 def CosSinAngle(U,V,G):

trunk/GSASIImath.py

@@ -2087 +2087 @@
 
 
-#  A schedule due to Lester Ingber
+#  A schedule due to Lester Ingber modified to use bounds - OK
 class fast_sa(base_schedule):
     def init(self, **options):
@@ -2101 +2101 @@
         u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
         T = self.T
-        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0
-        xc = y*(self.upper - self.lower)/2.0+self.lower
-        return xc
+        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
+        xnew = xc*(self.upper - self.lower)+self.lower
+        return xnew
 #        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
 #        xc = y*(self.upper - self.lower)
@@ -2114 +2114 @@
         return
 
-class cauchy_sa(base_schedule):
-#    def update_guess(self, x0):
-#        x0 = asarray(x0)
+class cauchy_sa(base_schedule):     #modified to use bounds - not good
+    def update_guess(self, x0):
+        x0 = asarray(x0)
+        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
+        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
+        xnew = xc*(self.upper - self.lower)+self.lower
+        return xnew
 #        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
 #        xc = self.learn_rate * self.T * tan(numbers)
@@ -2141 +2145 @@
         return
 
-class log_sa(base_schedule):
+class log_sa(base_schedule):        #OK
 
     def init(self,**options):
@@ -2361 +2365 @@
             if abs(af[-1]-best_state.cost) > feps*10:
                 retval = 5
-                print "Warning: Cooled to %f at %s but this is not" \
-                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
-                      + " the smallest point found."
+#                print "Warning: Cooled to %f at %s but this is not" \
+#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
+#                      + " the smallest point found."
             break
         if (Tf is not None) and (schedule.T < Tf):
@@ -2395 +2399 @@
     '''
     gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
+    
+    twopi = 2.0*np.pi
+    global tsum
+    tsum = 0.
     
     def getMDparms(item,pfx,parmDict,varyList):
@@ -2463 +2471 @@
     def GetAtomM(Xdata,SGData):
         Mdata = []
-        for xyz in Xdata.T:
+        for xyz in Xdata:
             Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
         return np.array(Mdata)
@@ -2485 +2493 @@
             if parmDict[pfx+'Type'] in ['Vector','Residue']:
                 if parmDict[pfx+'Type'] == 'Vector':
-                    RBId = parmDict[pfx+'RBId']
-                    RBRes = RBdata['Vector'][RBId]
+                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
                     aTypes = RBRes['rbTypes']
                     vecs = RBRes['rbVect']
@@ -2494 +2501 @@
                         Cart += vec*mag
                 elif parmDict[pfx+'Type'] == 'Residue':
-                    RBId = parmDict[pfx+'RBId']
-                    RBRes = RBdata['Residue'][RBId]
+                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
                     aTypes = RBRes['rbTypes']
                     Cart = np.array(RBRes['rbXYZ'])
@@ -2519 +2525 @@
                     if parm in parmDict:
                         keys[key][atNo] = parmDict[parm]
+                iatm += 1
             else:
                 continue        #skips March Dollase
-        return Tdata,Xdata
+        return Tdata,Xdata.T
     
     def calcMDcorr(MDval,MDaxis,Uniq,G):
@@ -2537 +2544 @@
 #            ParmDict:
 #        puts result F^2 in each ref[8] in refList
-#        '''        
-        twopi = 2.0*np.pi
+#        '''       
+        global tsum
         parmDict.update(dict(zip(varyList,values)))
         Tdata,Xdata = GetAtomTX(RBdata,parmDict)
@@ -2548 +2555 @@
         Srefs = np.zeros(len(refList))
         sumFcsq = 0.
+        t0 = time.time()
         for refl in refList:
-            fbs = 0.
-            H = refl[:3]
-            for i,El in enumerate(Tdata):
+            for i,El in enumerate(Tdata):       #NB: generator here is slower!
                 FF[i] = refl[7][El]
-            Uniq = refl[8]
-            phi = refl[9]
-            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
-            sinp = np.sin(phase)
-            cosp = np.cos(phase)
-            occ = Mdata/len(Uniq)
-            fa = np.array(FF*occ*cosp)
-            fas = np.sum(fa)
-            if not ifInv:
-                fb = np.array(FF*occ*sinp)
-                fbs = np.sum(fb)
-            fcsq = (fas**2+fbs**2)*refl[3]      #*calcMDcorr(MDval,MDaxis,Uniq,Gmat)
-            sumFcsq += fcsq
-            refl[5] = fcsq
+            FF *= Mdata/len(refl[8])            #FF*occ
+            phase = np.inner(refl[8],Xdata)     #hx+ky+lz
+            phase += refl[9][:,np.newaxis]      #hx+ky+lz+phi
+            cosp = np.cos(twopi*phase)
+            fas = np.sum(FF*cosp)
+            if ifInv:
+                fbs = 0.
+            else:
+                sinp = np.sin(twopi*phase)
+                fbs = np.sum(FF*sinp)
+            refl[5] = (fas**2+fbs**2)*refl[3]  #*calcMDcorr(MDval,MDaxis,Uniq,Gmat)
+            sumFcsq += refl[5]
+        tsum += (time.time()-t0)
         scale = (parmDict['sumFosq']/sumFcsq)
         for iref,refl in enumerate(refList):
@@ -2700 +2705 @@
         quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
         lower=lower, upper=upper, slope=MCSA['log slope'],dlg=pgbar)
-    return [False,results[1],results[2],results[0],varyList]
+    Result = [False,False,results[1],results[2],]+list(results[0])
+    Result.append(varyList)
+    return Result,tsum
 
         

trunk/GSASIIphsGUI.py

@@ -4208 +4208 @@
                     resultsGrid.ForceRefresh()
                     result = Results[r]
-                    for key,val in zip(result[4],result[3]):
+                    for key,val in zip(result[-1],result[4:-1]):
                         vals = key.split(':')
                         nObj,name = int(vals[0]),vals[1]
@@ -4234 +4234 @@
                     wx.CallAfter(UpdateMCSA)
                     G2plt.PlotStructure(G2frame,data)
+                elif c == 1:
+                    if Results[r][1]:
+                        Results[r][1] = False
+                    else:
+                        Results[r][1] = True
+                    resultsTable.SetValue(r,c,Results[r][1])
+                    resultsGrid.ForceRefresh()
+                        
                 
             resultsSizer = wx.BoxSizer(wx.VERTICAL)
@@ -4239 +4247 @@
             resultVals = []
             for result in Results:
-                maxVary = max(maxVary,len(result[3]))
-                resultVals.append(result[:3]+list(result[3]))
+                maxVary = max(maxVary,len(result[-1]))
+                resultVals.append(result[:-1])
             rowLabels = []
             for i in range(len(Results)): rowLabels.append(str(i))
-            colLabels = ['Select','Residual','Tmin',]
-            for item in result[4]: colLabels.append(item)
-#            for i in range(maxVary): colLabels.append('variable:'+str(i))
-            Types = [wg.GRID_VALUE_BOOL,wg.GRID_VALUE_FLOAT+':10,4',
+            colLabels = ['Select','Keep','Residual','Tmin',]
+            for item in result[-1]: colLabels.append(item)   #from last result from for loop above
+            Types = [wg.GRID_VALUE_BOOL,wg.GRID_VALUE_BOOL,wg.GRID_VALUE_FLOAT+':10,4',
                 wg.GRID_VALUE_FLOAT+':10,4',]+maxVary*[wg.GRID_VALUE_FLOAT+':10,5',]
             resultsTable = G2gd.Table(resultVals,rowLabels=rowLabels,colLabels=colLabels,types=Types)
@@ -4255 +4262 @@
             for r in range(resultsGrid.GetNumberRows()):
                 for c in range(resultsGrid.GetNumberCols()):
-                    if c == 0:
+                    if c in [0,1]:
                         resultsGrid.SetReadOnly(r,c,isReadOnly=False)
                     else:
@@ -4307 +4314 @@
         if G2frame.dataFrame.PhaseUserSize is None:
             mainSizer.FitInside(G2frame.dataFrame)
-            Size = mainSizer.GetMinSize()
+            Size = mainSizer.Fit()
             Size[0] += 40
-            Size[1] = max(Size[1],290) + 35
+            Size[1] = max(Size[1],350) + 35
             MCSA.SetSize(Size)
             MCSA.SetScrollbars(10,10,Size[0]/10-4,Size[1]/10-1)
@@ -4326 +4333 @@
         phaseName = generalData['Name']
         MCSAdata = data['MCSA']
+        saveResult = []
+        for result in MCSAdata['Results']:
+            if result[1]:       #keep?
+                saveResult.append(result)
+        MCSAdata['Results'] = saveResult           
         covData = {}
         if 'PWDR' in reflName:
@@ -4356 +4368 @@
         pgbar.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
         pgbar.SetSize(Size)
+        time1 = time.time()
         try:
+            tsf = 0.
             for i in range(mcsaControls['Cycles']):
-                MCSAdata['Results'].append(G2mth.mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar))
-                print ' MC/SA runs completed: ',i
+                Result,tsum = G2mth.mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar)
+                MCSAdata['Results'].append(Result)
+                print ' MC/SA runs completed: %d residual: %.3f%%'%(i,100*Result[2])
+                tsf += tsum
+            print ' MC/SA run time: %.2f'%(time.time()-time1)
+            print ' Structure factor time: %.2f'%(tsf)
         finally:
             pgbar.Destroy()
+        MCSAdata['Results'] = G2mth.sortArray(MCSAdata['Results'],2,reverse=False)
         UpdateMCSA()
         G2plt.PlotStructure(G2frame,data)

trunk/GSASIIplot.py

@@ -3247 +3247 @@
         
     def RenderPolyhedra(x,y,z,Faces,color):
+        glShadeModel(GL_FLAT)
         glPushMatrix()
         glTranslate(x,y,z)
@@ -3261 +3262 @@
             glEnd()
         glPopMatrix()
+        glShadeModel(GL_SMOOTH)
 
     def RenderMapPeak(x,y,z,color,den):

trunk/GSASIIpwdGUI.py

@@ -2116 +2116 @@
             data['Form Vol'] = max(10.0,SumElementVolumes())
             formVol.SetValue('%.2f'%(data['Form Vol']))
-            UpdatePDFGrid(G2frame,data)
+            wx.CallAfter(UpdatePDFGrid,G2frame,data)
             auxPlot = ComputePDF(data)
             G2plt.PlotISFG(G2frame,newPlot=True)        

trunk/GSASIIstrMath.py

@@ -58 +58 @@
         VRBData = RBData['Vector']
         for i,rbId in enumerate(VRBIds):
-            for j in range(len(VRBData[rbId]['VectMag'])):
-                name = '::RBV;'+str(j)+':'+str(i)
-                VRBData[rbId]['VectMag'][j] = parmDict[name]
+            if VRBData[rbId]['useCount']:
+                for j in range(len(VRBData[rbId]['VectMag'])):
+                    name = '::RBV;'+str(j)+':'+str(i)
+                    VRBData[rbId]['VectMag'][j] = parmDict[name]
     for phase in Phases:
         Phase = Phases[phase]