mean long axis with error bars

convergence_extension/average_over_clutches.py

+import numpy as np
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+sns.set_theme(style="white")
+# Number of samples in this clutch
+N1 = 10
+N2 = 11
+N3 = 4
+
+########## Load the data ###########################################################
+l_clutch2 = np.zeros(N2, dtype=object)  
+time2 = np.zeros(N2, dtype=object)   #
+for j in range(N1, N2+N1):
+    l_clutch2[j-N1-1] = np.load('long_axis_xy_coordinates_movie' + str(j) + '.npy')
+    time2[j-N1-1] = np.load('time_xy_coordinates_movie' + str(j) + '.npy')
+    
+#NOTE: Clutch1 is movie 1-10
+l_clutch1 = np.zeros(N1, dtype=object)  #l[i] is the long axis for movie i
+time1 = np.zeros(N1, dtype=object)   #time[i] is the time for movie i
+for i in range(1, N1+1):
+    l_clutch1[i-1] = np.load('long_axis_xy_coordinates_movie' + str(i) + '.npy')
+    time1[i-1] = np.load('time_xy_coordinates_movie' + str(i) + '.npy')
+
+#NOTE: Clutch2 is movie 11-21
+l_clutch2 = np.zeros(N2, dtype=object)  
+time2 = np.zeros(N2, dtype=object)   #
+for j in range(N1, N2+N1):
+    l_clutch2[j-N1-1] = np.load('long_axis_xy_coordinates_movie' + str(j) + '.npy')
+    time2[j-N1-1] = np.load('time_xy_coordinates_movie' + str(j) + '.npy')
+
+#NOTE: Clutch3 is movie 22, 23, 25, 26
+l_clutch3 = np.zeros(N3, dtype=object)
+time3 = np.zeros(N3, dtype=object)
+l_clutch3[0] = np.load('long_axis_xy_coordinates_movie22.npy')
+time3[0] = np.load('time_xy_coordinates_movie22.npy')
+l_clutch3[1] = np.load('long_axis_xy_coordinates_movie23.npy')
+time3[1] = np.load('time_xy_coordinates_movie23.npy')
+l_clutch3[2] = np.load('long_axis_xy_coordinates_movie25.npy')
+time3[2] = np.load('time_xy_coordinates_movie25.npy')
+l_clutch3[3] = np.load('long_axis_xy_coordinates_movie26.npy')
+time3[3] = np.load('time_xy_coordinates_movie26.npy')
+
+
+
+######### make all arrays the same length ##########################
+min_length1 = min([len(x) for x in l_clutch1])
+for i in range(N1):
+    l_clutch1[i] = l_clutch1[i][:min_length1]
+    time1[i] = time1[i][:min_length1]
+
+min_length2 = min([len(x) for x in l_clutch2])
+for i in range(N2):
+    l_clutch2[i] = l_clutch2[i][:min_length2]
+    time2[i] = time2[i][:min_length2]
+
+min_length3 = min([len(x) for x in l_clutch3])
+for i in range(N3):
+    l_clutch3[i] = l_clutch3[i][:min_length3]
+    time3[i] = time3[i][:min_length3]
+
+# All time arrays are the same. Randomly choose one to use as the time array
+time1 = time1[0]
+time2 = time2[0]
+time3 = time3[0]
+
+######### Plot the average over each clutch ############################
+avg_l1 = np.mean(l_clutch1, axis=0)
+std_l1 = np.std(l_clutch1, axis=0)
+avg_l2 = np.mean(l_clutch2, axis=0)
+std_l2 = np.std(l_clutch2, axis=0)
+avg_l3 = np.mean(l_clutch3, axis=0)
+std_l3 = np.std(l_clutch3, axis=0)
+plt.plot(time1, avg_l1, linewidth=2, color='#e07a5f', label='Clutch 1')
+plt.fill_between(time1, avg_l1-std_l1, avg_l1+std_l1, alpha=0.1, color='#e07a5f')
+plt.plot(time2, avg_l2, linewidth=2, color='#3d405b', label='Clutch 2')
+plt.fill_between(time2, avg_l2-std_l2, avg_l2+std_l2, alpha=0.1, color='#3d405b')
+plt.plot(time3, avg_l3, linewidth=2, color='#81b29a', label='Clutch 3')
+plt.fill_between(time3, avg_l3-std_l3, avg_l3+std_l3, alpha=0.1, color='#81b29a')
+plt.xlabel('Time (min)')
+plt.ylabel('Long axis ($\mu m$)')
+plt.legend()
+plt.savefig('Average_over_each_clutch.pdf')
+
+######## Average over 3 clutches #####
+min_length = min([len(avg_l1), len(avg_l2)])
+avg_l1 = avg_l1[:min_length]
+avg_l2 = avg_l2[:min_length]
+std_l1 = std_l1[:min_length]
+std_l2 = std_l2[:min_length]
+avg_over_clutches = np.mean(np.array([avg_l1, avg_l2]), axis=0)
+std_over_clutches = np.std(np.array([avg_l1, avg_l2]), axis=0)
+fig, ax = plt.subplots()
+ax.plot(time1[:min_length], avg_over_clutches, linewidth=2, color='black')
+ax.fill_between(time1[:min_length], avg_over_clutches-std_over_clutches, avg_over_clutches+std_over_clutches, alpha=0.1, color='black')
+ax.set_xlabel('Time (min)')
+ax.set_ylabel('Long axis ($\mu m$)')
+fig.savefig('Average_over_all_clutches.pdf')

convergence_extension/perimeter_long_axis.py

@@ -3,6 +3,9 @@ import matplotlib.pyplot as plt
 import pandas as pd
 from matplotlib.widgets import Slider
 from scipy.spatial import ConvexHull
+from matplotlib.animation import FuncAnimation
+import os
+from tempfile import TemporaryDirectory
 
 def perimeter_area_majoraxis(X, Y):
 
@@ -35,8 +38,8 @@ def perimeter_area_majoraxis(X, Y):
 # hull = ConvexHull(points)
 # print(perimeter_area_majoraxis(x, y))
 # raise SystemExit
-
-file = np.loadtxt('xy_coordinates_movie1.txt')
+filename = 'xy_coordinates_movie1.txt'
+file = np.loadtxt(filename)
 number_of_pixels = file.shape[0]
 time_points = file.shape[1]
 time_array = np.arange(0, 5*int(time_points/2), 5)
@@ -58,6 +61,10 @@ y = np.array(y.interpolate(method='linear', axis=0, limit_direction='both'))
 x = x.T
 y = y.T
 
+# NOTE: This is only for movie11 -- remaining NANs
+# x=x[1:]
+# y=y[1:]
+
 # make movie of the contours and convex hull
 fig, ax = plt.subplots(1, 1, sharex=True, figsize=(8,8))
 ax.set_xlabel(r'$x$')
@@ -66,9 +73,9 @@ ax.set_ylim(np.min(y), np.max(y))
 
 ax.set_ylabel(r'$y$')
 
-outlineplot, = ax.plot(x[0], y[0], marker = '.', linestyle='None')
+outlineplot, = ax.plot(x[0], y[0], marker = '.', color='darkgoldenrod', linestyle='None')
 hull = perimeter_area_majoraxis(x[0], y[0])[-1]
-hull_lines = [ax.plot(x[0][simplex], y[0][simplex], 'r-')[0] for simplex in hull.simplices]
+hull_lines = [ax.plot(x[0][simplex], y[0][simplex], '-', color='black')[0] for simplex in hull.simplices]
 
 def update(value):
     ti = np.abs(time_array - value).argmin()
@@ -83,7 +90,7 @@ def update(value):
             line.remove()
 
     hull = perimeter_area_majoraxis(x[ti], y[ti])[-1]
-    update.hull_lines = [ax.plot(x[ti][simplex], y[ti][simplex], 'r-')[0] for simplex in hull.simplices]
+    update.hull_lines = [ax.plot(x[ti][simplex], y[ti][simplex], '-', color='black')[0] for simplex in hull.simplices]
     plt.draw()
         
 sax = plt.axes([0.1, 0.92, 0.7, 0.02])
@@ -92,9 +99,28 @@ slider = Slider(sax, r'$t/\tau$', min(time_array), max(time_array),
             fc='#999999')
 slider.drawon = False
 slider.on_changed(update)
+
 plt.show()
 
+# save movie 
+print('Saving movie for the contour and hull..')
+FPS = 20
+fps = float(FPS)
+command = "ffmpeg -y -r"
+options = "-b:v 3600k -qscale:v 4 -vcodec mpeg4"
+tmp_dir = TemporaryDirectory()
+get_filename = lambda x: os.path.join(tmp_dir.name, x)
+for tt in range(len(time_array)):
+    slider.set_val(time_array[tt])
+    fr = get_filename("%03d.png" % tt)
+    fig.savefig(fr, facecolor=fig.get_facecolor(), dpi=100)
+os.system(command + " " + str(fps)
+            + " -i " + tmp_dir.name + os.sep
+            + "%03d.png " + options + " " + 'convex_hull_contour.mp4')
+tmp_dir.cleanup()
 # find the perimeter and area
+#NOTE: For movie 11, time_array_mov11 = time_array[1:] and reduce the length of the following arrays by 1
+
 perimeter_array = np.zeros(int(time_points/2))
 area_array = np.zeros(int(time_points/2))
 long_axis_array = np.zeros(int(time_points/2))
@@ -102,20 +128,29 @@ long_axis_array = np.zeros(int(time_points/2))
 for k in range(0,  int(time_points/2)):
     perimeter_array[k], area_array[k], long_axis_array[k], _ = perimeter_area_majoraxis(x[k], y[k])   
      
+np.save('long_axis_%s.npy' %filename, long_axis_array)
+np.save('time_%s.npy' %filename, time_array)
+
+
 plt.plot(time_array, perimeter_array, marker='o', linestyle='None')     
 plt.xlabel('Time (min)')
 plt.ylabel('Perimeter of the convex hull ($\mu m$)')
 plt.show()
+
+
 plt.plot(time_array, long_axis_array, marker='o', linestyle='None')     
 plt.xlabel('Time (min)')
 plt.ylabel('Long axis of the convex hull ($\mu m$)')
-coefficients = np.polyfit(time_array[25:-13], long_axis_array[25:-13], 1)
-polyomial = np.poly1d(coefficients)
-xline = time_array[25:-13]
-yline = polyomial(xline)
-plt.plot(xline, yline)
-print(coefficients[0])
+
+# coefficients = np.polyfit(time_array[25:-13], long_axis_array[25:-13], 1)
+# polyomial = np.poly1d(coefficients)
+# xline = time_array[25:-13]
+# yline = polyomial(xline)
+# plt.plot(xline, yline)
+# print(coefficients[0])
 plt.show()
+
+
 plt.plot(time_array, area_array, marker='o', linestyle='None')     
 plt.xlabel('Time (min)')
 plt.ylabel('Area of the convex hull ($\mu m^2$)')