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

Change make to use C++11 and added profiler to levmar in nonlin

parent 2ce55936
No related branches found
No related tags found
1 merge request!5Jesper development
......@@ -7,6 +7,8 @@ PROJNAME = miscmaths
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_BOOST} -I${INC_PROB}
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB}
USRCXXFLAGS = -std=c++11
OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o minimize.o cspline.o sparse_matrix.o sparsefn.o rungekutta.o nonlin.o bfmatrix.o Simplex.o SpMatMatrices.o
LIBS = -lutils -lnewmat -lprob -lm
......@@ -21,3 +23,6 @@ quick:${OBJS} quick.o
libmiscmaths.a: ${OBJS}
${AR} -r libmiscmaths.a ${OBJS}
%.o:%.cpp
$(CXX11) -c $(CPPFLAGS) $(CXXFLAGS) $< -o $@
......@@ -16,6 +16,7 @@
#include "nonlin.h"
#include "Simplex.h"
#include "utils/fsl_isfinite.h"
#include "utils/FSLProfiler.h"
using namespace std;
using namespace NEWMAT;
......@@ -338,8 +339,12 @@ NonlinOut nonlin(const NonlinParam& p, const NonlinCF& cfo)
NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
{
static Utilities::FSLProfiler prof("_"+string(__FILE__)+"_"+string(__func__));
double total_key = prof.StartEntry("Total");
// Calculate initial values
double cost_key = prof.StartEntry("Calculate cost");
p.SetCF(cfo.cf(p.Par())); // Cost-function evaluated at current parameters
prof.EndEntry(cost_key);
bool success = true; // True if last step decreased CF
double olambda = 0.0; // How much the diagonal of H was nudged last time
ColumnVector g; // Gradient
......@@ -347,9 +352,14 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
while (p.NextIter(success)) {
if (success) { // If last attempt decreased cost-function
double grad_key = prof.StartEntry("Calculate gradient");
g = cfo.grad(p.Par()); // Gradient evaluated at current parameters
prof.EndEntry(grad_key);
double hess_key = prof.StartEntry("Calculate hessian");
H = cfo.hess(p.Par(),H); // Hessian evaluated at current parameters
prof.EndEntry(hess_key);
}
double nudge_key = prof.StartEntry("Nudge hessian");
for (int i=1; i<=p.NPar(); i++) { // Nudge it
if (p.GaussNewtonType() == LM_LM) { // If Levenberg-Marquardt
// H->AddTo(i,i,(p.Lambda()-olambda)*H->Peek(i,i));
......@@ -359,16 +369,21 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
H->AddTo(i,i,p.Lambda()-olambda);
}
}
prof.EndEntry(nudge_key);
ColumnVector step;
double ncf = 0.0;
bool inv_fail = false; // Signals failure of equation solving
double solve_key = prof.StartEntry("Solve");
try {
step = -H->SolveForx(g,SYM_POSDEF,p.EquationSolverTol(),p.EquationSolverMaxIter());
ncf = cfo.cf(p.Par()+step);
}
catch(...) {
double fail_key = prof.StartEntry("Fail");
inv_fail = true;
prof.EndEntry(fail_key);
}
prof.EndEntry(solve_key);
if (!inv_fail && (success = (ncf < p.CF()))) { // If last step successful
olambda = 0.0; // Pristine Hessian, so no need to undo old lambda
p.SetPar(p.Par()+step); // Set attempt as new parameters
......@@ -393,6 +408,8 @@ NonlinOut levmar(const NonlinParam& p, const NonlinCF& cfo)
// Getting here means we did too many iterations
p.SetStatus(NL_MAXITER);
prof.EndEntry(total_key);
return(p.Status());
}
......
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