diff options
Diffstat (limited to 'macros/lsqnonlin.sci')
| -rw-r--r-- | macros/lsqnonlin.sci | 317 |
1 files changed, 317 insertions, 0 deletions
diff --git a/macros/lsqnonlin.sci b/macros/lsqnonlin.sci new file mode 100644 index 0000000..50a1a49 --- /dev/null +++ b/macros/lsqnonlin.sci @@ -0,0 +1,317 @@ +// Copyright (C) 2015 - IIT Bombay - FOSSEE +// +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Author: Harpreet Singh +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +function [xopt,resnorm,residual,exitflag,output,lambda,gradient] = lsqnonlin (varargin) + // Solves a non linear data fitting problems. + // + // Calling Sequence + // xopt = lsqnonlin(fun,x0) + // xopt = lsqnonlin(fun,x0,lb,ub) + // xopt = lsqnonlin(fun,x0,lb,ub,options) + // [xopt,resnorm] = lsqnonlin( ... ) + // [xopt,resnorm,residual] = lsqnonlin( ... ) + // [xopt,resnorm,residual,exitflag] = lsqnonlin( ... ) + // [xopt,resnorm,residual,exitflag,output,lambda,gradient] = lsqnonlin( ... ) + // + // Parameters + // fun : a function, representing the objective function and gradient (if given) of the problem + // x0 : a vector of double, contains initial guess of variables. + // lb : a vector of double, contains lower bounds of the variables. + // ub : a vector of double, contains upper bounds of the variables. + // options : a list containing the parameters to be set. + // xopt : a vector of double, the computed solution of the optimization problem. + // resnorm : a double, objective value returned as the scalar value i.e. sum(fun(x).^2). + // residual : a vector of double, solution of objective function i.e. fun(x). + // exitflag : The exit status. See below for details. + // output : The structure consist of statistics about the optimization. See below for details. + // lambda : The structure consist of the Lagrange multipliers at the solution of problem. See below for details. + // gradient : a vector of doubles, containing the Objective's gradient of the solution. + // + // Description + // Search the minimum of a constrained non-linear least square problem specified by : + // + // <latex> + // \begin{eqnarray} + // &\mbox{min}_{x} + // & (f_1(x)^2 + f_2(x)^2 + ... + f_n(x)^2) \\ + // & lb \leq x \leq ub \\ + // \end{eqnarray} + // </latex> + // + // The routine calls fmincon which calls Ipopt for solving the non-linear least square problem, Ipopt is a library written in C++. + // + // The options allows the user to set various parameters of the Optimization problem. + // It should be defined as type "list" and contains the following fields. + // <itemizedlist> + // <listitem>Syntax : options= list("MaxIter", [---], "CpuTime", [---],"GradObj", "on");</listitem> + // <listitem>MaxIter : a Scalar, containing the Maximum Number of Iteration that the solver should take.</listitem> + // <listitem>CpuTime : a Scalar, containing the Maximum amount of CPU Time that the solver should take.</listitem> + // <listitem>GradObj : a string, representing the gradient function is on or off.</listitem> + // <listitem>Default Values : options = list("MaxIter", [3000], "CpuTime", [600], "GradObj", "off");</listitem> + // </itemizedlist> + // + // The exitflag allows to know the status of the optimization which is given back by Ipopt. + // <itemizedlist> + // <listitem>exitflag=0 : Optimal Solution Found </listitem> + // <listitem>exitflag=1 : Maximum Number of Iterations Exceeded. Output may not be optimal.</listitem> + // <listitem>exitflag=2 : Maximum CPU Time exceeded. Output may not be optimal.</listitem> + // <listitem>exitflag=3 : Stop at Tiny Step.</listitem> + // <listitem>exitflag=4 : Solved To Acceptable Level.</listitem> + // <listitem>exitflag=5 : Converged to a point of local infeasibility.</listitem> + // </itemizedlist> + // + // For more details on exitflag see the ipopt documentation, go to http://www.coin-or.org/Ipopt/documentation/ + // + // The output data structure contains detailed informations about the optimization process. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>output.iterations: The number of iterations performed during the search</listitem> + // <listitem>output.constrviolation: The max-norm of the constraint violation.</listitem> + // </itemizedlist> + // + // The lambda data structure contains the Lagrange multipliers at the end + // of optimization. In the current version the values are returned only when the the solution is optimal. + // It has type "struct" and contains the following fields. + // <itemizedlist> + // <listitem>lambda.lower: The Lagrange multipliers for the lower bound constraints.</listitem> + // <listitem>lambda.upper: The Lagrange multipliers for the upper bound constraints.</listitem> + // </itemizedlist> + // + // Examples + // //A simple non-linear least square example taken from leastsq default present in scilab + // function y=yth(t, x) + // y = x(1)*exp(-x(2)*t) + // endfunction + // // we have the m measures (ti, yi): + // m = 10; + // tm = [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5]'; + // ym = [0.79, 0.59, 0.47, 0.36, 0.29, 0.23, 0.17, 0.15, 0.12, 0.08]'; + // // measure weights (here all equal to 1...) + // wm = ones(m,1); + // // and we want to find the parameters x such that the model fits the given + // // data in the least square sense: + // // + // // minimize f(x) = sum_i wm(i)^2 ( yth(tm(i),x) - ym(i) )^2 + // // initial parameters guess + // x0 = [1.5 ; 0.8]; + // // in the first examples, we define the function fun and dfun + // // in scilab language + // function y=myfun(x, tm, ym, wm) + // y = wm.*( yth(tm, x) - ym ) + // endfunction + // // the simplest call + // [xopt,resnorm,residual,exitflag,output,lambda,gradient] = lsqnonlin(myfun,x0) + // // Press ENTER to continue + // + // Examples + // //A basic example taken from leastsq default present in scilab with gradient + // function y=yth(t, x) + // y = x(1)*exp(-x(2)*t) + // endfunction + // // we have the m measures (ti, yi): + // m = 10; + // tm = [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5]'; + // ym = [0.79, 0.59, 0.47, 0.36, 0.29, 0.23, 0.17, 0.15, 0.12, 0.08]'; + // // measure weights (here all equal to 1...) + // wm = ones(m,1); + // // and we want to find the parameters x such that the model fits the given + // // data in the least square sense: + // // + // // minimize f(x) = sum_i wm(i)^2 ( yth(tm(i),x) - ym(i) )^2 + // // initial parameters guess + // x0 = [1.5 ; 0.8]; + // // in the first examples, we define the function fun and dfun + // // in scilab language + // function [y,dy]=myfun(x, tm, ym, wm) + // y = wm.*( yth(tm, x) - ym ) + // v = wm.*exp(-x(2)*tm) + // dy = [v , -x(1)*tm.*v] + // endfunction + // options = list("GradObj", "on") + // [xopt,resnorm,residual,exitflag,output,lambda,gradient] = lsqnonlin(myfun,x0,[],[],options) + // Authors + // Harpreet Singh + + + //To check the number of input and output argument + [lhs , rhs] = argn(); + + //To check the number of argument given by user + if ( rhs < 2 | rhs == 3 | rhs > 5 ) then + errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be in the set of [2 4 5]"), "lsqnonlin", rhs); + error(errmsg) + end + +// Initializing all the values to empty matrix + + + _fun = varargin(1); + x0 = varargin(2); + nbVar = size(x0,'*'); + + if(nbVar == 0) then + errmsg = msprintf(gettext("%s: Cannot determine the number of variables because input initial guess is empty"), "lsqcurvefit"); + error(errmsg); + end + + if (size(x0,2)~=1) then + errmsg = msprintf(gettext("%s: Initial Guess should be a column matrix"), "lsqcurvefit"); + error(errmsg); + end + + if ( rhs<3 ) then + lb = repmat(-%inf,nbVar,1); + ub = repmat(%inf,nbVar,1); + else + lb = varargin(3); + ub = varargin(4); + end + + if ( rhs<7 | size(varargin(5)) ==0 ) then + param = list(); + else + param =varargin(5); + end + + if (size(lb,2)==0) then + lb = repmat(-%inf,nbVar,1); + end + + if (size(ub,2)==0) then + ub = repmat(%inf,nbVar,1); + end + + if (type(param) ~= 15) then + errmsg = msprintf(gettext("%s: param should be a list "), "lsqnonlin"); + error(errmsg); + end + + //Check type of variables + Checktype("lsqnonlin", _fun, "fun", 1, "function") + Checktype("lsqnonlin", x0, "x0", 2, "constant") + Checktype("lsqnonlin", lb, "lb", 3, "constant") + Checktype("lsqnonlin", ub, "ub", 4, "constant") + + if (modulo(size(param),2)) then + errmsg = msprintf(gettext("%s: Size of parameters should be even"), "lsqnonlin"); + error(errmsg); + end + + options = list( "MaxIter" , [3000], ... + "CpuTime" , [600], ... + "GradObj" , ["off"]); + + for i = 1:(size(param))/2 + + select convstr(param(2*i-1),'l') + case "maxiter" then + options(2*i) = param(2*i); + Checktype("lsqcurvefit", param(2*i), "maxiter", i, "constant") + case "cputime" then + options(2*i) = param(2*i); + Checktype("lsqcurvefit", param(2*i), "cputime", i, "constant") + case "gradobj" then + if (convstr(param(2*i))=="on") then + options(2*i) = "on"; + else + errmsg = msprintf(gettext("%s: Unrecognized String [%s] entered for gradobj."), "lsqnonlin",param(2*i)); + error(errmsg); + end + else + errmsg = msprintf(gettext("%s: Unrecognized parameter name ''%s''."), "lsqnonlin", param(2*i-1)); + error(errmsg) + end + end + + // Check if the user gives row vector + // and Changing it to a column matrix + + if (size(lb,2)== [nbVar]) then + lb = lb(:); + end + + if (size(ub,2)== [nbVar]) then + ub = ub(:); + end + + //To check the match between fun (1st Parameter) and x0 (2nd Parameter) + if(execstr('init=_fun(x0)','errcatch')==21) then + errmsg = msprintf(gettext("%s: Objective function and x0 did not match"), "fmincon"); + error(errmsg); + end + + ierr = execstr('init=_fun(x0)', "errcatch") + if ierr <> 0 then + lamsg = lasterror(); + lclmsg = "%s: Error while evaluating the function: ""%s""\n"; + error(msprintf(gettext(lclmsg), "lsqnonlin", lamsg)); + end + + //Check the size of Lower Bound which should be equal to the number of variables + if ( size(lb,1) ~= nbVar) then + errmsg = msprintf(gettext("%s: The Lower Bound is not equal to the number of variables"), "lsqnonlin"); + error(errmsg); + end + + //Check the size of Upper Bound which should equal to the number of variables + if ( size(ub,1) ~= nbVar) then + errmsg = msprintf(gettext("%s: The Upper Bound is not equal to the number of variables"), "lsqnonlin"); + error(errmsg); + end + + //Check if the user gives a matrix instead of a vector + + if (size(lb,1)~=1)& (size(lb,2)~=1) then + errmsg = msprintf(gettext("%s: Lower Bound should be a vector"), "lsqnonlin"); + error(errmsg); + end + + if (size(ub,1)~=1)& (size(ub,2)~=1) then + errmsg = msprintf(gettext("%s: Upper Bound should be a vector"), "lsqnonlin"); + error(errmsg); + end + + for i = 1:nbVar + if(lb(i)>ub(i)) + errmsg = msprintf(gettext("%s: Problem has inconsistent variable bounds"), "lsqnonlin"); + error(errmsg); + end + end + + //Converting it into fmincon Problem + function [y] = __fun(x) + [xVal] = _fun(x); + y = sum(xVal.^2); + endfunction + + if (options(6) == "on") + ierr = execstr('[initx initdx]=_fun(x0)', "errcatch") + if ierr <> 0 then + lamsg = lasterror(); + lclmsg = "%s: Error while evaluating the function: ""%s""\n"; + error(msprintf(gettext(lclmsg), "lsqnonlin", lamsg)); + end + function [dy] = __fGrad(x) + [x_user,dx_user] = _fun(x) + dy = 2*dx_user*(x_user'); + endfunction + options(6) = __fGrad; + end + + + [xopt,fopt,exitflag,output,lambda_fmincon,gradient] = fmincon(__fun,x0,[],[],[],[],lb,ub,[],options); + + lambda = struct("lower", lambda_fmincon.lower,"upper",lambda_fmincon.upper); + + resnorm = sum(_fun(xopt).^2); + residual = _fun(xopt); + +endfunction |