Error using fittype for custom function: size and shape of input/output

Dear community, I am currently working on a code that could reproduce a so-called blaze function, a modulation of light peaks obtained in a spectrometer. Tha blaze function is a function of the wavelength, but must be contructed and is not a trivial mathematical expression, since it depends on angles an diffraction orders. In MATLAB, I am implementing a function that reproduces the profile of the blaze function, given as input:

• the wavelengths (vector, nm) on which the blaze function is modelled,
• parameters of the instrument (incident angle alpha, blaze angle of the facets phi and facet distance d)

the output profile of the function should have the same length as the input wavelength vector. The idea is to retrieve the parameters of the instrument alpha, phi and d by fitting the function to the experimental curve using fit() and fittype(). However, fittype keeps returning an error:

Error using fittype/testCustomModelEvaluation
Custom equations must produce an output vector,
matrix, or array that is the same size and shape as
the input data. This custom equation fails to meet
that requirement:

 blazeFunc(x,alpha,phi,d)

Error in fittype>iCreateFittype (line 373)
    testCustomModelEvaluation( obj );

Error in fittype (line 330)
                obj = iCreateFittype( obj, varargin{:} );

Error in blaze_function_manual (line 48)
fty = fittype('blazeFunc(x,alpha,phi,d)');

the function is the following

function profile = blazeFunc(wav,alpha,phi,d)
    
    %% Functions
    % These functions are derived from the physics of the instrument
    sinBetaBar = @(n,lam,d,alpha,phi) (n./d .* (lam*10^(-9)) - sin(alpha+phi)); % sin of the diffraction angle of order n and wavelength lam
    sinBetaBar0 = @(n,k,alpha,phi) (sin(alpha+phi) - 2*k*tan(phi)*cos(alpha+phi)./(n*cos(phi))); % sin of the diffraction angle at which the blaze peak of order n at wavelength lam is zero (cuts the function around the main peak)
    lam0 = @(n,k,d,alpha,phi) (d./n .* (sin(alpha+phi) + sinBetaBar0(n,k,alpha,phi)))*10^9; % wavelength at which the peak of order n is zero
    phi_ce = @(n,sBB,d,alpha,phi) ((n .* pi .* cos(phi) .* (sin(alpha+phi)-sBB))./ (2 .* tan(phi) .* cos(alpha+phi))); % phase difference to constuct the blaze function
    blazeFun = @(n,sBB,d,alpha,phi) (sin(phi_ce(n,sBB,d,alpha,phi)) ./ phi_ce(n,sBB,d,alpha,phi)).^2; % formula fo the blaze function of order n and angle beta
    
    orders = 68:246; % Key parameter form the instrument

    %% Shaping input
    wav = reshape(wav,1,[]); % Shapes input for correct computation (must be a row vector, input is a column vector)

    %% Construction of the function
    nOrd = length(orders); % Number of orders
    WavArray = repmat(wav,nOrd,1); % Repeats the wavelengths in an array for each order
    OrderArray = repmat(orders',1,length(WavArray)); % Repeats the order vector for each wavelength 
    
    % Initialization
    sBB = zeros(nOrd,length(WavArray));
    values = zeros(nOrd,length(WavArray));
    
    % Computation
    lambda0 = [lam0(orders,1,d,alpha,phi);lam0(orders,-1,d,alpha,phi)]; % Array in which each column is the limit in wavelength of each peak of each order
    OrderIndex = (WavArray>=lambda0(1,:)' & WavArray<=lambda0(2,:)'); % Define an index to compute only relevant elements
    sBB(OrderIndex) = sinBetaBar(OrderArray(OrderIndex),WavArray(OrderIndex),d,alpha,phi); % Compute for each wavelength (around the peak of each order) the sin of the diffraction angle
    values(OrderIndex) = blazeFun(OrderArray(OrderIndex),sBB(OrderIndex),d,alpha,phi); % Build the peaks of the blaze function for each order and sin(beta) range

    %% Output
    profile = max(values,[],1)'; % Peaks overlap, but the instrument removes automatically the overlapping, so the max at each wavelength is taken. The output is transposed as a column vector so it matches the input.    
end

The defined function is used here

%% Fit
range = wav.raw{1} > 200 & wav.raw{1} < 900; % Set the range of data to use (all) (in nm)
x = wav.raw{1}(range); % Is a column of values
y = int.norm{1}(range); % Is a column of intensities (spectrum)

fty = fittype('blazeFunc(x,alpha,phi,d)');...
    independent="x",...
    dependent="y",...
    coefficients={'alpha','phi','d'});
testFit = fit(x,y,fty);

From which the error occurs. I checked the output of blazeFunc(), and it has the same size as the input, so I really do not understand where the problem lies.

Note: for now the blazeFunc() output does not look like the original data (it´s just a series of peaks with height 1). As the blaze function si normalised, it can be multiplied by the spectrum to obtain a good match between the two and an accurate fit. For now, this operation is excluded for the sake of troubleshooting.

Thank you all for any help with this. I know that fittype() can behave wierdly with sizes, but I can't understand here where the issues is.

[Mathieu NOE]

would be nice to have some data to play with

Answer 1

[Shivansh]

Hi Francesco!

It seems like you are getting errors due to the "fittype" function in MATLAB. The error message suggests that the error occurs when the custom function used for fitting does not return an output array that matches the size of the input data array.

I have checked the function using sample data points and it is working as output and input size matches for R2023b.

% Example input
wav_example = linspace(200, 900, 700)'; % 700 points from 200nm to 900nm
alpha_example = 30; % Example value
phi_example = 45; % Example value
d_example = 1.5; % Example value
% Call the function
profile_example = blazeFunc(wav_example, alpha_example, phi_example, d_example);
% Check size
assert(length(profile_example) == length(wav_example), 'Output size mismatch.');
disp("Input and output size matches");

Input and output size matches

function profile = blazeFunc(wav,alpha,phi,d)

%% Functions
% These functions are derived from the physics of the instrument
sinBetaBar = @(n,lam,d,alpha,phi) (n./d .* (lam*10^(-9)) - sin(alpha+phi)); % sin of the diffraction angle of order n and wavelength lam
sinBetaBar0 = @(n,k,alpha,phi) (sin(alpha+phi) - 2*k*tan(phi)*cos(alpha+phi)./(n*cos(phi))); % sin of the diffraction angle at which the blaze peak of order n at wavelength lam is zero (cuts the function around the main peak)
lam0 = @(n,k,d,alpha,phi) (d./n .* (sin(alpha+phi) + sinBetaBar0(n,k,alpha,phi)))*10^9; % wavelength at which the peak of order n is zero
phi_ce = @(n,sBB,d,alpha,phi) ((n .* pi .* cos(phi) .* (sin(alpha+phi)-sBB))./ (2 .* tan(phi) .* cos(alpha+phi))); % phase difference to constuct the blaze function
blazeFun = @(n,sBB,d,alpha,phi) (sin(phi_ce(n,sBB,d,alpha,phi)) ./ phi_ce(n,sBB,d,alpha,phi)).^2; % formula fo the blaze function of order n and angle beta

orders = 68:246; % Key parameter form the instrument
%% Shaping input
wav = reshape(wav,1,[]); % Shapes input for correct computation (must be a row vector, input is a column vector)
%% Construction of the function
nOrd = length(orders); % Number of orders
WavArray = repmat(wav,nOrd,1); % Repeats the wavelengths in an array for each order
OrderArray = repmat(orders',1,length(WavArray)); % Repeats the order vector for each wavelength 

% Initialization
sBB = zeros(nOrd,length(WavArray));
values = zeros(nOrd,length(WavArray));

% Computation
lambda0 = [lam0(orders,1,d,alpha,phi);lam0(orders,-1,d,alpha,phi)]; % Array in which each column is the limit in wavelength of each peak of each order
OrderIndex = (WavArray>=lambda0(1,:)' & WavArray<=lambda0(2,:)'); % Define an index to compute only relevant elements
sBB(OrderIndex) = sinBetaBar(OrderArray(OrderIndex),WavArray(OrderIndex),d,alpha,phi); % Compute for each wavelength (around the peak of each order) the sin of the diffraction angle
values(OrderIndex) = blazeFun(OrderArray(OrderIndex),sBB(OrderIndex),d,alpha,phi); % Build the peaks of the blaze function for each order and sin(beta) range
%% Output
profile = max(values,[],1)'; % Peaks overlap, but the instrument removes automatically the overlapping, so the max at each wavelength is taken. The output is transposed as a column vector so it matches the input.    
end

Although the function works as expected for the example input, the fit may provide inputs during the fitting process that lead to unexpected behavior. You can check this by adding "assert" in "blazeFunc" to include input and output validation checks.

function profile = blazeFunc(wav,alpha,phi,d)
% Your existing function code here

% Add validation at the end of your function before returning the output
assert(isequal(size(profile), size(wav)), 'Output profile does not match input wav size.');
end

The above code will provide an error on encountering an input-output mismatch and point you towards the edge case.

You can also look for input and output data ranges and check for numerical instabilities leading to "inf" or "NaN" values.

if any(isnan(values(:))) || any(isinf(values(:)))
    warning('NaNs or Infs detected in calculations.');
end

You can refer to the following documentation for more information about the "fittype" function in MATLAB:

https://www.mathworks.com/help/curvefit/fittype.html.

I hope it helps!

[Francesco]

Thank you very much for your input. This will be very useful for future projects. Since, as you confirmed, fittype may be picky with data instabilities and the code revolving around this has expanded considerably, I opted for a different fitting algorithm (nlinfit) which seems more stable. I will be coming back to this answer in the future for sure!