Changeset 4578 for trunk

Timestamp:

Sep 28, 2020 1:23:10 PM (3 years ago)

Author:

toby

Message:

docs for freezing parameters started + docs cleanup; start scriptable for freezing params; record initial chi2; Show more post refinement info; noted but unfixed bkg GUI bug

Location:

trunk

Files:

• GSASIIdataGUI.py (modified) (1 diff)
• GSASIImath.py (modified) (8 diffs)
• GSASIIobj.py (modified) (1 diff)
• GSASIIplot.py (modified) (3 diffs)
• GSASIIpwdGUI.py (modified) (3 diffs)
• GSASIIscriptable.py (modified) (3 diffs)
• GSASIIstrMain.py (modified) (5 diffs)
• atmdata.py (modified) (1 diff)
• defaultIparms.py (modified) (1 diff)

trunk/GSASIIdataGUI.py

@@ -5037 +5037 @@
             lamMax = Msg.get('lamMax',0.001)
             lst = os.path.splitext(os.path.abspath(self.GSASprojectfile))[0]
-            text = 'Detailed results are in ' + lst + '.lst\n\nLoad new result?'
+            text = 'Detailed results are in ' + lst + '.lst\n'
+            if 'GOF0' in Msg and 'GOF' in Msg:
+                text += '\nReduced Chi^2 before refinement={:.3f} and after={:.3f}\n'.format(
+                    Msg['GOF0']**2,Msg['GOF']**2)
+            if 'Max shft/sig' in Msg:
+                text += '\nMax shift/sigma={:.3f}\n'.format(Msg['Max shft/sig'])
+            if 'msg' in Msg: text += '\n' + Msg['msg'] + '\n'
+            text += '\nLoad new result?'
             if lamMax >= 10.:
                 text += '\nWARNING: Steepest descents dominates;'+   \

trunk/GSASIImath.py

@@ -166 +166 @@
     Lam = np.zeros((n,n))
     Xvec = np.zeros(len(x0))
+    chisq00 = None
     while icycle < maxcyc:
         time0 = time.time()
@@ -172 +173 @@
         nfev += 1
         chisq0 = np.sum(M**2)
+        if chisq00 is None: chisq00 = chisq0
         Yvec,Amat = Hess(x0,*args)
         Adiag = np.sqrt(np.diag(Amat))
@@ -234 +236 @@
         Bmat = Bmat/Anorm
         return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],
-            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,'Xvec':Xvec}]
+            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,'Xvec':Xvec, 'chisq0':chisq00}]
     except nl.LinAlgError:
         G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error')
@@ -240 +242 @@
         if maxcyc:
             psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
-        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1,'Xvec':None}]          
+        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1,'Xvec':None, 'chisq0':chisq00}]
             
 def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
@@ -298 +300 @@
     if Print:
         G2fil.G2Print(' Hessian SVD refinement on %d variables:'%(n))
+    chisq00 = None
     while icycle < maxcyc:
         time0 = time.time()
@@ -303 +306 @@
         nfev += 1
         chisq0 = np.sum(M**2)
+        if chisq00 is None: chisq00 = chisq0
         Yvec,Amat = Hess(x0,*args)
         Adiag = np.sqrt(np.diag(Amat))
@@ -347 +351 @@
         Bmat = Bmat/Anorm
         return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':[],
-            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
+            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,
+                             'chisq0':chisq00}]
     except nl.LinAlgError:
         G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error')
@@ -353 +358 @@
         if maxcyc:
             psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
-        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]          
+        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1,
+                             'chisq0':chisq00}]
             
 def getVCov(varyNames,varyList,covMatrix):

trunk/GSASIIobj.py

@@ -1141 +1141 @@
 
 
+.. _ParameterLimits:
+
+.. index::
+   single: Parameter limits
+
+Parameter Limits
+------------------------------
+
+One of the most often requested "enhancements" for GSAS-II would be the inclusion 
+of constraints to force parameters such as occupancies or Uiso values to stay within 
+expected ranges. While it is possible for users to supply their own restraints that would
+perform this by supplying an appropriate expression with the "General" restraints, the 
+GSAS-II authors do not feel that use of restraints or constraints are a good solution for
+this common problem where parameters refine to non-physical values. This is because when 
+this occurs, most likely one of the following cases is occurring: 
+
+#. there is a significant problem 
+   with the model, for example for an x-ray fit if an O atom is placed where a S is actually 
+   present, the Uiso will refine artificially small or the occupancy much larger than unity 
+   to try to compensate for the missing electrons; or
+ 
+#. the data are simply insensitive 
+   to the parameter or combination of parameters, for example unless very high-Q data 
+   are included, the effects of a occupancy and Uiso value can have compensating effects,
+   so an assumption must be made; likewise, with neutron data natural-abundance V atoms 
+   are nearly invisible due to weak coherent scattering. No parameters can be fit for a 
+   V atom with neutrons.
+
+#. the parameter is non-physical (such as a negative Uiso value) but within 
+   two sigma (sigma = standard uncertainty, aka e.s.d.) of a reasonable value, 
+   in which case the
+   value is not problematic as it is experimentally indistinguishable from an 
+   expected value.  
+
+#. there is a systematic problem with the data (experimental error)
+
+In all these cases, this situation needs to be reviewed by a crystallographer to decide 
+how to best determine a structural model for these data. An implementation with a constraint
+or restraint is likely to simply hide the problem from the user, making it more probable 
+that a poor model choice is obtained. 
+
+What GSAS-II does implement is to allow users to specify ranges for parameters 
+that works by disabling 
+refinement of parameters that refine beyond either a lower limit or an upper limit, where 
+either or both may be optionally specified. Parameters limits are specified in the Controls 
+tree entry in dicts named as ``Controls['parmMaxDict']`` and ``Controls['parmMinDict']``, where 
+the keys are named using the standard GSAS-II variable 
+(see :func:`getVarDescr` and :func:`CompileVarDesc`) or with a variable name, where a 
+wildcard ('*') is used for histogram number or atom number (phase number is intentionally not 
+allowed as a wildcard as it makes little sense to group together different phases). Note
+that func:`prmLookup` is used to see if a name matches a wildcard. The upper or lower limit
+is placed into these dicts as a float value. These values can be edited using the window 
+created by the Calculate/"View LS parms" menu command or in scripting with the 
+:meth:`GSASIIscriptable.G2Project.set_Controls` function. 
+In the GUI, a checkbox labeled "match all histograms/atoms" is used to insert a wildcard
+into the appropriate part of the variable name.
+
+When a refinement is conducted, routine :func:`GSASIIstrMain.dropOOBvars` is used to 
+find parameters that have refined to values outside their limits. If this occurs, the parameter
+is set to the limiting value and the variable name is added to a list of frozen variables kept in the 
+``Controls['parmFrozen']`` dict. In a sequential refinement, this is kept as a list in 
+``Controls['parmFrozen'][histogram]`` (where the key is the histogram name) or in 
+``Controls['parmFrozen']['FrozenList']`` for a non-sequential fit. 
+This allows different variables
+to be frozen in each section of a sequential fit. 
+Frozen parameters are not included in refinements through removal from the 
+list of parameters to be refined (``varyList``) in :func:`GSASIIstrMain.Refine` or 
+:func:`GSASIIstrMain.SeqRefine`. Once a variable is frozen, it will not be refined in any 
+future refinements unless the the variable is removed (manually) from the list. This can also 
+be done with the Calculate/"View LS parms" menu window or ...
+
 *Classes and routines*
 ----------------------

trunk/GSASIIplot.py

@@ -209 +209 @@
 timeDebug = GSASIIpath.GetConfigValue('Show_timing',False)
 obsInCaption = True # include the observed, calc,... items in the plot caption (PlotPatterns)
-mplv = eval(mpl.__version__.replace('.',','))
 
 #matplotlib 2.0.x dumbed down Paired to 16 colors - 
@@ -333 +332 @@
         
     def ToolBarDraw(self):
+        mplv = eval(mpl.__version__.replace('.',','))
         if mplv[0] >= 3 and mplv[1] >= 3:
             self.toolbar.canvas.draw_idle()
@@ -10575 +10575 @@
 
 def SetupLegendPick(legend,new,delay=5):
+    mplv = eval(mpl.__version__.replace('.',','))
     legend.delay = delay*1000 # Hold time in ms for clear; 0 == forever
     for line in legend.get_lines():

trunk/GSASIIpwdGUI.py

@@ -1520 +1520 @@
 
         def OnPeaks(event):
+            'Respond to a change in the number of background peaks'
             data[1]['nPeaks'] = int(peaks.GetValue())
             M = len(data[1]['peaksList'])
@@ -1531 +1532 @@
             if N == 0:
                 CalcBack(G2frame.PatternId)
-                G2plt.PlotPatterns(G2frame,plotType='PWDR')                
+                G2plt.PlotPatterns(G2frame,plotType='PWDR')
+            # this callback is crashing wx when there is an open
+            # peaksGrid cell editor, at least on Mac. Code below
+            # should fix this, but it does not.
+            # https://stackoverflow.com/questions/64082199/wxpython-grid-destroy-with-open-celleditor-crashes-python-even-with-disablece
+            if peaksGrid and peaksGrid.IsCellEditControlEnabled():
+                # complete any grid edits in progress
+                #print('closing')
+                peaksGrid.HideCellEditControl()
+                peaksGrid.DisableCellEditControl()
+                #wx.CallLater(100,peaksGrid.Destroy) # crashes python
             wx.CallAfter(UpdateBackground,G2frame,data)
             
@@ -1567 +1578 @@
         peaksSizer.Add(topSizer)
         G2frame.dataWindow.currentGrids = []
+        peaksGrid = None
         if data[1]['nPeaks']:
             peaksSizer.Add(wx.StaticText(G2frame.dataWindow,-1,' Peak list:'),0,WACV)       

trunk/GSASIIscriptable.py

@@ -3360 +3360 @@
             raise Exception('{} is an invalid control value'.format(control))
         
-    def set_Controls(self, control, value):
+    def set_Controls(self, control, value, variable=None):
         '''Set project controls
 
         :param str control: the item to be set. See below for allowed values. 
         :param value: the value to be set.
+        :param str variable: used only with control set to "parmMin" or "parmMax"
 
         Allowed values for parameter control:
@@ -3378 +3379 @@
             * Reverse Seq: when True, sequential refinement is performed on the
               reversed list of histograms.
+            * parmMin & parmMax: set a maximum or minimum value for a refined 
+              parameter. Note that variable will be a GSAS-II variable name, 
+              optionally with * specified for a histogram or atom number and
+              value must be a float. 
+              (See :ref:`Parameter Limits<ParameterLimits>` description.)
 
         .. seealso::
             :meth:`get_Controls`
         '''
+        if control.startswith('parmM'):
+            if not variable:
+                raise Exception('set_Controls requires a value for variable for control=parmMin or parmMax')
+            for key in ('parmMinDict','parmMaxDict','parmFrozen'):
+                if key not in self.data['Controls']['data']: self.data['Controls']['data'][key] = {}
         if control == 'cycles':
             self.data['Controls']['data']['max cyc'] = int(value)
@@ -3398 +3409 @@
                                         .format(i,j))
             self.data['Controls']['data']['Seq Data'] = histlist
+        elif control == 'parmMin' or control == 'parmMax':
+            self.data['Controls']['data'][control+'Dict'][variable] = float(value)
         else:
             raise Exception('{} is an invalid control value'.format(control))

trunk/GSASIIstrMain.py

@@ -61 +61 @@
         ifPrint = False
     Rvals = {}
+    chisq0 = None
     while True:
         begin = time.time()
@@ -91 +92 @@
             Rvals['lamMax'] = result[2]['lamMax']
             Controls['Marquardt'] = -3  #reset to default
+            if 'chisq0' in result[2] and chisq0 is None:
+                chisq0 = result[2]['chisq0']
         elif 'Hessian SVD' in Controls['deriv type']:
             maxCyc = Controls['max cyc']
@@ -102 +105 @@
                 break
             ncyc = result[2]['num cyc']+1
+            if 'chisq0' in result[2] and chisq0 is None:
+                chisq0 = result[2]['chisq0']
         else:           #'numeric'
             result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
@@ -182 +187 @@
     if IfOK:
         G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
+    if chisq0 is not None:
+        Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
     return IfOK,Rvals,result,covMatrix,sig
 
@@ -293 +300 @@
             G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
             G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
-            if len(frozen) > 0:
-                G2fil.G2Print(
-                    ' {} variables were outside limits and were frozen (now {} frozen total)\n'
-                    .format(len(frozen),len(parmFrozenList)))
+            if len(frozen):
+                msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
+                    .format(len(frozen),len(parmFrozenList))
+                    )
+                G2fil.G2Print(msg)
+                Rvals['msg'] = msg
+            elif len(parmFrozenList):
+                msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
+                    .format(len(parmFrozenList))
+                    )
+                G2fil.G2Print('Note: ',msg)
+                Rvals['msg'] = msg
             G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
             printFile.close()

trunk/atmdata.py

@@ -1190 +1190 @@
         'N':[    1.8845,  9.161, 2.0746,  4.042,  -.1318,  1.723,  .0000,   .000,  .0020, 2.00]},
 }
-'''
-unused j<0>
+
+pass # placed here so that the following is not processed by Sphinx as a doc string
+''' # unused unused terms follow
+j<0>
 W_6s15d5    FFj0A=0.3811 FFj0a=62.707 FFj0B=0.7523 FFj0b=21.434 FFj0C=12.5449 FFj0c=2.702 FFj0D=12.4130 FFj0d=2.674 FFj0E=0.0023 0.0365
 W_6s25d4    FFj0A=0.3653 FFj0a=53.965 FFj0B=0.7926 FFj0b=20.078 FFj0C=0.8142 FFj0c=3.030 FFj0D=0.6581 FFj0d=2.476 FFj0E=0.0023 0.0247

trunk/defaultIparms.py

@@ -34 +34 @@
 defaultIparm_lbl = []
 defaultIparms = []
+
+pass # placed here so that the following is not processed by Sphinx as a doc string
 ''' the phrase 'lab data' in the following makes the default goniometer type='Bragg-Brentano'
     otherwise it is set to 'Debye-Scherrer'