fix horizontal fpca integration problem, adding plotting recipes for pca

src/fPCA.jl

@@ -88,32 +88,31 @@ Calculates horizontal functional principal component analysis on aligned data
 function horiz_fPCA(warp_data::fdawarp; no=1)
     gam = warp_data.gam;
     mu, gam_mu, psi, vec1 = sqrt_mean(gam);
-    tau = collect(1:6);
+    tau = collect(1:5);
     m = length(mu);
     n = length(tau);
-    TT = size(vec1,1)+1
+    TT = size(vec1,1)
 
     # TFPCA
     K = Statistics.covm(vec1, mean(vec1,dims=2), 2);
 
     U, s, V = svd(K);
     vm = mean(vec1, dims=2);
 
-    gam_pca = Array{Float64}(undef, (n, m+1, no));
-    psi_pca = Array{Float64}(undef, (n, m, no));
+    gam_pca = Array{Float64}(undef, (m, n, no));
+    psi_pca = Array{Float64}(undef, (m, n, no));
     for j in 1:no
         for k in tau
             v = (k-3)*sqrt(s[j]).*U[:,j];
             vn = norm(v) / sqrt(TT);
             if vn < 0.0001
-                psi_pca[k, :, j] = mu;
+                psi_pca[:, k, j] = mu;
             else
-                psi_pca[k, :, j] = cos(vn).*mu + sin(vn).*v/vn;
+                psi_pca[:, k, j] = cos(vn).*mu + sin(vn).*v/vn;
             end
 
-            tmp = zeros(TT);
-            tmp[2:TT] = cumsum(psi_pca[k,:,j] .* psi_pca[k,:,j],dims=2);
-            gam_pca[k,:,j] = (tmp .- tmp[1]) ./ (tmp[end] - tmp[1]);
+            gam0 = cumtrapz(collect(LinRange(0,1,TT)), psi_pca[:,k,j].*psi_pca[:,k,j]);
+            gam_pca[:,k,j] = norm_gam(gam0)
         end
     end
 

src/plots.jl

@@ -62,3 +62,79 @@
     end
 
 end
+
+@recipe function f(h::vfpca; no=1:3)
+    x = h.time
+    M, N = size(h.q_pca)
+
+    # set up the subplots
+    legend := false
+    layout := (2,3)
+
+    @series begin
+        subplot := 1
+        title := @sprintf "f domain: PD %d" no[1]
+        x, h.f_pca[:,:,no[1]]
+    end
+
+    @series begin
+        subplot := 2
+        title := @sprintf "f domain: PD %d" no[2]
+        x, h.f_pca[:,:,no[2]]
+    end
+
+    @series begin
+        subplot := 3
+        title := @sprintf "f domain: PD %d" no[3]
+        x, h.f_pca[:,:,no[3]]
+    end
+
+    @series begin
+        subplot := 4
+        title := @sprintf "q domain: PD %d" no[1]
+        x, h.q_pca[1:(M-1),:,no[1]]
+    end
+
+    @series begin
+        subplot := 5
+        title := @sprintf "q domain: PD %d" no[2]
+        x, h.q_pca[1:(M-1),:,no[2]]
+    end
+
+    @series begin
+        subplot := 6
+        title := @sprintf "q domain: PD %d" no[3]
+        x, h.q_pca[1:(M-1),:,no[3]]
+    end
+
+end
+
+@recipe function f(h::hfpca, no=1:3)
+    M, N, no1 = size(h.gam_pca)
+    x = LinRange(0,1,M)
+
+    # set up the subplots
+    legend := false
+    layout := (1,3)
+    aspect_ratio := :equal
+
+    @series begin
+        subplot := 1
+        title := @sprintf "gam: PD %d" no[1]
+        x, h.gam_pca[:,:,no[1]]
+    end
+
+    @series begin
+        subplot := 2
+        title := @sprintf "gam: PD %d" no[2]
+        x, h.gam_pca[:,:,no[2]]
+    end
+
+    @series begin
+        subplot := 3
+        title := @sprintf "gam: PD %d" no[3]
+        aspect_ratio := :equal
+        x, h.gam_pca[:,:,no[3]]
+    end
+
+end

test/runtests.jl

@@ -40,10 +40,10 @@ out = srsf_align(f, timet);
 out1 = align_fPCA(f, timet, MaxItr=2);
 
 # test vert_fPCA
-out1 = vert_fPCA(out);
+out1 = vert_fPCA(out, no=3);
 
 # test horiz_fPCA
-out1 = horiz_fPCA(out);
+out1 = horiz_fPCA(out, no=3);
 
 # test gauss_model
 out1 = gauss_model(out);