From e24e8c1f2f77ef32845c42af0a963d418765b766 Mon Sep 17 00:00:00 2001 From: Mladen Kolar Date: Sun, 25 Oct 2026 16:12:38 +0400 Subject: [PATCH 1/3] Standardize naming and use master checkouts for workspace dependencies --- .github/workflows/CI.yml | 16 +++--- .gitignore | 6 +++ CONTRIBUTING.md | 26 ++++++++++ NAMING_MIGRATION.md | 19 +++++++ Project.toml | 2 +- README.md | 18 ++++++- benchmark/sigma2/common.jl | 6 +-- docs/Project.toml | 2 +- docs/make.jl | 4 +- docs/src/bases.md | 4 +- docs/src/datasets.md | 4 +- docs/src/fdapace.md | 20 ++++---- docs/src/fpca.md | 22 ++++---- docs/src/index.md | 12 +++++ docs/src/m4_fgm.md | 6 +-- docs/src/m4_report.md | 8 +-- docs/src/mfpca.md | 4 +- docs/src/plotting.md | 34 ++++++------- docs/src/smoothing.md | 13 +++-- docs/src/uncertainty.md | 6 +-- scripts/check_naming.py | 85 +++++++++++++++++++++++++++++++ scripts/instantiate.jl | 6 +++ scripts/master_dependencies.jl | 71 ++++++++++++++++++++++++++ scripts/test_master.jl | 53 +++++++++++++++++++ src/FunctionalDataAnalysis.jl | 7 +-- src/bands.jl | 2 +- src/bases.jl | 10 ++-- src/compat_naming.jl | 8 +++ src/conversions.jl | 2 +- src/covariance.jl | 2 +- src/cross_structure.jl | 17 ++++--- src/derivatives.jl | 10 ++-- src/evaluation.jl | 2 +- src/face.jl | 18 +++---- src/face_sparse.jl | 20 ++++---- src/fpca.jl | 85 +++++++++++++++++-------------- src/fpcbasis.jl | 4 +- src/interface.jl | 2 +- src/mean.jl | 6 +-- src/mfpca.jl | 4 +- src/multivariate_fpca.jl | 27 +++++----- src/pace.jl | 30 +++++------ src/plot_data.jl | 2 +- src/psplines.jl | 8 +-- src/recipes.jl | 32 ++++++------ src/smoothed_fpca.jl | 31 +++++------ src/smoothing.jl | 4 +- src/smoothing2d.jl | 2 +- test/bands.jl | 10 ++-- test/bands_coverage.jl | 2 +- test/bases.jl | 34 ++++++------- test/covariance.jl | 8 +-- test/cross_covariance.jl | 6 +-- test/cross_structure.jl | 24 ++++----- test/data.jl | 10 ++-- test/derivatives.jl | 26 +++++----- test/face.jl | 10 ++-- test/face_sparse.jl | 6 +-- test/fpca.jl | 68 ++++++++++++------------- test/fpcbasis.jl | 14 ++--- test/kernelsmoothing_agreement.jl | 4 +- test/mean.jl | 4 +- test/mfpca.jl | 26 +++++----- test/multivariate_fpca.jl | 68 ++++++++++++------------- test/multivariate_fpca_golden.jl | 6 +-- test/naming.jl | 37 ++++++++++++++ test/pace.jl | 54 ++++++++++---------- test/pace_defaults.jl | 24 ++++----- test/pace_defaults_golden.jl | 6 +-- test/pace_golden.jl | 14 ++--- test/plotting.jl | 16 +++--- test/reference/mfpca_m4.jl | 4 +- test/runtests.jl | 1 + test/smoothed_fpca.jl | 28 +++++----- 74 files changed, 826 insertions(+), 466 deletions(-) create mode 100644 CONTRIBUTING.md create mode 100644 NAMING_MIGRATION.md create mode 100644 scripts/check_naming.py create mode 100644 scripts/instantiate.jl create mode 100644 scripts/master_dependencies.jl create mode 100644 scripts/test_master.jl create mode 100644 src/compat_naming.jl create mode 100644 test/naming.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index aaa3b46..c6bf950 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -2,7 +2,7 @@ name: CI on: push: - branches: [main] + branches: [master] tags: ['*'] pull_request: workflow_dispatch: @@ -29,11 +29,15 @@ jobs: with: version: ${{ matrix.julia-version }} - uses: julia-actions/cache@v3 - - name: Install unregistered KernelSmoothing on Julia 1.10 - if: matrix.julia-version == '1.10' - run: julia --project=. -e 'using Pkg; Pkg.add(url = "https://github.com/mlakolar/KernelSmoothing.jl")' + - name: Check naming conventions + run: ${{ runner.os == 'Windows' && 'python' || 'python3' }} scripts/check_naming.py + - name: Resolve dependencies from master checkouts + run: julia --color=yes scripts/instantiate.jl - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 + - name: Test against dependency master checkouts + run: julia --color=yes scripts/test_master.jl + env: + JULIA_TEST_COVERAGE: "true" - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v7 with: @@ -55,7 +59,7 @@ jobs: version: '1' - uses: julia-actions/cache@v3 - name: Install dependencies - run: julia --project=docs -e 'using Pkg; Pkg.add(url = "https://github.com/mlakolar/KernelSmoothing.jl"); Pkg.develop(path = "."); Pkg.instantiate()' + run: julia scripts/instantiate.jl docs - name: Build and deploy run: julia --project=docs docs/make.jl env: diff --git a/.gitignore b/.gitignore index cade24e..4d0553d 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,9 @@ Manifest.toml *.mem .DS_Store /data/build/ + +# Generated master dependency test environments +.master-test-env/ + +# Local Julia-minor environments are regenerated from master checkouts +/Manifest-v*.toml diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..1902767 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,26 @@ +# Naming and development dependencies + +Types and modules use `UpperCamelCase`. Descriptive functions, fields, keyword +arguments and private helpers use `snake_case`; mutating functions end in `!`. +Keep established numerical abbreviations (`nobs`, `atol`, `rtol`), mathematical +identifiers (`X`, `M`, `Σ`, `λ`) and names required by external interfaces. + +Develop and test against local `master` checkouts of workspace dependencies. +Retain package UUIDs and version numbers during this migration. Use canonical +names in new code; deprecated bindings remain available for downstream users. +Frozen reference data and recorded numerical results are not renamed. + +From the package root, run the complete suite with: + +```sh +julia scripts/test_master.jl +``` + +The runner prepares `.master-test-env/v.` with this package, +its test extras, and local master checkouts of workspace dependencies. Missing +workspace checkouts are cloned from master; existing checkouts must already be +on master. It supports Julia 1.10 and newer and keeps the package `Project.toml` +unchanged. CI uses the same runner. + +Use `julia scripts/test_master.jl --setup-only` to prepare the environment for +focused checks; the runner prints the environment path to use with `--project`. diff --git a/NAMING_MIGRATION.md b/NAMING_MIGRATION.md new file mode 100644 index 0000000..fa677f4 --- /dev/null +++ b/NAMING_MIGRATION.md @@ -0,0 +1,19 @@ +# Naming migration + +Versions and UUIDs are unchanged. The development branch is `master`; the +previous `main` branch remains available. Develop and test with the sibling +`KernelSmoothing.jl` master checkout, as specified in the project sources. + +| Previous public name | Canonical name | +| --- | --- | +| `basismatrix` | `basis_matrix` | +| `modesplot`, `screeplot` | `modes_plot`, `scree_plot` | +| `covplot`, `designplot`, `fitplot`, `bandplot` | `cov_plot`, `design_plot`, `fit_plot`, `band_plot` | +| FPCA keyword and stored key `maxK` | `max_k` | +| Stored key `requested_maxK` | `requested_max_k` | + +Old function names remain deprecated bindings to the same generic functions. +FPCA accepts the old `maxK` keyword and rejects simultaneous non-default old +and new values. Stored selection/settings named tuples retain the old keys +alongside canonical keys during the transition. Numerical defaults and +external StatsAPI, LinearAlgebra and Tables interfaces are unchanged. diff --git a/Project.toml b/Project.toml index 9294fe0..abbd3e9 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,7 @@ StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [sources] -KernelSmoothing = {url = "https://github.com/mlakolar/KernelSmoothing.jl", rev = "master"} +KernelSmoothing = {path = "../KernelSmoothing.jl"} [compat] Artifacts = "1" diff --git a/README.md b/README.md index ead7fd0..495ee2f 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,8 @@ # FunctionalDataAnalysis.jl -[![CI](https://github.com/mlakolar/FunctionalDataAnalysis.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/mlakolar/FunctionalDataAnalysis.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![CI](https://github.com/mlakolar/FunctionalDataAnalysis.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/mlakolar/FunctionalDataAnalysis.jl/actions/workflows/CI.yml?query=branch%3Amaster) [![Docs](https://img.shields.io/badge/docs-dev-blue.svg)](https://mlakolar.github.io/FunctionalDataAnalysis.jl/dev/) -[![codecov](https://codecov.io/gh/mlakolar/FunctionalDataAnalysis.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/mlakolar/FunctionalDataAnalysis.jl) +[![codecov](https://codecov.io/gh/mlakolar/FunctionalDataAnalysis.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/mlakolar/FunctionalDataAnalysis.jl) [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) Functional data analysis in Julia, with planned support for dense, irregular @@ -26,6 +26,20 @@ using Pkg Pkg.add(url = "https://github.com/mlakolar/FunctionalDataAnalysis.jl") ``` +## Development and testing + +Use a `master` checkout for this package and its workspace dependencies. From +the package root, run the complete test suite with: + +```sh +julia scripts/test_master.jl +``` + +The runner resolves sibling master checkouts and test extras in a dedicated +`.master-test-env` environment, leaving the package dependency declarations +unchanged. Missing workspace dependencies are cloned from their master branches. +See [CONTRIBUTING.md](CONTRIBUTING.md) for the naming and development policy. + ## Licence The package code is MIT-licensed (see `LICENSE`). The datasets loaded by diff --git a/benchmark/sigma2/common.jl b/benchmark/sigma2/common.jl index 005b0fe..5c53ea3 100644 --- a/benchmark/sigma2/common.jl +++ b/benchmark/sigma2/common.jl @@ -28,7 +28,7 @@ end # Same whole-curve CurveCV loss as the package, using sufficient statistics # at exact duplicate coordinates within each training/held-out fold. On dense # data this avoids predicting the same pair once per subject. No smoother is -# implemented here; every prediction uses KernelSmoothing.localpolynomial. +# implemented here; every prediction uses KernelSmoothing.local_polynomial. function grouped_cv(selector::CurveCV, fd, pool, kernel, kind) assignments = collect(1:nobs(fd)) if selector.folds != nobs(fd) @@ -52,7 +52,7 @@ function grouped_cv(selector::CurveCV, fd, pool, kernel, kind) end cv_predictions(train, held, h, kernel, ::Val{:covariance}) = - KernelSmoothing.localpolynomial((train.x, train.y), train.z, hcat(held.x, held.y); + KernelSmoothing.local_polynomial((train.x, train.y), train.z, hcat(held.x, held.y); bandwidth=h, kernel, weights=train.weights) cv_predictions(train, held, h, kernel, ::Val{:mean}) = FDA._mean_predict(train.x, train.z, train.weights, held.x, @@ -115,7 +115,7 @@ function estimate(prepared, estimator) # Benchmark-only FGM-inspired reference: use the package eigen step, # all positive components, and FGM's relative cutoff, not an FPCA K cap. system = FDA._fpca_eigen(surface, quadrature_weights(grid), mean.values; - maxK=length(grid), eigen_tol=100eps(Float64), sign_tol=1e-12) + max_k=length(grid), eigen_tol=100eps(Float64), sign_tol=1e-12) surface = (system.eigenfunctions .* system.eigenvalues') * system.eigenfunctions' surface = (surface + surface')/2 # FGM bsplinekit 6f0eac3, src/doubly.jl:223-226: arithmetic mean on diff --git a/docs/Project.toml b/docs/Project.toml index 8cb4d98..e91a610 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,7 +8,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [sources] -KernelSmoothing = {url = "https://github.com/mlakolar/KernelSmoothing.jl", rev = "master"} +KernelSmoothing = {path = "../../KernelSmoothing.jl"} FunctionalDataAnalysis = {path = ".."} [compat] diff --git a/docs/make.jl b/docs/make.jl index db7fbc4..7216245 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ import StatsAPI makedocs(; sitename = "FunctionalDataAnalysis.jl", modules = [FunctionalDataAnalysis], - format = Documenter.HTML(; edit_link = "main"), + format = Documenter.HTML(; edit_link = "master"), pages = [ "Home" => "index.md", "Interface" => "interface.md", @@ -32,5 +32,5 @@ makedocs(; deploydocs(; repo = "github.com/mlakolar/FunctionalDataAnalysis.jl.git", - devbranch = "main", + devbranch = "master", ) diff --git a/docs/src/bases.md b/docs/src/bases.md index b39fc36..806e048 100644 --- a/docs/src/bases.md +++ b/docs/src/bases.md @@ -8,7 +8,7 @@ julia> using FunctionalDataAnalysis julia> b = MonomialBasis((0, 1), 2); -julia> basismatrix(b, [0.0, 0.5, 1.0]) +julia> basis_matrix(b, [0.0, 0.5, 1.0]) 3×3 Matrix{Float64}: 1.0 0.0 0.0 1.0 0.5 0.25 @@ -34,7 +34,7 @@ BSplineBasis FourierBasis MonomialBasis nbasis -basismatrix +basis_matrix gram penalty ``` diff --git a/docs/src/datasets.md b/docs/src/datasets.md index bfb9079..cef0c8e 100644 --- a/docs/src/datasets.md +++ b/docs/src/datasets.md @@ -43,7 +43,7 @@ growth.covariates ``` ```@example datasets -model = StatsAPI.fit(FPCA, growth.data; ncomponents=2, maxK=2) +model = StatsAPI.fit(FPCA, growth.data; ncomponents=2, max_k=2) mu = mean_data(model) plot(mu.grid, mu.values; title="Mean height: four Berkeley sample curves", xlabel="Age (years)", ylabel="Height (cm)") @@ -53,7 +53,7 @@ savefig("datasets-growth-mean.svg"); nothing # hide ![](datasets-growth-mean.svg) ```@example datasets -modesplot(model, 1; xlabel="Age (years)", ylabel="Height (cm)") +modes_plot(model, 1; xlabel="Age (years)", ylabel="Height (cm)") savefig("datasets-growth-modes.svg"); nothing # hide ``` diff --git a/docs/src/fdapace.md b/docs/src/fdapace.md index b28424d..fd7e98c 100644 --- a/docs/src/fdapace.md +++ b/docs/src/fdapace.md @@ -73,7 +73,7 @@ identifies the matching conventions and the differences retained by design. | Choice | fdapace 0.6.0: executable source | FunctionalDataAnalysis.jl and reason | |:--|:--|:--| | Work grid / `nRegGrid` | `R/SetOptions.R`, `SetOptions`: Dense and DenseWithMV use the number of unique pooled times after `signif(t,6)`; Sparse and `p>>n` use 51. `R/FPCA.R`, `FPCA` then creates an equally spaced grid over the observed range. | Dense uses the supplied data grid and its stored weights, without rounding or resampling. It preserves the data representation. PACE uses 51 equally spaced points over the observed range unless a work grid is supplied, including when PACE is explicitly requested for common-grid data. | -| `maxK` | `SetOptions`: `min(nRegGrid-2,n-2)`, additionally capped at 20 **only for `methodMuCovEst="smooth"`**; values below 1 are reset to 1. An explicit `maxK` above `nRegGrid-2` is reduced with a message. | The dense default is `max(1,min(m-2,n-2))`; PACE uses `max(1,min(20,m-2,n-2))`, matching the respective audited routes. An explicit dense cap may use up to `min(m,n-1)` and an explicit PACE cap up to `m`; requested components above the cap or retained positive rank throw. | +| `max_k` | `SetOptions`: `min(nRegGrid-2,n-2)`, additionally capped at 20 **only for `methodMuCovEst="smooth"`**; values below 1 are reset to 1. An explicit `max_k` above `nRegGrid-2` is reduced with a message. | The dense default is `max(1,min(m-2,n-2))`; PACE uses `max(1,min(20,m-2,n-2))`, matching the respective audited routes. An explicit dense cap may use up to `min(m,n-1)` and an explicit PACE cap up to `m`; requested components above the cap or retained positive rank throw. | | `methodSelectK` | `SetOptions`: `"FVE"`. `R/SelectK.R`, `SelectK` also accepts `"AIC"`, `"BIC"`, or a positive integer. | `ncomponents=FVE(0.99)` or a positive integer. PACE also accepts `AIC()` and `BIC()`. | | `FVEthreshold` | `SetOptions`: 0.99 for `FPCA`. The standalone `SelectK` function separately defaults to 0.95. | `FVE(0.99)` follows the `FPCA` entry point. | | Kernel | `SetOptions` initially chooses `"epan"` for Dense and `"gauss"` for every other data type, including DenseWithMV. A later invalid-kernel recovery branch chooses Epanechnikov for DenseWithMV. | The dense sample mean and covariance require no kernel. Optional presmoothing records the supplied smoother and kernel. PACE defaults to Gaussian, following the audited non-dense default. Direct pooled-smoother calls retain their Epanechnikov default; every kernel remains overridable. | @@ -86,7 +86,7 @@ identifies the matching conventions and the differences retained by design. | Smoothed noise estimator | `PC_CovE` uses the rotated diagonal; `SetOptions` supplies the middle-half `rotationCut=c(.25,.75)`, while the internal formal default is `c(0,1)`. | `SmoothedCovariance` and PACE default to the rotated diagonal over the full observed domain, selected by the issue #17 Discovery benchmark. Explicit `noise_method=:surface_diagonal` and `noise_interval=:middle_half` retain the milestone 2 rule. | | `rho` | `SetOptions`: `methodRho="vanilla"`, `userRho=NULL`. This means no extra regularization; details are below. | `rho` regularization is not implemented. No ridge is silently added. | -The remembered claim that fdapace always uses `maxK=min(20,m-2)` is therefore +The remembered claim that fdapace always uses `max_k=min(20,m-2)` is therefore incorrect for 0.6.0. So is a description of GCV as its default bandwidth selector. The dense component cap follows the audited fdapace default, replacing the milestone plan's cap of 20. @@ -215,7 +215,7 @@ conditioned system is therefore not automatically regularized. An unused solver. The covariance used in that `Sigma_Y` is the **full positive-semidefinite -fitted covariance up to `maxK`**, not only the selected `K` components: +fitted covariance up to `max_k`**, not only the selected `K` components: `GetEigenAnalysisResults` constructs `eigObj$fittedCov` before selecting K, and `FPCA` passes that object to `GetCEScores`. `FVEfittedCov` controls the separate user-returned `fittedCovUser`; it does not change this score @@ -330,7 +330,7 @@ The Discovery benchmark led to `:rotated` with `:full` as the package defaults. `R/GetEigenAnalysisResults.R`, `GetEigenAnalysisResults` computes ordinary matrix eigenvectors of the covariance surface `G`. It keeps eigenvalues `d>=0`, including exact zeros, drops strictly negative ones, and truncates -the remaining descending list to `maxK`. It does not record a dropped +the remaining descending list to `max_k`. It does not record a dropped negative mass or apply a relative positive-eigenvalue cutoff. It reports eigenvalues `lambda=h*d`, with the first regular-grid spacing @@ -348,7 +348,7 @@ eigenvalues above `eigen_tol*maximum(abs,spectrum)` with default `eigen_tol=1e-10`, and records the number, absolute mass, and share of the strictly negative spectrum. The relative cutoff is an explicit numerical choice, rather than a statistical tuning rule from fdapace. Negative mass -is measured before the `maxK` cap. +is measured before the `max_k` cap. The sign agrees with fdapace when the mean overlap is nonzero. If its absolute value is at most @@ -358,7 +358,7 @@ This fallback follows the convention used by FunctionalGraphicalModels.jl. Both numerical tolerances are stored on the fit. They are explicit package choices that the milestone leaves to the implementation. -Both implementations compute FVE **after** the `maxK` truncation: +Both implementations compute FVE **after** the `max_k` truncation: `cumsum(lambda)/sum(lambda)` over that retained positive spectrum. Thus the last retained value is one even if the cap discarded positive components. `GetEigenAnalysisResults` first uses `>=` to find its threshold crossing; @@ -375,7 +375,7 @@ which is distinct from the full fitted covariance used by CE scoring. When `methodSelectK` is numeric, `SetOptions` first sets `FVEthreshold=1`, and `SelectK` then checks the fixed requested K against available components. -PACE stores the criterion for every candidate from 1 through `maxK`, then +PACE stores the criterion for every candidate from 1 through `max_k`, then applies the first-increase rule. Beyond the retained positive rank, extra zero-variance components leave the likelihood unchanged and increase only the penalty. It additionally records the @@ -433,7 +433,7 @@ estimation. These fixtures predate the rotated default; fdapace's enclosing The fixed-bandwidth run estimates its mean and covariance from the sparse fixture. It uses `userBwMu=0.30`, `userBwCov=0.40`, `kernel="epan"`, `useBinnedData="OFF"`, `useBinnedCov=FALSE`, observation weighting, -`methodXi="CE"`, `methodRho="vanilla"`, fixed `K=3`, and `maxK=5`. +`methodXi="CE"`, `methodRho="vanilla"`, fixed `K=3`, and `max_k=5`. Its 41-point work grid equals the observed-time union. Thus mean evaluation before forming residual products agrees in both implementations. Public Julia fits reproduce those mean and covariance values at relative error @@ -488,7 +488,7 @@ recreate `test/data/m3_dense_input.csv` and `test/data/m3_fdapace_dense.csv`. Each CSV begins with R and fdapace versions. The script has no random draws. It uses 18 rank-five curves on 21 grid points, `dataType="Dense"`, cross-sectional estimation, -`error=FALSE`, integration scores, `maxK=5`, no binning, and a fixed K of 3. +`error=FALSE`, integration scores, `max_k=5`, no binning, and a fixed K of 3. A second `FVE(0.99)` run selects K=4; a full-K run records the full FVE curve. The first three reference eigenvalues are approximately `4.448802620311068`, `1.165542244194415`, and `0.3479570319234146`. @@ -699,7 +699,7 @@ Two source issues require explicit corrections to exercise native `fitted.FPCA` on the milestone 3 supplied-input sparse fixture: 1. `SubsetFPCA` truncates `phi`, `xiEst`, and `lambda`, but leaves `xiVar` - at the full retained rank. With selected K smaller than `maxK`, native + at the full retained rank. With selected K smaller than `max_k`, native fitted bands fail with nonconformable matrices. The reference script verifies that failure and subsets only `xiVar` to its selected upper-left block. The CE denominator still contains all five retained components. diff --git a/docs/src/fpca.md b/docs/src/fpca.md index 74c3533..02fe7c2 100644 --- a/docs/src/fpca.md +++ b/docs/src/fpca.md @@ -68,8 +68,8 @@ some older help text. See [Comparison with fdapace](fdapace.md) for details. | `grid`, PACE | 51 equally spaced points over the observed time range. | fdapace [`nRegGrid=51` and work-grid construction](https://github.com/cran/fdapace/blob/0.6.0/R/FPCA.R), with the count set in [`SetOptions`](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass a strictly increasing `grid` covering all observations. | | `grid`, dense | The data grid and its stored quadrature weights. | Preserve the data representation. A new grid requires an explicit data conversion before fitting. | | `ncomponents` | `FVE(0.99)`. | fdapace [`FVEthreshold`](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass `FVE(threshold)` or a fixed integer; PACE also accepts `AIC()` and `BIC()` from [Yao, Müller, and Wang (2005), section 2.5](https://doi.org/10.1198/016214504000001745). | -| `maxK`, PACE | `max(1, min(20, g-2, n-2))`. | fdapace [`SetOptions`, smoothing route](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass an integer no larger than `g`. | -| `maxK`, dense | `max(1, min(g-2, n-2))`, without a cap of 20. | fdapace [`SetOptions`, cross-sectional route](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). An explicit cap may be as large as `min(g, n-1)`, the sample covariance rank bound. | +| `max_k`, PACE | `max(1, min(20, g-2, n-2))`. | fdapace [`SetOptions`, smoothing route](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass an integer no larger than `g`. | +| `max_k`, dense | `max(1, min(g-2, n-2))`, without a cap of 20. | fdapace [`SetOptions`, cross-sectional route](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). An explicit cap may be as large as `min(g, n-1)`, the sample covariance rank bound. | | `kernel`, PACE | `KernelSmoothing.GaussianKernel()`. | fdapace [`SetOptions`, non-dense data](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass another KernelSmoothing kernel. Unlike compact kernels, the Gaussian has no empty windows from finite support. Direct `SmoothedMean` and `SmoothedCovariance` calls retain their Epanechnikov default. | | `weighting`, PACE | `:obs`, equal weight per observation. | fdapace [default smoothing weighting](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass `:subj` for equal curve weights. Neither weighting scheme dominates in all designs; see [Zhang and Wang (2016)](https://doi.org/10.1214/16-AOS1446). | | `mean_bandwidth`, PACE | `0.05r`. | fdapace [`methodBwMu="Default"`, `userBwMu`](https://github.com/cran/fdapace/blob/0.6.0/R/SetOptions.R). Pass a number, `CurveCV`, or a supported KernelSmoothing regression bandwidth selector. | @@ -84,9 +84,9 @@ some older help text. See [Comparison with fdapace](fdapace.md) for details. | `presmooth`, dense | `nothing`. | Dense noise-free approximation in [Ramsay and Silverman (2005), chapter 8](https://doi.org/10.1007/b98888). Pass a `LocalPolynomial` or `PenalizedSpline` object, or use `method=:pace`. | All tuning settings are stored on the fit and displayed by `show`. The -`selection` record stores the rule, FVE threshold when applicable, `maxK`, +`selection` record stores the rule, FVE threshold when applicable, `max_k`, chosen component count, and rule scores. `all_eigenvalues` and `fve(model)` -refer to every retained positive eigenvalue up to `maxK`, before component +refer to every retained positive eigenvalue up to `max_k`, before component selection. Negative eigenvalues are removed before selection. Their count, absolute total, and share of the absolute spectrum are recorded in `dropped`. Values at or below the positive relative tolerance are excluded from the FVE @@ -501,7 +501,7 @@ selection metadata, the origin of every supplied or estimated input, and the noise truncation details are retained on the fit. Inputs are copied. Let ``\widetilde G`` denote the positive-semidefinite projection of the work-grid -surface over **all** positive components retained up to `maxK`. Negative +surface over **all** positive components retained up to `max_k`. Negative eigenvalues are dropped and reported in `model.dropped`. The selected component count ``K`` can be smaller than this retained rank. At subject ``i``'s observed times, let ``\Phi_i`` contain the selected eigenfunctions and @@ -596,7 +596,7 @@ import KernelSmoothing model = StatsAPI.fit(FPCA, simulation.data; method=:pace, grid=collect(range(extrema(simulation.data.times)...; length=61)), - ncomponents=2, maxK=6, + ncomponents=2, max_k=6, kernel=KernelSmoothing.GaussianKernel(), weighting=:subj, mean_bandwidth=0.15, covariance_bandwidth=0.25, diagonal_bandwidth=0.20, binning=nothing) @@ -698,8 +698,8 @@ Eigenanalysis, signs, and numerical tolerances are shared with the dense route. For `AIC()` and `BIC()`, each candidate ``k`` uses its own ``k``-component covariance plus ``\sigma^2 I`` in the Gaussian likelihood. The criterion is minus twice the log likelihood plus ``2k`` or ``k\log(n)``, respectively. -`selection.scores` stores the full criterion curve for `k=1:maxK`; -`selection.candidates` records those counts. When `maxK` exceeds the retained +`selection.scores` stores the full criterion curve for `k=1:max_k`; +`selection.candidates` records those counts. When `max_k` exceeds the retained positive rank, the additional components have zero variance, leave the likelihood unchanged, and increase only the penalty. Selection stops at the first increase, following fdapace. The chosen count @@ -781,7 +781,7 @@ node count throws. For example: ```julia model = StatsAPI.fit(FPCA, mfd; - basis=:marginal, ncomponents=[2, 3], maxK=[4, 5], + basis=:marginal, ncomponents=[2, 3], max_k=[4, 5], grid=[work_grid_node_1, work_grid_node_2]) ``` @@ -805,12 +805,12 @@ pooled surface and node scores are The weighted eigenstep and positive-eigenvalue cutoff are the same as for univariate FPCA. `ncomponents` is an integer or `FVE(threshold)`, defaulting -to `FVE(0.99)`; `maxK` follows +to `FVE(0.99)`; `max_k` follows the dense route. Selection uses the pooled spectrum. The sign rule uses the average of the node mean curves as its common reference and records `sign_reference=:average_node_mean`. Near-zero overlap uses the same largest-absolute-entry fallback as univariate FPCA. The pooled tuning -keywords are `ncomponents`, `maxK`, `eigen_tol`, and `sign_tol`; +keywords are `ncomponents`, `max_k`, `eigen_tol`, and `sign_tol`; `method=:auto` or `:dense` is accepted. Pooled presmoothing, custom work grids, and PACE smoothing options throw because no corresponding pooled method is implemented. Each node keeps its own sample mean, and diff --git a/docs/src/index.md b/docs/src/index.md index 03f5a31..0060d55 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -27,6 +27,18 @@ Pkg.add(url = "https://github.com/mlakolar/KernelSmoothing.jl") # required on Ju Pkg.add(url = "https://github.com/mlakolar/FunctionalDataAnalysis.jl") ``` +## Development tests + +From a master checkout's package root, run: + +```sh +julia scripts/test_master.jl +``` + +This prepares the dedicated master-based test environment, including test extras +and sibling workspace dependencies, and runs the complete suite. Package version +numbers and dependency declarations remain unchanged. + ## Quadrature Use [`quadrature_weights`](@ref) to integrate sampled curves on a finite, diff --git a/docs/src/m4_fgm.md b/docs/src/m4_fgm.md index 3bbc422..4c457c1 100644 --- a/docs/src/m4_fgm.md +++ b/docs/src/m4_fgm.md @@ -33,7 +33,7 @@ uses these explicit settings: | Binning | `binning=nothing` | No binning | | Measurement noise variances | `[0.04, 0.07]` | `[0.04, 0.07]` | | Selected components | `ncomponents=2` | `ncomponents=2` | -| Covariance rank for conditional scores | `maxK=21` | All retained positive eigenvalues | +| Covariance rank for conditional scores | `max_k=21` | All retained positive eigenvalues | | Relative eigenvalue cutoff | `eigen_tol=100eps(Float64)` | `100eps(Float64)` | The sparse test compares smoothed means, selected eigenvalues, @@ -45,8 +45,8 @@ correction. The fixture expresses FGM's settings explicitly rather than using FDA's univariate defaults. In particular, FDA's default relative eigenvalue cutoff is `1e-10`, whereas FGM uses `100eps(Float64)`. FDA also limits the covariance -rank through `maxK`; the comparison retains the complete positive spectrum -by setting `maxK` to the work-grid length. Both packages use all retained +rank through `max_k`; the comparison retains the complete positive spectrum +by setting `max_k` to the work-grid length. Both packages use all retained covariance components in the conditional-score denominator, even though only two scores are returned. FGM's noise-estimation rule differs from FDA's; the fixture supplies noise variances so that neither rule enters the diff --git a/docs/src/m4_report.md b/docs/src/m4_report.md index 95f17df..803b429 100644 --- a/docs/src/m4_report.md +++ b/docs/src/m4_report.md @@ -43,8 +43,8 @@ noise estimate, and smoothing provenance. | Pooled `method` | `:auto`, resolving to dense. | Pooled route requires complete common-grid data; `:dense` is also accepted. | | Marginal `grid` | Data grid for dense; 51 equally spaced observed-range points for PACE. | Dense data representation; fdapace 0.6.0 `nRegGrid` and `FPCA.R`. | | `ncomponents` | `FVE(0.99)` per node or on the pooled spectrum. | Task API and fdapace 0.6.0 `FVEthreshold`. | -| `maxK`, dense and pooled | `max(1,min(g-2,n-2))`; explicit caps no larger than `min(g,n-1)`. | Task requires the dense cap; fdapace 0.6.0 cross-sectional `SetOptions.R`. | -| `maxK`, marginal PACE | `max(1,min(20,g-2,n-2))`. | fdapace 0.6.0 smoothing-route `SetOptions.R`. | +| `max_k`, dense and pooled | `max(1,min(g-2,n-2))`; explicit caps no larger than `min(g,n-1)`. | Task requires the dense cap; fdapace 0.6.0 cross-sectional `SetOptions.R`. | +| `max_k`, marginal PACE | `max(1,min(20,g-2,n-2))`. | fdapace 0.6.0 smoothing-route `SetOptions.R`. | | Marginal PACE mean/covariance | Estimate missing inputs with existing local-linear smoothers. | Yao, Müller, and Wang (2005), inherited milestone 3 route. | | Marginal PACE `noise_variance` | Covariance smoother's nonnegative estimate; supplied surfaces require an explicit value. | Inherited milestone 3 contract and fdapace covariance smoothing. | | Marginal PACE `kernel` | Gaussian when a smoother runs. | fdapace 0.6.0 non-dense `SetOptions.R`. | @@ -71,7 +71,7 @@ remain separate from per-node score variances. equals the number of nodes. A vector of grids denotes per-node grids. Other vector-valued keyword containers must have one entry per node. Per-node `nothing` leaves that node's univariate default active. -- Pooled `maxK` retains the task's dense-route bound `min(g,n-1)`. The +- Pooled `max_k` retains the task's dense-route bound `min(g,n-1)`. The frozen FGM implementation allows `min(g,p*(n-1))` because averaging node covariances can increase rank. FDA follows the requested dense cap; this restriction is explicit and does not imply that the pooled @@ -83,7 +83,7 @@ remain separate from per-node score variances. training mean curves. Near-zero overlap uses the existing largest-entry fallback, with the first entry breaking ties. This choice affects only signs, not the pooled covariance or score variances. -- Pooled tuning accepts shared `ncomponents`, `maxK`, `eigen_tol`, and +- Pooled tuning accepts shared `ncomponents`, `max_k`, `eigen_tol`, and `sign_tol`, with `method=:auto` or `:dense`. Presmoothing, custom work grids, and PACE-specific options are rejected explicitly. The task did not define pooled presmoothing or a new pooled work-grid method. diff --git a/docs/src/mfpca.md b/docs/src/mfpca.md index 34859a2..10ccabe 100644 --- a/docs/src/mfpca.md +++ b/docs/src/mfpca.md @@ -198,7 +198,7 @@ strict given-expansion comparisons. | Eager fitting and bootstrap | Disabled, consistent with R defaults; reconstruction is an explicit operation. No bootstrap API. | The inherited first-step defaults include routing, work grid, quadrature, -`FVE(0.99)`, `maxK`, no dense presmoothing, sample dense mean/covariance, +`FVE(0.99)`, `max_k`, no dense presmoothing, sample dense mean/covariance, estimated PACE mean/covariance/noise, Gaussian kernel, observation weighting, 5%/10% range bandwidths, diagonal bandwidth inheritance, automatic binning, `eigen_tol=1e-10`, and `sign_tol=1e-12`. The linked univariate docstring/audit @@ -210,7 +210,7 @@ Recreation commands (R is only needed for regeneration): ```sh julia --project=. test/reference/mfpca_m4.jl Rscript test/reference/mfpca_m4.R -julia --project=. -e 'include("test/mfpca.jl")' +julia scripts/test_master.jl ``` Reference generation uses R **4.5.1**, MFPCA **1.3-11**, funData **1.3-9**, diff --git a/docs/src/plotting.md b/docs/src/plotting.md index ad28f47..84777a0 100644 --- a/docs/src/plotting.md +++ b/docs/src/plotting.md @@ -36,7 +36,7 @@ functions = fourier_eigenfunctions((0.0, 1.0), 3) dense = simulate_kl(Xoshiro(54), 60; mean=t -> 1 + t, eigenfunctions=functions, eigenvalues=[2.0, 0.7, 0.2], sampling=RegularSampling(grid)).data -model = StatsAPI.fit(FPCA, dense; ncomponents=2, maxK=3) +model = StatsAPI.fit(FPCA, dense; ncomponents=2, max_k=3) plot(dense; curves=1:8, legend=false) savefig("plotting-curves.svg"); nothing # hide ``` @@ -50,7 +50,7 @@ lines are a display convention, with no interpolated points in `curve_data`. sparse = simulate_kl(Xoshiro(55), 60; mean=t -> 1 + t, eigenfunctions=functions, eigenvalues=[2.0, 0.7, 0.2], noise_sd=0.2, sampling=SparseSampling((0.0, 1.0), 5:9)).data -pace = StatsAPI.fit(FPCA, sparse; grid, ncomponents=2, maxK=3, +pace = StatsAPI.fit(FPCA, sparse; grid, ncomponents=2, max_k=3, mean_bandwidth=0.3, covariance_bandwidth=0.4, noise_variance=0.04) plot(sparse; curves=1:5, legend=false) savefig("plotting-sparse.svg"); nothing # hide @@ -89,26 +89,26 @@ savefig("plotting-mean.svg"); nothing # hide ![](plotting-mean.svg) -`modesplot` shows the mean and the four default mode curves. The title gives +`modes_plot` shows the mean and the four default mode curves. The title gives the eigenvalue and the component's individual FVE. Joint `MFPCA` and marginal multivariate FPCA use one panel per node. Pooled multivariate modes use the shared pooled eigenvalue, rather than each node's score variance. ```@example plotting -modesplot(model, 1) +modes_plot(model, 1) savefig("plotting-modes.svg"); nothing # hide ``` ![](plotting-modes.svg) -`screeplot` shows every positive eigenvalue retained up to the fit's cap, +`scree_plot` shows every positive eigenvalue retained up to the fit's cap, including unselected components. FVE uses this full retained spectrum. The vertical line shows selected K. The horizontal line appears only when the fit used an FVE threshold. A marginal multivariate fit has per-node scree panels; pooled and joint fits have a shared spectrum. ```@example plotting -screeplot(model) +scree_plot(model) savefig("plotting-scree.svg"); nothing # hide ``` @@ -122,18 +122,18 @@ retained components, including those beyond selected K. Both surfaces are stored by FPCA; plotting does not rerun the eigendecomposition. ```@example plotting -covplot(model; surface=:fitted) +cov_plot(model; surface=:fitted) savefig("plotting-covariance.svg"); nothing # hide ``` ![](plotting-covariance.svg) -For a three-dimensional view, use `covplot(model; seriestype=:surface)`. +For a three-dimensional view, use `cov_plot(model; seriestype=:surface)`. The same recipe accepts `SmoothedCovariance`, `CrossCovariance`, and `SmoothedCrossCovariance`. ```@example plotting -designplot(sparse) +design_plot(sparse) savefig("plotting-design.svg"); nothing # hide ``` @@ -147,7 +147,7 @@ under the default and one pair with the diagonal enabled. Empty pair arrays are valid. No dense-grid approximation or time binning is used. ```@example plotting -fitplot(pace, 1) +fit_plot(pace, 1) savefig("plotting-fit.svg"); nothing # hide ``` @@ -155,7 +155,7 @@ savefig("plotting-fit.svg"); nothing # hide Fitted paths reuse training scores and `reconstruct`. The points are copied original observations, including when a dense fit used presmoothing. -`fitplot(model, i; grid=newgrid)` chooses a finite increasing reconstruction +`fit_plot(model, i; grid=newgrid)` chooses a finite increasing reconstruction grid inside the fitted domain. Joint MFPCA uses node panels, with either a shared grid or a vector of node grids. Plotting does not run a new score solve. @@ -196,7 +196,7 @@ displays their settings; no new fitted estimator is introduced. ## API Precomputed fitted-curve and mean uncertainty results also support -[`band_data`](@ref) and [`bandplot`](@ref). See [Uncertainty](uncertainty.md) +[`band_data`](@ref) and [`band_plot`](@ref). See [Uncertainty](uncertainty.md) for the covered targets and assumptions. Plotting bands never runs a new fit or bootstrap. Multivariate fitted bands use node panels, and `coverage` on mean-band plots selects among the types already computed. @@ -209,11 +209,11 @@ scree_data covariance_data design_data fitted_curve_data -modesplot -screeplot -covplot -designplot -fitplot +modes_plot +scree_plot +cov_plot +design_plot +fit_plot ``` Derivative fits support `curve_data(d; curves=:all, grid=nothing)` and diff --git a/docs/src/smoothing.md b/docs/src/smoothing.md index a9bd8ba..793be86 100644 --- a/docs/src/smoothing.md +++ b/docs/src/smoothing.md @@ -46,14 +46,17 @@ Pointwise confidence bands for the pooled mean are outside this API. They require cluster-robust variances, which KernelSmoothing.jl does not yet provide. Covariance and noise variance are estimated as described below. -KernelSmoothing.jl is currently unregistered. On Julia 1.10, install it before -instantiating this project (newer Julia versions use the root `[sources]`): +KernelSmoothing.jl is currently unregistered. From this package's development +checkout, prepare the sibling master dependencies on Julia 1.10 and newer with: -```julia -using Pkg -Pkg.add(url="https://github.com/mlakolar/KernelSmoothing.jl") +```sh +julia scripts/instantiate.jl ``` +The bootstrap develops the sibling KernelSmoothing master checkout and clones it +from master if missing. Run `julia scripts/test_master.jl` for the complete suite +in its dedicated environment with test extras. + KernelSmoothing must be registered before FunctionalDataAnalysis can be registered. ```@docs diff --git a/docs/src/uncertainty.md b/docs/src/uncertainty.md index 36dd751..df16952 100644 --- a/docs/src/uncertainty.md +++ b/docs/src/uncertainty.md @@ -10,7 +10,7 @@ independent curves. These have different targets and assumptions. bands = fitted_bands(pace; alpha=0.05, coverage=:interval, curves=[1, 3]) simultaneous = fitted_bands(pace; coverage=:band) # After loading Plots: -bandplot(bands) +band_plot(bands) ``` For curve i, the stored conditional score covariance is Ωᵢ. The variance @@ -69,7 +69,7 @@ uncertainty = mean_bands(fd, m; B=1000, rng=MersenneTwister(42)) # Recreate the bootstrap from its stored initial RNG state: repeated = mean_bands(fd, m; B=uncertainty.B, rng=uncertainty.rng_state) # After loading Plots: -bandplot(uncertainty; coverage=:band) +band_plot(uncertainty; coverage=:band) ``` The bootstrap samples whole curves with replacement, retaining all their @@ -176,5 +176,5 @@ MeanBands fitted_bands mean_bands band_data -bandplot +band_plot ``` diff --git a/scripts/check_naming.py b/scripts/check_naming.py new file mode 100644 index 0000000..e112dd8 --- /dev/null +++ b/scripts/check_naming.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +"""Check new Julia declarations; compatibility files and mathematical names are explicit exceptions.""" +import argparse +from pathlib import Path +import re + +# Matrix caches and scientific notation retain their conventional identifiers. +MATHEMATICAL = {'A','B','C','D','G','H','I','J','K','L','M','N','P','Q','R','S','T','U','V','W','X','Y','Z', + 'Ax','At','Bt','Xt','Yt','XtY','XtX','sUtY','Ws','SX','SY','Xn','Yn','SXn','SYn', + 'Xa','Ya','SXa','SYa','Xr','Yr','posA','posR','T0','F','Sx','Sy','Sxy','Sxx','Syy','AD','Vt','AU','D_calls','Dt_calls','Dx','Dxp','Dty','Dtyp','Bpath'} +FUNCTION = re.compile(r'^\s*(?:function\s+)([\w.!]+)\s*\(|^\s*([\w.!]+)\([^\n]*\)\s*(?:where\s+[^=]+)?=(?!=)') +FIELD = re.compile(r'^\s*(\w+)\s*::') +LEGACY_KEYWORDS = {'printEvery','normD','maxK','requested_maxK','epsTol','maxIter','checkEvery','bootSamples','estimSupport'} +FACTORIES = {'IndBallLinf','Constrained','_AtomIterate','CDLeastSquaresLoss','CDWeightedLSLoss','CDSqrtLassoLoss','CDSqrtQuadraticLoss','CDQuadraticLoss'} +TYPE = re.compile(r'^\s*(?:mutable\s+)?(?:struct|abstract\s+type|primitive\s+type)\s+(\w+)') + +def check(root): + errors=[] + types={m.group(1) for path in (root/'src').rglob('*.jl') + for line in path.read_text().splitlines() if (m := TYPE.match(line))} + types |= {m.group(1) for path in (root/'src').rglob('*.jl') for m in re.finditer(r'^\s*const\s+(\w+)\s*=',path.read_text(),re.M)} + for path in sorted((root/'src').rglob('*.jl')): + if 'deprecated' in path.stem: continue + in_doc=False; in_struct=False + lines=path.read_text().splitlines() + for number,line in enumerate(lines,1): + if line.count('"""') % 2: in_doc=not in_doc; continue + if in_doc: continue + line=line.split('#',1)[0] + typ=TYPE.match(line) + if typ: + name=typ.group(1) + if not name.lstrip('_')[0].isupper():errors.append((path,number,'type',name)) + in_struct='struct' in line + continue + if re.match(r'^\s*end\b',line):in_struct=False + function=FUNCTION.match(line) + if function: + name=function.group(1) or function.group(2) + # Qualified names implement external interfaces; constructors must name an actual type. + if '.' not in name and name not in types and name not in FACTORIES and name not in MATHEMATICAL and re.search('[A-Z]',name): + errors.append((path,number,'function',name)) + if ';' in line or line.lstrip().startswith('function '): + signature=line[line.find('(')+1:] + depth=1 + signature.count('(')-signature.count(')') + offset=number + while depth>0 and offset ("ProximalBase.jl", "ProximalBase.jl"), + "ProximalSolvers" => ("ProximalOPT.jl", "ProximalSolvers.jl"), + "HighDimensionalModels" => ("HighDimensionalModels.jl", "HighDimensionalModels.jl"), + "CovarianceSelection" => ("CovarianceSelection.jl", "CovarianceSelection.jl"), + "DiffPrecTest" => ("DiffPrecTest.jl", "DiffPrecTest.jl"), + "KernelSmoothing" => ("KernelSmoothing.jl", "KernelSmoothing.jl"), + "FunctionalDataAnalysis" => ("FunctionalDataAnalysis.jl", "FunctionalDataAnalysis.jl"), + "FunctionalGraphicalModels" => ("FunctionalGraphicalModels.jl", "FunctionalGraphicalModels.jl"), +) + +function instantiate_master(environment::AbstractString; + workspace::Union{Nothing,AbstractString}=nothing, + owner::AbstractString=environment) + environment = realpath(environment) + owner = realpath(owner) + workspace = workspace === nothing ? dirname(owner) : realpath(workspace) + owner_project = TOML.parsefile(joinpath(owner, "Project.toml")) + owner_name = get(owner_project, "name", "") + specs = Dict{String,Pkg.PackageSpec}() + visiting = Set{String}() + function visit(project_path) + project = TOML.parsefile(project_path) + deps = merge(get(project, "deps", Dict()), get(project, "extras", Dict())) + for name in sort!(collect(keys(deps))) + haskey(WORKSPACE_PACKAGES, name) || continue + name in visiting && continue + push!(visiting, name) + folder, repository = WORKSPACE_PACKAGES[name] + is_owner = name == owner_name + checkout = is_owner ? owner : joinpath(workspace, folder) + if !isdir(checkout) + url = "https://github.com/mlakolar/$repository" + run(`git clone --branch master --single-branch $url $checkout`) + end + branch = strip(read(`git -C $checkout branch --show-current`, String)) + (is_owner || branch == "master") || error("$name must use its master checkout: $checkout is on $branch") + dependency = TOML.parsefile(joinpath(checkout, "Project.toml")) + dependency["name"] == name || error("wrong package name at $checkout") + dependency["uuid"] == deps[name] || error("wrong UUID for $name at $checkout") + specs[name] = Pkg.PackageSpec(; path=checkout) + visit(joinpath(checkout, "Project.toml")) + end + end + project_path = joinpath(environment, "Project.toml") + original = read(project_path, String) + visit(project_path) + Pkg.activate(environment) + try + if isempty(specs) + Pkg.resolve() + else + Pkg.develop([specs[name] for name in sort!(collect(keys(specs)))]) + end + Pkg.instantiate() + finally + # Transitive master dependencies must not become stale direct dependencies. + write(project_path, original) + end + # Keep extra-only master paths in the manifest for Julia 1.10 test sandboxes. + # Resolving again after restoring Project.toml would prune those entries. + return nothing +end + +if abspath(PROGRAM_FILE) == @__FILE__ + isempty(ARGS) && error("usage: julia instantiate_master.jl PACKAGE_DIRECTORY [WORKSPACE_DIRECTORY]") + instantiate_master(ARGS[1]; workspace=length(ARGS) > 1 ? ARGS[2] : dirname(realpath(ARGS[1]))) +end diff --git a/scripts/test_master.jl b/scripts/test_master.jl new file mode 100644 index 0000000..464cd38 --- /dev/null +++ b/scripts/test_master.jl @@ -0,0 +1,53 @@ +# A dedicated test environment keeps extra-only master dependencies reachable on +# Julia 1.10 and leaves the package Project.toml unchanged for Aqua checks. +using Pkg, TOML +include(joinpath(@__DIR__, "master_dependencies.jl")) +package_root = dirname(@__DIR__) +project_path = joinpath(package_root, "Project.toml") +original_project = read(project_path, String) +project = TOML.parse(original_project) +deps = copy(get(project, "deps", Dict{String,Any}())) +extras = merge(get(project, "weakdeps", Dict()), get(project, "extras", Dict())) +for name in get(get(project, "targets", Dict()), "test", String[]) + deps[name] = extras[name] +end +deps[project["name"]] = project["uuid"] +compat = Dict(name => value for (name, value) in get(project, "compat", Dict()) + if name == "julia" || haskey(deps, name)) +sources = Dict{String,Any}() +for name in keys(deps) + if name == project["name"] + sources[name] = Dict("path" => package_root) + elseif haskey(WORKSPACE_PACKAGES, name) + folder, _ = WORKSPACE_PACKAGES[name] + sources[name] = Dict("path" => joinpath(dirname(package_root), folder)) + elseif haskey(get(project, "sources", Dict()), name) + source = copy(project["sources"][name]) + haskey(source, "path") && (source["path"] = abspath(joinpath(package_root, source["path"]))) + sources[name] = source + end +end +environment = joinpath(package_root, ".master-test-env", "v$(VERSION.major).$(VERSION.minor)") +mkpath(environment) +open(joinpath(environment, "Project.toml"), "w") do io + TOML.print(io, Dict("deps" => deps, "compat" => compat, "sources" => sources)) +end +instantiate_master(environment; workspace=dirname(package_root), owner=package_root) +read(project_path, String) == original_project || error("test setup changed package Project.toml") +installed = Pkg.dependencies() +for (name, uuid) in deps + if name == project["name"] || haskey(WORKSPACE_PACKAGES, name) + expected = name == project["name"] ? package_root : + joinpath(dirname(package_root), WORKSPACE_PACKAGES[name][1]) + realpath(installed[Base.UUID(uuid)].source) == realpath(expected) || + error("test environment did not retain the master checkout of $name") + end +end +println("Master test environment: ", environment) +if !("--setup-only" in ARGS) + depwarn = get(ENV, "JULIA_TEST_DEPWARN", "yes") + depwarn in ("yes", "no", "error") || error("invalid JULIA_TEST_DEPWARN") + coverage = get(ENV, "JULIA_TEST_COVERAGE", "false") == "true" ? "user" : "none" + test_file = joinpath(package_root, "test", "runtests.jl") + run(`$(Base.julia_cmd()) --startup-file=no --check-bounds=yes --depwarn=$depwarn --code-coverage=$coverage --project=$environment $test_file`) +end diff --git a/src/FunctionalDataAnalysis.jl b/src/FunctionalDataAnalysis.jl index 3fe9bd7..3316aab 100644 --- a/src/FunctionalDataAnalysis.jl +++ b/src/FunctionalDataAnalysis.jl @@ -14,7 +14,7 @@ import Random export AbstractFunctionalData, AbstractBasis, nobs, domain, evaluate, sampling_design, quadrature_weights -export nbasis, basismatrix, gram, penalty, BSplineBasis, FourierBasis, +export nbasis, basis_matrix, gram, penalty, BSplineBasis, FourierBasis, MonomialBasis, DerivativeBasis, derivative, FPCADerivative, MultivariateFPCADerivative export GridFunctionalData, IrregularFunctionalData, BasisFunctionalData, MultivariateFunctionalData, PiecewiseLinear, NoInterpolation, nobs_per_curve, @@ -31,8 +31,8 @@ export ids, RegularSampling, SparseSampling, simulate_kl, fourier_eigenfunctions, legendre_eigenfunctions export curve_data, mean_data, modes_of_variation, scree_data, covariance_data, design_data, fitted_curve_data -export modesplot, screeplot, covplot, designplot, fitplot -export FittedBands, MeanBands, fitted_bands, mean_bands, band_data, bandplot +export modes_plot, scree_plot, cov_plot, design_plot, fit_plot +export FittedBands, MeanBands, fitted_bands, mean_bands, band_data, band_plot export datasets, load_dataset include("interface.jl") @@ -66,5 +66,6 @@ include("derivatives.jl") include("bands.jl") include("plot_data.jl") include("recipes.jl") +include("compat_naming.jl") end diff --git a/src/bands.jl b/src/bands.jl index 2fdea63..034db6d 100644 --- a/src/bands.jl +++ b/src/bands.jl @@ -130,7 +130,7 @@ function fitted_bands(f::FPCA; alpha=0.05, coverage=:interval, grid=f.grid, curv grid isa AbstractVector{<:Real} || throw(ArgumentError("fitted bands grid must be a real vector")) index = _plot_curves(f, curves) recovered = reconstruct(f, f.scores[index, :]; grid) - functions = basismatrix(f.basis, recovered.grid) + functions = basis_matrix(f.basis, recovered.grid) variance = similar(recovered.values) for (row, i) in enumerate(index) omega = Symmetric(f.score_covariances[:, :, i]) diff --git a/src/bases.jl b/src/bases.jl index 816a257..e4beb2e 100644 --- a/src/bases.jl +++ b/src/bases.jl @@ -6,12 +6,12 @@ Number of functions in the basis `b`. function nbasis end """ - basismatrix(b, t; derivative=0) + basis_matrix(b, t; derivative=0) Evaluate a basis, or its nonnegative integer derivative, at `t`. Return a `length(t) × nbasis(b)` matrix. All points must be finite and in `domain(b)`. """ -function basismatrix end +function basis_matrix end """ gram(b) @@ -180,7 +180,7 @@ end Base.show(io::IO, ::MIME"text/plain", b::Union{BSplineBasis,FourierBasis,MonomialBasis}) = show(io, b) -function basismatrix(b::BSplineBasis, t::AbstractVector{<:Real}; derivative=0) +function basis_matrix(b::BSplineBasis, t::AbstractVector{<:Real}; derivative=0) m = _derivative_order(derivative) x = _basispoints(b, t) out = zeros(eltype(x), length(x), nbasis(b)) @@ -191,7 +191,7 @@ function basismatrix(b::BSplineBasis, t::AbstractVector{<:Real}; derivative=0) out end -function basismatrix(b::FourierBasis, t::AbstractVector{<:Real}; derivative=0) +function basis_matrix(b::FourierBasis, t::AbstractVector{<:Real}; derivative=0) m = _derivative_order(derivative) x = _basispoints(b, t) T = eltype(x) @@ -219,7 +219,7 @@ function _falling(p, m, ::Type{T}) where {T} value end -function basismatrix(b::MonomialBasis, t::AbstractVector{<:Real}; derivative=0) +function basis_matrix(b::MonomialBasis, t::AbstractVector{<:Real}; derivative=0) m = _derivative_order(derivative) x = _basispoints(b, t) T = eltype(x) diff --git a/src/compat_naming.jl b/src/compat_naming.jl new file mode 100644 index 0000000..e6036b2 --- /dev/null +++ b/src/compat_naming.jl @@ -0,0 +1,8 @@ +# Keep each old spelling bound to the canonical generic during the transition. +Base.@deprecate_binding basismatrix basis_matrix +Base.@deprecate_binding modesplot modes_plot +Base.@deprecate_binding screeplot scree_plot +Base.@deprecate_binding covplot cov_plot +Base.@deprecate_binding designplot design_plot +Base.@deprecate_binding fitplot fit_plot +Base.@deprecate_binding bandplot band_plot diff --git a/src/conversions.jl b/src/conversions.jl index 345a876..a4ee9ea 100644 --- a/src/conversions.jl +++ b/src/conversions.jl @@ -93,7 +93,7 @@ function to_basis(fd::Union{GridFunctionalData,IrregularFunctionalData}, length(times) >= k || lambda > 0 || throw(ArgumentError("curve $i has fewer observations than basis functions; lambda must be positive")) isempty(times) && throw(ArgumentError("curve $i has no observed values to fit")) - X = Matrix{T}(basismatrix(basis, times)) + X = Matrix{T}(basis_matrix(basis, times)) y = Vector{T}(values) if lambda == 0 factor = qr(X, ColumnNorm()) diff --git a/src/covariance.jl b/src/covariance.jl index 1cff709..9f8fa57 100644 --- a/src/covariance.jl +++ b/src/covariance.jl @@ -172,7 +172,7 @@ function _covariance_selection(selector::CurveCV, fd, pool, kernel) # Leaving out a fold of curves is leaving out a group of pairs. The within-cell # scatter of binned products adds the same constant to every candidate. hs = Float64.(selector.bandwidths) - cv = KernelSmoothing.selectbandwidth(KernelSmoothing.CrossValidation(), (pool.x, pool.y), pool.z; + cv = KernelSmoothing.select_bandwidth(KernelSmoothing.CrossValidation(), (pool.x, pool.y), pool.z; kernel, weights=pool.weights, groups=assignments[pool.subjects], bandwidths=hs) scatter = sum(pool.scatter) / sum(pool.weights) T = promote_type(eltype(pool.z), eltype(selector.bandwidths)) diff --git a/src/cross_structure.jl b/src/cross_structure.jl index eaac325..2abf340 100644 --- a/src/cross_structure.jl +++ b/src/cross_structure.jl @@ -1,5 +1,5 @@ """ - fit(FSVD, fd1, fd2; ncomponents=FVE(0.9999), maxK=nothing, + fit(FSVD, fd1, fd2; ncomponents=FVE(0.9999), max_k=nothing, maxK=nothing, bandwidth=nothing, kernel=KernelSmoothing.GaussianKernel(), fpca_options1=(;), fpca_options2=(;)) @@ -10,7 +10,7 @@ functions, scores, cumulative squared cross-covariance fractions, and settings are stored. Inputs are copied. Defaults follow fdapace 0.6.0 `FSVD` / `SetSVDOptions`: FVE 0.9999, -`maxK=min(20,n-2)`, Gaussian smoothing, scores enabled, and independent +`max_k=min(20,n-2)`, Gaussian smoothing, scores enabled, and independent cross-node measurement errors (equal-time products retained). Fixed component counts disable FVE selection. At least three subjects are required. Marginal FPCA options are named tuples passed to [`FPCA`](@ref); empty tuples use its @@ -106,8 +106,9 @@ end Base.@constprop :aggressive function fit(::Type{FSVD}, fd1::Union{GridFunctionalData,IrregularFunctionalData}, fd2::Union{GridFunctionalData,IrregularFunctionalData}; - ncomponents=FVE(0.9999), maxK=nothing, bandwidth=nothing, + ncomponents=FVE(0.9999), max_k=nothing, maxK=nothing, bandwidth=nothing, kernel=KernelSmoothing.GaussianKernel(), fpca_options1=(;), fpca_options2=(;)) + max_k = _max_k_keyword(max_k, maxK) _cross_subjects(fd1, fd2) n = nobs(fd1) n >= 3 || throw(ArgumentError("FSVD requires at least three subjects for its n-2 divisor and component cap")) @@ -153,19 +154,19 @@ Base.@constprop :aggressive function fit(::Type{FSVD}, fd1::Union{GridFunctional c = fit(SmoothedCrossCovariance, fd1, fd2, f1.mean, f2.mean; grid1=f1.grid, grid2=f2.grid, bandwidth, kernel) end - maxK === nothing && (maxK = min(20, n-2)) - maxK isa Integer && !(maxK isa Bool) && 1 <= maxK <= typemax(Int) || - throw(ArgumentError("FSVD maxK must be a positive integer representable as Int")) + max_k === nothing && (max_k = min(20, n-2)) + max_k isa Integer && !(max_k isa Bool) && 1 <= max_k <= typemax(Int) || + throw(ArgumentError("FSVD max_k must be a positive integer representable as Int")) root1, root2 = sqrt.(c.weights1), sqrt.(c.weights2) weighted = c.surface .* (root1 * transpose(root2)) all(isfinite, weighted) || throw(ArgumentError("FSVD weighted cross-covariance must be finite")) d = svd(weighted) - retained = min(Int(maxK), count(>(first(d.S)*_FPCA_EIGEN_TOL), d.S)) + retained = min(Int(max_k), count(>(first(d.S)*_FPCA_EIGEN_TOL), d.S)) retained > 0 || throw(ArgumentError("FSVD cross-covariance has no positive singular values above tolerance")) values = d.S[1:retained] fractions = cumsum((values ./ first(values)).^2) ./ sum((values ./ first(values)).^2) fractions[end] = one(eltype(fractions)) - selection = _fpca_selection(ncomponents, fractions, Int(maxK)) + selection = _fpca_selection(ncomponents, fractions, Int(max_k)) k = selection.chosen values = values[1:k] left, right = d.U[:,1:k] ./ root1, d.V[:,1:k] ./ root2 diff --git a/src/derivatives.jl b/src/derivatives.jl index 162b47c..bee8f71 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -3,7 +3,7 @@ The functions obtained by differentiating `basis` `order` times. The basis is copied; its domain and number of functions are retained. Evaluation delegates -to `basismatrix(basis, t; derivative=order+derivative)`. Gram and penalty matrices +to `basis_matrix(basis, t; derivative=order+derivative)`. Gram and penalty matrices are exact whenever the underlying basis supplies exact derivative penalties (B-splines through BSplineKit Galerkin, Fourier and monomial by closed forms). Zero functions and a singular Gram matrix are valid. Orders are nonnegative @@ -25,8 +25,8 @@ DerivativeBasis(basis::DerivativeBasis, order=1) = Base.eltype(::Type{DerivativeBasis{B}}) where {B} = eltype(B) nbasis(b::DerivativeBasis) = nbasis(b.basis) domain(b::DerivativeBasis) = domain(b.basis) -basismatrix(b::DerivativeBasis, t::AbstractVector{<:Real}; derivative=0) = - basismatrix(b.basis, t; derivative=_derivative_sum(b.order, derivative)) +basis_matrix(b::DerivativeBasis, t::AbstractVector{<:Real}; derivative=0) = + basis_matrix(b.basis, t; derivative=_derivative_sum(b.order, derivative)) penalty(b::DerivativeBasis, order) = penalty(b.basis, _derivative_sum(b.order, order)) Base.:(==)(a::DerivativeBasis, b::DerivativeBasis) = a.order == b.order && a.basis == b.basis Base.isequal(a::DerivativeBasis, b::DerivativeBasis) = @@ -150,7 +150,7 @@ end function _differentiate_columns(grid, values, method) out = similar(values, promote_type(eltype(values), eltype(grid), typeof(method.bandwidth))) # Every column is observed on the same grid, so one smoother matrix serves all. - L = KernelSmoothing.smoothermatrix(grid, grid; bandwidth=method.bandwidth, + L = KernelSmoothing.smoother_matrix(grid, grid; bandwidth=method.bandwidth, degree=method.degree, derivative=method.derivative, kernel=method.kernel) for (point, row) in zip(grid, eachrow(L)) all(isfinite, row) || throw(ArgumentError( @@ -260,7 +260,7 @@ function reconstruct(d::FPCADerivative, projected::AbstractMatrix{<:Real}=scores size(projected, 2) == ncomponents(d) || throw(DimensionMismatch("derivative score columns must equal ncomponents")) all(isfinite, projected) || throw(ArgumentError("derivative reconstruction scores must be finite")) quadrature_weights(grid) - functions = basismatrix(d.basis, grid) + functions = basis_matrix(d.basis, grid) mu = evaluate(d.mean, grid, PiecewiseLinear()) GridFunctionalData(projected * transpose(functions) .+ mu, grid; ids=projected === d.scores ? d.ids : nothing, diff --git a/src/evaluation.jl b/src/evaluation.jl index 996fb93..1bdcbae 100644 --- a/src/evaluation.jl +++ b/src/evaluation.jl @@ -27,7 +27,7 @@ function _evaluation_points(fd, t::AbstractVector{S}) where {S<:Real} end function evaluate(fd::BasisFunctionalData, t::AbstractVector{<:Real}; derivative=0) - fd.coefs * transpose(basismatrix(fd.basis, t; derivative)) + fd.coefs * transpose(basis_matrix(fd.basis, t; derivative)) end (fd::BasisFunctionalData)(t::AbstractVector{<:Real}; derivative=0) = evaluate(fd, t; derivative) diff --git a/src/face.jl b/src/face.jl index b8ad73e..2f49990 100644 --- a/src/face.jl +++ b/src/face.jl @@ -91,7 +91,7 @@ function _face_curve_ce(times, residuals, phi, values, grid, projected, projected_scores, (omega+transpose(omega))/2 end -Base.@constprop :aggressive function _fit_face(fd::GridFunctionalData{T}; ncomponents=FVE(0.99), maxK=nothing, +Base.@constprop :aggressive function _fit_face(fd::GridFunctionalData{T}; ncomponents=FVE(0.99), max_k=nothing, mean=nothing, eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL, knots=35, p=3, m=2, lambda=nothing, alpha=1, center=true, loglambda=collect(range(-20.0, 20.0; length=100)), mean_smoother=nothing, @@ -118,8 +118,8 @@ Base.@constprop :aggressive function _fit_face(fd::GridFunctionalData{T}; ncompo setting = _pspline_setting(fd.grid, knots, p, m) c = size(setting.A, 2) cap = min(n, c) - maxK = maxK === nothing ? cap : maxK - maxK isa Integer && !(maxK isa Bool) && 1 <= maxK <= cap || throw(ArgumentError("FACE maxK must be between 1 and min(nobs, knots+p)")) + max_k = max_k === nothing ? cap : max_k + max_k isa Integer && !(max_k isa Bool) && 1 <= max_k <= cap || throw(ArgumentError("FACE max_k must be between 1 and min(nobs, knots+p)")) mu, mean_settings = _face_mean(fd, mean, center, mean_smoother) Y, initialization = _face_initialize(fd, mu) all(isfinite, Y) || throw(ArgumentError("FACE centered observations must be finite; rescale the observations")) @@ -148,8 +148,8 @@ Base.@constprop :aggressive function _fit_face(fd::GridFunctionalData{T}; ncompo isfinite(sigma2) || throw(ArgumentError("FACE noise estimate must be finite; rescale the observations")) if observed !== nothing system = _fpca_eigen(smoothed*transpose(smoothed)/n, ones(T, c), zeros(T, c); - maxK, eigen_tol, sign_tol) - selection = _fpca_selection(ncomponents, system.fve, Int(maxK)) + max_k, eigen_tol, sign_tol) + selection = _fpca_selection(ncomponents, system.fve, Int(max_k)) U = system.eigenfunctions[:, 1:selection.chosen] d = system.eigenvalues[1:selection.chosen] completed = transpose(setting.A * (U * ((d ./ (d .+ sigma2)) .* (transpose(U)*transformed)))) @@ -169,15 +169,15 @@ Base.@constprop :aggressive function _fit_face(fd::GridFunctionalData{T}; ncompo R = Matrix(qr(sqrt.(fd.weights) .* setting.A).R) small = R * (smoothed*transpose(smoothed)/n) * transpose(R) small = (small+transpose(small))/2 - system = _fpca_eigen(small, ones(T, c), zeros(T, c); maxK, eigen_tol, sign_tol) + system = _fpca_eigen(small, ones(T, c), zeros(T, c); max_k, eigen_tol, sign_tol) functions = setting.A * (R \ system.eigenfunctions) _fpca_orient!(functions, mu, sign_tol) - selection = _fpca_selection(ncomponents, system.fve, Int(maxK)) + selection = _fpca_selection(ncomponents, system.fve, Int(max_k)) k = selection.chosen basis = FPCBasis(fd.grid, functions[:, 1:k]; weights=fd.weights) values = system.eigenvalues[1:k] # H' = Q*R implies H*H' = R'*R. Keep at most c columns, including all - # smoothed covariance directions regardless of maxK or eigen_tol. + # smoothed covariance directions regardless of max_k or eigen_tol. covariance_root = Matrix(qr(transpose(smoothed)).R) covariance = _LowRankCovariance(setting.A * transpose(covariance_root) / sqrt(T(n))) projected = _LowRankCovariance(functions .* transpose(sqrt.(system.eigenvalues))) @@ -187,7 +187,7 @@ Base.@constprop :aggressive function _fit_face(fd::GridFunctionalData{T}; ncompo scored = _pace_ce(completed_fd, fd.grid, mu, basis.values, values, projected, sigma2; curve_ce=_face_curve_ce) settings = (; requested_method=:face, grid_source=:data, ncomponents=deepcopy(ncomponents), - maxK=Int(maxK), eigen_tol, sign_tol, knots=Int(knots), p=Int(p), m=Int(m), + max_k=Int(max_k), maxK=Int(max_k), eigen_tol, sign_tol, knots=Int(knots), p=Int(p), m=Int(m), knot_sequence=copy(setting.knots), lambda=selected, lambda_source=lambda === nothing ? :GCV : :supplied, loglambda=T.(loglambda), gcv_lambdas=candidates, gcv_scores, alpha=T(alpha), center, mean=mean_settings, covariance_denominator=n, diff --git a/src/face_sparse.jl b/src/face_sparse.jl index 93b2e48..e1ad8d3 100644 --- a/src/face_sparse.jl +++ b/src/face_sparse.jl @@ -11,7 +11,7 @@ Base.getindex(c::_SplineCovariance, i::Int, j::Int) = dot(view(c.factor, i, :), Base.copy(c::_SplineCovariance) = _SplineCovariance(copy(c.factor), deepcopy(c.basis), copy(c.coefficient_root)) Base.Matrix(c::_SplineCovariance) = c.factor * transpose(c.factor) function _pace_observation_covariance(grid, c::_SplineCovariance, times) - root = basismatrix(c.basis, times)*c.coefficient_root + root = basis_matrix(c.basis, times)*c.coefficient_root root*transpose(root) end @@ -84,7 +84,7 @@ end _face_sparse_mean_values(mean, times) = mean.source === :zero ? zeros(eltype(times), length(times)) : mean.source === :supplied ? _pace_mean_values(mean.input, times) : - basismatrix(mean.basis, times)*mean.coefficients + basis_matrix(mean.basis, times)*mean.coefficients function _face_sparse_mean(fd, mean, center, knots, lambda_mean, mean_loglambda, grid) T = eltype(fd) @@ -127,7 +127,7 @@ function _face_sparse_mean(fd, mean, center, knots, lambda_mean, mean_loglambda, end Base.@constprop :aggressive function _fit_face_sparse(fd::IrregularFunctionalData{T}; grid=nothing, - ncomponents=FVE(.99), maxK=nothing, mean=nothing, + ncomponents=FVE(.99), max_k=nothing, mean=nothing, eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL, knots=7, p=3, m=2, lambda=nothing, center=true, loglambda=collect(range(-3.0, 10.0; length=14)), lambda_mean=nothing, @@ -154,8 +154,8 @@ Base.@constprop :aggressive function _fit_face_sparse(fd::IrregularFunctionalDat design = _pspline_design(fd.times, knots, p, m; interval=(a-T(.001)*width, b+T(.001)*width)) c = size(design.B, 2) cap = min(c, length(workgrid)) - maxK = maxK === nothing ? cap : maxK - maxK isa Integer && !(maxK isa Bool) && 1 <= maxK <= cap || throw(ArgumentError("sparse FACE maxK must be between 1 and min(knots+p, grid length)")) + max_k = max_k === nothing ? cap : max_k + max_k isa Integer && !(max_k isa Bool) && 1 <= max_k <= cap || throw(ArgumentError("sparse FACE max_k must be between 1 and min(knots+p, grid length)")) mean_settings = _face_sparse_mean(fd, mean, center, knots, lambda_mean, mean_loglambda, workgrid) mu = Vector{T}(_face_sparse_mean_values(mean_settings, workgrid)) residuals = fd.values - _face_sparse_mean_values(mean_settings, fd.times) @@ -202,7 +202,7 @@ Base.@constprop :aggressive function _fit_face_sparse(fd::IrregularFunctionalDat positive = findall(>(zero(T)), ce.values) isempty(positive) && throw(ArgumentError("sparse FACE covariance has no positive eigenvalues")) coefficient_root = ce.vectors[:,positive] .* transpose(sqrt.(ce.values[positive])) - Bnew = basismatrix(design.basis, workgrid) + Bnew = basis_matrix(design.basis, workgrid) # Transfer the stored quadrature into a small covariance eigensystem. gram = transpose(Bnew)*(weights .* Bnew) ge = eigen(Symmetric(gram)) @@ -211,7 +211,7 @@ Base.@constprop :aggressive function _fit_face_sparse(fd::IrregularFunctionalDat inverse_root = ge.vectors*Diagonal(inv.(sqrt.(ge.values)))*transpose(ge.vectors) smallroot = root*coefficient_root small = smallroot*transpose(smallroot) - system = _fpca_eigen((small+transpose(small))/2, ones(T,c), zeros(T,c); maxK, eigen_tol, sign_tol) + system = _fpca_eigen((small+transpose(small))/2, ones(T,c), zeros(T,c); max_k, eigen_tol, sign_tol) coefficients = inverse_root*system.eigenfunctions functions = Bnew*coefficients _fpca_orient!(functions, mu, sign_tol) @@ -219,12 +219,12 @@ Base.@constprop :aggressive function _fit_face_sparse(fd::IrregularFunctionalDat for j in axes(functions,2) dot(functions[:,j], Bnew*coefficients[:,j]) < 0 && (coefficients[:,j] .*= -1) end - selection = _fpca_selection(ncomponents, system.fve, Int(maxK)) + selection = _fpca_selection(ncomponents, system.fve, Int(max_k)) k = selection.chosen basis = FPCBasis(workgrid, functions[:,1:k]; weights) covariance = _SplineCovariance(Bnew*coefficient_root, design.basis, coefficient_root) settings = (; requested_method=:face, route=:sparse, grid_source, - ncomponents=deepcopy(ncomponents), maxK=Int(maxK), eigen_tol, sign_tol, + ncomponents=deepcopy(ncomponents), max_k=Int(max_k), maxK=Int(max_k), eigen_tol, sign_tol, knots=Int(knots), p=Int(p), m=Int(m), knot_sequence=design.knots, knot_range_expansion=T(.001), lambda=selected, lambda_source=lambda === nothing ? :subject_CV : :supplied, @@ -255,7 +255,7 @@ end function _face_sparse_scores(fd, grid, mu, functions, values, covariance, sigma2, settings) _pace_ce(fd, grid, mu, functions, values, covariance, sigma2; curve_ce=_face_sparse_curve_ce, - basis_evaluator=(g,f,t)->basismatrix(covariance.basis,t)*settings.eigenfunction_coefficients, + basis_evaluator=(g,f,t)->basis_matrix(covariance.basis,t)*settings.eigenfunction_coefficients, mean_evaluator=(g,m,t)->_face_sparse_mean_values(settings.mean,t)) end diff --git a/src/fpca.jl b/src/fpca.jl index 930b054..b7ec03e 100644 --- a/src/fpca.jl +++ b/src/fpca.jl @@ -1,5 +1,14 @@ abstract type _ComponentSelector end +function _max_k_keyword(max_k, maxK) + if maxK !== nothing + max_k === nothing || throw(ArgumentError("specify only max_k, not both keyword spellings")) + Base.depwarn("keyword maxK is deprecated; use max_k", :maxK) + return maxK + end + return max_k +end + const _FPCA_EIGEN_TOL = 1e-10 const _FPCA_SIGN_TOL = 1e-12 @@ -8,7 +17,7 @@ const _FPCA_SIGN_TOL = 1e-12 Select the smallest number of components whose cumulative fraction of variance explained is at least `threshold`, a finite number in `(0, 1]`. The denominator -is the sum of positive eigenvalues retained up to `maxK`. The default 0.99 is +is the sum of positive eigenvalues retained up to `max_k`. The default 0.99 is fdapace 0.6.0's `FVEthreshold` (`SetOptions.R`). """ struct FVE{T<:AbstractFloat} <: _ComponentSelector @@ -48,7 +57,7 @@ Base.show(io::IO, rule::FVE) = print(io, "FVE(", rule.threshold, ")") """ fit(FPCA, fd; method=:auto, grid=nothing, ncomponents=FVE(0.99), - maxK=nothing, presmooth=nothing, mean=nothing, covariance=nothing, + max_k=nothing, presmooth=nothing, mean=nothing, covariance=nothing, noise_variance=nothing, mean_bandwidth=nothing, covariance_bandwidth=nothing, diagonal_bandwidth=nothing, kernel=nothing, weighting=:obs, binning=:auto, face=(;), eigen_tol=1e-10, sign_tol=1e-12) @@ -76,7 +85,7 @@ the simple knot sequence extends `p` intervals beyond each endpoint. designs throw instead of refund's undocumented ridge. `lambda=nothing` minimizes refund's closed-form GCV over `exp.(loglambda)`; ties use the first candidate. A finite nonnegative fixed lambda is also accepted. The full -candidate grid and criterion are stored even for a fixed lambda. `maxK` +candidate grid and criterion are stored even for a fixed lambda. `max_k` defaults to `min(nobs, knots+p)`; FVE uses the retained positive spectrum. Fixed K, FVE, numerical eigenvalue cutoff and sign rule reuse the other routes; AIC and BIC throw. Refund uses strict `FVE > threshold`, while the package @@ -120,7 +129,7 @@ intervals (six strictly interior knots), `p=3`, `m=2`, `center=true`, `lambda=nothing`, and `loglambda=range(-3,10; length=14)`. Knots span the observed range expanded by 0.1% at each end. `grid=nothing` uses 100 equally spaced points across the observed range; an override must cover all observations -and stay within that range, with full spline rank. `maxK` defaults to `knots+p`. +and stay within that range, with full spline rank. `max_k` defaults to `knots+p`. The reference's `gcv` approximates leave-one-subject-out squared prediction error to first order, omitting the constant sum of squared raw products. It fits the covariance and noise jointly, including diagonal raw products. @@ -195,16 +204,16 @@ count are stored. With `P=sum(N_i*(N_i-1))` pairs and `M=length(grid)`, fixed covariance smoothing costs `O(M^2*P)` time. A 2D `CurveCV` with `H` candidates costs `O(H*P^2 + H*F*P)` for `F` folds, before the final fit; binning reduces `P` to the number of occupied per-curve cells (at most `n*binning^2`). -`maxK=nothing` resolves to `max(1,min(20,length(grid)-2,nobs(fd)-2))` for +`max_k=nothing` resolves to `max(1,min(20,length(grid)-2,nobs(fd)-2))` for PACE, following fdapace 0.6.0's smoothed route. An explicit PACE cap may be at most `length(grid)`. -The full positive-semidefinite covariance projection up to `maxK` is used +The full positive-semidefinite covariance projection up to `max_k` is used in each curve's Cholesky score solve. The selected eigenvalues and eigenfunctions define the score numerator and conditional covariance. PACE also supports [`AIC`](@ref) and [`BIC`](@ref). Their criteria use the K-component covariance, including the Gaussian normalization constant, for -every candidate `K=1:maxK`. Beyond the retained positive rank the covariance +every candidate `K=1:max_k`. Beyond the retained positive rank the covariance has zero additional variance, while the penalty still increases. Selection stops at the first increase as fdapace does; the global minimizer is retained too. No noise floor or ridge is added. A non-positive-definite observation covariance @@ -215,12 +224,12 @@ including on nonuniform grids. Eigenvalues must exceed `eigen_tol` times the largest absolute eigenvalue. The relative cutoff `1e-10` is a package numerical choice, adjustable here. All negative eigenvalues, including roundoff, are recorded in `dropped` as their `number`, absolute `total`, and `share` of the -absolute spectrum. Positive eigenvalues up to `maxK` define the FVE denominator. +absolute spectrum. Positive eigenvalues up to `max_k` define the FVE denominator. `ncomponents` is a positive integer or [`FVE`](@ref); the default threshold 0.99 -is fdapace's. The default `maxK=max(1, min(length(fd.grid)-2, nobs(fd)-2))` +is fdapace's. The default `max_k=max(1, min(length(fd.grid)-2, nobs(fd)-2))` follows fdapace 0.6.0's dense cross-sectional default (`SetOptions.R`), without -the cap of 20 used on its smoothing route. An explicit dense `maxK` cannot -exceed `min(length(fd.grid), nobs(fd)-1)`. A fixed K exceeding `maxK` +the cap of 20 used on its smoothing route. An explicit dense `max_k` cannot +exceed `min(length(fd.grid), nobs(fd)-1)`. A fixed K exceeding `max_k` or the retained positive rank throws. A zero-variance sample also throws. Signs make `sum(phi .* mean)` nonnegative, following fdapace. When its magnitude @@ -273,14 +282,14 @@ end # rectangular eigenvalue scaling, followed by trapezoidal column normalization. # Both weights coincide in every public fit. function _fpca_eigen(surface::AbstractMatrix{<:Real}, weights::AbstractVector{<:Real}, - mean::AbstractVector{<:Real}; maxK=length(weights), eigen_tol=_FPCA_EIGEN_TOL, + mean::AbstractVector{<:Real}; max_k=length(weights), eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL, normalization_weights=weights) Base.require_one_based_indexing(surface, weights, mean, normalization_weights) m = length(weights) size(surface) == (m, m) || throw(DimensionMismatch("covariance surface must be square with one row per weight")) length(mean) == m || throw(DimensionMismatch("mean length must match covariance grid")) length(normalization_weights) == m || throw(DimensionMismatch("normalization weights must match covariance grid")) - maxK isa Integer && 1 <= maxK <= m || throw(ArgumentError("maxK must be between 1 and the covariance grid length")) + max_k isa Integer && 1 <= max_k <= m || throw(ArgumentError("max_k must be between 1 and the covariance grid length")) all(isfinite, surface) && all(isfinite, mean) || throw(ArgumentError("covariance surface and mean must be finite")) issymmetric(surface) || throw(ArgumentError("covariance surface must be symmetric")) all(w -> isfinite(w) && w > 0, weights) && @@ -308,7 +317,7 @@ function _fpca_eigen(surface::AbstractMatrix{<:Real}, weights::AbstractVector{<: cutoff = eigen_tol * scale indices = reverse(findall(>(cutoff), spectrum)) isempty(indices) && throw(ArgumentError("FPCA covariance has no positive eigenvalues above eigen_tol")) - indices = indices[1:min(Int(maxK), length(indices))] + indices = indices[1:min(Int(max_k), length(indices))] values = spectrum[indices] functions = decomposition.vectors[:, indices] ./ root normalization_root = sqrt.(normalization_weights) @@ -339,20 +348,20 @@ function _fpca_orient!(functions, mean, sign_tol) functions end -function _fpca_selection(rule::Integer, curve, maxK) - 1 <= rule <= maxK || throw(ArgumentError("ncomponents must be between 1 and maxK ($maxK)")) +function _fpca_selection(rule::Integer, curve, max_k) + 1 <= rule <= max_k || throw(ArgumentError("ncomponents must be between 1 and max_k ($max_k)")) rule <= length(curve) || throw(ArgumentError("ncomponents exceeds the retained positive covariance rank ($(length(curve)))")) - (; rule=:fixed, threshold=nothing, maxK, chosen=Int(rule), scores=eltype(curve)[]) + (; rule=:fixed, threshold=nothing, max_k, maxK=max_k, chosen=Int(rule), scores=eltype(curve)[]) end -function _fpca_selection(rule::FVE, curve, maxK) +function _fpca_selection(rule::FVE, curve, max_k) chosen = findfirst(>=(rule.threshold), curve)::Int - (; rule=:FVE, threshold=rule.threshold, maxK, chosen, scores=copy(curve)) + (; rule=:FVE, threshold=rule.threshold, max_k, maxK=max_k, chosen, scores=copy(curve)) end -_fpca_selection(::_ComponentSelector, curve, maxK) = +_fpca_selection(::_ComponentSelector, curve, max_k) = throw(ArgumentError("AIC and BIC component selection require method=:pace")) -_fpca_selection(rule, curve, maxK) = +_fpca_selection(rule, curve, max_k) = throw(ArgumentError("ncomponents must be a positive integer, FVE(threshold), or (for PACE) AIC() or BIC()")) _fpca_preprocess(fd, ::Nothing) = (copy(fd.values), nothing) @@ -369,11 +378,12 @@ _fpca_preprocess(fd, method) = throw(ArgumentError("presmooth must be nothing, LocalPolynomial, or PenalizedSpline")) Base.@constprop :aggressive function fit(::Type{FPCA}, fd::GridFunctionalData; method=:auto, grid=nothing, - ncomponents=FVE(0.99), maxK=nothing, presmooth=nothing, + ncomponents=FVE(0.99), max_k=nothing, maxK=nothing, presmooth=nothing, mean=nothing, covariance=nothing, noise_variance=nothing, mean_bandwidth=nothing, covariance_bandwidth=nothing, diagonal_bandwidth=nothing, kernel=nothing, weighting=:obs, binning=:auto, face=(;), eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL) + max_k = _max_k_keyword(max_k, maxK) method in (:auto, :dense, :pace, :face) || throw(ArgumentError("FPCA method must be :auto, :dense, :pace, or :face")) face isa NamedTuple || throw(ArgumentError("face must be a named tuple of FACE settings")) if method == :face @@ -382,13 +392,13 @@ Base.@constprop :aggressive function fit(::Type{FPCA}, fd::GridFunctionalData; m covariance_bandwidth === nothing && diagonal_bandwidth === nothing && kernel === nothing && weighting == :obs && binning === :auto) || throw(ArgumentError("FACE uses the data grid and its own smoothing and noise estimates; grid, presmooth and PACE tuning are unsupported")) - return _fit_face(fd; ncomponents, maxK, mean, eigen_tol, sign_tol, face...) + return _fit_face(fd; ncomponents, max_k, mean, eigen_tol, sign_tol, face...) end isempty(face) || throw(ArgumentError("face settings require method=:face")) design = sampling_design(fd) if method == :pace || (method == :auto && design != :regular) presmooth === nothing || throw(ArgumentError("presmooth is only supported by dense FPCA; PACE smooths the pooled mean and covariance")) - return _fit_pace(fd; method, grid, ncomponents, maxK, mean, covariance, + return _fit_pace(fd; method, grid, ncomponents, max_k, mean, covariance, noise_variance, mean_bandwidth, covariance_bandwidth, diagonal_bandwidth, kernel, weighting, binning, eigen_tol, sign_tol) end @@ -400,9 +410,9 @@ Base.@constprop :aggressive function fit(::Type{FPCA}, fd::GridFunctionalData; m throw(ArgumentError("grid, mean, covariance, noise_variance and pooled smoothing options require method=:pace")) n, m = size(fd.values) n >= 2 || throw(ArgumentError("dense FPCA requires at least two curves")) - maxK = maxK === nothing ? max(1, min(m-2, n-2)) : maxK - maxK isa Integer && 1 <= maxK <= min(m, n-1) || throw(ArgumentError( - "dense FPCA maxK must be between 1 and min(grid length, nobs-1)")) + max_k = max_k === nothing ? max(1, min(m-2, n-2)) : max_k + max_k isa Integer && 1 <= max_k <= min(m, n-1) || throw(ArgumentError( + "dense FPCA max_k must be between 1 and min(grid length, nobs-1)")) _fpca_tolerance(eigen_tol, "eigen_tol") _fpca_tolerance(sign_tol, "sign_tol") values, preprocessing = _fpca_preprocess(fd, presmooth) @@ -413,8 +423,8 @@ Base.@constprop :aggressive function fit(::Type{FPCA}, fd::GridFunctionalData; m all(isfinite, covariance) && all(isfinite, mu) || throw(ArgumentError("dense FPCA mean and covariance must be finite; rescale the observations")) grid, weights = Vector{T}(fd.grid), Vector{T}(fd.weights) - system = _fpca_eigen(covariance, weights, mu; maxK, eigen_tol, sign_tol) - selection = _fpca_selection(ncomponents, system.fve, Int(maxK)) + system = _fpca_eigen(covariance, weights, mu; max_k, eigen_tol, sign_tol) + selection = _fpca_selection(ncomponents, system.fve, Int(max_k)) k = selection.chosen basis = FPCBasis(grid, system.eigenfunctions[:, 1:k]; weights) projected = residuals * (weights .* basis.values) @@ -424,7 +434,7 @@ Base.@constprop :aggressive function fit(::Type{FPCA}, fd::GridFunctionalData; m all(isfinite, full_projected_covariance) || throw(ArgumentError( "dense FPCA fitted covariance must be finite; rescale the observations")) settings = (; requested_method=method, grid_source=:data, ncomponents=deepcopy(ncomponents), - maxK=Int(maxK), eigen_tol, sign_tol, + max_k=Int(max_k), maxK=Int(max_k), eigen_tol, sign_tol, covariance_denominator=n-1, quadrature=:stored_trapezoidal, preprocessing, observations=:copied_original, full_projected_covariance) sign_rule = (; method=:mean_alignment, fallback=:largest_absolute_entry, tolerance=sign_tol) @@ -435,24 +445,25 @@ Base.@constprop :aggressive function fit(::Type{FPCA}, fd::GridFunctionalData; m end Base.@constprop :aggressive function fit(::Type{FPCA}, fd::IrregularFunctionalData; method=:auto, grid=nothing, - ncomponents=FVE(0.99), maxK=nothing, presmooth=nothing, + ncomponents=FVE(0.99), max_k=nothing, maxK=nothing, presmooth=nothing, mean=nothing, covariance=nothing, noise_variance=nothing, mean_bandwidth=nothing, covariance_bandwidth=nothing, diagonal_bandwidth=nothing, kernel=nothing, weighting=:obs, binning=:auto, face=(;), eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL) + max_k = _max_k_keyword(max_k, maxK) face isa NamedTuple || throw(ArgumentError("face must be a named tuple of FACE settings")) if method == :face (presmooth === nothing && covariance === nothing && noise_variance === nothing && mean_bandwidth === nothing && covariance_bandwidth === nothing && diagonal_bandwidth === nothing && kernel === nothing && weighting == :obs && binning === :auto) || throw(ArgumentError("sparse FACE uses its own smoothing and noise estimates; presmooth and PACE tuning are unsupported")) - return _fit_face_sparse(fd; grid, ncomponents, maxK, mean, eigen_tol, sign_tol, face...) + return _fit_face_sparse(fd; grid, ncomponents, max_k, mean, eigen_tol, sign_tol, face...) end isempty(face) || throw(ArgumentError("face settings require method=:face")) method in (:auto, :dense, :pace) || throw(ArgumentError("FPCA method must be :auto, :dense, :pace, or :face")) method == :dense && throw(ArgumentError("dense FPCA requires complete GridFunctionalData")) presmooth === nothing || throw(ArgumentError("presmooth is only supported by dense FPCA; PACE smooths the pooled mean and covariance")) - _fit_pace(fd; method, grid, ncomponents, maxK, mean, covariance, + _fit_pace(fd; method, grid, ncomponents, max_k, mean, covariance, noise_variance, mean_bandwidth, covariance_bandwidth, diagonal_bandwidth, kernel, weighting, binning, eigen_tol, sign_tol) end @@ -462,7 +473,7 @@ end Return the selected eigenfunctions as an [`FPCBasis`](@ref) for grid/PACE fits, or an exact linear combination of the input basis for basis-data fits. -Both implement `basismatrix`, `gram`, `penalty`, `domain`, and `nbasis`. +Both implement `basis_matrix`, `gram`, `penalty`, `domain`, and `nbasis`. """ eigenfunctions(f::FPCA) = f.basis @@ -496,7 +507,7 @@ score_covariances(f::FPCA) = f.score_covariances fve(f::FPCA) Return cumulative fractions of variance explained for all retained positive -eigenvalues up to `maxK`, including those beyond the selected components. +eigenvalues up to `max_k`, including those beyond the selected components. """ fve(f::FPCA) = f.fve @@ -535,7 +546,7 @@ MultivariateStats.jl's PCA projection interface. PACE accepts either representation and returns `(scores, covariances)` as a named tuple, with an `n × K` matrix and `K × K × n` array. It uses the training -mean, eigenfunctions, full projected covariance up to `maxK`, and noise +mean, eigenfunctions, full projected covariance up to `max_k`, and noise variance. Every observation time must lie inside the work grid. Predicting the training data reproduces [`scores`](@ref) and [`score_covariances`](@ref). @@ -595,7 +606,7 @@ function reconstruct(f::FPCA, projected::AbstractMatrix{<:Real}=scores(f); size(projected, 2) == ncomponents(f) || throw(DimensionMismatch("reconstruction score columns must equal ncomponents")) all(isfinite, projected) || throw(ArgumentError("reconstruction scores must be finite")) quadrature_weights(grid) - functions = basismatrix(f.basis, grid) + functions = basis_matrix(f.basis, grid) mu = _fpca_mean_values(f.mean, grid) reconstructed = projected * transpose(functions) .+ mu GridFunctionalData(reconstructed, grid; ids=projected === f.scores ? f.ids : nothing) diff --git a/src/fpcbasis.jl b/src/fpcbasis.jl index 7511642..01d6a86 100644 --- a/src/fpcbasis.jl +++ b/src/fpcbasis.jl @@ -7,7 +7,7 @@ must be finite and strictly increasing, with at least two points. At least one function is required; values must be finite and weights must be finite, positive and have the same length as the grid. -[`basismatrix`](@ref) evaluates the piecewise-linear interpolants inside the +[`basis_matrix`](@ref) evaluates the piecewise-linear interpolants inside the closed grid range and throws outside it. Only derivative zero is supported; [`penalty`](@ref) likewise supports only order zero. [`gram`](@ref) gives exact integrals of products of the piecewise-linear interpolants. @@ -64,7 +64,7 @@ end Base.show(io::IO, ::MIME"text/plain", b::FPCBasis) = show(io, b) -function basismatrix(b::FPCBasis, t::AbstractVector{<:Real}; derivative=0) +function basis_matrix(b::FPCBasis, t::AbstractVector{<:Real}; derivative=0) _derivative_order(derivative) == 0 || throw(ArgumentError("FPCBasis evaluation supports derivative=0 only")) points = _basispoints(b, t) diff --git a/src/interface.jl b/src/interface.jl index 73e81f8..78d7526 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -12,7 +12,7 @@ abstract type AbstractFunctionalData{T<:AbstractFloat} end AbstractBasis A basis for representing functions on a common domain. The interface -is `nbasis(b)`, `domain(b)`, `basismatrix(b, t; derivative)`, `gram(b)` and +is `nbasis(b)`, `domain(b)`, `basis_matrix(b, t; derivative)`, `gram(b)` and `penalty(b, order)`. """ abstract type AbstractBasis end diff --git a/src/mean.jl b/src/mean.jl index a0a9a5c..3edea28 100644 --- a/src/mean.jl +++ b/src/mean.jl @@ -63,10 +63,10 @@ _initial_bandwidth(cv::CurveCV) = cv.pilot === nothing ? first(cv.bandwidths) : function _regression_bandwidth(rule::KernelSmoothing.BandwidthSelector, t, y, w, kernel, degree) data = KernelSmoothing.SmoothingData(t, y; weights=w) - hasmethod(KernelSmoothing.selectbandwidth, + hasmethod(KernelSmoothing.select_bandwidth, Tuple{typeof(rule),typeof(data),typeof(kernel)}, (:degree,)) || throw(ArgumentError("$rule is density-only or does not support regression selection with degree; use a regression rule or CurveCV")) - selected = KernelSmoothing.selectbandwidth(rule, data, kernel; degree) + selected = KernelSmoothing.select_bandwidth(rule, data, kernel; degree) selected.h isa Real || throw(ArgumentError( "$rule returns pointwise bandwidths; supply a rule returning one global bandwidth or CurveCV")) h = LocalPolynomial(; bandwidth=selected.h, kernel, degree).bandwidth @@ -172,7 +172,7 @@ function _mean_pool(fd, weighting) end function _mean_predict(t, y, w, points, method; label="pooled mean") - values = KernelSmoothing.localpolynomial(t, y, points; bandwidth=method.bandwidth, + values = KernelSmoothing.local_polynomial(t, y, points; bandwidth=method.bandwidth, kernel=method.kernel, degree=method.degree, derivative=method.derivative, weights=w) for (point, value) in zip(points, values) isfinite(value) || throw(ArgumentError( diff --git a/src/mfpca.jl b/src/mfpca.jl index e806b99..8ab6749 100644 --- a/src/mfpca.jl +++ b/src/mfpca.jl @@ -102,7 +102,7 @@ function _mfpca(marginal::MultivariateFPCA{U}, weights::AbstractVector{W}, centered = scaled .- sum(scaled; dims=1) ./ nobs(marginal) covariance = transpose(centered) * centered system = _fpca_eigen(covariance, ones(T, total), zeros(T, total); - maxK=total, eigen_tol, sign_tol=0) + max_k=total, eigen_tol, sign_tol=0) selection = _fpca_selection(count, system.fve, total) c = system.eigenfunctions[:, 1:Int(count)] values = system.eigenvalues[1:Int(count)] @@ -226,7 +226,7 @@ function reconstruct(f::MFPCA{T}, projected::AbstractMatrix{S}=scores(f); for j in eachindex(f.eigenfunctions) points = grids[j] quadrature_weights(points) - functions = basismatrix(f.eigenfunctions[j], points) + functions = basis_matrix(f.eigenfunctions[j], points) mu = evaluate(f.univariate.means[j], points, PiecewiseLinear()) values = projected * transpose(functions) .+ mu push!(nodes, GridFunctionalData(Matrix{R}(values), Vector{R}(points); diff --git a/src/multivariate_fpca.jl b/src/multivariate_fpca.jl index 865e240..65ebeb2 100644 --- a/src/multivariate_fpca.jl +++ b/src/multivariate_fpca.jl @@ -24,11 +24,11 @@ The pooled spectrum, selection, dropped spectrum, sign rule, and settings are stored in `f.pooled`. The sign reference is the average node mean; the univariate largest-entry fallback resolves near-zero alignment. -Pooled defaults `ncomponents=FVE(0.99)`, `maxK=nothing`, `eigen_tol=1e-10`, +Pooled defaults `ncomponents=FVE(0.99)`, `max_k=nothing`, `eigen_tol=1e-10`, and `sign_tol=1e-12` follow the dense univariate defaults and have the sources given in [`FPCA`](@ref). The pooled covariance averages `p` node covariances of rank at most `nobs-1`, so its rank is at most `p*(nobs-1)`; -`maxK` therefore resolves to `max(1,min(grid_length-2,p*(nobs-1)-1))`, the +`max_k` therefore resolves to `max(1,min(grid_length-2,p*(nobs-1)-1))`, the dense rule with that rank, with explicit cap `min(grid_length,p*(nobs-1))`. Only shared values of these four keywords and `method=:auto`/`:dense` are supported for pooled fits. PACE tuning and `presmooth` throw instead of being ignored. Inputs and settings are copied. @@ -106,13 +106,14 @@ function _multivariate_marginal(mfd, options) end function _multivariate_pooled(mfd::MultivariateFunctionalData{T}; - ncomponents=FVE(0.99), maxK=nothing, eigen_tol=1e-10, sign_tol=1e-12, + ncomponents=FVE(0.99), max_k=nothing, maxK=nothing, eigen_tol=1e-10, sign_tol=1e-12, method=:auto, kwargs...) where {T} + max_k = _max_k_keyword(max_k, maxK) isempty(kwargs) || throw(ArgumentError( - "pooled FPCA supports only ncomponents, maxK, eigen_tol, sign_tol and method; PACE tuning and presmooth are unsupported")) + "pooled FPCA supports only ncomponents, max_k, eigen_tol, sign_tol and method; PACE tuning and presmooth are unsupported")) method in (:auto, :dense) || throw(ArgumentError( "pooled FPCA requires method=:auto or :dense; pooled PACE is not implemented")) - for (name, value) in pairs((; ncomponents, maxK, eigen_tol, sign_tol)) + for (name, value) in pairs((; ncomponents, max_k, eigen_tol, sign_tol)) value isa AbstractVector && throw(ArgumentError("pooled FPCA $name must be one shared value")) end all(c -> c isa GridFunctionalData && c.observed === nothing, mfd.components) || @@ -128,10 +129,10 @@ function _multivariate_pooled(mfd::MultivariateFunctionalData{T}; # The pooled covariance averages p rank-(n-1) node covariances, so its rank # is at most p(n-1); the defaults mirror the dense rule with that rank. rank = p * (n - 1) - requested_maxK = maxK - maxK = maxK === nothing ? max(1, min(g-2, rank-1)) : maxK - maxK isa Integer && 1 <= maxK <= min(g, rank) || throw(ArgumentError( - "pooled FPCA maxK must be between 1 and min(grid length, nodes*(nobs-1))")) + requested_max_k = max_k + max_k = max_k === nothing ? max(1, min(g-2, rank-1)) : max_k + max_k isa Integer && 1 <= max_k <= min(g, rank) || throw(ArgumentError( + "pooled FPCA max_k must be between 1 and min(grid length, nodes*(nobs-1))")) grid, weights = copy(firstnode.grid), copy(firstnode.weights) means = GridFunctionalData{T}[] residuals = Matrix{T}[] @@ -145,8 +146,8 @@ function _multivariate_pooled(mfd::MultivariateFunctionalData{T}; covariance .+= (transpose(centered) * centered / (n-1)) / p end reference_mean = vec(sum(mu.values ./ p for mu in means)) - system = _fpca_eigen(covariance, weights, reference_mean; maxK, eigen_tol, sign_tol) - selection = _fpca_selection(ncomponents, system.fve, Int(maxK)) + system = _fpca_eigen(covariance, weights, reference_mean; max_k, eigen_tol, sign_tol) + selection = _fpca_selection(ncomponents, system.fve, Int(max_k)) k = selection.chosen shared = FPCBasis(grid, system.eigenfunctions[:, 1:k]; weights) projected = Matrix{T}[centered * (weights .* shared.values) for centered in residuals] @@ -166,7 +167,7 @@ function _multivariate_pooled(mfd::MultivariateFunctionalData{T}; all_eigenvalues=system.eigenvalues, fve=system.fve, basis=shared, selection, dropped=system.dropped, sign_rule, reference_mean) settings = (; requested_method=method, ncomponents=deepcopy(ncomponents), - requested_maxK, maxK=Int(maxK), eigen_tol, sign_tol, + requested_max_k, requested_maxK=requested_max_k, max_k=Int(max_k), maxK=Int(max_k), eigen_tol, sign_tol, grid_source=:data, quadrature=:stored_trapezoidal, covariance_denominator=n-1, score_variance_denominator=n-1, covariance_pooling=:equal_nodes, sign_reference=:average_node_mean, @@ -364,7 +365,7 @@ function reconstruct(f::MultivariateFPCA{T}, projected::AbstractVector=scores(f) all(isfinite, s) || throw(ArgumentError("reconstruction scores at node $j must be finite")) points = grids[j] quadrature_weights(points) - functions = basismatrix(f.bases[j], points) + functions = basis_matrix(f.bases[j], points) mu = evaluate(f.means[j], points, PiecewiseLinear()) values = s * transpose(functions) .+ mu push!(nodes, GridFunctionalData(Matrix{R}(values), Vector{R}(points); diff --git a/src/pace.jl b/src/pace.jl index ce8615b..a2f464e 100644 --- a/src/pace.jl +++ b/src/pace.jl @@ -224,10 +224,10 @@ function _pace_curve_ce(times, residuals, phi, values, grid, projected::Abstract projected_scores, (omega + transpose(omega))/2 end -function _pace_likelihood_selection(rule::Union{AIC,BIC}, fd, grid, mu, system, maxK, sigma2) +function _pace_likelihood_selection(rule::Union{AIC,BIC}, fd, grid, mu, system, max_k, sigma2) values, functions = system.eigenvalues, system.eigenfunctions T = eltype(values) - candidates = collect(1:maxK) + candidates = collect(1:max_k) criterion = zeros(T, length(candidates)) for i in 1:nobs(fd) times, observations = fd[i] @@ -256,17 +256,17 @@ function _pace_likelihood_selection(rule::Union{AIC,BIC}, fd, grid, mu, system, "PACE component-selection criteria must be finite; rescale the observations")) increase = findfirst(>(zero(T)), diff(criterion)) chosen = increase === nothing ? length(candidates) : increase - (; rule=rule isa AIC ? :AIC : :BIC, threshold=nothing, maxK, + (; rule=rule isa AIC ? :AIC : :BIC, threshold=nothing, max_k, maxK=max_k, chosen, scores=criterion, candidates, global_minimizer=argmin(criterion), stopping=:first_increase) end -function _pace_selection(rule, fd, grid, mu, system, maxK, sigma2) - _fpca_selection(rule, system.fve, maxK) +function _pace_selection(rule, fd, grid, mu, system, max_k, sigma2) + _fpca_selection(rule, system.fve, max_k) end -function _pace_selection(rule::Union{AIC,BIC}, fd, grid, mu, system, maxK, sigma2) - _pace_likelihood_selection(rule, fd, grid, mu, system, maxK, sigma2) +function _pace_selection(rule::Union{AIC,BIC}, fd, grid, mu, system, max_k, sigma2) + _pace_likelihood_selection(rule, fd, grid, mu, system, max_k, sigma2) end _pace_mean_settings(mean::SmoothedMean) = (; bandwidth=mean.bandwidth, @@ -285,7 +285,7 @@ _pace_covariance_settings(covariance::SmoothedCovariance) = _pace_covariance_settings(covariance::NamedTuple) = (; interpolation=:bilinear, grid=copy(covariance.grid)) -function _fit_pace(fd; method=:auto, grid=nothing, ncomponents=FVE(0.99), maxK=nothing, +function _fit_pace(fd; method=:auto, grid=nothing, ncomponents=FVE(0.99), max_k=nothing, mean=nothing, covariance=nothing, noise_variance=nothing, mean_bandwidth=nothing, covariance_bandwidth=nothing, diagonal_bandwidth=nothing, kernel=nothing, weighting=:obs, binning=:auto, eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL) @@ -326,10 +326,10 @@ function _fit_pace(fd; method=:auto, grid=nothing, ncomponents=FVE(0.99), maxK=n covariance === nothing && covariance_bandwidth === nothing && (covariance_bandwidth = _pace_default_covariance_bandwidth(a, b)) _fpca_tolerance(eigen_tol, "eigen_tol") _fpca_tolerance(sign_tol, "sign_tol") - requested_maxK = maxK - maxK === nothing && (maxK = max(1, min(20, length(workgrid)-2, nobs(fd)-2))) - maxK isa Integer && 1 <= maxK <= length(workgrid) || - throw(ArgumentError("PACE maxK must be between 1 and the work-grid length")) + requested_max_k = max_k + max_k === nothing && (max_k = max(1, min(20, length(workgrid)-2, nobs(fd)-2))) + max_k isa Integer && 1 <= max_k <= length(workgrid) || + throw(ArgumentError("PACE max_k must be between 1 and the work-grid length")) resolved_mean, mean_source = _pace_mean(fd, mean, workgrid, mean_bandwidth, kernel, weighting) resolved_covariance, covariance_source = _pace_covariance(fd, resolved_mean, covariance, workgrid, covariance_bandwidth, diagonal_bandwidth, kernel, weighting, binning) @@ -341,15 +341,15 @@ function _fit_pace(fd; method=:auto, grid=nothing, ncomponents=FVE(0.99), maxK=n mu, surface, sigma2 = Vector{T}(mean_values), Matrix{T}(surface), T(sigma2) all(isfinite, mu) && all(isfinite, surface) && isfinite(sigma2) || throw(ArgumentError("PACE mean, covariance, and noise variance must remain finite after conversion")) - system = _fpca_eigen(surface, weights, mu; maxK, eigen_tol, sign_tol) - selection = _pace_selection(ncomponents, fd, workgrid, mu, system, Int(maxK), sigma2) + system = _fpca_eigen(surface, weights, mu; max_k, eigen_tol, sign_tol) + selection = _pace_selection(ncomponents, fd, workgrid, mu, system, Int(max_k), sigma2) k = selection.chosen projected = (system.eigenfunctions .* transpose(system.eigenvalues)) * transpose(system.eigenfunctions) projected = (projected + transpose(projected))/2 basis = FPCBasis(workgrid, system.eigenfunctions[:, 1:k]; weights) predicted = _pace_ce(fd, workgrid, mu, basis.values, system.eigenvalues[1:k], projected, sigma2) - settings = merge((; requested_method=method, grid_source, requested_maxK, maxK=Int(maxK), + settings = merge((; requested_method=method, grid_source, requested_max_k, requested_maxK=requested_max_k, max_k=Int(max_k), maxK=Int(max_k), ncomponents=deepcopy(ncomponents), observed_range=(a, b), observations=:copied_original, eigen_tol, sign_tol, quadrature=:work_grid_trapezoidal, mean_source, covariance_source, noise_variance_source=noise_source, diff --git a/src/plot_data.jl b/src/plot_data.jl index 87a7828..d4f3d7b 100644 --- a/src/plot_data.jl +++ b/src/plot_data.jl @@ -171,7 +171,7 @@ Accepts `SmoothedCovariance`, `SmoothedCrossCovariance`, or exact dense FPCA defaults to its stored estimated surface (`:stored`, a package display choice). `surface=:fitted` returns the positive-semidefinite projection over -all retained components up to `maxK`, including unselected components, as in +all retained components up to `max_k`, including unselected components, as in fdapace 0.6.0's full fitted covariance used for CE scores. It does not refit or diagonalize a surface. Other symbols throw. """ diff --git a/src/psplines.jl b/src/psplines.jl index 1b0c6b4..3e38a7b 100644 --- a/src/psplines.jl +++ b/src/psplines.jl @@ -120,8 +120,8 @@ end Base.eltype(::Type{<:_RInitializationBasis{T}}) where {T} = T nbasis(b::_RInitializationBasis) = nbasis(b.basis) domain(b::_RInitializationBasis) = domain(b.basis) -basismatrix(b::_RInitializationBasis, t::AbstractVector{<:Real}; derivative=0) = - basismatrix(b.basis, t; derivative) +basis_matrix(b::_RInitializationBasis, t::AbstractVector{<:Real}; derivative=0) = + basis_matrix(b.basis, t; derivative) penalty(b::_RInitializationBasis, order) = order == 2 ? copy(b.roughness) : penalty(b.basis, order) # Initialization only: reproduce stats::smooth.spline's default knot count and @@ -133,10 +133,10 @@ function _face_initial_curve(t, y, targets) nk = _spline_nknots(length(t)) indices = floor.(Int, range(1, length(t); length=nk)) cubic = BSplineBasis(4, x[indices]) - slopes = diff(basismatrix(cubic, cubic.breakpoints; derivative=2); dims=1) + slopes = diff(basis_matrix(cubic, cubic.breakpoints; derivative=2); dims=1) roughness = penalty(cubic, 2) - transpose(slopes)*(diff(cubic.breakpoints) .* slopes)/3000 basis = _RInitializationBasis(cubic, roughness) - X = basismatrix(basis, x) + X = basis_matrix(basis, x) P = penalty(basis, 2) XtX = transpose(X)*X inner = 3:(nbasis(basis)-3) diff --git a/src/recipes.jl b/src/recipes.jl index 58980d5..64daf16 100644 --- a/src/recipes.jl +++ b/src/recipes.jl @@ -60,14 +60,14 @@ end end """ - modesplot(fit, k; multipliers=(-2, -1, 1, 2), kwargs...) + modes_plot(fit, k; multipliers=(-2, -1, 1, 2), kwargs...) Plot the mean and mode-of-variation curves from [`modes_of_variation`](@ref). The default multipliers follow fdapace 0.6.0's envelope boundaries. Titles show the component, eigenvalue, and individual FVE. Multivariate fits have one panel per node. Load Plots.jl before calling; `kwargs` set plot attributes. """ -modesplot(fit, k; multipliers=(-2, -1, 1, 2), kwargs...) = +modes_plot(fit, k; multipliers=(-2, -1, 1, 2), kwargs...) = RecipesBase.plot(_plot_payload(Val(:modes), modes_of_variation(fit, k; multipliers)); kwargs...) @recipe function f(p::_FDAPlot{:modes}) @@ -92,14 +92,14 @@ modesplot(fit, k; multipliers=(-2, -1, 1, 2), kwargs...) = end """ - screeplot(fit; kwargs...) + scree_plot(fit; kwargs...) Plot individual FVE bars and cumulative FVE, following fdapace 0.6.0 `CreateScreePlot`, using [`scree_data`](@ref). Show selected K and the FVE threshold when applicable. All retained components are displayed. Load Plots.jl first. Plot attributes may be overridden with `kwargs`. """ -screeplot(fit; kwargs...) = RecipesBase.plot(_plot_payload(Val(:scree), scree_data(fit)); kwargs...) +scree_plot(fit; kwargs...) = RecipesBase.plot(_plot_payload(Val(:scree), scree_data(fit)); kwargs...) @recipe function f(p::_FDAPlot{:scree}) d = p.data @@ -136,16 +136,16 @@ screeplot(fit; kwargs...) = RecipesBase.plot(_plot_payload(Val(:scree), scree_da end """ - covplot(x; kwargs...) - covplot(f::FPCA; surface=:stored, kwargs...) + cov_plot(x; kwargs...) + cov_plot(f::FPCA; surface=:stored, kwargs...) Plot [`covariance_data`](@ref) as a heatmap, a package display default. Use `seriestype=:surface` for a surface plot. FPCA accepts `surface=:fitted` for the full retained positive-semidefinite projection. Load Plots.jl first; other keywords set plot attributes. """ -covplot(x; kwargs...) = RecipesBase.plot(_FDAPlot{:covariance}(covariance_data(x)); kwargs...) -covplot(f::FPCA; surface=:stored, kwargs...) = +cov_plot(x; kwargs...) = RecipesBase.plot(_FDAPlot{:covariance}(covariance_data(x)); kwargs...) +cov_plot(f::FPCA; surface=:stored, kwargs...) = RecipesBase.plot(_FDAPlot{:covariance}(covariance_data(f; surface)); kwargs...) @recipe function f(p::_FDAPlot{:covariance}) @@ -159,14 +159,14 @@ covplot(f::FPCA; surface=:stored, kwargs...) = end """ - designplot(fd; diagonal=false, kwargs...) + design_plot(fd; diagonal=false, kwargs...) Plot observed time pairs using [`design_data`](@ref). Diagonal pairs are excluded by default and color counts use 1, 2, 3, 4+, as in fdapace 0.6.0. Exact counts remain available in `design_data`. Multivariate data have one panel per node. Load Plots.jl first; `kwargs` override plot attributes. """ -designplot(fd; diagonal=false, kwargs...) = +design_plot(fd; diagonal=false, kwargs...) = RecipesBase.plot(_plot_payload(Val(:design), design_data(fd; diagonal)); kwargs...) @recipe function f(p::_FDAPlot{:design}) @@ -185,14 +185,14 @@ designplot(fd; diagonal=false, kwargs...) = end """ - fitplot(fit, i; grid=nothing, kwargs...) + fit_plot(fit, i; grid=nothing, kwargs...) Plot the original observations and selected-K reconstruction from [`fitted_curve_data`](@ref). `grid=nothing` uses the fit's stored grid(s), following fdapace 0.6.0 `CreatePathPlot`. MFPCA has one panel per node. Load Plots.jl first; `kwargs` override plot attributes. """ -function fitplot(fit, i; grid=nothing, kwargs...) +function fit_plot(fit, i; grid=nothing, kwargs...) data = grid === nothing ? fitted_curve_data(fit, i) : fitted_curve_data(fit, i; grid) RecipesBase.plot(_plot_payload(Val(:fit), data); kwargs...) end @@ -223,8 +223,8 @@ end end """ - bandplot(b; kwargs...) - bandplot(b::MeanBands; coverage=b.coverage, kwargs...) + band_plot(b; kwargs...) + band_plot(b::MeanBands; coverage=b.coverage, kwargs...) Plot precomputed [`fitted_bands`](@ref) or [`mean_bands`](@ref) results as ribbons around their fitted curves or means, using [`band_data`](@ref). @@ -233,9 +233,9 @@ computed types. Multivariate results use one panel per node. Load Plots.jl first; `kwargs` override plot attributes. Default ribbon opacity is 0.2, a package display choice. No bootstrap is run by plotting. """ -bandplot(b; kwargs...) = +band_plot(b; kwargs...) = RecipesBase.plot(_plot_payload(Val(:band), band_data(b)); kwargs...) -bandplot(b::MeanBands; coverage=b.coverage, kwargs...) = +band_plot(b::MeanBands; coverage=b.coverage, kwargs...) = RecipesBase.plot(_FDAPlot{:band}(band_data(b; coverage)); kwargs...) @recipe function f(p::_FDAPlot{:band}) diff --git a/src/smoothed_fpca.jl b/src/smoothed_fpca.jl index 186505d..032a53a 100644 --- a/src/smoothed_fpca.jl +++ b/src/smoothed_fpca.jl @@ -35,8 +35,8 @@ end Base.eltype(::Type{<:_BasisEigenfunctions{T}}) where {T} = T nbasis(b::_BasisEigenfunctions) = size(b.coefs, 2) domain(b::_BasisEigenfunctions) = domain(b.basis) -basismatrix(b::_BasisEigenfunctions, t::AbstractVector{<:Real}; derivative=0) = - basismatrix(b.basis, t; derivative) * b.coefs +basis_matrix(b::_BasisEigenfunctions, t::AbstractVector{<:Real}; derivative=0) = + basis_matrix(b.basis, t; derivative) * b.coefs penalty(b::_BasisEigenfunctions, order) = transpose(b.coefs) * penalty(b.basis, order) * b.coefs @@ -84,7 +84,7 @@ function _basis_cholesky(matrix, label) end # Reuse the common eigenvalue cutoff, spectrum accounting, and FVE machinery. -function _basis_pca(residuals::Matrix{T}, G, P, lambda, maxK, eigen_tol, +function _basis_pca(residuals::Matrix{T}, G, P, lambda, max_k, eigen_tol, sign_tol) where {T} n, k = size(residuals) U = _basis_cholesky(G + lambda * P, "FPCA penalized Gram matrix").U @@ -92,7 +92,7 @@ function _basis_pca(residuals::Matrix{T}, G, P, lambda, maxK, eigen_tol, covariance = transpose(residuals) * residuals / n transformed = J * covariance * transpose(J) transformed = Matrix(Symmetric((transformed + transpose(transformed)) / 2)) - system = _fpca_eigen(transformed, ones(T, k), zeros(T, k); maxK, + system = _fpca_eigen(transformed, ones(T, k), zeros(T, k); max_k, eigen_tol, sign_tol) raw = U \ system.eigenfunctions norms = sqrt.(vec(sum(raw .* (G * raw); dims=1))) @@ -108,11 +108,11 @@ function _basis_cv_error(err, lambda, i) end function _basis_pca_lambda(lambda::Real, residuals::Matrix{T}, G, P, - rule, maxK, eigen_tol, sign_tol) where {T} + rule, max_k, eigen_tol, sign_tol) where {T} _basis_lambda(lambda, T), nothing end function _basis_pca_lambda(cv::LambdaCV, residuals::Matrix{T}, G, P, - rule, maxK, eigen_tol, sign_tol) where {T} + rule, max_k, eigen_tol, sign_tol) where {T} rule isa Integer && !(rule isa Bool) && rule >= 1 || throw(ArgumentError( "FPCA LambdaCV requires a positive integer ncomponents so all folds and candidates use the same subspace count")) n = size(residuals, 1) @@ -125,8 +125,8 @@ function _basis_pca_lambda(cv::LambdaCV, residuals::Matrix{T}, G, P, for (c, lam) in enumerate(candidates), i in 1:n try train = residuals[[j for j in 1:n if j != i], :] - system = _basis_pca(train, G, P, lam, maxK, eigen_tol, sign_tol) - chosen = _fpca_selection(rule, system.fve, maxK).chosen + system = _basis_pca(train, G, P, lam, max_k, eigen_tol, sign_tol) + chosen = _fpca_selection(rule, system.fve, max_k).chosen for m in 1:chosen H = system.coefs[:, 1:m] rhs = transpose(H) * G * residuals[i, :] @@ -149,7 +149,7 @@ _basis_pca_lambda(lambda, args...) = throw(ArgumentError("lambda must be a nonne """ fit(FPCA, fd::BasisFunctionalData; lambda=0, penalty=nothing, - center=true, ncomponents=FVE(0.99), maxK=nothing, grid=nothing, + center=true, ncomponents=FVE(0.99), max_k=nothing, maxK=nothing, grid=nothing, eigen_tol=1e-10, sign_tol=1e-12) Roughness-penalized basis FPCA, following fda 6.3.0 `pca.fd`/`fdPar` and @@ -167,7 +167,7 @@ Each harmonic is divided by its exact L² norm q; scores are R's scores divided by q. Thus score variances (denominator n) are ν/q², and positive-lambda harmonics are individually L²-normalized but generally not L²-orthogonal. The package's `FVE(0.99)` default departs from R's `nharm=2`; a fixed integer -reproduces `nharm`. The default `maxK=min(nbasis,n-1)` (n when uncentered) +reproduces `nharm`. The default `max_k=min(nbasis,n-1)` (n when uncentered) retains the available rank, unlike the grid routes' caps. The same FPCA result supports exact basis prediction and reconstruction, @@ -178,15 +178,16 @@ No noise model is assumed; noise and score covariances are `nothing`. All settings and input copies are retained. See the FPCA page for conversions. """ function fit(::Type{FPCA}, fd::BasisFunctionalData{T}; lambda=0, penalty=nothing, - center=true, ncomponents=FVE(0.99), maxK=nothing, grid=nothing, + center=true, ncomponents=FVE(0.99), max_k=nothing, maxK=nothing, grid=nothing, eigen_tol=_FPCA_EIGEN_TOL, sign_tol=_FPCA_SIGN_TOL) where {T} + max_k = _max_k_keyword(max_k, maxK) center isa Bool || throw(ArgumentError("center must be Boolean")) n, k = size(fd.coefs) n >= 2 || throw(ArgumentError("basis FPCA requires at least two curves")) limit = min(k, center ? n-1 : n) - cap = maxK === nothing ? limit : maxK + cap = max_k === nothing ? limit : max_k cap isa Integer && !(cap isa Bool) && 1 <= cap <= limit || - throw(ArgumentError("basis FPCA maxK must be between 1 and min(nbasis, centered rank bound)")) + throw(ArgumentError("basis FPCA max_k must be between 1 and min(nbasis, centered rank bound)")) _fpca_tolerance(eigen_tol, "eigen_tol") _fpca_tolerance(sign_tol, "sign_tol") b = deepcopy(fd.basis) @@ -205,7 +206,7 @@ function fit(::Type{FPCA}, fd::BasisFunctionalData{T}; lambda=0, penalty=nothing points = grid === nothing ? collect(range(T.(domain(b))...; length=201)) : Vector{T}(_basispoints(b, grid)) weights = quadrature_weights(points) - X = Matrix{T}(basismatrix(b, points)) + X = Matrix{T}(basis_matrix(b, points)) values = X * system.coefs oriented = _fpca_orient!(copy(values), X * mu, sign_tol) H = copy(system.coefs) @@ -221,7 +222,7 @@ function fit(::Type{FPCA}, fd::BasisFunctionalData{T}; lambda=0, penalty=nothing all(isfinite, projected) && all(isfinite, covariance) && all(isfinite, fitted) || throw(ArgumentError("basis FPCA results must be finite; rescale data")) settings = (; lambda=lam, penalty=specification, penalty_matrix=P, center, - lambda_selection=tuning, ncomponents=deepcopy(ncomponents), maxK=Int(cap), + lambda_selection=tuning, ncomponents=deepcopy(ncomponents), max_k=Int(cap), maxK=Int(cap), eigen_tol, sign_tol, covariance_denominator=n, normalization=:exact_L2, original_harmonic_norms=system.normalization, eigenvalue_scale=:penalized_metric, grid_source=grid === nothing ? :display_201 : :user_display, diff --git a/src/smoothing.jl b/src/smoothing.jl index a918ee0..a55b9a3 100644 --- a/src/smoothing.jl +++ b/src/smoothing.jl @@ -112,7 +112,7 @@ function smooth_curves(fd::Union{GridFunctionalData,IrregularFunctionalData}, indices, L = get!(smoothers, t) do indices = findall(z -> first(t) <= z <= last(t), x) isempty(indices) && return indices, zeros(0, length(t)) - indices, KernelSmoothing.smoothermatrix(t, x[indices]; bandwidth=method.bandwidth, + indices, KernelSmoothing.smoother_matrix(t, x[indices]; bandwidth=method.bandwidth, degree=method.degree, derivative=method.derivative, kernel=method.kernel) end isempty(indices) && continue @@ -138,7 +138,7 @@ function _gcv_selection(fd, basis, order, candidates) for i in 1:nobs(fd) t, y = fd[i] isempty(t) && throw(ArgumentError("curve $i has no observed values for GCV")) - X = Matrix{T}(basismatrix(basis, t)) + X = Matrix{T}(basis_matrix(basis, t)) yy = Vector{T}(y) n = length(y) XtX = transpose(X) * X diff --git a/src/smoothing2d.jl b/src/smoothing2d.jl index 0f7c82b..0eeefa8 100644 --- a/src/smoothing2d.jl +++ b/src/smoothing2d.jl @@ -1,7 +1,7 @@ # Product-kernel local linear covariance surface on the grid g1 × g2, from # KernelSmoothing's bivariate smoother. Unidentified fits throw with their point. function _local_linear_surface(pool, g1, g2, h, kernel) - surface = KernelSmoothing.localpolynomial((pool.x, pool.y), pool.z, (g1, g2); + surface = KernelSmoothing.local_polynomial((pool.x, pool.y), pool.z, (g1, g2); bandwidth=h, degree=1, kernel, weights=pool.weights) for j in axes(surface, 2), i in axes(surface, 1) isfinite(surface[i, j]) || throw(ArgumentError( diff --git a/test/bands.jl b/test/bands.jl index 44e6816..4689326 100644 --- a/test/bands.jl +++ b/test/bands.jl @@ -15,7 +15,7 @@ function bands_fixture(; K=3, scalar=Float64) values = [input[input[:, 1] .== i, 3] for i in 1:n] fd = IrregularFunctionalData(times, values) model = StatsAPI.fit(FPCA, fd; grid, mean=(; grid, values=mean_input[:, 2]), - covariance=(; grid, surface), noise_variance=sigma2, ncomponents=K, maxK=5, + covariance=(; grid, surface), noise_variance=sigma2, ncomponents=K, max_k=5, eigen_tol=scalar === Float32 ? 1e-5 : 1e-10) (; fd, model, grid) end @@ -52,7 +52,7 @@ end mu = (; grid, values=1 .+ grid) covariance = (; grid, surface=fill(2.0, length(grid), length(grid))) f = StatsAPI.fit(FPCA, fd; grid, mean=mu, covariance, noise_variance=0.5, - ncomponents=1, maxK=1) + ncomponents=1, max_k=1) query = [0.0, 0.1, 0.8, 1.0] b = @inferred fitted_bands(f; grid=query, curves=[2, 1]) # A random intercept with prior variance 2 and two noisy measurements: @@ -74,7 +74,7 @@ end @test_throws ArgumentError fitted_bands(indefinite) mfd = MultivariateFunctionalData([fd, fd]) mf = StatsAPI.fit(FPCA, mfd; grid, mean=mu, covariance, noise_variance=0.5, - ncomponents=1, maxK=1) + ncomponents=1, max_k=1) mb = @inferred fitted_bands(mf; curves=2) @test length(mb.nodes) == 2 @test mb.nodes[1].variance ≈ fill(2/9, 1, 4) @@ -83,10 +83,10 @@ end @test off.nodes[2].grid == grid @test_throws DimensionMismatch fitted_bands(mf; grid=[grid]) dense = StatsAPI.fit(FPCA, GridFunctionalData([0.0 1 2 3; 1.0 2 3 4; 2.0 3 4 5], grid); - ncomponents=1, maxK=1) + ncomponents=1, max_k=1) @test_throws ArgumentError fitted_bands(dense) @test_throws ArgumentError fitted_bands(StatsAPI.fit(FPCA, - MultivariateFunctionalData([dense.training_data, dense.training_data]); ncomponents=1, maxK=1)) + MultivariateFunctionalData([dense.training_data, dense.training_data]); ncomponents=1, max_k=1)) for bad in (0, 1, -0.1, NaN, Inf, true, :auto) @test_throws ArgumentError fitted_bands(f; alpha=bad) end diff --git a/test/bands_coverage.jl b/test/bands_coverage.jl index c39bcbb..8498e09 100644 --- a/test/bands_coverage.jl +++ b/test/bands_coverage.jl @@ -25,7 +25,7 @@ function bands_coverage(; replicates=400, n=300, B=399, seed=20261025) fd = IrregularFunctionalData(times, values; domain=(0.0, 1.0)) model = StatsAPI.fit(FPCA, fd; method=:pace, grid, mean=(; grid, values=truth_mean), covariance=(; grid, surface), - noise_variance=0.2, ncomponents=2, maxK=2) + noise_variance=0.2, ncomponents=2, max_k=2) truth = transpose(truth_mean) .+ latent * phi' point = fitted_bands(model) band = fitted_bands(model; coverage=:band) diff --git a/test/bases.jl b/test/bases.jl index 0074fcf..88d01cb 100644 --- a/test/bases.jl +++ b/test/bases.jl @@ -9,22 +9,22 @@ @test domain(b) == (0.0, 1.0) @test nbasis(b) == 6 t = collect(range(0.0, 1.0; length=101)) - X = @inferred basismatrix(b, t) + X = @inferred basis_matrix(b, t) @test vec(sum(X; dims=2)) ≈ ones(length(t)) @test X[1, 1] == X[end, end] == 1 - @test vec(sum(basismatrix(b, t; derivative=1); dims=2)) ≈ + @test vec(sum(basis_matrix(b, t; derivative=1); dims=2)) ≈ zeros(length(t)) atol=1e-13 points = [0.13, 0.33, 0.72] h = 1e-6 for m in 1:2 - difference = (basismatrix(b, points .+ h; derivative=m - 1) - - basismatrix(b, points .- h; derivative=m - 1)) / (2h) - @test basismatrix(b, points; derivative=m) ≈ difference rtol=1e-8 + difference = (basis_matrix(b, points .+ h; derivative=m - 1) - + basis_matrix(b, points .- h; derivative=m - 1)) / (2h) + @test basis_matrix(b, points; derivative=m) ≈ difference rtol=1e-8 end - @test basismatrix(b, reverse(t)) ≈ reverse(X; dims=1) - @test basismatrix(b, t; derivative=4) == zeros(length(t), 6) + @test basis_matrix(b, reverse(t)) ≈ reverse(X; dims=1) + @test basis_matrix(b, t; derivative=4) == zeros(length(t), 6) @test penalty(b, 4) == zeros(6, 6) - @test size(basismatrix(b, Float64[])) == (0, 6) + @test size(basis_matrix(b, Float64[])) == (0, 6) @test nbasis(@inferred BSplineBasis(4, (0, 1), 7)) == 7 @test eltype(@inferred BSplineBasis(2, Float32[0, 0.5, 1])) == Float32 @test eltype(@inferred BSplineBasis(2, [0, 1])) == Float64 @@ -37,17 +37,17 @@ @test domain(f) == (0.0, 2.0) @test gram(f) ≈ Matrix{Float64}(I, 5, 5) @test diag(penalty(f, 1)) ≈ [0, pi^2, pi^2, (2pi)^2, (2pi)^2] - @test basismatrix(f, [0.0]) ≈ reshape([inv(sqrt(2)), 0, 1, 0, 1], 1, :) + @test basis_matrix(f, [0.0]) ≈ reshape([inv(sqrt(2)), 0, 1, 0, 1], 1, :) t = [0.17, 0.63, 1.21] h = 1e-6 - @test basismatrix(f, t; derivative=1) ≈ - (basismatrix(f, t .+ h) - basismatrix(f, t .- h)) / (2h) rtol=1e-8 + @test basis_matrix(f, t; derivative=1) ≈ + (basis_matrix(f, t .+ h) - basis_matrix(f, t .- h)) / (2h) rtol=1e-8 b = @inferred MonomialBasis((0, 1), 3) @test nbasis(b) == 4 @test gram(b) ≈ [1 / (i + j + 1) for i in 0:3, j in 0:3] @test penalty(b, 2) ≈ [0 0 0 0; 0 0 0 0; 0 0 4 6; 0 0 6 12] @test penalty(b, 4) == zeros(4, 4) - @test basismatrix(b, [0.0, 1.0]; derivative=2) == [0 0 2 0; 0 0 2 6] + @test basis_matrix(b, [0.0, 1.0]; derivative=2) == [0 0 2 0; 0 0 2 6] @test eltype(@inferred FourierBasis((0.0f0, 1.0f0), 3)) == Float32 @test eltype(@inferred MonomialBasis((0.0f0, 1.0f0), 2)) == Float32 end @@ -62,22 +62,22 @@ @test (@inferred domain(b)) == (-1.0, 2.0) @test (@inferred gram(b)) isa Matrix{Float64} for m in 0:2 - X = @inferred basismatrix(b, t; derivative=m) + X = @inferred basis_matrix(b, t; derivative=m) P = @inferred penalty(b, m) numerical = transpose(X) * (w .* X) @test norm(numerical - P) / norm(P) < 1e-6 end for badpoints in ([-1.01], [2.01], [NaN], [Inf]) - @test_throws ArgumentError basismatrix(b, badpoints) + @test_throws ArgumentError basis_matrix(b, badpoints) end - @test_throws ArgumentError basismatrix(b, [0.0]; derivative=-1) - @test_throws ArgumentError basismatrix(b, [0.0]; derivative=0.5) + @test_throws ArgumentError basis_matrix(b, [0.0]; derivative=-1) + @test_throws ArgumentError basis_matrix(b, [0.0]; derivative=0.5) @test_throws ArgumentError penalty(b, -1) @test_throws ArgumentError penalty(b, 0.5) end for b in (BSplineBasis(3, Float32[0, 0.5, 1]), FourierBasis((0.0f0, 1.0f0), 3), MonomialBasis((0.0f0, 1.0f0), 2)) - @test (@inferred basismatrix(b, Float32[0, 0.25, 1])) isa Matrix{Float32} + @test (@inferred basis_matrix(b, Float32[0, 0.25, 1])) isa Matrix{Float32} @test (@inferred gram(b)) isa Matrix{Float32} @test (@inferred penalty(b, 1)) isa Matrix{Float32} @test b == deepcopy(b) diff --git a/test/covariance.jl b/test/covariance.jl index 1cd2f33..158bd4f 100644 --- a/test/covariance.jl +++ b/test/covariance.jl @@ -293,7 +293,7 @@ end @test grouped.scatter[i] ≈ sum(w .* (z .- grouped.z[i]).^2) atol=1e-14 end fit = StatsAPI.fit(SmoothedCovariance, fd, m; diagonal_bandwidth=.65, grid, bandwidth=1.2, weighting, binning=5) - expected = KernelSmoothing.localpolynomial((x, y), raw.z, (grid, grid); + expected = KernelSmoothing.local_polynomial((x, y), raw.z, (grid, grid); bandwidth=1.2, kernel=fit.kernel, weights=raw.weights) @test fit.surface ≈ (expected + expected')/2 selector = CurveCV([.8, 1.2]; folds=2, rng=MersenneTwister(3)) @@ -304,7 +304,7 @@ end score = 0.0 for i in eachindex(raw.z) train = assignments[raw.subjects] .!= assignments[raw.subjects[i]] - prediction = only(KernelSmoothing.localpolynomial((x[train], y[train]), + prediction = only(KernelSmoothing.local_polynomial((x[train], y[train]), raw.z[train], [(x[i], y[i])]; bandwidth=h, kernel=fit.kernel, weights=raw.weights[train])) score += raw.weights[i]*abs2(raw.z[i]-prediction) end @@ -426,7 +426,7 @@ end x = binning === nothing ? original : round.(original .* 8, RoundNearestTiesUp) ./ 8 w = vcat([fill(weighting == :obs ? 1. : 1/length(t), length(t)) for t in times]...) for rule in (KernelSmoothing.RuleOfThumb(), KernelSmoothing.IMSEOptimal()) - direct = KernelSmoothing.selectbandwidth(rule, + direct = KernelSmoothing.select_bandwidth(rule, KernelSmoothing.SmoothingData(x, squares; weights=w), kernel; degree=1) # KernelSmoothing selectors' `details` types depend on the data (and on # the KernelSmoothing version), so selector fits are not inferred. @@ -443,7 +443,7 @@ end cv = CurveCV(; pilot, multipliers, folds=2, rng=MersenneTwister(52)) c = StatsAPI.fit(SmoothedCovariance, fd, m; grid, kernel, bandwidth=.5, diagonal_bandwidth=cv, weighting, binning) - h = KernelSmoothing.selectbandwidth(pilot, + h = KernelSmoothing.select_bandwidth(pilot, KernelSmoothing.SmoothingData(x, squares; weights=w), kernel; degree=1).h candidates = multipliers .* h explicit = StatsAPI.fit(SmoothedCovariance, fd, m; diff --git a/test/cross_covariance.jl b/test/cross_covariance.jl index 467fb39..d124832 100644 --- a/test/cross_covariance.jl +++ b/test/cross_covariance.jl @@ -207,7 +207,7 @@ end blocks[2,1].surface[1,1] += 1 @test blocks[1,2].surface != transpose(blocks[2,1].surface) for basis in (:marginal,:pooled) - fpca = StatsAPI.fit(FPCA,mfd; basis, ncomponents=1, maxK=2) + fpca = StatsAPI.fit(FPCA,mfd; basis, ncomponents=1, max_k=2) b = cross_covariances(mfd; fpca, bandwidth=2., grid=[grid, [0.,.5,1.]]) @test b[1,2].mean1.values == fpca.means[1].values @test b[1,2].mean2.values == fpca.means[2].values @@ -241,9 +241,9 @@ end end wrongids = MultivariateFunctionalData([GridFunctionalData(values,grid; ids=6:-1:1), GridFunctionalData(values,grid; ids=6:-1:1)]) - fpca = StatsAPI.fit(FPCA,wrongids; ncomponents=1,maxK=2) + fpca = StatsAPI.fit(FPCA,wrongids; ncomponents=1,max_k=2) @test_throws ArgumentError cross_covariances(mfd; fpca,bandwidth=2.) - fewer = StatsAPI.fit(FPCA,mfd[1:5]; ncomponents=1,maxK=2) + fewer = StatsAPI.fit(FPCA,mfd[1:5]; ncomponents=1,max_k=2) @test_throws DimensionMismatch cross_covariances(mfd; fpca=fewer,bandwidth=2.) @test size(cross_covariances(MultivariateFunctionalData([x]); bandwidth=2.)) == (1,1) end diff --git a/test/cross_structure.jl b/test/cross_structure.jl index 112bf7a..edcede7 100644 --- a/test/cross_structure.jl +++ b/test/cross_structure.jl @@ -9,7 +9,7 @@ _m6_relative(a, b) = norm(a-b)/norm(b) l, r = _m6_read("dense_left"), _m6_read("dense_right") x = GridFunctionalData(_m6_read("dense_input1"), l[:,1]; ids=1:18) y = GridFunctionalData(_m6_read("dense_input2"), r[:,1]; ids=1:18) - f = @inferred StatsAPI.fit(FSVD, x, y; ncomponents=2, maxK=2) + f = @inferred StatsAPI.fit(FSVD, x, y; ncomponents=2, max_k=2) expected = _m6_read("dense_values") @test _m6_relative(f.singular_values, expected[:,1]) < 1e-8 @test _m6_relative(f.left_functions, l[:,2:end]) < 1e-8 @@ -23,8 +23,8 @@ _m6_relative(a, b) = norm(a-b)/norm(b) @test f.left_functions' * Diagonal(x.weights) * f.left_functions ≈ I @test f.right_functions' * Diagonal(y.weights) * f.right_functions ≈ I @test f.settings.effective_fpca_options1 == f.settings.effective_fpca_options2 == (;method=:auto) - @test StatsAPI.fit(FSVD,x,y).selection.maxK == 16 - @test StatsAPI.fit(FSVD,x,y;ncomponents=FVE(.5),maxK=2).selection.chosen == 1 + @test StatsAPI.fit(FSVD,x,y).selection.max_k == 16 + @test StatsAPI.fit(FSVD,x,y;ncomponents=FVE(.5),max_k=2).selection.chosen == 1 @test occursin("denominator=16", sprint(show, f)) @test f.cross_covariance.ids == x.ids @test f.cross_covariance.ids !== x.ids @@ -47,7 +47,7 @@ _m6_relative(a, b) = norm(a-b)/norm(b) @test eltype((@inferred StatsAPI.fit(FSVD,xx,small;ncomponents=1)).singular_values) == Float64 @test_throws ArgumentError StatsAPI.fit(FSVD,xx,yy; fpca_options1=(;presmooth=LocalPolynomial(bandwidth=1.))) - for kw in ((;maxK=0), (;ncomponents=0), (;maxK=1,ncomponents=2), (;maxK=true), (;ncomponents=true), + for kw in ((;max_k=0), (;ncomponents=0), (;max_k=1,ncomponents=2), (;max_k=true), (;ncomponents=true), (;bandwidth=.5), (;fpca_options1=[]), (;kernel=:gauss)) @test_throws ArgumentError StatsAPI.fit(FSVD,xx,yy;kw...) end @@ -71,10 +71,10 @@ end phi = hcat(ones(51),g .- .5) cov1, cov2 = phi*Diagonal([.8,.3])*phi', phi*Diagonal([.7,.2])*phi' o1 = (;grid=g, mean=(;grid=g,values=1 .+ g), covariance=(;grid=g,surface=cov1), - noise_variance=.2, ncomponents=2,maxK=2) + noise_variance=.2, ncomponents=2,max_k=2) o2 = (;grid=g, mean=(;grid=g,values=-.2 .+ .3g), covariance=(;grid=g,surface=cov2), - noise_variance=.2,ncomponents=2,maxK=2) - f = @inferred StatsAPI.fit(FSVD,nodes...;ncomponents=2,maxK=2,bandwidth=(.55,.8), + noise_variance=.2,ncomponents=2,max_k=2) + f = @inferred StatsAPI.fit(FSVD,nodes...;ncomponents=2,max_k=2,bandwidth=(.55,.8), kernel=KernelSmoothing.EpanechnikovKernel(),fpca_options1=o1,fpca_options2=o2) @test _m6_relative(f.cross_covariance.surface,_m6_read("sparse_cross")) < 1e-8 @test _m6_relative(f.singular_values,_m6_read("sparse_values")[:,1]) < 1e-8 @@ -105,7 +105,7 @@ end [v + .4randn(rng, length(v)) for v in observations]; ids=1:18) # The complete grid has more points than the sample covariance's rank. @test length(g) > nobs(x)-1 - f = @inferred StatsAPI.fit(FSVD,x,y;bandwidth=.6,ncomponents=2,maxK=2) + f = @inferred StatsAPI.fit(FSVD,x,y;bandwidth=.6,ncomponents=2,max_k=2) @test f.fpca1.method == f.fpca2.method == :pace @test f.fpca1.noise_variance > 0 && f.fpca2.noise_variance > 0 @test length(f.fpca1.grid) == length(f.fpca2.grid) == 51 @@ -115,24 +115,24 @@ end @test f.settings.effective_fpca_options1 == f.settings.effective_fpca_options2 == (;method=:pace) @test occursin("effective_fpca_options1 = (method = :pace,)",sprint(show,f)) # Changing only the representation must preserve the conditional scores. - irregular = StatsAPI.fit(FSVD,to_irregular(x),y;bandwidth=.6,ncomponents=2,maxK=2) + irregular = StatsAPI.fit(FSVD,to_irregular(x),y;bandwidth=.6,ncomponents=2,max_k=2) @test f.singular_values ≈ irregular.singular_values @test f.left_functions ≈ irregular.left_functions @test f.right_functions ≈ irregular.right_functions @test f.left_scores ≈ irregular.left_scores @test f.right_scores ≈ irregular.right_scores - reversed = @inferred StatsAPI.fit(FSVD,y,x;bandwidth=.6,ncomponents=2,maxK=2) + reversed = @inferred StatsAPI.fit(FSVD,y,x;bandwidth=.6,ncomponents=2,max_k=2) @test reversed.fpca1.method == reversed.fpca2.method == :pace @test reversed.singular_values ≈ f.singular_values @test all(isfinite,reversed.left_scores) && all(isfinite,reversed.right_scores) masked = StatsAPI.fit(FSVD,x,to_grid(y,g,NoInterpolation()); - bandwidth=.6,ncomponents=2,maxK=2) + bandwidth=.6,ncomponents=2,max_k=2) @test masked.fpca1.method == masked.fpca2.method == :pace @test masked.left_scores ≈ f.left_scores @test masked.right_scores ≈ f.right_scores # Supplying noise alone must work without also supplying method=:pace. options = (;method=:auto,noise_variance=.25) - explicit = StatsAPI.fit(FSVD,x,y;bandwidth=.6,ncomponents=2,maxK=2, + explicit = StatsAPI.fit(FSVD,x,y;bandwidth=.6,ncomponents=2,max_k=2, fpca_options1=options) @test explicit.fpca1.noise_variance == .25 @test explicit.settings.fpca_options1 == options diff --git a/test/data.jl b/test/data.jl index ef0e39c..ca5447e 100644 --- a/test/data.jl +++ b/test/data.jl @@ -1,7 +1,7 @@ struct TestLinearBasis <: AbstractBasis end FunctionalDataAnalysis.domain(::TestLinearBasis) = (0.0, 1.0) FunctionalDataAnalysis.nbasis(::TestLinearBasis) = 2 -function FunctionalDataAnalysis.basismatrix(::TestLinearBasis, t; derivative=0) +function FunctionalDataAnalysis.basis_matrix(::TestLinearBasis, t; derivative=0) derivative == 0 && return hcat(ones(length(t)), t) derivative == 1 && return hcat(zeros(length(t)), ones(length(t))) zeros(length(t), 2) @@ -99,9 +99,9 @@ FunctionalDataAnalysis.penalty(b::TestLinearBasis, order) = @test domain(fd) == (0.0, 1.0) @test nobs_per_curve(fd) == [3, 3] t = [0.0, 0.25, 1.0] - @test (@inferred evaluate(fd, t)) == fd.coefs * transpose(basismatrix(b, t)) + @test (@inferred evaluate(fd, t)) == fd.coefs * transpose(basis_matrix(b, t)) @test (@inferred fd(t)) == evaluate(fd, t) - @test (@inferred evaluate(fd, t; derivative=1)) == fd.coefs * transpose(basismatrix(b, t; derivative=1)) + @test (@inferred evaluate(fd, t; derivative=1)) == fd.coefs * transpose(basis_matrix(b, t; derivative=1)) @test (@inferred BasisFunctionalData(Float32[1 2], MonomialBasis((0.0f0, 1.0f0), 1))) isa BasisFunctionalData{Float32} spline = BSplineBasis(3, [0.0, 0.5, 1.0]) sf = BasisFunctionalData(ones(1, nbasis(spline)), spline) @@ -170,7 +170,7 @@ FunctionalDataAnalysis.penalty(b::TestLinearBasis, order) = t = collect(range(0.0, 1.0; length=31)) for basis in (MonomialBasis((0, 1), 2), FourierBasis((0, 1), 5), BSplineBasis(4, (0, 1), 6)) c = reshape(collect(1.0:2nbasis(basis)), 2, :) - values = c * transpose(basismatrix(basis, t)) + values = c * transpose(basis_matrix(basis, t)) fd = GridFunctionalData(values, t) fit = @inferred to_basis(fd, basis) @test fit.coefs ≈ c rtol=1e-10 @@ -181,7 +181,7 @@ FunctionalDataAnalysis.penalty(b::TestLinearBasis, order) = b = MonomialBasis((0, 1), 2) values = reshape(sin.(3t), 1, :) raw = GridFunctionalData(values, t) - X = basismatrix(b, t) + X = basis_matrix(b, t) lambda = 0.3 expected = (transpose(X) * X + lambda * penalty(b, 2)) \ (transpose(X) * vec(values)) @test to_basis(raw, b; lambda).coefs ≈ transpose(expected) diff --git a/test/derivatives.jl b/test/derivatives.jl index f9bfe89..3dd419f 100644 --- a/test/derivatives.jl +++ b/test/derivatives.jl @@ -23,7 +23,7 @@ import RecipesBase @test gram(d.basis)[4,4] ≈ 9/5 @test penalty(d.basis,1)[4,4] ≈ 12 @test penalty(d.basis,3) == zeros(4,4) - @test basismatrix(d.basis,t;derivative=1) ≈ basismatrix(b,t;derivative=2) + @test basis_matrix(d.basis,t;derivative=1) ≈ basis_matrix(b,t;derivative=2) @test domain(d.basis) == (0.,1.) @test d.basis == deepcopy(d.basis) @test isequal(d.basis,deepcopy(d.basis)) @@ -56,7 +56,7 @@ import RecipesBase for nu in (-1,1.5,NaN,big(typemax(Int))+1) @test_throws ArgumentError derivative(fd,nu) end - @test_throws ArgumentError basismatrix(DerivativeBasis(b,typemax(Int)),t;derivative=1) + @test_throws ArgumentError basis_matrix(DerivativeBasis(b,typemax(Int)),t;derivative=1) @test_throws ArgumentError evaluate(d,[-.1]) end @@ -109,7 +109,7 @@ end covariance = (; grid, surface=ones(length(grid), length(grid))) mfd = MultivariateFunctionalData([raw, raw]) @test (@inferred StatsAPI.fit(FPCA, raw; method=:pace, grid, mean=m, - covariance, noise_variance=.1, ncomponents=1, maxK=1)).method === :pace + covariance, noise_variance=.1, ncomponents=1, max_k=1)).method === :pace @test (@inferred StatsAPI.fit(SmoothedCovariance, raw, m; grid, bandwidth=2.)) isa SmoothedCovariance @test (@inferred StatsAPI.fit(SmoothedCrossCovariance, raw, raw, m, m; @@ -119,9 +119,9 @@ end grid, bandwidth=2., derivative=nu, degree=3) err = ArgumentError("mean has derivative=$nu; a centering mean must have derivative=0") @test_throws err StatsAPI.fit(FPCA, raw; method=:pace, grid, mean=dm, - covariance, noise_variance=.1, ncomponents=1, maxK=1) + covariance, noise_variance=.1, ncomponents=1, max_k=1) @test_throws err StatsAPI.fit(FPCA, raw; method=:pace, grid, mean=dm, - covariance_bandwidth=2., noise_variance=.1, ncomponents=1, maxK=1) + covariance_bandwidth=2., noise_variance=.1, ncomponents=1, max_k=1) @test_throws err StatsAPI.fit(SmoothedCovariance, raw, dm; grid, bandwidth=2.) for means in ([dm, m], [m, dm]) @@ -145,13 +145,13 @@ function derivative_fixture(design) values = [Float64.(data[data[:,1].==i,3]) for i in 1:n] if design=="dense" fd = GridFunctionalData(permutedims(hcat(values...)),first(times);ids=collect(1:n)) - f = StatsAPI.fit(FPCA,fd;ncomponents=3,maxK=5) + f = StatsAPI.fit(FPCA,fd;ncomponents=3,max_k=5) else m = read("m3_sparse_mean.csv"); grid=Float64.(m[:,1]) fd = IrregularFunctionalData(times,values;ids=collect(1:n)) f = StatsAPI.fit(FPCA,fd;grid,mean=(;grid,values=m[:,2]+.2grid.^2), covariance=(;grid,surface=read("m3_sparse_surface.csv")), - noise_variance=.12,ncomponents=3,maxK=5) + noise_variance=.12,ncomponents=3,max_k=5) end fd,f end @@ -244,13 +244,13 @@ end t = [0.,.03,.1,.2,.4,.65,.8,.9,1.] v = collect(-3.:3.) .* transpose(t.^2 .+ t) .+ transpose(1 .+ t) fd = GridFunctionalData(v,t;ids=collect(1:7)) - f = StatsAPI.fit(FPCA,fd;ncomponents=1,maxK=1) + f = StatsAPI.fit(FPCA,fd;ncomponents=1,max_k=1) d = @inferred (f -> derivative(f;method=:DPC,bandwidth=2.,ncomponents=1))(f) # The quadratic local derivative is exact on this nonuniform grid. @test reconstruct(d).values ≈ collect(-3.:3.) .* transpose(2t .+ 1) .+ 1 atol=1e-10 mfd = MultivariateFunctionalData([fd,fd]) for mode in (:marginal,:pooled) - mf = StatsAPI.fit(FPCA,mfd;basis=mode,ncomponents=1,maxK=1) + mf = StatsAPI.fit(FPCA,mfd;basis=mode,ncomponents=1,max_k=1) md = derivative(mf;bandwidth=[.2,.3]) @test md isa MultivariateFPCADerivative @test nobs(md)==7 && ncomponents(md)==[1,1] @@ -278,7 +278,7 @@ end fd2 = GridFunctionalData(collect(-3.:3.) .* transpose(t.^3 .- .5t) .+ transpose(1 .+ t.^2), t; ids=collect(1:7)) mixed = MultivariateFunctionalData([fd, fd2]) - pooled = StatsAPI.fit(FPCA, mixed; basis=:pooled, ncomponents=1, maxK=2) + pooled = StatsAPI.fit(FPCA, mixed; basis=:pooled, ncomponents=1, max_k=2) dd = derivative(pooled; method=:DPC, bandwidth=2., degree=3, ncomponents=1) truth2 = collect(-3.:3.) .* transpose(3t.^2 .- .5) .+ transpose(2t) @test component(reconstruct(dd), 2).values ≈ truth2 atol=1e-10 @@ -291,7 +291,7 @@ end sparsefd, sparsefit = derivative_fixture("sparse") mixed_design = MultivariateFunctionalData([densefd, sparsefd[1:nobs(densefd)]]) sg = sparsefit.grid - mf = StatsAPI.fit(FPCA, mixed_design; ncomponents=3, maxK=5, + mf = StatsAPI.fit(FPCA, mixed_design; ncomponents=3, max_k=5, grid=[nothing, sg], mean=[nothing, (; grid=sg, values=vec(sparsefit.mean.values))], covariance=[nothing, (; grid=sg, surface=sparsefit.covariance)], noise_variance=[nothing, .12]) @@ -312,11 +312,11 @@ end sim = simulate_kl(MersenneTwister(512),250;mean=t->2+t,eigenfunctions=funcs,eigenvalues=lambda, sampling,noise_sd=design===:dense ? 0. : .05) if design===:dense - f = StatsAPI.fit(FPCA,sim.data;ncomponents=2,maxK=2) + f = StatsAPI.fit(FPCA,sim.data;ncomponents=2,max_k=2) else phi = hcat([fun.(grid) for fun in funcs]...) # Supplied covariance isolates derivative/CE recovery from smoothing bias. - f = StatsAPI.fit(FPCA,sim.data;grid,ncomponents=2,maxK=2, + f = StatsAPI.fit(FPCA,sim.data;grid,ncomponents=2,max_k=2, mean=(;grid,values=2 .+ grid),covariance=(;grid,surface=phi*Diagonal(lambda)*transpose(phi)), noise_variance=.05^2) end diff --git a/test/face.jl b/test/face.jl index dc48c9d..c100911 100644 --- a/test/face.jl +++ b/test/face.jl @@ -24,7 +24,7 @@ import RecipesBase @test setting.P == quantity(ps, "penalty") @test setting.knots ≈ vec(quantity(ps, "knots")) rtol=1e-14 @test nbasis(setting.basis) == 38 - @test basismatrix(setting.basis, grid) ≈ setting.B + @test basis_matrix(setting.basis, grid) ≈ setting.B @test domain(setting.basis) == (0.0, 2.0) @test FDA._difference_penalty(5, 0) == Matrix{Float64}(I, 5, 5) @test FDA._difference_penalty(5, 1) ≈ diff(Matrix{Float64}(I, 5, 5); dims=1)' * diff(Matrix{Float64}(I, 5, 5); dims=1) @@ -147,7 +147,7 @@ import RecipesBase pace = StatsAPI.fit(FPCA, newdata; method=:pace, grid, mean=(; grid, values=vec(model.mean.values)), covariance=(; grid, surface=selected_covariance), - noise_variance=noise_variance(model), ncomponents=3, maxK=3) + noise_variance=noise_variance(model), ncomponents=3, max_k=3) @test scores(pace) ≈ prediction.scores rtol=1e-10 @test score_covariances(pace) ≈ prediction.covariances rtol=1e-10 @test StatsAPI.predict(pace, newdata).scores ≈ scores(pace) @@ -225,7 +225,7 @@ import RecipesBase for selector in (AIC(), BIC()) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; method=:face, ncomponents=selector) end - @test_throws ArgumentError StatsAPI.fit(FPCA, fd; method=:face, maxK=39) + @test_throws ArgumentError StatsAPI.fit(FPCA, fd; method=:face, max_k=39) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; method=:face, ncomponents=100) for bad in (0, -1, 1.5, :invalid, true) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; method=:face, ncomponents=bad) @@ -275,13 +275,13 @@ end data[2, 20:24] .= missing end fd = GridFunctionalData(data, grid) - f = StatsAPI.fit(FPCA, fd; method=:face, ncomponents=1, maxK=4, + f = StatsAPI.fit(FPCA, fd; method=:face, ncomponents=1, max_k=4, mean=(; grid, values=mu), face=(; knots=8, lambda=0.5)) @test length(f.all_eigenvalues) == 4 > ncomponents(f) @test noise_variance(f) > 0 d = derivative(f; method=:DPC, bandwidth=0.4, ncomponents=2) # Observation-space matrices provide an independent check of both the - # cross-covariance numerator and the full maxK covariance denominator. + # cross-covariance numerator and the full max_k covariance denominator. G = Matrix(f.settings.full_projected_covariance) L = FunctionalDataAnalysis.KernelSmoothing.smoothermatrix(grid, grid; bandwidth=0.4, derivative=1, degree=2, diff --git a/test/face_sparse.jl b/test/face_sparse.jl index b658aad..31371bf 100644 --- a/test/face_sparse.jl +++ b/test/face_sparse.jl @@ -55,7 +55,7 @@ using DelimitedFiles predicted, variance = Float64[], Float64[] for i in subjects t,_ = fd[i] - phi = basismatrix(full.covariance.basis,t)*full.settings.eigenfunction_coefficients + phi = basis_matrix(full.covariance.basis,t)*full.settings.eigenfunction_coefficients append!(predicted, FDA._face_sparse_mean_values(full.settings.mean,t)+phi*scores(full)[i,:]) append!(variance,diag(phi*score_covariances(full)[:,:,i]*phi')) end @@ -77,7 +77,7 @@ using DelimitedFiles # Independent observation-space BLUP, with exact splines at new times. for i in 1:nobs(heldout) t,y = heldout[i] - B = basismatrix(f.covariance.basis,t) + B = basis_matrix(f.covariance.basis,t) Phi = B*f.settings.eigenfunction_coefficients theta = f.covariance.coefficient_root*f.covariance.coefficient_root' V = B*theta*B' + noise_variance(f)*I @@ -150,7 +150,7 @@ using DelimitedFiles for selector in (AIC(),BIC(),0,true,:invalid) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;method=:face,ncomponents=selector) end - @test_throws ArgumentError StatsAPI.fit(FPCA,fd;method=:face,maxK=11) + @test_throws ArgumentError StatsAPI.fit(FPCA,fd;method=:face,max_k=11) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;method=:face,ncomponents=11) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;method=:face,noise_variance=1) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;method=:face,presmooth=:invalid) diff --git a/test/fpca.jl b/test/fpca.jl index e6bc5bd..0d9c56c 100644 --- a/test/fpca.jl +++ b/test/fpca.jl @@ -9,7 +9,7 @@ using DelimitedFiles n = Int(maximum(input[:, 1])) values = reduce(vcat, [transpose(Float64.(input[input[:, 1] .== i, 3])) for i in 1:n]) fd = GridFunctionalData(values, grid) - model = StatsAPI.fit(FPCA, fd; ncomponents=3, maxK=5) + model = StatsAPI.fit(FPCA, fd; ncomponents=3, max_k=5) reference(quantity) = Float64.(golden[golden[:, 1] .== quantity, 4]) relative_error(actual, expected) = norm(actual - expected) / norm(expected) @test relative_error(vec(model.mean.values), reference("mean")) < 1e-8 @@ -21,7 +21,7 @@ using DelimitedFiles # Reproduce the audited R convention explicitly for golden comparison. rectangle = fill(grid[2] - grid[1], length(grid)) compatible = FunctionalDataAnalysis._fpca_eigen(model.covariance, - rectangle, vec(model.mean.values); maxK=5, + rectangle, vec(model.mean.values); max_k=5, normalization_weights=fd.weights) @test relative_error(compatible.eigenvalues[1:3], reference("eigenvalue")) < 1e-8 @test relative_error(compatible.eigenvalues, reference("all_eigenvalue")) < 1e-8 @@ -56,7 +56,7 @@ using DelimitedFiles labels = ["subject-$i" for i in 1:6] fd = GridFunctionalData(values, grid; ids=labels) model = @inferred (data -> StatsAPI.fit(FPCA, data; - method=:dense, ncomponents=2, maxK=3))(fd) + method=:dense, ncomponents=2, max_k=3))(fd) @test model.method == :dense @test model.design == :regular @test model.grid == grid @@ -75,13 +75,13 @@ using DelimitedFiles @test (@inferred ncomponents(model)) == 2 @test (@inferred nobs(model)) == 6 @test ids(model) == labels - phi = basismatrix(eigenfunctions(model), grid) + phi = basis_matrix(eigenfunctions(model), grid) @test phi ≈ truebasis atol=2e-14 @test transpose(phi) * (weights .* phi) ≈ Matrix{Float64}(I, 2, 2) @test transpose(scores(model)) * scores(model) / (nobs(model) - 1) ≈ Diagonal(eigvals(model)) @test model.selection.rule == :fixed @test model.selection.chosen == 2 - @test model.selection.maxK == 3 + @test model.selection.max_k == 3 @test model.selection.threshold === nothing @test isempty(model.selection.scores) @test model.presmooth === nothing @@ -105,23 +105,23 @@ using DelimitedFiles @test_throws ArgumentError StatsAPI.predict(model, GridFunctionalData(values, grid .+ 0.001)) @test_throws DimensionMismatch StatsAPI.predict(model, GridFunctionalData(values[:, 1:4], grid[1:4])) - @test StatsAPI.fit(FPCA, fd; ncomponents=FVE(0.6), maxK=3).selection.chosen == 1 + @test StatsAPI.fit(FPCA, fd; ncomponents=FVE(0.6), max_k=3).selection.chosen == 1 boundary = fve(model)[1] - atboundary = StatsAPI.fit(FPCA, fd; ncomponents=FVE(boundary), maxK=3) - aboveboundary = StatsAPI.fit(FPCA, fd; ncomponents=FVE(nextfloat(boundary)), maxK=3) + atboundary = StatsAPI.fit(FPCA, fd; ncomponents=FVE(boundary), max_k=3) + aboveboundary = StatsAPI.fit(FPCA, fd; ncomponents=FVE(nextfloat(boundary)), max_k=3) @test ncomponents(atboundary) == 1 @test ncomponents(aboveboundary) == 2 @test atboundary.selection.rule == :FVE @test atboundary.selection.threshold == boundary @test atboundary.selection.scores == fve(atboundary) - @test ncomponents(StatsAPI.fit(FPCA, fd; ncomponents=FVE(1), maxK=3)) == 2 - capped = StatsAPI.fit(FPCA, fd; ncomponents=FVE(0.99), maxK=1) + @test ncomponents(StatsAPI.fit(FPCA, fd; ncomponents=FVE(1), max_k=3)) == 2 + capped = StatsAPI.fit(FPCA, fd; ncomponents=FVE(0.99), max_k=1) @test ncomponents(capped) == 1 @test capped.all_eigenvalues ≈ [2.0] @test fve(capped) == [1.0] default = StatsAPI.fit(FPCA, fd) @test default.selection.threshold == 0.99 - @test default.selection.maxK == 3 + @test default.selection.max_k == 3 for field in ("method", "design", "selection", "dropped", "sign_rule", "presmooth") @test occursin(field, sprint(show, model)) end @@ -149,38 +149,38 @@ using DelimitedFiles true_scores = signs .* transpose(sqrt.(eigenvalues .* (31 / 32))) fd = GridFunctionalData(2 .+ true_scores * transpose(truebasis), grid) model = StatsAPI.fit(FPCA, fd) - @test model.selection.maxK == 24 # min(m-2, n-2), without a cap of 20. + @test model.selection.max_k == 24 # min(m-2, n-2), without a cap of 20. @test model.all_eigenvalues ≈ eigenvalues rtol=1e-12 @test fve(model) ≈ cumsum(eigenvalues) ./ sum(eigenvalues) rtol=1e-12 @test ncomponents(model) == 24 - @test ncomponents(StatsAPI.fit(FPCA, fd; maxK=20)) == 20 + @test ncomponents(StatsAPI.fit(FPCA, fd; max_k=20)) == 20 end @testset "Eigenvalue cutoff, sign rules, and quadrature convergence" begin step = FunctionalDataAnalysis._fpca_eigen - aligned = step(Matrix(Diagonal([3.0, 1.0])), [1.0, 1.0], [-1.0, 1.0]; maxK=2) + aligned = step(Matrix(Diagonal([3.0, 1.0])), [1.0, 1.0], [-1.0, 1.0]; max_k=2) @test aligned.eigenvalues == [3.0, 1.0] @test aligned.eigenfunctions ≈ [-1.0 0.0; 0.0 1.0] - fallback = step(Matrix(Diagonal([3.0, 1.0])), [1.0, 1.0], zeros(2); maxK=2) + fallback = step(Matrix(Diagonal([3.0, 1.0])), [1.0, 1.0], zeros(2); max_k=2) @test fallback.eigenfunctions ≈ Matrix{Float64}(I, 2, 2) # Mean inner products below the stated tolerance also use the fallback. tinymean = step(Matrix(Diagonal([3.0, 1.0])), [1.0, 1.0], [-1e-15, 1.0]; - maxK=2, sign_tol=1e-12) + max_k=2, sign_tol=1e-12) @test tinymean.eigenfunctions ≈ Matrix{Float64}(I, 2, 2) projected = step(Matrix(Diagonal([3.0, -2.0, 1e-12])), ones(3), ones(3); - maxK=3, eigen_tol=1e-10) + max_k=3, eigen_tol=1e-10) @test projected.eigenvalues == [3.0] @test projected.fve == [1.0] @test projected.dropped.number == 1 @test projected.dropped.total == 2.0 @test projected.dropped.share ≈ 2 / (5 + 1e-12) - @test_throws ArgumentError step(-Matrix{Float64}(I, 2, 2), ones(2), ones(2); maxK=2) + @test_throws ArgumentError step(-Matrix{Float64}(I, 2, 2), ones(2), ones(2); max_k=2) # Finite inputs can have a spectrum whose direct sum overflows. - large = step(Matrix(Diagonal([1e308, 9e307])), ones(2), ones(2); maxK=2) + large = step(Matrix(Diagonal([1e308, 9e307])), ones(2), ones(2); max_k=2) @test large.fve[1] ≈ 1 / 1.9 @test large.fve[end] == 1.0 - indefinite = step(Matrix(Diagonal([1e308, -1e308])), ones(2), ones(2); maxK=2) + indefinite = step(Matrix(Diagonal([1e308, -1e308])), ones(2), ones(2); max_k=2) @test indefinite.dropped.total == 1e308 @test indefinite.dropped.share ≈ 0.5 @@ -188,7 +188,7 @@ using DelimitedFiles # Their squares overflow, but their weighted norms remain exactly one. tinyweights = fill(1e-310, 2) small_weighted = step(Matrix(Diagonal([1e308, 5e307])), tinyweights, - ones(2); maxK=2) + ones(2); max_k=2) @test all(isfinite, small_weighted.eigenfunctions) @test all(k -> maximum(abs, small_weighted.eigenfunctions[:, k]) > 0, 1:2) rootweights = sqrt.(tinyweights) @@ -200,11 +200,11 @@ using DelimitedFiles # Relative mean alignment is unchanged by finite rescaling, even when # the unscaled alignment and product of norms would overflow. scaled_sign = step(Matrix(Diagonal([1e308, 5e307])), tinyweights, - [-1e308, 1e308]; maxK=2) + [-1e308, 1e308]; max_k=2) @test scaled_sign.eigenfunctions[1, 1] < 0 @test scaled_sign.eigenfunctions[2, 2] > 0 scaled_fallback = step(Matrix(Diagonal([1e308, 5e307])), tinyweights, - [-1e294, 1e308]; maxK=2, sign_tol=1e-12) + [-1e294, 1e308]; max_k=2, sign_tol=1e-12) @test scaled_fallback.eigenfunctions[1, 1] > 0 errors = Float64[] @@ -212,9 +212,9 @@ using DelimitedFiles grid = collect(range(0, 1; length=m)) trap = quadrature_weights(grid) surface = ones(m, m) # Kernel G(s,t)=1 has eigenvalue exactly 1. - ours = step(surface, trap, ones(m); maxK=1) + ours = step(surface, trap, ones(m); max_k=1) rectangle = fill(1 / (m - 1), m) - compatible = step(surface, rectangle, ones(m); maxK=1, normalization_weights=trap) + compatible = step(surface, rectangle, ones(m); max_k=1, normalization_weights=trap) @test only(ours.eigenvalues) ≈ 1.0 @test only(compatible.eigenvalues) ≈ 1 + 1 / (m - 1) @test compatible.eigenfunctions[:, 1] ≈ ones(m) @@ -236,7 +236,7 @@ using DelimitedFiles # Six score-variance standard errors; exact RNG values are not assumed. @test all(abs.(eigvals(model) - eigenvalues) .< 6sqrt(2 / (5000 - 1)) .* eigenvalues) - phi = basismatrix(eigenfunctions(model), grid) + phi = basis_matrix(eigenfunctions(model), grid) truephi = [f(t) for t in grid, f in functions] # The constant has positive population mean alignment; the sine has # negative alignment with 2+t, so the mean-alignment sign rule fixes both. @@ -309,10 +309,10 @@ using DelimitedFiles for count in (0, -1, 1.5, 4, AIC(), BIC()) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; ncomponents=count) end - @test_throws ArgumentError StatsAPI.fit(FPCA, fd; ncomponents=2, maxK=1) - @test_throws ArgumentError StatsAPI.fit(FPCA, fd; ncomponents=2, maxK=3) # Positive rank is one. + @test_throws ArgumentError StatsAPI.fit(FPCA, fd; ncomponents=2, max_k=1) + @test_throws ArgumentError StatsAPI.fit(FPCA, fd; ncomponents=2, max_k=3) # Positive rank is one. for cap in (0, -1, 1.5, 4) - @test_throws ArgumentError StatsAPI.fit(FPCA, fd; maxK=cap) + @test_throws ArgumentError StatsAPI.fit(FPCA, fd; max_k=cap) end for threshold in (0, -0.1, 1.1, NaN, Inf) @test_throws ArgumentError FVE(threshold) @@ -325,15 +325,15 @@ using DelimitedFiles @test_throws ArgumentError StatsAPI.fit(FPCA, GridFunctionalData(ones(3, 5), fd.grid)) twopoint = GridFunctionalData([1.0 2.0; 2.0 3.0; 3.0 4.0], [0.0, 1.0]) twopointfit = StatsAPI.fit(FPCA, twopoint) - @test twopointfit.selection.maxK == 1 # Clamp the m-2 default to one. + @test twopointfit.selection.max_k == 1 # Clamp the m-2 default to one. @test eigvals(twopointfit) ≈ [1.0] @test reconstruct(twopointfit).values ≈ twopoint.values twocurvefit = StatsAPI.fit(FPCA, fd[1:2]) - @test twocurvefit.selection.maxK == 1 # Clamp the n-2 default to one. + @test twocurvefit.selection.max_k == 1 # Clamp the n-2 default to one. @test ncomponents(twocurvefit) == 1 - @test StatsAPI.fit(FPCA, fd[1:3]).selection.maxK == 1 - # Explicit maxK can still use the sample covariance's n-1 rank cap. - @test StatsAPI.fit(FPCA, fd[1:3]; maxK=2).selection.maxK == 2 + @test StatsAPI.fit(FPCA, fd[1:3]).selection.max_k == 1 + # Explicit max_k can still use the sample covariance's n-1 rank cap. + @test StatsAPI.fit(FPCA, fd[1:3]; max_k=2).selection.max_k == 2 masked = GridFunctionalData([1.0 missing 3.0; 2.0 3.0 4.0], [0.0, 0.5, 1.0]) irregular = to_irregular(fd) for data in (fd, masked, irregular) diff --git a/test/fpcbasis.jl b/test/fpcbasis.jl index 82576c8..94fe00f 100644 --- a/test/fpcbasis.jl +++ b/test/fpcbasis.jl @@ -8,9 +8,9 @@ @test eltype(b) == Float64 @test b.weights == [0.125, 0.5, 0.375] points = [1.0, 0.125, 0.0, 0.625, 0.25] - @test (@inferred basismatrix(b, points)) == hcat(ones(5), points) - @test basismatrix(b, grid; derivative=0) == values - @test size(@inferred basismatrix(b, Float64[])) == (0, 2) + @test (@inferred basis_matrix(b, points)) == hcat(ones(5), points) + @test basis_matrix(b, grid; derivative=0) == values + @test size(@inferred basis_matrix(b, Float64[])) == (0, 2) @test (@inferred gram(b)) ≈ [1.0 0.5; 0.5 1 / 3] @test (@inferred penalty(b, 0)) == gram(b) @@ -31,11 +31,11 @@ single = @inferred FPCBasis(Float32[0, 0.25, 1], Float32[1 0; 1 0.25; 1 1]) @test eltype(single) == Float32 - @test (@inferred basismatrix(single, Float32[0, 0.5, 1])) isa Matrix{Float32} + @test (@inferred basis_matrix(single, Float32[0, 0.5, 1])) isa Matrix{Float32} @test (@inferred gram(single)) isa Matrix{Float32} @test (@inferred penalty(single, 0)) isa Matrix{Float32} @test eltype(@inferred FPCBasis([0, 1], [1 0; 1 1])) == Float64 - @test basismatrix(single, [0.5]) isa Matrix{Float64} + @test basis_matrix(single, [0.5]) isa Matrix{Float64} fd = BasisFunctionalData([2.0 3.0], b) @test evaluate(fd, [0.0, 0.4, 1.0]) ≈ [2.0 3.2 5.0] @@ -99,10 +99,10 @@ @test_throws ArgumentError FPCBasis(grid, invalid) end for points in ([-0.1], [1.1], [NaN], [Inf]) - @test_throws ArgumentError basismatrix(b, points) + @test_throws ArgumentError basis_matrix(b, points) end for order in (-1, 0.5, 1, 2) - @test_throws ArgumentError basismatrix(b, [0.5]; derivative=order) + @test_throws ArgumentError basis_matrix(b, [0.5]; derivative=order) @test_throws ArgumentError penalty(b, order) end end diff --git a/test/kernelsmoothing_agreement.jl b/test/kernelsmoothing_agreement.jl index ff26f84..debdcb3 100644 --- a/test/kernelsmoothing_agreement.jl +++ b/test/kernelsmoothing_agreement.jl @@ -2,8 +2,8 @@ import KernelSmoothing import StatsAPI const FDA = FunctionalDataAnalysis -# The smoothers FDA used before it switched to KernelSmoothing's localpolynomial, -# smoothermatrix and bivariate smoother (FDA 539b944). The new results must agree +# The smoothers FDA used before it switched to KernelSmoothing's local_polynomial, +# smoother_matrix and bivariate smoother (FDA 539b944). The new results must agree # with these to 1e-10 for the same kernel and bandwidth. module OldSmoothers using LinearAlgebra diff --git a/test/mean.jl b/test/mean.jl index 6dbc0e9..1ec3928 100644 --- a/test/mean.jl +++ b/test/mean.jl @@ -139,7 +139,7 @@ end for weighting in (:obs, :subj) w = vcat([fill(weighting == :obs ? 1. : 1/length(t), length(t)) for t in times]...) for rule in (KernelSmoothing.RuleOfThumb(), KernelSmoothing.IMSEOptimal()) - direct = KernelSmoothing.selectbandwidth(rule, + direct = KernelSmoothing.select_bandwidth(rule, KernelSmoothing.SmoothingData(x, y; weights=w), kernel; degree=1) # KernelSmoothing selectors' `details` types depend on the data (and on # the KernelSmoothing version), so selector fits are not inferred. @@ -157,7 +157,7 @@ end @test_throws ArgumentError("CurveCV multipliers must be finite and positive") CurveCV(; pilot, multipliers=[1., -2.]) multipliers[1] = 10. m = StatsAPI.fit(SmoothedMean, fd; grid, bandwidth=cv, kernel, weighting) - h = KernelSmoothing.selectbandwidth(pilot, + h = KernelSmoothing.select_bandwidth(pilot, KernelSmoothing.SmoothingData(x, y; weights=w), kernel; degree=1).h candidates = [1., 2., 4.] .* h explicit = StatsAPI.fit(SmoothedMean, fd; grid, kernel, weighting, diff --git a/test/mfpca.jl b/test/mfpca.jl index d3efc50..5dd3d5b 100644 --- a/test/mfpca.jl +++ b/test/mfpca.jl @@ -50,7 +50,7 @@ end for mode in ("given_unit", "given_weighted", "cross_method") weights = mode == "given_unit" ? ones(3) : [0.7, 1.5, 2.2] fitted = @inferred StatsAPI.fit(MFPCA, fixture.data; ncomponents=4, - univariate=(; ncomponents=2, maxK=2), weights) + univariate=(; ncomponents=2, max_k=2), weights) # Strict R values compare signed coefficients as well as functions # and scores. Cross-method signs align the whole multivariate function. signs = ones(4) @@ -103,14 +103,14 @@ end end for (j,g) in enumerate(grids)] data = MultivariateFunctionalData(nodes) one = StatsAPI.fit(MFPCA, MultivariateFunctionalData(nodes[1:1]); - ncomponents=2, univariate=(; ncomponents=2, maxK=2)) - uni = StatsAPI.fit(FPCA, nodes[1]; ncomponents=2, maxK=2) + ncomponents=2, univariate=(; ncomponents=2, max_k=2)) + uni = StatsAPI.fit(FPCA, nodes[1]; ncomponents=2, max_k=2) @test eigvals(one) ≈ eigvals(uni) rtol=1e-12 @test scores(one) ≈ scores(uni) atol=1e-13 @test eigenfunctions(one)[1].values ≈ eigenfunctions(uni).values atol=1e-13 fitted = @inferred StatsAPI.fit(MFPCA, data; ncomponents=4, - univariate=(; ncomponents=2, maxK=2)) + univariate=(; ncomponents=2, max_k=2)) @test (@inferred ncomponents(fitted)) == 4 @test (@inferred nobs(fitted)) == n @test (@inferred scores(fitted)) isa Matrix{Float64} @@ -127,7 +127,7 @@ end @test eigenfunctions(fitted)[j].grid == grids[j] end weighted = StatsAPI.fit(MFPCA, data; ncomponents=4, weights=[2.0, 3.0], - univariate=(; ncomponents=2, maxK=2)) + univariate=(; ncomponents=2, max_k=2)) @test eigvals(weighted) ≈ [18.0, 8.0, 6.0, 1.5] rtol=1e-12 geometry = sum(weighted.weights[j] .* (b.values' * (b.weights .* b.values)) for (j,b) in enumerate(eigenfunctions(weighted))) @@ -139,7 +139,7 @@ end end # Global weight scaling rescales eigenvalues and scores, preserving fits. scaled = StatsAPI.fit(MFPCA, data; ncomponents=4, weights=[14.0, 21.0], - univariate=(; ncomponents=2, maxK=2)) + univariate=(; ncomponents=2, max_k=2)) @test eigvals(scaled) ≈ 7eigvals(weighted) rtol=1e-12 @test scores(scaled) ≈ sqrt(7)*scores(weighted) atol=1e-12 for j in 1:2 @@ -150,7 +150,7 @@ end rescaled_data = MultivariateFunctionalData([nodes[1], GridFunctionalData(5nodes[2].values, grids[2]; ids=labels)]) units = StatsAPI.fit(MFPCA, rescaled_data; ncomponents=4, weights=[2.0, 3/25], - univariate=(; ncomponents=2, maxK=2)) + univariate=(; ncomponents=2, max_k=2)) @test eigvals(units) ≈ eigvals(weighted) rtol=1e-12 @test scores(units) ≈ scores(weighted) atol=1e-12 @@ -181,7 +181,7 @@ end sparse = IrregularFunctionalData(times, [fill(5+cos(i), length(t)) for (i,t) in enumerate(times)]; domain=(2,3), ids=ids(dense)) data = MultivariateFunctionalData([dense, sparse]) - options = (; ncomponents=[2, 1], maxK=[2, 1], grid=[nothing, g], + options = (; ncomponents=[2, 1], max_k=[2, 1], grid=[nothing, g], mean=[nothing, (; grid=g, values=zeros(5))], covariance=[nothing, (; grid=g, surface=ones(5, 5))], noise_variance=[nothing, 0.2]) @@ -235,7 +235,7 @@ end data64 = mfpca_fixture().data data = MultivariateFunctionalData([GridFunctionalData(Float32.(component(data64,j).values), Float32.(component(data64,j).grid)) for j in 1:3]) - fitted = StatsAPI.fit(MFPCA, data; ncomponents=3, univariate=(; ncomponents=2, maxK=2)) + fitted = StatsAPI.fit(MFPCA, data; ncomponents=3, univariate=(; ncomponents=2, max_k=2)) @test fitted isa MFPCA{Float32} @test scores(fitted) isa Matrix{Float32} @test eigenfunctions(fitted) isa Vector{FPCBasis{Float32}} @@ -258,14 +258,14 @@ end Float64.(component(data,j).grid)) for j in 1:3]) @test StatsAPI.predict(fitted, promoted) isa Matrix{Float64} weighted = StatsAPI.fit(MFPCA, data; ncomponents=3, weights=ones(Float64, 3), - univariate=(; ncomponents=2, maxK=2)) + univariate=(; ncomponents=2, max_k=2)) @test weighted isa MFPCA{Float64,Float32} @test eigenfunctions(weighted) isa Vector{FPCBasis{Float64}} end @testset "Validation and copied weights" begin data = mfpca_fixture().data - options = (; ncomponents=2, maxK=2) + options = (; ncomponents=2, max_k=2) @test_throws UndefKeywordError StatsAPI.fit(MFPCA, data) for count in (0, -1, 1.5, true, FVE(0.9)) @test_throws ArgumentError StatsAPI.fit(MFPCA, data; ncomponents=count) @@ -290,9 +290,9 @@ end # More requested components than the joint positive rank must throw. duplicate = MultivariateFunctionalData(fill(component(data, 1), 2)) @test_throws ArgumentError StatsAPI.fit(MFPCA, duplicate; ncomponents=3, - univariate=(; ncomponents=1, maxK=2)) + univariate=(; ncomponents=1, max_k=2)) @test_throws ArgumentError StatsAPI.fit(MFPCA, duplicate; ncomponents=2, - univariate=(; ncomponents=1, maxK=2)) + univariate=(; ncomponents=1, max_k=2)) weights = [0.7, 1.5, 2.2] fitted = StatsAPI.fit(MFPCA, data; ncomponents=3, univariate=options, weights) weights[1] = 10 diff --git a/test/multivariate_fpca.jl b/test/multivariate_fpca.jl index c0f34f8..e7e926b 100644 --- a/test/multivariate_fpca.jl +++ b/test/multivariate_fpca.jl @@ -22,7 +22,7 @@ end fixture = m4_dense_fixture() data, grid = fixture.data, fixture.grid model = @inferred (fd -> StatsAPI.fit(FPCA, fd; - basis=:marginal, ncomponents=2, maxK=2))(data) + basis=:marginal, ncomponents=2, max_k=2))(data) @test model isa MultivariateFPCA @test model.basis == :marginal @test (@inferred ncomponents(model)) == [2, 2] @@ -33,7 +33,7 @@ end @test (@inferred eigenfunctions(model)) isa Vector{FPCBasis{Float64}} @test (@inferred score_covariances(model)) == [nothing, nothing] for j in 1:2 - univariate = StatsAPI.fit(FPCA, component(data, j); ncomponents=2, maxK=2) + univariate = StatsAPI.fit(FPCA, component(data, j); ncomponents=2, max_k=2) node = component(model, j) @test node isa FPCA @test node.method == univariate.method == :dense @@ -85,11 +85,11 @@ end for j in 1:2] presmoothers = [LocalPolynomial(; bandwidth=0.4), LocalPolynomial(; bandwidth=0.6)] - presmoothed = StatsAPI.fit(FPCA, data; ncomponents=2, maxK=2, + presmoothed = StatsAPI.fit(FPCA, data; ncomponents=2, max_k=2, presmooth=presmoothers) for j in 1:2 reference = StatsAPI.fit(FPCA, component(data, j); - ncomponents=2, maxK=2, presmooth=presmoothers[j]) + ncomponents=2, max_k=2, presmooth=presmoothers[j]) @test scores(presmoothed, j) == scores(reference) @test eigvals(presmoothed)[j] == eigvals(reference) end @@ -102,7 +102,7 @@ end @test newprediction[j] ≈ StatsAPI.predict(component(model, j), component(shifted, j)) end - unequal = StatsAPI.fit(FPCA, data; ncomponents=[1, 2], maxK=2) + unequal = StatsAPI.fit(FPCA, data; ncomponents=[1, 2], max_k=2) @test ncomponents(unequal) == [1, 2] error = try score_array(unequal) @@ -158,7 +158,7 @@ end supplied_covariance = (; grid=copy(grid), surface=Matrix(Symmetric(functions * Diagonal(lambda) * functions'))) options = (; method=[:auto, :pace], grid=[nothing, grid], - ncomponents=[1, 2], maxK=[2, 2], mean=[nothing, supplied_mean], + ncomponents=[1, 2], max_k=[2, 2], mean=[nothing, supplied_mean], covariance=[nothing, supplied_covariance], noise_variance=[nothing, 0.13], eigen_tol=[1e-10, 1e-9], sign_tol=[1e-12, 2e-12]) model = StatsAPI.fit(FPCA, data; options...) @@ -169,9 +169,9 @@ end @test component(model, 2).grid == grid @test ids(model) == fixture.labels references = [StatsAPI.fit(FPCA, component(data, 1); ncomponents=1, - maxK=2, eigen_tol=1e-10, sign_tol=1e-12), + max_k=2, eigen_tol=1e-10, sign_tol=1e-12), StatsAPI.fit(FPCA, component(data, 2); method=:pace, grid, - ncomponents=2, maxK=2, mean=supplied_mean, + ncomponents=2, max_k=2, mean=supplied_mean, covariance=supplied_covariance, noise_variance=0.13, eigen_tol=1e-9, sign_tol=2e-12)] for j in 1:2 @@ -197,10 +197,10 @@ end # Shared values and explicit per-node values resolve identically. sparsepair = MultivariateFunctionalData([sparse, sparse]) - shared = StatsAPI.fit(FPCA, sparsepair; grid, ncomponents=2, maxK=2, + shared = StatsAPI.fit(FPCA, sparsepair; grid, ncomponents=2, max_k=2, mean=supplied_mean, covariance=supplied_covariance, noise_variance=0.13) separate = StatsAPI.fit(FPCA, sparsepair; grid=[grid, grid], - ncomponents=[2, 2], maxK=[2, 2], mean=[supplied_mean, supplied_mean], + ncomponents=[2, 2], max_k=[2, 2], mean=[supplied_mean, supplied_mean], covariance=[supplied_covariance, supplied_covariance], noise_variance=[0.13, 0.13]) @test scores(shared) == scores(separate) @test eigvals(shared) == eigvals(separate) @@ -211,13 +211,13 @@ end kernels = [KernelSmoothing.GaussianKernel(), KernelSmoothing.GaussianKernel()] weightings = [:obs, :subj] binnings = [nothing, 5] - smoothed = StatsAPI.fit(FPCA, sparsepair; grid, ncomponents=1, maxK=2, + smoothed = StatsAPI.fit(FPCA, sparsepair; grid, ncomponents=1, max_k=2, mean_bandwidth=mean_bandwidths, covariance_bandwidth=covariance_bandwidths, diagonal_bandwidth=diagonal_bandwidths, kernel=kernels, weighting=weightings, binning=binnings, noise_variance=0.13) for j in 1:2 reference = StatsAPI.fit(FPCA, sparse; - grid, ncomponents=1, maxK=2, mean_bandwidth=mean_bandwidths[j], + grid, ncomponents=1, max_k=2, mean_bandwidth=mean_bandwidths[j], covariance_bandwidth=covariance_bandwidths[j], diagonal_bandwidth=diagonal_bandwidths[j], kernel=kernels[j], weighting=weightings[j], binning=binnings[j], noise_variance=0.13) @@ -235,7 +235,7 @@ end fixture = m4_dense_fixture() data, grid = fixture.data, fixture.grid model = @inferred (fd -> StatsAPI.fit(FPCA, fd; - basis=:pooled, ncomponents=2, maxK=2))(data) + basis=:pooled, ncomponents=2, max_k=2))(data) @test model isa MultivariateFPCA @test model.basis == :pooled @test (@inferred ncomponents(model)) == [2, 2] @@ -244,7 +244,7 @@ end @test (@inferred score_covariances(model)) == [nothing, nothing] functions = @inferred eigenfunctions(model) @test functions[1] == functions[2] - @test basismatrix(functions[1], grid) ≈ fixture.functions atol=2e-14 + @test basis_matrix(functions[1], grid) ≈ fixture.functions atol=2e-14 @test functions[1].weights == fixture.weights @test functions[1].values' * (fixture.weights .* functions[1].values) ≈ I atol=2e-14 @test eigvals(model)[1] ≈ [2.0, 1.2] @@ -281,17 +281,17 @@ end end default = StatsAPI.fit(FPCA, data; basis=:pooled) @test ncomponents(default) == [2, 2] - selected = StatsAPI.fit(FPCA, data; basis=:pooled, ncomponents=FVE(0.85), maxK=2) + selected = StatsAPI.fit(FPCA, data; basis=:pooled, ncomponents=FVE(0.85), max_k=2) @test ncomponents(selected) == [1, 1] - capped = StatsAPI.fit(FPCA, data; basis=:pooled, maxK=1) - # Pooled rank is nodes*(n-1), so maxK may exceed n-1 when the grid allows. + capped = StatsAPI.fit(FPCA, data; basis=:pooled, max_k=1) + # Pooled rank is nodes*(n-1), so max_k may exceed n-1 when the grid allows. tgrid = collect(range(0, 1; length=12)) few = MultivariateFunctionalData([GridFunctionalData( permutedims(reduce(hcat, [sin.(k * (j + 1) .* tgrid) for k in 1:3])), tgrid) for j in 1:3]) wide = StatsAPI.fit(FPCA, few; basis=:pooled) - @test wide.settings.maxK == min(12 - 2, 3 * 2 - 1) - @test StatsAPI.fit(FPCA, few; basis=:pooled, maxK=6).settings.maxK == 6 - @test_throws ArgumentError StatsAPI.fit(FPCA, few; basis=:pooled, maxK=7) + @test wide.settings.max_k == min(12 - 2, 3 * 2 - 1) + @test StatsAPI.fit(FPCA, few; basis=:pooled, max_k=6).settings.max_k == 6 + @test_throws ArgumentError StatsAPI.fit(FPCA, few; basis=:pooled, max_k=7) @test ncomponents(capped) == [1, 1] error = try component(model, 1) @@ -307,7 +307,7 @@ end # Large differences in means must not create another pooled component. mean_shifted = MultivariateFunctionalData([component(data, 1), GridFunctionalData(fixture.values[2] .+ 1000, grid; ids=fixture.labels)]) - shifted_model = StatsAPI.fit(FPCA, mean_shifted; basis=:pooled, ncomponents=2, maxK=2) + shifted_model = StatsAPI.fit(FPCA, mean_shifted; basis=:pooled, ncomponents=2, max_k=2) @test eigvals(shifted_model) ≈ eigvals(model) rtol=1e-12 @test abs.(scores(shifted_model, 2)) ≈ abs.(scores(model, 2)) atol=3e-12 @@ -318,8 +318,8 @@ end GridFunctionalData(reversed_scores * fixture.functions' .+ fixture.means[2]', grid; ids=fixture.labels)]) reversed_model = StatsAPI.fit(FPCA, reversed_data; - basis=:pooled, ncomponents=2, maxK=2) - @test basismatrix(eigenfunctions(reversed_model)[1], grid) ≈ + basis=:pooled, ncomponents=2, max_k=2) + @test basis_matrix(eigenfunctions(reversed_model)[1], grid) ≈ fixture.functions[:, [2, 1]] atol=3e-14 @test eigvals(reversed_model)[1] ≈ [1.2, 2.0] @test eigvals(reversed_model)[2] ≈ [4.8, 0.125] @@ -346,8 +346,8 @@ end (randn(rng, n, 2) .* sqrt.(variances[j])') * functions' .+ (j .+ grid)', grid) for j in 1:2] model = StatsAPI.fit(FPCA, MultivariateFunctionalData(nodes); - basis=:pooled, ncomponents=2, maxK=2) - recovered = basismatrix(eigenfunctions(model)[1], grid) + basis=:pooled, ncomponents=2, max_k=2) + recovered = basis_matrix(eigenfunctions(model)[1], grid) # Six score-variance standard errors and a one-percent alignment # tolerance; neither check depends on exact Julia RNG output. expected = (variances[1] + variances[2]) / 2 @@ -361,7 +361,7 @@ end data = MultivariateFunctionalData([GridFunctionalData(Float32.(values), Float32.(fixture.grid)) for values in fixture.values]) for basis in (:marginal, :pooled) - model = StatsAPI.fit(FPCA, data; basis, ncomponents=2, maxK=2) + model = StatsAPI.fit(FPCA, data; basis, ncomponents=2, max_k=2) @test scores(model) isa Vector{Matrix{Float32}} @test eigvals(model) isa Vector{Vector{Float32}} @test eigenfunctions(model) isa Vector{FPCBasis{Float32}} @@ -408,12 +408,12 @@ end mean64 = (; grid=Float64.(pace_grid), values=1 .+ Float64.(pace_grid)) covariance64 = (; grid=Float64.(pace_grid), surface=ones(5, 5)) sigma2 = 0.123456789012345 - model = StatsAPI.fit(FPCA, mixed; ncomponents=1, maxK=[2, 1], + model = StatsAPI.fit(FPCA, mixed; ncomponents=1, max_k=[2, 1], grid=[nothing, pace_grid], mean=[nothing, mean64], covariance=[nothing, covariance64], noise_variance=[nothing, sigma2]) - references = [StatsAPI.fit(FPCA, component(mixed, 1); ncomponents=1, maxK=2), + references = [StatsAPI.fit(FPCA, component(mixed, 1); ncomponents=1, max_k=2), StatsAPI.fit(FPCA, component(mixed, 2); grid=pace_grid, - ncomponents=1, maxK=1, mean=mean64, + ncomponents=1, max_k=1, mean=mean64, covariance=covariance64, noise_variance=sigma2)] @test model isa MultivariateFPCA{Float64} @test scores(component(model, 1)) isa Matrix{Float32} @@ -440,7 +440,7 @@ end fixture = m4_dense_fixture() data, grid = fixture.data, fixture.grid @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:invalid) - for keyword in (:method, :ncomponents, :maxK, :presmooth, :mean, + for keyword in (:method, :ncomponents, :max_k, :presmooth, :mean, :covariance, :noise_variance, :mean_bandwidth, :covariance_bandwidth, :diagonal_bandwidth, :kernel, :weighting, :binning, :eigen_tol, :sign_tol) @@ -464,7 +464,7 @@ end end @test_throws ArgumentError StatsAPI.fit(FPCA, data[1:1]; basis=:pooled) @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, ncomponents=0) - @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, maxK=0) + @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, max_k=0) @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, eigen_tol=NaN) @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, sign_tol=-1) @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, @@ -477,7 +477,7 @@ end end @test pace_error isa ArgumentError @test occursin("pooled PACE", sprint(showerror, pace_error)) - for (name, value) in pairs((; ncomponents=[1, 1], maxK=[2, 2], + for (name, value) in pairs((; ncomponents=[1, 1], max_k=[2, 2], eigen_tol=[1e-10, 1e-10], sign_tol=[1e-12, 1e-12])) @test_throws ArgumentError StatsAPI.fit(FPCA, data; basis=:pooled, NamedTuple{(name,)}((value,))...) @@ -486,14 +486,14 @@ end ids=fixture.labels) basis_data = MultivariateFunctionalData([component(data, 1), basis_node]) basis_error = try - StatsAPI.fit(FPCA, basis_data; ncomponents=1, maxK=2) + StatsAPI.fit(FPCA, basis_data; ncomponents=1, max_k=2) catch exception exception end @test basis_error isa ArgumentError @test occursin("basis data", sprint(showerror, basis_error)) @test occursin("convert", sprint(showerror, basis_error)) - model = StatsAPI.fit(FPCA, data; basis=:pooled, ncomponents=2, maxK=2) + model = StatsAPI.fit(FPCA, data; basis=:pooled, ncomponents=2, max_k=2) @test_throws DimensionMismatch StatsAPI.predict(model, MultivariateFunctionalData([component(data, 1)])) @test_throws ArgumentError StatsAPI.predict(model, different_grid) diff --git a/test/multivariate_fpca_golden.jl b/test/multivariate_fpca_golden.jl index 3e4bcde..ffffcc6 100644 --- a/test/multivariate_fpca_golden.jl +++ b/test/multivariate_fpca_golden.jl @@ -34,7 +34,7 @@ using DelimitedFiles for j in 1:p expected_functions = reshape(reference(golden, string(basis), "eigenfunction", j), length(grid), 2) - actual_functions = basismatrix(eigenfunctions(fitted)[j], grid) + actual_functions = basis_matrix(eigenfunctions(fitted)[j], grid) signs = aligned_signs(actual_functions, expected_functions) @test all(abs.(signs) .== 1) @test relative_error(actual_functions .* transpose(signs), @@ -75,7 +75,7 @@ using DelimitedFiles noise = [only(reference(golden, "input", "noise_variance", j)) for j in 1:p] # Match FGM's full positive-rank denominator and relative cutoff. fitted = StatsAPI.fit(FPCA, MultivariateFunctionalData(nodes); - basis=:marginal, grid, ncomponents=2, maxK=length(grid), + basis=:marginal, grid, ncomponents=2, max_k=length(grid), mean_bandwidth=0.3, covariance_bandwidth=0.3, diagonal_bandwidth=0.3, noise_variance=noise, kernel=KernelSmoothing.GaussianKernel(), weighting=:subj, binning=nothing, eigen_tol=100eps(Float64)) @@ -88,7 +88,7 @@ using DelimitedFiles reference(golden, "marginal", "mean", j)) < 1e-8 expected_functions = reshape(reference(golden, "marginal", "eigenfunction", j), length(grid), 2) - actual_functions = basismatrix(eigenfunctions(fitted)[j], grid) + actual_functions = basis_matrix(eigenfunctions(fitted)[j], grid) signs = aligned_signs(actual_functions, expected_functions) @test all(abs.(signs) .== 1) @test relative_error(actual_functions .* transpose(signs), expected_functions) < 1e-8 diff --git a/test/naming.jl b/test/naming.jl new file mode 100644 index 0000000..a106f22 --- /dev/null +++ b/test/naming.jl @@ -0,0 +1,37 @@ +import StatsAPI + +@testset "canonical names and compatibility" begin + for (old, canonical) in ((:basismatrix, :basis_matrix), (:modesplot, :modes_plot), + (:screeplot, :scree_plot), (:covplot, :cov_plot), + (:designplot, :design_plot), (:fitplot, :fit_plot), + (:bandplot, :band_plot)) + @test getproperty(FunctionalDataAnalysis, old) === + getproperty(FunctionalDataAnalysis, canonical) + end + grid = collect(range(0, 1; length=7)) + values = [sin(i * t) + i * t for i in 1:6, t in grid] + fd = GridFunctionalData(values, grid) + canonical = StatsAPI.fit(FPCA, fd; ncomponents=1, max_k=2) + legacy = StatsAPI.fit(FPCA, fd; ncomponents=1, maxK=2) + @test scores(legacy) == scores(canonical) + @test eigvals(legacy) == eigvals(canonical) + @test legacy.settings.maxK == legacy.settings.max_k == 2 + @test legacy.selection.maxK == legacy.selection.max_k == 2 + @test_throws ArgumentError StatsAPI.fit(FPCA, fd; max_k=2, maxK=2) + b = FourierBasis((0, 1), 3) + basis_fd = BasisFunctionalData([(-1)^i * (i + j) for i in 1:6, j in 1:3], b) + basis_new = StatsAPI.fit(FPCA, basis_fd; ncomponents=1, max_k=2) + basis_old = StatsAPI.fit(FPCA, basis_fd; ncomponents=1, maxK=2) + @test scores(basis_old) == scores(basis_new) + @test basis_old.settings.maxK == basis_old.settings.max_k == 2 + @test_throws ArgumentError StatsAPI.fit(FPCA, basis_fd; max_k=2, maxK=2) + cross_new = StatsAPI.fit(FSVD, fd, fd; ncomponents=1, max_k=2) + cross_old = StatsAPI.fit(FSVD, fd, fd; ncomponents=1, maxK=2) + @test cross_old.left_scores == cross_new.left_scores + @test cross_old.selection.maxK == cross_old.selection.max_k == 2 + @test_throws ArgumentError StatsAPI.fit(FSVD, fd, fd; max_k=2, maxK=2) + mfd = MultivariateFunctionalData([fd, fd]) + pooled = StatsAPI.fit(FPCA, mfd; basis=:pooled, ncomponents=1, maxK=2) + @test pooled.settings.max_k == 2 + @test pooled.settings.maxK == 2 +end diff --git a/test/pace.jl b/test/pace.jl index b05a635..fe7fc20 100644 --- a/test/pace.jl +++ b/test/pace.jl @@ -21,7 +21,7 @@ import StatsAPI fd = IrregularFunctionalData(times, values; domain=(0.0, 1.0), ids=["curve-$i" for i in eachindex(times)]) sigma2 = 0.13 - model = @inferred StatsAPI.fit(FPCA, fd; grid, ncomponents=1, maxK=2, + model = @inferred StatsAPI.fit(FPCA, fd; grid, ncomponents=1, max_k=2, mean=supplied_mean, covariance=supplied_covariance, noise_variance=sigma2) @test model.method == :pace @test model.design == :sparse @@ -77,7 +77,7 @@ import StatsAPI expected_scores * ones(1, length(target)) .+ transpose(mu.(target)) displayed = sprint(show, model) - for choice in ("pace", "sparse", "selection", "maxK", "mean_source", + for choice in ("pace", "sparse", "selection", "max_k", "mean_source", "covariance_source", "noise_variance_source", "kernel", "weighting", "binning", "eigen_tol", "sign_tol", "dropped", "noise_variance") @test occursin(choice, displayed) @@ -85,7 +85,7 @@ import StatsAPI scale = 1e-6 scaled_data = IrregularFunctionalData(times, [scale .* y for y in values]; domain=(0.0, 1.0)) - scaled_model = StatsAPI.fit(FPCA, scaled_data; grid, ncomponents=1, maxK=2, + scaled_model = StatsAPI.fit(FPCA, scaled_data; grid, ncomponents=1, max_k=2, mean=(; grid=supplied_mean.grid, values=scale .* supplied_mean.values), covariance=(; grid=source_grid, surface=scale^2 .* surface), noise_variance=scale^2 * sigma2) @@ -112,7 +112,7 @@ import StatsAPI for (rule, penalty) in ((AIC(), 2.0), (BIC(), log(8))) model = @inferred StatsAPI.fit(FPCA, data; method=:pace, grid, mean=(; grid, values=zeros(3)), covariance=(; grid, surface), - noise_variance=sigma2, maxK=3, ncomponents=rule) + noise_variance=sigma2, max_k=3, ncomponents=rule) criteria = map(1:3) do k covariance = functions[:, 1:k] * Diagonal(lambda[1:k]) * functions[:, 1:k]' + sigma2 * I factor = cholesky(Symmetric(covariance)) @@ -132,7 +132,7 @@ import StatsAPI for (rule, penalty) in ((AIC(), 2.0), (BIC(), log(8))) model = StatsAPI.fit(FPCA, data; method=:pace, grid, mean=(; grid, values=zeros(3)), covariance=(; grid, surface=rank_one_surface), - noise_variance=sigma2, maxK=3, ncomponents=rule) + noise_variance=sigma2, max_k=3, ncomponents=rule) @test length(model.all_eigenvalues) == 1 @test model.selection.candidates == 1:3 @test diff(model.selection.scores) ≈ fill(penalty, 2) @@ -147,7 +147,7 @@ import StatsAPI covariance = (; grid, surface) fd = IrregularFunctionalData([[0.0, 0.3, 1.0], [0.1, 0.8]], [[1.0, 1.4, 2.0], [1.2, 1.6]]) - options = (; grid, mean=mu, covariance, noise_variance=0.2, maxK=1) + options = (; grid, mean=mu, covariance, noise_variance=0.2, max_k=1) @test StatsAPI.fit(FPCA, fd; options...).method == :pace dense_times = [collect(range(0.0, 1.0; length=m)) for m in (22, 23)] dense_irregular = IrregularFunctionalData(dense_times, [1 .+ t for t in dense_times]) @@ -157,46 +157,46 @@ import StatsAPI @test_throws ArgumentError StatsAPI.fit(FPCA, fd; method=:invalid, options...) @test_throws ArgumentError StatsAPI.fit(FPCA, fd[1:1]; options...) @test StatsAPI.fit(FPCA, fd; mean=mu, covariance, - noise_variance=0.2, maxK=1).settings.grid_source == :observed_range_51 + noise_variance=0.2, max_k=1).settings.grid_source == :observed_range_51 @test StatsAPI.fit(FPCA, fd; grid, covariance, - noise_variance=0.2, maxK=1).settings.mean_bandwidth_source == :default_fraction + noise_variance=0.2, max_k=1).settings.mean_bandwidth_source == :default_fraction @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=mu, - noise_variance=0.2, maxK=1) + noise_variance=0.2, max_k=1) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance_bandwidth=1.0, noise_variance=0.2, maxK=1) + covariance_bandwidth=1.0, noise_variance=0.2, max_k=1) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance, maxK=1) + covariance, max_k=1) for noise in (-1.0, Inf, NaN, :invalid) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance, noise_variance=noise, maxK=1) + covariance, noise_variance=noise, max_k=1) end short_mean = (; grid=[0.2, 0.8], values=[1.2, 1.8]) @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=short_mean, - covariance, noise_variance=0.2, maxK=1) + covariance, noise_variance=0.2, max_k=1) asymmetric = copy(surface) asymmetric[1, 2] += 0.1 @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance=(; grid, surface=asymmetric), noise_variance=0.2, maxK=1) + covariance=(; grid, surface=asymmetric), noise_variance=0.2, max_k=1) rounded = copy(surface) rounded[1, 2] *= 1 + 4eps() near = StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance=(; grid, surface=rounded), noise_variance=0.2, maxK=1) + covariance=(; grid, surface=rounded), noise_variance=0.2, max_k=1) exact = StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance=(; grid, surface), noise_variance=0.2, maxK=1) + covariance=(; grid, surface), noise_variance=0.2, max_k=1) @test eigvals(near) ≈ eigvals(exact) rtol=1e-12 @test_throws DimensionMismatch StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance=(; grid, surface=ones(4, 4)), noise_variance=0.2, maxK=1) + covariance=(; grid, surface=ones(4, 4)), noise_variance=0.2, max_k=1) @test_throws DimensionMismatch StatsAPI.fit(FPCA, fd; grid, - mean=(; grid, values=zeros(4)), covariance, noise_variance=0.2, maxK=1) + mean=(; grid, values=zeros(4)), covariance, noise_variance=0.2, max_k=1) nonfinite = copy(surface) nonfinite[2, 2] = NaN @test_throws ArgumentError StatsAPI.fit(FPCA, fd; grid, mean=mu, - covariance=(; grid, surface=nonfinite), noise_variance=0.2, maxK=1) + covariance=(; grid, surface=nonfinite), noise_variance=0.2, max_k=1) unordered = UInt[0, 2, 1, 3] for supplied in ((; mean=(; grid=unordered, values=[1.0, 3.0, 2.0, 4.0]), covariance), (; mean=mu, covariance=(; grid=unordered, surface=ones(4, 4)))) error = try - StatsAPI.fit(FPCA, fd; grid, supplied..., noise_variance=0.2, maxK=1) + StatsAPI.fit(FPCA, fd; grid, supplied..., noise_variance=0.2, max_k=1) catch exception exception end @@ -205,7 +205,7 @@ import StatsAPI end @test_throws ArgumentError StatsAPI.fit(FPCA, fd; options..., presmooth=LocalPolynomial(; bandwidth=0.5)) singular = try - StatsAPI.fit(FPCA, fd; grid, mean=mu, covariance, noise_variance=0.0, maxK=1) + StatsAPI.fit(FPCA, fd; grid, mean=mu, covariance, noise_variance=0.0, max_k=1) catch error error end @@ -232,7 +232,7 @@ import StatsAPI fit32 = @inferred StatsAPI.fit(FPCA, data32; method=:pace, grid=grid32, mean=(; grid=grid32, values=1f0 .+ grid32), covariance=(; grid=grid32, surface=ones(Float32, 5, 5)), - noise_variance=0.2f0, maxK=1) + noise_variance=0.2f0, max_k=1) @test eltype(scores(fit32)) == Float32 @test eltype(score_covariances(fit32)) == Float32 @test eltype(eigvals(fit32)) == Float32 @@ -250,11 +250,11 @@ import StatsAPI @test covariance.noise_details.noise_interval == :full # A positive explicit noise variance separates integration of the # smoothers from invertibility of their estimated latent covariance. - estimated = @inferred StatsAPI.fit(FPCA, fd; grid, ncomponents=1, maxK=2, + estimated = @inferred StatsAPI.fit(FPCA, fd; grid, ncomponents=1, max_k=2, mean_bandwidth=1.2, covariance_bandwidth=1.2, kernel=KernelSmoothing.EpanechnikovKernel(), noise_variance=0.2, weighting=:subj, binning=5) - supplied = @inferred StatsAPI.fit(FPCA, fd; grid, ncomponents=1, maxK=2, + supplied = @inferred StatsAPI.fit(FPCA, fd; grid, ncomponents=1, max_k=2, mean, covariance, noise_variance=0.2) @test vec(estimated.mean.values) ≈ mean.values @test estimated.covariance ≈ covariance.surface @@ -267,10 +267,10 @@ import StatsAPI @test supplied.settings.mean_input.bandwidth == 1.2 @test supplied.settings.covariance_input.binning == 5 @test covariance.noise_variance > 0 - inherited = StatsAPI.fit(FPCA, fd; grid, ncomponents=1, maxK=2, mean, covariance) + inherited = StatsAPI.fit(FPCA, fd; grid, ncomponents=1, max_k=2, mean, covariance) @test noise_variance(inherited) == covariance.noise_variance @test inherited.settings.noise_variance_source == :estimated - estimated_noise = StatsAPI.fit(FPCA, fd; grid, ncomponents=1, maxK=2, + estimated_noise = StatsAPI.fit(FPCA, fd; grid, ncomponents=1, max_k=2, mean_bandwidth=1.2, covariance_bandwidth=1.2, kernel=KernelSmoothing.EpanechnikovKernel(), weighting=:subj, binning=5) @test noise_variance(estimated_noise) == covariance.noise_variance @@ -289,7 +289,7 @@ import StatsAPI sampling=SparseSampling((0.0, 1.0), 3:8)) truth_functions = [f(t) for t in grid, f in functions] surface = Matrix(Symmetric(truth_functions * Diagonal(lambda) * truth_functions')) - model = StatsAPI.fit(FPCA, simulation.data; grid, ncomponents=2, maxK=2, + model = StatsAPI.fit(FPCA, simulation.data; grid, ncomponents=2, max_k=2, mean=(; grid, values=mu.(grid)), covariance=(; grid, surface), noise_variance=sigma2) # The true affine functions are exactly interpolated. Convert their diff --git a/test/pace_defaults.jl b/test/pace_defaults.jl index 7c83a30..72a4319 100644 --- a/test/pace_defaults.jl +++ b/test/pace_defaults.jl @@ -17,7 +17,7 @@ import RecipesBase @test model.grid == collect(range(extrema(small.times)...; length=51)) @test domain(small) == (-1.0, 2.0) @test extrema(model.grid) != domain(small) - @test model.settings.maxK == 20 + @test model.settings.max_k == 20 @test model.selection.rule == :FVE @test model.selection.threshold == 0.99 @test model.settings.kernel isa KernelSmoothing.GaussianKernel @@ -50,7 +50,7 @@ import RecipesBase @test length(RecipesBase.apply_recipe(Dict{Symbol,Any}(), payload)) == 2 for setting in ("grid_source", "GaussianKernel", "weighting", "obs", "mean_bandwidth", "covariance_bandwidth", "diagonal_bandwidth", - "requested_binning", "binning", "maxK", "0.99", + "requested_binning", "binning", "max_k", "0.99", "noise_variance", "mean_source", "covariance_source", "eigen_tol", "sign_tol", "regularization") @test occursin(setting, sprint(show, model)) @@ -77,16 +77,16 @@ import RecipesBase @test dense.method == :dense @test dense.grid == grid @test dense.grid !== regular.grid - @test dense.selection.maxK == 30 # The dense default has no cap of 20. + @test dense.selection.max_k == 30 # The dense default has no cap of 20. @test dense.selection.threshold == 0.99 @test noise_variance(dense) === nothing end @testset "Multivariate routes inherit rotated/full noise" begin data = MultivariateFunctionalData([small, small]) - marginal = StatsAPI.fit(FPCA, data; ncomponents=2, maxK=3) + marginal = StatsAPI.fit(FPCA, data; ncomponents=2, max_k=3) joint = StatsAPI.fit(MFPCA, data; ncomponents=2, - univariate=(; ncomponents=2, maxK=3)) + univariate=(; ncomponents=2, max_k=3)) for fitted in (marginal, joint), j in 1:2 node = component(fitted, j) details = node.settings.covariance_input.noise_details @@ -99,12 +99,12 @@ import RecipesBase @testset "Overrides, selectors, and errors" begin grid = collect(range(extrema(small.times)...; length=11)) - explicit = StatsAPI.fit(FPCA, small; grid, ncomponents=1, maxK=3, + explicit = StatsAPI.fit(FPCA, small; grid, ncomponents=1, max_k=3, mean_bandwidth=0.6, covariance_bandwidth=0.8, kernel=KernelSmoothing.EpanechnikovKernel(), weighting=:subj, binning=7, noise_variance=0.25, eigen_tol=1e-9, sign_tol=1e-11) @test explicit.grid == grid - @test explicit.settings.maxK == 3 + @test explicit.settings.max_k == 3 @test ncomponents(explicit) == 1 @test explicit.settings.mean_input.bandwidth == 0.6 @test explicit.settings.covariance_input.bandwidth == 0.8 @@ -166,7 +166,7 @@ import RecipesBase (; grid=fitted_covariance.grid, surface=fitted_covariance.surface)) for mean in means, covariance in covariances options = (; grid, mean, covariance, noise_variance=0.25, - ncomponents=1, maxK=3) + ncomponents=1, max_k=3) reused = @inferred StatsAPI.fit(FPCA, small; options...) @test eigvals(reused) ≈ eigvals(explicit) @test scores(reused) ≈ scores(explicit) @@ -189,7 +189,7 @@ import RecipesBase # still control the mean when it must be estimated. for covariance in covariances options = (; grid, covariance, noise_variance=0.25, - ncomponents=1, maxK=3, mean_bandwidth=0.6, + ncomponents=1, max_k=3, mean_bandwidth=0.6, kernel=KernelSmoothing.EpanechnikovKernel(), weighting=:subj) estimated_mean = @inferred StatsAPI.fit(FPCA, small; options...) @test estimated_mean.settings.mean_source === :estimated @@ -203,7 +203,7 @@ import RecipesBase end # Original fitted settings remain available when an input is reused. reused = StatsAPI.fit(FPCA, small; grid, mean=fitted_mean, - covariance=fitted_covariance, noise_variance=0.25, ncomponents=1, maxK=3) + covariance=fitted_covariance, noise_variance=0.25, ncomponents=1, max_k=3) @test reused.settings.mean_input.kernel isa KernelSmoothing.EpanechnikovKernel @test reused.settings.mean_input.weighting === :subj @test reused.settings.covariance_input.binning == 7 @@ -219,7 +219,7 @@ import RecipesBase covariance = (; grid, surface=ones(3, 3)) for mean in (nothing, (; grid, values=1.5 .+ grid)) model = StatsAPI.fit(FPCA, fd; grid, mean, covariance, - noise_variance=0.25, ncomponents=1, maxK=1) + noise_variance=0.25, ncomponents=1, max_k=1) @test model.settings.requested_binning === :auto @test model.settings.binning === nothing @test model.settings.binning_rule === :not_applicable @@ -366,7 +366,7 @@ import RecipesBase @test eigenvalue_relative_error < 0.15 agreements = Float64[] for k in 1:2 - dense_values = basismatrix(dense.basis, pace.grid)[:, k] + dense_values = basis_matrix(dense.basis, pace.grid)[:, k] dense_values ./= sqrt(dot(pace.weights, dense_values.^2)) agreement = abs(dot(pace.weights, dense_values .* pace.basis.values[:, k])) push!(agreements, agreement) diff --git a/test/pace_defaults_golden.jl b/test/pace_defaults_golden.jl index ce8e91c..f51ff3b 100644 --- a/test/pace_defaults_golden.jl +++ b/test/pace_defaults_golden.jl @@ -21,7 +21,7 @@ using DelimitedFiles k = Int(only(reference(table, prefix * "selected_k"))) step = model.grid[2] - model.grid[1] system = FunctionalDataAnalysis._fpca_eigen(model.covariance, fill(step, m), - vec(model.mean.values); maxK=model.selection.maxK, + vec(model.mean.values); max_k=model.selection.max_k, normalization_weights=model.weights) @test relative_error(system.eigenvalues[1:k], reference(table, prefix * "eigenvalue")) < 1e-8 @@ -51,7 +51,7 @@ using DelimitedFiles grid = reference(table, "grid") model = StatsAPI.fit(FPCA, fd; grid, mean_bandwidth=0.30, covariance_bandwidth=0.40, kernel=KernelSmoothing.EpanechnikovKernel(), - weighting=:obs, binning=nothing, ncomponents=3, maxK=5, + weighting=:obs, binning=nothing, ncomponents=3, max_k=5, noise_variance=sigma2) @test grid ≈ sort(unique(reduce(vcat, first.(collect(fd))))) @test relative_error(vec(model.mean.values), reference(table, "mean")) < 1e-8 @@ -76,7 +76,7 @@ using DelimitedFiles m = length(model.grid) @test model.grid ≈ reference(table, "native_grid") rtol=1e-14 @test m == 51 - @test model.selection.maxK == Int(only(reference(table, "default_maxK"))) + @test model.selection.max_k == Int(only(reference(table, "default_maxK"))) @test model.settings.mean_settings.bandwidth == reference(table, "default_bandwidth")[1] @test model.settings.covariance_settings.bandwidth == reference(table, "default_bandwidth")[2] @test model.settings.kernel isa KernelSmoothing.GaussianKernel diff --git a/test/pace_golden.jl b/test/pace_golden.jl index e090904..ecaeb98 100644 --- a/test/pace_golden.jl +++ b/test/pace_golden.jl @@ -22,7 +22,7 @@ using DelimitedFiles model = StatsAPI.fit(FPCA, fd; grid, mean=supplied_mean, covariance=supplied_covariance, noise_variance=sigma2, - ncomponents=3, maxK=5) + ncomponents=3, max_k=5) @test model.method == :pace @test nobs(model) == 40 @test ids(model) == labels @@ -35,13 +35,13 @@ using DelimitedFiles # trapezoidal eigenstep agrees with fdapace's rectangular operator. @test relative_error(eigvals(model), reference("eigenvalue")) < 1e-8 @test relative_error(model.all_eigenvalues, reference("all_eigenvalue")) < 1e-8 - @test relative_error(basismatrix(eigenfunctions(model), grid), + @test relative_error(basis_matrix(eigenfunctions(model), grid), reshape(reference("eigenfunction"), length(grid), 3)) < 1e-8 @test relative_error(fve(model), reference("fve")) < 1e-8 @test relative_error(scores(model), reshape(reference("score"), n, 3)) < 1e-8 for i in 1:n rows = (golden[:, 1] .== "score_covariance") .& (golden[:, 2] .== i) - # fdapace SubsetFPCA leaves xiVar at maxK even when selecting K=3; + # fdapace SubsetFPCA leaves xiVar at max_k even when selecting K=3; # the reference contains its selected upper-left 3-by-3 block. expected = reshape(Float64.(golden[rows, 4]), 3, 3) @test relative_error(score_covariances(model)[:, :, i], expected) < 1e-8 @@ -49,8 +49,8 @@ using DelimitedFiles @test relative_error(model.settings.full_projected_covariance, reshape(reference("fitted_covariance"), length(grid), length(grid))) < 1e-8 # The score denominator includes components 4 and 5 even though K=3. - selected_surface = basismatrix(eigenfunctions(model), grid) * - Diagonal(eigvals(model)) * transpose(basismatrix(eigenfunctions(model), grid)) + selected_surface = basis_matrix(eigenfunctions(model), grid) * + Diagonal(eigvals(model)) * transpose(basis_matrix(eigenfunctions(model), grid)) @test norm(model.settings.full_projected_covariance - selected_surface) > 0.1 training_prediction = StatsAPI.predict(model, fd) @test training_prediction.scores == scores(model) @@ -68,7 +68,7 @@ using DelimitedFiles for (j, (selector, quantity)) in enumerate(((AIC(), "aic"), (BIC(), "bic"))) chosen = StatsAPI.fit(FPCA, fd; grid, mean=supplied_mean, covariance=supplied_covariance, noise_variance=sigma2, - ncomponents=selector, maxK=5) + ncomponents=selector, max_k=5) @test ncomponents(chosen) == selected[j + 1] @test chosen.selection.global_minimizer == global_minimizers[j] @test relative_error(chosen.selection.scores, reference(quantity)) < 1e-8 @@ -78,6 +78,6 @@ using DelimitedFiles end fve_model = StatsAPI.fit(FPCA, fd; grid, mean=supplied_mean, covariance=supplied_covariance, noise_variance=sigma2, - ncomponents=FVE(0.99), maxK=5) + ncomponents=FVE(0.99), max_k=5) @test ncomponents(fve_model) == selected[4] end diff --git a/test/plotting.jl b/test/plotting.jl index 2e49a80..20ba46d 100644 --- a/test/plotting.jl +++ b/test/plotting.jl @@ -11,12 +11,12 @@ function plotting_fixture() latent = [sqrt(5.0) 0; -sqrt(5.0) 0; 0 sqrt(1.25); 0 -sqrt(1.25); 0 0; 0 0] fd = GridFunctionalData(1 .+ grid' .+ latent * phi', grid) model = @inferred (data -> StatsAPI.fit(FPCA, data; - method=:dense, ncomponents=FVE(0.7), maxK=2))(fd) + method=:dense, ncomponents=FVE(0.7), max_k=2))(fd) times = [[0.0, 0.1, 0.7], [0.1, 0.4, 1.0], [0.0, 0.4, 0.7, 1.0], [0.4]] values = [1 .+ t .+ i/4 for (i, t) in enumerate(times)] irregular = IrregularFunctionalData(times, values; domain=(0.0, 1.0)) surface = phi * Diagonal([2.0, 0.5]) * phi' - pace = @inferred StatsAPI.fit(FPCA, irregular; grid, ncomponents=1, maxK=2, + pace = @inferred StatsAPI.fit(FPCA, irregular; grid, ncomponents=1, max_k=2, mean=(; grid, values=1 .+ grid), covariance=(; grid, surface), noise_variance=0.1) (; grid, phi, latent, fd, model, times, values, irregular, pace, surface) end @@ -121,7 +121,7 @@ end @test_throws ArgumentError modes_of_variation(f.model, 1; multipliers=bad) end multi = MultivariateFunctionalData([f.fd, GridFunctionalData(2 .* f.fd.values, f.grid .+ 2)]) - joint = @inferred StatsAPI.fit(MFPCA, multi; ncomponents=2, univariate=(; ncomponents=2, maxK=2)) + joint = @inferred StatsAPI.fit(MFPCA, multi; ncomponents=2, univariate=(; ncomponents=2, max_k=2)) jm = @inferred modes_of_variation(joint, 1) @test length(jm.nodes) == 2 @test (@inferred mean_data(joint)).nodes[2].grid == f.grid .+ 2 @@ -137,7 +137,7 @@ end component(joint, 1).all_eigenvalues @test modes_of_variation(joint.univariate, 1).nodes[2].eigenvalue ≈ 8 samegrid = MultivariateFunctionalData([f.fd, GridFunctionalData(2 .* f.fd.values, f.grid)]) - pooled = StatsAPI.fit(FPCA, samegrid; basis=:pooled, ncomponents=1, maxK=2) + pooled = StatsAPI.fit(FPCA, samegrid; basis=:pooled, ncomponents=1, max_k=2) @test scree_data(pooled).eigenvalues ≈ [5.0, 1.25] @test modes_of_variation(pooled, 1).nodes[2].eigenvalue ≈ 5 @test length(mean_data(pooled).nodes) == 2 @@ -239,19 +239,19 @@ end f.fd.values[1, 1] = 800 @test fitted_curve_data(f.model, 1).observed.values[1] != 800 original = plotting_fixture().fd - smooth = StatsAPI.fit(FPCA, original; ncomponents=1, maxK=2, + smooth = StatsAPI.fit(FPCA, original; ncomponents=1, max_k=2, presmooth=LocalPolynomial(; bandwidth=2.0, degree=0)) @test fitted_curve_data(smooth, 1).observed.values == original.values[1, :] @test norm(fitted_curve_data(smooth, 1).values - original.values[1, :]) > 0.01 masked_values = Matrix{Union{Missing,Float64}}(original.values) masked_values[1, 2] = missing masked = GridFunctionalData(masked_values, f.grid) - mf = StatsAPI.fit(FPCA, masked; grid=f.grid, ncomponents=1, maxK=2, + mf = StatsAPI.fit(FPCA, masked; grid=f.grid, ncomponents=1, max_k=2, mean=(; grid=f.grid, values=1 .+ f.grid), covariance=(; grid=f.grid, surface=f.surface), noise_variance=0.1) @test fitted_curve_data(mf, 1).observed.times == f.grid[[1, 3, 4, 5]] multi = MultivariateFunctionalData([original, GridFunctionalData(2 .* original.values, f.grid .+ 2)]) - joint = StatsAPI.fit(MFPCA, multi; ncomponents=2, univariate=(; ncomponents=2, maxK=2)) + joint = StatsAPI.fit(MFPCA, multi; ncomponents=2, univariate=(; ncomponents=2, max_k=2)) jd = @inferred fitted_curve_data(joint, 2; grid=[grid, grid .+ 2]) @test length(jd.nodes) == 2 && jd.curve == 2 @test jd.nodes[2].observed.times == f.grid .+ 2 @@ -308,7 +308,7 @@ end @test [s.plotattributes[:subplot] for s in panels] == [1, 2] @test all(s -> s.plotattributes[:layout] == (2, 1), panels) @test occursin("Node 2", panels[2].plotattributes[:title]) - joint = StatsAPI.fit(MFPCA, multi; ncomponents=2, univariate=(; ncomponents=2, maxK=2)) + joint = StatsAPI.fit(MFPCA, multi; ncomponents=2, univariate=(; ncomponents=2, max_k=2)) modes = plotting_series(FunctionalDataAnalysis._plot_payload(Val(:modes), modes_of_variation(joint, 1))) @test length(modes) == 10 @test occursin("Node 2: Mode", modes[6].plotattributes[:title]) diff --git a/test/reference/mfpca_m4.jl b/test/reference/mfpca_m4.jl index 3352a6b..cdf4135 100644 --- a/test/reference/mfpca_m4.jl +++ b/test/reference/mfpca_m4.jl @@ -20,7 +20,7 @@ nodes = [begin GridFunctionalData(j .+ g' .+ joint[:, 2j-1:2j] * phi', g) end for (j,g) in enumerate(grids)] data = MultivariateFunctionalData(nodes) -fitted = StatsAPI.fit(FPCA, data; basis=:marginal, ncomponents=2, maxK=2) +fitted = StatsAPI.fit(FPCA, data; basis=:marginal, ncomponents=2, max_k=2) open(joinpath(@__DIR__, "..", "data", "m4_mfpca_input.csv"), "w") do io println(io, "# Julia ", VERSION, "; FunctionalDataAnalysis 0.1.0; deterministic 24 subjects, 3 nodes; quantity,node,row,column,value") function write_matrix(quantity, node, a) @@ -45,7 +45,7 @@ times = [g[[1, 2+i%2, 5]] for i in 1:8] sparse = IrregularFunctionalData(times, [fill(5+cos(i), length(t)) for (i,t) in enumerate(times)]; domain=(2,3)) pace_data = MultivariateFunctionalData([dense, sparse]) -options = (; ncomponents=[2, 1], maxK=[2, 1], grid=[nothing, g], +options = (; ncomponents=[2, 1], max_k=[2, 1], grid=[nothing, g], mean=[nothing, (; grid=g, values=zeros(5))], covariance=[nothing, (; grid=g, surface=ones(5, 5))], noise_variance=[nothing, 0.2]) diff --git a/test/runtests.jl b/test/runtests.jl index 999d121..271a3f0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ using Test @testset "FunctionalDataAnalysis" begin include("aqua.jl") + include("naming.jl") include("quadrature.jl") include("bases.jl") include("data.jl") diff --git a/test/smoothed_fpca.jl b/test/smoothed_fpca.jl index bf23e72..76c84d0 100644 --- a/test/smoothed_fpca.jl +++ b/test/smoothed_fpca.jl @@ -23,7 +23,7 @@ end @test penalty(FourierBasis((0, 2), 1), :harmonic_acceleration) == zeros(1,1) @test (@inferred penalty(FourierBasis((0.0f0, 2.0f0), 5), :harmonic_acceleration)) isa Matrix{Float32} points = collect(range(0, 2; length=401)) - L = basismatrix(b, points; derivative=3) + pi^2*basismatrix(b, points; derivative=1) + L = basis_matrix(b, points; derivative=3) + pi^2*basis_matrix(b, points; derivative=1) @test transpose(L)*(quadrature_weights(points).*L) ≈ P atol=1e-9 rtol=1e-12 end @@ -48,10 +48,10 @@ end @test_throws ArgumentError fitted_bands(model) points = [0.0, .073, .2, .497, .799, 1.0] @test reconstruct(model; grid=points).values ≈ evaluate(fd, points) - @test basismatrix(eigenfunctions(model), points; derivative=1) ≈ - basismatrix(b, points; derivative=1)*eigenfunctions(model).coefs - @test basismatrix(DerivativeBasis(eigenfunctions(model)), points) ≈ - basismatrix(eigenfunctions(model), points; derivative=1) + @test basis_matrix(eigenfunctions(model), points; derivative=1) ≈ + basis_matrix(b, points; derivative=1)*eigenfunctions(model).coefs + @test basis_matrix(DerivativeBasis(eigenfunctions(model)), points) ≈ + basis_matrix(eigenfunctions(model), points; derivative=1) @test_throws ArgumentError derivative(model) @test penalty(eigenfunctions(model), 2) ≈ penalty(b, 2) @test domain(eigenfunctions(model)) == domain(fd) @@ -135,7 +135,7 @@ end f = @inferred StatsAPI.fit(FPCA, left; lambda, ncomponents=2) @test eigvals(f) ≈ vec(ref["eigenvalues"])[1:2] rtol=1e-8 @test f.all_eigenvalues ≈ vec(ref["eigenvalues"]) rtol=1e-8 - @test basismatrix(eigenfunctions(f), points) ≈ ref["harmonics"] rtol=1e-8 + @test basis_matrix(eigenfunctions(f), points) ≈ ref["harmonics"] rtol=1e-8 @test scores(f) ≈ ref["scores"] rtol=1e-8 @test scree_data(f).individual[1:2] ≈ vec(ref["varprop"]) rtol=1e-8 @test f.settings.original_harmonic_norms[1:2] ≈ vec(ref["harmonic_norms"]) rtol=1e-8 @@ -143,7 +143,7 @@ end @test vec(sum(abs2, scores(f); dims=1))/nobs(left) ≈ eigvals(f) ./ f.settings.original_harmonic_norms[1:2].^2 rtol=1e-8 @test StatsAPI.predict(f,left) ≈ scores(f) - @test reconstruct(f; grid=points).values ≈ evaluate(f.mean,points) .+ scores(f)*transpose(basismatrix(eigenfunctions(f),points)) + @test reconstruct(f; grid=points).values ≈ evaluate(f.mean,points) .+ scores(f)*transpose(basis_matrix(eigenfunctions(f),points)) cc = @inferred StatsAPI.fit(FunctionalCCA, left, right; lambda) @test cc.correlations ≈ vec(ref["correlations"]) rtol=1e-8 @test cc.squared_correlations ≈ vec(ref["correlations"]).^2 rtol=1e-8 @@ -177,19 +177,19 @@ end values = 1 .+ transpose(grid) .+ latent*transpose(functions) fd = GridFunctionalData(values,grid) dense = @inferred (data -> StatsAPI.fit(FPCA,data; - method=:dense,ncomponents=2,maxK=2))(fd) + method=:dense,ncomponents=2,max_k=2))(fd) @test (@inferred StatsAPI.predict(dense,fd)) ≈ scores(dense) @test (@inferred reconstruct(dense)).values ≈ values @test fitted_curve_data(dense,1).observed.values == values[1,:] @test eigenfunctions(dense) isa FPCBasis - face = @inferred (data -> StatsAPI.fit(FPCA,data;method=:face,ncomponents=2,maxK=2, + face = @inferred (data -> StatsAPI.fit(FPCA,data;method=:face,ncomponents=2,max_k=2, mean=(;grid,values=1 .+ grid),face=(;knots=4,lambda=1.0)))(fd) @test StatsAPI.predict(face,fd).scores ≈ scores(face) @test eigenfunctions(face) isa FPCBasis @test size(reconstruct(face).values) == size(values) irregular = to_irregular(fd) surface = functions*Diagonal([2.,.5])*transpose(functions) - pace = @inferred StatsAPI.fit(FPCA,irregular;grid,ncomponents=2,maxK=2, + pace = @inferred StatsAPI.fit(FPCA,irregular;grid,ncomponents=2,max_k=2, mean=(;grid,values=1 .+ grid), covariance=(;grid,surface),noise_variance=.1) @test StatsAPI.predict(pace,irregular).scores ≈ scores(pace) @test eigenfunctions(pace) isa FPCBasis @@ -197,10 +197,10 @@ end # Basis FPCA must not widen the established multivariate fit APIs. multivariate = MultivariateFunctionalData([fd, fd]) marginal = @inferred (data -> StatsAPI.fit(FPCA,data; - ncomponents=2,maxK=2))(multivariate) + ncomponents=2,max_k=2))(multivariate) @test scores(marginal) == [scores(dense), scores(dense)] joint = @inferred StatsAPI.fit(MFPCA,multivariate;ncomponents=2, - univariate=(;ncomponents=2,maxK=2)) + univariate=(;ncomponents=2,max_k=2)) @test (@inferred StatsAPI.predict(joint,multivariate)) ≈ scores(joint) @test length((@inferred modes_of_variation(joint,1)).nodes) == 2 growth = load_dataset(:growth;sample=true).data @@ -246,8 +246,8 @@ end @test_throws DimensionMismatch StatsAPI.fit(FPCA,fd;penalty=ones(3,3)) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;penalty=-1) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;center=1) - @test_throws ArgumentError StatsAPI.fit(FPCA,fd;maxK=0) - @test_throws ArgumentError StatsAPI.fit(FPCA,fd;maxK=5) + @test_throws ArgumentError StatsAPI.fit(FPCA,fd;max_k=0) + @test_throws ArgumentError StatsAPI.fit(FPCA,fd;max_k=5) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;ncomponents=0) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;ncomponents=AIC()) @test_throws ArgumentError StatsAPI.fit(FPCA,fd;eigen_tol=NaN) From 7bb02d6162197fb796fb9d5bb7e0e286bb71c9c5 Mon Sep 17 00:00:00 2001 From: Mladen Kolar Date: Sun, 25 Oct 2026 16:27:24 +0400 Subject: [PATCH 2/3] Read Julia naming declarations explicitly as UTF-8 --- scripts/check_naming.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/check_naming.py b/scripts/check_naming.py index e112dd8..e4cc472 100644 --- a/scripts/check_naming.py +++ b/scripts/check_naming.py @@ -17,12 +17,12 @@ def check(root): errors=[] types={m.group(1) for path in (root/'src').rglob('*.jl') - for line in path.read_text().splitlines() if (m := TYPE.match(line))} - types |= {m.group(1) for path in (root/'src').rglob('*.jl') for m in re.finditer(r'^\s*const\s+(\w+)\s*=',path.read_text(),re.M)} + for line in path.read_text(encoding="utf-8").splitlines() if (m := TYPE.match(line))} + types |= {m.group(1) for path in (root/'src').rglob('*.jl') for m in re.finditer(r'^\s*const\s+(\w+)\s*=',path.read_text(encoding="utf-8"),re.M)} for path in sorted((root/'src').rglob('*.jl')): if 'deprecated' in path.stem: continue in_doc=False; in_struct=False - lines=path.read_text().splitlines() + lines=path.read_text(encoding="utf-8").splitlines() for number,line in enumerate(lines,1): if line.count('"""') % 2: in_doc=not in_doc; continue if in_doc: continue From cf99fe01db465406dd5c9386d2a4c2e728c351d3 Mon Sep 17 00:00:00 2001 From: Mladen Kolar Date: Sun, 25 Oct 2026 16:53:46 +0400 Subject: [PATCH 3/3] Exercise renamed plotting conveniences at the backend boundary --- test/plotting.jl | 51 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/test/plotting.jl b/test/plotting.jl index 20ba46d..816b01b 100644 --- a/test/plotting.jl +++ b/test/plotting.jl @@ -363,3 +363,54 @@ end @test !isdefined(Main, :Plots) end end + +# Convenience functions must prepare the correct recipe payload and forward +# plot attributes even when no plotting backend is loaded. Inspect the actual +# missing-backend call instead of installing test-only RecipesBase methods. +function plotting_boundary_payload(call) + failure = try + call() + nothing + catch error + error + end + @test failure isa MethodError + failure isa MethodError || return nothing + if failure.f === Core.kwcall + @test failure.args[1] == (; color=:blue) + @test failure.args[2] === RecipesBase.plot + else + @test failure.f === RecipesBase.plot + end + last(failure.args) +end + +@testset "plot convenience functions reach the backend boundary" begin + f = plotting_fixture() + cross = StatsAPI.fit(CrossCovariance, f.fd, f.fd) + bands = fitted_bands(f.pace; curves=[1]) + mean = StatsAPI.fit(SmoothedMean, f.fd; grid=f.grid, bandwidth=2.0) + mean_uncertainty = mean_bands(f.fd, mean; method=:sandwich) + custom_grid = [0.0, 0.5, 1.0] + cases = ( + (() -> modes_plot(f.model, 1; multipliers=[-1, 1], color=:blue), + :modes, modes_of_variation(f.model, 1; multipliers=[-1, 1])), + (() -> scree_plot(f.model; color=:blue), :scree, scree_data(f.model)), + (() -> cov_plot(cross; color=:blue), :covariance, covariance_data(cross)), + (() -> cov_plot(f.model; surface=:fitted, color=:blue), + :covariance, covariance_data(f.model; surface=:fitted)), + (() -> design_plot(f.irregular; diagonal=true, color=:blue), + :design, design_data(f.irregular; diagonal=true)), + (() -> fit_plot(f.model, 1; color=:blue), :fit, fitted_curve_data(f.model, 1)), + (() -> fit_plot(f.model, 1; grid=custom_grid, color=:blue), + :fit, fitted_curve_data(f.model, 1; grid=custom_grid)), + (() -> band_plot(bands; color=:blue), :band, band_data(bands)), + (() -> band_plot(mean_uncertainty; coverage=:interval, color=:blue), + :band, band_data(mean_uncertainty; coverage=:interval)), + ) + for (call, kind, expected) in cases + payload = plotting_boundary_payload(call) + @test payload isa FunctionalDataAnalysis._FDAPlot{kind} + @test isequal(payload.data, expected) + end +end