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