Changeset 5274

Timestamp:

May 12, 2022 2:57:13 PM (17 months ago)

Author:

toby

Message:

implement save of intensities by phase (partials)

Location:

trunk

Files:

• GSASIIdataGUI.py (modified) (4 diffs)
• GSASIIstrMath.py (modified) (5 diffs)
• help/gsasII.html (modified) (2 diffs)

trunk/GSASIIdataGUI.py

@@ -767 +767 @@
         self.Refine.append(item)
         self.Bind(wx.EVT_MENU, self.OnRefine, id=item.GetId())
-
+        # this could be disabled during sequential fits:
+        item = parent.Append(wx.ID_ANY,'Compute partials','Record the contribution from each phase')
+        self.Bind(wx.EVT_MENU, self.OnRefinePartials, id=item.GetId())
         item = parent.Append(wx.ID_ANY,'&Run Fprime','X-ray resonant scattering')
         self.Bind(wx.EVT_MENU, self.OnRunFprime, id=item.GetId())
@@ -5350 +5352 @@
         Called from the Calculate/Refine menu.
         '''
+        self._cleanPartials()  # phase partials invalid after a refinement
         if self.testSeqRefineMode():
             self.OnSeqRefine(event)
@@ -5448 +5451 @@
         either single or sequentially
         '''
+        self._cleanPartials()  # phase partials invalid after a refinement
         self.OnFileSave(event)
         item = GetGPXtreeItemId(self,self.root,'Covariance')
@@ -5502 +5506 @@
             return True
         return False
-
+    
+    def _cleanPartials(self):
+        '''Delete any partials created with :meth:`OnRefinePartials`; used
+        in GUI-based refinements, as the partials are no longer correct after 
+        any fit.
+        '''
+        PhasePartials = os.path.abspath(os.path.splitext(self.GSASprojectfile)[0])
+        histograms,phases = self.GetUsedHistogramsAndPhasesfromTree()
+        for h in histograms:
+            if 'PWDR' not in h[:4]: continue
+            hId = histograms[h]['hId']
+            phPartialFile = PhasePartials+'.partials'+str(hId)
+            if os.path.exists(phPartialFile):
+                os.remove(phPartialFile)
+                print('file deleted:',phPartialFile)
+
+
+    def OnRefinePartials(self,event):
+        '''Computes and saves the intensities from each phase for each powder 
+        histogram. These inten
+        Do a 0 cycle fit with no variables to pickle intensities for each 
+        phase into a file.
+        Not for sequential fits.
+        Sets Controls['PhasePartials'] to a file name to trigger save of
+        info in :meth:`GSASIIstrMath.getPowderProfile` and then clear that. 
+        '''
+        if self.testSeqRefineMode():
+            G2G.G2MessageBox(self,
+                'Phase partials cannot be computed for seqiential fits'
+                'not allowed')
+            return
+        Controls = self.GPXtree.GetItemPyData(GetGPXtreeItemId(self,self.root, 'Controls'))
+        savCyc,Controls['max cyc'] = Controls['max cyc'],0
+        Controls['PhasePartials'] = os.path.abspath(os.path.splitext(self.GSASprojectfile)[0])
+        self.OnFileSave(event)
+        dlg = G2G.RefinementProgress(parent=self)
+        self.SaveTreeSetting() # save the current tree selection
+        self.GPXtree.SaveExposedItems()             # save the exposed/hidden tree items
+        
+        import pickle
+        
+        try:
+            OK,Rvals = G2stMn.Refine(self.GSASprojectfile,dlg,refPlotUpdate=None,newLeBail=False)
+        except Exception as msg:
+            print('Refinement failed with message',msg)
+            return True
+        finally:
+            dlg.Update(101.) # forces the Auto_Hide; needed after move w/Win & wx3.0
+            # assemble the arrays by histogram
+            
+        histograms,phases = self.GetUsedHistogramsAndPhasesfromTree()
+        for h in histograms:
+            if 'PWDR' not in h[:4]: continue
+            hId = histograms[h]['hId']
+            hfx = ':%d:'%(hId)
+            Limits = histograms[h]['Limits'][1]
+            x = histograms[h]['Data'][0]
+            xB = np.searchsorted(x,Limits[0])
+            xF = np.searchsorted(x,Limits[1])+1
+            phPartialFile = Controls['PhasePartials']+'.partials'+str(hId)
+            fp = open(phPartialFile,'rb')
+            pickle.load(fp)   # skip over x values
+            valList = [x,histograms[h]['Data'][1],
+                           histograms[h]['Data'][2],histograms[h]['Data'][3],
+                           histograms[h]['Data'][4]]
+            lblList = ['x','obs','wgt','calc','background']
+            while True:
+                try:
+                    ph = pickle.load(fp)
+                    lblList.append(ph)
+                    ypartial = np.zeros_like(x)
+                    ypartial[xB:xF] = pickle.load(fp)
+                    valList.append(ypartial)
+                except EOFError:
+                    break
+            fp.close()
+            phPartialFile = Controls['PhasePartials']+'_part_'+str(histograms[h]['hId'])+'.csv'
+            fp = open(phPartialFile,'w')
+            fp.write(', '.join(lblList))
+            fp.write('\n')
+            for l in zip(*valList):
+                fp.write(', '.join([str(i) for i in l])+'\n')
+            fp.close()
+            print('File',phPartialFile,'written')
+        dlg.Destroy()
+        Controls['max cyc'] = savCyc
+        Controls['PhasePartials'] = None
+        self.OnFileSave(event)
+        return False
+    
     def reloadFromGPX(self,rtext=None):
         '''Deletes current data tree & reloads it from GPX file (after a 

trunk/GSASIIstrMath.py

@@ -19 +19 @@
 import scipy.stats as st
 import multiprocessing as mp
+import pickle
 import GSASIIpath
 GSASIIpath.SetVersionNumber("$Revision$")
@@ -3274 +3275 @@
         return alp,bet
 
+    # set up for save of phase partials if triggered in GSASIIdataGUI.OnRefinePartials
+    phasePartials = calcControls.get('PhasePartials',None)
+    def SavePartial(*args): pass
+    if phasePartials:
+        phPartialFile = phasePartials+'.partials'+str(Histogram['hId'])
+        phPartialFP = open(phPartialFile,'wb')  # create/overwrite a file
+        pickle.dump(x,phPartialFP)
+        phPartialFP.close()
+        print('Storing intensity by phase in',phPartialFile)
+        def SavePartial(phase,y):
+            phPartialFP = open(phPartialFile,'ab')  # append to file
+            pickle.dump(phase,phPartialFP)
+            pickle.dump(y,phPartialFP)
+            phPartialFP.close()
     hId = Histogram['hId']
     hfx = ':%d:'%(hId)
@@ -3345 +3360 @@
         useMP,ncores = G2mp.InitMP()
         if len(refDict['RefList']) < 100: useMP = False        
+        if phasePartials:
+            useMP = False
+            ypartial = np.zeros_like(yb)
         if useMP: # multiprocessing: create a set of initialized Python processes
             MPpool = mp.Pool(ncores,G2mp.InitPwdrProfGlobals,[im,shl,x])
@@ -3387 +3405 @@
                     fp = G2pwd.getFCJVoigt3(refl[5+im],refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))[0]
                     yc[iBeg:iFin] += refl[11+im]*refl[9+im]*fp   #>90% of time spent here
+                    if phasePartials: ypartial[iBeg:iFin] += refl[11+im]*refl[9+im]*fp
                 if Ka2:
                     pos2 = refl[5+im]+lamRatio*tand(refl[5+im]/2.0)       # + 360/pi * Dlam/lam * tan(th)
@@ -3403 +3422 @@
                         fp2 = G2pwd.getFCJVoigt3(pos2,refl[6+im],refl[7+im],shl,ma.getdata(x[iBeg:iFin]))[0]
                         yc[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*fp2       #and here
+                        if phasePartials: ypartial[iBeg:iFin] += refl[11+im]*refl[9+im]*kRatio*fp2
+            if phasePartials: SavePartial(phase,ypartial)
         elif 'E' in calcControls[hfx+'histType']:
             

trunk/help/gsasII.html

@@ -5525 +5525 @@
 style='font-family:Symbol;mso-fareast-font-family:Symbol;mso-bidi-font-family:
 Symbol'><span style='mso-list:Ignore'>·<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
+</span></span></span><![endif]><b>Compute partials - </b>This runs a
+zero-cycle refinement where the contributions from each phase (phase
+partial intensities) are written separately to a file for each histogram in the
+refinement. This information is written into files named
+<tt><I>project</I>.partials</tt><I>N</I> where <I>N</I> is the
+histogram number and <tt><I>project</I></tt> is the GSAS-II project (.gpx)
+name. From those file(s), files named
+<tt><I>project</I>_part_</tt><I>N</I><tt>.csv</tt> are also created.
+The first set of files are intended for future use internally by
+GSAS-II and the latter files are for use in external software. The
+former set will be deleted if a refinement is run. Use this menu
+command to recreate them if needed. 
+</P>
+
+<p class=MsoListParagraphCxSpMiddle style='text-indent:-.25in;mso-list:l2 level1 lfo5'><![if !supportLists]><span
+style='font-family:Symbol;mso-fareast-font-family:Symbol;mso-bidi-font-family:
+Symbol'><span style='mso-list:Ignore'>·<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
 </span></span></span><![endif]><b>Run <span class=SpellE>Fprime</span> - </b>This
 run the utility routine <span class=SpellE>Fprime</span> that displays real and
@@ -9746 +9763 @@
 </span></div>
 
-<p class=MsoNormal style='tab-stops:45.8pt 91.6pt 137.4pt 183.2pt 229.0pt 274.8pt 320.6pt 366.4pt 412.2pt 458.0pt 503.8pt 549.6pt 595.4pt 641.2pt 687.0pt 732.8pt'><!-- hhmts start -->Last
-modified: Wed Feb 23 16:21:54 CST 2022 <o:p></o:p></p>
-
-<!-- hhmts end --></div>
+<p class=MsoNormal style='tab-stops:45.8pt 91.6pt 137.4pt 183.2pt 229.0pt 274.8pt 320.6pt 366.4pt 412.2pt 458.0pt 503.8pt 549.6pt 595.4pt 641.2pt 687.0pt 732.8pt'><!-- hhmts start -->Last modified: Thu May 12 14:50:38 CDT 2022 <!-- hhmts end --></div>
 
 </body>