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)
%load_ext autoreload
%autoreload 2
from tqdm.notebook import tqdm as 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, signal
# 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 re
from scanf import scanf
import matplotlib
# matplotlib.use('agg')
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 tqdm import tqdm
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 [3]:
job_dir = 'ecoC01B05_passive_psi-0a'
table_name = 'ecoC01B05_tao1_wm0'
t_headle = '(.*?).pickle'
In [9]:
t_dir = os.path.join(PWD, job_dir)
data = spf_tb.load_table_data_pickle_dir(t_dir, t_headle)
lst_eta = data.lst_eta
theta_max_fre = data.theta_max_fre
phi_max_fre = data.phi_max_fre
psi_max_fre = data.psi_max_fre
eta_max_fre = data.eta_max_fre
data_idx = data.data_idx.fillna(-1).astype(int)
In [5]:
# check if the 3D parameter space is been traversed or not
do_check = False
if do_check:
t_headle = '(.*?).pickle'
t_path = os.listdir(os.path.join(PWD, job_dir))
filename_list = [filename for filename in os.listdir(os.path.join(PWD, job_dir))
if re.match(t_headle, filename) is not None]
ttheta = []
tphi = []
tpsi = []
for i0, tname in enumerate(tqdm_notebook(filename_list[:])):
tpath = os.path.join(PWD, job_dir, tname)
with open(tpath, 'rb') as handle:
tpick = pickle.load(handle)
ttheta.append(tpick['Table_theta'])
tphi.append(tpick['Table_phi'])
tpsi.append(tpick['Table_psi'])
ttheta, tphi, tpsi = np.hstack(ttheta), np.hstack(tphi), np.hstack(tpsi)
from evtk.hl import pointsToVTK
vtu_name = '/home/zhangji/check_th_ph_ps_%s' % job_dir
pointsToVTK(vtu_name,
ttheta, tphi, tpsi,
data={'th_ph_ps': ttheta})
print('save theta_phi_psi infomation to %s' % vtu_name)
In [6]:
# sort phase map of theta-phi using the name beging with pick frequience
importlib.reload(spf_tb)
do_sort = False
if do_sort:
# clear dir
dirpath = os.path.join(PWD, job_dir, 'th_ph_fft')
if os.path.exists(dirpath) and os.path.isdir(dirpath):
import shutil
shutil.rmtree(dirpath)
print('remove folder %s' % dirpath)
os.makedirs(dirpath)
print('make folder %s' % dirpath)
for tname in tqdm_notebook(filename_list[:]):
tpath = os.path.join(PWD, job_dir, tname)
with open(tpath, 'rb') as handle:
tpick = pickle.load(handle)
Table_t = tpick['Table_t']
Table_dt = tpick['Table_dt']
Table_X = tpick['Table_X']
Table_P = tpick['Table_P']
Table_P2 = tpick['Table_P2']
Table_theta = tpick['Table_theta']
Table_phi = tpick['Table_phi']
Table_psi = tpick['Table_psi']
Table_eta = tpick['Table_eta']
tmin = np.max((0, Table_t.max() - 1000))
idx = Table_t > tmin
freq_pk = spf_tb.get_major_fre(Table_t[idx], Table_theta[idx])
idx = Table_t > (Table_t.max() - 1 / freq_pk * 10)
freq_pk = spf_tb.get_major_fre(Table_t[idx], Table_theta[idx])
save_name = 'fre%.5f_%s.jpg' % (freq_pk, os.path.splitext(os.path.basename(tname))[0])
fig = spf_tb.light_save_theta_phi(os.path.join(dirpath, save_name),
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])
plt.close(fig)
In [5]:
importlib.reload(spf_tb)
theta, phi = 0.0, 1.714
tpick, _ = spf_tb.load_table_date_pickle(job_dir, theta, phi)
Table_t = tpick['Table_t']
Table_dt = tpick['Table_dt']
Table_X = tpick['Table_X']
Table_P = tpick['Table_P']
Table_P2 = tpick['Table_P2']
Table_theta = tpick['Table_theta']
Table_phi = tpick['Table_phi']
Table_psi = tpick['Table_psi']
Table_eta = tpick['Table_eta']
print('-ini_theta %f -ini_phi %f -ini_psi %f' %
(tpick['Table_theta'][0], tpick['Table_phi'][0], tpick['Table_psi'][0]))
freq_pk = spf_tb.get_major_fre(Table_t, Table_theta)
idx = Table_t > Table_t.max() - 1 / freq_pk * 10
# 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)
spf_tb.show_theta_phi(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], show_back_direction=False)
spf_tb.show_theta_phi_psi_eta(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])
spf_tb.show_center_X(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],
table_name=table_name)
Out[5]:
In [9]:
with pd.option_context('display.max_rows', 100, 'display.max_columns', 100):
display(data.theta_max_fre)
In [10]:
# sort all frequrents
with np.printoptions(precision=10, suppress=True, threshold=1e10):
print(np.flipud(np.sort(data.theta_max_fre.values.flatten())))
In [11]:
spf_tb.show_traj_phase_map_fre(theta_max_fre)
# spf_tb.show_traj_phase_map_fre(phi_max_fre)
# spf_tb.show_traj_phase_map_fre(psi_max_fre)
# spf_tb.show_traj_phase_map_fre(eta_max_fre)
Out[11]:
In [7]:
# put images with same frequence into a subdirect
importlib.reload(spf_tb)
def t_show_idx(iidx):
theta = type_fre.index.values[iidx[0][0]]
phi = type_fre.columns.values[iidx[1][0]]
print(theta, phi)
tipical_th_ph_list.append((theta, phi))
spf_tb.show_pickle_results(job_dir, theta, phi, table_name)
return True
def t_show_idx_list(iidx, nshow=5, Table_t_range1=np.array((0, np.inf)), Table_t_range2=np.array((0, np.inf))):
nshow = int(np.min((nshow, iidx[0].size)))
tidx = np.random.choice(iidx[0].size, nshow, replace=False)
theta = type_fre.index.values[iidx[0][tidx]]
phi = type_fre.columns.values[iidx[1][tidx]]
theta_phi_list = np.vstack((theta, phi)).T
# spf_tb.show_table_theta_phi_list(theta_phi_list, job_dir, Table_t_range=Table_t_range1, figsize=np.array((20, 20)), dpi=100)
spf_tb.show_table_result_list(theta_phi_list, job_dir, Table_t_range=Table_t_range2, figsize=np.array((16, 9)), dpi=100)
return True
tfre = theta_max_fre.copy()
check_fre_list = [0.0160, 0.0190, 0.0203, 0.0632]
atol_fre_list = [0.0001, 0.0001, 0.0008, 0.0008]
Table_t_range1 = np.array((0, 1000))
Table_t_range2 = np.array((4500, np.inf))
nshow = np.inf
tipical_th_ph_list = []
type_fre = tfre.copy()
type_fre.iloc[:, :] = len(check_fre_list)
for i0, (check_fre, atol_fre) in enumerate(zip(check_fre_list, atol_fre_list)):
use_idx = np.isclose(tfre, check_fre, rtol=0, atol=atol_fre)
type_fre.iloc[use_idx] = i0
iidx = np.where(use_idx)
# t_show_idx(iidx)
# t_show_idx_list(iidx, nshow=nshow, Table_t_range1=Table_t_range1, Table_t_range2=Table_t_range2)
# plot one of the remaind cases
if np.any(type_fre.values == len(check_fre_list)):
iidx = np.where(type_fre.values == len(check_fre_list))
# t_show_idx(iidx)
# t_show_idx_list(iidx, nshow=nshow, Table_t_range1=Table_t_range1, Table_t_range2=Table_t_range2)
spf_tb.show_traj_phase_map_type(type_fre)
# spf_tb.save_separate_angleList_fft(job_dir, tfre, check_fre_list, atol_fre_list)
Out[7]:
In [32]:
# create phase map
def tget_ax0():
n_xticks = 32
xticks = np.arange(n_xticks)
fig = plt.figure(figsize=(20, 20))
fig.patch.set_facecolor('white')
axs = []
axs.append(fig.add_subplot(221, polar=True))
axs.append(fig.add_subplot(222, polar=True))
axs.append(fig.add_subplot(223, polar=True))
axs.append(fig.add_subplot(224, polar=True))
for ax0 in axs:
ax0.set_xticks(xticks / n_xticks * 2 * np.pi)
ax0.set_xticklabels(['$\dfrac{%d}{%d}2\pi$' % (i0, n_xticks) for i0 in xticks])
ax0.set_yticklabels([])
ax0.set_ylim(0, np.pi)
plt.tight_layout()
return fig, axs
check_fre_list = [0.0160, 0.0190, 0.0203, 0.0632]
atol_fre_list = [0.0001, 0.0001, 0.0008, 0.0008]
color_list = [np.array(plt.get_cmap('tab20', len(check_fre_list))(i0)[:3]).reshape((1, 3))
for i0 in np.arange(len(check_fre_list))]
psi_lim_fct = 20
resampling_fct = 10
data0['use_max_fre'] = data0.theta_max_fre
case_path_list = spf_tb.separate_fre_path(check_fre_list, atol_fre_list, data0, pickle_path_list)
for idx, psi_lim1 in enumerate(np.linspace(0, 2 * np.pi, psi_lim_fct * 16,
endpoint=False)[::psi_lim_fct]):
fig, (ax0, ax1, ax2, ax3) = tget_ax0()
ax_list = [ax0, ax0, ax1, ax2, ax3]
psi_lim = (psi_lim1, psi_lim1 + 2 * np.pi / (psi_lim_fct * 16))
desc = '$\psi\in[%.3f\pi, %.3f\pi)$' % ((psi_lim[0] / np.pi), (psi_lim[1] / np.pi))
fig.suptitle(desc, fontsize=fontsize*0.8)
for check_fre, case_path, color, axi in zip(check_fre_list, case_path_list, color_list, ax_list):
thandle = '%f' % check_fre
spf_tb.draw_phase_map_theta(case_path, color, psi_lim, axs=(axi, ax_list[-1]), thandle=thandle,
resampling=True, resampling_fct=resampling_fct)
tdir = os.path.join(PWD, job_dir, 'phase_mape_fre')
if not os.path.exists(tdir):
os.makedirs(tdir)
figname = os.path.join(tdir, '%04d.png' % (idx))
fig.savefig(os.path.join(tdir, figname))
print('save to %s' % figname)
plt.close(fig)
In [ ]:
In [186]:
# show phase map of theta-phi, part 1
importlib.reload(spf_tb)
job_dir = 'ecoC01B05_T0.01_psi-0a'
table_name = 'ecoC01B05_T0.01'
t_headle = '(.*?).pickle'
t_path = os.listdir(os.path.join(PWD, job_dir))
filename_list = [filename for filename in os.listdir(os.path.join(PWD, job_dir))
if re.match(t_headle, filename) is not None]
ini_theta_list = []
ini_phi_list = []
lst_eta_list = []
pickle_path_list = []
idx_list = []
theta_primary_fre_list = []
phi_primary_fre_list = []
psi_primary_fre_list = []
for i0, tname in enumerate(tqdm_notebook(filename_list[:])):
tpath = os.path.join(PWD, job_dir, tname)
with open(tpath, 'rb') as handle:
tpick = pickle.load(handle)
ini_theta_list.append(tpick['ini_theta'])
ini_phi_list.append(tpick['ini_phi'])
lst_eta_list.append(tpick['Table_eta'][-1])
pickle_path_list.append(tpath)
idx_list.append(i0)
# fft rule
tx = tpick['Table_t']
tmin = np.max((0, tx.max() - 1000))
idx = tx > tmin
# the last frequence is the major frequence.
use_fft_number = 3
t1 = -use_fft_number - 1
theta_primary_fre_list.append(spf_tb.get_primary_fft_fre(tx[idx], tpick['Table_theta'][idx])[t1:-1])
phi_primary_fre_list.append(spf_tb.get_primary_fft_fre(tx[idx], tpick['Table_phi'][idx])[t1:-1])
psi_primary_fre_list.append(spf_tb.get_primary_fft_fre(tx[idx], tpick['Table_psi'][idx])[t1:-1])
In [216]:
# show phase map of theta-phi, part 2
def show_phase_map(tuse):
fig = plt.figure(figsize=(20, 12), dpi=300)
fig.patch.set_facecolor('white')
ax0 = fig.add_subplot(111, polar=True)
n_xticks = 32
xticks = np.arange(n_xticks)
ax0.set_xticks(xticks / n_xticks * 2 * np.pi)
ax0.set_xticklabels(['$\dfrac{%d}{%d}2\pi$' % (i0, n_xticks) for i0 in xticks])
ax0.set_yticklabels([])
ax0.set_ylim(0, np.pi)
tdata = tuse.values
im = ax0.pcolor(tuse.columns.values, tuse.index.values, tdata,
cmap=plt.get_cmap('Set2', np.nanmax(tdata)+1),
vmin=np.nanmin(tdata)-.5, vmax=np.nanmax(tdata)+.5)
ticks = np.arange(np.nanmin(tdata), np.nanmax(tdata)+1)
fig.colorbar(im, ax=ax0, orientation='vertical', ticks=ticks).ax.tick_params(labelsize=fontsize)
In [257]:
np.around(np.pi, 3)
Out[257]:
In [239]:
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
use_data = np.hstack((np.vstack(theta_primary_fre_list)[:, 1:],
np.vstack(phi_primary_fre_list)[:, 1:]))
#KMeans
km = KMeans(n_clusters=6, n_init=10, max_iter=1000, tol=1e-9, precompute_distances=True, n_jobs=-1, random_state=0)
km.fit(use_data)
km.predict(use_data)
tlabels = km.labels_
# # personal process
# tlabels[tlabels == 4] = 3
tdata1 = pd.DataFrame({'ini_theta': np.around(ini_theta_list, 3),
'ini_phi': np.around(ini_phi_list, 3),
'lst_eta': np.around(lst_eta_list, 3),
'use_fre': tlabels,
'data_idx': idx_list })
tdata1 = tdata1.pivot_table(index=['ini_theta'], columns=['ini_phi'])
show_phase_map(tdata1.use_fre)
In [250]:
np.array(ini_theta_list)[tlabels==show_tlabel]
Out[250]:
In [251]:
importlib.reload(spf_tb)
show_tlabel = 0
for theta, phi in zip(np.array(ini_theta_list)[tlabels==show_tlabel][:10],
np.array(ini_phi_list)[tlabels==show_tlabel][:10]):
spf_tb.show_pickle_results(job_dir, theta, phi, table_name, fast_mode=1)
In [241]:
importlib.reload(spf_tb)
for show_tlabel in np.arange(tlabels.max()+1):
theta = np.array(ini_theta_list)[tlabels==show_tlabel][0]
phi = np.array(ini_phi_list)[tlabels==show_tlabel][0]
spf_tb.show_pickle_results(job_dir, theta, phi, table_name, fast_mode=2)
In [ ]:
tlabels=