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

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.notebook import tqdm as tqdm_notebook
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 [4]:
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 [7]:
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)


-ini_theta 0.000000 -ini_phi 0.000000 -ini_psi 3.141593
Out[7]:
True

In [31]:
f


Out[31]:
array([0.00000000e+00, 1.71060853e-03, 3.42121706e-03, ...,
       7.00323133e+00, 7.00494194e+00, 7.00665255e+00])

In [39]:
figsize = np.array((16, 9)) * 1
dpi = 100

resamp_theta = spf_tb.get_continue_angle(Table_t, Table_theta)
fs = resamp_theta.size / Table_t.max() / 2
f, t, Zxx = signal.stft(np.cos(resamp_theta), fs=fs, nperseg=2**13)

fig, axs = plt.subplots(1, 1, figsize=figsize, dpi=dpi)
axi = axs
axi.pcolormesh(t, f, np.abs(Zxx))
axi.set_xlim(5000, 10000)
# axi.set_ylim(0, 0.1)
# axi.set_yscale('log')


Out[39]:
[<matplotlib.lines.Line2D at 0x7f7a9f115b38>]

In [102]:
14 / 0.016 * 4


Out[102]:
3500.0

In [143]:
figsize = np.array((16, 9)) * 0.5
dpi = 100
show_prim_freq = 3

use_t = np.linspace(Table_t.min(), Table_t.max(), 1 * Table_t.size)
# resamp_theta = spf_tb.get_continue_angle(Table_t, Table_theta, use_t)
taaa = np.random.sample(2)
print(taaa)
resamp_theta = np.cos(use_t * taaa[0] * 2 * np.pi) + np.cos(use_t * taaa[1] * 2 * np.pi)
ty = np.cos(resamp_theta)
# ty = resamp_theta

tmin = np.max((0, use_t.max() - 1000))
idx = use_t > tmin
freq_pk = spf_tb.get_major_fre(use_t[idx], ty[idx])
idx = use_t > (Table_t.max() - 1 / freq_pk * 30)
# tfft = np.fft.rfft(ty[idx])

fs = ty[idx].size / (use_t[idx].max() - use_t[idx].min())
nperseg = fs / freq_pk * 8
f2, Pxx_den = signal.welch(ty[idx], fs=fs, nperseg=nperseg)

fig, axs = plt.subplots(1, 2, figsize=figsize, dpi=dpi)
axi = axs[0]
axi.plot(use_t[idx], resamp_theta[idx], '.-')

axi = axs[1]
axi.loglog(f2, Pxx_den, '.')
tpk = signal.find_peaks(Pxx_den)[0]
fft_abs_pk = Pxx_den[tpk]
freq_pk = f2[tpk]
tidx = np.argsort(fft_abs_pk)[-show_prim_freq:]
axi.loglog(freq_pk[tidx], fft_abs_pk[tidx], '*')
t1 = 'starred freq: \n' + '\n'.join(['$%.5f$' % freq_pk[ti] for ti in tidx])
axi.text(axi.get_xlim()[0] * 1.1, axi.get_ylim()[0] * 2, t1)


[0.4617505  0.40895383]
Out[143]:
Text(0.005145975233880607, 4.935046782540978e-19, 'starred freq: \n$0.92107$\n$0.86806$\n$0.05301$')

In [131]:
0.949 / 0.15


Out[131]:
6.326666666666666

In [89]:
tpk = signal.find_peaks(ty[idx])[0]
fft_abs_pk = tfft_abs[tpk]
freq_pk = tfreq[tpk]
tidx = np.argsort(fft_abs_pk)[-show_prim_freq:]
ax1.loglog(freq_pk[tidx], fft_abs_pk[tidx], '*')


Out[89]:
array([ 380, 1276, 2172, 3068, 3964, 4860, 5756, 6652, 7548, 8444])

In [ ]:


In [5]:
# put images with same frequence into a subdirect
check_fre_list = [0.0160, 0.0190, 0.0203, 0.0632]
atol_fre_list =  [0.0001, 0.0001, 0.0008, 0.0008]
Table_t_range = np.array((4500, np.inf))

def sub_seperate_0(i1, type_fre, iidx, Table_t_range):
    t1 = np.zeros_like(iidx[0])
    return t1

def sub_seperate_1(i1, type_fre, iidx, Table_t_range):
    theta = type_fre.index.values[iidx[0]]
    phi = type_fre.columns.values[iidx[1]]
    theta_phi_list = np.vstack((theta, phi)).T
    tuse = []
    for theta, phi in tqdm_notebook(theta_phi_list, desc='No. %d  ' % i1):
        tpick, _ = spf_tb.load_table_date_pickle(job_dir, theta, phi)
        Table_t = tpick['Table_t']
        Table_theta = tpick['Table_theta']
        Table_phi = tpick['Table_phi']
        Table_psi = tpick['Table_psi']
        
        idx = np.logical_and(Table_t >= Table_t_range[0], Table_t <= Table_t_range[1])
        tuse.append(np.max(Table_theta[idx]))
    tuse = np.hstack(tuse)
    t1 = np.ones_like(tuse)
    t1[tuse < 1.6] = 0
    t1[tuse > 1.934] = 2
    return t1

def sub_seperate_2(i1, type_fre, iidx, Table_t_range):
    theta = type_fre.index.values[iidx[0]]
    phi = type_fre.columns.values[iidx[1]]
    theta_phi_list = np.vstack((theta, phi)).T
    tuse = []
    for theta, phi in tqdm_notebook(theta_phi_list, desc='No. %d  ' % i1):
        tpick, _ = spf_tb.load_table_date_pickle(job_dir, theta, phi)
        Table_t = tpick['Table_t']
        Table_theta = tpick['Table_theta']
        Table_phi = tpick['Table_phi']
        Table_psi = tpick['Table_psi']
        
        idx = np.logical_and(Table_t >= Table_t_range[0], Table_t <= Table_t_range[1])
        tuse.append(np.max(Table_phi[idx]))
    tuse = np.hstack(tuse)
    t1 = np.ones_like(tuse)
    t1[tuse < 4] = 0
    t1[tuse > 5.1] = 2
    return t1
    
sub_seperate = {0: sub_seperate_0, 
                1: sub_seperate_1, 
                2: sub_seperate_2, 
                3: sub_seperate_0, }

tfre = theta_max_fre.copy()
type_fre = tfre.copy()
type_fre.iloc[:, :] = -1
i0 = 0
process_total = 0
for i1, (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)
    iidx = np.where(use_idx)
    process_total = process_total + np.sum(use_idx)
    t1 = sub_seperate[i1](i1, type_fre, iidx, Table_t_range)
    type_fre.iloc[use_idx] = i0 + t1
    i0 = i0 + 1 + np.max(t1)

# plot one of the remaind cases
if np.any(type_fre.values == -1):
    pass

assert np.any(type_fre.values > 0)

type_fre2 = type_fre.copy()
type_fre2.iloc[type_fre.values == 0] = 0
type_fre2.iloc[type_fre.values == 1] = 1
type_fre2.iloc[type_fre.values == 2] = 2
type_fre2.iloc[type_fre.values == 3] = 3
type_fre2.iloc[type_fre.values == 4] = 1
type_fre2.iloc[type_fre.values == 5] = 1
type_fre2.iloc[type_fre.values == 6] = 4
type_fre2.iloc[type_fre.values == 7] = 1
spf_tb.show_traj_phase_map_type(type_fre2)
# spf_tb.save_separate_angleList_fft(job_dir, tfre, check_fre_list, atol_fre_list)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-8711920b5708> in <module>
     53                 3: sub_seperate_0, }
     54 
---> 55 tfre = theta_max_fre.copy()
     56 type_fre = tfre.copy()
     57 type_fre.iloc[:, :] = -1

NameError: name 'theta_max_fre' is not defined

In [7]:
t1 = []
for i0 in np.arange(type_fre2.values.max() + 1):
    t1.append(np.isclose(type_fre2.values, i0).sum())

print(np.sum(t1))
print(t1)
print(t1 / np.sum(t1))


1035
[275, 47, 232, 320, 161]
[0.26570048 0.04541063 0.22415459 0.30917874 0.15555556]

In [10]:
nshow = np.inf
# nshow = 2
fast_mode1 = 1
figsize = np.array((16, 9)) * 0.5
dpi = 200

tipical_th_ph_list = []
for i0 in np.arange(type_fre2.values.max() + 1)[:]:
    iidx = np.where(np.isclose(type_fre2.values, i0))
    spf_tb.phase_map_show_idx(type_fre2, tipical_th_ph_list, iidx, job_dir, table_name, fast_mode=fast_mode1, )


-ini_theta 0.000000 -ini_phi 1.714000
-ini_theta 0.000000 -ini_phi 0.000000 -ini_psi 3.141593
-ini_theta 0.000000 -ini_phi 0.000000
-ini_theta 0.000000 -ini_phi 0.000000 -ini_psi 0.000000
-ini_theta 0.143000 -ini_phi 0.000000
-ini_theta 0.142800 -ini_phi 6.283185 -ini_psi 0.000000
-ini_theta 0.143000 -ini_phi 0.428000
-ini_theta 0.142800 -ini_phi 0.428399 -ini_psi 0.000000
-ini_theta 0.143000 -ini_phi 0.571000
-ini_theta 0.142800 -ini_phi 0.571199 -ini_psi 0.000000