In [1]:
# import numpy as np
# # !/usr/bin/env python3
# # -*- coding: utf-8 -*-
# """
# Created on 20181219
# @author: zhangji
# Trajection of a ellipse, Jeffery equation.
# """
# %pylab inline
# pylab.rcParams['figure.figsize'] = (25, 11)
# fontsize = 40
# import numpy as np
# import scipy as sp
# from scipy.optimize import leastsq, curve_fit
# from scipy import interpolate
# from scipy.interpolate import interp1d
# from scipy.io import loadmat, savemat
# # import scipy.misc
# import matplotlib
# from matplotlib import pyplot as plt
# from matplotlib import animation, rc
# import matplotlib.ticker as mtick
# from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
# from mpl_toolkits.mplot3d import Axes3D, axes3d
# from sympy import symbols, simplify, series, exp
# from sympy.matrices import Matrix
# from sympy.solvers import solve
# from IPython.display import display, HTML
# from tqdm import tqdm_notebook as tqdm
# import pandas as pd
# import re
# from scanf import scanf
# import os
# import glob
# from codeStore import support_fun as spf
# from src.support_class import *
# from src import stokes_flow as sf
# rc('animation', html='html5')
# PWD = os.getcwd()
# font = {'size': 20}
# matplotlib.rc('font', **font)
# np.set_printoptions(linewidth=90, precision=5)
from tqdm import tqdm_notebook
import os
import glob
import natsort
import numpy as np
import scipy as sp
from scipy.optimize import leastsq, curve_fit
from scipy import interpolate, integrate
from scipy import spatial
# from scipy.interpolate import interp1d
from scipy.io import loadmat, savemat
# import scipy.misc
import importlib
from IPython.display import display, HTML
import pandas as pd
import pickle
import matplotlib
from matplotlib import pyplot as plt
import matplotlib.colors as colors
from matplotlib import animation, rc
import matplotlib.ticker as mtick
from mpl_toolkits.axes_grid1.inset_locator import inset_axes, zoomed_inset_axes
from mpl_toolkits.mplot3d import Axes3D, axes3d
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
from mpl_toolkits.mplot3d.art3d import Line3DCollection
from matplotlib import cm
from time import time
from src.support_class import *
from src import jeffery_model as jm
from codeStore import support_fun as spf
from codeStore import support_fun_table as spf_tb
# %matplotlib notebook
%matplotlib inline
rc('animation', html='html5')
fontsize = 40
PWD = os.getcwd()
In [2]:
fig = plt.figure(figsize=(2, 2))
fig.patch.set_facecolor('white')
ax0 = fig.add_subplot(1, 1, 1)
In [43]:
# active ecoli petsc family method
importlib.reload(spf_tb)
t0 = time()
t_theta, t_phi, t_psi = np.random.sample(3) * (1, 2, 2) * np.pi
# t_theta, t_phi, t_psi = 1.491744, 0.268875, 0.097399
# t_theta, t_phi, t_psi = 1.9643076668745707, 2.127873750065894, 1.86195827972309
# t_theta, t_phi, t_psi = 1.9182367797144864, 3.6042598037667752, 0.019126035606903825
t_theta, t_phi, t_psi = 0, 0, 0
max_t = 300
update_fun='3bs'
rtol=1e-6
atol=1e-9
eval_dt = 0.003
save_every = 1
# table_name='planeShearRatex_4c'
# omega_tail=48.41664953
# table_name='planeShearRatex_1d'
# omega_tail=193.66659814
# table_name='planeShearRatex_0.25c'
# omega_tail=774.66639255
# table_name='hlxC01_tau1a'
# omega_tail=0
# table_name='hlxB01_tau1a'
# omega_tail=0
# table_name='planeShearRatex_-1c'
# omega_tail=-193.66659814
# table_name='planeShearRatex_1d_passive'
# omega_tail=0
# table_name='ellipse_alpha3'
# omega_tail=0
table_name='ecoC01B01_tau1c'
# table_name='ecoC01B05_T193'
omega_tail=193.66659814
tnorm = np.array((np.sin(t_theta) * np.cos(t_phi), np.sin(t_theta) * np.sin(t_phi), np.cos(t_theta)))
Table_t, Table_dt, Table_X, Table_P, Table_P2, Table_theta, Table_phi, Table_psi, Table_eta \
= spf_tb.do_calculate_ecoli_Petsc4n(tnorm, t_psi, max_t, update_fun=update_fun, rtol=rtol, atol=atol,
eval_dt=eval_dt, save_every=save_every,
table_name=table_name, omega_tail=omega_tail)
t1 = time()
print('init norm: ', t_theta, t_phi, t_psi)
last_norm = np.array((np.sin(Table_theta[-1]) * np.cos(Table_phi[-1]),
np.sin(Table_theta[-1]) * np.sin(Table_phi[-1]),
np.cos(Table_theta[-1])))
print('last norm: ', last_norm, 'last psi:', Table_psi[-1])
print('%s: run %d loops/times using %fs' % ('do_calculate_ecoli_Petsc', max_t, (t1 - t0)))
print('%s_%s rt%.0e, at%.0e, dt%.0e %.1fs' % ('PETSC RK', update_fun, rtol, atol, eval_dt, (t1 - t0)))
spf_tb.show_theta_phi(Table_t, Table_dt, Table_X, Table_P, Table_P2,
Table_theta, Table_phi, Table_psi, Table_eta)
spf_tb.show_theta_phi_psi_eta(Table_t, Table_dt, Table_X, Table_P, Table_P2,
Table_theta, Table_phi, Table_psi, Table_eta, resampling_fct=2)
spf_tb.show_center_X(Table_t, Table_dt, Table_X, Table_P, Table_P2,
Table_theta, Table_phi, Table_psi, Table_eta,
table_name=table_name)
print(Table_t.size)
# idx = Table_t < 2000
# spf_tb.show_table_result(Table_t[idx], Table_dt[idx], Table_X[idx], Table_P[idx], Table_P2[idx],
# Table_theta[idx], Table_phi[idx], Table_psi[idx], Table_eta[idx], save_every)
In [51]:
print('%f, %f, %f' % (Table_theta[-1], Table_phi[-1], Table_psi[-1]))
In [57]:
# %matplotlib notebook
%matplotlib inline
pickle_name = table_name + '_kwargs'
problem_kwargs = spf_tb.load_problem_kwargs(pickle_name)
tnode1, tnode2 = get_ecoli_nodes_split_at(Table_theta[-1], Table_phi[-1], Table_psi[-1],
now_center=np.zeros(3), **problem_kwargs)
fig = plt.figure(figsize=(5, 5))
fig.patch.set_facecolor('white')
ax0 = fig.add_subplot(111, projection='3d')
for spine in ax0.spines.values():
spine.set_visible(False)
ax0.plot(tnode1[:, 0], tnode1[:, 1], tnode1[:, 2])[0]
ax0.plot(tnode2[:, 0], tnode2[:, 1], tnode2[:, 2])[0]
ax0.set_xlabel('X')
ax0.set_ylabel('Y')
ax0.set_zlabel('Z')
Out[57]:
In [58]:
# passive ecoli petsc family method
importlib.reload(spf_tb)
importlib.reload(jm)
t0 = time()
# t_theta, t_phi, t_psi = np.pi / 2, np.pi / 2, 0
t_theta, t_phi, t_psi = 0, 0, 0
max_t = 100
update_fun='5bs'
rtol=1e-12
atol=1e-15
eval_dt = 0.01
# save_every = np.ceil(1 / eval_dt / 100)
save_every = 1
# table_name='ecoC01_tau1c'
# omega_tail=193.66659814
t_theta, t_phi, t_psi = 1.5519607818186965 , 3.1554179700918383 , 1.1811466563914905
table_name='ecoC01B05_tau1c_passive'
omega_tail=0
tnorm = np.array((np.sin(t_theta) * np.cos(t_phi), np.sin(t_theta) * np.sin(t_phi), np.cos(t_theta)))
t_psi = np.ones(1) * t_psi
Table_t, Table_dt, Table_X, Table_P, Table_P2, Table_theta, Table_phi, Table_psi, Table_eta\
= spf_tb.do_calculate_ecoli_passive_Petsc4n(tnorm, t_psi, max_t, update_fun=update_fun, rtol=rtol, atol=atol,
eval_dt=eval_dt, save_every=save_every,
table_name=table_name, omega_tail=omega_tail)
t1 = time()
print('last norm: ', Table_theta[-1], ',', Table_phi[-1], ',', Table_psi[-1])
print('%s: run %d loops/times using %fs' % ('do_calculate_ecoli_Petsc4nPsi', max_t, (t1 - t0)))
print('%s_%s rt%.0e, at%.0e, dt%.0e %.1fs' % ('PETSC RK', update_fun, rtol, atol, eval_dt, (t1 - t0)))
spf_tb.show_table_result(Table_t, Table_dt, Table_X, Table_P, Table_P2,
Table_theta, Table_phi, Table_psi, Table_eta, save_every)
# t_pick = (t_theta, t_phi, t_psi, max_t, update_fun, rtol, atol, eval_dt,
# Table_t, Table_X, Table_P, Table_P2, Table_theta, Table_phi, Table_psi, Table_eta, save_every)
# idx = np.load('../motion_ecoliB01_t able/idx.npy')
# t_name = 'idx%03d_th%5.3f_ph%5.3f_ps%5.3f.pickle' % (idx, t_theta, t_phi, t_psi)
# np.save('../motion_ecoliB01_table/idx.npy', (idx + 1))
# with open('../motion_ecoliB01_table/%s' % t_name, 'wb') as handle:
# pickle.dump(t_pick, handle, protocol=pickle.HIGHEST_PROTOCOL)
# print('save to %s' % t_name)
Out[58]:
In [204]:
importlib.reload(spf_tb)
tidx = np.arange(Table_t.size)
tidx = np.arange(Table_t.size-10000, Table_t.size-2)
anim = spf_tb.make_table_video(Table_t[tidx], Table_X[tidx], Table_P[tidx], Table_P2[tidx],
Table_theta[tidx], Table_phi[tidx], Table_psi[tidx], Table_eta[tidx],
zm_fct=30, stp=1, interval=1, trange=70)
tname = '%s_th%5.3f_ph%5.3f_ps%5.3f.mp4' % (table_name, t_theta, t_phi, t_psi)
anim
# anim.save(tname, writer='ffmpeg', fps=15)
# print(tname, datetime.now())
Out[204]: