fixes #66 (#67)

* fixes #66 - there was a bug in `sportran-analysis` when the moving average window was not specified. * added test for cli when the user don't ask for the filtering of the data with the moving average approach * version updated to 1.0rc3 Co-authored-by: Loris Ercole <[email protected]>
sissaschool · Dec 7, 2022 · bc559d0 · bc559d0
1 parent 41005be
commit bc559d0
Show file tree

Hide file tree

Showing 12 changed files with 100,202 additions and 18 deletions.
diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS
@@ -1,4 +1,4 @@
 # The following owners will be the default owners for everything
 # in the repo unless a later match takes precedence.
-sportran/*      @lorisercole
-sportran_gui/*  @rikigigi
+sportran/      @lorisercole
+sportran_gui/  @rikigigi
diff --git a/MANIFEST.in b/MANIFEST.in
@@ -3,9 +3,9 @@ include LICENSE.txt
 include setup.json
 include sportran/README.md
 include sportran/metadata.json
-include sportran/utils/plot_style.mplstyle
-exclude sportran/utils/blocks.py
-exclude sportran/utils/blockanalysis.py
+exclude sportran/utils/obsolete/
 include sportran_gui/README_GUI.md
 include sportran_gui/assets/icon.gif
 include sportran_gui/assets/languages.json
+include sportran/plotter/styles/api_style.mplstyle
+include sportran/plotter/styles/cli_style.mplstyle
diff --git a/README.md b/README.md
@@ -9,7 +9,7 @@ A code to estimate transport coefficients from the cepstral analysis of a multi-
 https://sportran.readthedocs.io
 
 ### References
- - [Ercole L., Bertossa R., Bisacchi S., and Baroni S., "_SporTran: a code to estimate transport coefficients from the cepstral analysis of (multivariate) current time series_", *arXiv*:2202.11571 (2022)](https://arxiv.org/abs/2202.11571), submitted to *Comput. Phys. Commun.*
+ - [Ercole L., Bertossa R., Bisacchi S., and Baroni S., "_SporTran: a code to estimate transport coefficients from the cepstral analysis of (multivariate) current time series_", *Comput. Phys. Commun.*, 108470](https://doi.org/10.1016/j.cpc.2022.108470), [*arXiv*:2202.11571 (2022)](https://arxiv.org/abs/2202.11571)
  - (cepstral analysis) [Ercole, Marcolongo, Baroni, *Sci. Rep.* **7**, 15835 (2017)](https://doi.org/10.1038/s41598-017-15843-2)
  - (multicomponent systems) [Bertossa, Grasselli, Ercole, Baroni, *Phys. Rev. Lett.* **122**, 255901 (2019)](https://doi.org/10.1103/PhysRevLett.122.255901) ([arXiv](https://arxiv.org/abs/1808.03341))
  - (review) [Baroni, Bertossa, Ercole, Grasselli, Marcolongo, *Handbook of Materials Modeling* (2018)](https://doi.org/10.1007/978-3-319-50257-1_12-1) ([arXiv](https://arxiv.org/abs/1802.08006))

diff --git a/setup.json b/setup.json
@@ -1,6 +1,6 @@
 {
     "name": "sportran",
-    "version": "1.0.0rc2",
+    "version": "1.0.0rc3",
     "author": "Loris Ercole, Riccardo Bertossa, Sebastiano Bisacchi",
     "author_email": "[email protected]",
     "description": "Cepstral Data Analysis of current time series for Green-Kubo transport coefficients",

diff --git a/sportran/analysis.py b/sportran/analysis.py
@@ -211,6 +211,18 @@ def main():
     return 0
 
 
+def concatenate_if_not_none_with_labels(concat, labels=None):
+    out_arr = []
+    out_label = ''
+    if labels is None:
+        labels = ['' for i in concat]
+    for arr, label in zip(concat, labels):
+        if arr is not None:
+            out_arr.append(arr)
+            out_label += f' {label}'
+    return np.concatenate([out_arr], axis=1).transpose(), f'{out_label}\n'
+
+
 def run_analysis(args):
 
     inputfile = args.inputfile
@@ -472,8 +484,8 @@ def average(data, name, units=''):
 
     if not no_text_out:
         outfile_name = output + '.psd.dat'
-        outarray = np.c_[j.freqs_THz, j.psd, j.fpsd, j.logpsd, j.flogpsd]
-        outfile_header = 'freqs_THz  psd  fpsd  logpsd  flogpsd\n'
+        outarray, outfile_header = concatenate_if_not_none_with_labels(
+            [j.freqs_THz, j.psd, j.fpsd, j.logpsd, j.flogpsd], ['freqs_THz', 'psd', 'fpsd', 'logpsd', 'flogpsd'])
         np.savetxt(outfile_name, outarray, header=outfile_header, fmt=fmt)
         if j.MANY_CURRENTS:
             outfile_name = output + '.cospectrum.dat'
@@ -483,15 +495,17 @@ def average(data, name, units=''):
             np.savetxt(outfile_name, np.column_stack([outarray.real, outarray.imag]), fmt=fmt)
 
             outfile_name = output + '.cospectrum.filt.dat'
-            outarray = np.c_[j.freqs_THz,
-                             j.fcospectrum.reshape(
-                                 (j.fcospectrum.shape[0] * j.fcospectrum.shape[1], j.fcospectrum.shape[2])).transpose()]
-            np.savetxt(outfile_name, np.column_stack([outarray.real, outarray.imag]), fmt=fmt)
+            if j.fcospectrum is not None:
+                outarray = np.c_[j.freqs_THz,
+                                 j.fcospectrum.reshape((j.fcospectrum.shape[0] * j.fcospectrum.shape[1],
+                                                        j.fcospectrum.shape[2])).transpose()]
+                np.savetxt(outfile_name, np.column_stack([outarray.real, outarray.imag]), fmt=fmt)
 
         if resample:
             outfile_name = output + '.resampled_psd.dat'
-            outarray = np.c_[jf.freqs_THz, jf.psd, jf.fpsd, jf.logpsd, jf.flogpsd]
-            outfile_header = 'freqs_THz  psd  fpsd  logpsd  flogpsd\n'
+            outarray, outfile_header = concatenate_if_not_none_with_labels(
+                [jf.freqs_THz, jf.psd, jf.fpsd, jf.logpsd, jf.flogpsd],
+                ['freqs_THz', 'psd', 'fpsd', 'logpsd', 'flogpsd'])
             np.savetxt(outfile_name, outarray, header=outfile_header, fmt=fmt)
 
         outfile_name = output + '.cepstral.dat'
@@ -623,7 +637,8 @@ def write_old_binary(self, output):
         """Write old binary format."""
         opts = {'allow_pickle': False}
         optsa = {'axis': 1}
-        outarray = np.c_[self.j_freqs_THz, self.j_fpsd, self.j_flogpsd, self.j_psd, self.j_logpsd]
+        outarray, _ = concatenate_if_not_none_with_labels(
+            [self.j_freqs_THz, self.j_fpsd, self.j_flogpsd, self.j_psd, self.j_logpsd])
         np.save(output + '.psd.npy', outarray, **opts)
 
         if self.j_cospectrum is not None:
@@ -634,7 +649,8 @@ def write_old_binary(self, output):
             outarray = np.c_[self.j_freqs_THz, self.j_fcospectrum.reshape(-1, self.j_fcospectrum.shape[-1]).transpose()]
             np.save(output + '.cospectrum.filt.npy', outarray, **opts)
 
-        outarray = np.c_[self.jf_freqs_THz, self.jf_psd, self.jf_fpsd, self.jf_logpsd, self.jf_flogpsd]
+        outarray, _ = concatenate_if_not_none_with_labels(
+            [self.jf_freqs_THz, self.jf_psd, self.jf_fpsd, self.jf_logpsd, self.jf_flogpsd])
         np.save(output + '.resampled_psd.npy', outarray, **opts)
 
         outarray = np.c_[self.jf_cepf_logpsdK, self.jf_cepf_logpsdK_THEORY_std, self.jf_cepf_logtau,

diff --git a/sportran/md/mdsample.py b/sportran/md/mdsample.py
@@ -320,7 +320,7 @@ def filter_psd(self, PSD_FILTER_W=None, freq_units='THz', window_type='rectangul
         """
         Filter the periodogram with the given PSD_FILTER_W [freq_units].
           - PSD_FILTER_W  PSD filter window [freq_units]
-          - freq_units    frequency units   ['THz', 'red' (default)]
+          - freq_units    frequency units   ['THz' (default), 'red']
           - window_type   filtering window type ['rectangular']
         """
         if self.psd is None:

diff --git a/sportran/md/tools/filter.py b/sportran/md/tools/filter.py
@@ -13,5 +13,10 @@ def runavefilter(X, WF):
         WF = WF + 1
     W = int(WF / 2)
 
+    if W > X.shape[0] - 1:
+        W = X.shape[0] - 1
+        WF = 2 * W
+        print('Warning: reducing filtering window')
+
     Y = np.concatenate((X[W:0:-1], X, X[-2:-W - 2:-1]))
     return np.convolve(Y, np.array([1.0 / WF] * WF), 'valid')
diff --git a/tests/test_cli.py b/tests/test_cli.py
@@ -35,6 +35,27 @@ def test_cli_NaCl(tmpdir, run_cli, data_NaCl_path, num_regression, file_regressi
     num_regression.check(readed)
 
 
+def test_cli_NaCl_no_w(tmpdir, run_cli, data_NaCl_path, num_regression, file_regression):
+    fileout_pref = tmpdir.join('output_no_w')
+    output = run_cli(data_NaCl_path, '-k', 'flux', '-j', 'vcm[1]', '-t', '5.0', '--VOLUME', '65013.301261',
+                     '--param-from-input-file-column', 'Temp', 'TEMPERATURE', '--FSTAR', '14.0', '-r',
+                     '--test-suite-run', '-o', fileout_pref)
+    file_list = [
+        'output_no_w.cepstral.dat', 'output_no_w.cepstrumfiltered_psd.dat', 'output_no_w.cospectrum.dat',
+        'output_no_w.psd.dat', 'output_no_w.resampled_psd.dat'
+    ]
+    with open(str(tmpdir.join('output_no_w.log'))) as l:
+        file_regression.check(l.read(), basename='output_no_w.log')
+    file_regression.check(output.stdout.str(), basename='stdout_no_w')
+    file_regression.check(output.stderr.str(), basename='stderr_no_w')
+    readed = {}
+    for f in file_list:
+        file_out = tmpdir.join(f)
+        readed[f] = np.loadtxt(file_out).flatten()
+
+    num_regression.check(readed)
+
+
 def test_cli_bin_output_NaCl(tmpdir, run_cli, data_NaCl_path, num_regression, file_regression):
     fileout_pref = tmpdir.join('output')
     output = run_cli(data_NaCl_path, '-k', 'flux', '-j', 'vcm[1]', '-t', '5.0', '--VOLUME', '65013.301261',

diff --git a/tests/test_cli/output_no_w.log.txt b/tests/test_cli/output_no_w.log.txt
@@ -0,0 +1,66 @@
+ Units:      metal
+ Time step:      5.0 fs
+# Molten NaCl - Fumi-Tosi potential
+# https://journals.aps.org/prl/supplemental/10.1103/PhysRevLett.122.255901/Bertossa_et_al_Supplemental_Material.pdf
+# 1728 atoms, T~1400K
+# NVE, dt = 1.0 fs, 100 ps, print_step = 5.0 fs
+# Volume = 65013.301261 A^3
+# LAMMPS metal units
+Temp    c_flux[1] c_flux[2] c_flux[3] c_vcm[1][1] c_vcm[1][2] c_vcm[1][3]
+ #####################################
+  all_ckeys = [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))]
+ #####################################
+Data length = 20000
+  ckey = [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))]
+  ( 20000 ) steps read.
+VOLUME (input): 65013.301261
+Mean TEMPERATURE (computed): 1399.3477811999999 +/- 19.318785820942594
+ Time step (input):  5.0 fs
+  currents shape is (2, 20000, 3)
+snippet:
+[[[ 2.5086549e+02  2.0619423e+01  2.0011500e+02]
+  [ 1.9622265e+02  8.2667342e+01  2.8433250e+02]
+  [ 1.2639441e+02  1.6075472e+02  3.4036769e+02]
+  ...
+  [ 1.7991856e+02  1.8612706e+01 -1.3265623e+02]
+  [ 2.0471193e+02 -4.6643315e-01 -2.0401650e+02]
+  [ 2.4123318e+02 -1.8295461e+01 -2.5246475e+02]]
+
+ [[-1.5991832e-01 -7.1370426e-02  2.0687917e-02]
+  [-1.3755206e-01 -7.1002931e-02 -1.1279876e-02]
+  [-1.0615044e-01 -6.2381243e-02 -4.1568120e-02]
+  ...
+  [-9.1939899e-02 -8.4778292e-02  6.0011385e-02]
+  [-1.3384949e-01 -1.1474530e-01  8.9323546e-02]
+  [-1.8385053e-01 -1.3693430e-01  1.1434060e-01]]]
+Using multicomponent code.
+ Number of currents = 2
+ Number of equivalent components = 3
+ KAPPA_SCALE = 1.4604390788939313e-07
+ Nyquist_f   = 100.0  THz
+Using multicomponent code.
+-----------------------------------------------------
+  RESAMPLE TIME SERIES
+-----------------------------------------------------
+ Original Nyquist freq  f_Ny =     100.00000 THz
+ Resampling freq          f* =      14.28571 THz
+ Sampling time         TSKIP =             7 steps
+                             =        35.000 fs
+ Original  n. of frequencies =         10001
+ Resampled n. of frequencies =          1429
+ min(PSD)          (pre-filter&sample) =      0.31536
+ min(PSD)         (post-filter&sample) =  22166.11934
+ % of original PSD Power f<f* (pre-filter&sample)  = 96.679 %
+ fPSD not calculated before resampling
+ -----------------------------------------------------
+
+-----------------------------------------------------
+  CEPSTRAL ANALYSIS
+-----------------------------------------------------
+  cutoffK = (P*-1) = 3  (auto, AIC_Kmin = 3, corr_factor =  1.0)
+  L_0*   =          15.158757 +/-   0.056227
+  S_0*   =     6824108.702608 +/- 383697.095268
+-----------------------------------------------------
+  kappa* =           0.498310 +/-   0.028018  W/m/K
+-----------------------------------------------------
+
diff --git a/tests/test_cli/stderr_no_w.txt b/tests/test_cli/stderr_no_w.txt
diff --git a/tests/test_cli/stdout_no_w.txt b/tests/test_cli/stdout_no_w.txt
@@ -0,0 +1,65 @@
+ Units:      metal
+ Time step:      5.0 fs
+# Molten NaCl - Fumi-Tosi potential
+# https://journals.aps.org/prl/supplemental/10.1103/PhysRevLett.122.255901/Bertossa_et_al_Supplemental_Material.pdf
+# 1728 atoms, T~1400K
+# NVE, dt = 1.0 fs, 100 ps, print_step = 5.0 fs
+# Volume = 65013.301261 A^3
+# LAMMPS metal units
+Temp    c_flux[1] c_flux[2] c_flux[3] c_vcm[1][1] c_vcm[1][2] c_vcm[1][3]
+ #####################################
+  all_ckeys =  [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))]
+ #####################################
+Data length =  20000
+  ckey =  [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))]
+  ( 20000 ) steps read.
+VOLUME (input): 65013.301261
+Mean TEMPERATURE (computed): 1399.3477811999999 +/- 19.318785820942594
+ Time step (input):  5.0 fs
+  currents shape is (2, 20000, 3)
+snippet:
+[[[ 2.5086549e+02  2.0619423e+01  2.0011500e+02]
+  [ 1.9622265e+02  8.2667342e+01  2.8433250e+02]
+  [ 1.2639441e+02  1.6075472e+02  3.4036769e+02]
+  ...
+  [ 1.7991856e+02  1.8612706e+01 -1.3265623e+02]
+  [ 2.0471193e+02 -4.6643315e-01 -2.0401650e+02]
+  [ 2.4123318e+02 -1.8295461e+01 -2.5246475e+02]]
+
+ [[-1.5991832e-01 -7.1370426e-02  2.0687917e-02]
+  [-1.3755206e-01 -7.1002931e-02 -1.1279876e-02]
+  [-1.0615044e-01 -6.2381243e-02 -4.1568120e-02]
+  ...
+  [-9.1939899e-02 -8.4778292e-02  6.0011385e-02]
+  [-1.3384949e-01 -1.1474530e-01  8.9323546e-02]
+  [-1.8385053e-01 -1.3693430e-01  1.1434060e-01]]]
+Using multicomponent code.
+ Number of currents = 2
+ Number of equivalent components = 3
+ KAPPA_SCALE = 1.4604390788939313e-07
+ Nyquist_f   = 100.0  THz
+Using multicomponent code.
+-----------------------------------------------------
+  RESAMPLE TIME SERIES
+-----------------------------------------------------
+ Original Nyquist freq  f_Ny =     100.00000 THz
+ Resampling freq          f* =      14.28571 THz
+ Sampling time         TSKIP =             7 steps
+                             =        35.000 fs
+ Original  n. of frequencies =         10001
+ Resampled n. of frequencies =          1429
+ min(PSD)          (pre-filter&sample) =      0.31536
+ min(PSD)         (post-filter&sample) =  22166.11934
+ % of original PSD Power f<f* (pre-filter&sample)  = 96.679 %
+ fPSD not calculated before resampling
+ -----------------------------------------------------
+
+-----------------------------------------------------
+  CEPSTRAL ANALYSIS
+-----------------------------------------------------
+  cutoffK = (P*-1) = 3  (auto, AIC_Kmin = 3, corr_factor =  1.0)
+  L_0*   =          15.158757 +/-   0.056227
+  S_0*   =     6824108.702608 +/- 383697.095268
+-----------------------------------------------------
+  kappa* =           0.498310 +/-   0.028018  W/m/K
+-----------------------------------------------------