Changeset 2838


Ignore:
Timestamp:
May 24, 2017 2:32:22 PM (5 years ago)
Author:
vondreele
Message:

add a bit of printing to LS & show no. of penalty fxns

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r2775 r2838  
    120120        Yvec /= Adiag
    121121        Amat /= Anorm
     122        if Print:
     123            print 'initial chi^2 %.5g'%(chisq0)
    122124        while True:
    123125            Lam = np.eye(Amat.shape[0])*lam
     
    136138                lam *= 10.
    137139                if Print:
    138                     print 'matrix modification needed; lambda now %.1e'%(lam)
     140                    print 'new chi^2 %.5g; matrix modification needed; lambda now %.1e'%(chisq1,lam)
    139141            else:
    140142                x0 += Xvec
     
    160162    Anorm = np.outer(Adiag,Adiag)
    161163    Lam = np.eye(Amat.shape[0])*lam
    162     Amatlam = Amat/Anorm  #*(One+Lam)              #don't scale Amat to Marquardt array       
     164    Amatlam = Amat/Anorm       
     165#    Amatlam = Amat*(One+Lam)                #scale Amat to Marquardt array?       
    163166    try:
    164         Bmat = nl.inv(Amatlam)/Anorm  #*(One+Lam)      #don't rescale Bmat to Marquardt array
     167        Bmat = nl.inv(Amatlam)/Anorm
     168#        Bmat = Bmat*(One+Lam)               #rescale Bmat to Marquardt array?
    165169        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
    166170    except nl.LinAlgError:
  • trunk/GSASIIpwd.py

    r2808 r2838  
    954954def getFCJVoigt3(pos,sig,gam,shl,xdata):
    955955    'needs a doc string'
    956    
    957956    Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl)
    958957#    Df = pyd.pypsvfcjo(len(xdata),xdata-pos,pos,sig,gam,shl)
  • trunk/GSASIIstrMath.py

    r2837 r2838  
    306306    negWt = {}
    307307    pWsum = {}
     308    pWnum = {}
    308309    for phase in Phases:
    309310        pId = Phases[phase]['pId']
     
    325326        for name,rest in names:
    326327            pWsum[name] = 0.
     328            pWnum[name] = 0
    327329            itemRest = phaseRest[name]
    328330            if itemRest[rest] and itemRest['Use']:
     
    344346                        pWt.append(wt/esd**2)
    345347                        pWsum[name] += wt*((obs-calc)/esd)**2
     348                        pWnum[name] += 1
    346349                elif name in ['Torsion','Rama']:
    347350                    coeffDict = itemRest['Coeff']
     
    359362                        pWt.append(wt/esd**2)
    360363                        pWsum[name] += wt*(restr/esd)**2
     364                        pWnum[name] += 1
    361365                elif name == 'ChemComp':
    362366                    for i,[indx,factors,obs,esd] in enumerate(itemRest[rest]):
     
    368372                        pWt.append(wt/esd**2)                   
    369373                        pWsum[name] += wt*((obs-calc)/esd)**2
     374                        pWnum[name] += 1
    370375                elif name == 'Texture':
    371376                    SHkeys = textureData['SH Coeff'][1].keys()
     
    385390                            pWt.append(wt/esd1**2)
    386391                            pWsum[name] += wt*(-Z1[ind[0]][ind[1]]/esd1)**2
     392                            pWnum[name] += 1
    387393                        if ifesd2:
    388394                            Z2 = 1.-Z
     
    392398                                pWt.append(wt/esd2**2)
    393399                                pWsum[name] += wt*(Z2/esd2)**2
     400                                pWnum[name] += 1
    394401       
    395402    for phase in Phases:
     
    400407        cell = General['Cell'][1:7]
    401408        pWsum[name] = 0.0
     409        pWnum[name] = 0
    402410        for hist in Phases[phase]['Histograms']:
    403411            if not Phases[phase]['Histograms'][hist]['Use']:
     
    427435                            pWt.append(wt)
    428436                            pWsum[name] += wt*(Y[ind])**2
     437                            pWnum[name] += 1
    429438    pWsum['PWLref'] = 0.
     439    pWnum[name] = 0
    430440    for item in varyList:
    431441        if 'PWLref' in item and parmDict[item] < 0.:
     
    436446                pWt.append(negWt[pId])
    437447                pWsum['PWLref'] += negWt[pId]*(parmDict[item])**2
     448                pWnum[name] += 1
    438449    pVals = np.array(pVals)
    439450    pWt = np.array(pWt)         #should this be np.sqrt?
    440     return pNames,pVals,pWt,pWsum
     451    return pNames,pVals,pWt,pWsum,pWnum
    441452   
    442453def penaltyDeriv(pNames,pVal,HistoPhases,calcControls,parmDict,varyList):
     
    710721#reflection processing begins here - big arrays!
    711722    iBeg = 0
     723    time0 = time.time()
    712724    while iBeg < nRef:
    713725        iFin = min(iBeg+blkSize,nRef)
     
    795807                    refl.T[7] = np.copy(refl.T[9])               
    796808                    refl.T[10] = atan2d(fbs[0],fas[0])  #ignore f' & f"
    797 #        GSASIIpath.IPyBreak()
    798809#                refl.T[10] = atan2d(np.sum(fbs,axis=0),np.sum(fas,axis=0)) #include f' & f"
    799810        iBeg += blkSize
    800 #    print ' %d sf time %.4f\r'%(nRef,time.time()-time0)
     811#    print 'sf time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    801812   
    802813def StructureFactorDerv2(refDict,G,hfx,pfx,SGData,calcControls,parmDict):
     
    849860    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    850861        Flack = 1.-2.*parmDict[phfx+'Flack']
    851 #    time0 = time.time()
     862    time0 = time.time()
    852863#reflection processing begins here - big arrays!
    853864    iBeg = 0
     
    940951#        GSASIIpath.IPyBreak()
    941952        iBeg += blkSize
    942 #    print ' %d derivative time %.4f\r'%(nRef,time.time()-time0)
     953#    print 'derv time %.4f, nref %d, blkSize %d'%(time.time()-time0,nRef,blkSize)
    943954        #loop over atoms - each dict entry is list of derivatives for all the reflections
    944955    for i in range(len(Mdata)):
     
    30993110            else:
    31003111                StructureFactor2(refDict,G,hfx,pfx,SGData,calcControls,parmDict)
    3101 #            print 'sf calc time: %.3fs'%(time.time()-time0)
    31023112        badPeak = False
    31033113        for iref,refl in enumerate(refDict['RefList']):
     
    35723582                    for item in ['BabA','BabU']:
    35733583                        if phfx+item in varylist:
    3574                             dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
     3584                            dMdvh[varylist.index(phfx+item)][iref] = w*dFdvDict[phfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
    35753585                        elif phfx+item in dependentVars:
    3576                             depDerivDict[phfx+item][iref] = w*dFdvDict[pfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
     3586                            depDerivDict[phfx+item][iref] = w*dFdvDict[phfx+item][iref]*parmDict[phfx+'Scale']*ref[11+im]
    35773587    else:   #F refinement
    35783588        for iref,ref in enumerate(refDict['RefList']):
     
    36533663           
    36543664    GetFobsSq(Histograms,Phases,parmDict,calcControls)
    3655     pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
     3665    pNames,pVals,pWt,pWsum,pWnum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
    36563666    if np.any(pVals):
    36573667        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
     
    37303740            continue        #skip non-histogram entries
    37313741    GetFobsSq(Histograms,Phases,parmDict,calcControls)
    3732     pNames,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
     3742    pNames,pVals,pWt,pWsum,pWnum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
    37333743    if np.any(pVals):
    37343744        dpdv = penaltyDeriv(pNames,pVals,HistoPhases,calcControls,parmDict,varylist)
     
    39583968            dlg.Destroy()
    39593969            raise G2obj.G2Exception('User abort')         #Abort!!
    3960     pDict,pVals,pWt,pWsum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
     3970    pDict,pVals,pWt,pWsum,pWnum = penaltyFxn(HistoPhases,calcControls,parmDict,varylist)
    39613971    if len(pVals):
    39623972        pSum = np.sum(pWt*pVals**2)
    39633973        for name in pWsum:
    39643974            if pWsum[name]:
    3965                 print '  Penalty function for %8s = %12.5g'%(name,pWsum[name])
     3975                print '  Penalty function for %5d %8ss = %12.5g'%(pWnum[name],name,pWsum[name])
    39663976        print 'Total penalty function: %12.5g on %d terms'%(pSum,len(pVals))
    39673977        Nobs += len(pVals)
Note: See TracChangeset for help on using the changeset viewer.