%=========== Description: ===========
%
%            This is the main program to perform Functional Single Index 
%            Models (FSIM). Details can be found in Jiang and Wang (2011).
%	     The principal component scores are estimated through conditional 
%	     expectation as in PACE to accommodate sparse longitudinal data.
%
%====== Usage: ======
%
% function [no_opt,sigma,lambda,eigen,xi_est,xi_var,mu,bw_mu,xcov,bw_xcov,
%           xcovfit, AIC, BIC, FVE, y_pred, y_predOrig, out21]...
%        = FSIMs(y, z, t, bwmu, bwxcov, mfc, selection_k, FVE_threshold, 
%                error, ngrid)
%
%================ Input Arguments: ================
%
%Input y:          1*n cell array, y{i} is the vector of measurements for
%                  the ith subject, i=1,...,n.
%
%Input z:          n*1 cell array, z{i} is the p-dimensional longitudinal 
%                  covariates of the ith subject, i=1,...,n.
%
%Input t:          1*n cell array, t{i} is the vector of time points for
%                  the ith subject for which corresponding measurements y{i}
%                  are available, i=1,...,n.             
%
%Input bwmu:       1*2 vector, bandwidths for covariance surface used for
%                  smoothing of mu(t,z) bwmu(1)==0 & bwmu(2)==0: use
%                  leave-one-curve-out cross-validation (CV)  [Default]
%                  bwmu(1)>0  & bwmu(2)>0: user-specified bandwidths.
%                  bwmu(1): the bandwidth for time coordinate bwmu(2): the
%                  bandwidth for covariate coordinate
%
%Input bwxcov:     1*1 vector, bandwidth for covariance surface used for
%                  smoothing of the covariance G(t,s) 
%                  bwxcov(1)==0: use m-fold cross-validation (CV) [Default]
%                  bwxcov(1)>0 : user-specified bandwidths.
%
%Input mfc:        integer(>=2), number of fold used for CV when choosing
%                  bandwidths for smoothing the covariance surface. 
%
%Input selection_k: the method of choosing the number of principal
%                   components K.
%                   'AIC':  use AIC criterion with pseudo-likelihood of
%                           measurements (marginal likelihood).
%                   'BIC':  use BIC criterion with pseudo-likelihood of
%                           measurements (marginal likelihood).
%                                                               [Default]
%                   'FVE' (fraction of variance explained) : use scree plot
%                           approach to select number of principal
%                           components),
%
%                    positive integer K: user-specified number of principal 
%                                        components ( K <= ngrid1 )
%
%                   Note: BIC and FVE produce the most parsimonious models.
%
%Input FVE_threshold:  a positive number that is between 0 and 1 [Default
%                      is 0.80] It is used with the option selection_k = 'FVE' to
%                      select the number of principal components that
%                      explain at least "FVE_threshold" of total variation
%                      (the fraction of variance explained).
%
%Input error:     0, no additional measurement error assumed.
%                 1, with measurement error (constant variance)[Default] 
%                 2, with measurement error (variance varies with time)
%
%Input ngrid:     integer, number of support points in each direction of 
%                 covariance surface when performing principal component 
%                 analysis ( ngrid > K).                [Default is 51.]
%
%
%================= Output Arguments: =================
%
%Output noeig:    integer, automatically or subjectively selected value of
%                 K, the number of selected eigenfunctions.
%
%Output Esig:     a function measured at tout if error=2; 
%                 scalar, estimate of measurement error variance if
%                 error=1, while it is []; 
%                 0 if error=0.
%
%Output lambda:   1*K vector, estimated eigenvalues (variances of
%                 functional principal components scores).
%
%Output eigen:    ngrid*K matrix, estimated eigenfunctions,
%                 valued at out21, ngrid of the pooled distinct time points
%                 with ascending order of all subjects, phi is an
%                 interpolated version of eigen at out1
%
%Output xi_est:   1*n cell, predictions for random coeffecients (PC
%                 scores) for n subjects.
%
%Output beta:     1*p vector, coefficients of the single index
%
%Output muout:    ngrid*ngrid matrix, estimated mean functions valued
%
%Output bw_mu:    scalar(>0), automatically or subjectively selected
%                 bandwidth for smoothing mean curve.
%
%Output ecov:     ngrid*ngrid matrix, smoothed covariance surface
%                 (diagnal removed), corresponding to out21
%
%Output bw_xcov:  1*2 vector(>0), automatically or subjectively selected
%                 bandwidths for smoothing covariance surface.
%
%Output yfit:     1*n cell array
%
%Output tout:     1*ngrid vector, a grid of time points for which the
%                 smoothed covariance surface assumes values, i.e., ngrids
%                 from out1.
%