• Update GUI to work with the new core API.
• Add functionalities

LAMMPS input files

With reference to #18, LAMMPS input files would be useful.

• SiO2
• NaCl

NaCl

• The physical parameters of the system (time step, volume, number of atoms, potential, LAMMPS input file) are not reported in NaCl.dat example data.

Ok, I have a huge screen. In this screenshot the plot should occupy most of the space

After loading a saved .npy file, the volume and the timestep do not show up

in this folder there are 3 files: test.dat, test.dat.dat, test.dat.dat.dat, but only the first is displayed

first I put filter width 3THz, and I click 'Resample':

then, after clicking 'Resample'

What happened?

Then, both plots looks like they do not have the filter applied

Using matplotlib>=2.2 cause a warning to be issued when using the analysis program.

/home/lercole/.virtualenvs/thermotest2/local/lib/python2.7/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)

or

/home/lercole/.virtualenvs/thermotest3/lib/python3.6/site-packages/matplotlib/figure.py:98: MatplotlibDeprecationWarning: 
Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  "Adding an axes using the same arguments as a previous axes "

For matplotlib>=3.1 it becomes:

/home/lercole/.virtualenvs/thermotest3/lib/python3.6/site-packages/thermocepstrum/analysis.py:666: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.axes().xaxis.set_major_locator(MultipleLocator(dx1))
/home/lercole/.virtualenvs/thermotest3/lib/python3.6/site-packages/thermocepstrum/analysis.py:667: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.axes().xaxis.set_minor_locator(MultipleLocator(dx2))
/home/lercole/.virtualenvs/thermotest3/lib/python3.6/site-packages/thermocepstrum/analysis.py:668: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.axes().yaxis.set_major_locator(MultipleLocator(dy1))
/home/lercole/.virtualenvs/thermotest3/lib/python3.6/site-packages/thermocepstrum/analysis.py:669: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.

The problem seems to refer to these lines:
https://github.com/lorisercole/thermocepstrum/blob/4e89421c45d3eb95d445f1ed5d7f2e5e25532f22/thermocepstrum/analysis.py#L671-L674

I noticed that cepstral_log prints with float string format, which can be annoying when dealing with very big/small numbers:

In the examples of #62, this is the output, which is kind of hard to read

-----------------------------------------------------
  CEPSTRAL ANALYSIS
-----------------------------------------------------
  cutoffK = (P*-1) = 10  (auto, AIC_Kmin = 10, corr_factor =  1.0)
  L_0*   =          18.297017 +/-   0.057597
  S_0*   =   127801909.034038 +/- 7361041.782368
-----------------------------------------------------
  kappa* =           0.000229 +/-   0.000013  Pa*s
-----------------------------------------------------

I propose to use scientific notation whenever numbers are too big/small.

Hi, I want to spend two words on a possible source of problems in the multivariate cepstral analysis.
Using the multivariate methods, if you (wrongly) decorrelate from two linearly dependent currents, in some cases, it is possible that the code does not issue an error, and returns a result that seems reasonable, even though it is, of course, wrong.
The reason is that if the two currents are linearly dependent theoretically, but not exactly numerically, the resulting determinant of the covariate matrix is small but not exactly zero.

Working on branch fix_27_ci_setup, that was branched off develop, as this was the last stable version.

CI

• Setup GitHub Action CI workflow

Tests

• We should implement the dreaded fixtures to read the data files and do a better job.
  @rikigigi do you know how to do the magic? :)
• Tests should be extended to test more internal values, e.g. psd, fpsd, logpsd, cepstral coefficients, ...
• Tests should be extended to test the CLI and the GUI

Once these tests produce the expected results, we can merge this branch into core_redesign, and use them to test the new API.

Is there a summary of what needs to be done when switching from thermocepstrum-analysis to sportran-analysis? I had used thermocepstrum for a few years, and just upgraded machines and hence had to reinstall. I already noticed that I had to provide two new flags, --VOLUME and --TEMPERATURE, which appeared to require me to compute those quantities myself (they're columns in the input table, but before whatever processing was needed was getting done automatically). I'm getting strange
errors (which I'm trying to debug, and may open other issues), and I'm worried there's something else I need to change.

[edit] looks like --TEMPERATURE is supposed to be getting detected automatically (there's a table column with the correct name), and for some reason it isn't. I'll try to debug that. But if the answer to my original question is that there's not supposed to be any change required, then feel free to say so and close this issue. If there are changes, then documenting them would still be helpful.

If I try to run sportran-analysis without a -w flag, I get the following error

> sportran-analysis --param-from-input-file-column Temp TEMPERATURE --param-from-input-file-column Volume VOLUME -k flux -t 5 --input-format table -u metal -r -c 2 --TSKIP 2 --output TCEPS.out -- tmp2.txt
.
.
.
Traceback (most recent call last):
  File "/home/cluster2/bernstei/.local/bin/sportran-analysis", line 8, in <module>
    sys.exit(main())
  File "/home/cluster2/bernstei/.local/lib/python3.9/site-packages/sportran/analysis.py", line 210, in main
    run_analysis(args)
  File "/home/cluster2/bernstei/.local/lib/python3.9/site-packages/sportran/analysis.py", line 476, in run_analysis
    outarray = np.c_[j.freqs_THz, j.psd, j.fpsd, j.logpsd, j.flogpsd]
  File "/home/Software/python/system/extra/lib64/python3.9/site-packages/numpy/lib/index_tricks.py", line 412, in __getitem__
    res = self.concatenate(tuple(objs), axis=axis)
  File "<__array_function__ internals>", line 180, in concatenate
ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 5 and the array at index 2 has size 1

From what I can tell the issue is that two of the arrays being passed to np.c_, j.fpsd and j.flogpsd, are None, so the concatenation fails.

If I add a -w 1, which sets a filter width and therefore creates j.fpsd, it doesn't fail here. Note that with this example input it does fail later, but that's just because my example file is very short. With a reasonable length input it succeeds if and only if I pass a -w argument.

The example input file is attached.

tmp2.txt

https://github.com/lorisercole/thermocepstrum/blob/8a685b347947605c5d9016c2367db83541bdf64d/thermocepstrum/plotter/plotter.py#L487-L493

this line for example should stay in the GUI part

maybe some refactoring is needed

Copied from old wishlist.

THINGS TO FIX

schermata select a file

• BUG: se clicco sul link "README.md" nella finestra di Help mi dà un errore perché non trova il file:
  file:///home/lercole/software/README.md: Error when getting information for file “/home/lercole/software/README.md”: No such file or directory

• il cambiamento delle settings (lingua o font) avviene solo al riavvio del programma

• La schermata di selezione dei file confonde un po' le idee: sopra vengono elencati i file presenti nella cartella corrente, sotto c'è il preview del file selezionato, e in fondo nella casella "selected" è scritto tutto il percorso del file aperto. Il nome di questo file è spesso illeggibile perché troppo lungo.
  Per rendere tutto più lineare (dall'alto in basso) metterei:

• La casella per la selezione del file all'inizio.

• Non è chiaro perché ci sia una lista di files della cartella corrente, se il file si seleziona tramite menù nella casella precedente.

• Sarebbe bello che ci fosse la possibilità di navigare nelle cartelle direttamente in questa schermata (senza aprire una finestrella - tipo il Finder del Mac), e così si potrebbbe selezionare il file dalla lista di file nella cartella corrente.

• Se non si può: o eliminiamo la lista dei file, oppure la teniamo ma scriviamo qual è la cartella corrente ed evidenziamo il file selezionato.

• Vanno chiariti i tipi di file supportati... Ci dovremmo mettere un link al manuale (quando ci sarà un manuale), ma intanto si può scrivere una descrizione breve dei formati lì.

• C'è un modo per far sì che se chiudo e riapro il programma vada ad aprire direttamente la cartella dell'ultimo file usato, senza ripartire sempre dalla Home?

schermata define headers

• BUG: Se uso il file binario data/Silica.npy, setto Temperature: None, clicco next, poi back: Temperature sparisce :)

• BUG: il programma permette di settare più colonne Temperature, Volume, o DT_FS

• Vanno chiarite le unità di misura.

schermata set variable

• permetterei all'utente di cambiare i valori delle variabili anche se queste vengono lette da file. Il valore che compare dovrebbe essere quello calcolato

schermata $f^*$

• Manca un titolo/spiegazione.

• Il margine sx del grafico è troppo piccolo, i numeri sull'asse y mi vengono tagliati fuori.

• a cosa serve il pulsante bianco quadrato a sinistra di Reset View?

• farei lo sfondo azzurro della regione selezionata dallo slider un po' più chiaro, oppure giallino chiaro

• il sistema zoom-in reset-view può essere un po' macchinoso. Manca la possibilità di fare zoom-in più volte e zoom-out.

• BUG: se seleziono una zona e clicco zoom-in, poi clicco Back le variabili tornano a zero (se il file è di tipo Table). Oltre a questo, se ora imposto le variabili (o tengo quelle già presenti) e vado avanti, $f^$ rimane bloccata su zero, e lo slider non funziona. L'unico modo per sistemarlo è impostare manualmente un nuovo $f^$ nella casella.

• Se uno decide di utilizzare i bottoni del grafico per zommare/pannare, appena tocca lo slider lo zoom torna a quello di default. Questo comportamento potrebbe essere fastidioso... non so se si può risolvere (ad esempio forse si possono ottenere i nuovi limiti del grafico con plt.xlim()).

schermata $P^*$

• manca titolo/spiegazione

• cos'è la linea rossa a $f=0$?

• IMPORTANT: $P^$ a differenza di $f^$ non è un parametro che l'utente sceglie, di default. Esso viene calcolato automaticamente (è aic_Kmin + 1, mi sembra) tramite l'Akaike Information Criterion.

• Scriverei quindi chiaramente qual è il $P^*$ suggerito dall'AIC.

• Inserirei un bottone per resettarlo a questo valore.

• Se l'utente vuole modificarlo glielo lasciamo fare, ma chiamerei questa casella "$P^*$ correction"

• prevederei la possibilità di modificare il filtro del plot (Filter width) anche qui, così che uno possa modificarlo per analizzare meglio i risultati.

GENERAL WISH LIST

• Vanno inserite delle spiegazioni che aiutino a capire cosa rappresenta ogni schermata, senza aver letto approfonditamente il paper. Ad esempio spiegare come scegliere $f^*$.

• L'help per come è adesso è difficile da leggere. Avrebbe più senso preparare un help html più dettagliato, con degli screenshot di esempio. Tra le cose che vorrei implementare per la versione 1.0 di st c'è la documentazione su readthedocs, fatta con sphinx.

• si potrebbe fare che se il file letto ha estensione .npy o .npz viene riconosciuto come binario (dict) direttamente.

• Per tutti i plot: inserire un pulsante per visualizzare il log-periodogramma al posto del periodogramma.

• $P^*$ plot: inserire un pulsante per visualizzare i seguenti plot (in finestre esterne):

• coefficienti cepstrali (HeatCurrent.plot_ck())

• analisi di convergenza kappa(P*) (HeatCurrent.plot_kappa_Pstar() e/o HeatCurrent.plot_L0_Pstar()

• analisi di $f^$ -- più difficile, perché richiede di specificare una lista di $f^$ o TSKIP e può essere lenta. Dovremmo aggiungere delle funzioni che generino delle griglie lineari/logaritmiche, così da rendere la selezione più semplice. Per ora concentriamoci sulle altre cose...

packaging

• get rid of uncertainties package

• get rid of tk-html-widgets if possible

• ensure that the code is Python 2 compatible

if aic_Kmin=0 choosing different P* with the multiplication of Kmin_corrfactor is not possible. Moreover, the function plot_kappa_Pstar() brakes because it chooses the range of P* as 3 * self.dct.aic_Kmin e.g. line 170 of plotter.py

max_y = np.amax(self.kappa_scale * 0.5 *
                            (self.dct.tau + self.dct.tau_THEORY_std)[self.dct.aic_Kmin:3 * self.dct.aic_Kmin])

We are going to rename this repository to SporTran.
The new repository will be https://github.com/lorisercole/sportran. GitHub will redirect automatically the old address to the new one.
All issues, PRs, comments, commits, followers, settings, etc will be moved automatically.

• change all references in the code from thermocepstrum to sportran
• change all references to the repository from https://github.com/lorisercole/thermocepstrum to https://github.com/lorisercole/sportran
• add reference to the new article (in preparation)
• add notice in the README and in the repository's About
• update/add license header to all files where applicable

https://github.com/lorisercole/thermocepstrum/blob/8a685b347947605c5d9016c2367db83541bdf64d/thermocepstrum/plotter/plotter.py#L443

then after loading a new file, the messages are wrong

wrong reading of the temperature when opening a table file from scratch: the temperature field is editable, but if it is edited, it is silently ignored

The f* analysis is not well documented.
We should add it to the examples, and to the future documentation.

I'm trying to run sportran-analysis with a table that has a column for temperature, one for volume, and 48 for flux components (xyz x 16 independent runs). I get the following error

> sportran-analysis --VOLUME 392737.783943168 -k flux -t 5 --input-format table -u metal -r --FSTAR 1 -c 2 --output TCEPS.out -- tmp2.txt
.
.
Traceback (most recent call last):
  File "/home/cluster2/bernstei/.local/bin/sportran-analysis", line 8, in <module>
    sys.exit(main())
  File "/home/cluster2/bernstei/.local/lib/python3.9/site-packages/sportran/analysis.py", line 210, in main
    run_analysis(args)
  File "/home/cluster2/bernstei/.local/lib/python3.9/site-packages/sportran/analysis.py", line 393, in run_analysis
    j = st.current.all_currents[current_type][0](currents, DT_FS=DT_FS, UNITS=units, **parameters,
  File "/home/cluster2/bernstei/.local/lib/python3.9/site-packages/sportran/current/current.py", line 58, in __init__
    raise ValueError('The input parameters {} must be defined.'.format(self._input_parameters - keyset))
ValueError: The input parameters {'TEMPERATURE'} must be defined.

I see where the "Temp" field is detected (analysis.py line 317), but I don't see where value is supposed to be inserted into the parameters dict before it's passed in line 393.

tmp2.txt

With the term "multi-current" we refer to the a system described by N_CURRENTS currents associated to conserved quantities (e.g. energy, mass, charge, ...). For example a 2-component liquid will have 2 currents (energy and one of the particle currents). In the code, this case is signaled by the flag MANY_CURRENTS = True.

To avoid confusion, we dub as "multi-component" the case where we consider N_COMPONENTS (Cartesian) components for each current, e.g. the 3 Cartesian components of the energy current. This is flagged by MANY_COMPONENTS = True.
(see 385c0bf)

Issues to be addressed

• move the method MDSample.compute_kappa_multi to Current.
  https://github.com/lorisercole/thermocepstrum/blob/7a7650626a99f8fd354af31aef706ca1343db7d8/thermocepstrum/md/mdsample.py#L358-L359
• implement a function to check that the N_COMPONENTS is large enough wrt N_CURRENTS to comply with the constraints of the multi-current theory, for the system to be solvable.
• check that all the multi-current routines work properly.

current-specific units

Each Current subclass, such as HeatCurrent, can define its own units by placing them in the folder current/units in a file with the same name as the current's _current_type attribute (or the current module's file name) (e.g. current/units/heat.py).

We shall:

• find a way to automatically discover the units defined for a specific current
• define a class method to return the class' specific list of units

This machinery should be easily extendable by a user. For example, one would like to define new units by simply defining a new function/class in current/units/heat.py.
Perhaps a way to do this is by using a decorator, like @add_units to the function defining a type of units, similarly to what is done in AiiDA by the @process_handler decorator.

generic-current units

One could in principle also define a generic Current object (and not a specific subclass).
In this case the initialization would take one input parameter: the "kappa_scale", that is the constant that multiplies the Green-Kubo integral.

• implement this generic initialization for the base Current class.

Following the paper reviewers' comments:

• write a base template class with all the mandatory methods of a current raising NotImplementedError.
  This will insure that all the required methods would have to be implemented and not picked up from current.Current based class.

A new plotter module is needed.
This module should define all the needed plot functions to be called by the other code's modules.

An idea would be to define a general interface in a Plotter class. Each plot type will become a subclass of it.
We need to figure out a common, easy, and extendable way of passing the parameters. The user should be able to call the plotting methods and to add customizations to the plots.
We aim to have all the core, CLI, and GUI modules calling the same plotting methods. This shall avoid duplication of the same plotting routines in the 3 modules.

The test_cli.py pytest fails with the package pytest-regressions==2.3.0.
It works with version 2.2.0.
The error is

FAILED tests/test_cli.py::test_cli_NaCl - TypeError: unsupported operand type(s) for -: 'str' and 'str'
FAILED tests/test_cli.py::test_cli_bin_output_NaCl - TypeError: unsupported operand type(s) for -: 'str' and 'str'

Notice also the warnings appearing during tests:

================================================================================================ warnings summary ================================================================================================
tests/test_as_example.py: 22 warnings
tests/test_cli.py: 14 warnings
  /home/loris/.virtualenvs/pydevelop/lib/python3.8/site-packages/pytest_regressions/dataframe_regression.py:126: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
  Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
    if data_type in [float, np.float, np.float16, np.float32, np.float64]:

tests/test_notebooks.py: 62 warnings
  /home/loris/.virtualenvs/pydevelop/lib/python3.8/site-packages/nbformat/validator.py:315: DeprecationWarning: Passing a schema to Validator.iter_errors is deprecated and will be removed in a future release. Call validator.evolve(schema=new_schema).iter_errors(...) instead.
    for error in errors:

-- Docs: https://docs.pytest.org/en/stable/warnings.html
======================================================================================== 6 passed, 98 warnings in 26.52s =========================================================================================

https://www.python.org/download/mac/tcltk/

If you are using macOS 10.6 or later, the Apple-supplied Tcl/Tk 8.5 has serious bugs that can cause application crashes. If you wish to use IDLE or Tkinter, do not use the Apple-supplied Pythons. Instead, install and use a newer version of Python from python.org or a third-party distributor that supplies or links with a newer version of Tcl/Tk.

I think it would be useful for users who are trying to use this package for their own input to have the LAMMPS input files corresponding to the examples provided here. I'm trying to do exactly that, and getting terrible disagreement between my conventional autocorrelation integration analysis and the cepstrum-based output. I don't know if I'm somehow generating bad data, and if I had the LAMMPS structures and potentials (and things like the time step used for the MD), I could check if my problem is somehow inherent to my system, or merely the way I'm setting up the heat flux equilibration and sampling, for example.

The github CI image ubuntu-latest does not have 3.6 python version. We can do one of the following

• remove 3.6 (and add 3.10+ python version)
• modify the github CI workflow

I found some cases where AIC does not find a minimum, or it finds it at a P* that is clearly too high.
For example KCl, where the spectrum is essentially one peak at zero-frequency.
I suspect that this happens when dealing with very long trajectories.
The cepstral coefficients (dct.logpsdK) go to zero very quickly, but you can find a few high n quefrencies that are outliers for no reason.

This in turn causes the AIC (that has a minimum at low n, where it should be) to drop discontinuously at these points, therefore marking a new spurious global minimum.

We should figure out where this behavior comes from: is it a Fourier transform problem, a problem of the data, a problem with the filter&resample, ...?

The file list does not expand

• The units of kappa given in output should depend on the type of current (generic, heat, electric, etc...)
• physical constants should be moved to current/units/constants.py and loaded where needed
  @fgrassel

but I selected 'qepw' before saving it

We should add to the examples the calculation of the autocorrelation function and the Green-Kubo integral, for comparisons against the cepstral method.
This is useful to check the consistency of units.

• ElectricCurrent @fgrassel
• viscosity
• ...

As in the title I tried
pip install . with pip 22.3.1 and got the following error:

  Installing build dependencies ... done
  Getting requirements to build wheel ... error
  error: subprocess-exited-with-error
  
  × Getting requirements to build wheel did not run successfully.
  │ exit code: 1
  ╰─> [5 lines of output]
      /tmp/pip-build-env-22g7g9sb/overlay/lib/python3.9/site-packages/setuptools/_distutils/dist.py:265: UserWarning: Unknown distribution option: 'gui_version'
        warnings.warn(msg)
      /tmp/pip-build-env-22g7g9sb/overlay/lib/python3.9/site-packages/setuptools/_distutils/dist.py:265: UserWarning: Unknown distribution option: 'credits'
        warnings.warn(msg)
      error in sportran setup command: 'python_requires' must be a string containing valid version specifiers; Invalid specifier: '>=3.6.*'
      [end of output]
  

I solved it changing in the setup.json: '>=3.6.*' to '>=3.6'

@lorisercole
Right now, the default binary output is a pickle dumped blob that, for a first time user, I think it is difficult to understand. Its content is:

['KAPPA_SCALE',
 'TEMPERATURE',
 'TSKIP',
 'UNITS',
 'VOLUME',
 '__class__',
 '__delattr__',
 '__dict__',
 '__dir__',
 '__doc__',
 '__eq__',
 '__format__',
 '__ge__',
 '__getattribute__',
 '__gt__',
 '__hash__',
 '__init__',
 '__init_subclass__',
 '__le__',
 '__lt__',
 '__module__',
 '__ne__',
 '__new__',
 '__reduce__',
 '__reduce_ex__',
 '__repr__',
 '__setattr__',
 '__sizeof__',
 '__str__',
 '__subclasshook__',
 '__weakref__',
 'cepstral_log',
 'j_DT_FS',
 'j_Nyquist_f_THz',
 'j_PSD_FILTER_W_THz',
 'j_cospectrum',
 'j_fcospectrum',
 'j_flogpsd',
 'j_fpsd',
 'j_freqs_THz',
 'j_logpsd',
 'j_psd',
 'jf_DT_FS',
 'jf_Nyquist_f_THz',
 'jf_dct_Kmin_corrfactor',
 'jf_dct_aic_Kmin',
 'jf_dct_kappa',
 'jf_dct_kappa_THEORY_std',
 'jf_dct_logpsd',
 'jf_dct_logpsdK',
 'jf_dct_logpsdK_THEORY_std',
 'jf_dct_logtau',
 'jf_dct_logtau_THEORY_std',
 'jf_dct_psd',
 'jf_flogpsd',
 'jf_fpsd',
 'jf_freqs_THz',
 'jf_logpsd',
 'jf_psd',
 'jf_resample_log',
 'kappa_Kmin',
 'kappa_Kmin_std',
 'units',
 'write_old_binary']

Is it used by anyone or anywhere in the code? Is it safe to change the default binary output to the one equivalent to the human readable one but with numpy arrays?

Following PEP 621 it is now possible to fully define the build procedure of your package in pyproject.toml.
The setup.py will be deprecated.
MANIFEST.in is no longer necessary if we use flit.

Once the Plotter module is ready (#25), we should remove all traces of the old plot functions. In particular from

• thermocepstrum/md/mdsample.py
• thermocepstrum/md/resample.py
• thermocepstrum/analysis.py

Related to #44
We want to allow the possibility of manually select P*, either with the Kmin_corrfactor or by providing the value explicitly.

I propose:

• rename CosFilter to CepstralFilter
• aic_Kmin should be the minimum of AIC, it cannot be overwritten,
• introduce cutoffK, that is the P*-1 used in the cepstral filter. By default it is = aic_Kmin. It can differ if aic_Kmin_corrfactor is used, or a manual_cutoffK_flag is specified.
• in scan_filter_tau: K_PSD renamed cutoffK, if specified it is used to compute tau_Kmin (renamed tau_cutoffK), and aic_Kmin is not overwritten
• update current.py and all the other modules.
• This change of definition for aic_Kmin must be reflected onto ALL PLOTS: cutoffK should be used instead of aic_Kmin
• update GUI

• #51
• Conform plots to the common Plotting format #43
• Update docs
• Add support for binary format (see #37)

There is interest in performing (thermal) conductivity calculations for anisotropic systems.
In general, one wants to compute the full 3x3 conductivity matrix:

\begin{bmatrix} \kappa_{xx} & \kappa_{xy} & \kappa_{xz} \\ \kappa_{yx} & \kappa_{yy} & \kappa_{yz} \\ \kappa_{zx} & \kappa_{zy} & \kappa_{zz} \\ \end{bmatrix}
where
\kappa_{\alpha\beta} = \frac{1}{V k_B T} \int_0^\infty \langle J_{\alpha} (t) J_{\beta}(0) \rangle dt

Each Cartesian component (heat)-current component J_α (α=1,2,3) is not equivalent to the others, therefore one should define a HeatCurrent with one Cartesian component only.

Diagonal components

The diagonal components are readily estimated by defining a HeatCurrent for each Cartesian component, and performing cepstral analysis:

# Let's assume that a heat-current is contained in a (N, 3) numpy array (that here I define randomly for demonstration):
current = np.random.rand(N, 3)

# one should then create a HeatCurrent for each component, as they were different stochastic processes:
j1 = tc.HeatCurrent(current[:, 0], units, DT_FS, TEMPERATURE, VOLUME)
j2 = tc.HeatCurrent(current[:, 1], units, DT_FS, TEMPERATURE, VOLUME)
j3 = tc.HeatCurrent(current[:, 2], units, DT_FS, TEMPERATURE, VOLUME)
# and perform all the normal operations and cepstral analysis for each of them
j1f, ax1 = tc.heatcurrent.resample_current(j, fstar_THz=FSTAR_THZ, plot=True)
j2f, ax2 = tc.heatcurrent.resample_current(j, fstar_THz=FSTAR_THZ, plot=True)
j3f, ax3 = tc.heatcurrent.resample_current(j, fstar_THz=FSTAR_THZ, plot=True)
# ...
j1f.cepstral_analysis()
j2f.cepstral_analysis()
j3f.cepstral_analysis()

Off-diagonal components

Off-diagonal components are a more difficult subject. Similarly to the multi-component-fluid (aka multi-current) case, the distribution of the spectrum will not be a χ², but a Wishart distribution, and we have to deal with covariance matrices. This requires an extension of the theory...

Roadmap

• Design an interface to streamline the calculation of the diagonal components of the conductivity matrix for the typical 3D anisotropic case. This interface just defines 3 HeatCurrent objects and offers its main methods.
• Extend the theory to include off-diagonal components. The one-current case should be the easiest to deal with. The multi-current case might be quite convoluted.