Skip to content
Snippets Groups Projects
Commit efed6a9a authored by Elias Pipping's avatar Elias Pipping Committed by Elias Pipping
Browse files

Make descentDirection write to a reference

parent d9ee2155
No related branches found
No related tags found
No related merge requests found
...@@ -29,7 +29,7 @@ class SampleFunctional { ...@@ -29,7 +29,7 @@ class SampleFunctional {
return y * v + func_(v.two_norm()); // <1/2 Av - b,v> + H(|v|) return y * v + func_(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
} }
SmallVector descentDirection(const SmallVector x) const { void descentDirection(const SmallVector x, SmallVector &ret) const {
if (x == SmallVector(0.0)) { if (x == SmallVector(0.0)) {
SmallVector d = smoothGradient(x); SmallVector d = smoothGradient(x);
// Decline of the smooth part in the negative gradient direction // Decline of the smooth part in the negative gradient direction
...@@ -38,16 +38,17 @@ class SampleFunctional { ...@@ -38,16 +38,17 @@ class SampleFunctional {
func_.rightDifferential(0.0) * d.two_norm(); // TODO: is this correct? func_.rightDifferential(0.0) * d.two_norm(); // TODO: is this correct?
double combinedDecline = smoothDecline + nonlinearDecline; double combinedDecline = smoothDecline + nonlinearDecline;
return (combinedDecline < 0) ? d : SmallVector(0.0); ret = (combinedDecline < 0) ? d : SmallVector(0.0);
return;
} }
SmallVector const pg = upperGradient(x); SmallVector const pg = upperGradient(x);
SmallVector const mg = lowerGradient(x); SmallVector const mg = lowerGradient(x);
SmallVector ret;
// TODO: collinearity checks suck // TODO: collinearity checks suck
if (pg * x == pg.two_norm() * x.two_norm() && if (pg * x == pg.two_norm() * x.two_norm() &&
-(mg * x) == mg.two_norm() * x.two_norm()) { -(mg * x) == mg.two_norm() * x.two_norm()) {
return SmallVector(0); ret = SmallVector(0);
return;
} else if (pg * x >= 0 && mg * x >= 0) { } else if (pg * x >= 0 && mg * x >= 0) {
ret = pg; ret = pg;
} else if (pg * x <= 0 && mg * x <= 0) { } else if (pg * x <= 0 && mg * x <= 0) {
...@@ -56,7 +57,6 @@ class SampleFunctional { ...@@ -56,7 +57,6 @@ class SampleFunctional {
ret = project(smoothGradient(x), x); ret = project(smoothGradient(x), x);
} }
ret *= -1; ret *= -1;
return ret;
} }
SmallMatrix A; SmallMatrix A;
...@@ -97,7 +97,8 @@ template <class Functional> ...@@ -97,7 +97,8 @@ template <class Functional>
void minimise(const Functional J, const typename Functional::SmallVector x, void minimise(const Functional J, const typename Functional::SmallVector x,
typename Functional::SmallVector &corr) { typename Functional::SmallVector &corr) {
typedef typename Functional::SmallVector SmallVector; typedef typename Functional::SmallVector SmallVector;
SmallVector descDir = J.descentDirection(x); SmallVector descDir;
J.descentDirection(x, descDir);
if (descDir == SmallVector(0.0)) { if (descDir == SmallVector(0.0)) {
corr = SmallVector(0.0); corr = SmallVector(0.0);
......
...@@ -55,7 +55,8 @@ void testSampleFunction() { ...@@ -55,7 +55,8 @@ void testSampleFunction() {
SampleFunctional::SmallVector error; SampleFunctional::SmallVector error;
error[0] = -(102 - 1 + 2 / sqrt(5)); error[0] = -(102 - 1 + 2 / sqrt(5));
error[1] = -(161.5 - 2 + 4 / sqrt(5)); error[1] = -(161.5 - 2 + 4 / sqrt(5));
SampleFunctional::SmallVector returned = J.descentDirection(start); SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error -= returned; error -= returned;
assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it? assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
...@@ -97,7 +98,8 @@ void testSampleFunctionNonsmooth() { ...@@ -97,7 +98,8 @@ void testSampleFunctionNonsmooth() {
start /= (b.two_norm() + 1e-12); start /= (b.two_norm() + 1e-12);
assert(start.two_norm() < 1); assert(start.two_norm() < 1);
SampleFunctional::SmallVector returned = J.descentDirection(start); SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error[0] = -(7 / sqrt(5) - 1); error[0] = -(7 / sqrt(5) - 1);
error[1] = -(11.5 / sqrt(5) - 2); error[1] = -(11.5 / sqrt(5) - 2);
error -= returned; error -= returned;
...@@ -111,7 +113,8 @@ void testSampleFunctionNonsmooth() { ...@@ -111,7 +113,8 @@ void testSampleFunctionNonsmooth() {
start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1; start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
assert(start.two_norm() > 1); assert(start.two_norm() > 1);
SampleFunctional::SmallVector returned = J.descentDirection(start); SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error[0] = -(8 / sqrt(5) - 1); error[0] = -(8 / sqrt(5) - 1);
error[1] = -(13.5 / sqrt(5) - 2); error[1] = -(13.5 / sqrt(5) - 2);
error -= returned; error -= returned;
...@@ -148,7 +151,8 @@ void testTrivialFunction() { ...@@ -148,7 +151,8 @@ void testTrivialFunction() {
SampleFunctional::SmallVector error; SampleFunctional::SmallVector error;
error[0] = -101; error[0] = -101;
error[1] = -159.5; error[1] = -159.5;
SampleFunctional::SmallVector returned = J.descentDirection(start); SampleFunctional::SmallVector returned;
J.descentDirection(start, returned);
error -= returned; error -= returned;
assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it? assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
......
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