Changeset 4076 for trunk/GSASIIsasd.py
 Timestamp:
 Aug 2, 2019 1:33:41 PM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIsasd.py
r4072 r4076 1078 1078 print (' Final chi^2: %.3f'%(chisq)) 1079 1079 data['Size']['Distribution'] = [Bins,Dbins,BinMag/(2.*Dbins)] 1080 1081 ################################################################################ 1082 #### Modelling 1083 ################################################################################ 1084 1085 def 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 14001409 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))**26)/((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 14171428 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 14371448 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)**2dq**2) 1118 n += 1 1119 return POR 1120 1121 # def MooreINaught(cw,q,dmax): #lines 14571466 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)+IfbIo) 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.e8, #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/(IfinIbegN) 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((midBinBins)**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 1080 1196 1081 1197 ################################################################################ … … 1292 1408 except (ValueError,TypeError): #when bad LS refinement; covM missing or with nans 1293 1409 return False,0,0,0,0,0,0,Msg 1410 1411 def 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)+IfbIo) 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.e8, #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/(IfinIbeglen(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,IfinIbeg,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 1294 1498 1295 1499 def ModelFxn(Profile,ProfDict,Limits,Sample,sasdData):
Note: See TracChangeset
for help on using the changeset viewer.