_images/logo.png

Functional Alignment#

Group-wise function alignment using SRSF framework and Dynamic Programming

moduleauthor:: J. Derek Tucker <jdtuck@sandia.gov>

time_warping.align_fPCA(f, time, num_comp=3, showplot=True, smoothdata=False, cores=-1)[source]#

aligns a collection of functions while extracting principal components. The functions are aligned to the principal components

Parameters:

  • f (np.ndarray) – numpy ndarray of shape (M,N) of N functions with M samples

  • time (np.ndarray) – vector of size M describing the sample points

  • num_comp – number of fPCA components

  • showplot – Shows plots of results using matplotlib (default = T)

  • smooth_data (bool) – Smooth the data using a box filter (default = F)

  • cores – number of cores for parallel (default = -1 (all))

Return type:

tuple of numpy array

Return fn:

aligned functions - numpy ndarray of shape (M,N) of N functions with M samples

Return qn:

aligned srvfs - similar structure to fn

Return q0:

original srvf - similar structure to fn

Return mqn:

srvf mean or median - vector of length M

Return gam:

warping functions - similar structure to fn

Return q_pca:

srsf principal directions

Return f_pca:

functional principal directions

Return latent:

latent values

Return coef:

coefficients

Return U:

eigenvectors

Return orig_var:

Original Variance of Functions

Return amp_var:

Amplitude Variance

Return phase_var:

Phase Variance

time_warping.align_fPLS(f, g, time, comps=3, showplot=True, smoothdata=False, delta=0.01, max_itr=100)[source]#

This function aligns a collection of functions while performing principal least squares

Parameters:

  • f (np.ndarray) – numpy ndarray of shape (M,N) of N functions with M samples

  • g (np.ndarray) – numpy ndarray of shape (M,N) of N functions with M samples

  • time (np.ndarray) – vector of size M describing the sample points

  • comps – number of fPLS components

  • showplot – Shows plots of results using matplotlib (default = T)

  • smooth_data (bool) – Smooth the data using a box filter (default = F)

  • delta – gradient step size

  • max_itr – maximum number of iterations

Return type:

tuple of numpy array

Return fn:

aligned functions - numpy ndarray of shape (M,N) of N functions with M samples

Return gn:

aligned functions - numpy ndarray of shape (M,N) of N functions with M samples

Return qfn:

aligned srvfs - similar structure to fn

Return qgn:

aligned srvfs - similar structure to fn

Return qf0:

original srvf - similar structure to fn

Return qg0:

original srvf - similar structure to fn

Return gam:

warping functions - similar structure to fn

Return wqf:

srsf principal weight functions

Return wqg:

srsf principal weight functions

Return wf:

srsf principal weight functions

Return wg:

srsf principal weight functions

Return cost:

cost function value

class time_warping.fdawarp(f, time)[source]#

This class provides alignment methods for functional data using the SRVF framework

Usage: obj = fdawarp(f,t)

Parameters:

  • f – (M,N): matrix defining N functions of M samples

  • time – time vector of length M

  • fn – aligned functions

  • qn – aligned srvfs

  • q0 – initial srvfs

  • fmean – Karcher mean

  • mqn – mean srvf

  • gam – warping functions

  • psi – srvf of warping functions

  • stats – alignment statistics

  • qun – cost function

  • lambda – lambda

  • method – optimization method

  • gamI – inverse warping function

  • rsamps – random samples

  • fs – random aligned functions

  • gams – random warping functions

  • ft – random warped functions

  • qs – random aligned srvfs

  • type – alignment type

  • mcmc – mcmc output if bayesian

Author : J. D. Tucker (JDT) <jdtuck AT sandia.gov> Date : 15-Mar-2018

gauss_model(n=1, sort_samples=False)[source]#

This function models the functional data using a Gaussian model extracted from the principal components of the srvfs

Parameters:

  • n (integer) – number of random samples

  • sort_samples (bool) – sort samples (default = T)

joint_gauss_model(n=1, no=3)[source]#

This function models the functional data using a joint Gaussian model extracted from the principal components of the srsfs

Parameters:

  • n (integer) – number of random samples

  • no (integer) – number of principal components (default = 3)

multiple_align_functions(mu, omethod='DP2', smoothdata=False, parallel=False, lam=0.0, pen='l2psi', cores=-1, grid_dim=7)[source]#

This function aligns a collection of functions using the elastic square-root slope (srsf) framework.

Usage: obj.multiple_align_functions(mu)

Parameters:

  • mu – vector of function to align to

  • omethod – optimization method (DP, DP2, RBFGS, cRBFGS) (default = DP2)

  • smoothdata (bool) – Smooth the data using a box filter (default = F)

  • parallel – run in parallel (default = F)

  • lam (double) – controls the elasticity (default = 0)

  • penalty – penalty type (default=”l2psi”) options are “roughness”, “l2gam”, “l2psi”, “geodesic” and “none”. The penalty is weighted by lam, so it has no effect when lam is 0. “roughness” and “geodesic” have no dynamic programming counterpart, so with a nonzero lam they need omethod “RBFGS” or “cRBFGS”

  • cores – number of cores for parallel (default = -1 (all))

  • grid_dim – size of the grid, for the DP2 method only (default = 7)

plot()[source]#

plot functional alignment results

Usage: obj.plot()

ppd(max_lam=2, num_lam=10, pt=0.15, srvf=True, mu=None, method='mean', omethod='DP2', center=True, smoothdata=False, MaxItr=20, parallel=True, pen='l2psi', cores=-1, grid_dim=7, verbose=True)[source]#

This computes the peak persistance diagram over a range of lambda. This can help determine the proper elasticity (penalty). This can be slow and recommended to run in parallel.

Parameters:

  • max_lam – max regularization parameter (default 2)

  • num_lam – number of steps (default 10)

  • pt – the percentile of negative curvature of raw data (default .15)

  • srvf – perform in SRVF space (default=True)

  • mu – vector of function to align to (default None)

  • method – (string) warp calculate Karcher Mean or Median (options = “mean” or “median”) (default=”mean”)

  • omethod – optimization method (DP, DP2, RBFGS, cRBFGS) (default = DP2)

  • center – center warping functions (default = T)

  • smoothdata (bool) – Smooth the data using a box filter (default = F)

  • MaxItr – Maximum number of iterations (default = 20)

  • parallel – run in parallel (default = F)

  • penalty – penalty type (default=”l2psi”) options are “roughness”, “l2gam”, “l2psi”, “geodesic” and “none”. The penalty is weighted by lam, so it has no effect when lam is 0. “roughness” and “geodesic” have no dynamic programming counterpart, so with a nonzero lam they need omethod “RBFGS” or “cRBFGS”

  • cores – number of cores for parallel (default = -1 (all))

  • grid_dim – size of the grid, for the DP2 method only (default = 7)

  • verbose – print status output (default = T)

  • thresh – cost function threshold (default = .01)

srsf_align(method='mean', omethod='DP2', center=True, smoothdata=False, MaxItr=20, parallel=False, lam=0.0, pen='l2psi', cores=-1, grid_dim=7, verbose=True, thresh=0.01)[source]#

This function aligns a collection of functions using the elastic square-root slope (srsf) framework.

Parameters:

  • method – (string) warp calculate Karcher Mean or Median (options = “mean” or “median”) (default=”mean”)

  • omethod – optimization method (DP, DP2, RBFGS, cRBFGS) (default = DP2)

  • center – center warping functions (default = T)

  • smoothdata (bool) – Smooth the data using a box filter (default = F)

  • MaxItr – Maximum number of iterations (default = 20)

  • parallel – run in parallel (default = F)

  • lam (double) – controls the elasticity (default = 0)

  • penalty – penalty type (default=”l2psi”) options are “roughness”, “l2gam”, “l2psi”, “geodesic” and “none”. The penalty is weighted by lam, so it has no effect when lam is 0. “roughness” and “geodesic” have no dynamic programming counterpart, so with a nonzero lam they need omethod “RBFGS” or “cRBFGS”

  • cores – number of cores for parallel (default = -1 (all))

  • grid_dim – size of the grid, for the DP2 method only (default = 7)

  • verbose – print status output (default = T)

  • thresh – cost function threshold (default = .01)

Examples >>> import tables >>> fun=tables.open_file(“../Data/simu_data.h5”) >>> f = fun.root.f[:] >>> f = f.transpose() >>> time = fun.root.time[:] >>> obj = fs.fdawarp(f,time) >>> obj.srsf_align()

time_warping.pairwise_align_bayes(f1i, f2i, time, mcmcopts=None)[source]#

This function aligns two functions using Bayesian framework. It will align f2 to f1. It is based on mapping warping functions to a hypersphere, and a subsequent exponential mapping to a tangent space. In the tangent space, the Z-mixture pCN algorithm is used to explore both local and global structure in the posterior distribution.

The Z-mixture pCN algorithm uses a mixture distribution for the proposal distribution, controlled by input parameter zpcn. The zpcn$betas must be between 0 and 1, and are the coefficients of the mixture components, with larger coefficients corresponding to larger shifts in parameter space. The zpcn[“probs”] give the probability of each shift size.

Usage: out = pairwise_align_bayes(f1i, f2i, time)

Parameters:

  • f1i – vector defining M samples of function 1

  • f2i – vector defining M samples of function 2

  • time – time vector of length M

  • mcmopts – dict of mcmc parameters

default mcmc options: tmp = {“betas”:np.array([0.5,0.5,0.005,0.0001]),”probs”:np.array([0.1,0.1,0.7,0.1])} mcmcopts = {“iter”:2*(10**4) ,”burnin”:np.minimum(5*(10**3),2*(10**4)//2), “alpha0”:0.1, “beta0”:0.1,”zpcn”:tmp,”propvar”:1, “initcoef”:np.repeat(0,20), “npoints”:200, “extrainfo”:True}

:rtype collection containing :return f2_warped: aligned f2 :return gamma: warping function :return g_coef: final g_coef :return psi: final psi :return sigma1: final sigma

if extrainfo :return accept: accept of psi samples :return betas_ind :return logl: log likelihood :return gamma_mat: posterior gammas :return gamma_stats: posterior gamma stats :return xdist: phase distance posterior :return ydist: amplitude distance posterior)

time_warping.pairwise_align_bayes_infHMC(y1i, y2i, time, mcmcopts=None)[source]#

This function aligns two functions using Bayesian framework. It uses a hierarchical Bayesian framework assuming mearsurement error error It will align f2 to f1. It is based on mapping warping functions to a hypersphere, and a subsequent exponential mapping to a tangent space. In the tangent space, the infty-HMC algorithm is used to explore both local and global structure in the posterior distribution.

Usage: out = pairwise_align_bayes_infHMC(f1i, f2i, time)

Parameters:

  • y1i – vector defining M samples of function 1

  • y2i – vector defining M samples of function 2

  • time – time vector of length M

  • mcmopts – dict of mcmc parameters

default mcmc options: mcmcopts = {“iter”:1*(10**4), “nchains”:4, “vpriorvar”:1, “burnin”:np.minimum(5*(10**3),2*(10**4)//2), “alpha0”:0.1, “beta0”:0.1, “alpha”:1, “beta”:1, “h”:0.01, “L”:4, “f1propvar”:0.0001, “f2propvar”:0.0001, “L1propvar”:0.3, “L2propvar”:0.3, “npoints”:200, “thin”:1, “sampfreq”:1, “initcoef”:np.repeat(0,20), “nbasis”:10, “basis”:’fourier’, “extrainfo”:True}

Basis can be ‘fourier’ or ‘legendre’

:rtype collection containing :return f2_warped: aligned f2 :return gamma: warping function :return v_coef: final v_coef :return psi: final psi :return sigma1: final sigma

if extrainfo :return theta_accept: accept of psi samples :return f2_accept: accept of f2 samples :return SSE: SSE :return gamma_mat: posterior gammas :return gamma_stats: posterior gamma stats :return xdist: phase distance posterior :return ydist: amplitude distance posterior)

    1. Tucker, L. Shand, and K. Chowdhary. “Multimodal Bayesian Registration of Noisy Functions using Hamiltonian Monte Carlo”, Computational Statistics and Data Analysis, accepted, 2021.

time_warping.pairwise_align_functions(f1, f2, time, omethod='DP2', lam=0, pen='l2psi', grid_dim=7)[source]#
This function aligns f2 to f1 using the elastic square-root

slope (srsf) framework.

Usage: out = pairwise_align_functions(f1, f2, time)
out = pairwise_align_functions(f1, f2, time, omethod, lam,

grid_dim)

Parameters:

  • f1 – vector defining M samples of function 1

  • f2 – vector defining M samples of function 2

  • time – time vector of length M

  • omethod – optimization method (DP, DP2, RBFGS, cRBFGS) (default = DP)

  • lam – controls the elasticity (default = 0)

  • penalty – penalty type (default=”l2psi”) options are “roughness”, “l2gam”, “l2psi”, “geodesic” and “none”. The penalty is weighted by lam, so it has no effect when lam is 0. “roughness” and “geodesic” have no dynamic programming counterpart, so with a nonzero lam they need omethod “RBFGS” or “cRBFGS”

  • grid_dim – size of the grid, for the DP2 method only (default = 7)

:rtype list containing :return f2n: aligned f2 :return gam: warping function :return q2n: aligned q2 (srsf)