function [p, solverState] = lbfgsupdate(p, lossFcn, solverState, NameValueArgs)
    %lbfgsupdate   Update parameters using L-BFGS
    %
    %   [netUpdated,solverStateUpdated] = lbfgsupdate(net,lossFcn,solverState)
    %   updates the learnable parameters in the dlnetwork net using the
    %   L-BFGS algorithm with a line search and also returns the updated
    %   solver state. lossFcn is a function handle with syntax
    %   [loss, gradients] = f(net) that outputs the loss and the gradients
    %   with respect to the parameters in net, specified as a function
    %   handle. The state solverState must be a LBFGSState object or []. To
    %   create an LBFGSState object, use lbfgsState. If solverState is [],
    %   then the function uses solverState=lbfgsState.
    %
    %   Tip: To parameterize a model loss function that has additional
    %   inputs and at least one call to dlgradient, specify the loss
    %   function lossFcn=@(net)dlfeval(@modelLoss,net,arg1,...,argN) where
    %   modelLoss is a function that takes as input the network net and
    %   additional arguments arg1,...,argN and outputs the loss and gradients.
    % 
    %   [paramsUpdated,solverStateUpdated] = lbfgsupdate(params,lossFcn,solverState)
    %   updates the learnable parameters in the dlarray, numeric array,
    %   table, structure, or cell array params.
    %
    %   [___] = lbfgsupdate(___,Name=Value) specifies additional options
    %   using one or more name-value arguments:
    %       LineSearchMethod  - Method to find suitable learning rate, 
    %                           specified as one of these values:
    %                           "weak-wolfe"      - Search for a learning
    %                                               rate which satisfies
    %                                               the weak Wolfe conditions.
    %                                               This method maintains a
    %                                               positive definite
    %                                               approximation of the
    %                                               inverse Hessian.
    %                           "strong-wolfe"    - Search for a learning
    %                                               rate which satisfies
    %                                               the strong Wolfe conditions.
    %                                               This method maintains a
    %                                               positive definite
    %                                               approximation of the
    %                                               inverse Hessian.
    %                           "backtracking"    - Search for a learning
    %                                               rate which satisfies
    %                                               sufficient decrease
    %                                               conditions. This method
    %                                               does not maintain a
    %                                               positive definite
    %                                               approximation of the
    %                                               inverse Hessian.
    %                           The default value is "weak-wolfe".
    %
    %       MaxNumLineSearchIterations
    %                         - Maximum number of line search iterations,
    %                           specified as a positive integer. 
    %                           The default value is 20.
    %
    %       NumLossFunctionOutputs
    %                         - Number of outputs of lossFcn, specified as
    %                           an integer greater than or equal to 2. 
    %                           The default value is 2.
    %
    %   Example:
    %       % Perform a single epoch of training on a network that classifies 
    %       % iris flowers into three species.
    %
    %       % Load training data and one-hot encode the labels.
    %       [XTrain, TTrain] = iris_dataset;
    %       XTrain = dlarray(XTrain,'CB');
    %
    %       % Define the network.
    %       numInputFeatures = size(XTrain,1);
    %       numClasses = size(TTrain,1);
    %       numHiddenUnits = 32;
    %       layers = [
    %           featureInputLayer(numInputFeatures)
    %           fullyConnectedLayer(numHiddenUnits)
    %           reluLayer
    %           fullyConnectedLayer(numHiddenUnits)
    %           reluLayer
    %           fullyConnectedLayer(numClasses)
    %           softmaxLayer];
    %       net = dlnetwork(layers);
    %
    %       % Create loss function handle.
    %       lossFunction = @(net)dlfeval(@modelLoss, net, XTrain, TTrain);
    %
    %       % Initialize L-BFGS solver state with a maximum history size of
    %       % 3 and an initial inverse Hessian approximation factor with
    %       % value 1.1.
    %       solverState = lbfgsState(HistorySize=3, InitialInverseHessianFactor=1.1);
    %
    %       % Perform a single update on all the training data using a weak
    %       % Wolfe line search with a maximum of 10 line search iterations.
    %       [net, solverState] = lbfgsupdate(net, lossFunction, solverState, ...
    %           LineSearchMethod="weak-wolfe", MaxNumLineSearchIterations=10);
    %
    %       % Calculate the cross-entropy loss and the gradients of the
    %       % loss with respect to learnable parameters.
    %       function [loss, gradients] = modelLoss(net, X, T)
    %           loss = crossentropy(forward(net, X), T);
    %           gradients = dlgradient(loss,net.Learnables);
    %       end
    %
    %   See also lbfgsState, deep.LBFGSState, adamupdate, sgdmupdate, rmspropupdate.

    %   Copyright 2022 The MathWorks, Inc.

arguments
    p
    lossFcn                                         {iMustBeFunctionHandle}
    solverState                                     {iAssertValidSolverState}
    NameValueArgs.LineSearchMethod                  string {mustBeMember(NameValueArgs.LineSearchMethod ,["weak-wolfe","strong-wolfe","backtracking"])} = "weak-wolfe"
    NameValueArgs.MaxNumLineSearchIterations(1,1)   {mustBeNumeric, mustBeInteger, mustBePositive, mustBeFinite} = 20
    NameValueArgs.NumLossFunctionOutputs(1,1)       {mustBeNumeric, mustBeInteger, mustBeGreaterThanOrEqual(NameValueArgs.NumLossFunctionOutputs, 2), mustBeFinite} = 2
end

% Convert arguments to canonical form.
NameValueArgs = iConvertToCanonicalForm(NameValueArgs);

if isempty(solverState)
    % Initialize LBFGSState object with default values.
    solverState = lbfgsState();
end

[solverState, p] = solverState.step(p, lossFcn, NameValueArgs.LineSearchMethod, NameValueArgs.MaxNumLineSearchIterations, NameValueArgs.NumLossFunctionOutputs);
end

function iMustBeFunctionHandle(value)
if ~(isa(value, 'function_handle') || isa(value, 'deep.AcceleratedFunction'))
    error(message('deep:lbfgsupdate:funArgNotHandle'))
end
end

function iAssertValidSolverState(value)
if ~((isnumeric(value) && isempty(value)) || isa(value, 'deep.LBFGSState'))
    error(message('deep:lbfgsupdate:stateArgInvalid'))
end
end

function inputArguments = iConvertToCanonicalForm(params)
% Make sure values are converted to double
inputArguments = struct;
inputArguments.LineSearchMethod = params.LineSearchMethod;
inputArguments.MaxNumLineSearchIterations = double(params.MaxNumLineSearchIterations);
inputArguments.NumLossFunctionOutputs = double(params.NumLossFunctionOutputs);
end