Exiting: Maximum number of function evaluations has been exceeded

Hi,

I used ezyfit for gaussian functions. I got such errors.

Exiting: Maximum number of function evaluations has been exceeded

- increase MaxFunEvals option.

Current function value: 142.981920

Exiting: Maximum number of function evaluations has been exceeded

- increase MaxFunEvals option.

Current function value: 721.060509

I added "options = optimset('MaxFunEvals',1000);" but the error still exists.

Please suggest. thanks a lot!

load ('orientcurve.mat');
OSI=zeros(size(orientcurve,1),1);
tunewidth=zeros(size(orientcurve,1),1);
tunepref=NaN(size(orientcurve,1),1);
prefdF=zeros(size(orientcurve,1),1);
options = optimset('MaxFunEvals',1000);
for i=1:size(orientcurve,1)  
    orient=[-90 -45 0 45 90];     
    [~,I]=max(orientcurve(i,:));
    calc_pref_or=I;
    prefdF(i,3)=orientcurve(i,I);   
    if I==4
        I=1;
    else
        I=I+1;
    end
    prefdF(i,4)=orientcurve(i,I); 
    if I==4
        I=1;
    else
        I=I+1;
    end
    prefdF(i,5)=orientcurve(i,I); 
    prefdF(i,1)=orientcurve(i,I);
    if I==4
        I=1;
    else
        I=I+1;
    end
    prefdF(i,2)=orientcurve(i,I); 
    fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); x_0=0');  
    sigma=fit.m(2);
    if abs(sigma)<15 
        fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); sigma=15; x_0=0');
        sigma=fit.m(2);
    end   
end
    

[John D'Errico]

Note that even though you did this:

options = optimset('MaxFunEvals',1000);

that is irrelevant, IF you have not passed it into the code in question, and if the code would not use that anyway. ezfit is NOT part of the optimization toolbox, and it may not take an options structure.

Software cannot read your mind.

[Walter Roberson]

@John D'Errico

The third party ezyfit toolbox http://www.fast.u-psud.fr/ezyfit/html/ezfit.html is specifically documented as using fminsearch() internally.

However, it is not documented as accepting any fminsearch() options -- not unless its fitparam function http://www.fast.u-psud.fr/ezyfit/html/fitparam.html allows that.

[John D'Errico]

I'd looked at it, and saw no capability to pass in any parameters, so an edit would be needed, and an edit done to someone else's code is just dangerous.

Answer 1

[Walter Roberson]

https://www.mathworks.com/matlabcentral/fileexchange/10176-ezyfit-2-44 is not documented as accepting options. The function would have to be edited to permit options.

[HYZ]

could you help how to edit the function to permit options? thanks!

[Walter Roberson]

At the bottom of fitparam.m add

fp.options = [];

In ezfit.m near line 550, change

m=fminsearch(@fitlin, m0);

to

m=fminsearch(@fitlin, m0, fp.options);

and near line 562 change

m=fminsearch(@fitlog, m0);

to

m=fminsearch(@fitlog, m0, fp.options);

Once that is done, if you call ezfit() and pass in 'options', followed by the result of an optimset() call, then it should use the options. For example in your code

options = optimset('MaxFunEvals', 1000, 'MaxIters', 1000);    
fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); x_0=0', 'options', options);  

[HYZ]

I did what you suggested. I got this error.

Error using optimset
Unrecognized parameter name 'MaxIters'.  Please see the options table in the documentation for a list of acceptable option
parameters. Note that some parameters are only supported by OPTIMOPTIONS. Link to options table

Error in tuningcurve03 (line 38)
    options = optimset('MaxFunEvals', 1000, 'MaxIters', 1000);

[Walter Roberson]

options = optimset('MaxFunEvals', 1000, 'MaxIter', 1000);    
fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); x_0=0', 'options', options);  

[HYZ]

If I remove MaxIters in optimset, I still could not run as it showed "error during the fminsearch procedure". thank you!

[Walter Roberson]

What error does it show?

[HYZ]

this is after removing MaxIters.

Error using ezfit

Fit: error during the fminsearch procedure

Error in tuningcurve03 (line 39)

fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); x_0=0', 'options', options);

[HYZ]

so I guess we can't change the options in ezyfit toolbox. Is that correct?

[Walter Roberson]

I have a modified version of the ezyfit() function that permits the options. But because you do not want the warning to appear, I am testing to see how many iterations are needed to avoid the message.

So far one of the curves has needed 775 million iterations, which has taken about 1 1/2 days of execution time. At this rate it might well encounter the 1 billion iteration limit that I set, and at that point I would hae to re-test taking at least another 2 or 3 days of execution just for one of the curves (3 of the curves require more than 100 million iterations each.)

Might I suggest that it might be reasonable to have an iteration limit?

The algorithm used by minsearch() has a... challenge. Suppose you have a surface which has an asymptope slope away from a central peak, but inside the central peak is the minima. Consider for example if you had the equivalent of x^2 - 1 in a narrow section of the middle, and outside that section you had exp(-x^2) . Now suppose your sample points happen to be outside the catch basin for the (x^2-1) section -- your sample points happen to be on the exp(-x^2) slope. Then at each point locally on that slope, you can always go "downhill", so you move away from the central peak. You do not arrive at a minima for exp(-x^2) until x reaches +/- infinity --- there is always somewhere further down slope (at least until you get to the point in the computation where to the limits of floating point representation, you get 0 anyhow.)

The slope-descending algorithm never stops and says, "Wait, how do I know that there isn't any other section that I have not explored yet that might be better?": it just assumes that downhill is always better.

This is pretty much the situation in the search that I have running at the moment: it is moving along a slope in steps of roughly -0.004 at a time, and it has reached coordinate in the negative billion or so range, but is still searching.

[HYZ]

sorry ... I used fitting toolbox for the first time. could you suggest some other toolbox which can fulfil the purpose of my code? thank you.

[Walter Roberson]

After a few days of calculation, curve #62 reached the 1 billion iteration limit that I had set, at which time the fit.m was

m: [25.4045182948109 -22315434600485.8 -33293590339.937]

which is to say that the second and third model parameters were trying to escape towards negative infinity.

[Walter Roberson]

Here is adjusted code.

In the below code, usega controls whether to switch to ga() [requires the Global Optimization Toolbox] or to use fminsearch() [requires no additional toolboxes]. It is configured to use ga if you have the license for the toolbox, but you could hard-code it as false if you wanted to compare results between ga and fminsearch.

Using ga() is a little out of the spirit of the EZYFIT toolbox, which was designed to avoid having to use "expensive" optional toolboxes. However, I set the default to not use ga(), and in your script I only turn ga() on if you have a license for it.

Using ga() is sometimes much faster.

Note: with the 1e6 iteration limit I use here for the fminsearch case, you will still get messages about running out of iterations. Those are likely to be situations where fminsearch() is going in completely the wrong direction. If you are using fminsearch() then any of the curves that the modified gives you a "FYI" about using a large number of iterations, you really should test to see if the results are valid.

%your main script

load ('orientcurve.mat');
numcurve = size(orientcurve,1);
OSI = zeros(numcurve,1);
tunewidth = zeros(numcurve,1);
tunepref = NaN(numcurve,1);
prefdF = zeros(numcurve,1);
sigma = NaN(numcurve,1);
time_taken = NaN(numcurve,1);

%use ga if you have a license for Global Optimization Toolbox
usega = license('test', 'GADS_Toolbox');

if ~usega
    options = optimset('MaxFunEvals', 1e6, 'MaxIter', 1e6);
else
    options = optimoptions('ga', 'MaxGenerations', 5000, 'Display', 'none');
end

for i=1:numcurve
    orient=[-90 -45 0 45 90];     
    [~,I]=max(orientcurve(i,:));
    calc_pref_or=I;
    prefdF(i,3)=orientcurve(i,I);   
    if I==4
        I=1;
    else
        I=I+1;
    end
    prefdF(i,4)=orientcurve(i,I); 
    if I==4
        I=1;
    else
        I=I+1;
    end
    prefdF(i,5)=orientcurve(i,I); 
    prefdF(i,1)=orientcurve(i,I);
    if I==4
        I=1;
    else
        I=I+1;
    end
    prefdF(i,2)=orientcurve(i,I);
    fprintf('Curve #%d\n', i);
    starttime = tic;
    fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); x_0=0', 'options', options, 'usega', usega);  
    sigma(i)=fit.m(2);
    if abs(sigma(i))<15 
        fit=ezfit(orient,prefdF(i,:),'y(x) = a*exp(-((x-x_0)^2)/(2*sigma^2)); sigma=15; x_0=0', 'options', options, 'usega', usega);
        sigma(i)=fit.m(2);
    end
    time_taken(i) = toc(starttime);
end

plot(time_taken)

fitparam had to be modified slightly

function fp=fitparam
%FITPARAM  Default settings for the EzyFit toolbox
%   This M-File contains the default settings for the EzyFit Toolbox,
%   such as the extrapolation mode, fit colors, equation box location etc.
%   See the page 'Settings' in the help browser for the list of available
%   settings.
%
%   To change the default settings, edit this file (type: edit fitparam.m)
%   and follow the instructions. You may also choose the item
%   'Default Settings' in the EzyFit menu.
%
%   FP = FITPARAM returns the structure FP containing the settings.
%
%   See also SHOWFIT, EZFIT, SHOWEQBOX, DISPEQFIT, PICKDATA.

%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 1.21,  Date: 2007/09/18
%   This function is part of the EzyFit Toolbox




% Fit color mode:
%  - a [R G B] vector specifies a fixed color (eg, [0 0 0] is black)
%  - a string specifies one of eight predefined colors (full name: 'red',
%    'blue' etc, or single letters, 'r', 'b' etc). Type 'doc colorspec' for
%    details.
%  - a numeric value specifies that the color of the fit has the same color
%    as the data, but multiplied by the factor fitcolor (<1 for darker,
%    >1 for lighter)
fp.fitcolor = [1 0.5 0];


% Width and style of the fit (type 'doc linespec' to see the available
% Matlab line specifications):
fp.fitlinewidth = 2;
fp.fitlinestyle = '-';


% Extrapolation mode:
%  'fig'  = extrapolates the fit to the figure limits
%  'data' = extrapolates the fit to the whole data limits, even if only a
%           selection of the data is fitted.
%  'none' = no extrapolation
fp.extrapol = 'none';


% Number of X points used to compute and display the fitted curve:
% (typically 20 to 500, default=200)
fp.npt = 200;


% Display the equation in the command windows: 'on' or 'off'
fp.dispeqmode = 'on';


% Display the equation box in the figure: 'on' or 'off'
fp.dispeqboxmode = 'on';


% Display a legend in the figure: 'on' or 'off'
fp.dispfitlegend = 'off';


% Equation replace mode:
%  'on'    replace each parameter by its numerical value in the equation
%  'off'   keep the parameter names in the equation 
fp.eqreplacemode = 'off';


% Correlation coefficient mode:
%  'r'      display R
%  'r2'     display R^2
%  'none'   display nothing
fp.corrcoefmode = 'r';


% Lin/log mode display: 'on' or 'off'
% (this tells whether Y or LOG(Y) is fitted)
fp.linlogdisp = 'on';


% Open the Array Editor with the fit coefficients after showfit:
% 'on' or 'off'. See Editcoeff for details.
fp.editcoeffmode = 'off';


% Array for the Fit Coefficients: 'row' or 'line':
% (this option is useful for copy-paste the coefficients in Excel)
fp.coeffarray = 'row';


% Automatic call of MakeVarFit after each fit: 'on' or 'off'
% Setting to 'on' creates on the Matlab workspace the variables
% associated to the fit parameters (the variables will be overwritten
% if they already exist!). See Makevarfit for details
fp.automakevarfit = 'off';


% Maximum length of the equation string to be displayed in
% the equation box (longer strings are truncated).
% Set 'maxlength = inf' for no truncation:
fp.maxlengtheq = 35;
%fp.maxlengtheq = inf;


% Location and size of the equation box (in normalized units):
fp.boxlocation = [0.15 0.81 0.3 0.1]; % this is the top-left corner


% Which data to fit when several curves are present in the figure:
% 'first' or 'last'.
fp.whichpickdata = 'first';



% Name of the polynom coefficients. Using a '_' (underscore)
% at the end of the coeff name allows for using the coefficient order
% as a subscript (latex syntax)
% (default: 'a_')
fp.polynom_coeffname = 'a_';


% Number of digits for the coefficient values (default = 5)
fp.numberofdigit = 5;

% fminsearch or ga options
fp.options = [];
%use ga?
fp.usega = false;

end

ezfit() had to be modified more.

function f = ezfit(varargin)
%EZFIT   Fit data with arbitrary fitting function
%   EZFIT(FUN) fits the active curve with the function FUN. See below for
%   the syntax of FUN. If FUN is not specified, 'linear' is used.
%
%   By default, the first curve in the active figure is fitted - see
%   FITPARAM to change this default behavior. To fit another curve, select
%   it before calling EZFIT.  If some data are selected by the "Data
%   Brushing" tool (only for Matlab >= 7.6), only those data are fitted.
%
%   EZFIT(X,Y,FUN) or EZFIT(Y,FUN) fits the data (X,Y) (or Y) using the
%   function FUN. X and Y must be vectors of equal length. If X is not
%   specified, X=[1, 2, 3...] is assumed.
%
%   EZFIT(X,Y,FUN), where X is a 1-by-N vector and Y is a 2-by-N matrix,
%   also specifies the weights for Y(2,:) (error bars). By default, when
%   Y is a 1-by-N vector, all the weights are 1.
%
%   Note that EZFIT only computes the coefficients, but does not display the
%   fit. Use SHOWFIT to display the fit.
%
%   The function string FUN can be:
%      - the name of a predefined fitting function (see below).
%      - the name of a user-defined fitting function (see EDITFIT).
%      - an equation, in the form 'y(x)=...', where 'x' represents the
%        X-data, and all the other variables are parameters to be fitted
%        ('a', 'x_0', 'tau', ...). Example: 'y(x)=a*sin(b*x)'. If the
%        left-hand-side 'y(x)' is not specified, 'x' is taken for the
%        X-Data. All the parameter names are accepted, except Matlab
%        reserved strings ('sin', 'pi', ...)
%
%   The predefined fitting functions are:
%      - linear             y = m * x
%      - affine or poly1    y = a*x + b
%      - poly{n}            y = a0 + a1 * x + ... + an * x^n
%      - power              y = c*x^n
%      - sin                y = a * sin (b * x)
%      - cos                y = a * cos (b * x)
%      - exp                y = a * exp (b * x)
%      - log                y = a * log (b * x)
%      - cngauss            y = exp(-x^2/(2*s^2))/(2*pi*s^2)^(1/2)
%      - cfgauss            y = a*exp(-x^2/(2*s^2))
%      - ngauss             y = exp(-(x-x0)^2/(2*s^2))/(2*pi*s^2)^(1/2)
%      - gauss              y = a*exp(-(x-x0)^2/(2*s^2))
%   'ngauss' is a 2-parameters normalized Gaussian, and 'gauss' is a
%   3-parameters non-normalized (free) Gaussian. 'cngauss' and 'cfgauss'
%   are centered normalized and centered free Gaussian, respectively.
%
%   EZFIT is based on Matlab's built-in FMINSEARCH function (Nelder-Mead
%   method), which performs an unconstrained nonlinear minimization of
%   the SSR (sum of squared residuals) with respect to the various
%   parameters. 
%
%   The correlation coefficient R is defined as SSreg / SStot, where
%      SSreg = sum ((y_fit - mean(y)).^2)   % regression sum of squares
%      SStot = sum ((y     - mean(y)).^2)   % total sum of squares
%   (see https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient)
%   (NB: definition of R changed in Ezyfit v2.44)
%
%   Nonlinear minimization requires starting guesses (or starting estimates)
%   for the fit parameters. By default, all the starting guesses are taken
%   as 1, or, when using predefined fits (e.g., exp, gauss, power...), the
%   starting guesses are determined depending on the range of the data to
%   be fitted. However, in most cases, values closer to the expected
%   result should be specified to "help" the convergence. It is sufficient
%   to choose values that have the correct sign and correct order of
%   magnitude, e.g. 0.01, 1, 100...
%
%   The starting guesses for the parameters may be specified in two ways:
%     - directly in the string FUN, after the fit definition:
%          'c0 + a*sin(pi*x/lambda); c0=1; a=0.1; lambda=100'
%          ('!' or '$' may also be used instead of ';').
%     - by specifying them as an additional input argument for EZFIT:
%          EZFIT(x,y,'c0 + a*sin(pi*x/lambda)',[0.1 1 100]);
%       Note that in this case the parameters must be ordered alphabetically.
%   Note that if both methods are used, only the starting guesses given in
%   the string FUN are considered.
%
%   By default, Y is fitted in linear mode. If you want to fit LOG(Y)
%   instead, you must specify the option 'log' to the string FUN, separated
%   by the symbol ';' or '$' or '!' (eg, FUN='a*x^n;log'). This is
%   specially useful to fit power laws with equally weighted points in a
%   log scale.  If nothing specified, the option 'lin' is used.
%
%   Example:
%        plotsample('power')
%   and compare the results of:
%        ezfit('power;lin')
%        ezfit('power;log')
%
%   F = EZFIT(...) also returns a structure F having the following fields:
%      - name       name of the fit
%      - eq         equation of the fit
%      - param      cell array of strings: names of the parameters
%      - m          values of the coefficients
%      - m0         initial guess for the coefficients
%      - r          correlation coefficient R (Pearson's correlation)
%      - fitmode    'lin' (y is fitted) or 'log' (log(y) is fitted) mode
%
%   This structure F can be further used with SHOWFIT, DISPEQFIT,
%   SHOWEQBOX, MAKEVARFIT and EDITCOEFF.
%
%   From F, you can get the values of the fitted parameters. If you want to
%   create in the current Matlab workspace the variables associated to
%   these parameters, use MAKEVARFIT (or set the option 'automakevarfit'
%   to 'on' in FITPARAM).
%
%   Examples:
%     plotsample damposc
%     f = ezfit('u(t) = c + u_a * sin(2*pi*t/T) * exp(-t/tau); T=5; tau=20');
%     showfit(f);
%     
%     plotsample poly2
%     [x,y] = pickdata;
%     f = ezfit(x, y, 'z(v) = poly3');
%     editcoeff(f);
%
%     plotsample poly2
%     f = ezfit('beta(z) = poly2');
%     showfit(f, 'fitcolor', 'red', 'fitlinewidth', 2);
%
%   See also SHOWFIT, PLOTSAMPLE, DISPEQFIT, EDITCOEFF,
%   FMINSEARCH, MAKEVARFIT, FITPARAM.


%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 2.53,  Date: 2016/04/28
%   This function is part of the EzyFit Toolbox


% History:
% 2005/05/12: v1.00, first version.
% 2005/05/20: v1.10, Use 'eval', with generic functions
% 2005/05/21: v1.11, added the 'lin','log' options.
% 2005/05/24: v1.12, option 'log' by default for 'power' and 'exp'.
% 2005/07/27: v1.13, cosmetics.
% 2005/09/03: v1.14, check arg.
% 2005/10/07: v1.15, gaussian fits added (ngauss and fgauss, centered/not)
% 2005/10/20: v1.16, help text changed.
% 2005/10/31: v1.20, also returns R. cste and poly{n} fits added. Initial
%                    guess defined within the fitting function string. The
%                    order of the output parameters is changed.
% 2005/11/05: v1.21, evaluate strings for initial guess in FUN
% 2005/12/06: v1.22, opens a dialog box if the polynomial order is not
%                    specified.
% 2006/01/13: v1.24, check for negative data in log mode
% 2006/01/19: v1.25, bug fixed from 1.24
% 2006/01/25: v1.26, check the matlab version
% 2006/02/08: v2.00, new syntax. The output argument is now a fit
%                    structure, and the fitting equation string accepts
%                    arbitrary parameter names.
% 2006/02/14: v2.10, lhs 'y(x)=...' now accepted. Now case sensitive
% 2006/03/07: v2.11, bug fixed from v1.25
% 2006/03/09: v2.20, weighted chi-square criterion (ie: error bars accepted
%                    for y) - undocumented
% 2006/09/04: v2.21, '!' and '$' may be used instead of ';' in the FUN
%                    string: this allows to pass the argument like:
%                    fit power!log
% 2006/09/28: v2.22, 'x^1' replaced by 'x' for 1st order polynomial
% 2006/10/18: v2.30, input raw and column vectors accepted. accepts
%                    additional input parameters to change the default settings.
% 2007/04/17: v2.31, help text improved; standard error messages
% 2007/05/16: v2.40, weigthed fits (v2.20) now documented, with bug fixed.
% 2007/08/18: v2.41, help text improved.
% 2007/09/17: v2.50, guess the initial guess m0 for predefined fits
% 2011/11/04: v2.51, bug fixed for fits in log coordinates
% 2012/06/28: v2.52, check version number removed (unstable)
% 2016/04/28: v2.53, definitions of SSE and SSR exchanged (thanks Yoel!!)

% 2022/09/14: * add support for fminsearch options, especially so that the
%             number of iterations can be changed. Add warnings about
%             using many iterations (so that the user can adjust the
%             iteration limit more intelligently)
%             * Add support for using ga -- it often finds a solution much
%             much faster in cases where fminsearch would wander off in the
%             wrong direction.
%             * code cleanup to remove stray commas and semi-colons



% new v1.26, changed v2.11, removed v2.52
% if str2double(version('-release'))<14
%     error('Ezyfit:ezfit:compatibility','EzyFit requires Matlab 7 or higher.');
% end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fit parameters:  (new v2.30)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% loads the default fit parameters:
try
    fp=fitparam;
catch
    error('Ezyfit:ezfit:fitparamNotFound','''fitparam.m'' file not found.');
end


% change the default values of the fit parameters according to the
% additional input arguments:
for nopt=1:(nargin-1)
    if any(strcmp(varargin{nopt},fieldnames(fp))) % if the option string is one of the fit parameter
        fp.(varargin{nopt}) = varargin{nopt+1};
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input arguments:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==0 % if no input argument: do a linear fit of the current curve
    [x,y,h] = pickdata(fp);
    inputstr='linear';
elseif nargin==1
    if ischar(varargin{1})    %    EZFIT('a*x+b')
        inputstr=varargin{1};
        [x,y,h] = pickdata(fp);
    else                       %    EZFIT(y)
        y=varargin{1};
        x=1:length(y);
        inputstr='linear';
    end
else   % 2 or more input arguments
    if ~isnumeric(varargin{1})           % EZFIT('fun',...)
        inputstr=varargin{1};
        [x,y,h] = pickdata(fp);
        if isnumeric(varargin{2})    % EZFIT('fun',m0,...)
            m0=varargin{2};
        end
    elseif (isnumeric(varargin{1}) && ~isnumeric(varargin{2}))       % EZFIT(y,'a*x+b',...)
        [y,inputstr]=deal(varargin{1:2});
        x=1:length(y);
        if nargin>2
            if isnumeric(varargin{3})   % EZFIT(y,'a*x+b',m0,...)    
                m0=varargin{3};       
            end
        end
    elseif (isnumeric(varargin{1}) && isnumeric(varargin{2}))   % EZFIT(x,y,...)
        [x,y]=deal(varargin{1:2});      
        if nargin>2
            if ischar(varargin{3})
                inputstr=varargin{3};      % EZFIT(x,y,'fun',...)
                if nargin>3
                    if isnumeric(varargin{4})      % EZFIT(x,y,'fun',m0,...)
                        m0=varargin{4};
                    end
                end
            else
                error('Ezyfit:ezfit:syntaxError','Syntax error. 3rd paramater of EZFIT should be a string.');
            end
        end
    end
end               



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some checks about x and y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% turn all input vectors into row (1xN) vectors (new v2.30):
if size(x,1)>size(x,2)
    x=x';
end
if size(y,1)>size(y,2)
    y=y';
end


% check for error bars for y (new v2.20, fixed 2.40)
if size(y,1)>1
    dy = y(2,:);  % second line = error bars (1/weight)
    y = y(1,:);   % first line = data
else
    dy = ones(1,length(y));    % default error bars: 1
end

if length(x)~=length(y)
    error('Ezyfit:ezfit:dimagree','X and Y dimensions must agree.');
end







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% processing of the string FUN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% cleans the input string:
inputstr = strrep(inputstr,' ','');
inputstr = strrep(inputstr,'!',';');
inputstr = strrep(inputstr,'$',';');
inputstr = strrep([inputstr ';'],';;',';'); % ensures that fun is terminated by a ';'.

% the name of the fit is by default the first part of the input string
p=findstr(inputstr,';'); p=p(1);
f.name = inputstr(1:(p-1));

% separates the first part (fitting function itself) and the remainder:
p = strfind(inputstr,';'); p=p(1);
fun = inputstr(1:(p-1));
remfun = inputstr((p+1):end);


% search for predefined fit or user-defined fit
usepredefinedfit = '';
[defaultfit, userfit] = loadfit;
for i=1:length(defaultfit)
    if strcmp(fun, defaultfit(i).name)
        fun = defaultfit(i).eq;
        usepredefinedfit = defaultfit(i).name;   % new v2.50
    end
end
for i=1:length(userfit)
    if strcmp(fun, userfit(i).name)
        fun = userfit(i).eq;
    end
end


% separates again the first part (fitting function itself) and the remainder:
% (because the predefined/user-defined part may itself contain ';')
fun = [strrep(fun,' ','') ';' remfun];
p=strfind(fun,';'); p=p(1);
remfun = fun((p+1):end);
fun = fun(1:(p-1));


% recognize if a lhs is present
peq = strfind(fun,'=');
if ~isempty(peq)
    lhs = fun(1:(peq-1));    % left-hand side
    rhs = fun((peq+1):end);  % right-hand side
else
    lhs = '';
    rhs = fun;
end


% process the lhs (if present)
if ~isempty(lhs)
    pob = strfind(lhs,'('); % position of opening bracket
    pcb = strfind(lhs,')'); % position of closing bracket
    if ~isempty(pob)
        if pob==1
            f.yvar = 'y';
        else
            f.yvar = lhs(1:(pob-1));
        end
        f.xvar = lhs((pob+1):(pcb-1));
    else
        f.yvar = lhs;
        f.xvar = 'x';
    end
else % if no lhs present:
    f.xvar='x';
    f.yvar='y';
end


% process the 'poly' (polynomial fit) in the rhs
if strfind(rhs,'poly')   % polynomial fit:
    order=str2double(rhs(5:end));
    if isempty(order) % new v1.22
        str_ord=inputdlg('Order of the polynomial fit','Polynomial order',1,{'2'});
        if ~isempty(str_ord)
            order=str2double(str_ord{1});
            f.name=['poly' str_ord{1}];
        else
            clear f;
            return;
        end
    end
    if order>20
        error('Ezyfit:ezfit:invalidPolynomDegree','Invalid polynom degree.');
    end
    rhs = [fp.polynom_coeffname '0'];
    for i=1:order
        if i>1
            rhs = [rhs '+' fp.polynom_coeffname num2str(i) '*' f.xvar '^' num2str(i)];
        else
            rhs = [rhs '+' fp.polynom_coeffname num2str(i) '*' f.xvar];   % new v2.22
        end
    end
end


% search for option 'lin' or 'log':
% (if several are present, use the last one)
% (if none is present, check the Y-scale of the current figure)
% (if no figure present, take 'lin').
fun = strrep([rhs ';' remfun], ';;', ';');
f.fitmode='';
plin=strfind(fun,';lin');
if ~isempty(plin); plin=plin(end); else; plin=1; end
plog=strfind(fun,';log');
if ~isempty(plog); plog=plog(end); else; plog=1; end
plast=max(plin,plog);
if plast==1 % if nothing specified
    f.fitmode='lin'; % use 'lin'
else
    f.fitmode=fun((plast+1):(plast+3));
end
fun=strrep(fun,';lin','');
fun=strrep(fun,';log','');


% check for negative data in log mode (new v1.23, fixed 1.24):
if (strcmp(f.fitmode,'log') && sum(y<=0)>0)
    disp('Warning: Zero or negative data ignored');
    nonzero=find(y>0);
    x=x(nonzero);
    y=y(nonzero);
end


% separates again the first part (fitting function itself) and the
% remainder:
p = strfind(fun,';'); p=p(1);
f.eq = fun(1:(p-1));
remfun = fun((p+1):end);


 % convert the equation in the matlab syntax (parameters named m(1),
 % m(2)...)
[eqml,param] = eq2ml(f.eq, f.xvar);
maxm = length(param); % number of parameters
if maxm==0  % new v2.20
    error('Ezyfit:ezfit:noParameter','No parameter to be fitted.');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% processing the initial guess
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial guess for the m(i) by default (all m(i)=1):
if ~exist('m0','var')
    m0=ones(1,maxm);
    switch usepredefinedfit   % new v2.50 : "magic" initial guesses
        case 'linear'   % m*x
            m0(1) = (y(end)-y(1))/(x(end)-x(1));
        case 'affine'   % a*x+b
            m0(1) = (y(end)-y(1))/(x(end)-x(1));  % a
            m0(2) = mean(y-m0(1)*x);              % b
        case 'affineshift'   % a*(x-b)
            m0(1) = (y(end)-y(1))/(x(end)-x(1));  % a
            m0(2) = mean(x-y/m0(1));              % b
        case 'power'    % a*x^n;
            % if both x and y are of constant sign:
            if sum(diff(sign(x)))==0 || sum(diff(sign(y)))==0
                m0(2) = log(y(end)/y(1)) / log(x(end)/x(1));  % n
                m0(1) = mean(y./(x.^m0(2)));                  % a
            end
        case 'powerc'    % a*x^n+c;
            % if both x and y are of constant sign:
            if sum(diff(sign(x)))==0 || sum(diff(sign(y)))==0
                m0(3) = log(y(end)/y(1)) / log(x(end)/x(1));  % n
                m0(1) = mean(y./(x.^m0(3)));                  % a
            end
            m0(2) = mean(y);  % c
        case 'powershift'    % a*(x+b)^n;
            % if both x and y are of constant sign:
            m0(2) = mean(x);                              % b
            m0(3) = log(y(end)/y(1)) / log((x(end)+m0(2))/(x(1)+m0(2)));  % n
            m0(1) = mean(y./((x+m0(2)).^m0(3)));                  % a           
        case 'exp'    % a*exp(b*x)
            m0(2) = (log(y(end)/y(1))) / (x(end)-x(1));  % b
            m0(1) = mean(y./exp(m0(2)*x));               % a
        case 'expdiv'    % a*exp(x/b)
            m0(2) = (x(end)-x(1)) / (log(y(end)/y(1)));  % b
            m0(1) = mean(y./exp(x/m0(2)));               % a
        case 'explim'    % a*(1-exp(-x/b))
            m0(1) = max(y);                              % a
            m0(2) = mean(-x./(log(1-y/m0(1))));          % b
        case 'expc'    % a*exp(b*x)+c
            m0(3) = mean(y);                             % c
            m0(2) = (log((y(end)-m0(3))/(y(1)-m0(3)))) / (x(end)-x(1));  % b
            m0(1) = mean((y-m0(3))./exp(m0(2)*x));               % a
        case 'log'    % a*log(b*x)
            m0(1) = (y(end)-y(1))/(log(x(end)/x(1)));    % a
            m0(2) = mean(exp(y/m0(1))./x);               % b            
        case 'logc'    % a*log(x)+b
            m0(1) = (y(end)-y(1))/(log(x(end)/x(1)));    % a
            m0(2) = mean(y - m0(1)*log(x));              % b     
        case {'sin','cos'}    % a*sin(b*x), a*cos(b*x)
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
        case {'sinc','cosc'}    % a*sin(b*x)+c, a*cos(b*x)+c
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
            m0(3) = mean(y);           % c
        case 'sinphi'    % a*sin(b*x+phi)
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
            m0(3) = 1;                 % phi
        case 'sinphic'    % a*sin(b*x+phi)+c
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
            m0(3) = mean(y);           % c
            m0(4) = 1;                 % phi
        case 'cngauss'
            m0(1) = (mean(x.^2.*y)/mean(y))^(1/2);   % sigma
        case 'cfgauss'   % a*exp(-(x^2)/(2*sigma^2))
            m0(1) = max(y);      % a
            m0(2) = (mean(x.^2.*y)/mean(y))^(1/2);   % sigma
        case 'ngauss'    % exp(-((x-x_0)^2)/(2*sigma^2))/(2*pi*sigma^2)^(1/2)           
            m0(2) = mean(x.*y)/mean(y);             % x_0
            m0(1) = (mean((x-m0(2)).^2.*y)/mean(y))^(1/2);  % sigma
        case {'fgauss','gauss'}    % a*exp(-((x-x_0)^2)/(2*sigma^2))
            m0(1) = max(y);     % a
            m0(3) = mean(x.*y)/mean(y);             % x_0                   
            m0(2) = (mean((x-m0(3)).^2.*y)/mean(y))^(1/2);  % sigma
    end
    % if an initial guess is 0 or infinity or has imaginary part, set it to 1
    for np=1:length(m0)
        if m0(np)==0 || isinf(m0(np)) || imag(m0(np))~=0
            m0(np) = 1;
        end
    end
end
    

% search for initial guess defined into remfun
remfun = strrep([remfun ';'], ';;', ';'); % adds a final ';'
while strfind(remfun,'=')
    peq=strfind(remfun,'='); peq=peq(1);
    pc=strfind(remfun,';'); pc=pc(1);
    if pc>peq   % if ';' after '='
        par=remfun(1:(peq-1));
        for i=1:maxm
            if strcmp(param{i},par)
                m0(i)=eval(remfun((peq+1):(pc-1)));
            end
        end
    else
        error('Ezyfit:ezfit:syntaxError','Invalid syntax');
    end
    remfun=remfun((pc+1):end); % removes the processed i.g. and loops back
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fitting itself
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isfield(fp, 'options')
    fp.options = [];
    fp.usega = false;
end
if ~isfield(fp, 'usega');
    fp.usega = false;
end

switch f.fitmode
    case 'lin'
        x__ref = x;
        try
            if ~fp.usega
                [m, fval, exitflag, output] = fminsearch(@fitlin, m0, fp.options);
            else
                [m, fval, exitflag, output] = ga(@fitlin, length(m0), [], [], [], [], [], [], [],fp.options);
            end
            if isfield(output, 'funccount')
                fc = output.funccount;
            elseif isfield(output, 'funcCount')
                fc = output.funcCount;
            else
                fc = -inf;
            end
            if exitflag == 0
                fprintf(2, 'ran out of iterations at %d evaluations\n', fc);
            elseif fc > 10000
                fprintf(1, 'FYI, used %d function evaluations\n', fc);
            end
        catch ME
            error('Ezyfit:ezfit:searchError','Fit: error during the search procedure');
        end
        y_fit = eval(eqml);
        % new definition 2016:
        ssreg = sum(abs(y_fit-mean(y)).^2);
        sstot = sum(abs(y-mean(y)).^2);
        f.r = ssreg/sstot;
    case 'log'
        x__ref = x;
        try
            if ~fp.usega
                [m, fval, exitflag, output] = fminsearch(@fitlog, m0, fp.options);
            else
                [m, fval, exitflag, output] = ga(@fitlin, length(m0), [], [], [], [], [], [], [],fp.options);
            end
            if isfield(output, 'funccount')
                fc = output.funccount;
            elseif isfield(output, 'funcCount')
                fc = output.funcCount;
            else
                fc = -inf;
            end

            if exitflag == 0
                fprintf(2, 'ran out of iterations at %d evaluations\n', fc);
            elseif fc > 10000
                fprintf(1, 'FYI, used %d function evaluations\n', fc);
            end
        catch
            error('Ezyfit:ezfit:fminsearchError','Fit: error during the fminsearch procedure');
        end
        y_fit = eval(eqml);
        % new definition 2016:
        ssreg = sum(abs(log(y_fit)-mean(log(y))).^2);
        sstot = sum(abs(log(y)-mean(log(y))).^2);
        f.r = ssreg/sstot;
    otherwise
        error('Ezyfit:ezfit:unknownMode''Unknown fit mode');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% outputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% fills the output structure:
f.param = param;
f.m = m;
f.m0 = m0;
f.x = x;
f.y = y;
if sum(dy-1) % store the error bars only if defined
    f.dy = dy;
end
if exist('h','var'), f.hdata=h; end  % handle to the data

% stores the fit in the variable 'lastfit' in the 'base' workspace:
assignin('base','lastfit',f);

if strcmp(fp.automakevarfit,'on')
    makevarfit;
end

% ending displays (if no output argument):
if ~nargout
    if strcmp(fp.dispeqmode,'on') % new v2.30
        dispeqfit(f,fp);
    end
    clear f;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of the main function EZFIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Nested functions that evaluate the fit for prescribed parameters m(i),
% and return the chi2 (sum of the squared difference between the input
% curve and the fit), in lin or log mode:

    function chi2 = fitlin(m)
        y_fit = eval(eqml);
        chi2 = sum(((y_fit - y).^2)./(dy.^2));
    end

    function chi2 = fitlog(m)
        y_fit = eval(eqml);
        chi2 = sum((log(y_fit)-log(y)).^2);
    end
end

[HYZ]

thanks a lot! your modified codes work without showing error. I realised sigma values show everytime I ran the modified ezfit code as the number of function evaluations is also changed. Sigma changes in decimal points so it is not a lot. I guess the changes in function evaluations are unavoidable.

[Walter Roberson]

When you use fminsearch() then there is no randomness used at all.

When you use ga() then it uses a lot of random numbers.

So if you need to reproduce the results, then use rng() to set a random number generator before the call to ga()

Also note that ga() is not as good at finding the ultimate "best" location -- but it is pretty good at finding a "close" location.

You might want to consider editting the optimoptions() call to configure 'HybridFcn', 'fminsearch' to cause it to run an fminsearch pass after ga() in order to get even closer to the best value. https://www.mathworks.com/help/gads/genetic-algorithm-options.html#f17822

[Sim]

A great update! But, just a question.... Who added the following notes related to the version update?

% 2022/09/14: * add support for fminsearch options, especially so that the
%             number of iterations can be changed. Add warnings about
%             using many iterations (so that the user can adjust the
%             iteration limit more intelligently)
%             * Add support for using ga -- it often finds a solution much
%             much faster in cases where fminsearch would wander off in the
%             wrong direction.
%             * code cleanup to remove stray commas and semi-colons

Was it you @Walter Roberson?

In addition, I do not see this updated version on the Moisy and Ezfit-related websites, which look like unchanged since the "Version 2.45 23-Jan-2017" (http://www.fast.u-psud.fr/ezyfit/html/ezyfit.html):

• http://www.fast.u-psud.fr/~moisy/ml/
• http://www.fast.u-psud.fr/ezyfit/
• http://www.fast.u-psud.fr/ezyfit/html/ezyfit.html
• http://www.fast.u-psud.fr/ezyfit/html/ezyfit_releasenotes.html

Was there any contact between @Walter Roberson and Frederic Moisy (I found a contact here, even though it might be not working anymore: http://www.fast.u-psud.fr/ezyfit/html/ezyfit_knownbugs.html) in order to update Ezyfit in the corresponding websites where it would be possible to download it?

I think that this 2022 update would be very useful to many people! :-)

[Walter Roberson]

@Sim

Yes, those are my changes, and No, I did not contact Frederic Moisy

[Sim]

Thanks for the information :-)