nothing much

apply_gradient_and_lj_potential.py

@@ -19,7 +19,7 @@ from sklearn.neighbors import KDTree
 
 point_mass = 1
 epsilon = 1
-r_mean = 3
+r_mean = 1.5
 capmagnitude_map = 100
 capmagnitude_lj = 400
 scale_lj = 1
@@ -28,13 +28,14 @@ delta_time = 0.05
 iterrange = 120
 slice_index = 12
 num_points_fraction = 0.1
+lj_factor = 1.5
 #point1 = tuple(np.random.randint(5,45,(1,2))[0])
 #point2 = tuple(np.random.randint(5,45,(1,2))[0])
 #point1 = (25,34)
 #point2 = (25,32)
 g = 10
 
-friction = 10
+friction = 15
 
 fx = 4
 fy = 4
@@ -90,9 +91,9 @@ class PointClass:
         vy = self.velocity.y + self.acceleration.y*dt
        # print('velocity: '+str(tuple([vx,vy])))
         self.velocity = Vector((vx,vy))
-    def position_from_velocity(self,dt,voxelsize):
-        x = ((self.position.x + self.velocity.x*dt)*voxelsize)
-        y = ((self.position.y + self.velocity.y*dt)*voxelsize)
+    def position_from_velocity(self,dt):
+        x = ((self.position.x + self.velocity.x*dt))
+        y = ((self.position.y + self.velocity.y*dt))
       #  print('position: '+str(tuple([x,y])))
         self.position = Vector((x%50,y%50))
     
@@ -124,18 +125,18 @@ def sum_array_isfinite(a):
     return sum(list_a)
         
 
-def get_acceleration_from_lj_potential(targetpoint,allpoints,voxelsize):
+def get_acceleration_from_lj_potential(targetpoint,allpoints):
     
     k = 0
     points = [x for x in allpoints if x != targetpoint]
     r = np.zeros(len(points))
     angles = np.zeros(len(points))
     for point in points:
-        r[k] = targetpoint.distance_to(point)*voxelsize
+        r[k] = targetpoint.distance_to(point)
         angles[k] = targetpoint.angle_wrt_horizontal(point) #angles in radian
         k += 1
     eps = epsilon
-    rm = r_mean*1.5
+    rm = r_mean*lj_factor
     
     v_lj = eps * ((rm/r)**12 - 2*(rm/r)**6)
     
@@ -173,9 +174,9 @@ def plot_map_and_points(map_slice,points):
     nearest_neighbors = get_nearest_neighbor_distance(points)
     for point in points:
         if nearest_neighbors[k] > r_mean:
-            plt.plot(point.position.x,point.position.y,'w'+'o')
+            plt.plot(point.position.x,point.position.y,'w'+'.')
         else:
-            plt.plot(point.position.x,point.position.y,'r'+'o')
+            plt.plot(point.position.x,point.position.y,'r'+'.')
         k += 1
 
 def plot_point_series_and_map(map_slice,points):
@@ -228,7 +229,7 @@ def number_of_too_close_neighbors(points,threshold_distance):
 
 
 
-mrc = mrcfile.open('test_map.mrc')
+mrc = mrcfile.open('test_map1.mrc')
 emmap = mrc.data
 voxelsize = 1
 
@@ -246,8 +247,8 @@ all_inside_mask = np.asarray(np.where(mask_slice==1)).T.tolist()
 gy,gx = np.gradient(map_slice)
 
 num_points = int(len(all_inside_mask) * num_points_fraction)
-#points = [PointClass(tuple([x[1],x[0]])) for x in random.sample(all_inside_mask,num_points)]
-points = [PointClass(tuple([x[1],x[0]])) for x in np.load('random_selection_3.npy')]
+points = [PointClass(tuple([x[1]*voxelsize,x[0]*voxelsize])) for x in random.sample(all_inside_mask,num_points)]
+#points = [PointClass(tuple([x[1],x[0]])) for x in np.load('random_selection_3.npy')]
 #points = [PointClass((10,20)),PointClass((12,26)),PointClass((39,29)),PointClass((35,31)),PointClass((14,21)),PointClass((15,6)),PointClass((16,13)),PointClass((11,31))]
 #points = [PointClass(point1),PointClass(point2)]
 colors = ['b','g','r','c','m','y','k','w']
@@ -277,7 +278,7 @@ for iter in range(iterrange):
 
     for point in points:
         gradient_acceleration,map_potential = get_acceleration_from_gradient(map_slice, g, point=point, voxelsize=voxelsize)
-        lj_potential_acceleration,lj_potential = get_acceleration_from_lj_potential(point, points,voxelsize=voxelsize)
+        lj_potential_acceleration,lj_potential = get_acceleration_from_lj_potential(point, points)
         
         gradient_acceleration,lj_potential_acceleration = gradient_acceleration.scale(scale_map),lj_potential_acceleration.scale(scale_lj)
         acceleration = add_Vector(gradient_acceleration,lj_potential_acceleration)
@@ -297,7 +298,7 @@ for iter in range(iterrange):
     #plt.imshow(map_slice)
     #plot_point_series_and_map(map_slice,points)
     plot_map_and_points(map_slice,points)
-    '''
+    
     plt.title('Red: LJ, Black: Gradient, Yellow: Net, Iter = '+str(iter))
     for point in points:
         (X,Y) = point.position.get()
@@ -305,11 +306,11 @@ for iter in range(iterrange):
         (U_map,V_map) = point.map_acceleration.get()
         (U,V) = point.acceleration.get()
 
-        plt.quiver(X,Y,U_lj,-V_lj,scale=0.5,scale_units='x',color='r')
-        plt.quiver(X,Y,U_map,-V_map,scale=0.5,scale_units='x',color='k')
-        plt.quiver(X,Y,U,-V,scale=0.5,scale_units='x',color='y')
+        plt.quiver(X,Y,U_lj,-V_lj,scale=8,scale_units='x',color='r')
+       # plt.quiver(X,Y,U_map,-V_map,scale=2,scale_units='x',color='k')
+       # plt.quiver(X,Y,U,-V,scale=2,scale_units='x',color='y')
     #print(points[0].history['pos'])
-    '''
+    
     #plot_point_series_and_map(map_slice,points)    
     #plt.show()
     camera.snap()    
@@ -323,7 +324,7 @@ for iter in range(iterrange):
     for point in points:
         point.velocity_from_acceleration(dt)        
         #print(str(iter)+", "+str(k)+", "+str(point.lj_acceleration.magnitude())+", "+str(point.map_acceleration.magnitude())+", "+str(point.acceleration.magnitude()))
-        point.position_from_velocity(dt, voxelsize)
+        point.position_from_velocity(dt)
         point.update_history()
         
         # Following stuff for analysis
@@ -334,7 +335,7 @@ for iter in range(iterrange):
    
 df = pd.DataFrame(point_analysis)
 
-outfilename = 'gradient_and_lj_potential_3_int'
+outfilename = 'gradient_and_lj_potential'
 animation = camera.animate()
 animation.save(outfilename+'.gif')
 

study_iterations.py

@@ -19,7 +19,7 @@ from sklearn.neighbors import KDTree
 
 point_mass = 1
 epsilon = 1
-r_mean = 3
+r_mean = 1.5
 capmagnitude_map = 100
 capmagnitude_lj = 400
 scale_lj = 1
@@ -223,7 +223,7 @@ def number_of_too_close_neighbors(points,threshold_distance):
 
 mrc = mrcfile.open('test_map.mrc')
 emmap = mrc.data
-voxelsize = 1
+voxelsize = 1#mrc.voxel_size.x
 
 map_slice = normalize(emmap[:,:,slice_index],0,100)