addpath('/Volumes/iamdavecampbell/fdaM')
addpath('/Volumes/iamdavecampbell/fdaM/examples/lip')

%  Last modified 20 November 2006

%  -----------------------------------------------------------------------
%                       Lip Movement Data
%  -----------------------------------------------------------------------

%  ----------------  input the data  ------------------------

fid = fopen('lip.dat','rt');
lipmat = reshape(fscanf(fid,'%f'), [51, 20]);
nobs = size(lipmat,2);  %  number of replications

liptime  = linspace(0,1,51)';  %  sampling points for each curve

%  ----------  set up the b-spline basis object  ------------
%       use order 6 splines so we can look at acceleration

lipbasis = create_bspline_basis([0,1], 31, 6);

% we use 6 because we want acceleration (D^2) and so we will smooth the
% data with a penalty on two derivatives beyond (D^4) and we want that to
% be smooth so we use a 6th order spline


%  -----  create the fd object  ---------

%  set up the functional parameter object

Lfdobj   = int2Lfd(4);  %  penalize fourth derivative
lambda   = 1e-12;       %  use light smoothing


lipfdPar = fdPar(lipbasis, Lfdobj, lambda);

%  carry out the smoothing
 
[lipfd, df, gcv] = smooth_basis(liptime, lipmat, lipfdPar);
% note that the GCV is output from this function.  To minimize GCV for lambda in iPDA, we
% will need to collect the GCV values after iterating through the iPDA
% process, then we will get a value of GCV for each lambda


%  add names to dimensions

lipfd_fdnames{1} = 'Normalized time';
lipfd_fdnames{2} = 'Replications';
lipfd_fdnames{3} = 'mm';
lipfd = putnames(lipfd, lipfd_fdnames);

%  ------------  plot data and fit you will have to hit enter to see the next data fit added to the plot ----------------

subplot(1,1,1)
hold on
plotfit_fd(lipmat, liptime, lipfd)
hold off
%  ------------- plot residuals you will have to hit enter to see the next
%  data fit added to the plot ----------------
subplot(1,1,1)
hold on
plotfit_fd(lipmat, liptime, lipfd,  1)
hold off
%  ---------  summarize and plot the functions and their accelerations  -----
%  ---------  summarize and plot the functions and their accelerations  -----

for lp=0:5;
    subplot(3,2,lp+1);plot(lipfd,lp);
    title(strcat('\fontsize{5}derivative ',num2str(lp),' of the Lip functional observation'));
end

%NOTE: derivative 0 is the smooth fit to our data
% we are interested in modelling D^2Lip, the second derivative
% We smoothed with a penalty on the 4th derivative
% We are using a 6th order spline so the 5th derivative of the smooth is
% piecewise constant with jumps at the knot points.



%  -------------  principal differential analysis  -------------------

%  compute the weight functions by the basis expansion method

difeorder  = 2;  %  order of equation, so that the model that we will use 
% is D^2Y(t)= beta_0(t)Y(t)+beta_1(t)DY(t)

awtcell = {};% these are only used if we have a forcing function or other system input
ufdcell = {};% these are only used if we have a forcing function or other system input



nwbasis = 21;  %use about 1 per functional data observation
wbasis = create_bspline_basis([0,1], nwbasis);
wcoef0 = zeros(nwbasis,1);
betafd = fd(wcoef0, wbasis); 
betafdPar = fdPar(betafd);
bwtcell{1,1} = betafdPar;
bwtcell{1,2} = betafdPar;

xfdcell{1} = lipfd;

%  carry out principal differential analysis

[bfdcell, resfdcell, afdcell] = ...
    pda_fd(xfdcell, bwtcell, awtcell, ufdcell, difeorder);

%  plot the weight functions thiese are the resulting beta(t) functions
% NOTE1: to do this you either need to modify the FDA Matlab functions or
% get the file plot.m from the course file space and put it in the folder
% called @fdPar you will need to replace the existing file of that name.
% NOTE2: you will need to click the plot and/or hit enter to get the next
% plot to appear
figure(1)
for j=1:2
    subplot(2,1,j)
    plot(bfdcell{1,j});
    ylabel(['Weight function \beta',num2str(j-1),'(t)']);
end

%  set up a linear differential operator 

wcoef  = [getcoef(getfd(bfdcell{1})), getcoef(getfd(bfdcell{2}))];
wfd    = fd(wcoef,wbasis);
bwtcell = fd2cell(wfd);
Lfdobj = Lfd(difeorder, bwtcell);

% now we could re-smooth the data using this linear differential operator
% Lfdobj that contains our model and the model betas, and then update the
% betas.
% To do this we might take out some of the plots and run the above as a for
% loop where we say for lp=1:10 on line #35 and we say end on line #129.
% This would give us iterative PDA!




















%% Optional - only work on this as a hobby, not needed for the assignment

% The model we have have involves time varying coefficients, to compute an
% ode45 solution to this problem we have to use a give it our time varying
% coefficients, this is built into the function 'derivs' that we use below
%

%  solve equation for lip observation #1 and lip observation #2

ystart = eye(2);
[tp1, yp1] = ode45(@derivs, liptime, ystart(:,1), [], wfd);
[tp2, yp2] = ode45(@derivs, liptime, ystart(:,2), [], wfd);

%  plot the two solutions from ODE45

umat = [yp1(:,1),yp2(:,1)];
subplot(2,1,1)
plot(liptime, umat, liptime, zeros(51,1), ':'), title('2 lip Functions from using the ODE45 solver and our time varying coefficients');
Dumat = [yp1(:,2),yp2(:,2)];
subplot(2,1,2)
plot(liptime, Dumat, liptime, zeros(51,1), ':'), title('Derivatives of the above');

% do this part to see the residual to the model (the forcing function)


%  compute forcing functions

force = eval_fd(liptime, lipfd, lipLfd);

%  plot the forcing functions for each curve
subplot(1,1,1)
plot(liptime, force);
axis([0,1,-1e3,1e3]);
title('these are the forcing functions, they are like the residuals between the model and the smooth observations')
%  plot the mean forcing function along with second deriv.

forcemean = mean(force,2);
D2mean = eval_fd(mean(lipfd),liptime,int2Lfd(2));

subplot(1,1,1)
plot(liptime, forcemean, '-', liptime, D2mean, ':')
axis([0,1,-700,700])
title('The mean forcing function and the fitted values of D^2Lip')



