9bd398fac698eb7cd4b9ed36233400d3ba6452be
[scilab.git] / scilab / modules / optimization / macros / datafit.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12 //
13 function [p,err]=datafit(iprint,G,varargin)
14     //
15     //         [p,err]=datafit([iprint,] G [,DG],Z [,W],...)
16     //
17     //         Function used for fitting data to a model.
18     // For a given function G(p,z), this function finds the best vector
19     // of parameters p for approximating G(p,z_i)=0 for a set of measurement
20     // vectors z_i. Vector p is found by minimizing
21     //    G(p,z_1)'WG(p,z_1)+G(p,z_2)'WG(p,z_2)+...+G(p,z_n)'WG(p,z_n)
22     //
23     //      G: Scilab function (e=G(p,z), e: nex1, p: npx1, z: nzx1)
24     //     p0: initial guess (size npx1)
25     //      Z: matrix [z_1,z_2,...z_n] where z_i (nzx1) is the ith measurement
26     //      W: weighting matrix of size nexne (optional)
27     //     DG: partial of G wrt p (optional; S=DG(p,z), S: nexnp)
28     //
29     //                     Examples
30     //
31     //deff('y=FF(x)','y=a*(x-b)+c*x.*x')
32     //X=[];Y=[];
33     //a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
34     //Z=[Y;X];
35     //deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
36     //[p,err]=datafit(G,Z,[3;5;10])
37     //scf(0)
38     //clf();
39     //plot2d(X',Y',-1)
40     //plot2d(X',FF(X)',5,'002')
41     //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
42     //
43     //a=34;b=12;c=14;
44     //deff('s=DG(p,z)','y=z(1),x=z(2),s=-[x-p(2),-p(1),x*x]')
45     //[p,err]=datafit(G,DG,Z,[3;5;10])
46     //scf(1)
47     //clf();
48     //plot2d(X',Y',-1)
49     //plot2d(X',FF(X)',5,'002')
50     //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
51     //
52
53     [lhs,rhs]=argn(0)
54
55     if type(iprint)<>1 then
56         wflag=warning("query") // Disable warnings as the following lines may produce some
57         warning("off")
58         varargin(0)=G
59         G=iprint
60         iprint=0
61         warning(wflag)
62     end
63
64     if type(G)==15 then
65         Gparams=G;Gparams(1)=null();
66         G=G(1)
67     else
68         Gparams=list()
69     end
70
71
72     DG=varargin(1)
73     if type(DG)==10|type(DG)==13 then
74         GR=%t  //Jacobian provided
75         varargin(1)=null()
76     elseif type(DG)==15 then
77         error(msprintf(gettext("%s: Jacobian cannot be a list, parameters must be set in G."),"datafit"));
78     else
79         GR=%f
80     end
81
82     Z=varargin(1);
83     varargin(1)=null()
84     if type(Z)<>1 then
85         error(msprintf(gettext("%s: Wrong measurement matrix."),"datafit"));
86     end
87
88     nv=size(varargin)
89     if nv>=1 then
90         if size(varargin(1),2)==1 then // p0 ou 'b'
91             W=1
92         else
93             W=varargin(1);varargin(1)=null()
94             if size(W,1)~=size(W,2) then
95                 if size(W,1)==1 then
96                     error(msprintf(gettext("%s: Initial guess must be a column vector."),"datafit"));
97                 else
98                     error(msprintf(gettext("%s: Weighting matrix must be square."),"datafit"));
99                 end
100             end
101         end
102     end
103     if type(varargin(1))==1 then // p0
104         p0=varargin(1)
105     else
106         p0=varargin(4)
107     end
108
109     [mz,nz]=size(Z);np=size(p0,"*");
110
111     if type(G)==10 then  //form function to call hard coded external
112         if size(Gparams)==0 then
113             error(msprintf(gettext("%s: With hard coded function, user must give output size of G."),"datafit"));
114         end
115         m=Gparams(1);Gparams(1)=null()
116
117         // foo(m,np,p,mz,nz,Z,pars,f)
118         deff("f=G(p,Z)","f=call(''"+G+"'',"+..
119         "m,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'',"+..
120         "pars,7,''out'',["+string(m)+",1],8,''d'')")
121
122         pars=[];
123         for k=1:size(Gparams)
124             p=Gparams(k)
125             pars=[pars;p(:)]
126         end
127         Gparams=list()
128     end
129
130     if type(DG)==10 then //form function to call hard coded external
131         // dfoo(m,np,p,mz,nz,Z,pars,f)
132         deff("f=DG(p,Z)","f=call(''"+DG+"'',"+..
133         "m,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'',"+..
134         "pars,7,''out'',["+string(m)+","+string(np)+"],8,''d'')")
135     end
136
137
138     // form square cost gradient function DGG
139
140     if Gparams==list() then
141         GP   = "G(p,Z(:,i))"
142         GPPV = "G(p+v,Z(:,i))"
143         DGP  = "DG(p,Z(:,i))"
144     else
145         GP   = "G(p,Z(:,i),Gparams(:))"
146         GPPV = "G(p+v,Z(:,i),Gparams(:))"
147         DGP  = "DG(p,Z(:,i),Gparams(:))"
148     end
149
150     if ~GR then // finite difference
151         DGG=["g=0*p";
152         "pa=sqrt(%eps)*(1+1d-3*abs(p))"
153         "f=0;"
154         "for i=1:"+string(nz)
155         "  g1="+GP
156         "  f=f+g1''*W*g1"
157         "end"
158         "for j=1:"+string(np)
159         "  v=0*p;v(j)=pa(j),"
160         "  e=0;"
161         "  for i=1:"+string(nz)
162         "    g1="+GPPV
163         "    e=e+g1''*W*g1"
164         "  end"
165         "  g(j)=e-f;"
166         "end;"
167         "g=g./pa;"]
168     else // using Jacobian of G
169         DGG="g=0;for i=1:nz,g=g+2*"+DGP+"''*W*"+GP+",end"
170     end
171
172     // form cost function for optim
173     deff("[f,g,ind]=costf(p,ind)",[
174     "if ind==2|ind==4 then "
175     "  f=0;"
176     "   for i=1:"+string(nz)
177     "     g1="+GP
178     "     f=f+g1''*W*g1"
179     "   end"
180     "else "
181     "  f=0;"
182     "end";
183     "if ind==3|ind==4 then"
184     DGG
185     "else"
186     " g=0*p;"
187     "end"])
188
189     [err,p]=optim(costf,varargin(:),iprint=iprint)
190
191
192 endfunction