You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: fit_nlinOspreyMultiplex.m
    path of file: ./repos/osprey/libraries/FID-A/fitTools/fitModels/Osprey
    
    the code of the file until where you have to start completion: function F = fit_nlinOspreyMultiplex(x, dataToFitInput, basisSet, spl_pos, fitRangePPM)
%   This function receives the current set of non-linear parameters during
%   every iteration of the non-linear least-squares optimization. The
%   parameters are applied to the basis set. 
%
%   Subsequently, a fast non-negative linear least-squares (FNNLS) algorithm
%   is applied to determine the optimal amplitude to the current set of
%   non-linear parameters.
%
%   This function returns the difference between the input data and
%   the current optimized model. This difference is used by the NLLS
%   solver.
%
%   USAGE:
%       F           = fit_nlinOspreyMultiplex(x, dataToFit, basisSet, spl_pos, fitRangePPM);
%
%   INPUTS:
%       data        = Vector containing the concatenated real and imaginary
%                     values of an acquired signal.
%       basis       = Vector containing the concatenated real and imaginary
%                     values of a simulated basis set.
%       x           = Vector providing initial values to the non-linear

    
function F = fit_nlinOspreyMultiplex(x, dataToFitInput, basisSet, spl_pos, fitRangePPM)
%   This function receives the current set of non-linear parameters during
%   every iteration of the non-linear least-squares optimization. The
%   parameters are applied to the basis set. 
%
%   Subsequently, a fast non-negative linear least-squares (FNNLS) algorithm
%   is applied to determine the optimal amplitude to the current set of
%   non-linear parameters.
%
%   This function returns the difference between the input data and
%   the current optimized model. This difference is used by the NLLS
%   solver.
%
%   USAGE:
%       F           = fit_nlinOspreyMultiplex(x, dataToFit, basisSet, spl_pos, fitRangePPM);
%
%   INPUTS:
%       data        = Vector containing the concatenated real and imaginary
%                     values of an acquired signal.
%       basis       = Vector containing the concatenated real and imaginary
%                     values of a simulated basis set.
%       x           = Vector providing initial values to the non-linear
%                     least-squares solver.
%
%   OUTPUTS:
%       F           = Difference between data and model. This is the
%                     objective to be least-squares-optimized by the 
%                     non-linear solver.
%
%   AUTHOR:
%       Dr. Georg Oeltzschner (Johns Hopkins University, 2019-06-04)
%       goeltzs1@jhmi.edu
%
%   CREDITS:
%       This code uses the 'fnnls' implementation of a fast non-negative
%       least-squares optimization algorithm by Rasmus Bro.
%       https://www.mathworks.com/matlabcentral/fileexchange/3388-nnls-and-constrained-regression
%       Bro et al., Journal of Chemometrics 11:393-401 (1997)
%
%       This code is based on numerous functions from the FID-A toolbox by
%       Dr. Jamie Near (McGill University)
%       https://github.com/CIC-methods/FID-A
%       Simpson et al., Magn Reson Med 77:23-33 (2017)
%
%   HISTORY:
%       2019-06-04: First version of the code.


% 0. Loop over the number of data sets to be multiplex-fitted
nMultiplex = length(dataToFitInput);
nBasisFcts  = basisSet.nMets + basisSet.nMM; % number of basis functions

% 1. Extract parameters
zeroPhase   = x(1); % zero-order phase correction [rad]
firstPhase  = x(2); % first-order phase correction [rad]
gaussLB     = x(3); % Gaussian damping [Hz^2]
lorentzLB   = x(4:nBasisFcts+3); % Lorentzian damping [Hz] for each basis function
freqShift   = x(nBasisFcts+4:2*nBasisFcts+3); % Frequency shift [Hz] for each basis function
% We want to construct a cubic spline (i.e. of order 4); therefore a
% spline with k knots requires k-4 coefficients.
beta_j      = x(2*nBasisFcts+4:end