Rework octave bindings

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];

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;

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);
 }

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;