diff --git a/CMakeLists.txt b/CMakeLists.txt index 0528aacbe..ffba0d1b6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features project( EXP - VERSION "7.8.0" + VERSION "7.8.1" HOMEPAGE_URL https://github.com/EXP-code/EXP LANGUAGES C CXX Fortran) diff --git a/Paper/paper/paper.jats b/Paper/paper/paper.jats index 194915004..ca87725a4 100644 --- a/Paper/paper/paper.jats +++ b/Paper/paper/paper.jats @@ -82,33 +82,34 @@ a Creative Commons Attribution 4.0 International License (CC BY Summary -

Galaxies are nearly equilibrium ensembles of dark and baryonic - matter confined by their mutual gravitational attraction. Over many - decades of rich development, the dynamics of galactic evolution have - been studied with a combination of numerical simulations and - analytical methods +

Galaxies are ensembles of dark and baryonic matter confined by + their mutual gravitational attraction. The dynamics of galactic + evolution have been studied with a combination of numerical + simulations and analytical methods (Binney - & Tremaine, 2008). Recently, data describing the positions - and motions of billions of stars from the Gaia + & Tremaine, 2008) for decades. Recently, data describing + the positions and motions of billions of stars from the Gaia satellite (Gaia Collaboration, 2016, 2018) depict a Milky Way (our home galaxy) much farther from equilibrium - that we imagined and beyond the range of many analytic models. N-body - simulations, tracking the evolution of N bodies under their mutual - gravity, are capable of capturing such complexities, but robust links - to theoretical descriptions are still missing.

+ than we imagined and beyond the range of many analytic models. + Simulations that allow collections of bodies to evolve under their + mutual gravity are capable of reproducing such complexities but robust + links to fundamental theoretical explanations are still missing.

Basis Function Expansions (BFE) represent fields as a linear combination of orthogonal functions. BFEs are particularly well-suited - for studies of perturbations around equilibria, such as the evolution - of a galaxy. If the BFE is used to describe the evolution of a galaxy - with time, then one can also study the time series of function - coefficients. For any galaxy simulation, a biorthogonal BFE can fully + for studies of perturbations from equilibrium, such as the evolution + of a galaxy. For any galaxy simulation, a biorthogonal BFE can fully represent the density, potential and forces by time series of - coefficients. This results in huge compression of the information in - the dynamical fields; for example, 1.5 TB of phase space data becomes - 200 MB of coefficient data!

+ coefficients. The coefficients have physical meaning: they represent + the gravitational potential energy in a given function. The variation + the function coefficients in time encodes the dynamical evolution. The + representation of simulation data by BFE results in huge compression + of the information in the dynamical fields; for example, 1.5 TB of + phase space data enumerating the positions and velocities of millions + of particles becomes 200 MB of coefficient data!

For optimal representation, the lowest-order basis function should be similar to the mean or equilibrium profile of the galaxy. This allows for use of the fewest number of terms in the expansion. For @@ -120,10 +121,11 @@ a Creative Commons Attribution 4.0 International License (CC BY representing the cosmologically-motivated Navarro et al. (1997) profile. The EXP software package implements the - adaptive empirical orthogonal function (EOF) basis strategies - originally described in Weinberg + adaptive empirical orthogonal function (EOF; see + [fig:examplecylinder]) + basis strategy originally described in Weinberg (1999) - that match any physical system close to an equilibrium model. The + that matches any physical system close to an equilibrium model. The package includes both a high performance N-body simulation toolkit with computational effort scaling linearly with N (Petersen @@ -140,30 +142,31 @@ a Creative Commons Attribution 4.0 International License (CC BY computation to represent complete series of Basis Function Expansions that describe the variation of any field in space. In the context of galactic - dynamics, these fields may be density, potential, force, or even - velocity fields or chemical tags. By combining the information through - time using SSA - (Golyandina + dynamics, these fields may be density, potential, force, velocity + fields or any intrinsic field produced by simulations such as + chemistry data. By combining the coefficient information through time + using multichannel singular spectral analysis (mSSA, + Golyandina et al., 2001), a non-parametric spectral technique, EXP can deepen our understanding by discovering - the dynamics galaxy evolution directly from simulations and - observations.

-

EXP is a collection of object-oriented C++ + the dynamics of galaxy evolution directly from simulated, and by + analogy, observed data.

+

EXP decomposes a galaxy into multiple bases + for a variety of scales and geometries and is thus able to represent + arbitrarily complex simulation with many components (e.g., disk, + bulge, dark matter halo, satellites). EXP is + able to efficiently summarize the degree and nature of asymmetries + through coefficient amplitudes tracked through time and provide + details at multiple scales. The amplitudes themselves enable + ex-post-facto dynamical discovery. + EXP is a collection of object-oriented C++ libraries with an associated modular N-body code and a suite of - stand-alone analysis applications. The library includes many of the - published recursion relations for specific basis sets; we have - demonstrated that the adaptive method can reproduce all of them to - machine accuracy. By combining multiple bases for different scales and - geometries, EXP, can decompose a galaxy model - based on the geometry and symmetry for any number of components (disk, - bulge, dark-matter halo, satellites). The run of coefficient - amplitudes through time efficiently summarize the degree and nature of - asymmetries and provides detail at multiple scales and enable - ex-post-facto discovery.

-

pyEXP provides full Python interface, - implemented with pybind11 + stand-alone analysis applications.

+

pyEXP provides a full Python interface to + the EXP libraries, implemented with + pybind11 (Jakob - et al., 2017), that provides full interoperability with major + et al., 2017), which provides full interoperability with major astronomical packages including Astropy (Astropy Collaboration, 2013) and Gala @@ -172,31 +175,32 @@ a Creative Commons Attribution 4.0 International License (CC BY available and distributed as accompanying examples and tutorials. The examples and tutorials flatten the - learning curve for adopting BFE tools for generating and analyzing the - significance of coefficients and discovering dynamical relationships + learning curve for adopting BFE tools to generate and analyze the + significance of coefficients and discover dynamical relationships using time series analysis such as mSSA. We provide a full online manual hosted by ReadTheDocs.

The software package brings published – but difficult to implement – applied-math technologies into the astronomical mainstream. - EXP and particularly the Python interface + EXP and the associated Python interface pyEXP accomplish this by providing tools integrated with the Python ecosystem, and in particular are well-suited for interactive Python (Pérez & Granger, 2007) use through (e.g.) Jupyter notebooks (Kluyver - et al., 2016). We anticipate that EXP we - serve as the scaffolding for new imaginative applications in galactic - dynamics, providing a common dynamical language for simulations and - analytic theory.

+ et al., 2016). EXP serves as the + scaffolding for new imaginative applications in galactic dynamics, + providing a common dynamical language for simulations and analytic + theory.

Features and workflow

The core EXP library is built around methods to build the best basis function expansion for an arbitrary data set - in galactic dynamics. The table below lists available basis functions. - All computed bases and resulting coefficient data are stored in HDF5 + in galactic dynamics. The table below lists some of the available + basis functions. All computed bases and resulting coefficient data are + stored in HDF5 (The HDF Group, 2000-2010) format.

@@ -214,9 +218,13 @@ a Creative Commons Attribution 4.0 International License (CC BY sphereSL - Sturm-Liouville computed basis function solutions to - Poisson’s equation for any arbitrary input spherical - density + Sturm-Liouville basis function solutions to Poisson’s + equation for any arbitrary input spherical density + + + bessel + Basis constructed from eigenfunctions of the spherical + Laplacian cylinder @@ -244,22 +252,46 @@ a Creative Commons Attribution 4.0 International License (CC BY -

A picture of some basis functions will be trialled - here

+ +

Example cylinder basis functions, where the color + encodes the amplitude of the function, for an exponential disk with + a scalelength of 3 and a scaleheight of 0.3 in arbitrary units. We + select three functions at low, medium, and higher order + (corresponding to the number of nodes). The color scale has been + normalised such that the largest amplitude is unity in each panel. +

+ +
N-body simulation +

Computing the gravitational potential and forces from a + collection of N particles is typically an expensive endeavour. EXP + reduces the cost by using BFE to compute the potential and forces + such that computational effort scales with the number of particles. + Other modern N-body codes use direct summation + (Wang + et al., 2015) or tree-based solutions + (Springel + et al., 2021), which have computational effort that scales as + N + + 2 + and N log N, respectively. The trade off for BFE solutions comes in + the form of restricted degrees of freedom; for many problems in + near-equilibrium galactic dynamics this is not a problem, but rather + a feature.

Our design includes a wide choice of run-time summary diagnostics, phase-space output formats, dynamically loadable user libraries, and easy extensibility. Stand-alone routines include the - EOF and mSSA methods described above. EXP has - been used, enhanced, and tested for nearly twenty years. Its modular - software architecture makes it extensible and maintainable. This - code base has been tested rigorously over more than a decade and - described in published papers + EOF and mSSA methods described above, and the modular software + architecture of EXP enables users to easily + build and maintain extensions. The EXP code + base is described in published papers (Petersen et al., 2022; Weinberg, - 2023).

+ 2023) and has been used, enhanced, and rigorously tested for + nearly two decades.

The design and implementation of the N-body tools allows for execution on a wide variety of hardware, from personal laptops to high performance computing centers, with communication between @@ -269,14 +301,15 @@ a Creative Commons Attribution 4.0 International License (CC BY CUDA (NVIDIA et al., 2020). Owing to the linear scaling of computational - effort with N, the N-body methods in EXP - deliver performance in N-body simulations previously only accessible - with large dedicated CPU clusters.

+ effort with N and the novel GPU implementation, the N-body methods + in EXP deliver performance in collisionless + N-body simulations previously only accessible with large dedicated + CPU clusters.

The flexible N-body software design allows users to write their - own modules for on-the-fly execution during N-body integration. Such - modules can be powerful methods to design dynamical experiments in + own modules for on-the-fly execution during N-body integration. + Modules enable powerful and intricate dynamical experiments in N-body simulations, further reducing the gap between numerical - simulations and analytic deynamics. The package ships with several + simulations and analytic dynamics. The package ships with several examples, including imposed external potentials, as well as a basic example that can be extended by users.

@@ -293,19 +326,25 @@ a Creative Commons Attribution 4.0 International License (CC BY (Hunter, 2007) and Astropy. We include a verified set of stand-alone routines that read phase-space files from many major cosmological - tree codes and produce BFE-based analyses. The code suite includes - adapters for reading and writing phase space for many of the widely - used cosmology codes, with a base class for developing new ones. - There are multiple way to use the versatile and modular tools in - pyEXP, and we anticipate pipelines that we - have not yet imagined.

+ tree codes (for example, + Springel + et al., 2021) and produce BFE-based analyses. The code suite + includes adapters for reading and writing phase space for many of + the widely used cosmology codes, with a base class for developing + new ones. There are multiple ways to use the versatile and modular + tools in pyEXP, and we anticipate pipelines + that we have not yet imagined. The flexibility of the basis sets + available in EXP greatly enhances the number + of available basis sets implemented in Python (see, e.g. + Price-Whelan, + 2017).

Using pyEXP to analyze time series -

The EXP library includes multiple - different time series analysis tools, documented in the manual. - Here, we briefly highlight one technique that we have already used - in published work: mSSA +

The EXP library includes multiple time + series analysis tools, documented in the manual. Here, we briefly + highlight one technique that we have already used in published work: + mSSA (Johnson et al., 2023; Weinberg @@ -313,12 +352,13 @@ a Creative Commons Attribution 4.0 International License (CC BY the previous tools, mSSA summarizes signals in time that describes dynamically correlated responses and patterns. Essentially, this is BFE in time and space. These temporal and - spatial patterns allow users to better identify dynamical mechanism + spatial patterns allow users to better identify dynamical mechanisms and enable intercomparisons and filtering for features in simulation - suites; e.g. computing the fraction of grand design or bar - excitation. Options for eigenanalysis ensure that decomposition of - large data sets is possible. All mSSA decompositions are saved in - HDF5 format for reuse.

+ suites; e.g. computing the fraction galaxies with grand design + structure or hosting bars. Random-matrix techniques for + singular-value decomposition ensure that analyses of large data sets + is possible. All mSSA decompositions are saved in HDF5 format for + reuse.

@@ -345,6 +385,47 @@ a Creative Commons Attribution 4.0 International License (CC BY http://adsabs.harvard.edu/abs/2008gady.book.....B + + + + SpringelVolker + PakmorRüdiger + ZierOliver + ReineckeMartin + + Simulating cosmic structure formation with the GADGET-4 code + Monthly Notices of the RAS + 202109 + 506 + 2 + https://arxiv.org/abs/2010.03567 + 10.1093/mnras/stab1855 + 2871 + 2949 + + + + + + WangLong + SpurzemRainer + AarsethSverre + NitadoriKeigo + BerczikPeter + KouwenhovenM. B. N. + NaabThorsten + + NBODY6++GPU: ready for the gravitational million-body problem + Monthly Notices of the RAS + 201507 + 450 + 4 + https://arxiv.org/abs/1504.03687 + 10.1093/mnras/stv817 + 4070 + 4080 + + @@ -364,7 +445,7 @@ a Creative Commons Attribution 4.0 International License (CC BY Gaia Collaboration Gaia Data Release 2. Mapping the Milky Way disc kinematics - A&Ap + Astronomy and Astrophysics 201808 616 https://arxiv.org/abs/1804.09380 @@ -427,7 +508,7 @@ a Creative Commons Attribution 4.0 International License (CC BY WeinbergMartin D. New dipole instabilities in spherical stellar systems - MNRAS + Monthly Notices of the RAS 202311 525 4 @@ -444,7 +525,7 @@ a Creative Commons Attribution 4.0 International License (CC BY PetersenMichael S. Using multichannel singular spectrum analysis to study galaxy dynamics - MNRAS + Monthly Notices of the RAS 202103 501 4 @@ -462,7 +543,7 @@ a Creative Commons Attribution 4.0 International License (CC BY KatzNeal EXP: N-body integration using basis function expansions - + Monthly Notices of the RAS 202203 510 4 @@ -481,7 +562,7 @@ a Creative Commons Attribution 4.0 International License (CC BY WeinbergMartin D. Dynamical data mining captures disc-halo couplings that structure galaxies - MNRAS + Monthly Notices of the RAS 202305 521 2 @@ -567,7 +648,7 @@ a Creative Commons Attribution 4.0 International License (CC BY WeinbergMartin D. An Adaptive Algorithm for N-Body Field Expansions - + Astronomical Journal 199901 117 1 @@ -651,7 +732,7 @@ a Creative Commons Attribution 4.0 International License (CC BY HernquistLars An Analytical Model for Spherical Galaxies and Bulges - + Astrophysical Journal 199006 356 10.1086/168845 @@ -666,7 +747,7 @@ a Creative Commons Attribution 4.0 International License (CC BY OstrikerJeremiah P. A Self-consistent Field Method for Galactic Dynamics - + Astrophysical Journal 199202 386 10.1086/171025 @@ -682,7 +763,7 @@ a Creative Commons Attribution 4.0 International License (CC BY WhiteSimon D. M. A Universal Density Profile from Hierarchical Clustering - + Astrophysical Journal 199712 490 2 diff --git a/expui/expMSSA.H b/expui/expMSSA.H index 2fe4d0c3a..6e3563f4c 100644 --- a/expui/expMSSA.H +++ b/expui/expMSSA.H @@ -34,7 +34,10 @@ namespace MSSA void mssa_analysis(); //! MSSA control flags - bool computed, reconstructed, trajectory, useSignChoice; + bool computed, reconstructed, trajectory, useSignChoice, fullRecon; + + //! Last computed number of components + int nlast=-1; //! The reconstructed coefficients for each PC std::map RC; @@ -219,12 +222,14 @@ namespace MSSA @param name is the nmemonic name (see constructor) @param index is one of the specified keys (see constructor) */ - Eigen::MatrixXd wCorrKey(const Key& key); - Eigen::MatrixXd wCorr(const std::string& name, const Key& ckey); + Eigen::MatrixXd wCorrKey(const Key& key, + int nPC=std::numeric_limits::max()); + Eigen::MatrixXd wCorr(const std::string& name, const Key& ckey, + int nPC=std::numeric_limits::max()); //@} //! Compute w-correlation matrix for all channels - Eigen::MatrixXd wCorrAll(); + Eigen::MatrixXd wCorrAll(int nPC=std::numeric_limits::max()); /** Compute the diagnostic computation of contributions per channel In C++, consider using std::tie to get the unpacked values, e.g.: @@ -241,7 +246,7 @@ namespace MSSA void contributionsPNG(); //! Create wcorrlation matricies and output PNG - void wcorrPNG(); + void wcorrPNG(int nPC=std::numeric_limits::max()); //@{ /** diff --git a/expui/expMSSA.cc b/expui/expMSSA.cc index 106857ee2..176ae7715 100644 --- a/expui/expMSSA.cc +++ b/expui/expMSSA.cc @@ -65,12 +65,28 @@ namespace MSSA { (const Eigen::MatrixXd& X, const Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V); - Eigen::MatrixXd expMSSA::wCorrKey(const Key& key) + Eigen::MatrixXd expMSSA::wCorrKey(const Key& key, int nPC) { if (RC.find(key)==RC.end()) { throw std::runtime_error("expMSSA::wCorrKey: no such key"); } + if (nPC<2) { + throw std::runtime_error("expMSSA::wCorrKey: nPC must be >= 2 for a meaningful correlation"); + } + + // Get the number of components + int ncomp = std::min({numW, npc, nPC, static_cast(PC.cols())}); + + // Do the reconstruction + if (not fullRecon or ncomp>nlast) { + std::vector evlist(ncomp); + std::iota(evlist.begin(), evlist.end(), 0); + reconstruct(evlist); + fullRecon = true; + nlast = ncomp; + } + auto R = RC[key]; int numT = R.rows(); @@ -85,17 +101,19 @@ namespace MSSA { else return numT - i + 1; }; - Eigen::MatrixXd ret = Eigen::MatrixXd::Zero(numW, numW); - for (int m=0; m(nPC, numW); + + Eigen::MatrixXd ret = Eigen::MatrixXd::Zero(rank, rank); + for (int m=0; m0.0 and ret(n, n)>0.0) ret(m, n) /= sqrt(ret(m, m)*ret(n, n)); } @@ -103,11 +121,11 @@ namespace MSSA { // Unit diagonal // - for (int m=0; m= 2 for a meaningful correlation"); + } + + // Get the number of components + int ncomp = std::min({numW, npc, nPC, static_cast(PC.cols())}); + + // Do the reconstruction + if (not fullRecon or ncomp>nlast) { + std::vector evlist(ncomp); + std::iota(evlist.begin(), evlist.end(), 0); + reconstruct(evlist); + fullRecon = true; + nlast = ncomp; + } + int numT = RC.begin()->second.rows(); int numW = RC.begin()->second.cols(); int Lstar = std::min(numT - numW, numW); @@ -129,10 +163,12 @@ namespace MSSA { else return numT - i + 1; }; - Eigen::MatrixXd ret = Eigen::MatrixXd::Zero(numW, numW); + int rank = std::min(nPC, numW); + + Eigen::MatrixXd ret = Eigen::MatrixXd::Zero(rank, rank); for (auto R : RC) { - for (int m=0; m0.0 and ret(n, n)>0.0) ret(m, n) /= sqrt(ret(m, m)*ret(n, n)); } @@ -150,19 +186,24 @@ namespace MSSA { // Unit diagonal // - for (int m=0; m= 2 for a meaningful correlation"); + } + int indx = coefDB.index(name); if (indx<0) { std::cout << "No such name <" << name << ">" << std::endl @@ -179,11 +220,11 @@ namespace MSSA { ekey.push_back(0); ekey.push_back(indx); - auto mat = wCorrKey(ekey); + auto mat = wCorrKey(ekey, nPC); ekey[pos] = 1; if (RC.find(ekey)!=RC.end()) { - mat += wCorrKey(ekey); + mat += wCorrKey(ekey, nPC); mat *= 0.5; } return mat; @@ -191,7 +232,7 @@ namespace MSSA { else { auto ekey = key; ekey.push_back(indx); - return wCorrKey(ekey); + return wCorrKey(ekey, nPC); } } @@ -378,7 +419,8 @@ namespace MSSA { // PC = Y * U; - computed = true; + computed = true; + fullRecon = false; reconstructed = false; } @@ -533,6 +575,7 @@ namespace MSSA { } } + fullRecon = false; reconstructed = true; } @@ -983,7 +1026,7 @@ namespace MSSA { return {fw, pt}; } - void expMSSA::wcorrPNG() + void expMSSA::wcorrPNG(int nPC) { #ifdef HAVE_LIBPNGPP { @@ -997,7 +1040,7 @@ namespace MSSA { if (params["distance"]) use_dist = true; for (auto u : RC) { - Eigen::MatrixXd wc = wCorrKey(u.first); + Eigen::MatrixXd wc = wCorrKey(u.first, nPC); png::image< png::rgb_pixel > image(nDim*ndup, nDim*ndup); ColorGradient color; @@ -1027,7 +1070,7 @@ namespace MSSA { } { - Eigen::MatrixXd wc = wCorrAll(); + Eigen::MatrixXd wc = wCorrAll(nPC); png::image< png::rgb_pixel > image(nDim*ndup, nDim*ndup); ColorGradient color; @@ -1378,6 +1421,7 @@ namespace MSSA { // Compute flags // computed = false; + fullRecon = false; reconstructed = false; // Top level parameter flags diff --git a/pyEXP/MSSAWrappers.cc b/pyEXP/MSSAWrappers.cc index b6e556993..5846926df 100644 --- a/pyEXP/MSSAWrappers.cc +++ b/pyEXP/MSSAWrappers.cc @@ -346,6 +346,7 @@ void MSSAtoolkitClasses(py::module &m) { )"); f.def("wCorr", &expMSSA::wCorr, py::arg("name"), py::arg("key"), + py::arg("nPC") = std::numeric_limits::max(), R"( The w-correlation matrix for the selected component and channel key @@ -355,6 +356,8 @@ void MSSAtoolkitClasses(py::module &m) { The name of the selected component. key : tuple The identifier key of the selected channel. + nPC : int + The maximum rank for reconstruction Returns ------- @@ -363,6 +366,12 @@ void MSSAtoolkitClasses(py::module &m) { Notes ----- + The w-correlation matrix needs the reconstructed trajectory matrices for each + of the eigenvalue-PC pairs. Calling this method will recompute the reconstruction + for all eigenvalues up to 'nPC' and return an nPC x nPC matrix. If the 'nPC' + parameter is not specified, it will be set to `numpc` used to construct the + instance. Any previous reconstruction will be overwritten. + Returns the combined cosine+sine correlation for complex types for viewing (e.g., with 'imshow'). @@ -371,6 +380,7 @@ void MSSAtoolkitClasses(py::module &m) { )"); f.def("wCorrKey", &expMSSA::wCorrKey, py::arg("key"), + py::arg("nPC") = std::numeric_limits::max(), R"( Get the w-correlation matrix for the selected component and channel key @@ -386,9 +396,16 @@ void MSSAtoolkitClasses(py::module &m) { Notes ----- - The index key here is 'extended' by the prefixed component index + The index key here is 'extended' by the prefixed component index. - The rows and columns contain distinct cosine and sine indicies if the channel + Computation of the w-correlation matrix needs the reconstructed + trajectory matrices for each of the eigenvalue-PC pairs. Calling + this method will recompute the reconstruction for all eigenvalues up to + order 'npc' and return an (nPC x nPC) matrix. If the 'nPC' parameter is + not specified, it will be set to the `numpc` used in the original + construction. Any prior reconstruction will be overwritten. + + The rows and columns contain distinct cosine and sine indices if the channel is complex valued. This matrix can be visualized using 'imshow' for plotting. @@ -398,9 +415,16 @@ void MSSAtoolkitClasses(py::module &m) { )"); f.def("wCorrAll", &expMSSA::wCorrAll, + py::arg("nPC") = std::numeric_limits::max(), + R"( the w-correlation matrix for all channels in the reconstruction + Parameters + ---------- + nPC : int + The maximum rank for reconstruction + Returns ------- numpy.ndarray @@ -408,8 +432,15 @@ void MSSAtoolkitClasses(py::module &m) { Notes ----- - The w-correlation values range from 0 to 1, where a higher value corresponds to a - stronger correlation. + The w-correlation values range from 0 to 1, where a higher value + corresponds to a stronger correlation. + + Computation of the w-correlation matrix needs the reconstructed + trajectory matrices for each of the eigenvalue-PC pairs. Calling + this method will recompute the reconstruction for all eigenvalues up to + order 'npc' and return an (nPC x nPC) matrix. If the 'nPC' parameter is + not specified, it will be set to the `numpc` used in the original + construction. Any prior reconstruction will be overwritten. See also -------- @@ -419,13 +450,26 @@ void MSSAtoolkitClasses(py::module &m) { )"); f.def("wcorrPNG", &expMSSA::wcorrPNG, + py::arg("nPC") = std::numeric_limits::max(), R"( w-correlation matrices and output PNG image representations + Parameters + ---------- + nPC : int + The maximum rank for reconstruction + Notes ----- - The w-correlation values range from 0 to 1, where a higher value corresponds to a - stronger correlation. + The w-correlation values range from 0 to 1, where a higher value + corresponds to a stronger correlation. + + Computation of the w-correlation matrix needs the reconstructed + trajectory matrices for each of the eigenvalue-PC pairs. Calling + this method will recompute the reconstruction for all eigenvalues up to + order 'npc' and return an (nPC x nPC) matrix. If the 'nPC' parameter is + not specified, it will be set to the `numpc` used in the original + construction. Any prior reconstruction will be overwritten. See also -------- @@ -438,17 +482,18 @@ void MSSAtoolkitClasses(py::module &m) { py::arg("clusters") = 4, py::arg("stride") = 2, R"( - Do a k-means analysis on the reconstructed trajectory matrices for a single channel (specified key value) to - provide grouping insight. A vector of channel indices that identify clusters is return in a vector ordered by PC - index. + Do a k-means analysis on the reconstructed trajectory matrices for a + single channel (specified key value) to provide grouping insight. A + vector of channel indices that identify clusters is return in a vector + ordered by PC index. Parameters ---------- clusters : int, default=4 number of clusters for the k-means analysis stride : int, default=2 - if positive, the initial cluster centers are stride selected from the PC list. If zero, the centers are - selected randomly from the PC list + if positive, the initial cluster centers are stride selected from the PC list. + If zero, the centers are selected randomly from the PC list Returns ------- @@ -459,14 +504,18 @@ void MSSAtoolkitClasses(py::module &m) { Notes ----- - The k-means partitions n vector observations into k clusters in which each observation belongs to the cluster with - the nearest centers while minimizing the variance within each cluster. In this case, the vectors are the full - trajectory matrices and the distance is the distance between the trajectory matricies reconstructed from each - eigentriple from mSSA. The distance used here is the Frobenius distance or matrix norm distance: the square root - of the sum of squares of all elements in the difference between two matrices. - - This version does the analysis for all channels together, the most useful for estimating groups. For individual - contributions by channel, use kmeansChannel. + The k-means partitions n vector observations into k clusters in which + each observation belongs to the cluster with the nearest centers while + minimizing the variance within each cluster. In this case, the vectors + are the full trajectory matrices and the distance is the distance + between the trajectory matrices reconstructed from each eigentriples + from mSSA. The distance used here is the Frobenius distance or matrix + norm distance: the square root of the sum of squares of all elements in + the difference between two matrices. + + This version does the analysis for all channels together, the most + useful for estimating groups. For individual contributions by channel, + use kmeansChannel. )"); f.def("kmeansChannel", &expMSSA::kmeansChannel, @@ -474,9 +523,10 @@ void MSSAtoolkitClasses(py::module &m) { py::arg("clusters") = 4, py::arg("stride") = 2, R"( - Do a k-means analysis on the reconstructed trajectory matrices for a single channel (specified key value) to - provide grouping insight. In most cases, you will want to use the kmeans() version which analyzes all channels - together. + Do a k-means analysis on the reconstructed trajectory matrices for a + single channel (specified key value) to provide grouping insight. In + most cases, you will want to use the kmeans() version which analyzes all + channels together. Parameters ---------- @@ -512,18 +562,18 @@ void MSSAtoolkitClasses(py::module &m) { - F: Each PC's contribution to each channel. The columns are L2 normed. - G: Each channel's contribution to each PC. The rows are L2 normed. - By default, channels for non-zero 'm' are split into cosine and sine components - from the real+imaginary values. + By default, channels for non-zero 'm' are split into cosine and sine + components from the real+imaginary values. - The L2 norm, or Euclidean norm, computes the length of a vector in a multi-dimensional space. - For a vector v = [v1, v2, ..., vn], the L2 norm is calculated as sqrt(v1^2 + v2^2 + ... + vn^2). + The L2 norm, or Euclidean norm, computes the length of a vector in a + multi-dimensional space. For a vector v = [v1, v2, ..., vn], the L2 + norm is calculated as sqrt(v1^2 + v2^2 + ... + vn^2). - The L2 normed views provide a measure of the relative contribution of each PC to each channel - and the relative contribution of each channel to each PC. These contributions can be plotted - using 'imshow'. + The L2 normed views provide a measure of the relative contribution of + each PC to each channel and the relative contribution of each channel to + each PC. These contributions can be plotted using 'imshow'. )"); - f.def("saveState", &expMSSA::saveState, R"( Save the current MSSA state to an HDF5 file