9e68520052c39b8d464da1b33ef931f703746565
[scilab.git] / scilab / modules / optimization / macros / fit_dat.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 //
4 // This file must be used under the terms of the CeCILL.
5 // This source file is licensed as described in the file COPYING, which
6 // you should have received as part of this distribution.  The terms
7 // are also available at
8 // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
9 //
10 function [p,err]=fit_dat(G,p0,Z,W,pmin,pmax,DG)
11     //
12     //         [p,err]=fit_dat(G,p0,Z [,W] [,pmin,pmax] [,DG])
13     //
14     //         Function used for fitting data to a model.
15     // For a given function G(p,z), this function finds the best vector
16     // of parameters p for approximating G(p,z_i)=0 for a set of measurement
17     // vectors z_i. Vector p is found by minimizing
18     //    G(p,z_1)'WG(p,z_1)+G(p,z_2)'WG(p,z_2)+...+G(p,z_n)'WG(p,z_n)
19     //
20     //      G: Scilab function (e=G(p,z), e: nex1, p: npx1, z: nzx1)
21     //     p0: initial guess (size npx1)
22     //      Z: matrix [z_1,z_2,...z_n] where z_i (nzx1) is the ith measurement
23     //      W: weighting matrix of size nexne (optional; default 1)
24     //   pmin: lower bound on p (optional; size npx1)
25     //   pmax: upper bound on p (optional; size npx1)
26     //     DG: partial of G wrt p (optional; S=DG(p,z), S: nexnp)
27     //
28     //                     Examples
29     //
30     //deff('y=FF(x)','y=a*(x-b)+c*x.*x')
31     //X=[];Y=[];
32     //a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand-.5)];X=[X,x];end
33     //Z=[Y;X];
34     //deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
35     //[p,err]=fit_dat(G,[3;5;10],Z)
36     //xset('window',0)
37     //clf();
38     //plot2d(X',Y',-1)
39     //plot2d(X',FF(X)',5,'002')
40     //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
41     //
42     //a=34;b=12;c=14;
43     //deff('s=DG(p,z)','y=z(1),x=z(2),s=-[x-p(2),-p(1),x*x]')
44     //[p,err]=fit_dat(G,[3;5;10],Z,DG)
45     //xset('window',1)
46     //clf();
47     //plot2d(X',Y',-1)
48     //plot2d(X',FF(X)',5,'002')
49     //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
50     //
51     //
52     warnobsolete("datafit","5.5.0");
53
54     [lhs,rhs]=argn(0)
55     boun=%f
56     if rhs==3 then
57         W=1;GR=%f  //Jacobian not provided
58     elseif rhs==4 then
59         if type(W)==11|type(W)==13 then
60             DG=W;W=1;GR=%t  //Jacobian provided
61         else
62             GR=%f  //Jacobian not provided
63         end
64     elseif rhs==5 then
65         if type(pmin)==1 then
66             pmax=pmin,pmin=W,W=1,GR=%f,boun=%t
67         else
68             DG=pmin;GR=%t
69         end
70     elseif rhs==6 then
71         boun=%t
72         if type(pmax)==1 then
73             GR=%f
74         else
75             DG=pmax,pmax=pmin,pmin=W,W=1,GR=%t
76         end
77     elseif    rhs==7 then
78         GR=%t;
79     else
80         error(msprintf(gettext("%s: Wrong number of input argument(s).\n"),"fit_dat"));
81     end
82     if size(W,1)~=size(W,2) then
83         error(msprintf(gettext("%s: Weighting matrix must be square.\n"),"fit_dat"));
84     end
85     nz=size(Z,2);mz=size(p0,"*");
86     deff("e=GG(p,Z)","e=0;for i=1:nz,e=e+G(p,Z(:,i))''*W*G(p,Z(:,i)),end")
87     if ~GR then
88         deff("g=DGG(ps,Z)",["g=[]";
89         "nrm=norm(ps);if nrm<1 then nrm=1,end";
90         "[d,nu]=colcomp(rand(mz,mz));d=10*%eps*nrm*d";
91         "for v=d,e=GG(ps+v),g=[g,e-f],end;g=g/d;"])
92     else
93         deff("g=DGG(ps,Z)","g=0;for i=1:nz,g=g+2*DG(ps,Z(:,i))''*W*G(ps,Z(:,i)),end")
94     end
95     deff("[f,g,ind]=costf(p,ind)",["f=GG(p,Z);";
96     "g=DGG(p,Z);"])
97
98     if ~boun then
99         [err,p]=optim(costf,p0)
100     else
101         if or(p0<pmin)|or(p0>pmax) then
102             error(msprintf(gettext("%s: Initial guess not feasible.\n"),"fit_dat"));
103         end
104         [err,p]=optim(costf,"b",pmin,pmax,p0)
105     end
106
107 endfunction