classdef RegressionKernel < ...
        classreg.learning.regr.RegressionModel & classreg.learning.Kernel
%RegressionKernel Kernel model for regression.
%   RegressionKernel is a kernel model for regression for which the
%   predictors are explicitly mapped to a high dimensional space.
%
%   RegressionKernel properties:
%       BoxConstraint          - Box constraint for SVM.
%       Epsilon                - Half of the width of the epsilon-insensitive band for SVM.
%       ExpandedPredictorNames - Names of expanded predictors.
%       NumExpansionDimensions - Number of dimensions of expanded space.
%       FittedLoss             - Fitted loss function, 'epsiloninsensitive' or 'mse'.
%       KernelScale            - Kernel scale.
%       Lambda                 - Regularization strength.
%       Learner                - Learned model, 'svm' or 'leastsquares'.
%       PredictorNames         - Names of predictors used for this model.
%       ResponseName           - Name of the response variable.
%       ResponseTransform      - Transformation applied to predicted response values.
%       Mu                     - Predictor means.
%       Sigma                  - Predictor standard deviations.
%
%   RegressionKernel methods:
%       loss                  - Regression loss.
%       incrementalLearner    - Return an incremental regression kernel model.
%       partialDependence     - Partial Dependence values.
%       plotPartialDependence - Plot Partial Dependence.
%       predict               - Predicted response of this model.
%       resume                - Continue fitting this model.

%   Copyright 2017-2023 The MathWorks, Inc.

    properties(GetAccess=protected,SetAccess=protected,Hidden=true)
        FitInfo = [];
    end

    properties(GetAccess=public,SetAccess=protected)
        %Epsilon Epsilon-insensitive band
        %   The Epsilon property is a non negative scalar specifying half
        %   the width of the epsilon-insensitive band. Epsilon applies to
        %   SVM learners only.
        %
        %   See also RegressionKernel.
        Epsilon = [];

        %BoxConstraint Regularization strength.
        %   The BoxConstraint property is a non negative scalar used to
        %   control the strength of regularization. BoxConstraint applies
        %   to SVM learners only. BoxConstraint and Lambda are related by:
        %
        %               BoxConstraint * Lambda = 1 / N
        %
        %   where N is the number of observation in the training set.
        %
        %   See also ClassificationKernel, RegressionKernel.
        BoxConstraint = [];
    end

    methods(Hidden)
        function this = RegressionKernel(dataSummary,responseTransform,defaultResponse)
        % Protect against being called by the user
            if ~isstruct(dataSummary)
                error(message('stats:RegressionKernel:DoNotUseConstructor'));
            end

            this = this@classreg.learning.regr.RegressionModel(...
                dataSummary,responseTransform,defaultResponse);
            this = this@classreg.learning.Kernel;
            this.CategoricalVariableCoding = 'dummy';
            this.AllowedPredictArgs = struct('Name','PredictionForMissingValue','Value','Median');
        end

        function cmp = compact(this)
            cmp = this;
        end

        function varargout = gather(varargin)
            % Overload gather@statslib.internal.GatherableMixin to make
            % method hidden.
            [varargout{1:nargout}] = gather@statslib.internal.GatherableMixin(varargin{:});
        end
    end

    methods(Access=protected)
        function s = propsForDisp(this,s)
            s = propsForDisp@classreg.learning.regr.RegressionModel(this,s);
            s = propsForDisp@classreg.learning.Kernel(this,s);
            if strcmp(this.Learner,'svm')
                s.BoxConstraint  = this.BoxConstraint;
                s.Epsilon        = this.Epsilon;
            end
            s = rmfield(s,'CategoricalPredictors');
            if strcmp(s.ResponseTransform,'none')
                s = rmfield(s,'ResponseTransform');
            end
        end

        function S = response(this,X,varargin)

        % Implements abstract RegressionModel/response, assumes
        % observations in columns.

            [varargin{:}] = convertStringsToChars(varargin{:});

            args = {'predictionformissingvalue'};
            defs = {                   'median'};
            predMissing = internal.stats.parseArgs(args,defs,varargin{:});

            predMissing = classreg.learning.internal.validatePredictionForMissing(...
                predMissing,mfilename);

            internal.stats.checkSupportedNumeric('X',X)

            X = X'; %Because this.ObservationsInRows is false

            if this.ModelParameters.StandardizeData
                X = classreg.learning.Kernel.standardizePredict(X,this.Mu,this.Sigma);
            end
            
            [n,p] = size(X);
            if p~=this.FeatureMapper.d
                error(message('stats:RegressionKernel:XSizeMismatch',this.FeatureMapper.d))
            end

            maxChunkSize = RegressionKernel.estimateMaxChunkSize(this.FeatureMapper.n,this.ModelParams.BlockSize,false);
            numberChunks = ceil(n / maxChunkSize);

            if numberChunks<=1
                Xm = map(this.FeatureMapper,X,this.KernelScale);
                S = score(this.Impl,Xm,false,true); % doClass = false , obsInRows = true
            else
                % Block-wise
                S = cell(numberChunks,1);
                j = 1;
                for i = 1:numberChunks
                    k = min(n,j+maxChunkSize-1);
                    Xm = map(this.FeatureMapper,X(j:k,:),this.KernelScale);
                    S{i} = score(this.Impl,Xm,false,true); % doClass = false , obsInRows = true
                    j = j + maxChunkSize;
                end
                S = cell2mat(S);
            end

            nanX = any(isnan(X),2);
            if any(nanX)
                S(nanX,:) = classreg.learning.internal.predictionForMissingValue(...
                    predMissing,this.DefaultResponse);
            end
        end

        function n = getExpandedPredictorNames(this)
            n = classreg.learning.internal.expandPredictorNamesWithReference(...
                this.PredictorNames,this.CategoricalPredictors,...
                this.OrdinalPredictors,this.VariableRange,...
                this.ReferenceLevels);
        end
    end

    methods
        function  [model,fitinfo] = resume(this,X,Y,varargin)
        %RESUME Resume training of the model.
        %   MODEL_OUT = RESUME(MODEL,X,Y) continues training the same model
        %   starting at the current estimated parameters. X must be a table if
        %   MODEL was originally trained on a table and a numeric matrix if MODEL
        %   was originally trained on a matrix. If X is a table, it must contain
        %   all the predictors used for training the model. If X is a matrix, it
        %   must have size N-by-P, where P is the number of predictors used for
        %   training. Y must be a vector of N floating-point numbers, where N is
        %   the number of observations (rows) in X. You can omit Y if it appears in
        %   the table X. To effectively continue the training, X and Y should be
        %   the same arrays used to initially fit the model MODEL.
        %
        %   [MODEL_OUT,FITINFO] = RESUME(...) also returns FITINFO, a struct
        %   containing fit information.
        %
        %   [...] = RESUME(MODEL,X,Y,'PARAM1',val1,'PARAM2',val2,...) specifies
        %   optional name-value parameter pairs:
        %
        %            'Weights'  - Vector of observation weights, one weight per
        %                         observation. RESUME normalizes the weights to
        %                         add up to one. Default is equal weights.
        %      'BetaTolerance'  - Relative tolerance on linear coefficients and
        %                         bias term. A non-negative scalar. Default
        %                         is same BetaTolerance value used to train
        %                         MODEL.
        %   'GradientTolerance' - Absolute gradient tolerance. A non-negative
        %                         scalar. Default is the same GradientTolerance
        %                         value used to train MODEL.
        %     'IterationLimit'  - Maximal number of extra optimization iterations.
        %                         A positive integer. Default is 1000 if the
        %                         transformed data fits in memory, 100 if the
        %                         fitting algorithm switches to a block-wise
        %                         strategy, or 20 if inputs are tall.
        %
        %   See also fitrkernel, RegressionKernel.

            [varargin{:}] = convertStringsToChars(varargin{:});
            args = {'weights'  'betatolerance'                    'gradienttolerance'                    'iterationlimit'};
            defs = {[]         this.ModelParameters.BetaTolerance this.ModelParameters.GradientTolerance []};
            [W,betaTolerance,gradientTolerance,iterationLimit] = internal.stats.parseArgs(args,defs,varargin{:});
            temp = RegressionKernel.template('Weights',W,...
                                             'BetaTolerance',betaTolerance,...
                                             'GradientTolerance',gradientTolerance,...
                                             'Epsilon',this.Epsilon,...  % Taken from object as ModelParams may have 'auto'
                                             'FitBias',this.ModelParameters.FitBias,...
                                             'HessianHistorySize',this.ModelParameters.HessianHistorySize,...
                                             'Learner',this.ModelParameters.Learner,...
                                             'Lambda',this.Lambda,... % Taken from object as ModelParams may have 'auto'
                                             'LineSearch',this.ModelParameters.LineSearch,...
                                             'LossFunction',this.ModelParameters.LossFunction,...
                                             'Regularization',this.ModelParameters.Regularization,...
                                             'Solver',this.ModelParameters.Solver,...
                                             'RandomStream',this.ModelParameters.Stream,...
                                             'Verbose',this.ModelParameters.VerbosityLevel,...
                                             'NumExpansionDimensions',this.NumExpansionDimensions,... % Taken from object as ModelParams may have 'auto'
                                             'KernelScale',this.KernelScale,... % Taken from object as ModelParams may have 'auto'
                                             'Transformation',this.ModelParameters.Transformation,...
                                             'BlockSize',this.ModelParameters.BlockSize,...
                                             'ADMMIterationLimit',0,... % No initialization needed when resuming
                                             'Consensus',this.ModelParameters.Consensus,...
                                             'InitialStepSize',this.ModelParameters.InitialStepSize,...
                                             'Standardize',this.ModelParameters.StandardizeData);

            % IterationLimit is handled differently as two parameters
            % depend on it (IterationLimit and IterationLimitBlockWise),
            % if these parameters were defaulted in the original fit, they have
            % *different* values and we need to set them individually in the new
            % template. We cannot set them with the template constructor (as all
            % other parameters) because they would be set to the same value.
            if ~isempty(iterationLimit)
                tempDummy = RegressionKernel.template('iterationLimit',iterationLimit);
                temp.ModelParams.IterationLimit = tempDummy.ModelParams.IterationLimit;
                temp.ModelParams.IterationLimitBlockWise = tempDummy.ModelParams.IterationLimitBlockWise;
            end

            %Starting point
            temp.ModelParams.InitialBeta = this.Beta;
            temp.ModelParams.InitialBias = this.Bias;
            temp.ModelParams.FeatureMapper = this.FeatureMapper;


            if istall(X) || istall(Y) || istall(W)
                if getNumberPredictors(this) ~= gather(size(X,2))
                    error(message('stats:RegressionKernel:XSizeMismatch',this.FeatureMapper.d))
                end
                [model,fitinfo] = fitrkernel(X,Y,'Template',temp);
            elseif ~istable(X)
                if getNumberPredictors(this) ~= size(X,2)
                    error(message('stats:RegressionKernel:XSizeMismatch',this.FeatureMapper.d))
                end
                if any(this.CategoricalPredictors)
                    % Check that the input data encodes similarly to model,
                    % otherwise fitted coefficients will not match.
                    D =getNumberPredictors(this);
                    foundvrange = cell(1,D);
                    iscat = false(1,D);
                    iscat(this.CategoricalPredictors) = true;
                    for k=1:D
                        if iscat(k)
                            x = X(:,k);
                            vrk = unique(x);
                            foundvrange{k} = vrk(~isnan(vrk));
                        end
                    end
                    if ~isequal(foundvrange,this.VariableRange)
                        error(message('stats:RegressionKernel:CategoricalEncoding'))
                    end
                end
                temp = setBaseArg(temp,'CategoricalPredictors',this.CategoricalPredictors);
                [model,fitinfo] = fit(temp,X,Y);
            else
                if internal.stats.isString(Y)
                    Y = this.ResponseName;
                end

                % Check the input table has the required Predictors, also
                % checks their type is appropriate as it is trying to
                % encode as prescribed by the model
                [~,Yt] = classreg.learning.internal.table2PredictMatrix(X,Y,W,...
                                                                           this.VariableRange,this.CategoricalPredictors,...
                                                                           this.PredictorNames,true); % ObsInRows is true

                % Filter potential extra predictors in data before encoding
                Xt = X(:,this.PredictorNames);
                if internal.stats.isString(Y)
                    Xt.(this.ResponseName) = Yt;
                end

                % Check that the input data encodes similarly to model,
                % otherwise fitted coefficients will not match.
                [~,~,foundvrange] = classreg.learning.internal.table2FitMatrix(Xt,Y,...
                                                                               'CategoricalPredictors',this.CategoricalPredictors);

                if ~isequal(foundvrange,this.VariableRange)
                    error(message('stats:RegressionKernel:CategoricalEncoding'))
                end

                temp = setBaseArg(temp,'CategoricalPredictors',this.CategoricalPredictors);
                temp = setBaseArg(temp,'PredictorNames',this.PredictorNames);
                temp = setBaseArg(temp,'OrdinalPredictors',this.OrdinalPredictors);
                temp = setBaseArg(temp,'ResponseName',this.ResponseName);

                [model,fitinfo] = fit(temp,Xt,Y);
            end

            % Check if it is needed to concatenate the History of a
            % previous fit:
            if ~isempty(this.FitInfo.History)
                if isempty(fitinfo.History)
                    % Current fit did not return History, so only the
                    % previous History is copied:
                    fitinfo.History = this.FitInfo.History;
                else
                    % Concatenate both History structures
                    fns = fieldnames(fitinfo.History);
                    for i = 1:numel(fns)
                        fitinfo.History.(fns{i}) = [this.FitInfo.History.(fns{i});fitinfo.History.(fns{i})];
                    end
                end
                model.FitInfo.History = fitinfo.History;
            end

        end

        function [varargout] = predict(this,X,varargin)
        %PREDICT Predict responses of the model.
        %   YFIT=PREDICT(MODEL,X) returns predicted responses for the regression
        %   kernel model MODEL and  predictors X. X must be a table if MODEL was
        %   originally trained on a table and a numeric matrix if MODEL was
        %   originally trained on a matrix. If X is a table, it must contain all
        %   the predictors used for training the model.  If X is a matrix, it must
        %   have size N-by-P, where P is the number of predictors used for training.
        %   YFIT is a N-by-1 array with the same type data as the response data
        %   used for training.
        %
        %   YFIT = PREDICT(MODEL,X,'NAME1',VALUE1,'NAME2',VALUE2,...)
        %   specifies optional name/value arguments:
        %       'PredictionForMissingValue' - Numeric scalar, character vector, or string scalar. 
        %                                     If a character vector or string scalar, must be one of: 
        %                                     'median' or 'mean'. If a numeric scalar, PREDICT uses
        %                                     this response for observations with missing predictor 
        %                                     values. If 'mean' or 'median', PREDICT uses the mean 
        %                                     or median of observed response values in the training data. 
        %                                     Default: 'median'
        %
        %   See also RegressionKernel.
            [varargout{1:nargout}] = predict@classreg.learning.regr.RegressionModel(this,X,varargin{:});
        end

        function l = loss(this,X,varargin)
        %LOSS Regression error.
        %   ERR=LOSS(MODEL,X,Y) returns the mean squared error for the model
        %   MODEL using predictors X and observed response Y. X must be a table if
        %   MODEL was originally trained on a table and a numeric matrix if MODEL
        %   was originally trained on a matrix. If X is a table, it must contain
        %   all the predictors used for training the model. If X is a matrix, it
        %   must have size N-by-P, where P is the number of predictors used for
        %   training. Y must be a vector of N floating-point numbers, where N is
        %   the number of observations (rows) in X. You can omit Y if it appears in
        %   the table X.
        %
        %   ERR=LOSS(MODEL,X,Y,'PARAM1',val1,'PARAM2',val2,...) specifies optional
        %   name-value parameter pairs:
        %       'LossFun'          - Function handle for loss, or string
        %                            representing a built-in loss function.
        %                            Available loss functions for regression: 'mse'
        %                            and 'epsiloninsensitive'. If you pass a
        %                            function handle FUN, LOSS calls it as shown
        %                            below:
        %                               FUN(Y,Yfit,W)
        %                            where Y, Yfit and W are numeric vectors of
        %                            length N. Y is the observed response, Yfit is
        %                            the predicted response, and W is the
        %                            observation weights. Default: 'mse'
        %       'Weights'          - Vector of observation weights. The length of
        %                            this vector must be equal to the number of
        %                            rows in X. Default is equal weights.
        %       'PredictionForMissingValue' - Numeric scalar, character vector, or
        %                                     string scalar. If a character vector or
        %                                     string scalar, must be one of:
        %                                     'median', 'mean', or 'omitted'. If a
        %                                     numeric scalar, LOSS uses this response 
        %                                     for observations with missing predictor 
        %                                     values. If 'mean' or 'median', LOSS uses 
        %                                     the mean or median of observed response 
        %                                     values in the training data. If 'omitted',
        %                                     LOSS excludes observations with missing 
        %                                     predictor values from the computation. 
        %                                     Default: 'median'
        %
        %   See also RegressionKernel, predict.

        % Overrides RegressionModel/loss because tables are not
        % supported and because it adds 'epsiloninsensitive' loss
            [varargin{:}] = convertStringsToChars(varargin{:});

            % Handle input data such as "tall" requiring a special adapter
            adapter = classreg.learning.internal.makeClassificationModelAdapter(this,X,varargin{:});
            if ~isempty(adapter)
                l = loss(adapter,X,varargin{:});
                return
            end

            [Y,varargin] = classreg.learning.internal.inferResponse(this.ResponseName,X,varargin{:});

            obsInRows = true;
            N = size(X,1);

            args = {                  'lossfun'  'weights' 'predictionformissingvalue'};
            defs = {@classreg.learning.loss.mse  ones(N,1)                    'median'};
            [funloss,W,predMissing,~,extraArgs] = internal.stats.parseArgs(args,defs,varargin{:});

            % Validate predMissing
            predMissing = classreg.learning.internal.validatePredictionForMissing(predMissing,mfilename,'loss');
            if strcmp(predMissing,'omitted')
                cleanNaNsInX = true;
            else
                cleanNaNsInX = false;
                extraArgs = [extraArgs {'predictionformissingvalue' predMissing}];
            end

            % Prepare data
            [X,Y,W] = prepareDataForLoss(this,X,Y,W,this.VariableRange,false,obsInRows,cleanNaNsInX);
            internal.stats.checkSupportedNumeric('Weights',W,true);

            % Adjust loss function in case of 'epsiloninsensitive'
            if strncmpi(funloss,'epsiloninsensitive',length(funloss))
                if isempty(this.Epsilon) % this is empty for models other than SVM
                    error(message('stats:RegressionKernel:UseEpsilonInsensitiveForSVM'));
                end
                funloss = @(Y,Yfit,W) classreg.learning.loss.epsiloninsensitive(...
                    Y,Yfit,W,this.Epsilon);
            end
            funloss = classreg.learning.internal.lossCheck(funloss,'regression');

            % Get predictions
            Yfit = predict(this,X,extraArgs{:});

            % Check
            classreg.learning.internal.regrCheck(Y,Yfit(:,1),W);

            % Get loss
            l = funloss(Y,Yfit,W);
        end

        function incrmdl = incrementalLearner(this,varargin)
        %INCREMENTALLEARNER Convert regression kernel model to incremental learner.
        %   INCREMENTALMODEL=INCREMENTALLEARNER(MODEL) returns an
        %   incremental regression kernel model INCREMENTALMODEL initialized
        %   using the parameters provided in MODEL.
        %
        %   INCREMENTALLEARNER initializes these INCREMENTALMODEL properties by
        %   using the properties of MODEL:
        %
        %       'Epsilon'                 - If MODEL.Learner is 'svm',
        %                                   MODEL.Epsilon is preserved.
        %
        %       'FitBias'
        %       'Learner'
        %       'KernelScale'
        %       'Mu'                      - Mu is copied to INCREMENTALMODEL if Standardize is true.
        %       'Sigma'                   - Sigma is copied to INCREMENTALMODEL if Standardize is true.
        %       'NumExpansionDimensions'
        %       'RandomStream'
        %
        %   INCREMENTALMODEL=INCREMENTALLEARNER(MODEL,'PARAM1',val1,'PARAM2',val2,...)
        %   specifies optional name/value arguments:
        %       'BatchSize'               - Mini-batch size for the stochastic
        %                                   solvers, specified as a positive
        %                                   scalar.
        %                                   Default is 10 for stochastic solvers.
        %                                   'BatchSize' is not supported
        %                                   for scale-invariant solver.
        %       'EstimationPeriod'        - Nonnegative integer specifying the number of
        %                                   observations that are used for
        %                                   estimation of hyperparameters prior to
        %                                   training. For a positive integer, if
        %                                   model has unknown hyperparameters
        %                                   that require estimation,
        %                                   hyperparameters are estimated.
        %                                   Default value is 1000 when estimation is needed.
        %                                   Otherwise, default value is 0.
        %                                   Refer to documentation for more information.
        %       'Lambda'                  - Ridge (L2) regularization strength specified as a
        %                                   nonnegative scalar.
        %                                   'Lambda' is not supported
        %                                   for scale-invariant solver.
        %                                   Default: 1e-5 for stochastic solvers.
        %       'LearnRate'               - Learning rate for the stochastic
        %                                   solvers, specified either as 'auto' or as a positive
        %                                   scalar. If 'Solver' is specified
        %                                   as 'sgd' or 'asgd', default is 'auto'.
        %                                   'LearnRate' is not supported for scale-invariant solver.
        %                                   Refer to documentation for info
        %                                   on 'auto' behavior.
        %       'LearnRateSchedule'       - String, one of 'decaying' or 'constant'.
        %                                   For stochastic solvers, if 'decaying',
        %                                   learning rate decreases over
        %                                   iterations. If 'constant', learning
        %                                   rate stays constant. Default is 'decaying'.
        %                                   'LearnRateSchedule' is not
        %                                   supported for scale-invariant solver.
        %       'Metrics'                 - String vector representing built-in
        %                                   metrics, cell array of function handles
        %                                   or a struct of function handles specifying
        %                                   metrics that are tracked during incremental
        %                                   training. 'Metrics' can be a
        %                                   combination of strings, function
        %                                   handles or struct of function handles
        %                                   passed as a cell array.
        %                                   Built-in metrics for regression:
        %                                       - 'epsiloninsensitive'
        %                                       - 'mse'
        %                                   'epsiloninsensitive' can be selected
        %                                   only when MODEL.Learner is 'svm'.
        %                                   If you pass a function handle FUN, it
        %                                   is called as shown below:
        %                                       FUN(Y,Yfit)
        %                                   where Y and Yfit are numeric vectors of
        %                                   length N. Y is observed response, Yfit is
        %                                   predicted response.
        %                                   Metric name will be extracted from provided function
        %                                   handle if it's not an anonymous
        %                                   function handle. If it is an anonymous
        %                                   function handle, metric will be named
        %                                   as 'CustomMetric' and it will be post-fixed
        %                                   with an underscore and an integer
        %                                   identifying its relative position in
        %                                   'Metrics' property.
        %                                   If you pass a struct, fieldnames will
        %                                   be treated as the metric names and
        %                                   values must be function handles.
        %                                   Incremental models track two values for
        %                                   each metric; a 'Cumulative' value and a
        %                                   'Window' value. 'MetricsWindowSize' determines
        %                                   the update frequency of
        %                                   'Window' metrics.
        %                                   Default is 'mse' if 'Learner' is
        %                                   'leastsquares'. Otherwise, default is
        %                                   'epsiloninsensitive'.
        %       'MetricsWarmupPeriod'     - Nonnegative integer specifying the number of
        %                                   observations that will be used for
        %                                   training the model prior to metrics
        %                                   update. Default is 0.
        %       'MetricsWindowSize'       - Positive integer specifying the window
        %                                   size that will be used for metrics
        %                                   update. Default value is 200.
        %       'Shuffle'                 - Logical value specifying
        %                                   whether observations
        %                                   are shuffled at each step.
        %                                   Default: true.
        %                                   'Shuffle' is not supported for
        %                                   stochastic solvers.
        %                                   Observations are always
        %                                   shuffled for stochastic
        %                                   solvers.
        %       'Solver'                  - String, one of: 'sgd', 'asgd' or 'scale-invariant'.
        %                                   If 'sgd', stochastic gradient descent
        %                                   solver is used. If 'asgd', average
        %                                   stochastic gradient descent solver is
        %                                   used. If 'scale-invariant', adaptive
        %                                   scale-invariant solver is used.
        %                                   Default is 'scale-invariant'.
        %
        %   Example: Initialize an incremental kernel model using a traditionally trained model
        %   and continue incremental training.
        %
        %       ttModel = fitrkernel(Xoffline,Yoffline);
        %       mdl = incrementalLearner(ttModel);
        %       mdl = updateMetricsAndFit(mdl,Xnew,Ynew,'Weights',Wnew);
        %
        %   See also incrementalRegressionKernel.
            [args,additionalargs] = incremental.utils.offlineLinearParameters(this,false,varargin{:});
            [args,additionalargs] = incremental.utils.offlineRegressionParameters(this,args,additionalargs);
            mergedargs = incremental.utils.mergeOfflineParameters(args,additionalargs);
            % set parameters exclusive to kernel models : KernelScale,Mu,Sigma
            % RandomStream, and NumExpansionDimensions
            mergedargs.KernelScale = this.KernelScale;
            mergedargs.Mu = this.Mu;
            mergedargs.Sigma = this.Sigma;
            mergedargs.NumExpansionDimensions = this.NumExpansionDimensions;
            mergedargs.RandomStream = this.ModelParameters.Stream; % user passed stream
            mergedargs.NumPredictors = numel(this.PredictorNames); % explicitly pass NumPredictors because they are not the same size as Beta
            if this.Learner=="svm"
                mergedargs.Epsilon = this.Epsilon;
            end
            incrmdl = incrementalRegressionKernel(mergedargs);
            incrmdl.Impl  = setFeatureMapper(incrmdl.Impl,this.FeatureMapper); % use the same feature mapper
        end
    end

    methods(Static,Hidden)

        function [lambda,bc] = resolveAutoLambdaOrEmptyBoxConstraint(lambda,bc,n)
            if strcmp(lambda,'auto')
                lambda = 1/bc/n;
            elseif isempty(bc)
                bc = 1/lambda/n;
            end
        end

        function numExpansionDimensions = resolveAutoNumExpansionDimensions(numExpansionDimensions,d)
            if strcmp(numExpansionDimensions,'auto')
                numExpansionDimensions = 2.^ceil(min(log2(d)+5,15));
            end
        end

        function epsilon = resolveAutoEpsilon(epsilon,lossFunction,Y)
            if strcmp(epsilon,'auto')
                switch lossFunction
                  case 'epsiloninsensitive'
                      iqrange = iqr(Y);
                      if iqrange == 0
                          epsilon = 0.1;
                      else
                          epsilon = iqrange./13.49;
                      end
                  case 'mse'
                    epsilon = [];
                end
            end
        end

        function kernelScale = resolveAutoKernelScale(kernelScale,X,Y,type)
            if strcmp(kernelScale,'auto')
                kernelScale = classreg.learning.svmutils.optimalKernelScale(X,Y,type);
            end
        end

        function featureMapper = resolveEmptyFeatureMapper(featureMapper,d,expDim,transformation,stream)
            if isempty(featureMapper)
                featureMapper = classreg.learning.rkeutils.featureMapper(stream,d,expDim,transformation);
            end
        end

        function initialBeta = resolveEmptyInitialBeta(initialBeta,expDim)
            if isempty(initialBeta)
                initialBeta = zeros(expDim,1);
            end
        end

        function initialBias = resolveEmptyInitialBias(initialBias,lossFunction,Y,W,~)
            if isempty(initialBias)
                switch lossFunction
                  case 'epsiloninsensitive'
                    % Epsilon is not used for initial Bias to avoid cases that could
                    % lead to zero objective gradient...
                    initialBias =  classreg.learning.linearutils.fitbias(...
                        'epsiloninsensitive',Y,zeros(numel(Y),1,'like',Y),W,0);
                  case 'mse'
                    initialBias = sum(W'*Y)/sum(W);
                end
            end
        end

        function maxChunkSize = estimateMaxChunkSize(numExpansionDimensions,blockSize,isFitting)
            memoryPerExpandedRow = numExpansionDimensions * 8;
            memoryBlock = blockSize * 1e6;
            maxChunkSize = floor(memoryBlock ./ memoryPerExpandedRow); %RowsPerBlock (Maximum)
            if isFitting && maxChunkSize<1000 %Only check for few rows per chunk when we are fitting
                error(message('stats:RegressionKernel:FewObservationsPerBlock'))
            end
        end

        function temp = template(varargin)
            classreg.learning.FitTemplate.catchType(varargin{:}); %You are not allowed to pass 'type' argument.
            temp = classreg.learning.FitTemplate.make('Kernel',...
                                                      'type','regression',varargin{:});
        end

        function [varargout] = fit(X,Y,varargin)
            temp = RegressionKernel.template(varargin{:});
            [varargout{1:nargout}] = fit(temp,X,Y);
        end

        function [X,Y,W,dataSummary,responseTransform] = ...
                prepareData(X,Y,varargin)
        % Reuse RegressionLinear.prepareData
            [X,Y,W,dataSummary,responseTransform] = ...
                RegressionLinear.prepareData(X,Y,varargin{:});
            % ... but we also restrict sparse data
            internal.stats.checkSupportedNumeric('X',X,false,false,false)
        end

        function obj = makeNoFit(beta,bias,lambda,featureMapper,...
                kernelScale,boxConstraint,epsilon,modelParams,dataSummary,...
                fitinfo,responseTransform,Mu,Sigma,defaultResponse)
        % Makes a RegressionKernel object without fitting by
        % feeding model parameters (beta and bias). input model
        % parameters (modelParams) and featureMapper

            if isempty(responseTransform)
                responseTransform = @classreg.learning.transform.identity;
            end

            obj = RegressionKernel(dataSummary,responseTransform,defaultResponse);
            obj.DefaultLoss = @classreg.learning.loss.mse;

            switch lower(modelParams.LossFunction)
              case 'epsiloninsensitive'
                obj.Learner = 'svm';
              case 'mse'
                obj.Learner = 'leastsquares';
            end

            obj.Epsilon = epsilon;
            obj.KernelScale = kernelScale;
            obj.FeatureMapper = featureMapper;
            obj.BoxConstraint = boxConstraint;

            modelParams.Lambda = lambda; % Impl keeps the real lambda in case modelParams is still 'auto'
            obj.Impl = classreg.learning.impl.LinearImpl.makeNoFit(modelParams,beta,bias,fitinfo);

            obj.ModelParams = modelParams;

            % save fitInfo into a private property, used only by resume method to
            % concatenate history information. User should only obtain
            % fitInfo from the second output argument.
            obj.FitInfo = fitinfo;
            obj.Mu=Mu;
            obj.Sigma=Sigma;
        end

        function [obj,fitInfo] = fitRegressionKernel(...
            X,Y,W,modelParams,dataSummary,responseTransform)

            defaultResponse.Median = classreg.learning.internal.wnanmedian(Y,W);
            defaultResponse.Mean   = classreg.learning.internal.wnanmean(Y,W);

            [X,Y,W,dataSummary] = RegressionLinear.removeNaNsFromX(X,Y,W,dataSummary);

            modelParams = fillIfNeeded(modelParams,X,Y,W,dataSummary,[]);
            
            obj = RegressionKernel(dataSummary,responseTransform,defaultResponse);

            switch lower(modelParams.LossFunction)
              case 'epsiloninsensitive'
                obj.Learner = 'svm';
              case 'mse'
                obj.Learner = 'leastsquares';
            end

            % Expand categorical predictors if necessary
            if any(obj.CategoricalPredictors)
                [X,catcols] = classreg.learning.internal.expandCategorical(X,...
                                                                 obj.CategoricalPredictors,obj.OrdinalPredictors,...
                                                                 obj.VariableRange,obj.ReferenceLevels,false);
                if isempty(X)                
                    error(message('stats:classreg:learning:impl:SVMImpl:make:OneOrdinalOneLvl'));                
                end
            else
                catcols = false(1,size(X,1));
            end
            
            if modelParams.StandardizeData
                [X,obj.Mu,obj.Sigma] = classreg.learning.Kernel.standardize(X,W,catcols);
            end

            [d,n] = size(X); % Number of predictors and observations
            if d>=4096  %2^12
                        % fitrkernel shall not be used in high dimensional datasets
                error(message('stats:RegressionKernel:InvalidInputDimensions'))
            end

            % Resolve parameters that may still be set to 'auto' or empty
            % (Order matters!). All these modifications are not placed back
            % into modelParams; the actual lambda is saved in obj.impl after
            % fitting, and the actual numExpansionDimensions is saved in
            % obj.FeatureMapper after fitting.
            [lambda,obj.BoxConstraint] = RegressionKernel.resolveAutoLambdaOrEmptyBoxConstraint(modelParams.Lambda,modelParams.BoxConstraint,n);
            numExpansionDimensions =  RegressionKernel.resolveAutoNumExpansionDimensions(modelParams.NumExpansionDimensions,d);
            if numExpansionDimensions<= d
                error(message('stats:RegressionKernel:InvalidNumExpansionDimensions'))
            end
            obj.FeatureMapper = RegressionKernel.resolveEmptyFeatureMapper(modelParams.FeatureMapper,d,numExpansionDimensions,modelParams.Transformation,modelParams.Stream);
            obj.KernelScale = RegressionKernel.resolveAutoKernelScale(modelParams.KernelScale,X',Y,0);
            obj.Epsilon = RegressionKernel.resolveAutoEpsilon(modelParams.Epsilon,modelParams.LossFunction,Y);


            initialBeta = RegressionKernel.resolveEmptyInitialBeta(modelParams.InitialBeta,numExpansionDimensions);
            initialBias = RegressionKernel.resolveEmptyInitialBias(modelParams.InitialBias,modelParams.LossFunction,Y,W,obj.Epsilon);

            % Figure out maxChunkSize from BlockSize
            maxChunkSize = RegressionKernel.estimateMaxChunkSize(numExpansionDimensions,modelParams.BlockSize,true);
            numberChunks = ceil(n / maxChunkSize);

            if numberChunks == 1
                %%%%%%%% FIT THE DATA IN ONE MEMORY BLOCK %%%%%%%%%%%%%%
                if modelParams.IterationLimit == 0
                    error(message('stats:RegressionKernel:InvalidIterationLimit'))
                end

                Xm = map(obj.FeatureMapper,X',obj.KernelScale)';

                obj.Impl = classreg.learning.impl.LinearImpl.make(0,...
                                                                  initialBeta,initialBias,...
                                                                  Xm,Y,W,...
                                                                  modelParams.LossFunction,...
                                                                  strcmp(modelParams.Regularization,'ridge'),...
                                                                  lambda,...
                                                                  [],...  % PassLimit
                                                                  [],...  % BatchLimit
                                                                  [],...  % NumCheckConvergence
                                                                  [],...  % BatchIndex
                                                                  [],...  % BatchSize
                                                                  modelParams.Solver,...
                                                                  modelParams.BetaTolerance,...
                                                                  modelParams.GradientTolerance,...
                                                                  1e-6,... % DeltaGradientTolerance
                                                                  [],...   % LearnRate
                                                                  [],...   % OptimizeLearnRate
                                                                  [],...   % ValidationX
                                                                  [],...   % ValdiationY
                                                                  [],...   % ValidationW
                                                                  modelParams.IterationLimit,...
                                                                  [],...   % TruncationPeriod,...
                                                                  modelParams.FitBias,...
                                                                  modelParams.PostFitBias,...
                                                                  obj.Epsilon,... % Epsilon
                                                                  modelParams.HessianHistorySize,...
                                                                  modelParams.LineSearch,...
                                                                  0, ... % Consensus
                                                                  [],... % Stream
                                                                  modelParams.VerbosityLevel);

                fiHis = obj.Impl.FitInfo.History;
                if isempty(fiHis)
                    fiHisNew = [];
                else
                    % Adapt the history of the Fast Linear Solver to
                    % RegressionKernel's history
                    fiHisNew = struct('ObjectiveValue',fiHis.Objective,...
                                      'GradientMagnitude',fiHis.Gradient,...
                                      'Solver',categorical(repmat({'LBFGS-fast'},numel(fiHis.Solver),1)),...
                                      'IterationNumber',fiHis.NumIterations,...
                                      'DataPass',cumsum(fiHis.NumPasses),...
                                      'RelativeChangeInBeta',fiHis.RelativeChangeInBeta,...
                                      'ElapsedTime',[nan(numel(fiHis.Solver)-1,1);obj.Impl.FitInfo.FitTime]);
                end

                fitInfo = struct('Solver','LBFGS-fast',...
                                 'LossFunction',modelParams.LossFunction,...
                                 'Lambda',lambda,...
                                 'BetaTolerance', modelParams.BetaTolerance,...
                                 'GradientTolerance',modelParams.GradientTolerance,...
                                 'ObjectiveValue',obj.Impl.FitInfo.Objective,...
                                 'GradientMagnitude',obj.Impl.FitInfo.GradientNorm,...
                                 'RelativeChangeInBeta',obj.Impl.FitInfo.RelativeChangeInBeta,...
                                 'FitTime',obj.Impl.FitInfo.FitTime,...
                                 'History',fiHisNew);

            else
                if (modelParams.IterationLimit + modelParams.ADMMIterationLimit) == 0
                    error(message('stats:RegressionKernel:InvalidIterationLimitADMM'))
                end
                %%%%%%%% FIT THE DATA WITH A BLOCK-WISE SOLVER %%%%%%%%%%%%%%
                X = X';
                if modelParams.VerbosityLevel>0
                    fprintf(getString(message('stats:RegressionKernel:FoundNBlocks',numberChunks)))
                end

                Beta = [initialBias;initialBeta]; %including Bias

                doridge = strcmp(modelParams.Regularization,'ridge');
                rho = modelParams.Consensus;

                Wk = zeros(numberChunks,1);
                X_ = cell(numberChunks,1);
                Y_ = cell(numberChunks,1);
                W_ = cell(numberChunks,1);
                h = 0;
                for i = 1:maxChunkSize:n
                    h = h + 1;
                    j = min(i+maxChunkSize-1,n);
                    Wk(h) = sum(W(i:j,:),1);
                    X_{h} = X(i:j,:);
                    Y_{h} = Y(i:j,:);
                    W_{h} = W(i:j,:);
                end
                hfixedchunkfun = @fixedchunkfun;
                hmarkforreuse = @(x)x; % Doesn't matter because X is not tall
                hupdateforreuse = @(y,x)x; % Doesn't matter because X is not tall

                betaTol = modelParams.BetaTolerance;
                gradTol = modelParams.GradientTolerance;
                admmIterationLimit = modelParams.ADMMIterationLimit;
                tallPassLimit = intmax;
                verbose = modelParams.VerbosityLevel;

                doBias = modelParams.FitBias;
                lossfun = modelParams.LossFunction;

                % Parameters for the LBFGS inside ADMM iterations:
                iterationlimit_ADMMLBFGS = [modelParams.WarmStartIterationLimit modelParams.ADMMUpdateIterationLimit];
                hessianHistorySize_ADMMLBFGS = 15;  % fixed, cannot be controlled by user
                dowolfe_ADMMLBFGS = true;           % fixed, cannot be controlled by user
                doBias_ADMMLBFGS  = doBias;         % already restricted to always true above
                gradTol_ADMMLBFGS = gradTol;        % same as outer loop
                betaTol_ADMMLBFGS = betaTol;        % same as outer loop
                doridge_ADMMLBFGS = doridge;        % already restricted to always true above
                lossfun_ADMMLBFGS = lossfun;        % already set to 'mse' or 'epsiloninsensitive' above

                % Make the objective function
                expType = modelParams.Transformation;
                objgraF = makeobjgradF(X_,Y_,W_,lossfun,lambda./numberChunks,obj.FeatureMapper,obj.KernelScale,obj.Epsilon);
                % Wrap objective function with verbosity and progress tracking behavior
                hclientfun = @(x)(x); % Doesn't matter because X is not tall
                keepHist = true;
                progressF = classreg.learning.linearutils.linearSolverProgressFunction(objgraF,false,admmIterationLimit>1,verbose>0,keepHist,hclientfun);

                %%%%%%%%%%%%%%%%%%%%%%%%%%% ADMM FIT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                Beta = classreg.learning.linearutils.ADMMimpl(X_,Y_,W_,Beta,rho,lambda,doridge,...
                                                              numberChunks,Wk,hmarkforreuse,hupdateforreuse,hfixedchunkfun,...
                                                              betaTol,gradTol,admmIterationLimit,tallPassLimit,progressF,verbose,...
                                                              {lossfun_ADMMLBFGS,obj.Epsilon},betaTol_ADMMLBFGS,gradTol_ADMMLBFGS,...
                                                              iterationlimit_ADMMLBFGS, doridge_ADMMLBFGS,...
                                                              hessianHistorySize_ADMMLBFGS,dowolfe_ADMMLBFGS,doBias_ADMMLBFGS,...
                                                              obj.FeatureMapper,expType,obj.KernelScale,false);

                %%%%%%%%%%%%%%%%%%%%%%%%%%% CENTRALIZED LBFGS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                hessianHistorySize = modelParams.HessianHistorySize;
                lineSearch = modelParams.LineSearch;
                iterationlimit = modelParams.IterationLimitBlockWise;
                initialStepSize = modelParams.InitialStepSize;

                Beta = classreg.learning.linearutils.LBFGSimpl(Beta,progressF,verbose,...
                                                               betaTol,gradTol,iterationlimit,tallPassLimit,...
                                                               hessianHistorySize,lineSearch,initialStepSize);

                beta = Beta(2:end);
                bias = Beta(1);

                if verbose>0
                    printLastLine(progressF);
                end

                modelParams.Lambda = lambda; % Impl keeps the real lambda in case modelParams is still 'auto'

                obj.Impl = classreg.learning.impl.LinearImpl.makeNoFit(modelParams,beta,bias,progressF.History);

                fiHis = obj.Impl.FitInfo;
                if verbose>0
                    % Adapt the history of linearSolverProgressFunction to
                    % RegressionKernel's history
                    fiHisNew = struct('ObjectiveValue',fiHis.ObjectiveValue(:),...
                                      'GradientMagnitude',fiHis.GradientMagnitude(:),...
                                      'Solver',categorical(repmat({'LBFGS-blockwise'},numel(fiHis.Solver),1)),...
                                      'IterationNumber',fiHis.IterationNumber(:),...
                                      'DataPass',fiHis.DataPass(:),...
                                      'RelativeChangeInBeta',fiHis.RelativeChangeBeta(:),...
                                      'ElapsedTime',fiHis.ElapsedTime(:));
                else
                    fiHisNew = [];
                end

                fitInfo = struct('Solver','LBFGS-blockwise',...
                                 'LossFunction',modelParams.LossFunction,...
                                 'Lambda',modelParams.Lambda,...
                                 'BetaTolerance', modelParams.BetaTolerance,...
                                 'GradientTolerance',modelParams.GradientTolerance,...
                                 'ObjectiveValue',obj.Impl.FitInfo.ObjectiveValue(end),...
                                 'GradientMagnitude',obj.Impl.FitInfo.GradientMagnitude(end),...
                                 'RelativeChangeInBeta',obj.Impl.FitInfo.RelativeChangeBeta(end),...
                                 'FitTime',obj.Impl.FitInfo.ElapsedTime(end),...
                                 'History',fiHisNew);


            end

            obj.ModelParams = modelParams;

            % save fitInfo into a private property, used only by resume method to
            % concatenate history information. User should only obtain
            % fitInfo from the second output argument.
            obj.FitInfo = fitInfo;

        end

        function obj = fromStruct(s)
            % Make a Kernel object from a codegen struct.

            s.ResponseTransform = s.ResponseTransformFull;

            % Remove fields added for codegen
            s.DataSummary = rmfield(s.DataSummary,{'NumExpandedPredictors','NumCategories'});
            s = rmfield(s,'MaximumChunkSize');

            modelTrainedWithTallArray = s.DataSummary.ModelTrainedWithTallArray;
            s.DataSummary = rmfield(s.DataSummary,'ModelTrainedWithTallArray');
            s = classreg.learning.coderutils.structToRegr(s);

            if isfield(s.DataSummary,'RemoveObservationsWithMissingValues')
                s.DataSummary = rmfield(s.DataSummary,'RemoveObservationsWithMissingValues');
            end

            if modelTrainedWithTallArray
                % Model trained with non-table data but has PredictorNames,
                % must have been trained with tall data.
                s.DataSummary.VariableRange = {};
                s.DataSummary = rmfield(s.DataSummary,'BinnedX');
                s.DataSummary = rmfield(s.DataSummary,'BinEdges');
            end            

            beta = s.Impl.Beta;
            bias = s.Impl.Bias;
            lambda = s.ModelParams.Lambda;
            featureMapper = classreg.learning.rkeutils.featureMapper.fromStruct(s.FeatureMapper);
            kernelScale = s.ModelParams.KernelScale;
            boxConstraint = s.ModelParams.BoxConstraint;
            modelParams = classreg.learning.coderutils.coderStructToKernelParams(s.ModelParams,modelTrainedWithTallArray);
            dataSummary = s.DataSummary;
            fitinfo = s.Impl.FitInfo;
            responseTransform = s.ResponseTransform;
            if isfield(s,'Mu') 
                mu=s.Mu;
                sigma=s.Sigma;
            else  % for 23a or earlier releases
                mu=[];
                sigma=[];
            end
            %handle tall case separately
            if modelTrainedWithTallArray &&  strcmpi(s.Impl.LossFunction,'mse')
                epsilon = [];
            else
                epsilon = s.Impl.Epsilon;
            end

            % Make an object
            obj = RegressionKernel.makeNoFit(beta,bias,lambda,featureMapper,...
                kernelScale,boxConstraint,epsilon,...
                modelParams,dataSummary,fitinfo,responseTransform,mu,sigma,s.DefaultResponse);

            obj.FitInfo = s.FitInfo;
            % Override Impl using fromStruct
            try
                obj.Impl = classreg.learning.impl.LinearImpl.fromStruct(s.Impl);
            catch ME
                if contains(ME.identifier,'stats:classreg:learning:impl:LinearImpl:make:BadStatus') % reset termination code MaxIter reached
                    s.Impl.FitInfo.TerminationCode = 0;
                    obj.Impl = classreg.learning.impl.LinearImpl.fromStruct(s.Impl);
                else
                    rethrow(ME);
                end
            end 
        end 

        function name = matlabCodegenRedirect(~)
            name = 'classreg.learning.coder.regr.RegressionKernel';
        end 
    end

    methods (Hidden)
        function s = toStruct(this,varargin)
            % Convert to a struct for codegen.

            warnState  = warning('query','all');
            warning('off','MATLAB:structOnObject');
            cleanupObj = onCleanup(@() warning(warnState));

            % convert common properties to struct
            s = classreg.learning.coderutils.regrToStruct(this);

            modelTrainedWithTallArray = false;
            if ischar(s.DataSummary.PredictorNames) && isfield(s.DataSummary,'TableInput') && ~s.DataSummary.TableInput
                modelTrainedWithTallArray = true;
            end    
            s.DataSummary.ModelTrainedWithTallArray = modelTrainedWithTallArray;             
    
            s.Impl = toStruct(this.Impl);
            s.FeatureMapper = toStruct(this.FeatureMapper);

            % model params
            s.ModelParams = classreg.learning.coderutils.kernelParamsToCoderStruct(...
                this.ModelParams);            

            % save the path to the fromStruct method
            s.FromStructFcn = 'RegressionKernel.fromStruct';

            s.DataSummary.NumExpandedPredictors = numel(getExpandedPredictorNames(this));
            s.DataSummary.NumCategories = cellfun(@numel,this.VariableRange);
            s.MaximumChunkSize = RegressionKernel.estimateMaxChunkSize(this.FeatureMapper.n,this.ModelParams.BlockSize,false);
            s.FitInfo = this.FitInfo;
            s.Mu = this.Mu;
            s.Sigma = this.Sigma;
        end
    end

    methods(Static,Hidden)        
        function this = loadobj(obj)
            this = loadobj@classreg.learning.Kernel(obj);
            if isempty(this.AllowedPredictArgs)
                this.AllowedPredictArgs = struct('Name','PredictionForMissingValue','Value','Median');
            end
        end
    end
end

function objgraF = makeobjgradF(X,Y,W,lossfun,lambdaK,FM,kernelScale,epsilon)
    objgraF = @fcn;
    function [obj,gra] = fcn(Beta)
        [obj,gra] = fixedchunkfun(1, @(x,y,w) ...
                                  chunkObjGraFun(Beta,x,y,w,lossfun,lambdaK,FM,kernelScale,epsilon), ...
                                  X,Y,W);
        obj = sum(obj,1);
        gra = sum(gra,1);
    end
end

function  varargout = fixedchunkfun(FixedNumSlices,fcn,varargin)
% Replicates tall.fixedchunkfun on non-tall inputs so we can use the same
% implementation of the solvers. Assumes all variables in varargin will be
% partitioned (i.e. no broadcasted inputs). This functions is called in two
% places: 1) from the objgra functor (above), and 2) from the Beta update
% in the ADMM algorithm.

    n = size(varargin{1},1); % height
    nia = numel(varargin);
    varargout = cell(1,nargout);
    h = 0;
    for i = 1:FixedNumSlices:n
        h = h +1;
        j = min(i+FixedNumSlices-1,n);
        invar = cell(1,nia);
        outvar = cell(1,nargout);
        for k = 1:nia
            invar{k} = varargin{k}(i:j,:);
        end
        [outvar{1:nargout}] = fcn(invar{:});
        for k = 1:nargout
            varargout{k}(h,:) = outvar{k};
        end
    end

end

function [obj,gra] = chunkObjGraFun(Beta,x,y,w,lossfun,lambda,FM,kernelScale,epsilon)
    x = x{1};
    y = y{1};
    w = w{1};
    useBias = true;
    doridge = true;
    if isempty(x)
        obj = (Beta(2:end)'*Beta(2:end))*lambda/2;
        gra = [0;Beta(2:end)*lambda]';
    else
        xm = map(FM,x,kernelScale);
        [obj,gra] = classreg.learning.linearutils.objgrad( Beta(2:end),Beta(1),xm',y,w,...
                                                           lossfun,doridge,lambda,epsilon,useBias,...
                                                           zeros(numel(Beta)-1,1,'like',Beta),0,0);
        gra = [gra(end);gra(1:end-1)]';
        % obj = sum(max(0,1-(x*Beta(2:end)+Beta(1)).*y).*w)+(Beta(2:end)'*Beta(2:end))*lambda/2;
        % gra = (sum((y.*(x*Beta(2:end)+Beta(1))<1) .* (-y.*[ones(numel(w),1) x].*w))' + [0;Beta(2:end)*lambda])';
    end
end

% LocalWords:  INCREMENTALLEARNER INCREMENTALMODEL Warmup tt Xoffline Yoffline Xnew Ynew Wnew
% LocalWords:  NBlocks fixedchunkfun objgra gra predictionformissingvalue Valdiation
