Re: Algorithmic Differentiation in Octave

[Brad Bell]

On 02/04/2017 03:29 PM, Olaf Till wrote:
... snip ...
> > Cppad is a header only package.
> Ok, I didn't expect this. But still, although a file in cppad_swig
> included 'cppad.hpp', an oct-file was compiled even at the time when
> no cppad.hpp was present on my system (in the last compilations it was
> present, of course).

cppad_swig.hpp defines an interface to
    CppAD::AD<double>
without having to include the CppAD include files.

The cppad_swig object library (*.cpp files) includes CppAD, but does not 
even know swig exists. The greatly reduces the complexity of the files 
that swig needs to process.

>
... snip ...
> No, first 'make auto_lib_octave', then 'make
> check_swig_xam_octave'. Now I've run 'make check_lib_octave' also.
>
> For compilation to succeed, I put the lines
>
> #include <octave/oct-alloc.h>
>
> #undef HAVE_HDF5
>
> #include <limits.h>
>
> into the generated 'cppad_swigOCTAVE_wrap.cxx'.
>
> > This should build all the example *.m files, the octave swig interface, and
> > run the tests.
> Yes. But the tests make Octave segfault, just as my previous manual
> tests.

This needs to be fixed, but perhaps your problems below make it not 
worth it.

> >
... snip ...
> > I thought that Octave would automatically generate matrix operations in C++
> > that used the overloaded classes (as numpy does for python) ?
> That is a misunderstanding. (In jwe's example, in the line
> 'numeric_result = builtin ('mtimes', x.data, y.data);', x.data and
> y.data refer to Octaves _native_ matrix type.)
>
> Even only to _generate_ a matrix type which contains elements of a
> type corresponding to 'a_double', one has to define the matrix type
> manually in Octaves C++ API. And, after the Octave interpreter
> commands:
>
>    ax = m_cppad.independent(x);
>
> , for the concatenation in the command:
>
>    A = [ ax(0), ax(1) ; ax(2), ax(3) ];
>
> to work, the a_double type would have to have certain methods defined
> -- without looking into the source now I think it was at least the
> resize() method. The resize() method would have to produce an object
> of our new matrix type (which also must have the resize() method,
> BTW).
>
> More directly, our newly defined matrix type could be generated from
> 'a_double_vec' by a function (written in C++), say 'm_a_double':
>
>    m_ax = m_a_double (m_cppad.independent (x));
>
>    A = [ m_ax(0), m_ax(1) ; m_ax(2), m_ax(3) ];
>
> And a method for matrix multiplication, e.g., must also be implemented
> manually for the new matrix type.
... snip ...

This is unfortunate. I had just assumed that octave had a similar 
feature to how numpy works; see
    https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html
It is my understanding that numpy implements its matrix operations using 
loops in C or C++, and not in python. That way one does not need special 
code for each type.

Perhaps Octave has a C++ template library that works for an arbitrary 
type, but it would still require, as you suggest, special C or C++ code 
for the special type.