Changeset 4076 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Aug 2, 2019 1:33:41 PM (2 years ago)
Author:
vondreele
Message:

fix axes amnesia for PDF plots
add plotting of P(R) pair distribution for small angle data
force element selection if none for 1st pdf calc
Size & pair dist plots now hollow bars
implement Moore method for making pair distribution P(R) from protein small angle data

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r4072 r4076  
    10781078    print (' Final chi^2: %.3f'%(chisq))
    10791079    data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     1080   
     1081################################################################################
     1082#### Modelling
     1083################################################################################
     1084
     1085def PairDistFxn(Profile,ProfDict,Limits,Sample,data):
     1086    print('create P(R) fit to data - TBD')
     1087   
     1088    def CalcMoore():
     1089       
     1090#        def MooreRg(cw,q):      #lines 1400-1409
     1091#            n = 0
     1092#            nmax = len(cw)-5
     1093#            POR = 0.
     1094#            dmax = cw[nmax+2]
     1095#            while n < nmax:
     1096#                POR += cw[n]*((-1**n)*((np.pi*(n+1))**2-6)/((n+1)**3))
     1097#                n += 1
     1098#            POR *= 4*(dmax**4)/((np.pi**2)*cw[nmax+3])
     1099#            return POR
     1100       
     1101        def MoorePOR(cw,r,dmax):    #lines 1417-1428
     1102            n = 0
     1103            nmax = len(cw)
     1104            POR = np.zeros(len(r))
     1105            while n < nmax:
     1106                POR += 4.*r*cw[n]*np.sin((n+1.)*np.pi*r/dmax)
     1107                n += 1
     1108            return POR
     1109       
     1110        def MooreIOREFF(cw,q,dmax):      #lines 1437-1448
     1111            n = 0
     1112            nmax = len(cw)
     1113            POR = np.zeros(len(q))
     1114            dq = dmax*q
     1115            fpd = 8.*(np.pi**2)*dmax*np.sin(dq)/q
     1116            while n < nmax:
     1117                POR += cw[n]*(n+1.)*((-1)**n)*fpd/(((n+1.)*np.pi)**2-dq**2)
     1118                n += 1
     1119            return POR
     1120       
     1121#        def MooreINaught(cw,q,dmax):         #lines 1457-1466
     1122#            n = 0
     1123#            nmax = len(cw)
     1124#            POR = 0.
     1125#            while n < nmax:
     1126#                POR += cw[n]*8*(dmax**2)*((-1.)**n)/(n+1)
     1127#                n += 1
     1128#            return POR
     1129       
     1130        def calcSASD(values,Q,Io,wt,Ifb,dmax,ifBack):
     1131            if ifBack:
     1132                M = np.sqrt(wt)*(MooreIOREFF(values[:-1],Q,dmax)+Ifb+values[-1]-Io)
     1133            else:
     1134                M = np.sqrt(wt)*(MooreIOREFF(values,Q,dmax)+Ifb-Io)
     1135            return M
     1136
     1137        Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     1138        Qmin = Limits[1][0]
     1139        Qmax = Limits[1][1]
     1140        wtFactor = ProfDict['wtFactor']
     1141        Back,ifBack = data['Back']
     1142        Ibeg = np.searchsorted(Q,Qmin)
     1143        Ifin = np.searchsorted(Q,Qmax)+1    #include last point
     1144        Ic[Ibeg:Ifin] = 0
     1145        Bins = np.linspace(0.,pairData['MaxRadius'],pairData['NBins']+1,True)
     1146        Dbins = np.diff(Bins)
     1147        Bins = Bins[:-1]+Dbins/2.
     1148        N = pairData['Moore']
     1149        if ifBack:
     1150            N += 1
     1151        MPV = np.zeros(N)
     1152        MPS = np.zeros(N)
     1153        MPV[0] = Q[Ibeg]
     1154        dmax = pairData['MaxRadius']
     1155        result = so.leastsq(calcSASD,MPV,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     1156            args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],dmax,ifBack))
     1157        if ifBack:
     1158            MPVR = result[0][:-1]
     1159            data['Back'][0] = result[0][-1]
     1160            Back = data['Back'][0]
     1161        else:       
     1162            MPVR = result[0]
     1163            Back = 0.
     1164        chisq = np.sum(result[2]['fvec']**2)
     1165        Ic[Ibeg:Ifin] = MooreIOREFF(MPVR,Q[Ibeg:Ifin],dmax)+Ifb+Back
     1166        covM = result[1]
     1167        GOF = chisq/(Ifin-Ibeg-N)
     1168        MPS = np.sqrt(np.diag(covM)*GOF)
     1169        BinMag = MoorePOR(MPVR,Bins,dmax)/2.
     1170        return Bins,Dbins,BinMag
     1171   
     1172    def CalcRegular():
     1173        Bins = np.linspace(0.,pairData['MaxRadius'],pairData['NBins']+1,True)
     1174        Dbins = np.diff(Bins)
     1175        Bins = Bins[:-1]+Dbins/2.
     1176
     1177#test       
     1178        midBin = pairData['MaxRadius']/4.
     1179        wid = midBin/4.
     1180        BinMag  = 1./np.sqrt(2.*np.pi*wid*2)*np.exp(-(midBin-Bins)**2/(2.*wid**2))
     1181        return Bins,Dbins,BinMag
     1182               
     1183    pairData = data['Pair']
     1184   
     1185    if pairData['Method'] == 'Regularization':
     1186        print('Regularization calc; dummy Gaussian - TBD')
     1187        Bins,Dbins,BinMag = CalcRegular()
     1188       
     1189       
     1190    elif pairData['Method'] == 'Moore':
     1191        Bins,Dbins,BinMag = CalcMoore()       
     1192   
     1193    data['Pair']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)]
     1194   
     1195   
    10801196       
    10811197################################################################################
     
    12921408    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
    12931409        return False,0,0,0,0,0,0,Msg
     1410
     1411def RgFit(Profile,ProfDict,Limits,Sample,Model):
     1412    print('unified fit single Rg to data to estimate a size')
     1413   
     1414    def GetModelParms():
     1415        parmDict = {}
     1416        varyList = []
     1417        values = []
     1418        Back = Model['Back']
     1419        if Back[1]:
     1420            varyList += ['Back',]
     1421            values.append(Back[0])
     1422        parmDict['Back'] = Back[0]
     1423        pairData = Model['Pair']
     1424        parmDict['G'] = pairData.get('Dist G',100.)
     1425        parmDict['Rg'] = pairData['MaxRadius']/2.5
     1426        varyList  += ['G','Rg',]
     1427        values += [parmDict['G'],parmDict['Rg'],]
     1428        return parmDict,varyList,values
     1429
     1430    def calcSASD(values,Q,Io,wt,Ifb,parmDict,varyList):
     1431        parmDict.update(zip(varyList,values))
     1432        M = np.sqrt(wt)*(getSASD(Q,parmDict)+Ifb-Io)
     1433        return M
     1434   
     1435    def getSASD(Q,parmDict):
     1436        Ic = np.zeros_like(Q)
     1437        Rg,G = parmDict['Rg'],parmDict['G']
     1438        Guin = G*np.exp(-(Q*Rg)**2/3)
     1439        Ic += Guin
     1440        Ic += parmDict['Back']  #/parmDict['Scale']
     1441        return Ic
     1442       
     1443    def SetModelParms():
     1444        print (' Refined parameters: ')
     1445        if 'Back' in varyList:
     1446            Model['Back'][0] = parmDict['Back']
     1447            print ('  %15s %15.4f esd: %15.4g'%('Background:',parmDict['Back'],sigDict['Back']))
     1448        pairData = Model['Pair']
     1449        pairData['Dist G'] = parmDict['G']
     1450        pairData['MaxRadius'] = parmDict['Rg']*2.5
     1451        for item in varyList:
     1452            print (' %15s: %15.4g esd: %15.4g'%(item,parmDict[item],sigDict[item]))
     1453
     1454    Q,Io,wt,Ic,Ib,Ifb = Profile[:6]
     1455    Qmin = Limits[1][0]
     1456    Qmax = Limits[1][1]
     1457    wtFactor = ProfDict['wtFactor']
     1458    Ibeg = np.searchsorted(Q,Qmin)
     1459    Ifin = np.searchsorted(Q,Qmax)+1    #include last point
     1460    Ic[:] = 0
     1461    parmDict,varyList,values = GetModelParms()
     1462    result = so.leastsq(calcSASD,values,full_output=True,epsfcn=1.e-8,   #ftol=Ftol,
     1463        args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Ifb[Ibeg:Ifin],parmDict,varyList))
     1464    parmDict.update(zip(varyList,result[0]))
     1465    chisq = np.sum(result[2]['fvec']**2)
     1466    ncalc = result[2]['nfev']
     1467    covM = result[1]
     1468    Rvals = {}
     1469    Rvals['Rwp'] = np.sqrt(chisq/np.sum(wt[Ibeg:Ifin]*Io[Ibeg:Ifin]**2))*100.      #to %
     1470    Rvals['GOF'] = chisq/(Ifin-Ibeg-len(varyList))       #reduced chi^2
     1471    Ic[Ibeg:Ifin] = getSASD(Q[Ibeg:Ifin],parmDict)
     1472    Msg = 'Failed to converge'
     1473    try:
     1474        Nans = np.isnan(result[0])
     1475        if np.any(Nans):
     1476            name = varyList[Nans.nonzero(True)[0]]
     1477            Msg = 'Nan result for '+name+'!'
     1478            raise ValueError
     1479        Negs = np.less_equal(result[0],0.)
     1480        if np.any(Negs):
     1481            name = varyList[Negs.nonzero(True)[0]]
     1482            Msg = 'negative coefficient for '+name+'!'
     1483            raise ValueError
     1484        if len(covM):
     1485            sig = np.sqrt(np.diag(covM)*Rvals['GOF'])
     1486            sigDict = dict(zip(varyList,sig))
     1487        print (' Results of Rg fit:')
     1488        print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(ncalc,Ifin-Ibeg,len(varyList)))
     1489        print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF']))
     1490        SetModelParms()
     1491        covMatrix = covM*Rvals['GOF']
     1492        return True,result,varyList,sig,Rvals,covMatrix,parmDict,''
     1493    except (ValueError,TypeError):      #when bad LS refinement; covM missing or with nans
     1494        return False,0,0,0,0,0,0,Msg
     1495
     1496    return [None,]
     1497
    12941498   
    12951499def ModelFxn(Profile,ProfDict,Limits,Sample,sasdData):
Note: See TracChangeset for help on using the changeset viewer.