Commit 1ff5e4cf authored by amichaut's avatar amichaut
Browse files

added resolution_factor + notebook cleaning + synthetic data fixing + voronoi line_width

parent e56e5a47
This source diff could not be displayed because it is too large. You can view the blob instead.
...@@ -847,7 +847,8 @@ def fit_msd(msd, mean_vel=None, traj=None, dim=2, model_bounds={'P': [0, 300], ' ...@@ -847,7 +847,8 @@ def fit_msd(msd, mean_vel=None, traj=None, dim=2, model_bounds={'P': [0, 300], '
raise Exception("Error: mean_vel is required to fit with {}. Aborting...".format(model)) raise Exception("Error: mean_vel is required to fit with {}. Aborting...".format(model))
if model == 'PRW': if model == 'PRW':
param_list = ['P'] param_list = ['P']
func = lambda t, P: 2 * (mean_vel ** 2) * P * (t - P * (1 - np.exp(-t / P))) # mistake in Liu et al 2015, use Stokes 1991 definition instead
func = lambda t, P: (mean_vel ** 2) * P * (t - P * (1 - np.exp(-t / P)))
func_model = Model(func) func_model = Model(func)
p = func_model.make_params(P=10) p = func_model.make_params(P=10)
p['P'].set(min=model_bounds['P'][0], max=model_bounds['P'][1]) p['P'].set(min=model_bounds['P'][0], max=model_bounds['P'][1])
......
...@@ -51,6 +51,7 @@ def make_plot_config(data_dir=None, export_config=False): ...@@ -51,6 +51,7 @@ def make_plot_config(data_dir=None, export_config=False):
# Plotting config # Plotting config
plot_config = {'figsize': (5, 5), plot_config = {'figsize': (5, 5),
'dpi': 300, 'dpi': 300,
'figsize_factor':1, # factor to modulate figsize on traj and map plots
'color_list': color_list, 'color_list': color_list,
'format': '.png', # plot image format 'format': '.png', # plot image format
'despine': True, # use seaborn despine function 'despine': True, # use seaborn despine function
...@@ -124,6 +125,7 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim' ...@@ -124,6 +125,7 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim'
# unpack config # unpack config
color_list = plot_config['color_list'] color_list = plot_config['color_list']
figsize_factor = plot_config['figsize_factor']
show_tail = traj_parameters['show_tail'] show_tail = traj_parameters['show_tail']
color_code = traj_parameters['color_code'] color_code = traj_parameters['color_code']
cmap = traj_parameters['cmap'] cmap = traj_parameters['cmap']
...@@ -164,7 +166,8 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim' ...@@ -164,7 +166,8 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim'
# else: # else:
# Get background image # Get background image
bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg, bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg,
image_size=image_size,axis_on=show_axis, dpi=plot_config['dpi']) image_size=image_size,axis_on=show_axis, dpi=plot_config['dpi'],
figsize_factor=figsize_factor)
fig = bkgd['fig'] fig = bkgd['fig']
ax = bkgd['ax'] ax = bkgd['ax']
xmin = bkgd['xmin'] xmin = bkgd['xmin']
...@@ -273,7 +276,7 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim' ...@@ -273,7 +276,7 @@ def plot_traj(df, frame, data_dir, groups=None, image={'image_fn': None, 't_dim'
return fig return fig
else: else:
filename = osp.join(plot_dir, '{:04d}{}'.format(int(frame), plot_config['format'])) filename = osp.join(plot_dir, '{:04d}{}'.format(int(frame), plot_config['format']))
fig.savefig(filename, dpi=plot_config['dpi']) fig.savefig(filename)
plt.close(fig) plt.close(fig)
...@@ -294,6 +297,7 @@ def plot_scalar_field(data_dir, df, data, field, frame, image={'image_fn': None, ...@@ -294,6 +297,7 @@ def plot_scalar_field(data_dir, df, data, field, frame, image={'image_fn': None,
plot_config = make_plot_config() if plot_config is None else plot_config plot_config = make_plot_config() if plot_config is None else plot_config
save_as_stack = plot_config['save_as_stack'] save_as_stack = plot_config['save_as_stack']
figsize_factor = plot_config['figsize_factor']
no_bkg = map_param['no_bkg'] no_bkg = map_param['no_bkg']
show_axis = map_param['show_axis'] show_axis = map_param['show_axis']
cmap = map_param['cmap'] cmap = map_param['cmap']
...@@ -320,7 +324,8 @@ def plot_scalar_field(data_dir, df, data, field, frame, image={'image_fn': None, ...@@ -320,7 +324,8 @@ def plot_scalar_field(data_dir, df, data, field, frame, image={'image_fn': None,
no_bkg = True no_bkg = True
bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg, bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg,
image_size=image_size, axis_on=show_axis,dpi=plot_config['dpi']) image_size=image_size, axis_on=show_axis,dpi=plot_config['dpi'],
figsize_factor=figsize_factor)
fig = bkgd['fig'] fig = bkgd['fig']
ax = bkgd['ax'] ax = bkgd['ax']
xmin = bkgd['xmin'] xmin = bkgd['xmin']
...@@ -349,7 +354,7 @@ def plot_scalar_field(data_dir, df, data, field, frame, image={'image_fn': None, ...@@ -349,7 +354,7 @@ def plot_scalar_field(data_dir, df, data, field, frame, image={'image_fn': None,
return fig return fig
else: else:
filename = osp.join(plot_dir, field + '_{:04d}.png'.format(int(frame))) filename = osp.join(plot_dir, field + '_{:04d}.png'.format(int(frame)))
fig.savefig(filename, dpi=plot_config['dpi']) fig.savefig(filename)
plt.close(fig) plt.close(fig)
...@@ -371,6 +376,7 @@ def plot_vector_field(data_dir, df, data, field, frame, plot_on_field=None, dim= ...@@ -371,6 +376,7 @@ def plot_vector_field(data_dir, df, data, field, frame, plot_on_field=None, dim=
# misc param # misc param
plot_config = make_plot_config() if plot_config is None else plot_config plot_config = make_plot_config() if plot_config is None else plot_config
save_as_stack = plot_config['save_as_stack'] save_as_stack = plot_config['save_as_stack']
figsize_factor = plot_config['figsize_factor']
no_bkg = map_param['no_bkg'] no_bkg = map_param['no_bkg']
show_axis = map_param['show_axis'] show_axis = map_param['show_axis']
# cmap=map_param['cmap'] # cmap=map_param['cmap']
...@@ -408,8 +414,8 @@ def plot_vector_field(data_dir, df, data, field, frame, plot_on_field=None, dim= ...@@ -408,8 +414,8 @@ def plot_vector_field(data_dir, df, data, field, frame, plot_on_field=None, dim=
if no_plot_on_field: if no_plot_on_field:
bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg, bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg,
image_size=image_size, axis_on=show_axis, image_size=image_size, axis_on=show_axis,
dpi=plot_config['dpi']) dpi=plot_config['dpi'],figsize_factor=figsize_factor)
fig = bkgd['fig'] fig = bkgd['fig']
ax = bkgd['ax'] ax = bkgd['ax']
xmin = bkgd['xmin'] xmin = bkgd['xmin']
...@@ -440,13 +446,13 @@ def plot_vector_field(data_dir, df, data, field, frame, plot_on_field=None, dim= ...@@ -440,13 +446,13 @@ def plot_vector_field(data_dir, df, data, field, frame, plot_on_field=None, dim=
return fig return fig
else: else:
filename = osp.join(plot_dir, 'vector_' + field + '_{:04d}.png'.format(int(frame))) filename = osp.join(plot_dir, 'vector_' + field + '_{:04d}.png'.format(int(frame)))
fig.savefig(filename, dpi=plot_config['dpi']) fig.savefig(filename)
plt.close(fig) plt.close(fig)
def plot_Voronoi(data_dir, df, frame, data, show_local_area=True, def plot_Voronoi(data_dir, df, frame, data, show_local_area=True,
image={'image_fn': None, 't_dim': None, 'z_dim': None}, image={'image_fn': None, 't_dim': None, 'z_dim': None},
map_param={'no_bkg': False, 'vlim': None, 'show_axis': False, 'cmap': 'plasma', 'size_factor': 1}, map_param={'no_bkg': False, 'vlim': None, 'show_axis': False,'cmap':'plasma','size_factor': 1,'line_width':1.},
plot_dir=None, plot_config=None, dont_print_count=False): plot_dir=None, plot_config=None, dont_print_count=False):
""" """
Plot Voronoi tesselation and local area in 2D only. Plot Voronoi tesselation and local area in 2D only.
...@@ -474,10 +480,12 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True, ...@@ -474,10 +480,12 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True,
# misc param # misc param
plot_config = make_plot_config() if plot_config is None else plot_config plot_config = make_plot_config() if plot_config is None else plot_config
save_as_stack = plot_config['save_as_stack'] save_as_stack = plot_config['save_as_stack']
figsize_factor = plot_config['figsize_factor']
no_bkg = map_param['no_bkg'] no_bkg = map_param['no_bkg']
show_axis = map_param['show_axis'] show_axis = map_param['show_axis']
cmap = map_param['cmap'] cmap = map_param['cmap']
vlim = map_param['vlim'] vlim = map_param['vlim']
line_width = map_param['line_width']
info = tpr.get_info(data_dir) info = tpr.get_info(data_dir)
if 'image_width' not in info.keys() or 'image_height' not in info.keys(): if 'image_width' not in info.keys() or 'image_height' not in info.keys():
image_size = None image_size = None
...@@ -490,7 +498,8 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True, ...@@ -490,7 +498,8 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True,
no_bkg = True no_bkg = True
bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg, bkgd = tpr.get_background(image=image, frame=frame, df=df, no_bkg=no_bkg,
image_size=image_size, axis_on=show_axis, dpi=plot_config['dpi']) image_size=image_size, axis_on=show_axis, dpi=plot_config['dpi'],
figsize_factor=figsize_factor)
fig = bkgd['fig'] fig = bkgd['fig']
ax = bkgd['ax'] ax = bkgd['ax']
xmin = bkgd['xmin'] xmin = bkgd['xmin']
...@@ -502,7 +511,7 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True, ...@@ -502,7 +511,7 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True,
# plot tesselation # plot tesselation
vor = data[frame]['vor'] vor = data[frame]['vor']
if vor is not None: if vor is not None:
voronoi_plot_2d(vor, show_points=False, show_vertices=False, ax=ax) voronoi_plot_2d(vor, show_points=False, show_vertices=False, ax=ax, line_width=line_width)
# plot local area on top # plot local area on top
if show_local_area: if show_local_area:
...@@ -513,7 +522,7 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True, ...@@ -513,7 +522,7 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True,
area = areas[pt_id] area = areas[pt_id]
if not np.isnan(area): if not np.isnan(area):
color = tpr.get_cmap_color(area, cmap, vmin=vlim[0], vmax=vlim[1]) color = tpr.get_cmap_color(area, cmap, vmin=vlim[0], vmax=vlim[1])
ax.fill(*zip(*vor.vertices[indices]), color=color, alpha=0.5) ax.fill(*zip(*vor.vertices[indices]), color=color, alpha=0.5, linewidth=0)
# ensure axis limits are constant # ensure axis limits are constant
ax.set_xlim(xmin,xmax) ax.set_xlim(xmin,xmax)
...@@ -532,7 +541,7 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True, ...@@ -532,7 +541,7 @@ def plot_Voronoi(data_dir, df, frame, data, show_local_area=True,
return fig return fig
else: else:
filename = osp.join(plot_dir, 'voronoi_{:04d}.png'.format(int(frame))) filename = osp.join(plot_dir, 'voronoi_{:04d}.png'.format(int(frame)))
fig.savefig(filename, dpi=plot_config['dpi']) fig.savefig(filename)
plt.close(fig) plt.close(fig)
...@@ -1390,14 +1399,14 @@ def plot_total_traj(data_dir, df, dim=3, plot_dir=None, plot_fn=None, plot_confi ...@@ -1390,14 +1399,14 @@ def plot_total_traj(data_dir, df, dim=3, plot_dir=None, plot_fn=None, plot_confi
cmap_lim = [df['z_scaled'].min(), df['z_scaled'].max()] cmap_lim = [df['z_scaled'].min(), df['z_scaled'].max()]
if len(cmap_lim) == 2: if len(cmap_lim) == 2:
plot_cmap(plot_dir, tpr.make_param_label('z', l_unit=info["length_unit"]), cmap, cmap_lim[0], cmap_lim[1]) plot_cmap(plot_dir, tpr.make_param_label('z', l_unit=info["length_unit"]), cmap, cmap_lim[0], cmap_lim[1],suffix='_total_traj')
elif color_code == 't': elif color_code == 't':
if cmap_lim is None: if cmap_lim is None:
cmap_lim = [df['t'].min(), df['t'].max()] cmap_lim = [df['t'].min(), df['t'].max()]
if len(cmap_lim) == 2: if len(cmap_lim) == 2:
plot_cmap(plot_dir, tpr.make_param_label('t', l_unit=info["time_unit"]), cmap, cmap_lim[0], cmap_lim[1]) plot_cmap(plot_dir, tpr.make_param_label('t', l_unit=info["time_unit"]), cmap, cmap_lim[0], cmap_lim[1],suffix='_total_traj')
elif color_code == 'group': elif color_code == 'group':
if 'subset' in df.columns: if 'subset' in df.columns:
......
...@@ -1005,7 +1005,7 @@ def select_sub_data(df, filters=[]): ...@@ -1005,7 +1005,7 @@ def select_sub_data(df, filters=[]):
def get_background(image=None, frame=None, df=None, no_bkg=False, image_size=None, orig=None, axis_on=False, def get_background(image=None, frame=None, df=None, no_bkg=False, image_size=None, orig=None, axis_on=False,
dpi=plot_param['dpi'], figsize=(5, 5)): dpi=plot_param['dpi'], figsize=(5, 5), figsize_factor=1):
"""Get image background or create white background if no_bkg. The image can be a time nd stack or a single image.""" """Get image background or create white background if no_bkg. The image can be a time nd stack or a single image."""
if orig is None: if orig is None:
# orig = 'lower' if image_dir is None else 'upper' #trick to plot for the first time only inverting Yaxis: not very elegant... # orig = 'lower' if image_dir is None else 'upper' #trick to plot for the first time only inverting Yaxis: not very elegant...
...@@ -1032,7 +1032,7 @@ def get_background(image=None, frame=None, df=None, no_bkg=False, image_size=Non ...@@ -1032,7 +1032,7 @@ def get_background(image=None, frame=None, df=None, no_bkg=False, image_size=Non
ymax = df['y'].max() + 0.05 * (df['y'].max() - df['y'].min()) ymax = df['y'].max() + 0.05 * (df['y'].max() - df['y'].min())
figsize = ((xmax - xmin) / dpi, (ymax - ymin) / dpi) figsize = ((xmax - xmin) / dpi, (ymax - ymin) / dpi)
fig = plt.figure(frameon=False, figsize=figsize, dpi=dpi) fig = plt.figure(frameon=False, figsize=figsize, dpi=dpi*figsize_factor)
ax = fig.add_axes([0, 0, 1, 1]) ax = fig.add_axes([0, 0, 1, 1])
ax.set_aspect('equal') ax.set_aspect('equal')
...@@ -1054,7 +1054,7 @@ def get_background(image=None, frame=None, df=None, no_bkg=False, image_size=Non ...@@ -1054,7 +1054,7 @@ def get_background(image=None, frame=None, df=None, no_bkg=False, image_size=Non
im = img_as_ubyte(im) # 8bit conversion im = img_as_ubyte(im) # 8bit conversion
n = im.shape[0] n = im.shape[0]
m = im.shape[1] m = im.shape[1]
fig = plt.figure(frameon=False, figsize=(m / dpi, n / dpi), dpi=dpi) fig = plt.figure(frameon=False, figsize=(m / dpi, n / dpi), dpi=dpi*figsize_factor)
ax = fig.add_axes([0, 0, 1, 1]) ax = fig.add_axes([0, 0, 1, 1])
ax.imshow(im, aspect='equal', origin=orig, cmap='gray', vmin=0, vmax=255) ax.imshow(im, aspect='equal', origin=orig, cmap='gray', vmin=0, vmax=255)
xmin, xmax, ymin, ymax = ax.axis() xmin, xmax, ymin, ymax = ax.axis()
...@@ -1306,7 +1306,8 @@ def make_traj_config(data_dir=None, export_config=True): ...@@ -1306,7 +1306,8 @@ def make_traj_config(data_dir=None, export_config=True):
'area_threshold': 3, # exclude areas above this multiple of area median 'area_threshold': 3, # exclude areas above this multiple of area median
'no_bkg': False, # don't show background image if an image is passed 'no_bkg': False, # don't show background image if an image is passed
'size_factor': 1., # to multiply the default size of lines !! Not implemented yet !! 'size_factor': 1., # to multiply the default size of lines !! Not implemented yet !!
'show_axis': False, # to show the plot axes (by default just image) 'show_axis': False, # to show the plot axes (by default just image)
'line_width': 1., # diagram line width
} }
# package all in a dict # package all in a dict
......
...@@ -75,8 +75,14 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa ...@@ -75,8 +75,14 @@ def traj_analysis(data_dir, data=None, image=None, refresh=False, parallelize=Fa
# check that all configs are in traj_confign, if not load default # check that all configs are in traj_confign, if not load default
for key in ["traj_config_", "MSD_config", "scatter_config", "hist_config", "total_traj_config", "voronoi_config"]: for key in ["traj_config_", "MSD_config", "scatter_config", "hist_config", "total_traj_config", "voronoi_config"]:
traj_config_default_key = traj_config_default[key]
# if don't exist, add
if key not in traj_config.keys(): if key not in traj_config.keys():
traj_config[key] = traj_config_default[key] traj_config[key] = traj_config_default_key
# check all defaut param are present
for k in traj_config_default_key.keys():
if k not in traj_config[key].keys():
traj_config[key][k] = traj_config_default_key[k]
traj_config_ = traj_config["traj_config_"] traj_config_ = traj_config["traj_config_"]
MSD_config = traj_config["MSD_config"] MSD_config = traj_config["MSD_config"]
......
...@@ -23,180 +23,183 @@ ...@@ -23,180 +23,183 @@
# If not, see <https://www.gnu.org/licenses/>. # # If not, see <https://www.gnu.org/licenses/>. #
########################################################################## ##########################################################################
import sys
import os.path as osp import os.path as osp
import datetime
import csv
import shutil import shutil
import itertools
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import pandas as pd import pandas as pd
from skimage import io from skimage import io
from skimage.color import rgb2gray
from skimage.util import img_as_ubyte from skimage.util import img_as_ubyte
import scipy.interpolate as sci
import pickle
import seaborn as sns import seaborn as sns
from lmfit import Parameters, Model
from scipy.optimize import curve_fit
from scipy import stats
import matplotlib.tri as tri
import napari
import tifffile as tifff import tifffile as tifff
from track_analyzer import prepare as tpr
# Plotting parameters # Plotting parameters
color_list=[c['color'] for c in list(plt.rcParams['axes.prop_cycle'])]+sns.color_palette("Set1", n_colors=9, desat=.5) color_list = [c['color'] for c in list(plt.rcParams['axes.prop_cycle'])] + sns.color_palette("Set1", n_colors=9,
plot_param={'figsize':(5,5),'dpi':300,'color_list':color_list,'format':'.png','despine':True,'logx':False,'logy':False,'invert_yaxis':True,'export_data_pts':False} desat=.5)
plot_param = {'figsize': (5, 5), 'dpi': 300, 'color_list': color_list, 'format': '.png', 'despine': True, 'logx': False,
'logy': False, 'invert_yaxis': True, 'export_data_pts': False}
def make_diff_traj(part_index=0, grid_size=[500, 500, 500], dim=3, tmax=10, periodic=True, noise_amp=10,
def make_diff_traj(part_index=0,grid_size=[500,500,500],dim=3,tmax=10,periodic=True,noise_amp=10,x0=[250,250,250],bias=[0,0,0]): x0=[250, 250, 250], bias=[0, 0, 0]):
"""Generate a trajectory with a diffusive trajectory, with a bias. bias gives the amplitude at each step along each dimension.""" """
#time and index Generate a trajectory with a diffusive trajectory, with a bias.
t = arange(tmax) bias gives the amplitude at each step along each dimension.
index = ones(tmax)*part_index """
#displacement # time and index
displacement=pd.DataFrame(np.random.randn(tmax,dim),columns=list('xyz')[0:dim]) t = np.arange(tmax)
displacement['r2']=0 index = np.ones(tmax) * part_index
# displacement
displacement = pd.DataFrame(np.random.randn(tmax, dim), columns=list('xyz')[0:dim])
displacement['r2'] = 0
for i in range(dim): for i in range(dim):
displacement['r2']+=displacement[list('xyz')[i]]**2 displacement['r2'] += displacement[list('xyz')[i]] ** 2
displacement['r']=np.sqrt(displacement['r2']) displacement['r'] = np.sqrt(displacement['r2'])
for i in range(dim): for i in range(dim):
displacement[list('xyz')[i]]/=displacement['r'] #normalize raw displacement displacement[list('xyz')[i]] /= displacement['r'] # normalize raw displacement
displacement[list('xyz')[i]]*=noise_amp #amply amplitude displacement[list('xyz')[i]] *= noise_amp # amply amplitude
displacement[list('xyz')[i]]+=bias[i] #add bias displacement[list('xyz')[i]] += bias[i] # add bias
displacement=displacement[list('xyz')[0:dim]].values displacement = displacement[list('xyz')[0:dim]].values
#traj # traj
traj=np.zeros((tmax,dim)) traj = np.zeros((tmax, dim))
for i in range(dim): for i in range(dim):
traj[:,i]=np.cumsum(displacement[:,i])+x0[i] traj[:, i] = np.cumsum(displacement[:, i]) + x0[i]
if periodic: if periodic:
traj[:,i]=np.remainder(traj[:,i],grid_size[i]) traj[:, i] = np.remainder(traj[:, i], grid_size[i])
return pd.DataFrame(np.concatenate([index[:,None],t[:,None],traj],axis=1),columns=['traj','frame']+list('xyz')[0:dim]) return pd.DataFrame(np.concatenate([index[:, None], t[:, None], traj], axis=1),
columns=['traj', 'frame'] + list('xyz')[0:dim])
def make_spatial_gradient(part_num=100,grid_size=[500,500,500],dim=3,tmax=10,periodic=True,noise_amp=10,bias_basis=[0,0,0],
diff_grad={'min':0,'max':10},bias_grad={'min':0,'max':10,'dim':0},grad={'step_num':4,'dim':0}, def make_spatial_gradient(part_num=100, grid_size=[500, 500, 500], dim=3, tmax=10, periodic=True,
x0_range={'x':[0.1,0.9],'y':[0.1,0.9],'z':[0.1,0.9]},dt=1): bias_basis=[0, 0, 0],
diff_grad={'min': 0, 'max': 10}, bias_grad={'min': 0, 'max': 10, 'dim': 0},
grad={'step_num': 4, 'dim': 0},
x0_range={'x': [0.1, 0.9], 'y': [0.1, 0.9], 'z': [0.1, 0.9]}, dt=1):
"""Make a spatial gradient (number of steps on the gradient given by grad['step_num'})in diffusion or bias, along a specific dimension, given by grad['dim']. """Make a spatial gradient (number of steps on the gradient given by grad['step_num'})in diffusion or bias, along a specific dimension, given by grad['dim'].
The gradient can be in diffusion with diff_grad or bias_grad. min and max give the extrema of the gradient, and bias_grad['dim'] give the dimension along the gradient in bias is applied. The gradient can be in diffusion with diff_grad or bias_grad. min and max give the extrema of the gradient, and bias_grad['dim'] give the dimension along the gradient in bias is applied.
An overall constant bias can be passed by bias_basis. An overall constant bias can be passed by bias_basis.
""" """
df=pd.DataFrame([],columns=['traj','frame']+list('xyz')[0:dim]) df = pd.DataFrame([], columns=['traj', 'frame'] + list('xyz')[0:dim])
df_param=pd.DataFrame([],columns=['traj','v','D']) df_param = pd.DataFrame([], columns=['traj', 'v', 'D'])
diff_grad_=np.linspace(diff_grad['min'],diff_grad['max'],grad['step_num']) diff_grad_ = np.linspace(diff_grad['min'], diff_grad['max'], grad['step_num'])
bias_grad_=np.linspace(bias_grad['min'],bias_grad['max'],grad['step_num']) bias_grad_ = np.linspace(bias_grad['min'], bias_grad['max'], grad['step_num'])
#spatial boundaries of the regions of particles # spatial boundaries of the regions of particles
lims=[[x0_range['x'][0]*grid_size[0],x0_range['x'][1]*grid_size[0]], lims = [[x0_range['x'][0] * grid_size[0], x0_range['x'][1] * grid_size[0]],
[x0_range['y'][0]*grid_size[1],x0_range['y'][1]*grid_size[1]], [x0_range['y'][0] * grid_size[1], x0_range['y'][1] * grid_size[1]],
[x0_range['z'][0]*grid_size[2],x0_range['z'][1]*grid_size[2]]] [x0_range['z'][0] * grid_size[2], x0_range['z'][1] * grid_size[2]]]
part_count=0 part_count = 0
for i in range(grad['step_num']): for i in range(grad['step_num']):
grad_increment=(lims[grad['dim']][1]-lims[grad['dim']][0])/grad['step_num'] grad_increment = (lims[grad['dim']][1] - lims[grad['dim']][0]) / grad['step_num']
lims_=lims[:] lims_ = lims[:]
lims_[grad['dim']]=[lims_[grad['dim']][0]+i*grad_increment,lims_[grad['dim']][0]+(i+1)*grad_increment] lims_[grad['dim']] = [lims_[grad['dim']][0] + i * grad_increment,
noise_amp=diff_grad_[i] lims_[grad['dim']][0] + (i + 1) * grad_increment]
bias=bias_basis[:] noise_amp = diff_grad_[i]
bias[bias_grad['dim']]=bias_grad_[i] bias = bias_basis[:]
bias_ampl=0 bias[bias_grad['dim']] = bias_grad_[i]
bias_ampl = 0
for k in range(dim): for k in range(dim):
bias_ampl+=bias[k]**2 bias_ampl += bias[k] ** 2
bias_ampl=np.sqrt(bias_ampl) bias_ampl = np.sqrt(bias_ampl)
for j in range(int(part_num/grad['step_num'])): for j in range(int(part_num / grad['step_num'])):
x0=[np.random.uniform(lims_[0][0],lims_[0][1]), x0 = [np.random.uniform(lims_[0][0], lims_[0][1]),
np.random.uniform(lims_[1][0],lims_[1][1]), np.random.uniform(lims_[1][0], lims_[1][1]),
np.random.uniform(lims_[2][0],lims_[2][1])] np.random.uniform(lims_[2][0], lims_[2][1])]
traj=make_diff_traj(part_index=part_count,noise_amp=noise_amp,x0=x0,bias=bias,tmax=tmax,periodic=periodic,dim=dim) traj = make_diff_traj(part_index=part_count, noise_amp=noise_amp, x0=x0, bias=bias, tmax=tmax,
df=pd.concat([df,traj]) periodic=periodic, dim=dim)
v=bias_ampl/dt df = pd.concat([df, traj])
D=noise_amp**2/(2.*dim*dt) v = bias_ampl / dt
df_param.loc[part_count,:]=[part_count,v,D] D = noise_amp ** 2 / (2. * dim * dt)
df_param.loc[part_count, :] = [part_count, v, D]
part_count+=1
return df,df_param part_count += 1
return df, df_param
def make_attraction_node(part_num=100,grid_size=[500,500,500],dim=3,tmax=10,periodic=True,noise_amp=10,bias_basis=[0,0,0],
attraction_ampl=10,node=None,x0_range={'x':[0.1,0.9],'y':[0.1,0.9],'z':[0.1,0.9]},dt=1):
def make_attraction_node(part_num=100, grid_size=[500, 500, 500], dim=3, tmax=10, periodic=True, noise_amp=10,
bias_basis=[0, 0, 0],
attraction_ampl=10, node=None, x0_range={'x': [0.1, 0.9], 'y': [0.1, 0.9], 'z': [0.1, 0.9]},
dt=1):
"""Make array of diffusive particles biased toward a node (or away if attraction_ampl is negative)""" """Make array of diffusive particles biased toward a node (or away if attraction_ampl is negative)"""
df=pd.DataFrame([],columns=['traj','frame']+list('xyz')[0:dim]) df = pd.DataFrame([], columns=['traj', 'frame'] + list('xyz')[0:dim])
df_param=pd.DataFrame([],columns=['traj','v','D']) df_param = pd.DataFrame([], columns=['traj', 'v', 'D'])
if node is None: if node is None:
node=[grid_size[d]/2 for d in range(dim)] #by default center node = [grid_size[d] / 2 for d in range(dim)] # by default center
# spatial boundaries of the regions of particles
lims = [[x0_range['x'][0] * grid_size[0], x0_range['x'][1] * grid_size[0]],
[x0_range['y'][0] * grid_size[1], x0_range['y'][1] * grid_size[1]],
[x0_range['z'][0] * grid_size[2], x0_range['z'][1] * grid_size[2]]]
#spatial boundaries of the regions of particles
lims=[[x0_range['x'][0]*grid_size[0],x0_range['x'][1]*grid_size[0]],
[x0_range['y'][0]*grid_size[1],x0_range['y'][1]*grid_size[1]],
[x0_range['z'][0]*grid_size[2],x0_range['z'][1]*grid_size[2]]]
for i in range(part_num): for i in range(part_num):
x0=[np.random.uniform(lims[0][0],lims[0][1]), x0 = [np.random.uniform(lims[0][0], lims[0][1]),
np.random.uniform(lims[1][0],lims[1][1]), np.random.uniform(lims[1][0], lims[1][1]),