1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function [f,x,g]=leastsq(imp,fun,varargin)
// n p
// min sum (fun(x).^2) where fun is a function from R to R
// x
//
// [f]=fun(x) computes the value f of the function at the point x
// and the gradient g of f at x g(i,j)=Dfi/dxj
if type(imp)<>1 then
varargin(0)=fun
fun=imp
imp=0
end
if type(fun)==15 then
fn=fun(1);params=fun;params(1)=null()
else
fn=fun;params=list()
end
Dfun=varargin(1),kr=1;
if type(Dfun)==10 then //the 'b' keyword or the jacobian entry point name
if Dfun=="b" & size(varargin) >= 3 then
if type(varargin(2))==1 & type(varargin(3))==1 then
J=%f //bounds specification
else
J=%t //jacobian
end
else
J=%t //jacobian
end
elseif type(Dfun)==11|type(Dfun)==13 then
J=%t //Jacobian provided
elseif type(Dfun)==15 then
error(msprintf(gettext("%s: Jacobian cannot be a list, parameters must be set in fun."),"leastsq"));
else
J=%f;
end
if J then, varargin(1)=null(), end // to correct bug 1219 (bruno, 22 feb 2005)
kr=1
if varargin(kr)=="b" then kr=kr+3,end
x0=varargin(kr)
if type(fn)==10 then //hard coded function given by its name
if size(params)==0 then
error(msprintf(gettext("%s: With hard coded function, user must give output size of fun."),"leastsq"));
end
m=params(1);params(1)=null()
n=size(x0,"*")
// foo(m,nx,x,params,f)
deff("f=fn(x)","f=call(''"+fn+"'',"+..
"m,1,''i'',n,2,''i'',x,3,''d'',"+..
"pars,4,''d'',''out'',["+string(m)+",1],5,''d'')")
pars=[];
for k=1:size(params)
p=params(k)
pars=[pars;p(:)]
end
params=list()
end
if J then //jacobian given
if type(Dfun)==10 then //form function to call hard coded external
// dfoo(m,nx,x,params,g)
deff("g=Dfun(x)","g=call(''"+Dfun+"'',"+..
"m,1,''i'',n,2,''i'',x,3,''d'',"+..
"pars,4,''d'',''out'',["+string(m)+","+string(n)+"],5,''d'')")
end
else
if params==list() then
deff("g=Dfun(x)","g=numderivative(fn,x)")
else
deff("g=Dfun(x,varargin)","g=numderivative(list(fn,varargin(:)),x)")
end
end
if params==list() then
deff("[f,g,ind]=%opt(x,ind)",[
"ff=fn(x);gf=Dfun(x)"
"f=sum(ff.^2)"
"g=2*(gf''*ff(:))"])
else
deff("[f,g,ind]=%opt(x,ind)",[
"ff=fn(x,params(:));gf=Dfun(x,params(:))"
"f=sum(ff.^2)"
"g=2*(gf''*ff(:))"])
end
[f,x,g]=optim(%opt,varargin(:),imp=imp)
endfunction
|
