Changeset 4099 for trunk/G2shapes.py


Ignore:
Timestamp:
Aug 18, 2019 6:06:56 PM (2 years ago)
Author:
vondreele
Message:

more speedup of SHAPES

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/G2shapes.py

    r4098 r4099  
    800800        return dr;
    801801   
     802    # Return non-zero distances one coordinate & the rest
     803   
     804    def get_drs(xyz,XYZ):
     805       
     806        return  np.sqrt(np.sum((XYZ-xyz)**2,axis=1))
     807   
    802808    # Find center of beads within a radii
    803809   
     
    806812        num_beads = len(aList_beads_x)
    807813   
    808         xsum = 0.0
    809         ysum = 0.0
    810         zsum = 0.0
    811         count_beads = 0.0
    812    
    813         i = 0
    814         while i < num_beads:
     814#        xsum = 0.0
     815#        ysum = 0.0
     816#        zsum = 0.0
     817#        count_beads = 0.0
     818#   
     819#        i = 0
     820#        while i < num_beads:
     821#           
     822#            dr = get_dr(x,y,z,aList_beads_x[i],aList_beads_y[i],aList_beads_z[i])
     823#   
     824#            if dr > radii_1 and dr < radii_2:
     825#                count_beads = count_beads + 1.0
     826#                xsum = xsum + float(aList_beads_x[i])
     827#                ysum = ysum + float(aList_beads_y[i])
     828#                zsum = zsum + float(aList_beads_z[i])           
     829#   
     830#            i = i + 1
     831#   
     832#        if count_beads > 0.0:
     833#            xsum = xsum/count_beads
     834#            ysum = ysum/count_beads
     835#            zsum = zsum/count_beads
     836#            delta = (xsum - x)**2 + (ysum - y)**2 + (zsum - z)**2
     837#            delta = math.sqrt(delta)
     838#        else:
     839#            delta = 0.0
    815840           
    816             dr = get_dr(x,y,z,aList_beads_x[i],aList_beads_y[i],aList_beads_z[i])
    817    
    818             if dr > radii_1 and dr < radii_2:
    819                 count_beads = count_beads + 1.0
    820                 xsum = xsum + float(aList_beads_x[i])
    821                 ysum = ysum + float(aList_beads_y[i])
    822                 zsum = zsum + float(aList_beads_z[i])           
    823    
    824             i = i + 1
    825    
    826         if count_beads > 0.0:
    827             xsum = xsum/count_beads
    828             ysum = ysum/count_beads
    829             zsum = zsum/count_beads
    830             delta = (xsum - x)**2 + (ysum - y)**2 + (zsum - z)**2
    831             delta = math.sqrt(delta)
    832         else:
    833             delta = 0.0
    834        
     841        XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T
     842        xyz = np.array([x,y,z])
     843        drs = get_drs(xyz,XYZ)
     844        sumXYZ = np.array([XYZ[i] for i in range(num_beads) if radii_1<drs[i]<radii_2])
     845        count_beads = sumXYZ.shape[0]
     846       
     847        delta = 0.0
     848        if count_beads:
     849            delta = np.sqrt(np.sum(((np.sum(sumXYZ,axis=0)/count_beads)-xyz)**2))
     850           
    835851        return delta;
    836852   
     
    10271043        num_hist = len(aList_r)
    10281044        max_dr = (float(num_hist)-1.0)*hist_grid
    1029         num_beads = len(aList_beads_x)
    1030    
    1031         i = 0
    1032         while i < num_hist:
    1033             aList_pr_model_test2[i] = 0.0
    1034             i = i + 1
    1035    
    1036         i = 0
    1037         while i < num_beads:
    1038    
    1039             if i != ii:
    1040                 x2 = float(aList_beads_x[i])
    1041                 y2 = float(aList_beads_y[i])
    1042                 z2 = float(aList_beads_z[i])
    1043                 dr = get_dr(x1,y1,z1,x2,y2,z2)
    1044                 dr = min(dr,max_dr)
    1045                 dr_grid = dr/hist_grid
    1046                 int_dr_grid = int(dr_grid)
    1047                 int_dr_grid = max(int_dr_grid,0)
    1048                 int_dr_grid = min(int_dr_grid,num_hist-2)               
    1049                 ip_low = int_dr_grid
    1050                 ip_high = ip_low + 1           
    1051                 ip_high_frac = dr_grid - float(int_dr_grid)
    1052                 ip_low_frac = 1.0 - ip_high_frac
    1053                 aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac
    1054                 aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac
    1055    
    1056             i = i + 1
    1057    
     1045#        num_beads = len(aList_beads_x)
     1046       
     1047        aList_pr_model_test2 = num_hist*[0.0,]
     1048       
     1049        XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T
     1050        xyz = np.array([x1,y1,z1])
     1051        drs = get_drs(xyz,XYZ)
     1052        drs_grid = np.where(drs>max_dr,max_dr,drs)/hist_grid
     1053        int_drs_grid = np.array(drs_grid,dtype=np.int)
     1054        int_drs_grid = np.where(int_drs_grid>num_hist-2,num_hist-2,int_drs_grid)               
     1055        ip_lows = int_drs_grid
     1056        ip_highs = ip_lows + 1           
     1057        ip_high_fracs = drs_grid - int_drs_grid
     1058        ip_low_fracs = 1.0 - ip_high_fracs
     1059        for ip_low in ip_lows:
     1060            aList_pr_model_test2[ip_low] += ip_low_fracs[ip_low]
     1061        for ip_high in ip_highs:
     1062            aList_pr_model_test2[ip_high] += ip_high_fracs[ip_high]
     1063   
     1064#        i = 0
     1065#        while i < num_hist:
     1066#            aList_pr_model_test2[i] = 0.0
     1067#            i = i + 1
     1068#   
     1069#        i = 0
     1070#        while i < num_beads:
     1071#   
     1072#            if i != ii:
     1073#                x2 = float(aList_beads_x[i])
     1074#                y2 = float(aList_beads_y[i])
     1075#                z2 = float(aList_beads_z[i])
     1076#                dr = get_dr(x1,y1,z1,x2,y2,z2)
     1077#                dr = min(dr,max_dr)
     1078#                dr_grid = dr/hist_grid
     1079#                int_dr_grid = int(dr_grid)
     1080#                int_dr_grid = max(int_dr_grid,0)
     1081#                int_dr_grid = min(int_dr_grid,num_hist-2)               
     1082#                ip_low = int_dr_grid
     1083#                ip_high = ip_low + 1           
     1084#                ip_high_frac = dr_grid - float(int_dr_grid)
     1085#                ip_low_frac = 1.0 - ip_high_frac
     1086#                aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac
     1087#                aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac
     1088#   
     1089#            i = i + 1
     1090#   
    10581091        return;
    10591092   
Note: See TracChangeset for help on using the changeset viewer.