Changeset 3896 for trunk/GSASIImath.py


Ignore:
Timestamp:
Apr 14, 2019 9:23:46 AM (3 years ago)
Author:
vondreele
Message:

enormous speed improvement for contour plot - not dependent on map size

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIImath.py

    r3894 r3896  
    36573657    :returns: density at xyz
    36583658    '''
     3659    def getBoxes(rho,I):
     3660        Rhos = np.zeros((2,2,2))
     3661        Mx,My,Mz = rho.shape
     3662        Ix,Iy,Iz = I
     3663        Rhos = np.array([[[rho[Ix%Mx,Iy%My,Iz%Mz],rho[Ix%Mx,Iy%My,(Iz+1)%Mz]],
     3664                          [rho[Ix%Mx,(Iy+1)%My,Iz%Mz],rho[Ix%Mx,(Iy+1)%My,(Iz+1)%Mz]]],
     3665                        [[rho[(Ix+1)%Mx,Iy%My,Iz%Mz],rho[(Ix+1)%Mx,Iy%My,(Iz+1)%Mz]],
     3666                          [rho[(Ix+1)%Mx,(Iy+1)%My,Iz%Mz],rho[(Ix+1)%Mx,(Iy+1)%My,(Iz+1)%Mz]]]])
     3667        return Rhos
     3668       
    36593669    Blk = 8     #8 seems optimal
    36603670    nBlk = len(XYZ)//Blk        #select Blk so this is an exact divide
    3661     rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)[:2,:2,:2]
    36623671    mapShape = np.array(rho.shape)
    36633672    mapStep = 1./mapShape
     
    36653674    iBeg = 0
    36663675    R = np.zeros(len(XYZ))
    3667 #actually a bit faster!
     3676#actually a lot faster!
    36683677    for iblk in range(nBlk):
    36693678        iFin = iBeg+Blk
    36703679        Xs = X[iBeg:iFin]
    36713680        I = np.array(np.rint(Xs*mapShape),dtype='int')
    3672         Rhos = np.array([rollMap(rho,-i) for i in I])
     3681        Rhos = np.array([getBoxes(rho,i) for i in I])
    36733682        Ds = Xs-I*mapStep
    36743683        RIJs = Rhos[:,0,:2,:2]*(1.-Ds[:,0][:,nxs,nxs])
     
    36763685        R[iBeg:iFin] = RIs[:,0]*(1.-Ds[:,2])+RIs[:,1]*Ds[:,2]
    36773686        iBeg += Blk
    3678 #one at a time
    3679 #    for i,x in enumerate(X):
    3680 #        D = x-I[i]*mapStep         #position inside map cell
    3681 #        Rho = rollMap(rho,-I[i])    #shifts map so point is in corner
    3682 #        RIJ = Rho[0,:2,:2]*(1.-D[0])
    3683 #        RI = RIJ[0]*(1.-D[1])+RIJ[1]*D[1]
    3684 #        R[i] = RI[0]*(1.-D[2])+RI[1]*D[2]
    36853687    return R
    36863688       
Note: See TracChangeset for help on using the changeset viewer.