* Bugs 5824, 7732, 15344 fixed: datafit() upgraded + page overhauled
[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 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 // Copyright (C) 2018 - Samuel GOUGEON
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, dmin, status] = datafit(iprint, G, varargin)
14     //
15     //         [p,err]=datafit([iprint,] G [,DG], Data [,Wd] [,Wg],...)
16     //
17     // Data: matrix where Data(:,i) are the nzc coordinates of the ith measurement point
18     //       nz  = size(Data,"c") : number of measured available points.
19     //       nzc = size(Data,"r") : number of coordinates for each Data point.
20     //    G: Scilab function (e = G(p,Data) (e: ne x nz)
21     //   p0: initial guess of parameters values (vector or matrix)
22     //       np = length(p0) = length(p) : number of fitting parameters
23     //   Wd: Data weights (optionnal). Row of length nz.
24     //   Wg: weighting matrix of gaps definitions, of size ng x ng (optional)
25     //   DG: partial of G wrt p (optional; S=DG(p,Data), S: ne x np x nz)
26     //
27     //      p: vector of best parameters values (minimizing dmin)
28     // status: integer as returned by optim() when the default "qn" algo is used.
29     //   dmin: scalar = average distance of data points to model.
30
31     [lhs,rhs] = argn(0)
32
33     // Processing the 2 first argin imp and G:
34     if type(iprint)<>1 then      // 1st argin iprint is actually G (no iprint directive)
35         wflag = warning("query") // Disable warnings as the following lines may produce some
36         warning("off")
37         varargin(0) = G     // so the argin G is the first arg after the true G
38         G = iprint          // G gets its true value
39         iprint = 0          // iprint gets its default value
40         warning(wflag)
41     end
42
43     if type(G)==15 then     // G = list(G, params)
44         Gparams = G;
45         Gparams(1) = null();
46         G=G(1)
47     else
48         Gparams = list()    // no explicit parameter
49     end
50
51     // Next argins: optional dG (Jacobian), or mandatory points Data
52     DG = varargin(1)
53     if type(DG)==10 | type(DG)==13 then   // is a function
54         GR = %t      // Jacobian provided
55         varargin(1) = null()
56     elseif type(DG)==15 then
57         msg = gettext("%s: Jacobian cannot be a list, parameters must be set in G.")
58         error(msprintf(msg, "datafit"));
59     else
60         GR = %f  // Jacobian not provided
61     end
62
63     Data = varargin(1);          // Data points to fit
64     varargin(1) = null()
65     if type(Data) <>1 | ~isreal(Data,0) then
66         msg = gettext("%s: Data points must be real numbers.");
67         error(msprintf(msg, "datafit"));
68     end
69     [mz,nz] = size(Data);    // nz: Number of points
70
71     // Data weights: Must be a line vector of length nz (Wg can _neither_ be so long)
72     if size(varargin(1),2)==nz
73         Wd = varargin(1)
74         varargin(1) = null()
75     else
76         Wd = ones(1,nz);
77     end
78     // Wg assignment / initialization  (from [,Wg][,'b', pinf, psup], p0 [,algo])
79     nv = size(varargin);
80     if nv==0 then
81         msg =  gettext("%s: Initial guess p0 is missing.");
82         error(msprintf(msg, "datafit"));
83     elseif nv==1 then   // p0 expected
84         Wg = 1;
85     else
86         v1 = varargin(1);
87         v2 = varargin(2);
88         if v1=="b"| (nv>=5 & or(varargin(5)==["qn" "gc" "nd" "ar" "in"]))
89             Wg = 1
90         elseif (nv>1 & v2=="b") | ..
91                (nv>=6 & or(varargin(6)==["qn" "gc" "nd" "ar" "in"])) | ..
92                (type(v1)==1 & issquare(v1) & nv>1 & type(v2)==1 & (size(v2,"*")>1 | nv==2))
93             Wg = varargin(1)
94             if size(Wg,1) ~= size(Wg,2) then
95                 msg = gettext("%s: Weighting matrix must be square.");
96                 error(msprintf(msg, "datafit"));
97             end
98             varargin(1) = null()
99             nv = nv-1
100         else
101             Wg = 1;
102         end
103     end
104     // next: [,'b', pinf, psup], p0 [,algo])
105     if varargin(1)=="b" then
106         if nv<4
107             msg = gettext("%s: error with parameters bounds or initial guess\n");
108             error(msprintf(msg, "datafit"));
109         end
110         p0 = varargin(4);
111     else
112         p0 = varargin(1);
113     end
114     sp0 = size(p0);
115     np = size(p0,"*");    // Number of parameters
116     _format = format()
117     format("v",23)        // for string()
118
119     if type(G)==10 then   // form function to call hard coded external
120         if size(Gparams)==0 then
121             format(_format(2),_format(1))
122             msg = gettext("%s: With hard coded function, user must give output size of G.")
123             error(msprintf(msg, "datafit"));
124         end
125         m = Gparams(1);
126         Gparams(1) = null();
127
128         // foo(m,np,p,mz,nz,Data[,Wd],pars,f)
129         tmp = "f=call(''" + G + "''," + ..
130               "m,1,''i'', np,2,''i'', p,3,''d'', mz,4,''i'', nz,5,''i''," + ..
131               "Data,6,''d'', pars,7,''out'',[" + string(m) + ",1],8,''d'')"
132         deff("f = G(p,Data)", tmp)
133
134         pars = [];
135         for k = 1:size(Gparams)
136             p = Gparams(k)
137             pars = [pars;p(:)]
138         end
139         Gparams = list()
140     end
141
142     if type(DG)==10 then // form function to call hard coded external
143         // dfoo(m,np,p,mz,nz,Data[,Wd],pars,f)
144         tmp = "f = call(''" + DG + "''," + ..
145           "m,1,''i'', np,2,''i'', p,3,''d'', mz,4,''i'', nz,5,''i'',Data,6,''d''," + ..
146           "pars,7,''out'',[" + string(m) + "," + string(np) + "],8,''d'')"
147         deff("f = DG(p,Data)", tmp)
148     end
149
150     // Is the G function vectorized?
151     try
152         if Gparams==list() then
153             e = G(p0,[Data(:,1) Data(:,1)]);
154         else
155             e = G(p0,[Data(:,1) Data(:,1)], Gparams(:));
156         end
157         isGvec = size(e,2)==2
158     catch
159         isGvec = %f;
160     end
161
162     // form square cost gradient function DGG
163     if Gparams==list() then
164         GP   = "G(p,  Data(:,i))"
165         GPPV = "G(p+v,Data(:,i))"
166         DGP  = "DG(p, Data(:,i))"
167     else
168         GP   = "G(p,  Data(:,i), Gparams(:))"
169         GPPV = "G(p+v,Data(:,i), Gparams(:))"
170         DGP  = "DG(p, Data(:,i), Gparams(:))"
171     end
172
173     if ~GR then // finite difference
174         // Version vectorized over Data:
175         if isGvec then
176             DGG = [
177             "g  = 0*p"
178             "pa = sqrt(%eps)*(1+1d-3*abs(p))"
179             "i = : "
180             "g1 = " + GP
181             "f = sum(sum((g1''*Wg) .* g1'',""c"").*Wd'');"
182             "for j = 1:" + msprintf("%d\n",np)
183             "    v = 0*p;"
184             "    v(j) = pa(j);"
185             "    g1 = " + GPPV
186             "    e = sum(sum((g1''*Wg) .* g1'',""c"").*Wd'');"
187             "    g(j) = e - f;"
188             "end"
189             "g = g ./ pa;"
190             ];
191         else
192             DGG = [
193             "g  = 0*p"
194             "pa = sqrt(%eps)*(1+1d-3*abs(p))"
195             "f = 0;"
196             "for i = 1:" + msprintf("%d\n", nz)
197             "    g1 = " + GP
198             "    f = f + g1''*Wg*g1 * Wd(i)"
199             "end"
200             "for j = 1:" + msprintf("%d\n", np)
201             "    v = 0*p;"
202             "    v(j) = pa(j),"
203             "    e = 0;"
204             "    for i = 1:" + msprintf("%d\n", nz)
205             "        g1 = " + GPPV
206             "        e = e + g1''*Wg*g1 * Wd(i)"
207             "    end"
208             "    g(j) = e - f;"
209             "end"
210             "g = g ./ pa;"
211             ];
212         end
213     else // using Jacobian of G  (UNVECTORIZED)
214         DGG = "g=0; for i=1:nz, g = g + 2*"+DGP+"''*Wg*"+GP+", end"
215     end
216
217     // Defining the costf() function for optim:
218     if isGvec then
219         tmp = [
220             "    i = : "
221             "    g1 = " + GP
222             "    f = sum(sum((g1''*Wg) .* g1'',""c"") .* Wd'');"
223             ]
224     else
225         tmp = [
226             "    f = 0;"
227             "    for i = 1:" + msprintf("%d\n", nz)
228             "        g1 = " + GP
229             "        f = f + g1''*Wg*g1 * Wd(i)"
230             "    end"
231             ]
232     end
233     deff("[f,g,ind] = costf(p,ind)",[
234         "if ind==2 | ind==4 then "
235         tmp
236         "else"
237         "    f = 0;"
238         "end";
239         "if ind==3 | ind==4 then"
240         DGG
241         "else"
242         "    g = 0*p;"
243         "end"
244         ])
245
246     format(_format(2), _format(1))
247     qn = %t;
248     for v = varargin
249         if or(v==['gc' 'nd'])
250             qn = %f;
251             break
252         end
253     end
254     if qn then
255         [fmin, p, t,t,t,t, status] = optim(costf, varargin(:), iprint=iprint)
256     else
257         [fmin, p] = optim(costf, varargin(:), iprint=iprint)
258         status = %nan;
259     end
260
261     p = matrix(p, sp0);
262
263     // fmin => dmin: average Model-to-data distance
264     ng = size(G(p, Data(:,1)),1);
265     dmin = sqrt(fmin/ng/sum(Wd));
266 endfunction