% Example for WFPCA.m

%%%%%%%%%%%%%%%%% Generate data %%%%%%%%%%%%%%%%%%%%%%%%
% balanced data with N=50 curves and each with M=20 measurements
% the model is Y_i(t)=X_i(h_i^{-1}(t))+error 
% X_i(t)=mu(t) are the same for all i.
% mu(t)=exp(-10*(t-0.5).^2) with one peak
% the inverse warping functions h_i^{-1} are generated by the area under random linear splines
% noises from N(0,0.1) are added at each measurement
N = 50;
eps = 0.1;
M = 20;
me = @(t)exp(-10*(t-0.5).^2);
T = 0:(1/(M-1)):1;
Treal = [];
Y = [];
for i=1:N
    tempT = spapi(2,[0:0.2:1],rand(1,6));
    tempT = fnint(tempT);
    Treal(i,:) = fnval(tempT,T)/fnval(tempT,1);
    Y(i,:) = me(Treal(i,:))+eps*randn(1,M);
end

%%%%%%%%%%%%%%%% use WFPCA function %%%%%%%%%%%%%%%%%%%%%%%%%%%
d = path;
isExist = regexp(d, 'PACE');
if isempty(isExist) == 1
addpath(genpath('../PACE/'));
end
[X,aligned,h,hinv]=WFPCA(Y,T);


%%%%%%%%%%%%%%%% use WFPCA with Eigen option %%%%%%%%%%%%%%%%%%
%% Please note that Eigen library must be installed for this option,
%% which yields faster computation (see help file)
[X,aligned,h,hinv]=WFPCA(Y,T,[],[],[],[],[],[],1);


%%%%%%%%%%%%% analysis and plots based on the output %%%%%%%%%%%%%%%%%%%%
% plot aligned Y_i(t) (with noises) comparing to the simulated X_i(t) (without noises) 
figure 
subplot(1,2,1)
hold on
for i=1:N
    plot(T,Y(i,:),'-r')
end
xlabel('t')
ylabel('Y_i(t)')
title('simulated Y_i(t) with warping and noises')
hold off
subplot(1,2,2)
hold on
for i=1:size(aligned,1)
    plot(T,aligned(i,:),'-r')
end
xlabel('t')
ylabel('Y*_i(t)')
title('aligned Y_i(t)')
hold off
% plot the cross-sectional mean of Y_i(t), cross-sectional mean of aligned Y_i(t) and the true mean 
figure
hold on
plot(T,mean(Y),'-r')
plot(T,mean(aligned),'-b')
plot(T,me(T),'-m')
legend('mean of Y_i','mean of aligned Y_i','true mean','location','best')
xlabel('t')
title('cross-sectional means')
hold off
% plot the estimated inverse warping functions comparing to the simulated ones
figure
subplot(1,2,1)
hold on
for i=1:size(Treal,1)
    plot(T,Treal(i,:),'-r')
end
hold off
xlabel('t')
ylabel('Y(t)')
title('simulated inverse warping functions')
subplot(1,2,2)
hold on
for i=1:size(hinv,1)
    plot(T,hinv(i,:),'-b')
end
hold off
xlabel('t')
ylabel('Y(t)')
title('estimated inverse warping functions')
