In [1]:
# !/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 20181217

@author: zhangji
"""

%pylab inline
pylab.rcParams['figure.figsize'] = (25, 11)
fontsize = 40

from time import time
import os
import glob
import numpy as np
import scipy as sp
import scipy.misc
import pandas as pd
import re
from scanf import scanf
from matplotlib import pyplot as plt
import matplotlib.ticker as mtick
import matplotlib
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
from mpl_toolkits.mplot3d import Axes3D, axes3d
from mpl_toolkits.mplot3d.art3d import Line3DCollection
from scipy.optimize import leastsq, curve_fit
from scipy.interpolate import interp1d
from scipy.io import loadmat, savemat
from IPython.display import display, HTML
from scipy import interpolate
from matplotlib import animation, rc
from IPython.display import HTML
rc('animation', html='html5')
from tqdm import tqdm_notebook as tqdm
import natsort 
from shutil import copyfile

from src import jeffery_model as jm
from codeStore import support_fun as spf
from src.support_class import *

PWD = os.getcwd()
font = {'size': 20}
matplotlib.rc('font', **font)
np.set_printoptions(linewidth=90, precision=5)


Populating the interactive namespace from numpy and matplotlib

In [2]:
def set_axes_equal(ax):
    if  ax.name == "3d":
        '''Make axes of 3D plot have equal scale so that spheres appear as spheres,
        cubes as cubes, etc..  This is one possible solution to Matplotlib's
        ax.set_aspect('equal') and ax.axis('equal') not working for 3D.

        Input
          ax: a matplotlib axis, e.g., as output from plt.gca().
        '''

        limits = np.array([
            ax.get_xlim3d(),
            ax.get_ylim3d(),
            ax.get_zlim3d(),
        ])

        origin = np.mean(limits, axis=1)
        radius = 0.6 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
        ax.set_xlim3d([origin[0] - radius, origin[0] + radius])
        ax.set_ylim3d([origin[1] - radius, origin[1] + radius])
        ax.set_zlim3d([origin[2] - radius, origin[2] + radius])
    else:
        limits = np.array([
            ax.get_xlim(),
            ax.get_ylim(),
        ])

        origin = np.mean(limits, axis=1)
        radius = 0.6 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
        ax.set_xlim([origin[0] - radius, origin[0] + radius])
        ax.set_ylim([origin[1] - radius, origin[1] + radius])

In [3]:
def read_helix_data(dir_name, file_handle, readLength=None):
    t_dir = os.path.join(PWD, dir_name, file_handle)
    mat_names = natsort.natsorted(glob.glob('%s/%s_*.mat' % (t_dir, file_handle)))
    if readLength is None:
        readLength = len(mat_names)
    ecoli_U_list = []
    ecoli_norm_list = []
    ecoli_center_list = []
    ecoli_nodes_list = []
    ecoli_u_list = []
    ecoli_f_list = []
    ecoli_t_list = []
    for mati in mat_names[:readLength]:
        mat_contents = loadmat(mati)
        ecoli_U = mat_contents['ecoli_U'].flatten()
        ecoli_norm = mat_contents['ecoli_norm'].flatten()
        ecoli_center = mat_contents['ecoli_center'].flatten()
        ecoli_t = mat_contents['ti'].flatten()
        ecoli_nodes = mat_contents['ecoli_nodes']
        ecoli_u = mat_contents['ecoli_u']
        ecoli_f = mat_contents['ecoli_f']
        ecoli_U_list.append(ecoli_U)
        ecoli_norm_list.append(ecoli_norm)
        ecoli_center_list.append(ecoli_center)
        ecoli_nodes_list.append(ecoli_nodes)
        ecoli_t_list.append(ecoli_t)
    #     ecoli_u_list.append(ecoli_u)
    #     ecoli_f_list.append(ecoli_f)
    #     set_axes_equal(ax)

    ecoli_U = np.vstack(ecoli_U_list)
    ecoli_norm = np.vstack(ecoli_norm_list)
    ecoli_center = np.vstack(ecoli_center_list)
    ecoli_t = np.hstack(ecoli_t_list)
    print(dir_name, file_handle, ecoli_t.size)
    return ecoli_U, ecoli_norm, ecoli_center, ecoli_t, ecoli_nodes_list

In [32]:
# dir_name = 'motion_helix_passive1'
# file_handle = 'hlx_tau01.00'
# anim_step = 1

# dir_name = 'motion_ecoli_torque6'
# file_handle = 'passive_tau01.00'
# anim_step = 1

# dir_name = 'motion_helix_passive2b'
# file_handle = 'hlx_ntail1'

dir_name = 'motion_helix_passive1'
file_handle = 'hlx_prl_ch15'

# dir_name = 'motion_helix_passive2c'
# file_handle = 'hlx_ntail1_theta0.15_phi4.34_t0.1'

anim_step = 1
tmp = scanf('hlx_ntail%d', file_handle)
n_tail = 1 if tmp is None else tmp[0]
ecoli_U, ecoli_norm, ecoli_center, ecoli_t, ecoli_nodes_list = read_helix_data(dir_name, file_handle)
ecoli_theta = np.arccos(ecoli_norm[:, 2] / np.linalg.norm(ecoli_norm, axis=1))
t_phi = np.arctan2(ecoli_norm[:, 1], ecoli_norm[:, 0])
ecoli_phi = np.hstack([t1 + 2 * np.pi if t1 < 0 else t1 for t1 in t_phi])

fig = plt.figure(figsize=(20, 15))
fig.patch.set_facecolor('white')
ax = axes3d.Axes3D(fig)
for spine in ax.spines.values():
    spine.set_visible(False)

Writer = animation.writers['ffmpeg']
writer = Writer(fps=25, bitrate=1800, metadata=dict(artist='ZhangJi', title=file_handle, copyright='MIT'))
t0 = np.split(ecoli_nodes_list[0], 2 * n_tail)
t1 = np.vstack(t0[0::2])
t2 = np.vstack(t0[1::2])
tl1 = ax.plot(t1[:, 0], t1[:, 1], t1[:, 2], color='b')[0]
tl2 = ax.plot(t2[:, 0], t2[:, 1], t2[:, 2], color='r')[0]
ax.set_xlim3d([0, 4])
ax.set_xlabel('X')
ax.set_ylim3d([-2, 2])
ax.set_ylabel('Y')
ax.set_zlim3d([-2, 2])
ax.set_zlabel('Z')
ax.set_title(file_handle, size=fontsize)

def update_fun(num, tl1, tl2, ecoli_nodes_list):
    num = num * anim_step
    t0 = np.split(ecoli_nodes_list[num], 2 * n_tail)
    t1 = np.vstack(t0[0::2])
    t2 = np.vstack(t0[1::2])
    tl1.set_data(t1[:, 0], t1[:, 1])
    tl1.set_3d_properties(t1[:, 2])
    tl2.set_data(t2[:, 0], t2[:, 1])
    tl2.set_3d_properties(t2[:, 2])
    return tl1, tl2, 

# line_ani = animation.FuncAnimation(fig, update_fun, len(ecoli_nodes_list) // anim_step, interval=50, 
#                                    blit=False, fargs=(tl1, tl2, ecoli_nodes_list, ),)
# line_ani.save(os.path.join(PWD, dir_name, file_handle + '.mp4'), writer=writer)
# line_ani
pass


motion_helix_passive1 hlx_prl_ch15 2740

In [33]:
print(ecoli_theta[-1], ecoli_phi[-1])


1.56958829509 1.57486732341

In [34]:
# calculate psi and lateral norm. 
# def get_rot_matrix(norm=np.array([0, 0, 1]), theta=0):
#     norm = np.array(norm).reshape((3,))
#     theta = -1 * float(theta)
#     if np.linalg.norm(norm) > 0:
#         norm = norm / np.linalg.norm(norm)
#     a = norm[0]
#     b = norm[1]
#     c = norm[2]
#     rotation = np.array([
#         [a ** 2 + (1 - a ** 2) * np.cos(theta),
#          a * b * (1 - np.cos(theta)) + c * np.sin(theta),
#          a * c * (1 - np.cos(theta)) - b * np.sin(theta)],
#         [a * b * (1 - np.cos(theta)) - c * np.sin(theta),
#          b ** 2 + (1 - b ** 2) * np.cos(theta),
#          b * c * (1 - np.cos(theta)) + a * np.sin(theta)],
#         [a * c * (1 - np.cos(theta)) + b * np.sin(theta),
#          b * c * (1 - np.cos(theta)) - a * np.sin(theta),
#          c ** 2 + (1 - c ** 2) * np.cos(theta)]])
#     return rotation

# def vector_rotation(P2, norm=np.array([0, 0, 1]), theta=0, rotation_origin=np.zeros(3)):
#     rotation = get_rot_matrix(norm, theta)
#     P20 = np.dot(rotation, (P2 - rotation_origin)) + rotation_origin
#     P20 = P20 / np.linalg.norm(P20)
#     return P20

ecoli_lateral_norm = []
ecoli_psi = []
idx = np.random.randint(0, ecoli_nodes_list[0].shape[0], 1)[0]
for t_nodes, t_center, t_norm in zip(ecoli_nodes_list, ecoli_center, ecoli_norm):
    r0 = t_nodes[idx] - t_center
    n0 = np.dot(r0, t_norm) * t_norm / np.dot(t_norm, t_norm)
    t0 = r0 - n0
    ecoli_lateral_norm.append(t0 / np.linalg.norm(t0))
ecoli_lateral_norm = np.vstack(ecoli_lateral_norm)
t_lateral_norm0 = ecoli_lateral_norm[0]
t_norm0 = ecoli_norm[0]
# t_lateral_norm0 = vector_rotation(t_lateral_norm0, norm=np.array((0, 0, 1)), theta=-ecoli_phi[0])
# t_lateral_norm0 = vector_rotation(t_lateral_norm0, norm=np.array((0, 1, 0)), theta=-ecoli_theta[0])
# t_norm0 = np.array((0, 0, 1))
for t_lateral_norm, t_theta, t_phi in zip(ecoli_lateral_norm, ecoli_theta, ecoli_phi):
    t_lateral_norm = vector_rotation(t_lateral_norm, norm=np.array((0, 0, 1)), theta=-t_phi)
    t_lateral_norm = vector_rotation(t_lateral_norm, norm=np.array((0, 1, 0)), theta=ecoli_theta[0]-t_theta)
    t_lateral_norm = vector_rotation(t_lateral_norm, norm=np.array((0, 0, 1)), theta=ecoli_phi[0])
    sign = np.sign(np.dot(t_norm0, np.cross(t_lateral_norm0, t_lateral_norm)))
    t_psi = sign * np.arccos(np.clip(np.dot(t_lateral_norm0, t_lateral_norm) 
                                     / np.linalg.norm(t_lateral_norm) / np.linalg.norm(t_lateral_norm0),
                                     -1, 1))
    t_psi = t_psi + 2 * np.pi if t_psi < 0 else t_psi  # (-pi,pi) -> (0, 2pi)
    ecoli_psi.append(t_psi)
ecoli_psi = np.hstack(ecoli_psi).flatten()
ecoli_eta = np.arccos(np.sin(ecoli_theta) * np.sin(ecoli_phi))

fig = plt.figure(figsize=(15, 10))
fig.patch.set_facecolor('white')
ax0, ax1, ax2, ax3 = fig.subplots(nrows=4, ncols=1)
ax0.plot(ecoli_t[1:], ecoli_theta[1:] / np.pi, '-*')
ax1.plot(ecoli_t[1:], ecoli_phi[1:] / np.pi, '-*')
ax2.plot(ecoli_t[1:], ecoli_psi[1:] / np.pi, '-*')
ax3.plot(ecoli_t[1:], ecoli_eta[1:] / np.pi, '-*')
ax0.set_ylabel('theta')
ax1.set_ylabel('phi')
ax2.set_ylabel('psi')
ax3.set_ylabel('eta')
# ax0.set_xlim(0,50)
# ax1.set_xlim(0,50)
# ax2.set_xlim(0,50)
pass
print(ecoli_norm[1])


[  2.93700e-01   5.01406e-05   9.55898e-01]

In [ ]:


In [35]:
# compare simulation and table results. 
ini_idx = 0
plot_length = ecoli_U.shape[0] - ini_idx
# ini_idx = np.random.randint(0, ecoli_t.shape[0]-plot_length, 1)[0]

# calculate Table result
importlib.reload(jm)
planeShearRate = np.array((1, 0, 0))
tcenter = ecoli_center[ini_idx]
tnorm = ecoli_norm[ini_idx]
ini_psi = ecoli_psi[ini_idx]
tlateral_norm = ecoli_lateral_norm[ini_idx]
eval_dt = np.mean(np.diff(ecoli_t))
max_iter = plot_length
update_order = 1
helix_kwargs = {'name':         'helix',
                'center':       tcenter,
                'norm':         tnorm / np.linalg.norm(tnorm),
                'lateral_norm': tlateral_norm / np.linalg.norm(tlateral_norm),
                'speed':        0,
                'lbd':          np.nan, 
                'ini_psi':      ini_psi, 
                'table_name':   'hlxB01_tau1a', }
fileHandle = 'ShearTableProblem'
helix_obj = jm.TableObj(**helix_kwargs)
helix_obj.set_update_para(fix_x=False, fix_y=False, fix_z=False, update_order=update_order)
problem = jm.ShearTableProblem(name=fileHandle, planeShearRate=planeShearRate)
problem.add_obj(helix_obj)
t0 = time()
for idx in range(1, max_iter + 1):
    problem.update_location(eval_dt, print_handle='%d / %d' % (idx, max_iter))
t1 = time()
Table_X = np.vstack(helix_obj.center_hist)
Table_U = np.vstack(helix_obj.U_hist)
Table_P = np.vstack(helix_obj.norm_hist)
Table_P2 = np.vstack(helix_obj.lateral_norm_hist)
Table_t = np.arange(max_iter) * eval_dt + eval_dt
Table_theta, Table_phi, Table_psi = helix_obj.theta_phi_psi
print('%s: run %d loops using %f' % (fileHandle, max_iter, (t1 - t0)))

theta_lim = [np.min(np.hstack((ecoli_theta[ini_idx:ini_idx+plot_length], Table_theta))) / np.pi, 
             np.max(np.hstack((ecoli_theta[ini_idx:ini_idx+plot_length], Table_theta))) / np.pi]
phi_lim = [np.min(np.hstack((ecoli_phi[ini_idx:ini_idx+plot_length], Table_phi))) / np.pi, 
             np.max(np.hstack((ecoli_phi[ini_idx:ini_idx+plot_length], Table_phi))) / np.pi]
psi_lim = [np.min(np.hstack((ecoli_psi[ini_idx:ini_idx+plot_length], Table_psi))) / np.pi, 
             np.max(np.hstack((ecoli_psi[ini_idx:ini_idx+plot_length], Table_psi))) / np.pi]
########################################################################################

# show simulation results. 
fig = plt.figure(figsize=(20, 8))
fig.patch.set_facecolor('white')
ax0 = plt.subplot2grid((3, 2), (0, 0), rowspan=3)
ax1 = plt.subplot2grid((3, 2), (0, 1), )
ax2 = plt.subplot2grid((3, 2), (1, 1), )
ax3 = plt.subplot2grid((3, 2), (2, 1), )
norm=plt.Normalize(ecoli_t[ini_idx:ini_idx+plot_length].min(), ecoli_t[ini_idx:ini_idx+plot_length].max())
cmap=plt.get_cmap('jet')
ax0.plot(ecoli_phi[ini_idx:ini_idx+plot_length] / np.pi, ecoli_theta[ini_idx:ini_idx+plot_length] / np.pi, ' ')
# ax0.plot(ecoli_phi / np.pi, ecoli_theta / np.pi, '*', ms=fontsize*0.5)
lc = spf.colorline(ecoli_phi[ini_idx:ini_idx+plot_length] / np.pi, 
                   ecoli_theta[ini_idx:ini_idx+plot_length] / np.pi, 
                   ecoli_t[ini_idx:ini_idx+plot_length], 
                   ax=ax0, cmap=cmap, norm=norm, linewidth=3)
clb = fig.colorbar(lc, ax=ax0, orientation="vertical")
ax0.set_xlabel('$\\phi / \pi$', size=fontsize)
ax0.set_ylabel('$\\theta / \pi$', size=fontsize)
plt.sca(ax0)
plt.xticks(fontsize=fontsize*0.5)
plt.yticks(fontsize=fontsize*0.5)
# ax1.plot(Jeffery_t, Jeffery_theta / np.pi, label='Jeffery')
# ax2.plot(Jeffery_t, Jeffery_phi / np.pi, label='Jeffery')
ax1.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_theta[ini_idx:ini_idx+plot_length] / np.pi, '-*', label='Sim')
ax2.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_phi[ini_idx:ini_idx+plot_length] / np.pi, '-*', label='Sim')
ax3.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_psi[ini_idx:ini_idx+plot_length] / np.pi, '-*', label='Sim')
for axi, axyi in zip((ax1, ax2, ax3), ('$\\theta / \pi$', '$\\phi / \pi$', '$\\psi / \pi$')):
    plt.sca(axi)
    axi.set_xlabel('t', size=fontsize)
    axi.set_ylabel('%s' % axyi, size=fontsize)
    axi.legend()
    plt.xticks(fontsize=fontsize*0.5)
    plt.yticks(fontsize=fontsize*0.5)
plt.tight_layout()
ax0.set_xlim(phi_lim)
ax0.set_ylim(theta_lim)
ax1.set_ylim(theta_lim)
ax2.set_ylim(phi_lim)
ax3.set_ylim(psi_lim)
print('ini_idx', ini_idx)
print('ecoli_norm[ini_idx]', ecoli_norm[ini_idx])

# velocity of simulation results
fig, [(ax0, ax3), (ax1, ax4), (ax2, ax5), (ax6, ax7)] = plt.subplots(nrows=4, ncols=2, figsize=(20, 20))
fig.patch.set_facecolor('white')
ax0.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_U[ini_idx:ini_idx+plot_length, 0], '-*')
ax1.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_U[ini_idx:ini_idx+plot_length, 1], '-*')
ax2.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_U[ini_idx:ini_idx+plot_length, 2], '-*')
ax3.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_U[ini_idx:ini_idx+plot_length, 3], '-*')
ax4.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_U[ini_idx:ini_idx+plot_length, 4], '-*')
ax5.plot(ecoli_t[ini_idx:ini_idx+plot_length], ecoli_U[ini_idx:ini_idx+plot_length, 5], '-*')
ax6.plot(ecoli_t[ini_idx:ini_idx+plot_length], np.linalg.norm(ecoli_U[ini_idx:ini_idx+plot_length, :3], axis=1), '-*')
ax7.plot(ecoli_t[ini_idx:ini_idx+plot_length], np.linalg.norm(ecoli_U[ini_idx:ini_idx+plot_length, 3:], axis=1), '-*')
for axi, axyi in zip((ax0, ax1, ax2, ax3, ax4, ax5, ax6, ax7),
                     ('ecoli_u_x', 'ecoli_u_y', 'ecoli_u_z', 
                      'ecoli_w_x', 'ecoli_w_y', 'ecoli_w_z', 
                      '|ecoli_U|', '|ecoli_W|')):
    plt.sca(axi)
    axi.set_xlabel('time', size=fontsize)
    axi.set_ylabel('%s' % axyi, size=fontsize)
    plt.xticks(fontsize=fontsize*0.8)
    plt.yticks(fontsize=fontsize*0.8)
plt.tight_layout()

###################################################################################
# show table results. 
fig = plt.figure(figsize=(20, 8))
fig.patch.set_facecolor('white')
ax0 = plt.subplot2grid((3, 2), (0, 0), rowspan=3)
ax1 = plt.subplot2grid((3, 2), (0, 1), )
ax2 = plt.subplot2grid((3, 2), (1, 1), )
ax3 = plt.subplot2grid((3, 2), (2, 1), )
norm=plt.Normalize(Table_t.min(), Table_t.max())
cmap=plt.get_cmap('jet')
ax0.plot(Table_phi / np.pi, Table_theta / np.pi, ' ')
lc = spf.colorline(Table_phi / np.pi, Table_theta / np.pi, Table_t, 
                   ax=ax0, cmap=cmap, norm=norm, linewidth=3)
clb = fig.colorbar(lc, ax=ax0, orientation="vertical")
ax0.set_xlabel('$\\phi / \pi$', size=fontsize)
ax0.set_ylabel('$\\theta / \pi$', size=fontsize)
plt.sca(ax0)
plt.xticks(fontsize=fontsize*0.5)
plt.yticks(fontsize=fontsize*0.5)
ax1.plot(Table_t, Table_theta / np.pi, '-*', label='Table')
ax2.plot(Table_t, Table_phi / np.pi, '-*', label='Table')
ax3.plot(Table_t, Table_psi / np.pi, '-*', label='Table')
for axi, axyi in zip((ax1, ax2, ax3), ('$\\theta / \pi$', '$\\phi / \pi$', '$\\psi / \pi$')):
    plt.sca(axi)
    axi.set_xlabel('t', size=fontsize)
    axi.set_ylabel('%s' % axyi, size=fontsize)
    axi.legend()
    plt.xticks(fontsize=fontsize*0.5)
    plt.yticks(fontsize=fontsize*0.5)
plt.tight_layout()
ax0.set_xlim(phi_lim)
ax0.set_ylim(theta_lim)
ax1.set_ylim(theta_lim)
ax2.set_ylim(phi_lim)
ax3.set_ylim(psi_lim)
print('ini_psi', ini_psi)

# velocity of table results
fig, [(ax0, ax3), (ax1, ax4), (ax2, ax5), (ax6, ax7)] = plt.subplots(nrows=4, ncols=2, figsize=(20, 20))
fig.patch.set_facecolor('white')
ax0.plot(Table_t[:], Table_U[:, 0], '-*')
ax1.plot(Table_t[:], Table_U[:, 1], '-*')
ax2.plot(Table_t[:], Table_U[:, 2], '-*')
ax3.plot(Table_t[:], Table_U[:, 3], '-*')
ax4.plot(Table_t[:], Table_U[:, 4], '-*')
ax5.plot(Table_t[:], Table_U[:, 5], '-*')
ax6.plot(Table_t[:], np.linalg.norm(Table_U[:, :3], axis=1), '-*')
ax7.plot(Table_t[:], np.linalg.norm(Table_U[:, 3:], axis=1), '-*')
for axi, axyi in zip((ax0, ax1, ax2, ax3, ax4, ax5, ax6, ax7),
                     ('Table_u_x', 'Table_u_y', 'Table_u_z', 
                      'Table_w_x', 'Table_w_y', 'Table_w_z', 
                      '|Table_U|', '|Table_W|')):
    plt.sca(axi)
    axi.set_xlabel('time', size=fontsize)
    axi.set_ylabel('%s' % axyi, size=fontsize)
    plt.xticks(fontsize=fontsize*0.8)
    plt.yticks(fontsize=fontsize*0.8)
plt.tight_layout()


ShearTableProblem: run 2740 loops using 3.853821
ini_idx 0
ecoli_norm[ini_idx] [ 0.  0.  1.]
ini_psi 0.0