Skip to content

Commit

Permalink
Fixed PSD_FILTER_W definition in resample_current. Fixed minor
Browse files Browse the repository at this point in the history 
copy-paste errors with write_log functions.
  • Loading branch information
@lorisercole
lorisercole committed Jul 19, 2019
1 parent 05b1c17 commit 6f76ef6
Show file tree
Hide file tree
Showing 10 changed files with 80 additions and 65 deletions.
Empty file modified examples/example_commandline_NaCl/command.sh
100644 → 100755
View file
Empty file.
9 changes: 5 additions & 4 deletions thermocepstrum/analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,11 @@
except ImportError:
raise ImportError('Cannot locate thermocepstrum.')

# import log-print method
from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()
log.set_method('bash')

#try to import plotstyle (not fundamental)
try:
import pkg_resources
Expand All @@ -39,10 +44,6 @@
plt.style.use(pltstyle_filename)
except:
pass
from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()

log.set_method('bash')

c = plt.rcParams['axes.prop_cycle'].by_key()['color']

Expand Down
71 changes: 38 additions & 33 deletions thermocepstrum/heatcurrent.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,16 +313,16 @@ def plot_cepstral_spectrum(self, freq_units='thz', freq_scale=1.0, axes=None, ka
axes[1].grid()
return axes

def resample_current(self, TSKIP=None, fstar_THz=None, FILTER_W=None, plot=True, PSD_FILTER_W=None, freq_units='thz',
FIGSIZE=None):
return resample_current(self, TSKIP=TSKIP, fstar_THz=fstar_THz, FILTER_W=FILTER_W, plot=plot, PSD_FILTER_W=PSD_FILTER_W,
freq_units=freq_units)
def resample_current(self, TSKIP=None, fstar_THz=None, FILTER_W=None, plot=True, PSD_FILTER_W=None,
req_units='thz', FIGSIZE=None): # yapf: disable
return resample_current(self, TSKIP=TSKIP, fstar_THz=fstar_THz, FILTER_W=FILTER_W, plot=plot,
PSD_FILTER_W=PSD_FILTER_W, freq_units=freq_units) # yapf: disable


#is this function needed?
# def compute_kappa_multi(self, others, FILTER_WINDOW_WIDTH=None):
# """Multi-component kappa calculation."""
# log.write_log "HeatCurrent.compute_kappa_multi"
# log.write_log("HeatCurrent.compute_kappa_multi")
# if FILTER_WINDOW_WIDTH is None:
# FILTER_WINDOW_WIDTH = self.FILTER_WINDOW_WIDTH
# multi_mdsample = super(HeatCurrent, self).compute_kappa_multi(others, FILTER_WINDOW_WIDTH, DT_FS=self.DT_FS, call_other=True)
Expand Down Expand Up @@ -375,32 +375,34 @@ def resample_current(x, TSKIP=None, fstar_THz=None, FILTER_W=None, plot=True, PS
TSKIP = int(round(x.Nyquist_f_THz / fstar_THz))
if plot:
figure, axes = plt.subplots(2, sharex=True, figsize=FIGSIZE)
axes = x.plot_periodogram(PSD_FILTER_W, freq_units, 1.0, axes=axes)
axes = x.plot_periodogram(PSD_FILTER_W, freq_units, 1.0, axes=axes) # this also updates x.FILTER_WINDOW_WIDTH

fstar_THz = x.Nyquist_f_THz / TSKIP
fstar_idx = np.argmin(x.freqs_THz < fstar_THz)

# filter and sample
if FILTER_W is None:
FILTER_W = TSKIP
trajf = md.tools.filter_and_sample(x.traj, FILTER_W, TSKIP, 'rectangular')

# resample filtering window width in order to use the same filtering frequency window in the plot
# if PSD_FILTER_W was specified, then x.FILTER_WINDOW_WIDTH was updated by the previous plot function
if x.FILTER_WINDOW_WIDTH is not None:
PSD_FILTER_W = x.FILTER_WINDOW_WIDTH * TSKIP

# define new HeatCurrent
if not x.multicomponent:
if x.FILTER_WINDOW_WIDTH is not None:
xf = HeatCurrent(trajf, x.units, x.DT_FS * TSKIP, x.TEMPERATURE, x.VOLUME, PSD_FILTER_W= PSD_FILTER_W, freq_units=freq_units)
else:
xf = HeatCurrent(trajf, x.units, x.DT_FS * TSKIP, x.TEMPERATURE, x.VOLUME)
xf = HeatCurrent(trajf, x.units, x.DT_FS * TSKIP, x.TEMPERATURE, x.VOLUME, PSD_FILTER_W, freq_units)
else:
if x.otherMD is None:
raise ValueError('x.otherMD cannot be none (wrong/missing initialization?)')
raise RuntimeError('x.otherMD cannot be none (wrong/missing initialization?)')
# filter_and_sample also other trajectories
yf = []
yf.append(trajf)
for y in x.otherMD:
tmp = md.tools.filter_and_sample(y.traj, FILTER_W, TSKIP, 'rectangular')
yf.append(tmp)
if x.FILTER_WINDOW_WIDTH is not None:
xf = HeatCurrent(yf, x.units, x.DT_FS * TSKIP, x.TEMPERATURE, x.VOLUME, PSD_FILTER_W= PSD_FILTER_W, freq_units=freq_units)
else:
xf = HeatCurrent(yf, x.units, x.DT_FS * TSKIP, x.TEMPERATURE, x.VOLUME)
xf = HeatCurrent(yf, x.units, x.DT_FS * TSKIP, x.TEMPERATURE, x.VOLUME, PSD_FILTER_W, freq_units)
if plot:
if (freq_units == 'thz') or (freq_units == 'THz'):
xf.plot_periodogram(x.FILTER_WINDOW_WIDTH * 1000. / x.DT_FS, 'thz', TSKIP, axes=axes)
Expand All @@ -409,26 +411,29 @@ def resample_current(x, TSKIP=None, fstar_THz=None, FILTER_W=None, plot=True, PS
log.write_log(x.FILTER_WINDOW_WIDTH)
xf.plot_periodogram(x.FILTER_WINDOW_WIDTH * TSKIP, 'red', TSKIP, axes=axes)

xf.resample_log = '-----------------------------------------------------\n' +\
' RESAMPLE TIME SERIES\n' +\
'-----------------------------------------------------\n' +\
' Original Nyquist freq f_Ny = {:12.5f} THz\n'.format(x.Nyquist_f_THz) +\
' Resampling freq f* = {:12.5f} THz\n'.format(fstar_THz) +\
' Sampling time TSKIP = {:12d} steps\n'.format(TSKIP) +\
' = {:12.3f} fs\n'.format(TSKIP * x.DT_FS) +\
' Original n. of frequencies = {:12d}\n'.format(x.Nfreqs) +\
' Resampled n. of frequencies = {:12d}\n'.format(xf.Nfreqs)
# write log
xf.resample_log = \
'-----------------------------------------------------\n' +\
' RESAMPLE TIME SERIES\n' +\
'-----------------------------------------------------\n' +\
' Original Nyquist freq f_Ny = {:12.5f} THz\n'.format(x.Nyquist_f_THz) +\
' Resampling freq f* = {:12.5f} THz\n'.format(fstar_THz) +\
' Sampling time TSKIP = {:12d} steps\n'.format(TSKIP) +\
' = {:12.3f} fs\n'.format(TSKIP * x.DT_FS) +\
' Original n. of frequencies = {:12d}\n'.format(x.Nfreqs) +\
' Resampled n. of frequencies = {:12d}\n'.format(xf.Nfreqs)
if x.fpsd is not None:
xf.resample_log = xf.resample_log + ' PSD @cutoff (pre-filter) = {:12.5f}\n'.format(x.fpsd[fstar_idx]) +\
' (post-filter) = {:12.5f}\n'.format(xf.fpsd[-1]) +\
' log(PSD) @cutoff (pre-filter) = {:12.5f}\n'.format(x.flogpsd[fstar_idx]) +\
' (post-filter) = {:12.5f}\n'.format(xf.flogpsd[-1]) +\
' min(PSD) (pre-filter) = {:12.5f}\n'.format(x.psd_min) +\
' min(PSD) (post-filter) = {:12.5f}\n'.format(xf.psd_min) +\
' % of original PSD Power f<f* (pre-filter) = {:5f}\n'.format(np.trapz(x.psd[:fstar_idx+1]) / x.psd_power * 100.)
xf.resample_log += \
' PSD @cutoff (pre-filter) = {:12.5f}\n'.format(x.fpsd[fstar_idx]) +\
' (post-filter) = {:12.5f}\n'.format(xf.fpsd[-1]) +\
' log(PSD) @cutoff (pre-filter) = {:12.5f}\n'.format(x.flogpsd[fstar_idx]) +\
' (post-filter) = {:12.5f}\n'.format(xf.flogpsd[-1]) +\
' min(PSD) (pre-filter) = {:12.5f}\n'.format(x.psd_min) +\
' min(PSD) (post-filter) = {:12.5f}\n'.format(xf.psd_min) +\
' % of original PSD Power f<f* (pre-filter) = {:5f}\n'.format(np.trapz(x.psd[:fstar_idx+1]) / x.psd_power * 100.)
else:
xf.resample_log = xf.resample_log + ' fPSD not calculated before resampling!\n '
xf.resample_log = xf.resample_log + '-----------------------------------------------------\n'
xf.resample_log += ' fPSD not calculated before resampling!\n '
xf.resample_log += '-----------------------------------------------------\n'
log.write_log(xf.resample_log)

if plot:
Expand Down
7 changes: 4 additions & 3 deletions thermocepstrum/i_o/read_lammps_dump.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()


def is_string(string):
try:
float(string)
Expand Down Expand Up @@ -78,7 +79,7 @@ class LAMMPS_Dump(object):
traj.read_timesteps(10, select_ckeys=['id,xu,yu,vu']) -->> Read the next 10 timesteps, only the specified columns (DELTA_TIMESTEP is assumed)
traj.read_timesteps((10,30)) -->> Read from TIMESTEP 10 to 30
traj.read_timesteps((10,30,2)) -->> Read every 2 steps from TIMESTEP 10 to 30
log.write_log traj.data
print(traj.data)
"""

def __init__(self, *args, **kwargs):
Expand All @@ -104,7 +105,7 @@ def __init__(self, *args, **kwargs):
self._read_ckeys(group_vectors, preload_timesteps)
self.ckey = None
#self.MAX_NSTEPS = data_length(self.filename)
#log.write_log "Data length = ", self.MAX_NSTEPS
#log.write_log("Data length = ", self.MAX_NSTEPS)
return

def __repr__(self):
Expand Down Expand Up @@ -409,7 +410,7 @@ def read_timesteps(self, selection, start_step=-1, select_ckeys=None, fast_check
progbar.description = '%g %%' % progbar.value
else:
log.write_log(' step = {:9d} - {:6.2f}% completed'.format(istep + 1,
float(istep + 1) / self.nsteps * 100.))
float(istep + 1) / self.nsteps * 100.))
if self._GUI:
progbar.close()
# check number of steps read, keep an even number of steps
Expand Down
10 changes: 6 additions & 4 deletions thermocepstrum/i_o/read_lammps_log.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
### Example script:
### data = LAMMPSLogFile('lammps.log', run_keyword='PRODUCTION RUN')
### data.read_datalines(NSTEPS=100, start_step=0, select_ckeys=['Step', 'Temp', 'flux'])
### log.write_log data.data
### print(data.data)
###
### # to save data into a Numpy file:
### save_hc_npz(data, ['flux'], 'lammps.data', 'flux.npz')
Expand All @@ -51,6 +51,7 @@
from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()


def is_string(string):
try:
float(string)
Expand Down Expand Up @@ -120,7 +121,7 @@ class LAMMPSLogFile(object):
Example script:
data = LAMMPSLogFile(filename, run_keyword='PRODUCTION RUN')
data.read_datalines(NSTEPS=100, start_step=0, select_ckeys=['Step', 'Temp', 'flux'])
log.write_log data.data
print(data.data)
# to save data into a Numpyz file:
save_hc_npz(data, ['flux'], 'lammps.data', 'flux.npz')
Expand Down Expand Up @@ -326,7 +327,8 @@ def read_datalines(self, NSTEPS=0, start_step=-1, select_ckeys=None, max_vector_
progbar.value = float(step + 1) / NSTEPS * 100.
progbar.description = '{:6.2f}%'.format(progbar.value)
else:
log.write_log(' step = {:9d} - {:6.2f}% completed'.format(step + 1, float(step + 1) / NSTEPS * 100.))
log.write_log(' step = {:9d} - {:6.2f}% completed'.format(step + 1,
float(step + 1) / NSTEPS * 100.))

if self._GUI:
progbar.close()
Expand Down Expand Up @@ -400,7 +402,7 @@ def get_box(filename):
if 'Time' in lammpslogfile.ckey:
dic['DT_TIMEUNITS'] = lammpslogfile.data['Time'][1, 0] - lammpslogfile.data['Time'][0, 0]

log.write_log("These keys will be saved in file \"{:}\" :".format(outfilename))
log.write_log('These keys will be saved in file \"{:}\" :'.format(outfilename))
log.write_log(' ', list(dic.keys()))
np.savez(outfilename, **dic)
return
Expand Down
8 changes: 5 additions & 3 deletions thermocepstrum/i_o/read_tablefile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,14 +30,15 @@
### example:
### current = TableFile(filename)
### current.read_datalines(NSTEPS=100, start_step=0, select_ckeys=['Step', 'Temp', 'flux'])
### log.write_log current.data
### print(current.data)
################################################################################

import numpy as np
from time import time
from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()


def is_string(string):
try:
float(string)
Expand Down Expand Up @@ -78,7 +79,7 @@ class TableFile(object):
Example:
current = TableFile(filename)
current.read_datalines(NSTEPS=100, select_ckeys=['Step', 'Temp', 'flux'])
log.write_log current.data
print(current.data)
Variables (columns) are organized into a dictionary according to the column headers.
LAMMPS-style vector variables header are grouped together (only if group_vector = True).
Expand Down Expand Up @@ -270,7 +271,8 @@ def read_datalines(self, NSTEPS=0, start_step=-1, select_ckeys=None, max_vector_
progbar.value = float(step + 1) / NSTEPS * 100.
progbar.description = '{:6.2f}%'.format(progbar.value)
else:
log.write_log(' step = {:9d} - {:6.2f}% completed'.format(step + 1, float(step + 1) / NSTEPS * 100.))
log.write_log(' step = {:9d} - {:6.2f}% completed'.format(step + 1,
float(step + 1) / NSTEPS * 100.))

if self._GUI:
progbar.close()
Expand Down
2 changes: 1 addition & 1 deletion thermocepstrum/md/cepstral.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -292,5 +292,5 @@ def compute_logtau_density(self,
# else:
# self.optimalK_idx = np.NaN
# self.optimalK = np.NaN
# log.write_log 'Warning: optimal cutoff K NOT FOUND.'
# log.write_log(Warning: optimal cutoff K NOT FOUND.')
# return
2 changes: 1 addition & 1 deletion thermocepstrum/md/tools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ def resample_psd(freqs, psd, cutfrequency):
newpsd = psd.copy()[:cutidx + 1]
newpsd = newpsd + psd[:-cutidx - 2:-1]
newfreqs = freqs[:cutidx + 1]
#log.write_log cutidx, DT, freqs[cutidx], newpsd.size
#log.write_log(cutidx, DT, freqs[cutidx], newpsd.size)
else:
raise NotImplementedError('Not implemented.')
return newfreqs, newpsd
12 changes: 7 additions & 5 deletions thermocepstrum/utils/blockanalysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -395,22 +395,24 @@ def compute_averages(self):
self.flogtau_density_XAVE, self.flogtau_density_XSTD,
np.mean(np.mean(list(self.flogtau_density_xstd()))), np.std(list(self.flogtau_density_xstd()))))
log.write_log(' flogtau[ave AIC_w] = {:12f} +/- {:8f}'.format(self.flogtau_avedensity_XAVE,
self.flogtau_avedensity_XSTD))
self.flogtau_avedensity_XSTD))

self.ftau_density_XAVE = np.mean(list(self.ftau_density_xave()))
self.ftau_density_XSTD = np.std(list(self.ftau_density_xave()))
#self.ftau_avedensity = np.mean(list(self.ftau_density()), axis=1)
#self.ftau_avedensity_XAVE, self.ftau_avedensity_XSTD = ta.grid_statistics(self.density_grid, self.ftau_avedensity)
log.write_log(' ftau[AIC_w] = {:12f} +/- {:8f}'.format(self.ftau_density_XAVE, self.ftau_density_XSTD))
#log.write_log ' ftau[ave AIC_w] = {:12f} +/- {:8f}'.format(self.ftau_avedensity_XAVE, self.ftau_avedensity_XSTD)
log.write_log(' ftau[AIC_w] = {:12f} +/- {:8f}'.format(self.ftau_density_XAVE,
self.ftau_density_XSTD))
#log.write_log(' ftau[ave AIC_w] = {:12f} +/- {:8f}'.format(self.ftau_avedensity_XAVE, self.ftau_avedensity_XSTD))

self.FTAU_density_XAVE = self.ftau_density_XAVE * 0.5 * self.tau_scale
self.FTAU_density_XSTD = self.ftau_density_XSTD * 0.5 * self.tau_scale
#self.FTAU_avedensity = self.ftau_avedensity*0.5*self.tau_scale
#self.FTAU_avedensity_XAVE = self.ftau_avedensity_XAVE*0.5*self.tau_scale
#self.FTAU_avedensity_XSTD = self.ftau_avedensity_XSTD*0.5*self.tau_scale
log.write_log(' FTAU[AIC_w] = {:12f} +/- {:8f}'.format(self.FTAU_density_XAVE, self.FTAU_density_XSTD))
#log.write_log ' FTAU[ave AIC_w] = {:12f} +/- {:8f}'.format(self.FTAU_avedensity_XAVE, self.FTAU_avedensity_XSTD)
log.write_log(' FTAU[AIC_w] = {:12f} +/- {:8f}'.format(self.FTAU_density_XAVE,
self.FTAU_density_XSTD))
#log.write_log(' FTAU[ave AIC_w] = {:12f} +/- {:8f}'.format(self.FTAU_avedensity_XAVE, self.FTAU_avedensity_XSTD))
log.write_log()
return

Expand Down
24 changes: 13 additions & 11 deletions thermocepstrum/utils/blocks.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,13 @@
#!/usr/bin/env python

from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()
import sys
import os
import math
import os.path
import thermocepstrum as tc

abs_path = os.path.abspath(sys.argv[0])
tc_path = abs_path[:abs_path.rfind('/')]
tc_path = tc_path[:tc_path.rfind('/')]
Expand All @@ -24,16 +30,11 @@
#plt.rc('text', usetex=True)
c = plt.rcParams['axes.prop_cycle'].by_key()['color']
from matplotlib.ticker import MultipleLocator
import math
import os.path
import thermocepstrum as tc

try:
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
except:
log.write_log('Error: cannot import inset_axes (will not be able to plot some parts of the plots)')
from thermocepstrum.utils.utils import PrintMethod
log = PrintMethod()


def main():
Expand Down Expand Up @@ -91,7 +92,7 @@ def main():
periodograms.append(np.loadtxt(fname + '.psd.dat', usecols=(3, 4), unpack=True))
if periodograms[-1].shape != periodograms[0].shape:
log.write_log(fname, ' not used (inconsistent shape with firts element)', periodograms[-1].shape,
periodograms[0].shape)
periodograms[0].shape)
del periodograms[-1]
continue
if os.path.isfile(fname + '.cepstral.npy'):
Expand Down Expand Up @@ -131,7 +132,7 @@ def main():
cepstrals[i].resize(cepstrals[0].shape)
# *this does not work when using a lot of data
#import pdb; pdb.set_trace()
#log.write_log periodograms[0].shape
#log.write_log(periodograms[0].shape)
#for i in range(1,len(periodograms)):
# log.write_log i
# periodograms[i]=np.resize(periodograms[i],periodograms[0].shape)
Expand Down Expand Up @@ -199,8 +200,9 @@ def main():
for i in range(all_cut):
if mean_periodogram[0, i] > zero:
selection_not_zero.append(i)
log.write_log('Number of components > {}: {}. Last is {}'.format(zero, len(selection_not_zero), selection_not_zero[-1]))
# log.write_log selection_not_zero
log.write_log('Number of components > {}: {}. Last is {}'.format(zero, len(selection_not_zero),
selection_not_zero[-1]))
# log.write_log(selection_not_zero)

log.write_log('Index = {} , {} THz'.format(p95, freqs[p95]))

Expand Down Expand Up @@ -375,7 +377,7 @@ def ffpsd(self, w, single=False):
plt.close()

# np.savetxt(output+'.histogram_all',np.c_[(intervals[1:]+intervals[:-1])/2.0,histogram/np.sum(histogram)])
# log.write_log 'Histogram bin width: {}'.format(intervals[1]-intervals[0])
# log.write_log('Histogram bin width: {}'.format(intervals[1]-intervals[0]))

# np.savetxt(output+'.kolmogorov_smirnov',[ks_0,ks_1,ks_all])
log.write_log('Statistical test results (psd(0), psd(1), psd(all but 0)): {}'.format([ks__0, ks__1, ks_all]))
Expand Down Expand Up @@ -536,7 +538,7 @@ def n_tick_in_range(beg, end, n, n_c=1, nit=0):
if cifre == 0:
cifre = 1.0
delta = cifre * e / 10**(n_c)
#log.write_log "n=",n, " dx0=",dx0," e=",e," m=" ,m," cifre=",cifre
#log.write_log("n=",n, " dx0=",dx0," e=",e," m=" ,m," cifre=", cifre)
if nit < 30:
if delta >= size:
return n_tick_in_range(beg, end, n + 1, n_c, nit + 1)
Expand Down

0 comments on commit 6f76ef6

Please sign in to comment.