Changeset 2123


Ignore:
Timestamp:
Jan 13, 2016 1:19:43 PM (6 years ago)
Author:
vondreele
Message:

code cleanup & simplification in StructureFactor? & Dervs
block reflections for twin SF derv function

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIstrMath.py

    r2122 r2123  
    10181018    dFdbab = np.zeros((nRef,2))
    10191019    dFdfl = np.zeros((nRef))
    1020     dFdtw = np.zeros((nRef))
    10211020    Flack = 1.0
    10221021    if not SGData['SGInv'] and 'S' in calcControls[hfx+'histType'] and phfx+'Flack' in parmDict:
    … …  
    11001099        SB = fbs[0]+fbs[1]
    11011100        if 'P' in calcControls[hfx+'histType']: #checked perfect for centro & noncentro
    1102             dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
    1103                 2.*(fbs[0,:,nxs]*dfbdfr[0]-fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
    1104             dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
    1105                 2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
    1106             dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
    1107                 2.*(fbs[0,:,nxs]*dfbdui[0]-fbs[1,:,nxs]*dfbdui[1])
    1108             dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
    1109                 2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
     1101            dFdfr[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadfr+fbs[:,:,nxs]*dfbdfr,axis=0)*Mdata/len(SGMT)
     1102#            dFdfr[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadfr[0]+fas[1,:,nxs]*dfadfr[1])*Mdata/len(SGMT)+   \
     1103#                2.*(fbs[0,:,nxs]*dfbdfr[0]+fbs[1,:,nxs]*dfbdfr[1])*Mdata/len(SGMT)
     1104            dFdx[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadx+fbs[:,:,nxs,nxs]*dfbdx,axis=0)
     1105#            dFdx[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadx[0]+fas[1,:,nxs,nxs]*dfadx[1])+  \
     1106#                2.*(fbs[0,:,nxs,nxs]*dfbdx[0]+fbs[1,:,nxs,nxs]*dfbdx[1])
     1107            dFdui[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs]*dfadui+fbs[:,:,nxs]*dfbdui,axis=0)
     1108#            dFdui[iBeg:iFin] = 2.*(fas[0,:,nxs]*dfadui[0]+fas[1,:,nxs]*dfadui[1])+   \
     1109#                2.*(fbs[0,:,nxs]*dfbdui[0]+fbs[1,:,nxs]*dfbdui[1])
     1110            dFdua[iBeg:iFin] = 2.*np.sum(fas[:,:,nxs,nxs]*dfadua+fbs[:,:,nxs,nxs]*dfbdua,axis=0)
     1111#            dFdua[iBeg:iFin] = 2.*(fas[0,:,nxs,nxs]*dfadua[0]+fas[1,:,nxs,nxs]*dfadua[1])+   \
     1112#                2.*(fbs[0,:,nxs,nxs]*dfbdua[0]+fbs[1,:,nxs,nxs]*dfbdua[1])
    11101113        else:
    11111114            dFdfr[iBeg:iFin] = (2.*SA[:,nxs]*(dfadfr[0]+dfadfr[1])+2.*SB[:,nxs]*(dfbdfr[0]+dfbdfr[1]))*Mdata/len(SGMT)
    … …  
    11911194#            FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
    11921195#        H = np.array(refl[:3])
    1193 #        H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
     1196#        H = np.inner(H.T,TwinLaw)   #maybe array(3,nTwins) or (3)
    11941197#        TwMask = np.any(H,axis=-1)
    11951198#        if iref in TwDict:
    … …  
    12141217#        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    12151218#        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    1216 #        Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
     1219#        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6))
    12171220#        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    12181221#        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
    … …  
    13321335    dFdtw = np.zeros((nRef,nTwin))
    13331336    time0 = time.time()
    1334     nref = len(refDict['RefList'])/100   
    1335     for iref,refl in enumerate(refDict['RefList']):
     1337#reflection processing begins here - big arrays!
     1338    iBeg = 0
     1339    blkSize = 16       #no. of reflections in a block - optimized for speed
     1340    while iBeg < nRef:
     1341        iFin = min(iBeg+blkSize,nRef)
     1342        refl = refDict['RefList'][iBeg:iFin]    #array(blkSize,nItems)
     1343        H = refl.T[:3]
     1344        H = np.inner(H.T,TwinLaw)   #array(3,nTwins)
     1345        TwMask = np.any(H,axis=-1)
     1346        for ir in range(blkSize):
     1347            iref = ir+iBeg
     1348            if iref in TwDict:
     1349                for i in TwDict[iref]:
     1350                    for n in range(NTL):
     1351                        H[ir][i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
     1352        TwMask = np.any(H,axis=-1)
     1353        SQ = 1./(2.*refl.T[4])**2             # or (sin(theta)/lambda)**2
     1354        SQfactor = 8.0*SQ*np.pi**2
    13361355        if 'T' in calcControls[hfx+'histType']:
    1337             FP,FPP = G2el.BlenResCW(Tdata,BLtables,refl.T[12])
    1338         H = np.array(refl[:3])
    1339         H = np.squeeze(np.inner(H.T,TwinLaw))   #maybe array(3,nTwins) or (3)
    1340         TwMask = np.any(H,axis=-1)
    1341         if iref in TwDict:
    1342             for i in TwDict[iref]:
    1343                 for n in range(NTL):
    1344                     H[i+n*NM] = np.inner(TwinLaw[n*NM],np.array(TwDict[iref][i])*TwinInv[i+n*NM])
    1345             TwMask = np.any(H,axis=-1)
    1346         SQ = 1./(2.*refl[4])**2             # or (sin(theta)/lambda)**2
    1347         SQfactor = 8.0*SQ*np.pi**2
     1356            if 'P' in calcControls[hfx+'histType']:
     1357                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[14])
     1358            else:
     1359                FP,FPP = G2el.BlenResTOF(Tdata,BLtables,refl.T[12])
     1360            FP = np.repeat(FP.T,len(SGT)*len(TwinLaw),axis=0)
     1361            FPP = np.repeat(FPP.T,len(SGT)*len(TwinLaw),axis=0)
    13481362        dBabdA = np.exp(-parmDict[phfx+'BabU']*SQfactor)
    1349         Bab = parmDict[phfx+'BabA']*dBabdA
     1363        Bab = np.repeat(parmDict[phfx+'BabA']*dBabdA,len(SGT)*nTwin)
    13501364        Tindx = np.array([refDict['FF']['El'].index(El) for El in Tdata])
    1351         FF = refDict['FF']['FF'][iref].T[Tindx].T
    1352         Uniq = np.inner(H,SGMT)             # array(nSGOp,3) or (nTwin,nSGOp,3)
     1365        FF = np.repeat(refDict['FF']['FF'][iBeg:iFin].T[Tindx].T,len(SGT)*len(TwinLaw),axis=0)
     1366        Uniq = np.inner(H,SGMT)             # (nTwin,nSGOp,3)
    13531367        Phi = np.inner(H,SGT)
    13541368        phase = twopi*(np.inner(Uniq,(dXdata+Xdata).T).T+Phi.T).T
    … …  
    13601374        HbH = -np.sum(Uniq.T*np.swapaxes(np.inner(bij,Uniq),2,-1),axis=1)
    13611375        Hij = np.array([Mast*np.multiply.outer(U,U) for U in np.reshape(Uniq,(-1,3))])
    1362         Hij = np.squeeze(np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(nTwin,-1,6)))
     1376        Hij = np.reshape(np.array([G2lat.UijtoU6(Uij) for Uij in Hij]),(-1,nTwin,len(SGT),6))
    13631377        Tuij = np.where(HbH<1.,np.exp(HbH),1.0)
    1364         Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*occ
    1365         fot = (FF+FP-Bab)*Tcorr
     1378        Tcorr = (np.reshape(Tiso,Tuij.shape)*Tuij).T*Mdata*Fdata/len(SGMT)
     1379        fot = np.reshape(((FF+FP).T-Bab).T,cosp.shape)*Tcorr
    13661380        fotp = FPP*Tcorr       
    1367         fa = np.array([((FF+FP).T-Bab).T*cosp*Tcorr,-FPP*sinp*Tcorr])
    1368         fb = np.array([((FF+FP).T-Bab).T*sinp*Tcorr,FPP*cosp*Tcorr])
    1369 #        GSASIIpath.IPyBreak()
     1381        if 'T' in calcControls[hfx+'histType']: #fa,fb are 2 X blkSize X nTwin X nOps x nAtoms
     1382            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-np.reshape(FPP,sinp.shape)*sinp*Tcorr])
     1383            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,np.reshape(FPP,cosp.shape)*cosp*Tcorr])
     1384        else:
     1385            fa = np.array([np.reshape(((FF+FP).T-Bab).T,cosp.shape)*cosp*Tcorr,-FPP*sinp*Tcorr])
     1386            fb = np.array([np.reshape(((FF+FP).T-Bab).T,sinp.shape)*sinp*Tcorr,FPP*cosp*Tcorr])
    13701387        fas = np.sum(np.sum(fa,axis=-1),axis=-1)      #real sum over atoms & unique hkl array(2,nTwins)
    13711388        fbs = np.sum(np.sum(fb,axis=-1),axis=-1)      #imag sum over atoms & uniq hkl
    … …  
    13731390            fbs[0] *= 0.
    13741391            fas[1] *= 0.
    1375         fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,ntwi,nEqv,nAtoms)
     1392        fax = np.array([-fot*sinp,-fotp*cosp])   #positions array(2,nRef,ntwi,nEqv,nAtoms)
    13761393        fbx = np.array([fot*cosp,-fotp*sinp])
    13771394        #sum below is over Uniq
    1378         dfadfr = np.sum(fa/occ,axis=-2)        #array(2,ntwin,nAtom) Fdata != 0 avoids /0. problem
     1395        dfadfr = np.sum(np.sum(fa/occ,axis=-2),axis=0)        #array(2,nRef,ntwin,nAtom) Fdata != 0 avoids /0. problem
    13791396        dfadba = np.sum(-cosp*Tcorr[:,nxs],axis=1)
    1380         dfadui = np.sum(-SQfactor*fa,axis=-2)
    1381         dfadx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fax,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    1382         dfadua = np.array([np.sum(-Hij[it]*np.swapaxes(fa,-2,-1)[:,it,:,:,nxs],axis=-2) for it in range(nTwin)])
    1383         # array(nTwin,2,nAtom,3) & array(nTwin,2,nAtom,6)
     1397        dfadui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fa,axis=-2),axis=0)           
     1398        dfadx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fax[:,:,:,:,:,nxs],axis=-3),axis=0) # nRef x nTwin x nAtoms x xyz; sum on ops & A,A'
     1399        dfadua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fa[:,:,:,:,:,nxs],axis=-3),axis=0)
    13841400        if not SGData['SGInv']:
    1385             dfbdfr = np.sum(fb/occ,axis=-2)        #Fdata != 0 avoids /0. problem
     1401            dfbdfr = np.sum(np.sum(fb/occ,axis=-2),axis=0)        #Fdata != 0 avoids /0. problem
    13861402            dfadba /= 2.
    13871403            dfbdba = np.sum(-sinp*Tcorr[:,nxs],axis=1)/2.
    1388             dfbdui = np.sum(-SQfactor*fb,axis=-2)
    1389             dfbdx = np.array([np.sum(twopi*Uniq[it]*np.swapaxes(fbx,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])           
    1390             dfbdua = np.array([np.sum(-Hij[it]*np.swapaxes(fb,-2,-1)[:,it,:,:,nxs],axis=2) for it in range(nTwin)])
     1404            dfbdui = np.sum(np.sum(-SQfactor[nxs,:,nxs,nxs,nxs]*fb,axis=-2),axis=0)
     1405            dfbdx = np.sum(np.sum(twopi*Uniq[nxs,:,:,:,nxs,:]*fbx[:,:,:,:,:,nxs],axis=-3),axis=0)
     1406            dfbdua = np.sum(np.sum(-Hij[nxs,:,:,:,nxs,:]*fb[:,:,:,:,:,nxs],axis=-3),axis=0)
    13911407        else:
    13921408            dfbdfr = np.zeros_like(dfadfr)
    … …  
    13971413        SA = fas[0]+fas[1]
    13981414        SB = fbs[0]+fbs[1]
    1399         dFdfr[iref] = [2.*TwMask[it]*(SA[it]*dfadfr[0,it]+SA[it]*dfadfr[1,it]+SB[it]*dfbdfr[0,it]+SB[it]*dfbdfr[1,it])*Mdata/len(Uniq[it]) for it in range(nTwin)]
    1400         dFdx[iref] = [2.*TwMask[it]*(SA[it]*dfadx[it,0]+SA[it]*dfadx[it,1]+SB[it]*dfbdx[it,0]+SB[it]*dfbdx[it,1]) for it in range(nTwin)]
    1401         dFdui[iref] = [2.*TwMask[it]*(SA[it]*dfadui[0,it]+SA[it]*dfadui[1,it]+SB[it]*dfbdui[0,it]+SB[it]*dfbdui[1,it]) for it in range(nTwin)]
    1402         dFdua[iref] = [2.*TwMask[it]*(SA[it]*dfadua[it,0]+SA[it]*dfadua[it,1]+SB[it]*dfbdua[it,0]+SB[it]*dfbdua[it,1]) for it in range(nTwin)]
     1415#        GSASIIpath.IPyBreak()
     1416        dFdfr[iBeg:iFin] = ((2.*TwMask*SA)[:,:,nxs]*dfadfr+(2.*TwMask*SB)[:,:,nxs]*dfbdfr)*Mdata[nxs,nxs,:]/len(SGMT)
     1417        dFdx[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadx+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdx
     1418        dFdui[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs]*dfadui+(2.*TwMask*SB)[:,:,nxs]*dfbdui
     1419        dFdua[iBeg:iFin] = (2.*TwMask*SA)[:,:,nxs,nxs]*dfadua+(2.*TwMask*SB)[:,:,nxs,nxs]*dfbdua
    14031420        if SGData['SGInv']: #centrosymmetric; B=0
    1404             dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
     1421            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2
    14051422        else:               
    1406             dFdtw[iref] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
    1407         dFdbab[iref] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
    1408             fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1423            dFdtw[iBeg:iFin] = np.sum(TwMask[nxs,:]*fas,axis=0)**2+np.sum(TwMask[nxs,:]*fbs,axis=0)**2
     1424#        dFdbab[iBeg:iFin] = fas[0,:,nxs]*np.array([np.sum(dfadba*dBabdA),np.sum(-dfadba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T+ \
     1425#            fbs[0,:,nxs]*np.array([np.sum(dfbdba*dBabdA),np.sum(-dfbdba*parmDict[phfx+'BabA']*SQfactor*dBabdA)]).T
     1426        iBeg += blkSize
    14091427#        GSASIIpath.IPyBreak()
    14101428    print ' %d derivative time %.4f\r'%(len(refDict['RefList']),time.time()-time0)
    14111429    #loop over atoms - each dict entry is list of derivatives for all the reflections
    1412     for i in range(len(Mdata)):     #these all OK?
     1430    for i in range(len(Mdata)):     #these all OK
    14131431        dFdvDict[pfx+'Afrac:'+str(i)] = np.sum(dFdfr.T[i]*TwinFr[:,nxs],axis=0)
    14141432        dFdvDict[pfx+'dAx:'+str(i)] = np.sum(dFdx.T[0][i]*TwinFr[:,nxs],axis=0)
    … …  
    36673685    if im: # split to nontwin/twin versions
    36683686        if len(TwinLaw) > 1:
    3669             dFdvDict = SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3687            dFdvDict = SStructureFactorDervTw(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)    #?
    36703688        else:
    3671             dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)
     3689            dFdvDict = SStructureFactorDerv(refDict,im,G,hfx,pfx,SGData,SSGData,calcControls,parmDict)  #OK
    36723690    else:
    36733691        if len(TwinLaw) > 1:
Note: See TracChangeset for help on using the changeset viewer.