Added CeCILL headers
[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 // 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-en.txt
9 //
10 function [p,err]=datafit(imp,G,varargin)
11 //
12 //         [p,err]=datafit([imp,] G [,DG],Z [,W],...)
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)
24 //     DG: partial of G wrt p (optional; S=DG(p,z), S: nexnp)
25 //
26 //                     Examples
27 //
28 //deff('y=FF(x)','y=a*(x-b)+c*x.*x')
29 //X=[];Y=[];
30 //a=34;b=12;c=14;for x=0:.1:3, Y=[Y,FF(x)+100*(rand()-.5)];X=[X,x];end
31 //Z=[Y;X];
32 //deff('e=G(p,z)','a=p(1),b=p(2),c=p(3),y=z(1),x=z(2),e=y-FF(x)')
33 //[p,err]=datafit(G,Z,[3;5;10])
34 //xset('window',0)
35 //xbasc();
36 //plot2d(X',Y',-1) 
37 //plot2d(X',FF(X)',5,'002')
38 //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
39 //
40 //a=34;b=12;c=14;
41 //deff('s=DG(p,z)','y=z(1),x=z(2),s=-[x-p(2),-p(1),x*x]')
42 //[p,err]=datafit(G,DG,Z,[3;5;10])
43 //xset('window',1)
44 //xbasc(); 
45 //plot2d(X',Y',-1) 
46 //plot2d(X',FF(X)',5,'002')
47 //a=p(1),b=p(2),c=p(3);plot2d(X',FF(X)',12,'002')
48 //
49
50 [lhs,rhs]=argn(0)
51
52 if type(imp)<>1 then
53   varargin(0)=G
54   G=imp
55   imp=0
56 end
57
58 if type(G)==15 then
59   Gparams=G;Gparams(1)=null();
60   G=G(1)
61 else
62   Gparams=list()
63 end
64
65   
66 DG=varargin(1)
67 if type(DG)==10|type(DG)==11|type(DG)==13 then 
68   GR=%t  //Jacobian provided
69   varargin(1)=null()
70 elseif type(DG)==15 then 
71   error('Jacobian cannot be a list, parameters must be set in G')
72 else
73   GR=%f
74 end
75
76 Z=varargin(1);
77 varargin(1)=null()
78 if type(Z)<>1 then 
79   error('datafit : wrong measurement matrix')
80 end
81
82 nv=size(varargin)
83 if nv>=1 then
84   if size(varargin(1),2)==1 then // p0 ou 'b'
85     W=1
86   else
87     W=varargin(1);varargin(1)=null()
88     if size(W,1)~=size(W,2) then 
89       error('Weighting matrix must be square');
90     end
91   end
92 end  
93 if type(varargin(1))==1 then // p0
94   p0=varargin(1)
95 else
96   p0=varargin(4)
97 end
98
99 [mz,nz]=size(Z);np=size(p0,'*');
100
101 if type(G)==10 then  //form function to call hard coded external
102   if size(Gparams)==0 then 
103     error('With hard coded function, user must give output size of G'),
104   end
105   m=Gparams(1);Gparams(1)=null()
106   
107   // foo(m,np,p,mz,nz,Z,pars,f)
108   deff('f=G(p,Z)','f=call('''+G+''','+..
109       'm,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'','+..
110       'pars,7,''out'',['+string(m)+',1],8,''d'')')
111   
112   pars=[];
113   for k=1:size(Gparams)
114     p=Gparams(k)
115     pars=[pars;p(:)]
116   end
117   Gparams=list()
118 end
119
120 if type(DG)==10 then //form function to call hard coded external
121   // dfoo(m,np,p,mz,nz,Z,pars,f)
122   deff('f=DG(p,Z)','f=call('''+DG+''','+..
123       'm,1,''i'',np,2,''i'',p,3,''d'',mz,4,''i'',nz,5,''i'',Z,6,''d'','+..
124       'pars,7,''out'',['+string(m)+','+string(np)+'],8,''d'')')
125 end
126
127
128 // form square cost gradient function DGG
129
130 if Gparams==list() then
131   GP   = 'G(p,Z(:,i))'
132   GPPV = 'G(p+v,Z(:,i))'
133   DGP  = 'DG(p,Z(:,i))'
134 else
135   GP   = 'G(p,Z(:,i),Gparams(:))'
136   GPPV = 'G(p+v,Z(:,i),Gparams(:))'
137   DGP  = 'DG(p,Z(:,i),Gparams(:))'  
138 end
139   
140 if ~GR then // finite difference
141   DGG=['g=0*p';
142        'pa=sqrt(%eps)*(1+1d-3*abs(p))'
143        'f=0;'
144        'for i=1:'+string(nz)
145        '  g1='+GP
146        '  f=f+g1''*W*g1'
147        'end'
148        'for j=1:'+string(np)
149        '  v=0*p;v(j)=pa(j),'
150        '  e=0;'
151        '  for i=1:'+string(nz)
152        '    g1='+GPPV
153        '    e=e+g1''*W*g1'
154        '  end'
155        '  g(j)=e-f;'
156        'end;'
157        'g=g./pa;']
158 else // using Jacobian of G
159   DGG='g=0;for i=1:nz,g=g+2*'+DGP+'''*W*'+GP+',end'
160 end
161
162 // form cost function for optim
163 deff('[f,g,ind]=costf(p,ind)',[
164     'if ind==2|ind==4 then '
165     '  f=0;'
166     '   for i=1:'+string(nz)
167     '     g1='+GP
168     '     f=f+g1''*W*g1'
169     '   end'
170     'else '
171     '  f=0;'
172     'end';
173     'if ind==3|ind==4 then' 
174       DGG
175     'else' 
176     ' g=0*p;'
177     'end'])
178
179 [err,p]=optim(costf,varargin(:),imp=imp)
180
181  
182 endfunction