expose lambda to curve functions

jdtuck · Aug 21, 2023 · 591406f · 591406f
1 parent 832207f
commit 591406f
Show file tree

Hide file tree

Showing 5 changed files with 80 additions and 13 deletions.
diff --git a/Find_Rotation_and_Seed_unique.m b/Find_Rotation_and_Seed_unique.m
@@ -1,4 +1,4 @@
-function [q2best,Rbest,gamIbest] = Find_Rotation_and_Seed_unique(q1,q2,reparamFlag,rotFlag,closed,method)
+function [q2best,Rbest,gamIbest] = Find_Rotation_and_Seed_unique(q1,q2,reparamFlag,rotFlag,closed,lam,method)
 % FIND_ROTATION_AND_SEED_UNIQUE find roation and seed of two curves
 % -------------------------------------------------------------------------
 % find roation and seed of closed curve
@@ -11,6 +11,7 @@
 % q2: matrix (n,T) defining T points on n dimensional SRVF
 % reparamFlag: flag to calculate reparametrization (default F)
 % closed: flag if closed curve (default F)
+% lam: penalty (default 0.0)
 % method: controls which optimization method (default="DP") options are
 % Dynamic Programming ("DP") and Riemannian BFGS
 % ("RBFGSM")
@@ -24,15 +25,21 @@
     reparamFlag = true;
     rotFlag = true;
     closed = false;
+    lam = 0.0;
     method = 'DP';
 elseif nargin < 4
     rotFlag = true;
     closed = false;
+    lam = 0.0;
     method = 'DP';
 elseif nargin < 5
     closed = false;
+    lam = 0.0;
     method = 'DP';
 elseif nargin < 6
+    lam = 0.0;
+    method = 'DP';
+elseif nargin < 7
     method = 'DP';
 end
 
@@ -62,7 +69,7 @@
     if(reparamFlag)
 
         if norm(q1-q2n,'fro') > 0.0001
-            gam = optimum_reparam_curve(q2,q1,method);
+            gam = optimum_reparam_curve(q2,q1,lam,method);
             gamI = invertGamma(gam);
             gamI = (gamI-gamI(1))/(gamI(end)-gamI(1));
             p2n = q_to_curve(q2n);

diff --git a/fdacurve.m b/fdacurve.m
@@ -3,7 +3,59 @@
     % ---------------------------------------------------------------------
     % This class provides alignment methods for curves in R^n using the
     % SRVF framework
-
+    %
+    % Usage:  bj = fdacurve(beta, closed, N, scale, center)
+    %
+    % where:
+    %   beta: (n,T,K) matrix defining n dimensional curve on T samples with K curves
+    %   closed: true or false if closed curve
+    %   N: resample curve to N points (default = T)
+    %   scale: include scale (true/false (default))
+    %   center: center curve (true (default)/false)
+    %
+    %
+    % fdacurve Properties:
+    %   beta            % (n,T,K) matrix defining n dimensional curve on T samples with K curves
+    %   q               % (n,T,K) matrix defining n dimensional srvf on T samples with K srvfs
+    %   betan           % aligned curves
+    %   qn              % aligned srvfs
+    %   basis           % calculated basis
+    %   beta_mean       % karcher mean curve
+    %   q_mean          % karcher mean srvf
+    %   gams            % warping functions
+    %   v               % shooting vectors
+    %   C               % karcher covariance
+    %   s               % pca singular values
+    %   U               % pca singular vectors
+    %   pca             % principal directions
+    %   coef            % pca coefficients
+    %   closed          % closed curve if true
+    %   lambda          % warping penalty
+    %   qun             % cost function
+    %   samples         % random samples
+    %   gamr            % random warping functions
+    %   cent            % center
+    %   scale           % scale
+    %   E               % energy
+    %   len             % length of curves
+    %   len_q           % length of SRVFs
+    %   mean_scale      % mean length
+    %   mean_scale_q    % mean length SRVF
+    %   center          % centering of curves done
+    %
+    %
+    % fdacurve Methods:
+    %   fdacurve - class constructor
+    %   karcher_mean - find karcher mean
+    %   karcher_cov - find karcher covariance
+    %   shape_pca - compute shape pca
+    %   sample_shapes - sample shapes from generative model
+    %   plot - plot results and functions in object
+    %   plot_pca - plot shape pca
+    %
+    %
+    % Author :  J. D. Tucker (JDT) <jdtuck AT sandia.gov>
+    % Date   :  15-Mar-2020
     properties
         beta            % (n,T,K) matrix defining n dimensional curve on T samples with K curves
         q               % (n,T,K) matrix defining n dimensional srvf on T samples with K srvfs
@@ -20,6 +72,7 @@
         pca             % principal directions
         coef            % pca coefficients
         closed          % closed curve if true
+        lambda          % warping penalty
         qun             % cost function
         samples         % random samples
         gamr            % random warping functions
@@ -100,6 +153,7 @@
             % default options
             % option.reparam = true; % computes optimal reparamertization
             % option.rotation = true; % computes optimal rotation
+            % option.lambda = 0.0;  % penalty
             % option.parallel = 0; % turns on MATLAB parallel processing (need
             % parallel processing toolbox)
             % option.closepool = 0; % determines wether to close matlabpool
@@ -115,6 +169,7 @@
             if nargin < 2
                 option.reparam = true;
                 option.rotation = true;
+                option.lambda = 0.0;
                 option.parallel = 0;
                 option.closepool = 0;
                 option.MaxItr = 20;
@@ -176,7 +231,7 @@
                         q1=obj.q(:,:,i);
 
                         % Compute shooting vector from mu to q_i
-                        [qn_t,~,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.method);
+                        [qn_t,~,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.lambda,option.method);
                         qn_t = qn_t/sqrt(InnerProd_Q(qn_t,qn_t));
 
                         gamma(:,i) = gamI;
@@ -213,7 +268,7 @@
                         q1=obj.q(:,:,i);
 
                         % Compute shooting vector from mu to q_i
-                        [qn_t,~,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.method);
+                        [qn_t,~,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.lambda,option.method);
                         qn_t = qn_t/sqrt(InnerProd_Q(qn_t,qn_t));
 
                         gamma(:,i) = gamI;
@@ -297,7 +352,7 @@
                     beta1 = betan1(:,:,i);
 
                     % Compute shooting vector from mu to q_i
-                    [~,R,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.method);
+                    [~,R,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.lambda,option.method);
                     beta1 = R*beta1;
                     beta1n = warp_curve_gamma(beta1,gamI);
                     q1n = curve_to_q(beta1n);
@@ -316,7 +371,7 @@
                     beta1 = obj.beta(:,:,i);
 
                     % Compute shooting vector from mu to q_i
-                    [~,R,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.method);
+                    [~,R,gamI] = Find_Rotation_and_Seed_unique(mu,q1,option.reparam,option.rotation,obj.closed,option.lambda,option.method);
                     beta1 = R*beta1;
                     beta1n = warp_curve_gamma(beta1,gamI);
                     q1n = curve_to_q(beta1n);
@@ -336,6 +391,7 @@
             obj.beta_mean = betamean;
             obj.q_mean = mu;
             obj.gams = gamma;
+            obj.lambda = lambda;
             obj.v = v1;
             obj.qun = sumd(1:iter);
             obj.E = normvbar(1:iter-1);

diff --git a/fdawarp.m b/fdawarp.m
@@ -23,7 +23,7 @@
     %    psi - srvf of warping functions
     %    stats - alignment statistics
     %    qun - cost function
-    %    lambda - lambda
+    %    lambda - warping penalty
     %    method - optimization method
     %    gamI - invserse warping function
     %    rsamps - random samples

diff --git a/optimum_reparam_curve.m b/optimum_reparam_curve.m
@@ -1,4 +1,4 @@
-function gam = optimum_reparam_curve(q1,q2,method)
+function gam = optimum_reparam_curve(q1,q2,lam,method)
 % OPTIMUM_REPARAM_CURVE Calculates Warping for two SRVFs
 % -------------------------------------------------------------------------%
 % This function aligns two SRVF functions using Dynamic Programming
@@ -9,13 +9,17 @@
 % Input:
 % q1: srvf of function 1
 % q2: srvf of function 2
+% lam: penalty (default 0.0)
 % method: controls which optimization method (default="DP") options are
 % Dynamic Programming ("DP") and Riemannian BFGS
 % ("RBFGSM")
 %
 % Output:
 % gam: warping function
 if nargin < 3
+    lam = 0.0;
+    method = 'DP';
+elseif nargin < 4
     method = 'DP';
 end
 
@@ -24,7 +28,7 @@
 
 switch upper(method)
     case 'DP'
-        gam0 = DynamicProgrammingQ(q2,q1,0,0);
+        gam0 = DynamicProgrammingQ(q2,q1,lam,0);
     case 'RBFGSM'
         t1 = linspace(0,1,size(q1,2));
         options.verbosity=0;

diff --git a/pairwise_align.m b/pairwise_align.m
@@ -1,4 +1,4 @@
-function [f2_align, gam] = pairwise_align(f1,f2,time,lam)
+function [f2_align, gam] = pairwise_align(f1,f2,time,lambda)
 % PAIRWISE_ALIGN Align two functions
 % -------------------------------------------------------------------------
 % This function aligns two functions using SRSF framework. It will align f2
@@ -19,11 +19,11 @@
     f1 double
     f2 double
     time double
-    lam = 0
+    lambda = 0
 end
 q1 = f_to_srvf(f1,time);
 q2 = f_to_srvf(f2,time);
 
-gam = optimum_reparam(q1,q2,time,lam,'DP1',0.0,0.0,0.0);
+gam = optimum_reparam(q1,q2,time,lambda,'DP1',0.0,0.0,0.0);
 
 f2_align = warp_f_gamma(f2,gam,time);