In [10]:
"""
This notebook will be used as a template for each dataset giving the key outputs I need for presentations
and publications
"""
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
sys.path.append(module_path)
from trajectory_visualization import plot_trajectory, sidebyside, shift_trajectory, overlay, shift_trajectory3D
from trajectory_visualization import plot_trajectories3D, plot_3Doverlay, plot_MSDorDeff, plot_MeanMSDorDeff, randtraj, multrandtraj
from trajectory_visualization import randtraj2, plot_Mean2DMSDsorDeff, plot_MSDorDeffLR, LRfor3D2D, fillin, prettify
from methods_of_Dcalc import Dpointder, Dlinalg, Dlinalgw
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.optimize as opt
from mpl_toolkits.mplot3d import Axes3D
from operator import itemgetter
import random
import numpy as np
import numpy.linalg as la
from mpl_toolkits.mplot3d import Axes3D
pi = np.pi
sin = np.sin
cos = np.cos
In [16]:
# Type in datasets you would like to combine. Adjust totvids to reflect the total number of videos. Also type in the
# umppx and fps in the conversion dictionary.
logplot = 'gopher_S1_R1_logplot'
Mplot = 'gopher_S1_R1_Mplot'
Dplot = 'gopher_S1_R1_Dplot'
Hplot = 'gopher_S1_R1_Hplot'
Hlogplot = 'gopher_S1_R1_Hlogplot'
# Chad sets up conversion factors here, uses these to convert distance from pixels to microns and frames to seconds
conversion = dict() #First element is the umppx, second is fps
conversion[1] = (0.11, 4, 1)
conversion[2] = (0.11, 4, 1)
conversion[3] = (0.11, 4, 1)
conversion[4] = (0.11, 4, 1)
# Importing all trajectories obtained from all videos taken on this specific sample, the size of this
# trajectory dictionary will fluctuate depending on how many vids are taken on a given sample
trajectory = dict()
trajectory[1] = np.genfromtxt('../../sample_data/Regional_Data/gopher_S1_R1_V1_10NPs_traj.csv',
delimiter =",")
trajectory[2] = np.genfromtxt('../../sample_data/Regional_Data/gopher_S1_R1_V2_10NPs_traj.csv',
delimiter =",")
trajectory[3] = np.genfromtxt('../../sample_data/Regional_Data/gopher_S1_R1_V3_10NPs_traj.csv',
delimiter =",")
trajectory[4] = np.genfromtxt('../../sample_data/Regional_Data/gopher_S1_R1_V4_10NPs_traj.csv',
delimiter =",")
totvids = 4
# for loop that goes through all data tables in all dictionary entries and removes the titles of each column
# and removes the labels of each row. This leaves only the hard data.
for num in range(1, totvids + 1):
#Remove titles from columns
trajectory[num]=np.delete(trajectory[num], 0,0)
#Remove the number column from dataset
trajectory[num]=np.delete(trajectory[num],0,1)
#Get rid of trajectories that are too short, as determined by prettify
final = dict()
tots = dict()
for num in range(1, totvids + 1):
(final[num], tots[num]) = prettify(trajectory[num], 40, 40, conversion[num][0], conversion[num][1], conversion[num][2])
#Collect trajectories of separate videos into one array
parts = dict()
lehi = dict()
for num in range(1, totvids + 1):
parts[num] = tots[num]
lehi[num] = final[num][1]
counter = 1
while counter < parts[num]:
counter = counter + 1
lehi[num] = np.append(lehi[num], final[num][counter], axis=0)
for num in range(2, totvids + 1):
lehi[num][:,0] = lehi[num][:,0] + max(lehi[num-1][:, 0])
nephi = lehi[1]
counter = 1
while counter < totvids:
counter = counter + 1
nephi = np.append(nephi, lehi[counter], axis=0)
In [3]:
D = dict()
SD = dict()
total = dict()
t = dict()
D['pd2D'], SD['pd2D'], total['pd2D'], t['pd2D'] = Dpointder(nephi, 0, 15, 17, 1, '2D')
D['pd1Dx'], SD['pd1Dx'], total['pd1Dx'], t['pd1Dx'] = Dpointder(nephi, 0, 15, 17, 1, '1Dx')
D['pd1Dy'], SD['pd1Dy'], total['pd1Dy'], t['pd1Dy'] = Dpointder(nephi, 0, 15, 17, 1, '1Dy')
print(D['pd2D'], SD['pd2D'])
print(D['pd1Dx'], SD['pd1Dx'])
print(D['pd1Dy'], SD['pd1Dy'])
print(total['pd2D'])
In [18]:
def plot_Mean1DMSDsorDeffLOG(traj, n1, n2, n3, dec1, dec2, datatype, filename, limit1, limit2, tick1, tick2):
"""
Plots the MSDs or Deffs from a trajectory dataset.
n1: particle numbers (typically 0)
n2: time (typically 1)
n3: MSDs or Deffs (11 or 15 typically)
"""
# Creates an array 'particles' that contains the particle number at each frame.
particles = traj[:, n1]
total = int(max(particles))
total1 = total + 1
rawtime = traj[:, n2]
bow = traj.shape[0]
raw2DMSDs = np.zeros((bow, 4))
raw2DMSDs[:, 0] = traj[:, n3]
raw2DMSDs[:, 1:4] = traj[:, n3 + 3:n3 + 6]
MSD = dict()
time = dict()
# Creates an array for each trajectory containing all xyz data
for num in range(1, total1):
hold = np.where(particles == num)
itindex = hold[0]
min1 = min(itindex)
max1 = max(itindex)
MSD[num] = (raw2DMSDs[min1+2:max1, :])
time[num] = (rawtime[min1+2:max1])
MMSD = MSD[1]
for num in range(2, total1):
MMSD = MMSD + MSD[num]
MMSD = MMSD/total1
SD = np.zeros(np.shape(MMSD))
#Now to calculate the standard dev at each point:
for num in range (1, total1):
SDunit = (MSD[num] - MMSD)**2
SD = SD + SDunit
SD = np.sqrt(SD/total1)
SE = SD/np.sqrt(total1)
#Linear algebra to find Deff:
t = time[1][:]
w = dict()
line = dict()
A = np.ones((np.shape(t)[0], 2))
A[:, 0] = t
w[0] = np.linalg.lstsq(A, MMSD[:, 0])[0]
w[1] = np.linalg.lstsq(A, MMSD[:, 1])[0]
w[2] = np.linalg.lstsq(A, MMSD[:, 2])[0]
line[0] = w[0][0]*t + w[0][1]
line[1] = w[1][0]*t + w[1][1]
line[2] = w[2][0]*t + w[2][1]
#Linear algebra for fit on log plot:
wl = dict()
linel = dict()
lt = np.log(t)
lA = np.ones((np.shape(t)[0], 2))
lA[:, 0] = lt
lM = np.log(MMSD)
wl[0] = np.linalg.lstsq(lA, lM[:, 0])[0]
wl[1] = np.linalg.lstsq(lA, lM[:, 1])[0]
wl[2] = np.linalg.lstsq(lA, lM[:, 2])[0]
linel[0] = np.exp(wl[0][0]*lt + wl[0][1])
linel[1] = np.exp(wl[1][0]*lt + wl[1][1])
linel[2] = np.exp(wl[2][0]*lt + wl[2][1])
# Creates figure
fig = plt.figure(figsize=(24, 18), dpi=80)
ax = fig.add_subplot(111)
# ax.set_title('Particle Trajectories', x=0.5, y=1.15)
ax.plot(t, MMSD[:, 0], linestyle='', linewidth=10, label='2D', marker='o', ms=10, color='blue')
ax.plot(t, MMSD[:, 1], linestyle='', linewidth=10, label='1D x', marker='o', ms=10, color='red')
ax.plot(t, MMSD[:, 2], linestyle='', linewidth=10, label='1D x', marker='o', ms=10, color='green')
ax.errorbar(t, MMSD[:, 0], yerr=SE[:, 0], fmt='', linestyle='', capsize=7, capthick=2, elinewidth=2, color='blue')
ax.errorbar(t, MMSD[:, 1], yerr=SE[:, 1], fmt='', linestyle='', capsize=7, capthick=2, elinewidth=2, color='red')
ax.errorbar(t, MMSD[:, 2], yerr=SE[:, 2], fmt='', linestyle='', capsize=7, capthick=2, elinewidth=2, color='green')
#ax.plot(t, line[0], linewidth=3, color='blue')
#ax.plot(t, line[1], linewidth=3, color='red')
#ax.plot(t, line[2], linewidth=3, color='green')
ax.plot(t, linel[0], linewidth=3, color='blue')
ax.plot(t, linel[1], linewidth=3, color='red')
ax.plot(t, linel[2], linewidth=3, color='green')
hfont = {'fontname':'Arial'}
#mpl.rc('font', **{'family':'serif','serif':['monbaiti']})
# A few adjustments to prettify the graph
for item in ([ax.xaxis.label, ax.yaxis.label] +
ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(90)
xmajor_ticks = np.arange(0, limit1, tick1)
ymajor_ticks = np.arange(0, limit2, tick2)
ax.set_xticks(xmajor_ticks)
ax.set_yticks(ymajor_ticks)
ax.title.set_fontsize(70)
ax.set_xlabel('Time (s)',fontsize=115, **hfont)
ax.set_ylabel(r'MSD ($\mu$m$^2$)', fontsize=115, **hfont)
ax.tick_params(direction='out', pad=16)
ax.legend(loc=(0.05, 0.60), prop={'size': 65})
plt.gca().xaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%.{}f'.format(dec1)))
plt.gca().yaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%.{}f'.format(dec2)))
plt.yscale('log')
plt.xscale('log')
plt.gca().set_xlim([0.1, limit1])
plt.gca().set_ylim([0.01, limit2])
# Save your figure
plt.savefig('{}.png'.format(filename), bbox_inches='tight')
return MMSD, total1
In [19]:
one1, total1 = plot_Mean1DMSDsorDeffLOG(nephi, 0, 15, 9, 1, 1, 'MSD (um^2)', logplot, 10.1, 1000.1, 1, 1)
plt.show()
In [20]:
nephi
Out[20]:
In [6]:
def plot_Mean1DMSDsorDeff(traj, n1, n2, n3, dec1, dec2, datatype, filename, limit1, limit2, tick1, tick2):
"""
Plots the MSDs or Deffs from a trajectory dataset.
n1: particle numbers (typically 0)
n2: time (typically 1)
n3: MSDs or Deffs (11 or 15 typically)
"""
# Creates an array 'particles' that contains the particle number at each frame.
particles = traj[:, n1]
total = int(max(particles))
total1 = total + 1
rawtime = traj[:, n2]
bow = traj.shape[0]
raw2DMSDs = np.zeros((bow, 4))
raw2DMSDs[:, 0] = traj[:, n3]
raw2DMSDs[:, 1:4] = traj[:, n3 + 3:n3 + 6]
MSD = dict()
time = dict()
# Creates an array for each trajectory containing all xyz data
for num in range(1, total1):
hold = np.where(particles == num)
itindex = hold[0]
min1 = min(itindex)
max1 = max(itindex)
MSD[num] = (raw2DMSDs[min1+2:max1, :])
time[num] = (rawtime[min1+2:max1])
MMSD = MSD[1]
for num in range(2, total1):
MMSD = MMSD + MSD[num]
MMSD = MMSD/total1
SD = np.zeros(np.shape(MMSD))
#Now to calculate the standard dev at each point:
for num in range (1, total1):
SDunit = (MSD[num] - MMSD)**2
SD = SD + SDunit
SD = np.sqrt(SD/total1)
SE = SD/np.sqrt(total1)
#Linear algebra to find Deff:
t = time[1][:]
w = dict()
line = dict()
A = np.ones((np.shape(t)[0], 2))
A[:, 0] = t
w[0] = np.linalg.lstsq(A, MMSD[:, 0])[0]
w[1] = np.linalg.lstsq(A, MMSD[:, 1])[0]
w[2] = np.linalg.lstsq(A, MMSD[:, 2])[0]
line[0] = w[0][0]*t + w[0][1]
line[1] = w[1][0]*t + w[1][1]
line[2] = w[2][0]*t + w[2][1]
#Linear algebra for fit on log plot:
wl = dict()
linel = dict()
lt = np.log(t)
lA = np.ones((np.shape(t)[0], 2))
lA[:, 0] = lt
lM = np.log(MMSD)
wl[0] = np.linalg.lstsq(lA, lM[:, 0])[0]
wl[1] = np.linalg.lstsq(lA, lM[:, 1])[0]
wl[2] = np.linalg.lstsq(lA, lM[:, 2])[0]
linel[0] = np.exp(wl[0][0]*lt + wl[0][1])
linel[1] = np.exp(wl[1][0]*lt + wl[1][1])
linel[2] = np.exp(wl[2][0]*lt + wl[2][1])
# Creates figure
fig = plt.figure(figsize=(24, 18), dpi=80)
ax = fig.add_subplot(111)
# ax.set_title('Particle Trajectories', x=0.5, y=1.15)
ax.plot(t, MMSD[:, 0], linestyle='', linewidth=20, marker='o', ms=20, color='blue')
ax.plot(t, MMSD[:, 1], linestyle='', linewidth=20, marker='o', ms=20, color='red')
ax.plot(t, MMSD[:, 2], linestyle='', linewidth=20, marker='o', ms=20, color='green')
ax.errorbar(t, MMSD[:, 0], yerr=SE[:, 0], fmt='', linestyle='', capsize=7, capthick=5, elinewidth=5, color='blue')
ax.errorbar(t, MMSD[:, 1], yerr=SE[:, 1], fmt='', linestyle='', capsize=7, capthick=5, elinewidth=5, color='red')
ax.errorbar(t, MMSD[:, 2], yerr=SE[:, 2], fmt='', linestyle='', capsize=7, capthick=5, elinewidth=5, color='green')
ax.plot(t, line[0], linewidth=6, label='2D', color='blue')
ax.plot(t, line[1], linewidth=6, label='1D x', color='red')
ax.plot(t, line[2], linewidth=6, label='1D y', color='green')
#ax.plot(t, linel[0], linewidth=3, color='blue')
#ax.plot(t, linel[1], linewidth=3, color='red')
#ax.plot(t, linel[2], linewidth=3, color='green')
# A few adjustments to prettify the graph
for item in ([ax.xaxis.label, ax.yaxis.label] +
ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(70)
xmajor_ticks = np.arange(0, limit1, tick1)
ymajor_ticks = np.arange(0, limit2, tick2)
hfont = {'fontname':'Arial'}
ax.set_xticks(xmajor_ticks)
ax.set_yticks(ymajor_ticks)
ax.title.set_fontsize(70)
ax.set_xlabel('Time (s)', fontsize=125, **hfont)
ax.set_ylabel(r'MSD ($\mu$m$^2$)', fontsize=125, **hfont)
ax.tick_params(direction='out', pad=16)
ax.legend(loc=(0.05, 0.66), prop={'size': 70}, frameon=False)
plt.gca().xaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%.{}f'.format(dec1)))
plt.gca().yaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%.{}f'.format(dec2)))
#plt.yscale('log')
#plt.xscale('log')
plt.gca().set_xlim([0, limit1])
plt.gca().set_ylim([0, limit2])
# Save your figure
plt.savefig('{}.png'.format(filename), bbox_inches='tight')
return MMSD, total1
In [7]:
one1, total1= plot_Mean1DMSDsorDeff(nephi, 0, 15, 9, 0, 0, 'MSD (um^2)', Mplot, 10.1, 500.01, 2, 100)
plt.show()
In [8]:
def plot_Mean1DMSDsorDeffold(traj, n1, n2, n3, dec1, dec2, datatype, filename, limit1, limit2, tick1, tick2):
"""
Plots the MSDs or Deffs from a trajectory dataset.
n1: particle numbers (typically 0)
n2: time (typically 1)
n3: MSDs or Deffs (11 or 15 typically)
"""
# Creates an array 'particles' that contains the particle number at each frame.
particles = traj[:, n1]
total = int(max(particles))
total1 = total + 1
rawtime = traj[:, n2]
bow = traj.shape[0]
raw2DMSDs = np.zeros((bow, 4))
raw2DMSDs[:, 0] = traj[:, n3]
raw2DMSDs[:, 1:4] = traj[:, n3 + 3:n3 + 6]
MSD = dict()
time = dict()
# Creates an array for each trajectory containing all xyz data
for num in range(1, total1):
hold = np.where(particles == num)
itindex = hold[0]
min1 = min(itindex)
max1 = max(itindex)
MSD[num] = (raw2DMSDs[min1+2:max1, :])
time[num] = (rawtime[min1+2:max1])
MMSD = MSD[1]
for num in range(2, total1):
MMSD = MMSD + MSD[num]
MMSD = MMSD/total1
# Creates figure
fig = plt.figure(figsize=(24, 18), dpi=80)
ax = fig.add_subplot(111)
# ax.set_title('Particle Trajectories', x=0.5, y=1.15)
ax.plot(time[1][:], MMSD[:, 0], linewidth=10, label='2D')
ax.plot(time[1][:], MMSD[:, 1], linewidth=10, label='1D x')
ax.plot(time[1][:], MMSD[:, 2], linewidth=10, label='1D y')
# A few adjustments to prettify the graph
for item in ([ax.xaxis.label, ax.yaxis.label] +
ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(90)
xmajor_ticks = np.arange(0, limit1, tick1)
ymajor_ticks = np.arange(0, limit2, tick2)
hfont = {'fontname':'Arial'}
ax.set_xticks(xmajor_ticks)
ax.set_yticks(ymajor_ticks)
ax.title.set_fontsize(70)
ax.set_xlabel('Time (s)', fontsize=115, **hfont)
ax.set_ylabel(r'D ($\mu$m$^2$/s)', fontsize=115, **hfont)
ax.tick_params(direction='out', pad=16)
ax.legend(loc=(0.65, 0.20), prop={'size': 70}, frameon=False)
plt.gca().xaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%.{}f'.format(dec1)))
plt.gca().yaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%.{}f'.format(dec2)))
plt.gca().set_xlim([0, limit1])
plt.gca().set_ylim([0, limit2])
# Save your figure
plt.savefig('{}.png'.format(filename), bbox_inches='tight')
return MMSD
In [9]:
one1 = plot_Mean1DMSDsorDeffold(nephi, 0, 15, 17, 0, 0, 'Deff (um^2/s)', Dplot, 10.1, 100.1, 2, 20)
plt.show()
In [ ]: