Skip to content
Snippets Groups Projects
Commit d4686fb3 authored by Jesper Andersson's avatar Jesper Andersson
Browse files

Explicitly prefixed all sqrt with std::

parent 7f320a67
No related branches found
No related tags found
No related merge requests found
// Definitions for module nonlin
#include <ctime>
#include <iostream>
#include <fstream>
#include <iomanip>
......@@ -263,6 +264,7 @@ ColumnVector operator*(const VarmetMatrix& m, const ColumnVector& v)
}
}
// Gateway function to routines for non-linear optimisation
NonlinOut nonlin(const NonlinParam& p, const NonlinCF& cfo)
......@@ -480,8 +482,9 @@ NonlinOut sccngr(const NonlinParam& np, const NonlinCF& cfo)
while (np.NextIter()) {
double p2 = DotProduct(p,p); // p'*p, Temporary variable to save some time
if (success == true) { // If last step led to reduction of cos-function
double sigma_k = sigma/sqrt(p2); // Normalised step-length when estimating H*p
if (success == true) { // If last step led to reduction of cost-function
double sigma_k = sigma/std::sqrt(p2); // Normalised step-length when estimating H*p
// cout << "np.NIter() = " << np.NIter() << ", p2 = " << p2 << ", sigma_k = " << sigma_k << endl;
s = (cfo.grad(np.Par()+sigma_k*p) + r) / sigma_k; // Approximation to H*p
delta = DotProduct(p,s); // Approximation to p'*H*p
}
......@@ -497,12 +500,16 @@ NonlinOut sccngr(const NonlinParam& np, const NonlinCF& cfo)
double alpha = mu/delta; // Step size in direction p
double tmp_cf = cfo.cf(np.Par()+alpha*p); // Value of cost-function at attempted new point
// cout << "np.NIter() " << np.NIter() << ", delta = " << delta << ", mu = " << mu << ", alpha = " << alpha << endl;
/*
char fname[100];
sprintf(fname,"/Users/jesper/Desktop/gradient_%02d.txt",np.NIter());
sprintf(fname,"scg_debug_gradient_%02d.txt",np.NIter());
print_newmat(r,fname);
sprintf(fname,"/Users/jesper/Desktop/step_%02d.txt",np.NIter());
sprintf(fname,"scg_debug_step_%02d.txt",np.NIter());
ColumnVector step(p); step *= alpha;
print_newmat(step,fname);
*/
double Delta = 2.0*delta*(np.CF()-tmp_cf) / (mu*mu); // > 0 means attempted step reduced cost-function
......@@ -511,7 +518,7 @@ NonlinOut sccngr(const NonlinParam& np, const NonlinCF& cfo)
np.SetPar(np.Par() + alpha*p); // Update best set of parameters
lambda_bar = 0.0;
success = true;
if ((np.NIter()%np.NPar()) == 0) { // If npar iterations since last resetting of directions
if ((np.NIter()%np.NPar()) == 0) { // If npar iterations since last resetting of directions
r = -cfo.grad(np.Par()); // Reset search direction to negative gradient
p = r;
}
......@@ -519,6 +526,7 @@ NonlinOut sccngr(const NonlinParam& np, const NonlinCF& cfo)
ColumnVector oldr = r;
r = -cfo.grad(np.Par());
double beta = (DotProduct(r,r)-DotProduct(oldr,r)) / mu;
// cout << "np.NIter() = " << np.NIter() << ", beta = " << beta << endl;
p = r + beta*p; // New search direction
}
if (Delta > 0.75) { // If attempted step was \emph{REALLY} good
......@@ -633,7 +641,7 @@ LinOut linsrch(// Input
// First make sure that the step-length suggested
// by pdir isn't completely unreasonable.
double totstep=sqrt(DotProduct(dir,dir));
double totstep=std::sqrt(DotProduct(dir,dir));
ColumnVector pdir(dir);
if (totstep > sm) {pdir *= sm/totstep;}
......@@ -690,7 +698,7 @@ LinOut linsrch(// Input
y << f1-fp0*l1-f0 << f2-fp0*l2-f0;
ColumnVector b = X.i()*y;
// Find value for lambda that yield minimum of cubic
*lambda = (-b.element(1) + sqrt(std::pow(b.element(1),2.0) - 3.0*b.element(0)*fp0)) / (3.0*b.element(0));
*lambda = (-b.element(1) + std::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
*lambda = std::max(lmin*l1,*lambda);
*lambda = std::min(lmax*l1,*lambda);
......@@ -786,7 +794,7 @@ LinOut linmin(// Input
return(LM_CONV);
}
// Try parabolic fit, but not before third iteration
double tmp = 10.0*sqrt(MISCMATHS::EPS);
double tmp = 10.0*std::sqrt(MISCMATHS::EPS);
if (std::abs(ostep) > tol/2.0 && // If second to last step big enough
std::abs(x->first-w.first) > tmp &&
std::abs(x->first-v.first) > tmp &&
......@@ -903,7 +911,7 @@ pair<double,double> bracket(// Input
return(p_l);
}
// Let's see if a parabolic might help us
if (std::abs(l2-l1) > 10.0*sqrt(MISCMATHS::EPS)) {
if (std::abs(l2-l1) > 10.0*std::sqrt(MISCMATHS::EPS)) {
X << std::pow(l1,2.0) << l1 << std::pow(l2,2.0) << l2;
y << cf1 << cf2;
ColumnVector b = X.i()*y;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment