Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
M
miscmaths
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Container Registry
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
FSL
miscmaths
Commits
78c36d65
Commit
78c36d65
authored
17 years ago
by
Jesper Andersson
Browse files
Options
Downloads
Patches
Plain Diff
Qualified abs, max, min and pow with std:: to avoid problems with overloaded functions in miacmaths
parent
38118c40
No related branches found
Branches containing commit
No related tags found
Tags containing commit
No related merge requests found
Changes
1
Hide whitespace changes
Inline
Side-by-side
Showing
1 changed file
nonlin.cpp
+26
-26
26 additions, 26 deletions
nonlin.cpp
with
26 additions
and
26 deletions
nonlin.cpp
+
26
−
26
View file @
78c36d65
...
@@ -111,7 +111,7 @@ ReturnMatrix NonlinCF::grad(const ColumnVector& p) const
...
@@ -111,7 +111,7 @@ ReturnMatrix NonlinCF::grad(const ColumnVector& p) const
double
tiny
=
1e-8
;
double
tiny
=
1e-8
;
double
cf0
=
cf
(
tmpp
);
double
cf0
=
cf
(
tmpp
);
for
(
int
i
=
0
;
i
<
p
.
Nrows
();
i
++
)
{
for
(
int
i
=
0
;
i
<
p
.
Nrows
();
i
++
)
{
double
step
=
tiny
*
max
(
tmpp
.
element
(
i
),
1.0
);
double
step
=
tiny
*
std
::
max
(
tmpp
.
element
(
i
),
1.0
);
tmpp
.
element
(
i
)
+=
step
;
tmpp
.
element
(
i
)
+=
step
;
gradv
.
element
(
i
)
=
(
cf
(
tmpp
)
-
cf0
)
/
step
;
gradv
.
element
(
i
)
=
(
cf
(
tmpp
)
-
cf0
)
/
step
;
tmpp
.
element
(
i
)
-=
step
;
tmpp
.
element
(
i
)
-=
step
;
...
@@ -146,7 +146,7 @@ boost::shared_ptr<BFMatrix> NonlinCF::hess(const ColumnVector& p,
...
@@ -146,7 +146,7 @@ boost::shared_ptr<BFMatrix> NonlinCF::hess(const ColumnVector& p,
// First calculate all f(x+dx_i) values
// First calculate all f(x+dx_i) values
for
(
int
i
=
0
;
i
<
p
.
Nrows
();
i
++
)
{
for
(
int
i
=
0
;
i
<
p
.
Nrows
();
i
++
)
{
step
.
element
(
i
)
=
tiny
*
max
(
tmpp
.
element
(
i
),
1.0
);
step
.
element
(
i
)
=
tiny
*
std
::
max
(
tmpp
.
element
(
i
),
1.0
);
tmpp
.
element
(
i
)
+=
step
.
element
(
i
);
tmpp
.
element
(
i
)
+=
step
.
element
(
i
);
fdx
.
element
(
i
)
=
cf
(
tmpp
);
fdx
.
element
(
i
)
=
cf
(
tmpp
);
tmpp
.
element
(
i
)
-=
step
.
element
(
i
);
tmpp
.
element
(
i
)
-=
step
.
element
(
i
);
...
@@ -647,7 +647,7 @@ LinOut linsrch(// Input
...
@@ -647,7 +647,7 @@ LinOut linsrch(// Input
double
almin
=
0.0
;
double
almin
=
0.0
;
for
(
int
i
=
0
;
i
<
p0
.
Nrows
();
i
++
)
{
for
(
int
i
=
0
;
i
<
p0
.
Nrows
();
i
++
)
{
almin
=
max
(
almin
,
abs
(
pdir
.
element
(
i
))
/
max
(
abs
(
p0
.
element
(
i
)),
1.0
));
almin
=
std
::
max
(
almin
,
std
::
abs
(
pdir
.
element
(
i
))
/
std
::
max
(
std
::
abs
(
p0
.
element
(
i
)),
1.0
));
}
}
almin
=
ptol
/
almin
;
almin
=
ptol
/
almin
;
...
@@ -665,8 +665,8 @@ LinOut linsrch(// Input
...
@@ -665,8 +665,8 @@ LinOut linsrch(// Input
*
lambda
=
-
fp0
/
(
2.0
*
(
f2
-
f0
-
fp0
));
// Minumum of f(lambda)
*
lambda
=
-
fp0
/
(
2.0
*
(
f2
-
f0
-
fp0
));
// Minumum of f(lambda)
// Make sure new lambda is 0.1*old_l < lambda < 0.5*old_l
// Make sure new lambda is 0.1*old_l < lambda < 0.5*old_l
*
lambda
=
max
(
lmin
,
*
lambda
);
*
lambda
=
std
::
max
(
lmin
,
*
lambda
);
*
lambda
=
min
(
lmax
,
*
lambda
);
*
lambda
=
std
::
min
(
lmax
,
*
lambda
);
(
*
np
)
=
p0
+
(
*
lambda
)
*
pdir
;
// Second set of new parameters to try
(
*
np
)
=
p0
+
(
*
lambda
)
*
pdir
;
// Second set of new parameters to try
double
f1
=
sf
*
cfo
.
cf
(
*
np
);
// Cost-function value for par
double
f1
=
sf
*
cfo
.
cf
(
*
np
);
// Cost-function value for par
...
@@ -686,14 +686,14 @@ LinOut linsrch(// Input
...
@@ -686,14 +686,14 @@ LinOut linsrch(// Input
// See if present value is acceptable
// See if present value is acceptable
if
(
f1
<
f0
+
alpha
*
(
*
lambda
)
*
DotProduct
(
grad
,(
*
np
)
-
p0
))
{
*
of
=
f1
;
return
(
LM_CONV
);}
if
(
f1
<
f0
+
alpha
*
(
*
lambda
)
*
DotProduct
(
grad
,(
*
np
)
-
p0
))
{
*
of
=
f1
;
return
(
LM_CONV
);}
// Find parameter values for cubic and square on lambda
// Find parameter values for cubic and square on lambda
X
<<
pow
(
l1
,
3
)
<<
pow
(
l1
,
2
)
<<
pow
(
l2
,
3
)
<<
pow
(
l2
,
2
);
X
<<
std
::
pow
(
l1
,
3
.0
)
<<
std
::
pow
(
l1
,
2
.0
)
<<
std
::
pow
(
l2
,
3
.0
)
<<
std
::
pow
(
l2
,
2
.0
);
y
<<
f1
-
fp0
*
l1
-
f0
<<
f2
-
fp0
*
l2
-
f0
;
y
<<
f1
-
fp0
*
l1
-
f0
<<
f2
-
fp0
*
l2
-
f0
;
ColumnVector
b
=
X
.
i
()
*
y
;
ColumnVector
b
=
X
.
i
()
*
y
;
// Find value for lambda that yield minimum of cubic
// Find value for lambda that yield minimum of cubic
*
lambda
=
(
-
b
.
element
(
1
)
+
sqrt
(
pow
(
b
.
element
(
1
),
2
)
-
3.0
*
b
.
element
(
0
)
*
fp0
))
/
(
3.0
*
b
.
element
(
0
));
*
lambda
=
(
-
b
.
element
(
1
)
+
sqrt
(
std
::
pow
(
b
.
element
(
1
),
2
.0
)
-
3.0
*
b
.
element
(
0
)
*
fp0
))
/
(
3.0
*
b
.
element
(
0
));
// Make sure new lambda is 0.1*old_l < lambda < 0.5*old_l
// Make sure new lambda is 0.1*old_l < lambda < 0.5*old_l
*
lambda
=
max
(
lmin
*
l1
,
*
lambda
);
*
lambda
=
std
::
max
(
lmin
*
l1
,
*
lambda
);
*
lambda
=
min
(
lmax
*
l1
,
*
lambda
);
*
lambda
=
std
::
min
(
lmax
*
l1
,
*
lambda
);
// Get new function value and update parameters
// Get new function value and update parameters
f2
=
f1
;
f2
=
f1
;
(
*
np
)
=
p0
+
(
*
lambda
)
*
pdir
;
(
*
np
)
=
p0
+
(
*
lambda
)
*
pdir
;
...
@@ -781,27 +781,27 @@ LinOut linmin(// Input
...
@@ -781,27 +781,27 @@ LinOut linmin(// Input
for
(
int
i
=
0
;
i
<
maxiter
;
i
++
)
{
for
(
int
i
=
0
;
i
<
maxiter
;
i
++
)
{
double
midp
=
(
rp
.
first
+
lp
.
first
)
/
2.0
;
// Midpoint of bracketing points
double
midp
=
(
rp
.
first
+
lp
.
first
)
/
2.0
;
// Midpoint of bracketing points
double
tol
=
2.0
*
ftol
*
abs
(
x
->
first
)
+
MISCMATHS
::
EPS
;
// Absolute tolerance
double
tol
=
2.0
*
ftol
*
std
::
abs
(
x
->
first
)
+
MISCMATHS
::
EPS
;
//
Std::
Absolute tolerance
if
(
abs
(
x
->
first
-
midp
)
<=
(
tol
-
0.5
*
(
rp
.
first
-
lp
.
first
)))
{
// Convergence check
if
(
std
::
abs
(
x
->
first
-
midp
)
<=
(
tol
-
0.5
*
(
rp
.
first
-
lp
.
first
)))
{
// Convergence check
return
(
LM_CONV
);
return
(
LM_CONV
);
}
}
// Try parabolic fit, but not before third iteration
// Try parabolic fit, but not before third iteration
double
tmp
=
10.0
*
sqrt
(
MISCMATHS
::
EPS
);
double
tmp
=
10.0
*
sqrt
(
MISCMATHS
::
EPS
);
if
(
abs
(
ostep
)
>
tol
/
2.0
&&
// If second to last step big enough
if
(
std
::
abs
(
ostep
)
>
tol
/
2.0
&&
// If second to last step big enough
abs
(
x
->
first
-
w
.
first
)
>
tmp
&&
std
::
abs
(
x
->
first
-
w
.
first
)
>
tmp
&&
abs
(
x
->
first
-
v
.
first
)
>
tmp
&&
std
::
abs
(
x
->
first
-
v
.
first
)
>
tmp
&&
abs
(
w
.
first
-
v
.
first
)
>
tmp
)
{
// And points not degenerate
std
::
abs
(
w
.
first
-
v
.
first
)
>
tmp
)
{
// And points not degenerate
step
=
ostep
;
step
=
ostep
;
ostep
=
d
;
ostep
=
d
;
y
<<
x
->
second
<<
w
.
second
<<
v
.
second
;
y
<<
x
->
second
<<
w
.
second
<<
v
.
second
;
X
<<
pow
(
x
->
first
,
2.0
)
<<
x
->
first
<<
1.0
<<
X
<<
std
::
pow
(
x
->
first
,
2.0
)
<<
x
->
first
<<
1.0
<<
pow
(
w
.
first
,
2.0
)
<<
w
.
first
<<
1.0
<<
std
::
pow
(
w
.
first
,
2.0
)
<<
w
.
first
<<
1.0
<<
pow
(
v
.
first
,
2.0
)
<<
v
.
first
<<
1.0
;
std
::
pow
(
v
.
first
,
2.0
)
<<
v
.
first
<<
1.0
;
ColumnVector
b
=
X
.
i
()
*
y
;
ColumnVector
b
=
X
.
i
()
*
y
;
if
(
b
.
element
(
0
)
<
4
*
MISCMATHS
::
EPS
||
// If on line or going for maximum
if
(
b
.
element
(
0
)
<
4
*
MISCMATHS
::
EPS
||
// If on line or going for maximum
(
test
.
first
=
-
b
.
element
(
1
)
/
(
2.0
*
b
.
element
(
0
)))
<=
lp
.
first
(
test
.
first
=
-
b
.
element
(
1
)
/
(
2.0
*
b
.
element
(
0
)))
<=
lp
.
first
||
test
.
first
>=
rp
.
first
||
// If outside bracketed interval
||
test
.
first
>=
rp
.
first
||
// If outside bracketed interval
abs
(
test
.
first
-
x
->
first
)
>
0.5
*
step
)
{
// Or if step too big (indicates oscillation)
std
::
abs
(
test
.
first
-
x
->
first
)
>
0.5
*
step
)
{
// Or if step too big (indicates oscillation)
// Take golden step into larger interval
// Take golden step into larger interval
if
(
rp
.
first
-
x
->
first
>
x
->
first
-
lp
.
first
)
{
// If right interval larger
if
(
rp
.
first
-
x
->
first
>
x
->
first
-
lp
.
first
)
{
// If right interval larger
test
.
first
=
x
->
first
+
gold
*
(
rp
.
first
-
x
->
first
);
test
.
first
=
x
->
first
+
gold
*
(
rp
.
first
-
x
->
first
);
...
@@ -869,7 +869,7 @@ pair<double,double> bracket(// Input
...
@@ -869,7 +869,7 @@ pair<double,double> bracket(// Input
// Find maximum relative component of search direction
// Find maximum relative component of search direction
double
test
=
0.0
;
double
test
=
0.0
;
for
(
int
i
=
0
;
i
<
pdir
.
Nrows
();
i
++
)
{
test
=
max
(
test
,
abs
(
pdir
.
element
(
i
))
/
max
(
p
.
element
(
i
),
1.0
));}
for
(
int
i
=
0
;
i
<
pdir
.
Nrows
();
i
++
)
{
test
=
std
::
max
(
test
,
std
::
abs
(
pdir
.
element
(
i
))
/
std
::
max
(
p
.
element
(
i
),
1.0
));}
// Do a crude initial search for order of magnitude
// Do a crude initial search for order of magnitude
...
@@ -903,8 +903,8 @@ pair<double,double> bracket(// Input
...
@@ -903,8 +903,8 @@ pair<double,double> bracket(// Input
return
(
p_l
);
return
(
p_l
);
}
}
// Let's see if a parabolic might help us
// Let's see if a parabolic might help us
if
(
abs
(
l2
-
l1
)
>
10.0
*
sqrt
(
MISCMATHS
::
EPS
))
{
if
(
std
::
abs
(
l2
-
l1
)
>
10.0
*
sqrt
(
MISCMATHS
::
EPS
))
{
X
<<
pow
(
l1
,
2.0
)
<<
l1
<<
pow
(
l2
,
2.0
)
<<
l2
;
X
<<
std
::
pow
(
l1
,
2.0
)
<<
l1
<<
std
::
pow
(
l2
,
2.0
)
<<
l2
;
y
<<
cf1
<<
cf2
;
y
<<
cf1
<<
cf2
;
ColumnVector
b
=
X
.
i
()
*
y
;
ColumnVector
b
=
X
.
i
()
*
y
;
if
(
b
.
element
(
0
)
>
4.0
*
MISCMATHS
::
EPS
)
{
// Check they are not on a line and not for maximum
if
(
b
.
element
(
0
)
>
4.0
*
MISCMATHS
::
EPS
)
{
// Check they are not on a line and not for maximum
...
@@ -960,9 +960,9 @@ bool zero_grad_conv(const ColumnVector& par,
...
@@ -960,9 +960,9 @@ bool zero_grad_conv(const ColumnVector& par,
{
{
double
test
=
0.0
;
// test will be largest relative component of gradient
double
test
=
0.0
;
// test will be largest relative component of gradient
for
(
int
i
=
0
;
i
<
par
.
Nrows
();
i
++
)
{
for
(
int
i
=
0
;
i
<
par
.
Nrows
();
i
++
)
{
test
=
max
(
test
,
abs
(
grad
.
element
(
i
))
*
max
(
abs
(
par
.
element
(
i
)),
1.0
));
test
=
std
::
max
(
test
,
std
::
abs
(
grad
.
element
(
i
))
*
std
::
max
(
std
::
abs
(
par
.
element
(
i
)),
1.0
));
}
}
test
/=
max
(
cf
,
1.0
);
// Protect against near-zero values for cost-function
test
/=
std
::
max
(
cf
,
1.0
);
// Protect against near-zero values for cost-function
return
(
test
<
gtol
);
return
(
test
<
gtol
);
}
}
...
@@ -973,7 +973,7 @@ bool zero_cf_diff_conv(double cfo,
...
@@ -973,7 +973,7 @@ bool zero_cf_diff_conv(double cfo,
double
cfn
,
double
cfn
,
double
cftol
)
double
cftol
)
{
{
return
(
2.0
*
abs
(
cfo
-
cfn
)
<=
cftol
*
(
abs
(
cfo
)
+
abs
(
cfn
)
+
MISCMATHS
::
EPS
));
return
(
2.0
*
std
::
abs
(
cfo
-
cfn
)
<=
cftol
*
(
std
::
abs
(
cfo
)
+
std
::
abs
(
cfn
)
+
MISCMATHS
::
EPS
));
}
}
// Based on zero (neglible) step in parameter space
// Based on zero (neglible) step in parameter space
...
@@ -984,7 +984,7 @@ bool zero_par_step_conv(const ColumnVector& par,
...
@@ -984,7 +984,7 @@ bool zero_par_step_conv(const ColumnVector& par,
{
{
double
test
=
0.0
;
double
test
=
0.0
;
for
(
int
i
=
0
;
i
<
par
.
Nrows
();
i
++
)
{
for
(
int
i
=
0
;
i
<
par
.
Nrows
();
i
++
)
{
test
=
max
(
test
,
abs
(
step
.
element
(
i
))
/
max
(
abs
(
par
.
element
(
i
)),
1.0
));
test
=
std
::
max
(
test
,
std
::
abs
(
step
.
element
(
i
))
/
std
::
max
(
std
::
abs
(
par
.
element
(
i
)),
1.0
));
}
}
return
(
test
<
ptol
);
return
(
test
<
ptol
);
}
}
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment