Changeset 4094 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Aug 15, 2019 2:54:58 PM (2 years ago)
Author:
vondreele
Message:

SHAPES now complete, but slow.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r4077 r4094  
    666666    DISTANCE_LIMIT_FACTOR = 0.1                 # limitation on df to constrain runaways
    667667   
    668     MAX_MOVE_LOOPS    = 500                     # for no solution in routine: move,
     668    MAX_MOVE_LOOPS    = 5000                     # for no solution in routine: move,
    669669    MOVE_PASSES       = 0.001                   # convergence test in routine: move
    670670
     
    11281128        if ifBack:
    11291129            N += 1
    1130         MPV = np.zeros(N)
    1131         MPV[0] = Q[Ibeg]
     1130        MPV = np.ones(N)*1.e-5
     1131#        MPV[0] = Q[Ibeg]
    11321132        dmax = pairData['MaxRadius']
    11331133        result = so.leastsq(calcSASD,MPV,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     
    11411141            Back = 0.
    11421142        chisq = np.sum(result[2]['fvec']**2)
    1143         Ic[Ibeg:Ifin] = MooreIOREFF(MPVR,Q[Ibeg:Ifin],dmax)+Ifb+Back
     1143        covM = result[1]
     1144        Ic[Ibeg:Ifin] = MooreIOREFF(MPVR,Q[Ibeg:Ifin],dmax)+Ifb[Ibeg:Ifin]+Back
    11441145        ncalc = result[2]['nfev']
    11451146        GOF = chisq/(Ifin-Ibeg-N)
     
    11481149        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,N))
    11491150        print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF))
     1151        if len(covM):
     1152            sig = np.sqrt(np.diag(covM)*GOF)
     1153            for val,esd in zip(result[0],sig):
     1154                print(' parameter: %.4g esd: %.4g'%(val,esd))           
    11501155        BinMag = MoorePOR(MPVR,Bins,dmax)/2.
    11511156        return Bins,Dbins,BinMag
    11521157   
    1153     def CalcRegular():
    1154         Bins = np.linspace(0.,pairData['MaxRadius'],pairData['NBins']+1,True)
    1155         Dbins = np.diff(Bins)
    1156         Bins = Bins[:-1]+Dbins/2.
    1157 
    1158 #test       
    1159         midBin = pairData['MaxRadius']/4.
    1160         wid = midBin/4.
    1161         BinMag  = 1./np.sqrt(2.*np.pi*wid*2)*np.exp(-(midBin-Bins)**2/(2.*wid**2))
    1162         return Bins,Dbins,BinMag
    1163                
    11641158    pairData = data['Pair']
    11651159   
    1166     if pairData['Method'] == 'Regularization':
     1160    if pairData['Method'] == 'Regularization':      #not used
    11671161        print('Regularization calc; dummy Gaussian - TBD')
    1168         Bins,Dbins,BinMag = CalcRegular()
     1162        pairData['Method'] = 'Moore'
    11691163       
    11701164       
    11711165    elif pairData['Method'] == 'Moore':
    1172         Bins,Dbins,BinMag = CalcMoore()       
    1173    
    1174     data['Pair']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     1166        Bins,Dbins,BinMag = CalcMoore()
     1167        BinSum = np.sum(BinMag*Dbins)
     1168        BinMag /= BinSum       
     1169   
     1170    data['Pair']['Distribution'] = [Bins,Dbins,BinMag]      #/(2.*Dbins)]
     1171    if 'Pair Calc' in data['Pair']:
     1172        del data['Pair']['Pair Calc']
    11751173   
    11761174   
     
    13901388        return False,0,0,0,0,0,0,Msg
    13911389
     1390def getSASDRg(Q,parmDict):
     1391    Ic = np.zeros_like(Q)
     1392    Rg,G,B = parmDict['Rg'],parmDict['G'],parmDict['B']
     1393    Qstar = Q/(scsp.erf(Q*Rg/np.sqrt(6)))**3
     1394    Guin = G*np.exp(-(Q*Rg)**2/3)
     1395    Porod = (B/Qstar**4)        #*np.exp(-(Q*B)**2/3)
     1396    Ic += Guin+Porod+parmDict['Back']
     1397    return Ic
     1398       
    13921399def RgFit(Profile,ProfDict,Limits,Sample,Model):
    13931400    print('unified fit single Rg to data to estimate a size')
     
    14031410        parmDict['Back'] = Back[0]
    14041411        pairData = Model['Pair']
    1405         parmDict['G'] = pairData.get('Dist G',100.)
     1412        parmDict['G'] = pairData.get('Dist G',Io[Ibeg])
    14061413        parmDict['Rg'] = pairData['MaxRadius']/2.5
    1407         varyList  += ['G','Rg',]
    1408         values += [parmDict['G'],parmDict['Rg'],]
     1414        parmDict['B'] = pairData.get('Dist B',Io[Ifin-1]*Q[Ifin-1]**4)
     1415        varyList  += ['G','Rg','B']
     1416        values += [parmDict['G'],parmDict['Rg'],parmDict['B']]
    14091417        return parmDict,varyList,values
    14101418
    14111419    def calcSASD(values,Q,Io,wt,Ifb,parmDict,varyList):
    14121420        parmDict.update(zip(varyList,values))
    1413         M = np.sqrt(wt)*(getSASD(Q,parmDict)+Ifb-Io)
     1421        M = np.sqrt(wt)*(getSASDRg(Q,parmDict)-Io)
    14141422        return M
    14151423   
    1416     def getSASD(Q,parmDict):
    1417         Ic = np.zeros_like(Q)
    1418         Rg,G = parmDict['Rg'],parmDict['G']
    1419         Guin = G*np.exp(-(Q*Rg)**2/3)
    1420         Ic += Guin
    1421         Ic += parmDict['Back']  #/parmDict['Scale']
    1422         return Ic
    1423        
    14241424    def SetModelParms():
    14251425        print (' Refined parameters: ')
     
    14301430        pairData['Dist G'] = parmDict['G']
    14311431        pairData['MaxRadius'] = parmDict['Rg']*2.5
     1432        pairData['Dist B'] = parmDict['B']
    14321433        for item in varyList:
    14331434            print (' %15s: %15.4g esd: %15.4g'%(item,parmDict[item],sigDict[item]))
     
    14411442    Ic[:] = 0
    14421443    parmDict,varyList,values = GetModelParms()
    1443     result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     1444    result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-12,factor=0.1,  #ftol=Ftol,
    14441445        args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],parmDict,varyList))
    1445     parmDict.update(zip(varyList,result[0]))
     1446    parmDict.update(dict(zip(varyList,result[0])))
    14461447    chisq = np.sum(result[2]['fvec']**2)
    14471448    ncalc = result[2]['nfev']
     
    14501451    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
    14511452    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
    1452     Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],parmDict)
     1453    Ic[Ibeg:Ifin] = getSASDRg(Q[Ibeg:Ifin],parmDict)
    14531454    Msg = 'Failed to converge'
    14541455    try:
Note: See TracChangeset for help on using the changeset viewer.