% pop_fda_5 - Perform ERP analysis using Function-on-Function Regression (FOF).  
%             This function smooths EEG curves using B-spline and Fourier 
%             basis functions, selects the optimal smoothing penalty using 
%             Generalized Cross-Validation (GCV), and computes FOF regression 
%             coefficients.  
%             It also estimates and plots bootstrap-based p-values and 
%             confidence intervals for statistical inference. 
%
% Usage:
%   pop_fda_5(EEG);     % Opens a GUI prompting the user to select analysis options
%
% Inputs:
%   EEG - EEGLAB dataset containing epoched data with ICA weights computed.
%
% Outputs:
%   In version 0.1, the function does not return any output, except for a  
%   command string appended to the EEGLAB history.
%
% Version History:
%   Version 0.1 (2025)
%       - Initial release.
%       - GUI for Function-on-Function Regression added.
%       - Integrated smoothing, FOF regression, and bootstrap computation  
%         routines from the fdaM toolbox.
%       - Added visualization of bootstrap p-values and confidence intervals.
%
% Author:
%   Mohammad Fayaz. FDA plug-in GUI code.
%   E-mail: Mohammad.Fayaz.89@gmail.com
%
%   This function calls routines (with permission) from the fdaM toolbox:
%
%       J.O. Ramsay, Giles Hooker, and Spencer Graves,
%       "Functional Data Analysis with R and MATLAB",
%       Springer, 2009.
%
%       Repository (MATLAB): 
%           https://www.psych.mcgill.ca/misc/fda/downloads/FDAfuns/Matlab/fdaM.zip
%
%       R alternative: The 'fda' package on CRAN 
%           https://cran.r-project.org/web/packages/fda/index.html
%           License: GPL-2 | GPL-3.
%
% Example:
%   % Invoke Function-on-Function Regression (GUI mode):
%   pop_fda_5(EEG);
%
% Copyright:
%   © 2025 Mohammad Fayaz – Version 0.1
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.



function com = pop_fda_5(EEG, varargin)
%%%% Tools -> Regreesion -> FoF
fprintf('Regression Analysis \n');
fprintf('Function on Function Regression (FoF) \n');


com = '';
if nargin < 1
    help pop_sourcereconstruction;
    return
end

if nargin < 2
    %%%% GUI 
    g = [1 0.5 0.5];
    geometry = { g g g g g g g g g g g g g g g g g g g g g g g g g g g};
    uilist = { ...
      { 'Style', 'text', 'string', 'Channel Selection 1 ', 'fontweight', 'bold' } {} {} ...                       % 1
      { 'Style', 'text', 'string', '   Channel number', 'FontAngle', 'italic'  } ...                             % 2
      { 'Style', 'edit', 'string', '7' 'tag' 'SelectedChanIndex_1'} {} ...
      ...
      { 'Style', 'text', 'string', '   Time limits [min max] (msec)', 'FontAngle', 'italic' } ...                % 3
      { 'Style', 'edit', 'string', '-0.2 0.800 ' 'tag' 'TIS_1' } {} ...
      ...
      { 'Style', 'text', 'string', 'Channel Selection 2 ', 'fontweight', 'bold' } {} {} ...                       % 4
      { 'Style', 'text', 'string', '   Channel number', 'FontAngle', 'italic'  } ...                             % 5
      { 'Style', 'edit', 'string', '25' 'tag' 'SelectedChanIndex_2'} {} ...
      ...
      { 'Style', 'text', 'string', '   Time limits [min max] (msec)', 'FontAngle', 'italic' } ...                % 6         
      { 'Style', 'edit', 'string', '-0.2 0.800 ' 'tag' 'TIS_2' } {} ...
      ...    
      { 'Style', 'text', 'string', 'Preprocessing (Channel 1) ', 'fontweight', 'bold' } {} {} ...                % 7               
      { 'Style', 'text', 'string', '   Choose Basis Function (Default: B-Spline)', 'FontAngle', 'italic' } ...   % 8
      { 'Style', 'popupmenu', 'string', 'B-Spline|Fourier' , 'tag' 'TypeBasisFunc_1' } {} ...
      { 'Style', 'text', 'string', '   Number of Basis', 'FontAngle', 'italic' } ...                             % 9
      { 'Style', 'edit', 'string', '120' 'tag' 'NB_1'  } {} ...
      ...
      { 'Style', 'text', 'string', '   Order of Basis','FontAngle', 'italic' } ...                               % 10
      { 'Style', 'edit', 'string', '6' 'tag' 'NORD_1' } {} ...
      { 'Style', 'text', 'string', '   Parameter Estimation', 'FontAngle', 'italic' } ...                        % 11
      { 'Style', 'checkbox', 'string' 'GCV' 'value' 1 'tag' 'TickGCV_1' } ...
      { 'Style', 'checkbox', 'string' 'Plot GCV' 'value' 1 'tag' 'PlotGCV_1' } ...
      { 'Style', 'text', 'string', 'Preprocessing (Channel 2) ', 'fontweight', 'bold' } {} {} ...                % 12
      { 'Style', 'text', 'string', '   Choose Basis Function (Default: B-Spline)', 'FontAngle', 'italic' } ...   % 13
      { 'Style', 'popupmenu', 'string', 'B-Spline|Fourier' , 'tag' 'TypeBasisFunc_2' } {} ...
      { 'Style', 'text', 'string', '   Number of Basis', 'FontAngle', 'italic' } ...                             % 14
      { 'Style', 'edit', 'string', '120' 'tag' 'NB_2'  } {} ...
      ...
      { 'Style', 'text', 'string', '   Order of Basis','FontAngle', 'italic' } ...                               % 15
      { 'Style', 'edit', 'string', '6' 'tag' 'NORD_2' } {} ...
      { 'Style', 'text', 'string', '   Parameter Estimation', 'FontAngle', 'italic' } ...                        % 16
      { 'Style', 'checkbox', 'string' 'GCV' 'value' 1 'tag' 'TickGCV_2' } ...
      { 'Style', 'checkbox', 'string' 'Plot GCV' 'value' 1 'tag' 'PlotGCV_2' } ...
      { 'Style', 'text', 'string', 'Function on Function Regression  (FOF)', 'fontweight', 'bold' } {} {} ...      % 17
      { 'Style', 'text', 'string', '   Set Up Coefficient Functions',  'FontAngle', 'italic' } ...
      { 'Style', 'checkbox', 'string' 'Diagnostic Plots' 'value' 1 'tag' 'FOF_DiagPlot ' }  {} ...
      { 'Style', 'text', 'string', '   Choose Basis Function (Default: B-Spline)', 'FontAngle', 'italic' } ...     % 19
      { 'Style', 'popupmenu', 'string', 'B-Spline|Fourier' , 'tag' 'TypeBasisFunc_3' } {} ...
      { 'Style', 'text', 'string', '   Number of Basis', 'FontAngle', 'italic' } ...                               % 20
      { 'Style', 'edit', 'string', '8' 'tag' 'NB_3'  } {} ...
      { 'Style', 'text', 'string', '   Order of Basis','FontAngle', 'italic' } ...                               % 21
      { 'Style', 'edit', 'string', '6' 'tag' 'NORD_3' } {} ...
      { 'Style', 'text', 'string', '   Plot Coefficients', 'FontAngle', 'italic' } ...
      { 'Style', 'checkbox', 'string' 'Intercept' 'value' 1 'tag' 'FOF_Beta0' }  ...
      { 'Style', 'checkbox', 'string' 'Beta' 'value' 1 'tag' 'FOF_Beta1' }  ...
      { 'Style', 'text', 'string', '   Bootstrap (default: 200 resamples) ', 'FontAngle', 'italic' } ...
      { 'Style', 'checkbox', 'string' 'Yes' 'value' 1 'tag' 'FOF_BS' } {} ...
      { 'Style', 'text', 'string', '   Bootstrap Iterations', 'FontAngle', 'italic' } ...
      { 'Style', 'edit', 'string', '200' 'tag' 'nBoot'  } {} ...
      { 'Style', 'text', 'string', '   Alpha', 'FontAngle', 'italic' } ...
      { 'Style', 'edit', 'string', '0.05' 'tag' 'ValueAlpha'  } {} ...
      { 'Style', 'text', 'string', 'General Setting', 'fontweight', 'bold' } {} {} ...                                       % 21
      { 'Style', 'text', 'string', '   Plot Status',  'FontAngle', 'italic' } ...                                            % 18
      { 'Style', 'checkbox', 'string' 'Plot All Smooth Curve' 'value' 0 'tag' 'AllSmoothPlot' } {} ...
      };

      %%%% Calling (GUI)
   [ tmp1 tmp2 strhalt structout ] = inputgui( geometry, uilist, ...
           'pophelp(''pop_newtimef'');', 'ERP Analysis | function-on-function Regression (FOF) for EEGLAB - Version 0.1');
    
    % Check if user canceled
    if isempty(structout)
        disp('User canceled the dialog. Exiting function.');
        return;  % stop the function safely
end  


%%%%%%%%% Bootstrap --> Residual Bootstrap or model based bootstrapp --->
%%%%%%%%% Improve it by Block Bootstrap



      %%% Import GUI items into the variables.
      %%%%% General Settings

            nChannel = EEG.nbchan;           % Number of Channenls
            nTrials = EEG.trials;            % Number of Trials
            ChanLocs = EEG.chanlocs;         % Channels Location
            Events = EEG.event;              % Events
            Epochs = EEG.epoch;              % Epochs
            Times = EEG.times;               % EEG Times
            EData = EEG.data;                % EEG Dataset
            ICAWeights = EEG.icawinv;        % ICA weights
            ICAData = EEG.icaact;            % ICA Source activation

            %%% Input by User (General)
            %%% Channel (Electrode) Selection
            % K = 3;
            %%% Channel 1 
            SelectedChanIndex_1 = str2num(structout.SelectedChanIndex_1);
            SelectedChan_1 = ChanLocs(:,SelectedChanIndex_1) ;
            NSelectedChan_1 = length(SelectedChanIndex_1);
            
            %%% Channel 2
            SelectedChanIndex_2 = str2num(structout.SelectedChanIndex_2);
            SelectedChan_2 = ChanLocs(:,SelectedChanIndex_2) ;
            NSelectedChan_2 = length(SelectedChanIndex_2);

            %%% Event Selection
            % SelectedEventIndex = M;
            % SelectedEvent = Events(:,SelectedEventIndex);

            %%% Time Interval Selection 
            %%%% Channel 1
            TIS_1 = str2num(structout.TIS_1) ;
            TimeIntervalSele_1 = TIS_1; %%% Time Interval Selection (defualt: -200, 800 ms)
            TimeLO_1 = min(TimeIntervalSele_1) * 1.0e+03 ;
            TimeUP_1 = max(TimeIntervalSele_1) * 1.0e+03 ;
            TimeIntervalSeleInput_1      = Times(Times(1,:) >= TimeLO_1 & Times(1,:) <= TimeUP_1);
            TimeIntervalSeleIndexInput_1 = Times(1,:) >= TimeLO_1 & Times(1,:) <= TimeUP_1;

            %%%% Channel 2
            TIS_2 = str2num(structout.TIS_2) ;
            TimeIntervalSele_2 = TIS_2; %%% Time Interval Selection (defualt: -200, 800 ms)
            TimeLO_2 = min(TimeIntervalSele_2) * 1.0e+03 ;
            TimeUP_2 = max(TimeIntervalSele_2) * 1.0e+03 ;
            TimeIntervalSeleInput_2      = Times(Times(1,:) >= TimeLO_2 & Times(1,:) <= TimeUP_2);
            TimeIntervalSeleIndexInput_2 = Times(1,:) >= TimeLO_2 & Times(1,:) <= TimeUP_2;
            
            %%% Analysis Data Set with selected time intervals
            ADS_EData_1 = EData(:,TimeIntervalSeleIndexInput_1(1,:),:);
            ADS_EData_2 = EData(:,TimeIntervalSeleIndexInput_2(1,:),:);

            TypeBasisFunc_1 = structout.TypeBasisFunc_1;
            TypeBasisFunc_2 = structout.TypeBasisFunc_2;
            
            %%%% FOF Setting
            %%%%%%% Coefficients 
            FOF_TypeBasisFunc  = structout.TypeBasisFunc_3     ; % B-Spline|Fourier for coefficient
            FOF_NB             = str2num(structout.NB_3)       ; % Number of basics for coefficient
            FOF_ORDER          = str2num(structout.NORD_3)     ; % Order for coefficient
            FOF_Plot_Beta0     = structout.FOF_Beta0           ; % Plot Beta0 
            FOF_Plot_Beta1     = structout.FOF_Beta1           ; % Plot Beta1
            FOF_BootStrap      = structout.FOF_BS              ; % Bootstrap 
            FOF_BootStrap_N    = str2num(structout.nBoot)      ; % Number of Boot
            FOF_BootStrap_Alpa = str2num(structout.ValueAlpha) ; % Alpha
            Diagnostic_Plot    = structout.FOF_DiagPlot        ; %Diagnostic_Plot for FOF

            %%% Analysis 
            %%%% Preprocessing
            %%% B-Spline
            %% Input
            %% Channel 1
            NumBasis_1 = str2num(structout.NB_1); % Number of Splines
            norder_1   = str2num(structout.NORD_1); % Number of Order
            TickGCV_1  = structout.TickGCV_1; % If GCV is ticked (TRUE-Default), the parameters are optimised with GCV, O.W. it is FALSE.
            PlotGCV_1 = structout.PlotGCV_1; % Plot GCV against Lambda Status (Deafult=0), If PlotGCV=1, it plots
            rng_1 = [min(TimeIntervalSeleInput_1),max(TimeIntervalSeleInput_1)];  % Range of B-Spline
            ADS_EData_Bspline_1_out_1=cell(NSelectedChan_1,1,nTrials);

            %% Channel 2
            NumBasis_2 = str2num(structout.NB_2); % Number of Splines
            norder_2   = str2num(structout.NORD_2); % Number of Order
            TickGCV_2  = structout.TickGCV_2; % If GCV is ticked (TRUE-Default), the parameters are optimised with GCV, O.W. it is FALSE.
            PlotGCV_2 = structout.PlotGCV_2; % Plot GCV against Lambda Status (Deafult=0), If PlotGCV=1, it plots
            rng_2 = [min(TimeIntervalSeleInput_2),max(TimeIntervalSeleInput_1)];  % Range of B-Spline
            ADS_EData_Bspline_1_out_2=cell(NSelectedChan_2,1,nTrials);

                               

            Bootstrap_Yes = structout.FOF_BS;
            AllSmoothPlot = structout.AllSmoothPlot; % Plot all smoothed curves (default = 0 ) , if AllSmoothPlot=1, plot all smoothed curves. 
            %% Anlaysis Code
            %% for nch = 1:NSelectedChan
            %% 
            
                SelectedChnInd_1 = SelectedChanIndex_1(NSelectedChan_1);
                SelectedChanLabels_1 = SelectedChan_1(NSelectedChan_1).labels;
                ADS_EData_Bspline_1 = ADS_EData_1(SelectedChnInd_1,:,:);

                if TypeBasisFunc_1 == 1
                 %%% Defining B-Spline Basis  
                 %%% Preparing Data
                 wbasis_1 = create_bspline_basis(rng_1, NumBasis_1, norder_1);
                 ADS_EData_Bspline_1 = double(squeeze(ADS_EData_Bspline_1));
                 ADS_EData_Bspline_1 = ADS_EData_Bspline_1';
                 ADS_EData_Bspline_1 = double(ADS_EData_Bspline_1');
                 TimeIntervalSeleInput_BSpline_1 = double(TimeIntervalSeleInput_1');
                end 
                if TypeBasisFunc_1 == 2
                 %%% Defining Fourier Basis  
                 %%% Preparing Data
                wbasis_1 = create_fourier_basis(rng_1, NumBasis_1);
                 ADS_EData_Bspline_1 = double(squeeze(ADS_EData_Bspline_1));
                 ADS_EData_Bspline_1 = ADS_EData_Bspline_1';
                 ADS_EData_Bspline_1 = double(ADS_EData_Bspline_1');
                 TimeIntervalSeleInput_BSpline_1 = double(TimeIntervalSeleInput_1');
                end
               
                 if TickGCV_1 == 1
                     fprintf('The GCV is selected for the first channel.\n')

                           n_1 = length(ADS_EData_Bspline_1(1,:));
                      loglam_1 = (-9:0.25:9)'; %  set up the range of log lambda values
                        nlam_1 = length(loglam_1);
                         Lfd_1 = 4;

                    dfsave_1  = zeros(nlam_1,1);
                    gcvsave_1 = zeros(nlam_1,n_1)';
                    MSEsave_1 = zeros(nlam_1,n_1)';

                    %  loop through the log lambda values (Modify Later this part for parameters)
                    for ilam_1=1:nlam_1
                        lambda_1 = 10^loglam_1(ilam_1);
                        hgtfdPar_1 = fdPar(wbasis_1, Lfd_1, lambda_1);
                        [hgtfd_1, df_1, gcv_1, coef_1,SSE_1] = smooth_basis(TimeIntervalSeleInput_BSpline_1,ADS_EData_Bspline_1,hgtfdPar_1);
                        accest_1 = eval_fd(TimeIntervalSeleInput_BSpline_1, hgtfd_1, 2);
                        dfsave_1(ilam_1)    = df_1;
                        gcvsave_1(:,ilam_1)   = gcv_1;  % check later for mean 
                        MSEsave_1(:,ilam_1)   = SSE_1;  % check later for mean 
                        fprintf(['Channel Number 1 : ', num2str(SelectedChnInd_1),' , ', num2str(ilam_1),'-Log Lambda : ',num2str(loglam_1(ilam_1)), ' , Mean GCV of Epochs : ',num2str(round(mean(gcvsave_1(:,ilam_1)),4)),' , Mean SSE of Epochs: ',num2str(round(mean(MSEsave_1(:,ilam_1)),4))  ,'\n']);                       

                    end

                    [MVALVE_1,MININDEX_1] =  min(mean(gcvsave_1));
                    LOGLAMMIN_1 = loglam_1(MININDEX_1);

                    lambdaSelected_1  = 10^loglam_1(MININDEX_1);
                    hgtfdPar_1 = fdPar(wbasis_1, Lfd_1, lambdaSelected_1);
                    ADS_EData_Bspline_1_1 = smooth_basis(TimeIntervalSeleInput_BSpline_1,ADS_EData_Bspline_1,hgtfdPar_1);
                    fprintf(['*** Channel Number 1: ', num2str(SelectedChnInd_1),' , Minimum log Lambda 1: ',num2str(LOGLAMMIN_1), ' , Minimum GCV 1: ', num2str(MVALVE_1)]);

                    if PlotGCV_1 == 1
                    fprintf('The plot GCV is selected for the first channel.\n')
                    %%% Plot GCV (Modify)
                    figure(NSelectedChan_1)
                    plot(loglam_1,mean(gcvsave_1)' )        
                    axis([ (min(loglam_1)-3) (max(loglam_1)+3) (min(mean(gcvsave_1)-5)) (max(mean(gcvsave_1)+5)) ])
                    subtitle(['Channel Name: ',SelectedChanLabels_1])
                    title(['The minimum GCV (', num2str(MVALVE_1) ,') at Lambda:' num2str(loglam_1(MININDEX_1)) ])
                    yy_1 = MVALVE_1;
                    xx_1 = loglam_1(MININDEX_1);
                   line([xx_1, xx_1], [min(mean(gcvsave_1))-20 ,max(mean(gcvsave_1))+20],'Color','red','LineStyle','--')
                   line([min(loglam_1)-10, max(loglam_1)+10], [yy_1, yy_1],'Color','red','LineStyle','--')
                   ylabel('GCV') 
                   xlabel('Logarithm Lambda') 
                   end

                    %%% Plot MSE
                    %%%% plot(loglam,mean(MSEsave)) --> Based on Values
                    %%%% plot(mean(MSEsave))  %% --> Based on Index

                 elseif TickGCV_1 == 0
                      fprintf('The GCV is not selected for the first channel.\n')
                    %%% Smoothing Data
                    [ADS_EData_Bspline_1_1,df_1_1, gcv_1_1,coef_1_1,SSE_1_1] = smooth_basis(TimeIntervalSeleInput_BSpline_1,ADS_EData_Bspline_1,wbasis_1);
                 LOGLAMMIN_1 = 0;
                 end

                    %%% Output
                    EEGDatafd_names_1{1} = 'Time (ms)';
                    EEGDatafd_names_1{2} = SelectedChanLabels_1 ;
                    EEGDatafd_names_1{3} = '\mu. volt';
                    ADS_EData_Bspline_1_1 = putnames(ADS_EData_Bspline_1_1, EEGDatafd_names_1);

                                      %%% Output
                   %% Save all fd Object 
                   ADS_EData_Bspline_1_out_1(NSelectedChan_1,1,:) = fd2cell(ADS_EData_Bspline_1_1);


            %% Channel 2
           
                SelectedChnInd_2 = SelectedChanIndex_2(NSelectedChan_2);
                SelectedChanLabels_2 = SelectedChan_2(NSelectedChan_2).labels;
                ADS_EData_Bspline_2 = ADS_EData_2(SelectedChnInd_2,:,:);

                if TypeBasisFunc_2 == 1
                 %%% Defining B-Spline Basis  
                 %%% Preparing Data
                 wbasis_2 = create_bspline_basis(rng_2, NumBasis_2, norder_2);
                 ADS_EData_Bspline_2 = double(squeeze(ADS_EData_Bspline_2));
                 ADS_EData_Bspline_2 = ADS_EData_Bspline_2';
                 ADS_EData_Bspline_2 = double(ADS_EData_Bspline_2');
                 TimeIntervalSeleInput_BSpline_2 = double(TimeIntervalSeleInput_2'); 
                end 
                if TypeBasisFunc_2 == 2
                 %%% Defining Fourier Basis  
                 %%% Preparing Data
                 wbasis_2 = create_fourier_basis(rng_2, NumBasis_2);
                 ADS_EData_Bspline_2 = double(squeeze(ADS_EData_Bspline_2));
                 ADS_EData_Bspline_2 = ADS_EData_Bspline_2';
                 ADS_EData_Bspline_2 = double(ADS_EData_Bspline_2');
                 TimeIntervalSeleInput_BSpline_2 = double(TimeIntervalSeleInput_2'); 
                end


                 if TickGCV_2 == 1
                      fprintf('The GCV is selected for the second channel.\n')

                           n_2 = length(ADS_EData_Bspline_2(1,:));
                      loglam_2 = (-9:0.25:9)'; %  set up the range of log lambda values
                        nlam_2 = length(loglam_2);
                         Lfd_2 = 4;

                    dfsave_2  = zeros(nlam_2,1);
                    gcvsave_2 = zeros(nlam_2,n_2)';
                    MSEsave_2 = zeros(nlam_2,n_2)';

                    %  loop through the log lambda values (Modify Later this part for parameters)
                    for ilam_2=1:nlam_2
                        lambda_2 = 10^loglam_2(ilam_2);
                        hgtfdPar_2 = fdPar(wbasis_2, Lfd_2, lambda_2);
                        [hgtfd_2, df_2, gcv_2, coef_2,SSE_2] = smooth_basis(TimeIntervalSeleInput_BSpline_2,ADS_EData_Bspline_2,hgtfdPar_2);
                        accest_2 = eval_fd(TimeIntervalSeleInput_BSpline_2, hgtfd_2, 2);
                        dfsave_2(ilam_2)    = df_2;
                        gcvsave_2(:,ilam_2)   = gcv_2;  % check later for mean 
                        MSEsave_2(:,ilam_2)   = SSE_2;  % check later for mean 
                        fprintf(['Channel Number 2 : ', num2str(SelectedChnInd_2),' , ', num2str(ilam_2),'-Log Lambda : ',num2str(loglam_2(ilam_2)), ' , Mean GCV of Epochs : ',num2str(round(mean(gcvsave_2(:,ilam_2)),4)),' , Mean SSE of Epochs: ',num2str(round(mean(MSEsave_2(:,ilam_2)),4))  ,'\n']);                       
                    end

                    [MVALVE_2,MININDEX_2] =  min(mean(gcvsave_2));
                    LOGLAMMIN_2 = loglam_2(MININDEX_2);

                    lambdaSelected_2  = 10^loglam_2(MININDEX_2);
                    hgtfdPar_2 = fdPar(wbasis_2, Lfd_2, lambdaSelected_2);
                    ADS_EData_Bspline_1_2 = smooth_basis(TimeIntervalSeleInput_BSpline_2,ADS_EData_Bspline_2,hgtfdPar_2);
                    fprintf(['*** Channel Number 2: ', num2str(SelectedChnInd_2),' , Minimum log Lambda 2: ',num2str(LOGLAMMIN_2), ' , Minimum GCV 2: ', num2str(MVALVE_2)]);

                    
                    
                    if PlotGCV_2 == 1
                    %%% Plot GCV (Modify)
                    fprintf('The GCV plot is selected for the second channel.\n')
                    figure(NSelectedChan_2)
                    plot(loglam_2,mean(gcvsave_2)' )        
                    axis([ (min(loglam_2)-3) (max(loglam_2)+3) (min(mean(gcvsave_2)-5)) (max(mean(gcvsave_2)+5)) ])
                    subtitle(['Channel Name: ',SelectedChanLabels_2])
                    title(['The minimum GCV (', num2str(MVALVE_2) ,') at Lambda:' num2str(loglam_2(MININDEX_2)) ])
                    yy_2 = MVALVE_2;
                    xx_2 = loglam_2(MININDEX_2);
                   line([xx_2, xx_2], [min(mean(gcvsave_2))-20 ,max(mean(gcvsave_2))+20],'Color','red','LineStyle','--')
                   line([min(loglam_2)-10, max(loglam_2)+10], [yy_2, yy_2],'Color','red','LineStyle','--')
                   ylabel('GCV') 
                   xlabel('Logarithm Lambda') 
                   end

                    %%% Plot MSE
                    %%%% plot(loglam,mean(MSEsave)) --> Based on Values
                    %%%% plot(mean(MSEsave))  %% --> Based on Index

                 elseif TickGCV_2 == 0
                        fprintf('The GCV is not selected for second channel.\n')
                    %%% Smoothing Data
                    [ADS_EData_Bspline_1_2,df_1_2, gcv_1_2,coef_1_2,SSE_1_2] = smooth_basis(TimeIntervalSeleInput_BSpline_2,ADS_EData_Bspline_2,wbasis_2);
                 LOGLAMMIN_2 = 0;
                 end

                    %%% Output
                    EEGDatafd_names_2{1} = 'Time (ms)';
                    EEGDatafd_names_2{2} = SelectedChanLabels_2 ;
                    EEGDatafd_names_2{3} = '\mu. volt';
                    ADS_EData_Bspline_1_2 = putnames(ADS_EData_Bspline_1_2, EEGDatafd_names_2);

                   %% Save all fd Object 
                   ADS_EData_Bspline_1_out_2(NSelectedChan_2,1,:) = fd2cell(ADS_EData_Bspline_1_2);

                   
                   %%%% Plot each Curve Separatley  
                   %%% plotfit_fd(ADS_EData_Bspline, TimeIntervalSeleInput_BSpline, ADS_EData_Bspline_1)
                   if AllSmoothPlot == 1
                   fprintf('Plot all smoothed function is selected. \n');
                   %%% Curve 1     
                   AllSmoothPlot_1 = getcoef(ADS_EData_Bspline_1_1);
                   AllSmoothPlot_range_1 = getbasisrange(getbasis(ADS_EData_Bspline_1_1));
                   AllSmoothPlot_nbasis_1 = getnbasis(getbasis(ADS_EData_Bspline_1_1));
                   AllSmoothPlot_nx_1 = max([501, 10*AllSmoothPlot_nbasis_1+1]);
                   AllSmoothPlot_x_1  = linspace(AllSmoothPlot_range_1(1),AllSmoothPlot_range_1(2),AllSmoothPlot_nx_1)';
                   AllSmoothPlot_Lfdobj_1 = int2Lfd(int2Lfd(0));                           
                   AllSmoothPlot_fdmat_1 = eval_fd(AllSmoothPlot_x_1,ADS_EData_Bspline_1_1, AllSmoothPlot_Lfdobj_1);
                   AllSmoothPlot_max_1 = max(AllSmoothPlot_fdmat_1);
                   AllSmoothPlot_min_1 = min(AllSmoothPlot_fdmat_1);
                   
                   MaxS_1 = max(AllSmoothPlot_max_1);
                   MinS_1 = min(AllSmoothPlot_min_1);
                   
                   if TypeBasisFunc_1 == 1 
                    MName_1 = 'B-Spline';
                   end
                   
                   if TypeBasisFunc_1 == 2 
                    MName_1 = 'Fourier';
                   end

                   MinT_1 = min(TimeIntervalSeleInput_BSpline_1);
                   MaxT_1 = max(TimeIntervalSeleInput_BSpline_1);
                   
                  %%% Curve 2  
                   AllSmoothPlot_2 = getcoef(ADS_EData_Bspline_1_2);
                   AllSmoothPlot_range_2 = getbasisrange(getbasis(ADS_EData_Bspline_1_2));
                   AllSmoothPlot_nbasis_2 = getnbasis(getbasis(ADS_EData_Bspline_1_2));
                   AllSmoothPlot_nx_2 = max([501, 10*AllSmoothPlot_nbasis_2+1]);
                   AllSmoothPlot_x_2  = linspace(AllSmoothPlot_range_2(1),AllSmoothPlot_range_2(2),AllSmoothPlot_nx_2)';
                   AllSmoothPlot_Lfdobj_2 = int2Lfd(int2Lfd(0));                           
                   AllSmoothPlot_fdmat_2 = eval_fd(AllSmoothPlot_x_2,ADS_EData_Bspline_1_2, AllSmoothPlot_Lfdobj_2);
                   AllSmoothPlot_max_2 = max(AllSmoothPlot_fdmat_2);
                   AllSmoothPlot_min_2 = min(AllSmoothPlot_fdmat_2);
                   
                   MaxS_2 = max(AllSmoothPlot_max_2);
                   MinS_2 = min(AllSmoothPlot_min_2);
                   
                   if TypeBasisFunc_2 == 1 
                    MName_2 = 'B-Spline';
                   end
                   
                   if TypeBasisFunc_2 == 2 
                    MName_2 = 'Fourier';
                   end

                   MinT_2 = min(TimeIntervalSeleInput_BSpline_2);
                   MaxT_2 = max(TimeIntervalSeleInput_BSpline_2);
                   
                   figure;
                   plot(AllSmoothPlot_x_1 ,AllSmoothPlot_fdmat_1)
                   title(['The Smoothing Functions of ', SelectedChanLabels_1 ,' is:  ', MName_1 ] )
                   line([MinT_1,MaxT_1], [0,0],'Color',[0 0 0]+0.05,'LineStyle','--')
                   subtitle(['Channel Name: ',SelectedChanLabels_1,' (Best Log Lambda:' num2str(LOGLAMMIN_1),')'])
                   axis([MinT_1-10 MaxT_1+10 MinS_1-25 MaxS_1+25])
                   xlabel('Time (ms)') 
                   ylabel('Mean')      
                   
                   figure;
                   plot(AllSmoothPlot_x_2 ,AllSmoothPlot_fdmat_2)
                   title(['The Smoothing Functions of ', SelectedChanLabels_2 ,' is:  ', MName_2 ] )
                   line([MinT_2,MaxT_2], [0,0],'Color',[0 0 0]+0.05,'LineStyle','--')
                   subtitle(['Channel Name: ',SelectedChanLabels_2,' (Best Log Lambda:' num2str(LOGLAMMIN_2),')'])
                   axis([MinT_2-10 MaxT_2+10 MinS_2-25 MaxS_2+25])
                   xlabel('Time (ms)') 
                   ylabel('Mean') 
                   end

                   
                   
                   %%% Function on Function Regression 
                   fprintf('The Function on Function Analysis Starts. \n');
                   % Set Up Coefficient Functions
                   if FOF_TypeBasisFunc == 1 
                     %  FOF_NB = 8;
                     %  FOF_ORDER = 6;
                       wbasis_FOF = create_bspline_basis(rng_2, FOF_NB, FOF_ORDER);
                       FOF_Beta0fdPar = fdPar(wbasis_FOF, 2, 1e-5);
                       FOF_Beta1fd    = bifd(zeros(FOF_NB), wbasis_FOF, wbasis_FOF);
                       FOF_Beta1fdPar = bifdPar(FOF_Beta1fd, 2, 1e3);
                       p = 2;
                       FOF_BetaCell = cell(p, 1);
                       FOF_BetaCell{1} = FOF_Beta0fdPar;
                       FOF_BetaCell{2} = FOF_Beta1fdPar;
                   end 

                    if FOF_TypeBasisFunc == 2
                       wbasis_FOF = create_fourier_basis(rng_2, FOF_NB);
                       FOF_Beta0fdPar = fdPar(wbasis_FOF, 2, 1e-5);
                       FOF_Beta1fd    = bifd(zeros(FOF_NB), wbasis_FOF, wbasis_FOF);
                       FOF_Beta1fdPar = bifdPar(FOF_Beta1fd, 2, 1e3);
                       p = 2;
                       FOF_BetaCell = cell(p, 1);
                       FOF_BetaCell{1} = FOF_Beta0fdPar;
                       FOF_BetaCell{2} = FOF_Beta1fdPar;
                   end 
             
                   %%%% XFD = ADS_EData_Bspline_1_1 , YFD = ADS_EData_Bspline_1_2
                   FOF_linmodStr_1 = linmod( ADS_EData_Bspline_1_1,  ADS_EData_Bspline_1_2, FOF_BetaCell);
                   %%%%  XFD = ADS_EData_Bspline_1_2 , YFD = ADS_EData_Bspline_1_1
                   FOF_linmodStr_2 = linmod(ADS_EData_Bspline_1_2, ADS_EData_Bspline_1_1, FOF_BetaCell);
                   fprintf('The Function on Function Regression is done ! \n');


                   FOF_beta1fd_1 = eval_bifd(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, ...
                          FOF_linmodStr_1.beta);

                   FOF_beta1fd_2 = eval_bifd(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, ...
                          FOF_linmodStr_2.beta);



                   if Diagnostic_Plot == 1

                % Get evaluation grid from linmod output
                tfine_1 = TimeIntervalSeleInput_BSpline_1;
                
                Y_obs_1 = eval_fd(ADS_EData_Bspline_1_1, tfine_1);    % [nfine x ncurves]
                Y_pred_1 = eval_fd(FOF_linmodStr_1.yhat, tfine_1);    % [nfine x ncurves]
                
                sqerr_1 = (Y_obs_1 - Y_pred_1).^2;
                RMSE_per_curve_1 = sqrt(mean(sqerr_1,1));          % per curve
                RMSE_overall_1   = sqrt(mean(sqerr_1(:)));         % overall
                
                % Compute R^2 for each curve
                varY_1 = var(Y_obs_1,0,1);
                R2_per_curve_1 = 1 - mean(sqerr_1,1)./varY_1;
                R2_overall_1   = 1 - mean(sqerr_1(:))/mean(varY_1);
                
                
                % Residuals
                residuals_1 = Y_obs_1 - Y_pred_1;
                
                % Create figure
                figure('Name', 'Residual Diagnostics', 'NumberTitle', 'off', 'Position', [100 100 1400 800]);
                
                %% 1. Residuals over time (all curves)
                subplot(2,3,1);
                plot(tfine_1, residuals_1, 'k');
                hold on;
                plot(tfine_1, mean(residuals_1, 2), 'r', 'LineWidth', 2);
                xlabel('Time');
                ylabel('Residual');
                title('Residuals over Time');
                legend('Residuals (curves)', 'Mean residual', 'Location', 'best');
                grid on;
                
                %% 2. Residual variance over time
                subplot(2,3,2);
                plot(tfine_1, var(residuals_1, 0, 2), 'b', 'LineWidth', 2);
                xlabel('Time');
                ylabel('Variance');
                title('Residual Variance over Time');
                grid on;
                
                %% 3. Histogram of residuals
                subplot(2,3,3);
                histogram(residuals_1(:), 30, 'FaceColor', [0.2 0.6 0.8]);
                xlabel('Residual');
                ylabel('Frequency');
                title('Histogram of Residuals');
                grid on;
                
                %% 4. Q-Q plot
                subplot(2,3,4);
                qqplot(residuals_1(:));
                title('Q-Q Plot of Residuals');
                
                %% 5. Residuals vs Predicted values (Homoscedasticity check)
                subplot(2,3,5);
                scatter(Y_pred_1(:), residuals_1(:), 15, 'filled', 'MarkerFaceAlpha', 0.6);
                xlabel('Predicted Values');
                ylabel('Residuals');
                title('Residuals vs Predicted');
                refline(0,0);
                grid on;
                
                %% 6. Autocorrelation of residuals
                subplot(2,3,6);
                My_Autocorr(residuals_1(:), 50);
                title('Residual Autocorrelation');
                
                sgtitle(['Residual Diagnostics for Functional Linear Model 1: ', ...
         '$\mathbf{y}(s)=$ ', SelectedChanLabels_2, ...
         ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_1], ...
         'Interpreter', 'latex', 'FontSize', 13, 'FontWeight', 'bold');
                
           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Each plot 
                        ncurves_1 = size(Y_obs_1, 2);
                        
                        % Fixed grid size: 5 x 5
                        nrows_1 = 5;
                        ncols_1 = 5;
                        plots_per_page = nrows_1 * ncols_1;
                        
                        % Number of pages needed
                        nPages = ceil(ncurves_1 / plots_per_page);
                        
                        for page = 1:nPages
                            figure('Name', sprintf('Observed vs Predicted (Page %d)', page), ...
                                   'NumberTitle', 'off', 'Position', [100 100 1400 800]);
                            
                            start_idx = (page-1) * plots_per_page + 1;
                            end_idx = min(page * plots_per_page, ncurves_1);
                            
                            for i = start_idx:end_idx
                                subplot(nrows_1, ncols_1, i - start_idx + 1);
                                plot(tfine_1, Y_obs_1(:,i), 'b', 'LineWidth', 1.2); hold on;
                                plot(tfine_1, Y_pred_1(:,i), 'r--', 'LineWidth', 1.2);
                                xlabel('Time');
                                ylabel('Amplitude');
                                title(sprintf('Curve %d', i), 'FontSize', 8);
                                legend('Obs', 'Pred', 'Location', 'best', 'FontSize', 6);
                                grid on;
                            end
                            
                            sgtitle(['Observed vs Predicted for All Curves in Model 1: ', ...
                                      '$\mathbf{y}(s)=$ ', SelectedChanLabels_2, ...
                                      ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_1, ...
                                      sprintf(' (Page %d of %d)', page, nPages)], ...
                                      'Interpreter', 'latex', 'FontSize', 13, 'FontWeight', 'bold');
                        end 


                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    
                    % Get evaluation grid from linmod output
                    tfine_2 = TimeIntervalSeleInput_BSpline_2;
                    
                    Y_obs_2 = eval_fd(ADS_EData_Bspline_1_2, tfine_2);    % [nfine x ncurves]
                    Y_pred_2 = eval_fd(FOF_linmodStr_2.yhat, tfine_2);    % [nfine x ncurves]
                    
                    sqerr_2 = (Y_obs_2 - Y_pred_2).^2;
                    RMSE_per_curve_2 = sqrt(mean(sqerr_2,1));          % per curve
                    RMSE_overall_2   = sqrt(mean(sqerr_2(:)));         % overall
                    
                    % Compute R^2 for each curve
                    varY_2 = var(Y_obs_2,0,1);
                    R2_per_curve_2 = 1 - mean(sqerr_2,1)./varY_2;
                    R2_overall_2   = 1 - mean(sqerr_2(:))/mean(varY_2);
                    
                    
                    % Residuals
                    residuals_2 = Y_obs_2 - Y_pred_2;
                    
                    % Create figure
                    figure('Name', 'Residual Diagnostics', 'NumberTitle', 'off', 'Position', [100 100 1400 800]);
                    
                    %% 1. Residuals over time (all curves)
                    subplot(2,3,1);
                    plot(tfine_2, residuals_2, 'k');
                    hold on;
                    plot(tfine_2, mean(residuals_2, 2), 'r', 'LineWidth', 2);
                    xlabel('Time');
                    ylabel('Residual');
                    title('Residuals over Time');
                    legend('Residuals (curves)', 'Mean residual', 'Location', 'best');
                    grid on;
                    
                    %% 2. Residual variance over time
                    subplot(2,3,2);
                    plot(tfine_2, var(residuals_2, 0, 2), 'b', 'LineWidth', 2);
                    xlabel('Time');
                    ylabel('Variance');
                    title('Residual Variance over Time');
                    grid on;
                    
                    %% 3. Histogram of residuals
                    subplot(2,3,3);
                    histogram(residuals_2(:), 30, 'FaceColor', [0.2 0.6 0.8]);
                    xlabel('Residual');
                    ylabel('Frequency');
                    title('Histogram of Residuals');
                    grid on;
                    
                    %% 4. Q-Q plot
                    subplot(2,3,4);
                    qqplot(residuals_2(:));
                    title('Q-Q Plot of Residuals');
                    
                    %% 5. Residuals vs Predicted values (Homoscedasticity check)
                    subplot(2,3,5);
                    scatter(Y_pred_2(:), residuals_2(:), 15, 'filled', 'MarkerFaceAlpha', 0.6);
                    xlabel('Predicted Values');
                    ylabel('Residuals');
                    title('Residuals vs Predicted');
                    refline(0,0);
                    grid on;
                    
                    %% 6. Autocorrelation of residuals
                    subplot(2,3,6);
                    My_Autocorr(residuals_2(:), 50);
                    title('Residual Autocorrelation');
                    
                sgtitle(['Residual Diagnostics for Functional Linear Model 2: ', ...
         '$\mathbf{y}(s)=$ ', SelectedChanLabels_1, ...
         ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_2], ...
         'Interpreter', 'latex', 'FontSize', 13, 'FontWeight', 'bold');



                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Each plot 
                        ncurves_2 = size(Y_obs_2, 2);
                        
                        % Fixed grid size: 5 x 5
                        nrows_2 = 5;
                        ncols_2 = 5;
                        plots_per_page = nrows_2 * ncols_2;
                        
                        % Number of pages needed
                        nPages = ceil(ncurves_2 / plots_per_page);
                        
                        for page = 1:nPages
                            figure('Name', sprintf('Observed vs Predicted (Page %d)', page), ...
                                   'NumberTitle', 'off', 'Position', [100 100 1400 800]);
                            
                            start_idx = (page-1) * plots_per_page + 1;
                            end_idx = min(page * plots_per_page, ncurves_2);
                            
                            for i = start_idx:end_idx
                                subplot(nrows_2, ncols_2, i - start_idx + 1);
                                plot(tfine_2, Y_obs_2(:,i), 'b', 'LineWidth', 1.2); hold on;
                                plot(tfine_2, Y_pred_2(:,i), 'r--', 'LineWidth', 1.2);
                                xlabel('Time');
                                ylabel('Amplitude');
                                title(sprintf('Curve %d', i), 'FontSize', 8);
                                legend('Obs', 'Pred', 'Location', 'best', 'FontSize', 6);
                                grid on;
                            end
                            
                            sgtitle(['Observed vs Predicted for All Curves in Model 2: ', ...
                                      '$\mathbf{y}(s)=$ ', SelectedChanLabels_1, ...
                                      ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_2, ...
                                      sprintf(' (Page %d of %d)', page, nPages)], ...
                                      'Interpreter', 'latex', 'FontSize', 13, 'FontWeight', 'bold');
                        end 
                                 
                        
                        % ----- Functional accuracy for Model 1 -----
                        res_1  = Y_obs_1 - Y_pred_1;                            % residuals
                        
                        % Timewise mean curve (baseline for total variation)
                        ybar_t_1 = mean(Y_obs_1, 2);                        % [nt x 1]
                        
                        % Time-integrated sums of squares (trapezoidal rule)
                        SS_res_1 = 0; 
                        SS_tot_1 = 0;
                        for i = 1:size(Y_obs_1,2)
                            SS_res_1 = SS_res_1 + trapz(tfine_1(:), (res_1(:,i)).^2);
                            SS_tot_1 = SS_tot_1 + trapz(tfine_1(:), (Y_obs_1(:,i) - ybar_t_1).^2);
                        end
                        
                        R2_func_1 = 1 - SS_res_1/SS_tot_1;
                        
                        % Integrated MSE and functional RMSE
                        T_range_1   = tfine_1(end) - tfine_1(1);                     % |T|
                        IMSE_1    = SS_res_1 / (size(Y_obs_1,2) * T_range_1);
                        RMSE_func_1 = sqrt(IMSE_1);
                        
                        
                        % ----- Functional accuracy for Model 2 -----
                        res_2  = Y_obs_2 - Y_pred_2;                            % residuals
                        
                        % Timewise mean curve (baseline for total variation)
                        ybar_t_2 = mean(Y_obs_2, 2);                        % [nt x 1]
                        
                        % Time-integrated sums of squares (trapezoidal rule)
                        SS_res_2 = 0; 
                        SS_tot_2 = 0;
                        for i = 1:size(Y_obs_2,2)
                            SS_res_2 = SS_res_2 + trapz(tfine_2(:), (res_2(:,i)).^2);
                            SS_tot_2 = SS_tot_2 + trapz(tfine_2(:), (Y_obs_2(:,i) - ybar_t_2).^2);
                        end
                        
                        R2_func_2 = 1 - SS_res_2/SS_tot_2;
                        
                        % Integrated MSE and functional RMSE
                        T_range_2   = tfine_2(end) - tfine_2(1);                     % |T|
                        IMSE_2    = SS_res_2 / (size(Y_obs_2,2) * T_range_2);
                        RMSE_func_2 = sqrt(IMSE_2);
                        
                        
                    
                % --- After computing RMSE_overall_1 and R2_overall_1 ---
                fprintf('\n===== Model 1: y(s) = %s on x(t) = %s =====\n', SelectedChanLabels_2, SelectedChanLabels_1);
                fprintf('Total RMSE: %.4f\n', RMSE_overall_1);
                fprintf('Total R^2: %.4f\n', R2_overall_1);
                fprintf('\n===== Model 1 (functional metrics) =====\n');
                fprintf('Functional RMSE: %.4f\n', RMSE_func_1);
                fprintf('Functional R^2 : %.4f\n', R2_func_1);
                fprintf('\n');
                fprintf('\n');
            
                % --- After computing RMSE_overall_2 and R2_overall_2 ---
                fprintf('\n===== Model 2: y(s) = %s on x(t) = %s =====\n', SelectedChanLabels_1, SelectedChanLabels_2);
                fprintf('Total RMSE: %.4f\n', RMSE_overall_2);
                fprintf('Total R^2: %.4f\n', R2_overall_2);
                fprintf('\n===== Model 2 (functional metrics) =====\n');
                fprintf('Functional RMSE: %.4f\n', RMSE_func_2);
                fprintf('Functional R^2 : %.4f\n', R2_func_2);
                fprintf('\n');
                fprintf('\n');




                   end 



                   if Bootstrap_Yes == 0
                  fprintf('No Bootstrap Selected. \n');


                   %%%%%%% Surface 3D Plot
                   figure;
                     % Plot the transparent surface
                        %% --- Top left: 3D surface (y1 on x2) ---
                        subplot(2,2,1)
                        hold on
                        hSurf_1_2 = surf(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, FOF_beta1fd_1);
                        set(hSurf_1_2, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
                        contour3(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, FOF_beta1fd_1, 20, 'k');
                        shading interp; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        zlabel('$\beta(s,t)$', 'Interpreter', 'latex', 'FontSize', 12);
                        title(['3D: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_2, ...
                               ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 13);
                        view(45,30); camlight left; lighting gouraud;
                        hold off
                        
                        %% --- Top right: 3D surface (y2 on x1) ---
                        subplot(2,2,2)
                        hold on
                        hSurf_2_1 = surf(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, FOF_beta1fd_2);
                        set(hSurf_2_1, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
                        contour3(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, FOF_beta1fd_2, 20, 'k');
                        shading interp; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        zlabel('$\beta(s,t)$', 'Interpreter', 'latex', 'FontSize', 12);
                        title(['3D: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_1, ...
                               ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 13);
                        view(45,30); camlight left; lighting gouraud;
                        hold off
                        
                        %% --- Bottom left: Heatmap (y1 on x2) ---
                        subplot(2,2,3)
                        imagesc(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, FOF_beta1fd_1)
                        axis xy; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        title(['Heatmap: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 13);
                        
                        %% --- Bottom right: Heatmap (y2 on x1) ---
                        subplot(2,2,4)
                        imagesc(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, FOF_beta1fd_2)
                        axis xy; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        title(['Heatmap: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 13);

                     fprintf('The Estimated Coefficinets are plotted. \n');

                      %%%%%%%%%%%%%%%%%%%% Intercept Plots
                      %% --- Row 3: Intercept curves ---
                        figure;
                        plot(FOF_linmodStr_1.alpha)
                        xlabel('Time', 'FontSize', 12)
                        ylabel('$\alpha(t)$', 'Interpreter', 'latex', 'FontSize', 12)
                        title(['Intercept $\alpha(t)$ for $\mathbf{y}(s)=$ ', SelectedChanLabels_2], ...
                            'Interpreter', 'latex', 'FontSize', 13)
                        grid on

                        figure;
                        plot(FOF_linmodStr_2.alpha)
                        xlabel('Time', 'FontSize', 12)
                        ylabel('$\alpha(t)$', 'Interpreter', 'latex', 'FontSize', 12)
                        title(['Intercept $\alpha(t)$ for $\mathbf{y}(s)=$ ', SelectedChanLabels_1], ...
                            'Interpreter', 'latex', 'FontSize', 13)
                        grid on
                        fprintf('The Estimated Intercept are plotted. \n');

                   end

                   if Bootstrap_Yes == 1 
                       fprintf('Bootstrap Selected. \n');                      
                       %% Bootstrap setup
                        B = FOF_BootStrap_N ;        % number of bootstrap resamples
                        alpha = FOF_BootStrap_Alpa ;   % 95% CI
                        CI_level = (1 - alpha) * 100;  % converts to percentage


                        fprintf('Bootstrap Selected.\nNumber of Bootstrap: %d\nAlpha Level: %.2f\n', B, alpha);


                        % Evaluate beta(s,t) on a grid
                        svals = TimeIntervalSeleInput_BSpline_1; % same s grid as you used for eval_bifd
                        tvals = TimeIntervalSeleInput_BSpline_2; % same t grid
                        beta_hat_1 = eval_bifd(FOF_linmodStr_1.beta, svals, tvals);
                        beta_hat_2 = eval_bifd(FOF_linmodStr_2.beta, svals, tvals);

                        beta_boot_1 = zeros([size(beta_hat_1), B]);
                        beta_boot_2 = zeros([size(beta_hat_2), B]);

                        
                        % Get coefficients for original data
                        coefx = getcoef(ADS_EData_Bspline_1_1);
                        coefy = getcoef(ADS_EData_Bspline_1_2);
                        ncurves = size(coefx, 2);


                        %% Bootstrap loop
                        for b = 1:B
                            % Resample curve indices
                            idx = randsample(ncurves, ncurves, true);
                            
                            % Bootstrap functional data objects
                            xfd_b = fd(coefx(:, idx), getbasis(ADS_EData_Bspline_1_1));
                            yfd_b = fd(coefy(:, idx), getbasis(ADS_EData_Bspline_1_2));
                            
                            % Fit model to bootstrap sample
                            fit_b_1 = linmod(xfd_b, yfd_b, FOF_BetaCell);
                            fit_b_2 = linmod(yfd_b, xfd_b, FOF_BetaCell);
                            
                            % Store beta(s,t) on same grid
                            beta_boot_1(:,:,b) = eval_bifd(fit_b_1.beta, svals, tvals);
                            beta_boot_2(:,:,b) = eval_bifd(fit_b_2.beta, svals, tvals);

                            if mod(b,10) == 0
                                 fprintf('Bootstrap iteration %d of %d completed.\n', b, B);
                            end
                        end
                                 fprintf('Bootstrap is Completed.\n ');


                     %% Basic bootstrap CI
                    beta_1_ci_low  = 2*beta_hat_1  - prctile(beta_boot_1, (1-alpha/2)*100, 3);
                    beta_1_ci_high = 2*beta_hat_1  - prctile(beta_boot_1, (alpha/2)*100, 3);

                    beta_2_ci_low  = 2*beta_hat_2 - prctile(beta_boot_2, (1-alpha/2)*100, 3);
                    beta_2_ci_high = 2*beta_hat_2 - prctile(beta_boot_2, (alpha/2)*100, 3);


                    
                    %% Pointwise p-values
                    pvals_1 = zeros(size(beta_hat_1));
                    for si = 1:size(beta_hat_1,1)       % loop over s
                        for ti = 1:size(beta_hat_1,2)   % loop over t
                            less0 = mean(beta_boot_1(si,ti,:) <= 0);
                            more0 = mean(beta_boot_1(si,ti,:) >= 0);
                            pvals_1(si,ti) = 2 * min(less0, more0);
                        end
                    end

                    pvals_2 = zeros(size(beta_hat_2));
                    for si = 1:size(beta_hat_2,1)       % loop over s
                        for ti = 1:size(beta_hat_2,2)   % loop over t
                            less0 = mean(beta_boot_2(si,ti,:) <= 0);
                            more0 = mean(beta_boot_2(si,ti,:) >= 0);
                            pvals_2(si,ti) = 2 * min(less0, more0);
                        end
                    end




                    %% Plot p-values heatmap
                    figure;
                    subplot(2,2,1)
                    imagesc(svals, tvals, pvals_1');
                    colorbar;
                    axis xy;
                    title('Pointwise bootstrap p-values for \beta(s,t)');
                    xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                    ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
   
                    hold on;
                    % Contours at p=0.05 and p=0.1
                    [C,h] = contour(svals, tvals, pvals_1', [0.05 0.1], 'k', 'LineWidth', 1.5);
                    clabel(C,h, 'FontSize', 10, 'Color', 'k');
                    hold off;

                    subplot(2,2,2)
                    imagesc(svals, tvals, pvals_2');
                    colorbar;
                    axis xy;
                    xlabel('s'); ylabel('t');
                    title('Pointwise bootstrap p-values for \beta(s,t)');
                    xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                    ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                    hold on;
                    % Contours at p=0.05 and p=0.1
                    [C,h] = contour(svals, tvals, pvals_2', [0.05 0.1], 'k', 'LineWidth', 1.5);
                    clabel(C,h, 'FontSize', 10, 'Color', 'k');
                    hold off;



                    %% Example: CI for middle t
                    subplot(2,2,3)
                    t_index = ceil(length(tvals)/2);
                    plot(svals, beta_hat_1(:,t_index), 'k', 'LineWidth', 2); hold on;
                    plot(svals, beta_1_ci_low(:,t_index), 'r--');
                    plot(svals, beta_1_ci_high(:,t_index), 'r--');
                    xlabel('s'); ylabel('\beta(s,t_{mid})');
                    title(sprintf('Bootstrap %.0f%% CI for Intercept (pointwise)', CI_level));
                    legend('Estimate','Lower CI','Upper CI');

                    subplot(2,2,4)
                    t_index = ceil(length(tvals)/2);
                    plot(svals, beta_hat_2(:,t_index), 'k', 'LineWidth', 2); hold on;
                    plot(svals, beta_2_ci_low(:,t_index), 'r--');
                    plot(svals, beta_2_ci_high(:,t_index), 'r--');
                    xlabel('s'); ylabel('\beta(s,t_{mid})');
                    title(sprintf('Bootstrap %.0f%% CI for Intercept (pointwise)', CI_level));
                    legend('Estimate','Lower CI','Upper CI');

                     fprintf('Plot I: Pointwise P-Value for Coefficients are done ! \n\n\n');


                     figure;
                     % Plot the transparent surface
                        %% --- Top left: 3D surface (y1 on x2) ---
                        subplot(2,2,1)
                        hold on
                        hSurf_1_2 = surf(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, FOF_beta1fd_1);
                        set(hSurf_1_2, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
                        contour3(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, FOF_beta1fd_1, 20, 'k');
                        shading interp; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        zlabel('$\beta(s,t)$', 'Interpreter', 'latex', 'FontSize', 12);
                        title(['3D: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_2, ...
                               ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 13);
                        view(45,30); camlight left; lighting gouraud;
                        hold off
                        
                        %% --- Top right: 3D surface (y2 on x1) ---
                        subplot(2,2,2)
                        hold on
                        hSurf_2_1 = surf(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, FOF_beta1fd_2);
                        set(hSurf_2_1, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
                        contour3(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, FOF_beta1fd_2, 20, 'k');
                        shading interp; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        zlabel('$\beta(s,t)$', 'Interpreter', 'latex', 'FontSize', 12);
                        title(['3D: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_1, ...
                               ' on $\mathbf{x}(t)=$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 13);
                        view(45,30); camlight left; lighting gouraud;
                        hold off
                        
                        %% --- Bottom left: Heatmap (y1 on x2) ---
                        subplot(2,2,3)
                        imagesc(TimeIntervalSeleInput_BSpline_1, TimeIntervalSeleInput_BSpline_2, FOF_beta1fd_1)
                        axis xy; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        title(['Heatmap: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 13);
                        subtitle('P-Value Contour Map', 'Interpreter', 'latex', 'FontSize', 13);
                        hold on;
                        % Contours at p=0.05 and p=0.1
                        [C,h] = contour(svals, tvals, pvals_1', [0.05 0.1], 'k', 'LineWidth', 1.5);
                        clabel(C,h, 'FontSize', 10, 'Color', 'k');
                        hold off;


                        %% --- Bottom right: Heatmap (y2 on x1) ---
                        subplot(2,2,4)
                        imagesc(TimeIntervalSeleInput_BSpline_2, TimeIntervalSeleInput_BSpline_1, FOF_beta1fd_2)
                        axis xy; colormap(jet); colorbar;
                        xlabel(['$\mathbf{x}(t) =$ ', SelectedChanLabels_2], 'Interpreter', 'latex', 'FontSize', 12);
                        ylabel(['$\mathbf{y}(s) =$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 12);
                        title(['Heatmap: $\beta(s,t)$, $\mathbf{y}(s)=$ ', SelectedChanLabels_1], 'Interpreter', 'latex', 'FontSize', 13);
                        subtitle('P-Value Contour Map', 'Interpreter', 'latex', 'FontSize', 13);
                        hold on;
                        % Contours at p=0.05 and p=0.1
                        [C,h] = contour(svals, tvals, pvals_2', [0.05 0.1], 'k', 'LineWidth', 1.5);
                        clabel(C,h, 'FontSize', 10, 'Color', 'k');
                        hold off;
                         fprintf('Plot II: Estimated Coefficients and P-Values Contour Map are done ! \n');



                  end
                   
end 