Changeset 4098 for trunk/G2shapes.py


Ignore:
Timestamp:
Aug 16, 2019 3:57:43 PM (2 years ago)
Author:
vondreele
Message:

speedup of SHAPES

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/G2shapes.py

    r4097 r4098  
    386386   
    387387        dmax_over2 = dmax/2.0
    388         num_box = int(dmax/box_step)+1
     388        num_box = int(dmax/box_step)
     389        search_sq = rsearch**2
    389390       
    390         box_grid = np.zeros((num_box,num_box,num_box),dtype=int)
    391         aList_beads_xyz_loc = np.array([(np.array(aList_beads_x)+dmax_over2)/box_step,
    392             (np.array(aList_beads_y)+dmax_over2)/box_step,
    393             (np.array(aList_beads_z)+dmax_over2)/box_step],dtype=int)
    394    
    395         for ix,iy,iz in aList_beads_xyz_loc.T:
    396             box_grid[ix,iy,iz] += 1
    397            
    398         non_zero = np.argwhere(box_grid)
    399         for ix,iy,iz in non_zero:
    400             aList_box_x_all.append((ix*box_step)-dmax_over2)
    401             aList_box_y_all.append((iy*box_step)-dmax_over2)
    402             aList_box_z_all.append((iz*box_step)-dmax_over2)
    403             aList_box_score.append(box_grid[ix,iy,iz])
    404    
     391        XYZ = np.meshgrid(np.linspace(-dmax_over2,dmax_over2,num_box),
     392            np.linspace(-dmax_over2,dmax_over2,num_box),
     393            np.linspace(-dmax_over2,dmax_over2,num_box))
     394        XYZ = np.array([XYZ[0].flatten(),XYZ[1].flatten(),XYZ[2].flatten()]).T
     395        xyz = np.array((aList_beads_y,aList_beads_x,aList_beads_z)).T
     396        for XYZi in XYZ:
     397            dsq = np.sum((xyz-XYZi)**2,axis=1)
     398            count = int(np.sum(np.array([1 for dist in dsq if dist < search_sq])))
     399            if count>1:
     400                aList_box_x_all.append(XYZi[0])
     401                aList_box_y_all.append(XYZi[1])
     402                aList_box_z_all.append(XYZi[2])
     403                aList_box_score.append(count)
    405404        return;
    406405   
Note: See TracChangeset for help on using the changeset viewer.