Changeset 4578


Ignore:
Timestamp:
Sep 28, 2020 1:23:10 PM (12 months 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:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r4573 r4578  
    50375037            lamMax = Msg.get('lamMax',0.001)
    50385038            lst = os.path.splitext(os.path.abspath(self.GSASprojectfile))[0]
    5039             text = 'Detailed results are in ' + lst + '.lst\n\nLoad new result?'
     5039            text = 'Detailed results are in ' + lst + '.lst\n'
     5040            if 'GOF0' in Msg and 'GOF' in Msg:
     5041                text += '\nReduced Chi^2 before refinement={:.3f} and after={:.3f}\n'.format(
     5042                    Msg['GOF0']**2,Msg['GOF']**2)
     5043            if 'Max shft/sig' in Msg:
     5044                text += '\nMax shift/sigma={:.3f}\n'.format(Msg['Max shft/sig'])
     5045            if 'msg' in Msg: text += '\n' + Msg['msg'] + '\n'
     5046            text += '\nLoad new result?'
    50405047            if lamMax >= 10.:
    50415048                text += '\nWARNING: Steepest descents dominates;'+   \
  • trunk/GSASIImath.py

    r4560 r4578  
    166166    Lam = np.zeros((n,n))
    167167    Xvec = np.zeros(len(x0))
     168    chisq00 = None
    168169    while icycle < maxcyc:
    169170        time0 = time.time()
     
    172173        nfev += 1
    173174        chisq0 = np.sum(M**2)
     175        if chisq00 is None: chisq00 = chisq0
    174176        Yvec,Amat = Hess(x0,*args)
    175177        Adiag = np.sqrt(np.diag(Amat))
     
    234236        Bmat = Bmat/Anorm
    235237        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[],
    236             'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,'Xvec':Xvec}]
     238            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,'Xvec':Xvec, 'chisq0':chisq00}]
    237239    except nl.LinAlgError:
    238240        G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error')
     
    240242        if maxcyc:
    241243            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    242         return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1,'Xvec':None}]         
     244        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing,'SVD0':-1,'Xvec':None, 'chisq0':chisq00}]
    243245           
    244246def HessianSVD(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.e-6, maxcyc=0,lamda=-3,Print=False,refPlotUpdate=None):
     
    298300    if Print:
    299301        G2fil.G2Print(' Hessian SVD refinement on %d variables:'%(n))
     302    chisq00 = None
    300303    while icycle < maxcyc:
    301304        time0 = time.time()
     
    303306        nfev += 1
    304307        chisq0 = np.sum(M**2)
     308        if chisq00 is None: chisq00 = chisq0
    305309        Yvec,Amat = Hess(x0,*args)
    306310        Adiag = np.sqrt(np.diag(Amat))
     
    347351        Bmat = Bmat/Anorm
    348352        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':[],
    349             'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2}]
     353            'SVD0':Nzero,'Converged': ifConverged, 'DelChi2':deltaChi2,
     354                             'chisq0':chisq00}]
    350355    except nl.LinAlgError:
    351356        G2fil.G2Print('ouch #2 linear algebra error in making v-cov matrix', mode='error')
     
    353358        if maxcyc:
    354359            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
    355         return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1}]         
     360        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':0.,'psing':psing,'SVD0':-1,
     361                             'chisq0':chisq00}]
    356362           
    357363def getVCov(varyNames,varyList,covMatrix):
  • trunk/GSASIIobj.py

    r4534 r4578  
    11411141
    11421142
     1143.. _ParameterLimits:
     1144
     1145.. index::
     1146   single: Parameter limits
     1147
     1148Parameter Limits
     1149------------------------------
     1150
     1151One of the most often requested "enhancements" for GSAS-II would be the inclusion
     1152of constraints to force parameters such as occupancies or Uiso values to stay within
     1153expected ranges. While it is possible for users to supply their own restraints that would
     1154perform this by supplying an appropriate expression with the "General" restraints, the
     1155GSAS-II authors do not feel that use of restraints or constraints are a good solution for
     1156this common problem where parameters refine to non-physical values. This is because when
     1157this occurs, most likely one of the following cases is occurring:
     1158
     1159#. there is a significant problem
     1160   with the model, for example for an x-ray fit if an O atom is placed where a S is actually
     1161   present, the Uiso will refine artificially small or the occupancy much larger than unity
     1162   to try to compensate for the missing electrons; or
     1163 
     1164#. the data are simply insensitive
     1165   to the parameter or combination of parameters, for example unless very high-Q data
     1166   are included, the effects of a occupancy and Uiso value can have compensating effects,
     1167   so an assumption must be made; likewise, with neutron data natural-abundance V atoms
     1168   are nearly invisible due to weak coherent scattering. No parameters can be fit for a
     1169   V atom with neutrons.
     1170
     1171#. the parameter is non-physical (such as a negative Uiso value) but within
     1172   two sigma (sigma = standard uncertainty, aka e.s.d.) of a reasonable value,
     1173   in which case the
     1174   value is not problematic as it is experimentally indistinguishable from an
     1175   expected value. 
     1176
     1177#. there is a systematic problem with the data (experimental error)
     1178
     1179In all these cases, this situation needs to be reviewed by a crystallographer to decide
     1180how to best determine a structural model for these data. An implementation with a constraint
     1181or restraint is likely to simply hide the problem from the user, making it more probable
     1182that a poor model choice is obtained.
     1183
     1184What GSAS-II does implement is to allow users to specify ranges for parameters
     1185that works by disabling
     1186refinement of parameters that refine beyond either a lower limit or an upper limit, where
     1187either or both may be optionally specified. Parameters limits are specified in the Controls
     1188tree entry in dicts named as ``Controls['parmMaxDict']`` and ``Controls['parmMinDict']``, where
     1189the keys are named using the standard GSAS-II variable
     1190(see :func:`getVarDescr` and :func:`CompileVarDesc`) or with a variable name, where a
     1191wildcard ('*') is used for histogram number or atom number (phase number is intentionally not
     1192allowed as a wildcard as it makes little sense to group together different phases). Note
     1193that func:`prmLookup` is used to see if a name matches a wildcard. The upper or lower limit
     1194is placed into these dicts as a float value. These values can be edited using the window
     1195created by the Calculate/"View LS parms" menu command or in scripting with the
     1196:meth:`GSASIIscriptable.G2Project.set_Controls` function.
     1197In the GUI, a checkbox labeled "match all histograms/atoms" is used to insert a wildcard
     1198into the appropriate part of the variable name.
     1199
     1200When a refinement is conducted, routine :func:`GSASIIstrMain.dropOOBvars` is used to
     1201find parameters that have refined to values outside their limits. If this occurs, the parameter
     1202is set to the limiting value and the variable name is added to a list of frozen variables kept in the
     1203``Controls['parmFrozen']`` dict. In a sequential refinement, this is kept as a list in
     1204``Controls['parmFrozen'][histogram]`` (where the key is the histogram name) or in
     1205``Controls['parmFrozen']['FrozenList']`` for a non-sequential fit.
     1206This allows different variables
     1207to be frozen in each section of a sequential fit.
     1208Frozen parameters are not included in refinements through removal from the
     1209list of parameters to be refined (``varyList``) in :func:`GSASIIstrMain.Refine` or
     1210:func:`GSASIIstrMain.SeqRefine`. Once a variable is frozen, it will not be refined in any
     1211future refinements unless the the variable is removed (manually) from the list. This can also
     1212be done with the Calculate/"View LS parms" menu window or ...
     1213
    11431214*Classes and routines*
    11441215----------------------
  • trunk/GSASIIplot.py

    r4573 r4578  
    209209timeDebug = GSASIIpath.GetConfigValue('Show_timing',False)
    210210obsInCaption = True # include the observed, calc,... items in the plot caption (PlotPatterns)
    211 mplv = eval(mpl.__version__.replace('.',','))
    212211
    213212#matplotlib 2.0.x dumbed down Paired to 16 colors -
     
    333332       
    334333    def ToolBarDraw(self):
     334        mplv = eval(mpl.__version__.replace('.',','))
    335335        if mplv[0] >= 3 and mplv[1] >= 3:
    336336            self.toolbar.canvas.draw_idle()
     
    1057510575
    1057610576def SetupLegendPick(legend,new,delay=5):
     10577    mplv = eval(mpl.__version__.replace('.',','))
    1057710578    legend.delay = delay*1000 # Hold time in ms for clear; 0 == forever
    1057810579    for line in legend.get_lines():
  • trunk/GSASIIpwdGUI.py

    r4567 r4578  
    15201520
    15211521        def OnPeaks(event):
     1522            'Respond to a change in the number of background peaks'
    15221523            data[1]['nPeaks'] = int(peaks.GetValue())
    15231524            M = len(data[1]['peaksList'])
     
    15311532            if N == 0:
    15321533                CalcBack(G2frame.PatternId)
    1533                 G2plt.PlotPatterns(G2frame,plotType='PWDR')               
     1534                G2plt.PlotPatterns(G2frame,plotType='PWDR')
     1535            # this callback is crashing wx when there is an open
     1536            # peaksGrid cell editor, at least on Mac. Code below
     1537            # should fix this, but it does not.
     1538            # https://stackoverflow.com/questions/64082199/wxpython-grid-destroy-with-open-celleditor-crashes-python-even-with-disablece
     1539            if peaksGrid and peaksGrid.IsCellEditControlEnabled():
     1540                # complete any grid edits in progress
     1541                #print('closing')
     1542                peaksGrid.HideCellEditControl()
     1543                peaksGrid.DisableCellEditControl()
     1544                #wx.CallLater(100,peaksGrid.Destroy) # crashes python
    15341545            wx.CallAfter(UpdateBackground,G2frame,data)
    15351546           
     
    15671578        peaksSizer.Add(topSizer)
    15681579        G2frame.dataWindow.currentGrids = []
     1580        peaksGrid = None
    15691581        if data[1]['nPeaks']:
    15701582            peaksSizer.Add(wx.StaticText(G2frame.dataWindow,-1,' Peak list:'),0,WACV)       
  • trunk/GSASIIscriptable.py

    r4551 r4578  
    33603360            raise Exception('{} is an invalid control value'.format(control))
    33613361       
    3362     def set_Controls(self, control, value):
     3362    def set_Controls(self, control, value, variable=None):
    33633363        '''Set project controls
    33643364
    33653365        :param str control: the item to be set. See below for allowed values.
    33663366        :param value: the value to be set.
     3367        :param str variable: used only with control set to "parmMin" or "parmMax"
    33673368
    33683369        Allowed values for parameter control:
     
    33783379            * Reverse Seq: when True, sequential refinement is performed on the
    33793380              reversed list of histograms.
     3381            * parmMin & parmMax: set a maximum or minimum value for a refined
     3382              parameter. Note that variable will be a GSAS-II variable name,
     3383              optionally with * specified for a histogram or atom number and
     3384              value must be a float.
     3385              (See :ref:`Parameter Limits<ParameterLimits>` description.)
    33803386
    33813387        .. seealso::
    33823388            :meth:`get_Controls`
    33833389        '''
     3390        if control.startswith('parmM'):
     3391            if not variable:
     3392                raise Exception('set_Controls requires a value for variable for control=parmMin or parmMax')
     3393            for key in ('parmMinDict','parmMaxDict','parmFrozen'):
     3394                if key not in self.data['Controls']['data']: self.data['Controls']['data'][key] = {}
    33843395        if control == 'cycles':
    33853396            self.data['Controls']['data']['max cyc'] = int(value)
     
    33983409                                        .format(i,j))
    33993410            self.data['Controls']['data']['Seq Data'] = histlist
     3411        elif control == 'parmMin' or control == 'parmMax':
     3412            self.data['Controls']['data'][control+'Dict'][variable] = float(value)
    34003413        else:
    34013414            raise Exception('{} is an invalid control value'.format(control))
  • trunk/GSASIIstrMain.py

    r4534 r4578  
    6161        ifPrint = False
    6262    Rvals = {}
     63    chisq0 = None
    6364    while True:
    6465        begin = time.time()
     
    9192            Rvals['lamMax'] = result[2]['lamMax']
    9293            Controls['Marquardt'] = -3  #reset to default
     94            if 'chisq0' in result[2] and chisq0 is None:
     95                chisq0 = result[2]['chisq0']
    9396        elif 'Hessian SVD' in Controls['deriv type']:
    9497            maxCyc = Controls['max cyc']
     
    102105                break
    103106            ncyc = result[2]['num cyc']+1
     107            if 'chisq0' in result[2] and chisq0 is None:
     108                chisq0 = result[2]['chisq0']
    104109        else:           #'numeric'
    105110            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
     
    182187    if IfOK:
    183188        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
     189    if chisq0 is not None:
     190        Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
    184191    return IfOK,Rvals,result,covMatrix,sig
    185192
     
    293300            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
    294301            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
    295             if len(frozen) > 0:
    296                 G2fil.G2Print(
    297                     ' {} variables were outside limits and were frozen (now {} frozen total)\n'
    298                     .format(len(frozen),len(parmFrozenList)))
     302            if len(frozen):
     303                msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
     304                    .format(len(frozen),len(parmFrozenList))
     305                    )
     306                G2fil.G2Print(msg)
     307                Rvals['msg'] = msg
     308            elif len(parmFrozenList):
     309                msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
     310                    .format(len(parmFrozenList))
     311                    )
     312                G2fil.G2Print('Note: ',msg)
     313                Rvals['msg'] = msg
    299314            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
    300315            printFile.close()
  • trunk/atmdata.py

    r4337 r4578  
    11901190        'N':[    1.8845,  9.161, 2.0746,  4.042,  -.1318,  1.723,  .0000,   .000,  .0020, 2.00]},
    11911191}
    1192 '''
    1193 unused j<0>
     1192
     1193pass # placed here so that the following is not processed by Sphinx as a doc string
     1194''' # unused unused terms follow
     1195j<0>
    11941196W_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
    11951197W_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

    r3157 r4578  
    3434defaultIparm_lbl = []
    3535defaultIparms = []
     36
     37pass # placed here so that the following is not processed by Sphinx as a doc string
    3638''' the phrase 'lab data' in the following makes the default goniometer type='Bragg-Brentano'
    3739    otherwise it is set to 'Debye-Scherrer'
Note: See TracChangeset for help on using the changeset viewer.