f65bb587457a8846ae1aab1984731265770384b0
[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(imp,G,varargin)
14     //
15     //         [p,err]=datafit([imp,] 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     //xset('window',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     //xset('window',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(imp)<>1 then
56         varargin(0)=G
57         G=imp
58         imp=0
59     end
60
61     if type(G)==15 then
62         Gparams=G;Gparams(1)=null();
63         G=G(1)
64     else
65         Gparams=list()
66     end
67
68
69     DG=varargin(1)
70     if type(DG)==10|type(DG)==11|type(DG)==13 then
71         GR=%t  //Jacobian provided
72         varargin(1)=null()
73     elseif type(DG)==15 then
74         error(msprintf(gettext("%s: Jacobian cannot be a list, parameters must be set in G."),"datafit"));
75     else
76         GR=%f
77     end
78
79     Z=varargin(1);
80     varargin(1)=null()
81     if type(Z)<>1 then
82         error(msprintf(gettext("%s: Wrong measurement matrix."),"datafit"));
83     end
84
85     nv=size(varargin)
86     if nv>=1 then
87         if size(varargin(1),2)==1 then // p0 ou 'b'
88             W=1
89         else
90             W=varargin(1);varargin(1)=null()
91             if size(W,1)~=size(W,2) then
92                 if size(W,1)==1 then
93                     error(msprintf(gettext("%s: Initial guess must be a column vector."),"datafit"));
94                 else
95                     error(msprintf(gettext("%s: Weighting matrix must be square."),"datafit"));
96                 end
97             end
98         end
99     end
100     if type(varargin(1))==1 then // p0
101         p0=varargin(1)
102     else
103         p0=varargin(4)
104     end
105
106     [mz,nz]=size(Z);np=size(p0,"*");
107
108     if type(G)==10 then  //form function to call hard coded external
109         if size(Gparams)==0 then
110             error(msprintf(gettext("%s: With hard coded function, user must give output size of G."),"datafit"));
111         end
112         m=Gparams(1);Gparams(1)=null()
113
114         // foo(m,np,p,mz,nz,Z,pars,f)
115         deff("f=G(p,Z)","f=call(''"+G+"'',"+..
116         "m,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'',"+..
117         "pars,7,''out'',["+string(m)+",1],8,''d'')")
118
119         pars=[];
120         for k=1:size(Gparams)
121             p=Gparams(k)
122             pars=[pars;p(:)]
123         end
124         Gparams=list()
125     end
126
127     if type(DG)==10 then //form function to call hard coded external
128         // dfoo(m,np,p,mz,nz,Z,pars,f)
129         deff("f=DG(p,Z)","f=call(''"+DG+"'',"+..
130         "m,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'',"+..
131         "pars,7,''out'',["+string(m)+","+string(np)+"],8,''d'')")
132     end
133
134
135     // form square cost gradient function DGG
136
137     if Gparams==list() then
138         GP   = "G(p,Z(:,i))"
139         GPPV = "G(p+v,Z(:,i))"
140         DGP  = "DG(p,Z(:,i))"
141     else
142         GP   = "G(p,Z(:,i),Gparams(:))"
143         GPPV = "G(p+v,Z(:,i),Gparams(:))"
144         DGP  = "DG(p,Z(:,i),Gparams(:))"
145     end
146
147     if ~GR then // finite difference
148         DGG=["g=0*p";
149         "pa=sqrt(%eps)*(1+1d-3*abs(p))"
150         "f=0;"
151         "for i=1:"+string(nz)
152         "  g1="+GP
153         "  f=f+g1''*W*g1"
154         "end"
155         "for j=1:"+string(np)
156         "  v=0*p;v(j)=pa(j),"
157         "  e=0;"
158         "  for i=1:"+string(nz)
159         "    g1="+GPPV
160         "    e=e+g1''*W*g1"
161         "  end"
162         "  g(j)=e-f;"
163         "end;"
164         "g=g./pa;"]
165     else // using Jacobian of G
166         DGG="g=0;for i=1:nz,g=g+2*"+DGP+"''*W*"+GP+",end"
167     end
168
169     // form cost function for optim
170     deff("[f,g,ind]=costf(p,ind)",[
171     "if ind==2|ind==4 then "
172     "  f=0;"
173     "   for i=1:"+string(nz)
174     "     g1="+GP
175     "     f=f+g1''*W*g1"
176     "   end"
177     "else "
178     "  f=0;"
179     "end";
180     "if ind==3|ind==4 then"
181     DGG
182     "else"
183     " g=0*p;"
184     "end"])
185
186     [err,p]=optim(costf,varargin(:),imp=imp)
187
188
189 endfunction