optimum.reparam(method = "RBFGS") is now
as accurate as "DP": on random function pairs its median amplitude
distance is 0.94-0.99x "DP"'s and at most 1.05x (it stopped in local
minima, up to 3x), at about "DP"'s run time for 101 points and 3x and 9x
less for 201 and 401 points. It starts from a "DP" warp on a coarse grid
of at most 101 points and refines it on the SRVFs resampled scl times
finer (scl now applies to "RBFGS" as well as "DPo"); with
lambda > 0 it also starts from the identity and keeps the lower cost, and
with lambda = 0 it keeps the "DP" start when the refined warp fits worse
on the original grid. MATLAB's coarse grid of 51 points and keeping the
refined warp left the worst pairs at up to 1.2x "DP"'s distance. In the RBFGS solver the SRVFs
are scaled to unit L2 rather than Euclidean norm, so the gradient-norm
stopping test no longer loosens with the grid size (lam is rescaled so it
weighs the penalty as before), and q2 is warped from the original by the
accumulated warping each step, so interpolation error no longer accumulates.
The internal rlbfgs() takes an optional starting warp gam0 and returns
a list with the warp gam and info = (iterations, stop reason, cost,
gradient norm)q2 is a natural cubic spline and the gradient
differentiates the discretised cost exactly (the old continuum gradient was
not a descent direction at warps with kinks, so the solver stalled). It
stops when the relative change of the cost stays below 1e-3 for more than
5 iterations after at least 10 (new stop reason 4 in rlbfgs()'s info).
optimum.reparam(method = "RBFGS") now refines the coarse "DP" warp on
SRVFs resampled linearly rather than by a spline, starting, as MATLAB and
the paper do, from a "DP" warp on a coarse grid of at most 51 points
(was 101). On "DP"'s own objective its median cost is 0.98x "DP"'s on
random function pairs and 0.86x on consecutive growth_vel pairs, at about
0.5x "DP"'s run time for 101 points and 0.16x for 201. From the coarse
start it can land in a different minimum: up to 1.4x "DP"'s cost on
random three-bump pairs (a 101 point start kept that to 1.03x at 2.3x the
run time)classify() for ampbox and phbox
objects aligns new functions to the fitted template (new align_new() for
fdawarp objects) and reports their distance to the median, boxplot region
(50%, 1 - alpha, within whiskers, outlier) and quartile side. The boxplot
statistics now also store alpha, the whisker scale (ka/kp), the region
edges r50, ralpha and whisker_dis, and the alignment templatef_to_srvf() gains a smooth argument (default FALSE; TRUE uses a
smoothing spline) and optimum.reparam() gains a scl argument for
method = "DPo" (see below)optimum.reparam(method = "DPo")
(DPQ, also used for curves) resamples the SRVFs three times finer for its
cost (new scl argument, default 3; scl = 1 uses the original grid and
matches MATLAB DP1 and Python DP to machine precision) and interpolates
q2 linearly instead of rounding to the nearest grid point, so warps are
less jagged and fit better. The previous hard-coded 5x upsampling is goneoptimum.reparam(method = "DP") (DPQ2) now integrates the SRVFs as linear
between samples (previously piecewise constant, last sample ignored) and
nbhd_dim defaults to 10 (was 7), here and in the curve alignment code;
it now fits as well as "DPo"lambda >= 0); results
are unchanged and the DP runs about 1.5-2.9x fastertime_warping() now finishes like MATLAB's fdawarp.time_warping: the
final matching step's warpings are centered and the original data are warped
once by them, so fn is exactly f composed with warping_functions (and
mqn is the mean of qn for the mean). Previously fn/qn kept the
previous iteration's alignment and warping_functions came from an extra
alignment pass, so they disagreed by a few percent. qun now holds one
relative template change per iteration. Results match MATLAB and Python
fdasrsf on simu_data to about 1e-8SqrtMeanInverse now starts its Karcher-mean iteration from the psi
closest to the pointwise mean, which lies on the unit sphere, rather than
from the pointwise mean itself, matching the Python and MATLAB
implementations (results change only at the level of the 1e-8 stopping
tolerance)simu_warp and simu_warp_median datasets are regenerated with the
updated time_warping() and DP solversf_to_srvf now differentiates the interpolating cubic spline and
srvf_to_f integrates q|q| with a spline-corrected trapezoid rule instead
of plain cumtrapz, so the round trip f -> q -> f is now O(h^4)
accurate (previously O(h^2) and poor for oscillatory functions)gradient.spline with smooth_data = FALSE) and reconstructs functions
from SRVFs (cumtrapzmid, mean/median function, boxplot quantiles) with the
same cubic-spline derivative/integral, so the two directions are consistentgam_to_psi (and so
gam_to_v) rescales psi to unit norm so it stays on the Hilbert sphere;
inv_exp_map and elastic.distance compute the angle with atan2 instead
of acos, which was only accurate to ~1e-8 near zero. gam_to_psi,
gam_to_v and gam_to_h now default to smooth = FALSE, and gam_to_h
clamps gamdot at eps before the logQ1a/Q3a when the
quantile points away from its quartile, as the phase boxplot already didkmeans_align(rotation = TRUE) now applies the estimated rotation to the
aligned curves instead of discarding itelastic.mlpcr.regression now uses the standard softmax in its
multinomial logistic loss and gradient, so fitting, warping and prediction
agree, and predict for these fits picks the class with the highest
probability rather than the lowestfind_basis_normal), return NaN distances when the inner product falls
slightly below -1 (inverse_exp, inverse_exp_coord), or fail on
(nearly) antipodal curves when parallel transportinggroup_action_by_gamma uses the correct grid spacing, and square roots of
warping-function derivatives (e.g. in SqrtMean and the phase boxplot) are
clamped at zero so small negative numerical derivatives no longer give
NaNconfigure script probes the toolchain and builds
without OpenMP if no usable runtime is found (src/Makevars is now
generated from src/Makevars.in)interparc function for down-sampling closed curvesoptimum.reparam now applies the selected penalty ("none", "roughness",
"l2gam", "l2psi", "geodesic") for the "DP" method as welltime_warping and ppd now accept every optimum.reparam penalty
("roughness", "l2gam", "l2psi", "geodesic", "none");
penalty_method = "norm", which always failed with "invalid penalty
selection", is now an alias for "l2gam"pair_align_functions, multiple_align_functions, elastic.distance and
elastic.depth now also accept pen = "norm" as an alias for "l2gam"
instead of failing with "invalid penalty selection", and document every
optimum.reparam penaltytime_warping and ppd no longer list optim_method = "DP2": its
coordinate-descent solver was removed from optimum.reparam in 2.1.2, so
the option always failed partway through with "'arg' should be one of ...".
It is now rejected up front; use "DP", "DPo" or "RBFGS". The
pair_align_functions and multiple_align_functions docs no longer
mention "DP2" eitherpredict for elastic.pcr.regression, elastic.lpcr.regression and
elastic.mlpcr.regression fits no longer fails when newdata is supplied;
new functions are now aligned with the lambda, penalty and optimization
method used by time_warping, which the vertFPCA, horizFPCA and
jointFPCA predict methods now also honorpredict for jointFPCAh fits no longer fails with "requires
numeric/complex matrix/vector arguments": it now projects onto the
horizontal basis stored in U1 and computes h without smoothing, as
jointFPCAh doesget_distance_matrix(scale = TRUE) now compares the scaled curves; it
ignored the scaling and returned all-zero distancesmultivariate_karcher_mean no longer fails with "object 'v' not found" when
it converges on its first iteration, and the returned qn, gamma, R and
v now all come from the same alignment steppredict for multivariate_pca fits no longer fails; curves are aligned
and mapped to the tangent space the same way multivariate_karcher_mean
does for the fitted curvesv_to_curve accepts a matrix of shooting vectors, one per columnppd chooses lambda from the persistent peaks again: the peak clustering
always came back empty and the exact-match test never succeededbootTB computes the phase tolerance bounds from the bootstrap samples
rather than from the original warping functionsgauss_model(sort_samples = TRUE) no longer fails with "the condition has
length > 1", and gauss_model stores gams with one warping function per
column, like warping_functionskmeans_align no longer fails when a cluster holds a single curve, and no
longer swaps the curve and SRVF of multivariate medoid templateselastic_ph_change_ff uses d Monte Carlo draws as documented,
elastic_amp_change_ff and elastic_ph_change_ff center the curve at the
changepoint with the "before" mean, and elastic_change_fpca no longer fails
when a single principal component is retainedjointFPCA(srvf = FALSE) chooses C by comparing its reconstructions with
the original functions rather than with their SRSFsreparam_curve(method = "DPo") returns the optimal rotation instead of
NULL, and the "DP" method honours modepair_align_image reads the width of the second image correctly and takes
image gradients of non-square images without failing, using each axis's own
spacingcurve_boxplot no longer fails with "argument is of length zero" on
multivariate_karcher_mean results, and with scale = FALSE it measures
distances to the lower whisker from the lower whisker"s" basis uses the right frequency and axis, composing warps
without the interp package no longer returns zeros, and
reparam_image records the previous energy rather than an unset one when
it rejects a step"SIMUL" method of optimum.reparam interpolates correctly across
several flat segments, and "DPo" falls back to "DP" whenever the two
grids differ anywheretime_warping now stops once the template changes by less than 1%, as the
Python and MATLAB implementations do, instead of as soon as the change
shrinks; qun now includes the last value computedalign_fPCA returns warping functions that match the aligned functionsfunction_group_warp_bayes now iterates the Karcher mean of the warping
functions instead of taking a single stepsample_shapes follows the whole geodesic for closed curves, supports
curves in more than two dimensions when rotation = FALSE, and keeps the
sign of sampled rotation angleswarp_q_gamma is defined once, fixing the regression
helpers that passed its arguments in the other order; duplicate definitions
of gam_to_psi, psi_to_gam, l2_norm, inner_product and Enorm were
removed; several unexported helpers (calc_j, calculate_variance,
inverse_exp, curve_align_sub, karcher_calc, elastic.regression,
elastic.logistic) were corrected; a memory leak in the multinomial
logistic warping code was fixedpair_align_image and reparam_image work again, for square and
non-square images. With interp 1.1.6, whose bicubic() is an
irregular-grid workaround, warping failed with "length of y0 and x0
differs!" or returned NAs, and every warp other than the identity paired
the wrong row and column coordinates. Images and warps are now
interpolated with the package's own tensor-product cubic spline (as in the
C++ q-map code): the identity warp reproduces the image exactly, and warps
are no longer upsampled 8 times, so iterations take a fraction of a second
instead of minutes. interp and fields are no longer used.
reparam_image now also checks the updated warp, rather than the previous
one, for foldingfasatPNSs2eppdhorizFPNS functioncurve_karcher_mean and curve_srvf_aligncurve_karcher_cov to multivariate_karcher_covcurve_principal_direcions to multivariate_pcamultivariate_karcher_meanmultivariate_pcamultivariate_karcher_mean (#59)curve_karcher_meancurve_karcher_mean and curve_srvf_alignsample_shapes and created curveboxplot (#39)elastic_amp_change_ff, elastic_ph_change_ff, and elastic_change_fpcacurve_depth function (#31)calc_shape_dist for different pre-shape spacesSqrtMeanInverse and inv_exp_map to global (#29)Added functionalities to kmeans_align() (#24). Specifically:
use_verbose which allows the user to suppress
information displayed in the console;Further small improvements, documentation updates and improved code coverage (#23).
bootTB() (#19);gfortran issuesgropt library