diff --git a/dune/tectonic/octave/duneoctave.hh b/dune/tectonic/octave/duneoctave.hh
index 902306a7f9ea5fab76fa4254aec27abb6b2cbbd6..223881ac635a8399c49be71ba82b19e406d597db 100644
--- a/dune/tectonic/octave/duneoctave.hh
+++ b/dune/tectonic/octave/duneoctave.hh
@@ -11,17 +11,29 @@
namespace Dune {
template <int dimension>
-void octaveToDune(ColumnVector const &from,
+void octaveToDune(Array<double> const &from,
typename Dune::SampleFunctional<dimension>::SmallVector &to) {
assert(from.length() == dimension);
for (size_t i = 0; i < dimension; ++i)
to[i] = from(i);
}
+template <int dimension>
+void octaveToDune(Matrix const &from,
+ typename Dune::SampleFunctional<dimension>::SmallMatrix &to) {
+ dim_vector dims = from.dims();
+ assert(dims.length() == 2);
+ assert(dims(0) == 2);
+ assert(dims(1) == 2);
+ for (size_t i = 0; i < dimension; ++i)
+ for (size_t j = 0; j < dimension; ++j)
+ to[i][j] = from(i, j);
+}
+
template <int dimension>
void duneToOctave(
typename Dune::SampleFunctional<dimension>::SmallVector const &from,
- ColumnVector &to) {
+ Array<double> &to) {
assert(to.length() == dimension);
for (size_t i = 0; i < dimension; ++i)
to(i) = from[i];
diff --git a/src/duneevaluate.cc b/src/duneevaluate.cc
index 9db4a2f59dfa84537d970bc95585cd5f76d18ccf..dd28b626b8e56c47240f8808cfb87d8e691ea12f 100644
--- a/src/duneevaluate.cc
+++ b/src/duneevaluate.cc
@@ -11,30 +11,41 @@
#include <cassert>
-DEFUN_DLD(duneevaluate, args, nargout, "duneevaluate(A)\n\
+DEFUN_DLD(duneevaluate, args, nargout, "duneevaluate(A,b,points)\n\
\n\
Evaluate x -> 1/2<Ax,x> - <b,x> + H(|x|) at each point y that is a column vector of A using DUNE\n") {
- assert(args.length() == 1);
- assert(nargout <= 1);
-
int const dim = 2;
- // {{{ The usual setup
+ assert(args.length() == 4);
+ assert(nargout <= 1);
+
typedef Dune::SampleFunctional<dim> Functional;
+
Functional::SmallMatrix A;
- A[0][0] = 3;
- A[0][1] = 1.5;
- A[1][0] = 1.5;
- A[1][1] = 4;
- Functional::SmallVector b;
- b[0] = 1;
- b[1] = 2;
+ Dune::octaveToDune<dim>(args(0).matrix_value(), A);
- Dune::SampleFunction f;
- Functional J(A, b, Dune::MyNonlinearity<dim>(f));
- // }}}
+ Functional::SmallVector b;
+ Dune::octaveToDune<dim>(args(1).vector_value(), b);
+
+ Dune::SteepFunctionNonsmooth f_steep;
+ Dune::SampleFunction f_sample;
+
+ charNDArray bar = args(2).char_array_value();
+ Dune::NiceFunction *foo;
+ switch (bar(0)) {
+ case 'S':
+ foo = &f_steep;
+ break;
+ case 's':
+ foo = &f_sample;
+ break;
+ default:
+ assert(false);
+ }
+ Dune::MyNonlinearity<dim> phi(*foo);
+ Functional J(A, b, phi);
- Matrix points(args(0).matrix_value());
+ Matrix points(args(3).matrix_value());
size_t const columns = points.columns();
RowVector ret(columns);
Functional::SmallVector our_point;
diff --git a/src/duneminimise.cc b/src/duneminimise.cc
index 39ab3dd0d07fc6faaa6eed06c364574365d93d00..fb66823b67c22ae38f15ce680c5f084f2a1ce2b4 100644
--- a/src/duneminimise.cc
+++ b/src/duneminimise.cc
@@ -11,36 +11,46 @@
#include <cassert>
-DEFUN_DLD(duneminimise, args, nargout, "duneminimise(y)\n\
+DEFUN_DLD(duneminimise, args, nargout, "duneminimise(A,b,f,y)\n\
\n\
Make a minimisation step of x -> 1/2<Ax,x> - <b,x> + H(|x|) using DUNE starting at x\n") {
- assert(args.length() == 1);
- assert(nargout <= 1);
- ColumnVector current(args(0).vector_value());
-
int const dim = 2;
- // {{{ The usual setup
+ assert(args.length() == 4);
+ assert(nargout <= 1);
+
typedef Dune::SampleFunctional<dim> Functional;
- Functional::SmallMatrix A;
- A[0][0] = 3;
- A[0][1] = 1.5;
- A[1][0] = 1.5;
- A[1][1] = 4;
- Functional::SmallVector b;
- b[0] = 1;
- b[1] = 2;
- Dune::SampleFunction f;
- Functional J(A, b, Dune::MyNonlinearity<dim>(f));
- // }}}
+ Functional::SmallMatrix A;
+ Dune::octaveToDune<dim>(args(0).matrix_value(), A);
+ Functional::SmallVector b;
+ Dune::octaveToDune<dim>(args(1).vector_value(), b);
+
+ Dune::SteepFunctionNonsmooth f_steep;
+ Dune::SampleFunction f_sample;
+
+ charNDArray bar = args(2).char_array_value();
+ Dune::NiceFunction *foo;
+ switch (bar(0)) {
+ case 'S':
+ foo = &f_steep;
+ break;
+ case 's':
+ foo = &f_sample;
+ break;
+ default:
+ assert(false);
+ }
+ Dune::MyNonlinearity<dim> phi(*foo);
+ Functional J(A, b, phi);
+
+ ColumnVector start_octave(args(3).vector_value());
Functional::SmallVector start;
- Dune::octaveToDune<dim>(current, start);
+ Dune::octaveToDune<dim>(start_octave, start);
Dune::minimise(J, start);
- Dune::duneToOctave<dim>(start, current);
-
- return octave_value(current);
+ Dune::duneToOctave<dim>(start, start_octave);
+ return octave_value(start_octave);
}
diff --git a/src/foo.m b/src/foo.m
index cf5f4d047fdcf6a085f6fefc7869b1d818b73160..121210074656cf623a5c30b8c0583f2cf4bfc980 100644
--- a/src/foo.m
+++ b/src/foo.m
@@ -5,6 +5,15 @@ close all
if exist('graphics_toolkit','file') graphics_toolkit('fltk') end
+A=[
+ 3 1.5
+ 1.5 4
+];
+b=[
+ 1
+ 2
+];
+
x = -50:1:300;
y = -125:1:50;
[X, Y] = meshgrid(x,y);
@@ -14,7 +23,7 @@ ylabel('y')
tic
for i=1:length(y)
in = [ X(i,:); Y(i,:) ];
- f(i,:) = duneevaluate(in);
+ f(i,:) = duneevaluate(A,b, 's', in);
end
clear X Y;
toc
@@ -22,23 +31,21 @@ toc
surf(x, y, f)
hold on;
-oldvec = [279; 0]; # Something random that takes a couple of iterations
-olddiff = [ 0; 0];
-
-for i = 1:10 # Something random
- newvec = duneminimise(oldvec);
- newdiff = newvec - oldvec;
- line([oldvec(1) newvec(1)], ...
- [oldvec(2) newvec(2)], ...
- [duneevaluate(oldvec) duneevaluate(newvec)], ...
- 'color', 'r');
- printf("Norm of step: %f\n", norm(newdiff));
- if (i != 1)
- printf("Angle between steps: %f degrees\n", ...
- acos( dot(olddiff,newdiff) / (norm(olddiff) * norm(newdiff)) ) / pi * 180);
- end
- oldvec=newvec;
- olddiff=newdiff;
+steps = 10;
+vecs = zeros(steps + 1, 2);
+diffs = zeros(steps, 2);
+
+vecs(1,:) = [ 0; 0]; % So that we can always calculate a difference
+vecs(2,:) = [279; 0]; % Start
+
+for j = 2:steps+1 % Something random
+ % One is code for 'SampleFunction'
+ vecs(j+1,:) = duneminimise(A,b, 's', vecs(j,:));
+ diffs(j,:) = vecs(j+1,:) - vecs(j,:);
+ line([vecs(j,1) vecs(j+1,1)], ...
+ [vecs(j,2) vecs(j+1,2)], ...
+ [duneevaluate(A,b, 's', vecs(j,:)') duneevaluate(A,b, 's', vecs(j+1,:)')], ...
+ 'color', 'y');
end
axis tight;