Quantcast

Re: Algorithmic Differentiation in Octave

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
22 messages Options
12
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
On Tue, Jan 24, 2017 at 12:38:55PM -0700, Brad Bell wrote:

> On 01/24/2017 12:32 PM, Olaf Till wrote:
> ... snip ...
> >It seems you've been able to create Octave plugins with swig. There
> >have been problems in the past with swig and Octave. I don't know the
> >current state. But it makes me a bit uneasy to depend on swig being
> >up-to-date with Octave, since AFAIK the Octave part of swig is
> >maintained independently from Octave. But maybe we could use the
> >generated code directly for adaptions to future Octave versions, if
> >swig should not be up-to-date anymore.
> >
> >Olaf
> >
> I have kept to a very minimal swig interface to make it more robust. In
> fact, I have implemented my own version of exception handling
>     http://www.seanet.com/~bradbell/cppad_swig/error_message.htm
> because I could not get the swig version to work reliably. See
>     http://www.seanet.com/~bradbell/cppad_swig/swig_xam.i.htm
> For some swig examples I used to test what worked.
Ok, sounds good again.

Unrelated: Your license is AGPL3. Combined with Octave, in particular
with further Octave plugins used, your code may get _indirectly_
accessible over network. Does this mean that your code, if modified by
us, should always give a 'prominent notice' to the user on where the
source code is located?

Olaf

PS: We should better keep the list CCed, but I've switched to the
maintainers list now.

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Richard Crozier
On 24/01/17 20:13, Olaf Till wrote:

> On Tue, Jan 24, 2017 at 12:38:55PM -0700, Brad Bell wrote:
>> On 01/24/2017 12:32 PM, Olaf Till wrote:
>> ... snip ...
>>> It seems you've been able to create Octave plugins with swig. There
>>> have been problems in the past with swig and Octave. I don't know the
>>> current state. But it makes me a bit uneasy to depend on swig being
>>> up-to-date with Octave, since AFAIK the Octave part of swig is
>>> maintained independently from Octave. But maybe we could use the
>>> generated code directly for adaptions to future Octave versions, if
>>> swig should not be up-to-date anymore.
>>>
>>> Olaf
>>>
>> I have kept to a very minimal swig interface to make it more robust. In
>> fact, I have implemented my own version of exception handling
>>     http://www.seanet.com/~bradbell/cppad_swig/error_message.htm
>> because I could not get the swig version to work reliably. See
>>     http://www.seanet.com/~bradbell/cppad_swig/swig_xam.i.htm
>> For some swig examples I used to test what worked.
>
> Ok, sounds good again.
>
> Unrelated: Your license is AGPL3. Combined with Octave, in particular
> with further Octave plugins used, your code may get _indirectly_
> accessible over network. Does this mean that your code, if modified by
> us, should always give a 'prominent notice' to the user on where the
> source code is located?
>
> Olaf
>
> PS: We should better keep the list CCed, but I've switched to the
> maintainers list now.
>

The dev version of SWIG has a new interface method with a classdef based
wrapper which is the same for both Matlab and Octave. It is conceptually
similar to:

https://uk.mathworks.com/matlabcentral/fileexchange/38964-example-matlab-class-wrapper-for-a-c++-class

I haven't used SWIG for this, but I have used this method for wrapping
C++ classes for Octave/Matlab.

Richard

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
On 01/24/2017 01:37 PM, Richard Crozier wrote:

> On 24/01/17 20:13, Olaf Till wrote:
>> On Tue, Jan 24, 2017 at 12:38:55PM -0700, Brad Bell wrote:
>>> On 01/24/2017 12:32 PM, Olaf Till wrote:
>>> ... snip ...
>>>> It seems you've been able to create Octave plugins with swig. There
>>>> have been problems in the past with swig and Octave. I don't know the
>>>> current state. But it makes me a bit uneasy to depend on swig being
>>>> up-to-date with Octave, since AFAIK the Octave part of swig is
>>>> maintained independently from Octave. But maybe we could use the
>>>> generated code directly for adaptions to future Octave versions, if
>>>> swig should not be up-to-date anymore.
>>>>
>>>> Olaf
>>>>
>>> I have kept to a very minimal swig interface to make it more robust. In
>>> fact, I have implemented my own version of exception handling
>>> http://www.seanet.com/~bradbell/cppad_swig/error_message.htm
>>> because I could not get the swig version to work reliably. See
>>>     http://www.seanet.com/~bradbell/cppad_swig/swig_xam.i.htm
>>> For some swig examples I used to test what worked.
>>
>> Ok, sounds good again.
>>
>> Unrelated: Your license is AGPL3. Combined with Octave, in particular
>> with further Octave plugins used, your code may get _indirectly_
>> accessible over network. Does this mean that your code, if modified by
>> us, should always give a 'prominent notice' to the user on where the
>> source code is located?
>>
>> Olaf
>>
>> PS: We should better keep the list CCed, but I've switched to the
>> maintainers list now.

Is some other license preferred by Octave ?

>>
>
> The dev version of SWIG has a new interface method with a classdef
> based wrapper which is the same for both Matlab and Octave. It is
> conceptually similar to:
>
> https://uk.mathworks.com/matlabcentral/fileexchange/38964-example-matlab-class-wrapper-for-a-c++-class 
>
>
> I haven't used SWIG for this, but I have used this method for wrapping
> C++ classes for Octave/Matlab.
>
> Richard
>
Does it look the same on the Swig side. My goal to support a very
light-weight Swig interface to works for all scripting languages. Then
each individual language might have a specific package that wraps the
interface to something more convent for that language.

I have take care to make the interface efficient, and to minimize the
amount of C++ that gets processed by swig (long story with details if
you would like more).

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
On Tue, Jan 24, 2017 at 01:59:15PM -0700, Brad Bell wrote:
> Is some other license preferred by Octave ?

I'm just not familiar enough with the AGPL3 to know its implications
if the code gets _indirectly_ (via different modules) remotely
accessible. If we'd have to print messages containing the source code
location at each usage, this would be awkward. In this case, the plain
GPL3 would probably be preferrable.

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
In reply to this post by Brad Bell
After a look at the source, it seems to me that the interface to cppad
only overloads scalar operations, not matrix multiplication and not
element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
Octave). This would be a rather severe limitation for application in
Octave...

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
In reply to this post by Olaf Till-2
That does sound very awkard.  Looking on the GNU page about this, on
     /licenses/gpl-faq.en.html#AGPLProxy

The AGPL says you must make the offer to "all users". If you know that a
certain user has already been shown the offer, for the current version
of the software, you don't have to repeat it to that user again.

This does seem awkward. I will switch the license to GPL3.

On 01/24/2017 11:58 PM, Olaf Till wrote:

> On Tue, Jan 24, 2017 at 01:59:15PM -0700, Brad Bell wrote:
>> Is some other license preferred by Octave ?
> I'm just not familiar enough with the AGPL3 to know its implications
> if the code gets _indirectly_ (via different modules) remotely
> accessible. If we'd have to print messages containing the source code
> location at each usage, this would be awkward. In this case, the plain
> GPL3 would probably be preferrable.
>
> Olaf
>


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
In reply to this post by Olaf Till-2
On 01/25/2017 01:59 AM, Olaf Till wrote:
> After a look at the source, it seems to me that the interface to cppad
> only overloads scalar operations, not matrix multiplication and not
> element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
> Octave). This would be a rather severe limitation for application in
> Octave...
>
> Olaf
>
There are lots of ways to proceed here. This simplest thing to do would
be to define the matrix operations in octave and have them just be element
by element operators. (I assume this is possible in the Octave
language). Recording the function would be slow, but once it is
recorded, playback would be fast.

Perhaps a better approach, in the long run, might be to include Cppad
atomic operations in the Cppad Swig interface. It is not obvious (to me)
how to do this. For example, see the atomic matrix multiply at
     https://www.coin-or.org/CppAD/Doc/atomic_mat_mul.cpp.xml
For a paper on matrix AD calculations; see
     https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf

Note that, one only need to define the parts of the atomic operation
they are using. For example, if you do not compute the sparsity patterns
for Hessians, you do not need to include it. If you do not use forward
mode for derivative of order higher than 3, you do not need to include them.



Brad



Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

AD matrix operations

Brad Bell
In python, numpy will provide matrix operations for overloaded types. This worked for AD matrix operations; see
    http://www.seanet.com/~bradbell/pycppad/pycppad.htm
Is there a package, similar to numpy, that works for Octave ?

On 01/25/2017 06:04 AM, Brad Bell wrote:
On 01/25/2017 01:59 AM, Olaf Till wrote:
After a look at the source, it seems to me that the interface to cppad
only overloads scalar operations, not matrix multiplication and not
element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
Octave). This would be a rather severe limitation for application in
Octave...

Olaf

There are lots of ways to proceed here. This simplest thing to do would be to define the matrix operations in octave and have them just be element
by element operators. (I assume this is possible in the Octave language). Recording the function would be slow, but once it is recorded, playback would be fast.

Perhaps a better approach, in the long run, might be to include Cppad atomic operations in the Cppad Swig interface. It is not obvious (to me) how to do this. For example, see the atomic matrix multiply at
    https://www.coin-or.org/CppAD/Doc/atomic_mat_mul.cpp.xml
For a paper on matrix AD calculations; see
    https://people.maths.ox.ac.uk/gilesm/files/NA-08-01.pdf

Note that, one only need to define the parts of the atomic operation they are using. For example, if you do not compute the sparsity patterns for Hessians, you do not need to include it. If you do not use forward mode for derivative of order higher than 3, you do not need to include them.



Brad





Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: AD matrix operations

John W. Eaton
Administrator
On 01/26/2017 08:50 AM, Brad Bell wrote:
> In python, numpy will provide matrix operations for overloaded types.
> This worked for AD matrix operations; see
>     http://www.seanet.com/~bradbell/pycppad/pycppad.htm
> Is there a package, similar to numpy, that works for Octave ?

In Matlab, you can create a function like this (the @double directory
needs to be somewhere in your MATLABPATH):

   @double/mtimes.m:
   function r = mtimes (x, y)
     fprintf ('my mtimes!\n');
     r = builtin ('mtimes', x, y);
   end

to overload the builtin mtimes function that is called for the syntax "x
* y" when x and y are double precision numeric objects (that includes
sparse and complex).  You'd need a separate definition for single.

Is that what you are looking for?

The overloading feature is not yet implemented in Octave for built-in
types like double and single.

jwe


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: AD matrix operations

Brad Bell
The type a_double in Octave (and other languages)
     http://www.seanet.com/~bradbell/cppad_swig/a_double_ctor.htm
is a scalar type that will record floating point operations with the
purpose of creating a corresponding function object
     http://www.seanet.com/~bradbell/cppad_swig/a_fun_ctor.htm

In python, you can create a numpy matrix of type a_double and perform
operations, like matrix addition, multiplication, etc on the objects.
The actual matrix loops are done in C (not in python). Can one do  a
similar thing in Octave ?



On 01/26/2017 07:20 AM, John W. Eaton wrote:

> On 01/26/2017 08:50 AM, Brad Bell wrote:
>> In python, numpy will provide matrix operations for overloaded types.
>> This worked for AD matrix operations; see
>>     http://www.seanet.com/~bradbell/pycppad/pycppad.htm
>> Is there a package, similar to numpy, that works for Octave ?
>
> In Matlab, you can create a function like this (the @double directory
> needs to be somewhere in your MATLABPATH):
>
>   @double/mtimes.m:
>   function r = mtimes (x, y)
>     fprintf ('my mtimes!\n');
>     r = builtin ('mtimes', x, y);
>   end
>
> to overload the builtin mtimes function that is called for the syntax
> "x * y" when x and y are double precision numeric objects (that
> includes sparse and complex).  You'd need a separate definition for
> single.
>
> Is that what you are looking for?
>
> The overloading feature is not yet implemented in Octave for built-in
> types like double and single.
>
> jwe
>
>


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: AD matrix operations

John W. Eaton
Administrator
On 01/26/2017 09:30 AM, Brad Bell wrote:

> The type a_double in Octave (and other languages)
>     http://www.seanet.com/~bradbell/cppad_swig/a_double_ctor.htm
> is a scalar type that will record floating point operations with the
> purpose of creating a corresponding function object
>     http://www.seanet.com/~bradbell/cppad_swig/a_fun_ctor.htm
>
> In python, you can create a numpy matrix of type a_double and perform
> operations, like matrix addition, multiplication, etc on the objects.
> The actual matrix loops are done in C (not in python). Can one do  a
> similar thing in Octave ?

If you are already creating a custom object that wraps a matrix, then
you can just do something like this

   @a_double/mtimes.m
   function r = mtimes (x, y)
     ## compute the normal numeric result using Octave's builtin
     ## function for matrix multiplication.
     ## assumes "data" is the data member that contains
     ## the matrix that a_dobule wraps.
     numeric_result = builtin ('mtimes', x.data, y.data);
     ... do the other stuff that your a_double type needs to do ...
     ## package the result in an a_double object:
     r = a_double (...);
   endfunction

The above is the old-style class syntax, but you should also be able to
do this with a classdef object and methods.  If there are problems,
report the bugs.

jwe


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
In reply to this post by Olaf Till-2
On 01/24/2017 11:58 PM, Olaf Till wrote:

> On Tue, Jan 24, 2017 at 01:59:15PM -0700, Brad Bell wrote:
>> Is some other license preferred by Octave ?
> I'm just not familiar enough with the AGPL3 to know its implications
> if the code gets _indirectly_ (via different modules) remotely
> accessible. If we'd have to print messages containing the source code
> location at each usage, this would be awkward. In this case, the plain
> GPL3 would probably be preferrable.
>
> Olaf
>
I have changed the license to GPL3; see the heading Jan 26, 2017 on
     https://github.com/bradbell/cppad_swig/commits/master

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
In reply to this post by Olaf Till-2
On 01/25/2017 01:59 AM, Olaf Till wrote:
> After a look at the source, it seems to me that the interface to cppad
> only overloads scalar operations, not matrix multiplication and not
> element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
> Octave). This would be a rather severe limitation for application in
> Octave...
>
> Olaf
>
Does John's method for extending to matrix operations; see
http://lists.gnu.org/archive/html/octave-maintainers/2017-01/msg00235.html
solve your objection above ?



Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
On Fri, Jan 27, 2017 at 04:46:22AM -0700, Brad Bell wrote:

> On 01/25/2017 01:59 AM, Olaf Till wrote:
> >After a look at the source, it seems to me that the interface to cppad
> >only overloads scalar operations, not matrix multiplication and not
> >element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
> >Octave). This would be a rather severe limitation for application in
> >Octave...
> >
> >Olaf
> >
> Does John's method for extending to matrix operations; see
> http://lists.gnu.org/archive/html/octave-maintainers/2017-01/msg00235.html
> solve your objection above ?
Maybe I see it wrong, but I had seen JWEs comment more as an
explanation to you of the way(s) used in Octave to overload e.g. a
matrix multiplication, not as a complete sketch of extending your
interface to matrix operations...

(To avoid misunderstandings: Octave performs matrix multiplications
('a * b') or element-wise matrix operations (as 'a .* b') already for
its builtin matrix types, without any overloading.)

Your interface obviously already creates the type a_double (in a
different way?). Do you mean to change your code to make swig produce
m-files similar the one shown by JWE?

The core problem still is how to make your code record matrix
multiplications and element-wise matrix operations. Even if it should
be slow only at recording time, I wouldn't like a solution which
involves splitting into element operations at the Octave interpreter
level. As for the other solution mentioned by you (at the swig code
level with C-code) the details are not obvious to me either. But I
never used swig and am not familiar with it.

In any case, whether the splitting into element operations takes place
in Octave or in the swig interface, the result, after recording and
running, must be translated back into Jakobians and Hessians in the
matrix- or array-representatin of Octave. This is an additional
problem, which AFAICS can't be reasonably solved. So maybe it would be
more reasonable to support matrix operations directly in the cppad
library?

But all these matrix issues could be implemented in a simple forward
method in pure interpreted Octave code in a rather straightforward
way... and consider that probably in most large user functions the
path of computation depends on the parameters...

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
On 01/27/2017 08:41 AM, Olaf Till wrote:

> On Fri, Jan 27, 2017 at 04:46:22AM -0700, Brad Bell wrote:
>> On 01/25/2017 01:59 AM, Olaf Till wrote:
>>> After a look at the source, it seems to me that the interface to cppad
>>> only overloads scalar operations, not matrix multiplication and not
>>> element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
>>> Octave). This would be a rather severe limitation for application in
>>> Octave...
>>>
>>> Olaf
>>>
>> Does John's method for extending to matrix operations; see
>> http://lists.gnu.org/archive/html/octave-maintainers/2017-01/msg00235.html
>> solve your objection above ?
> Maybe I see it wrong, but I had seen JWEs comment more as an
> explanation to you of the way(s) used in Octave to overload e.g. a
> matrix multiplication, not as a complete sketch of extending your
> interface to matrix operations...
>
> (To avoid misunderstandings: Octave performs matrix multiplications
> ('a * b') or element-wise matrix operations (as 'a .* b') already for
> its builtin matrix types, without any overloading.)
>
> Your interface obviously already creates the type a_double (in a
> different way?). Do you mean to change your code to make swig produce
> m-files similar the one shown by JWE?
It is my intention to have a developer, for each target language, who
creates a natural interface for that target language. See the heading
Purpose on
     http://www.seanet.com/~bradbell/cppad_swig/cppad_swig.htm
I would provide any needed support and additional connections to the
Cppad in order to make this possible and fast.
>
> The core problem still is how to make your code record matrix
> multiplications and element-wise matrix operations. Even if it should
> be slow only at recording time, I wouldn't like a solution which
> involves splitting into element operations at the Octave interpreter
> level. As for the other solution mentioned by you (at the swig code
> level with C-code) the details are not obvious to me either. But I
> never used swig and am not familiar with it.
Given that Octave, like python with numpy, has a way to create matrix
operations from scalar ones, I do not think the recording will be slow.

The solutions, with atomic functions, and matrix AD operations are much
more complicated. They may give some benefit in the long run, but I do
not think they are a good idea for a first implementation.

>
> In any case, whether the splitting into element operations takes place
> in Octave or in the swig interface, the result, after recording and
> running, must be translated back into Jakobians and Hessians in the
> matrix- or array-representatin of Octave. This is an additional
> problem, which AFAICS can't be reasonably solved. So maybe it would be
> more reasonable to support matrix operations directly in the cppad
> library?
I do not see any problem here. The Jacobians and Hessians are returned
as simple vectors and special purpose *.m files  could put them in any
desired form.

>
> But all these matrix issues could be implemented in a simple forward
> method in pure interpreted Octave code in a rather straightforward
> way... and consider that probably in most large user functions the
> path of computation depends on the parameters...
>
> Olaf
Forward, tapeless mode, has advantages. In the simple case, of only
needing first order derivatives, the complex step method provides this
with no extra work (in Octave) and is probably as fast as possible.

On the other hand, optimizing the tape, computing sparsity patterns and
sparse derivatives, reverse mode, and just function evaluation with C++
speed, has advantages also. For example, reverse mode can compute the
derivative of a scalar value function with respect to an arbitrary
number of variables in a small multiple of the number of floating
operations for the original function.



Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
On Fri, Jan 27, 2017 at 09:05:35AM -0700, Brad Bell wrote:
> Given that Octave, like python with numpy, has a way to create matrix
> operations from scalar ones, I do not think the recording will be slow.

I do not understand, could you explain? What do you mean by creating
matrix operations from scalar ones? What is the relation to the speed
of recording? (I know what Octave can do, but I don't know what you
mean.)

> >In any case, whether the splitting into element operations takes place
> >in Octave or in the swig interface, the result, after recording and
> >running, must be translated back into Jakobians and Hessians in the
> >matrix- or array-representatin of Octave. This is an additional
> >problem, which AFAICS can't be reasonably solved. So maybe it would be
> >more reasonable to support matrix operations directly in the cppad
> >library?
> I do not see any problem here. The Jacobians and Hessians are returned as
> simple vectors and special purpose *.m files  could put them in any desired
> form.
Still can't see it... Say we have matrices A:=[a11, a12; a21, a22] and
B:=[b11, b12; b21, b22], both A and B somehow depending on a vector of
a_double. If we overload matrix multiplication (normally done symply
by typing A*B in Octave), to make it recordable, into scalar
operations, we compute C:=A*B as c11=a11*b11+a12*b21,
c12=a11*b12+a12*b22, c21=..., c22=... . We have to return the zero
order result, so we concatenate c11, ..., c22 into a Matrix C
e.g. with the Octave command C=[c11, c12; c21, c22], and let the
overloaded method return C. Since all computed c.. should be of type
a_vector, the operation should have been recorded by cppad. But if we
run the record to get the derivatives, in which form will they be
returned? cppad has only records for the four single values c11, ...,
c22, it doesn't know that they belong together ... or what?

If you correct the above, please be detailed, so that I can understand
it although I'm not familiar with cppad.

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell
On 01/27/2017 01:20 PM, Olaf Till wrote:
... snip ...

> Still can't see it... Say we have matrices A:=[a11, a12; a21, a22] and
> B:=[b11, b12; b21, b22], both A and B somehow depending on a vector of
> a_double. If we overload matrix multiplication (normally done symply
> by typing A*B in Octave), to make it recordable, into scalar
> operations, we compute C:=A*B as c11=a11*b11+a12*b21,
> c12=a11*b12+a12*b22, c21=..., c22=... . We have to return the zero
> order result, so we concatenate c11, ..., c22 into a Matrix C
> e.g. with the Octave command C=[c11, c12; c21, c22], and let the
> overloaded method return C. Since all computed c.. should be of type
> a_vector, the operation should have been recorded by cppad. But if we
> run the record to get the derivatives, in which form will they be
> returned? cppad has only records for the four single values c11, ...,
> c22, it doesn't know that they belong together ... or what?
>
> If you correct the above, please be detailed, so that I can understand
> it although I'm not familiar with cppad.
>
> Olaf
>

There are two separate steps. One is the recording of function
evaluation to create an AD function object.
The second is using the function object to evaluate derivatives (of any
order including order zero which is function values).

The function gets recorded by:
1. first declaring the independent variables, say the vector ax. I use
ax to denote the fact that its elements have type a_double.
2. Performing Operations to eventually get the vector of dependent
varialbes ay.
3. Stopping the recording and declaring that the function maps ax -> ay.
For example, suppose that we have overloaded matrix multiplication for
a_double and we enter

ax = m_cppad.independent(x)
A = [ ax(0), ax(1) ; ax(2), ax(3) ]
B = A * A
ay=m_cppad.vec_a_double(4)
ay(0)=B(1,1)
ay(1)=B(1,2)
ay(2)=B(2,1)
ay(3)=B(2,2)
af = m_cppad.a_fun(ax, ay)

Then the function mapping x to y (which uses matrix multiplication in
its definition) would have been recorded in the function object af.
Does this help ?


Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
On Fri, Jan 27, 2017 at 01:47:49PM -0700, Brad Bell wrote:

> <...snip...>
> The function gets recorded by:
> 1. first declaring the independent variables, say the vector ax. I use ax to
> denote the fact that its elements have type a_double.
> 2. Performing Operations to eventually get the vector of dependent varialbes
> ay.
> 3. Stopping the recording and declaring that the function maps ax -> ay.
> For example, suppose that we have overloaded matrix multiplication for
> a_double and we enter
>
> ax = m_cppad.independent(x)
> A = [ ax(0), ax(1) ; ax(2), ax(3) ]
> B = A * A
> ay=m_cppad.vec_a_double(4)
> ay(0)=B(1,1)
> ay(1)=B(1,2)
> ay(2)=B(2,1)
> ay(3)=B(2,2)
> af = m_cppad.a_fun(ax, ay)
>
> Then the function mapping x to y (which uses matrix multiplication in its
> definition) would have been recorded in the function object af.
> Does this help ?
Yes, think I got it, thanks. The declaration of mapping at the end was
the clue.

For finding an overall way to fit this into Octave, I still have to
find out about some issues in your interface which are not clear to
me. But before asking, I'd like to give it a try, making it compile
and experimenting with it a bit. Since I have currently some other
things on my schedule, I guess it'll take a week.

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Olaf Till-2
On Sat, Jan 28, 2017 at 09:42:43AM +0100, Olaf Till wrote:
> For finding an overall way to fit this into Octave, I still have to
> find out about some issues in your interface which are not clear to
> me. But before asking, I'd like to give it a try, making it compile
> and experimenting with it a bit. Since I have currently some other
> things on my schedule, I guess it'll take a week.

I havn't come very far. The cppad library, compiled according to the
instructions at

https://www.coin-or.org/CppAD/Doc/cmake.htm

with

cd build
cmake -G "Unix Makefiles" -D cppad_prefix=/usr/local ..

with 'make install' only installs headers, no library. And I find no
library even under the build/directory.

swig-cppad compiled an oct-file, after fiddling with Octaves include
files. (BTW: Instead of `octave-config -p INCLUDEDIR`, use `mkoctfile
-p INCFLAGS`, yielding 2 include dirs, with '-I' already prepended.)

swig-cppad compiled this although no cppad library was present
anywhere, and not even the cppad headers were present at this time
(and I don't find them in swig-cppad either)! The resulting oct-file,
according to 'ldd', referenced nothing which looked like a cppad
library, of course. And it segfaulted at the time of loading into
Octave.

So I can only give 'theoretical' notes now.

Let's assume swig-cppad worked as expected. Then, if we took it
unchanged, we'd have to define a wrapper class around it. Doing this
in iterpreted Octave code would yield matrix operations which are too
slow.

I think the best bet would be to generate an Octave type at the c++
level (similar to what swig-cppad currently does) which interfaces to
cppad and already contains overloading of matrix operations.

If the latter is done, I would not like to do it with swig, but rather
try it directly with Octaves API. And the build system of swig-cppad
would have to be changed anyway, IMHO, to make it distributable as an
Octave package.

This would be rather a long term project for me (but maybe others
would be faster). I have difficulties in deducing from swig-cppad the
interaction with the cppad interface. I'd have to learn how to
interface with cppad from the cppad documentation.

Problems are that cppad currently installs no library for me (see
above), and that cppad doesn't seem to be available in current Debian
(i.e. only the headers are available there).

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc (836 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Algorithmic Differentiation in Octave

Brad Bell

Did your run the following command ?
     make check_lib_octave
(see same question below)

On 02/04/2017 05:48 AM, Olaf Till wrote:

> On Sat, Jan 28, 2017 at 09:42:43AM +0100, Olaf Till wrote:
>> For finding an overall way to fit this into Octave, I still have to
>> find out about some issues in your interface which are not clear to
>> me. But before asking, I'd like to give it a try, making it compile
>> and experimenting with it a bit. Since I have currently some other
>> things on my schedule, I guess it'll take a week.
> I havn't come very far. The cppad library, compiled according to the
> instructions at
>
> https://www.coin-or.org/CppAD/Doc/cmake.htm
>
> with
>
> cd build
> cmake -G "Unix Makefiles" -D cppad_prefix=/usr/local ..
>
> with 'make install' only installs headers, no library. And I find no
> library even under the build/directory.
Cppad is a header only package. There is are a few library's that will
install if you choose certain options, but the headers are all that is
needed by cppad_swig. You can run a check of CppAD for your particular
machine with the following command
     make check
see
     https://www.coin-or.org/CppAD/Doc/cmake_check.xml

Note that the cppad_swig package does provide a C++ AD object library;
e.g., see
http://www.seanet.com/~bradbell/cppad_swig/a_fun_jacobian_xam.cpp.htm

>
> swig-cppad compiled an oct-file, after fiddling with Octaves include
> files. (BTW: Instead of `octave-config -p INCLUDEDIR`, use `mkoctfile
> -p INCFLAGS`, yielding 2 include dirs, with '-I' already prepended.
Thanks, I will give this a try.
>
> swig-cppad compiled this although no cppad library was present
> anywhere, and not even the cppad headers were present at this time
> (and I don't find them in swig-cppad either)! The resulting oct-file,
> according to 'ldd', referenced nothing which looked like a cppad
> library, of course. And it segfaulted at the time of loading into
> Octave.
See the heading Testing on
     http://www.seanet.com/~bradbell/cppad_swig/cppad_swig.htm
When you run
     bin/run_cmake.sh
there will a list of available make commands (for your configuration).

Did your run the following command ?
     make check_lib_octave
This should build all the example *.m files, the octave swig interface,
and run the tests.
Setting verbose to true in bin/run_cmake.sh (before running it) may help.

No install is yet available because I figure that should be wrapped with
the special code (like matrix multiply in Octaves case) for each language.
>
> So I can only give 'theoretical' notes now.
>
> Let's assume swig-cppad worked as expected. Then, if we took it
> unchanged, we'd have to define a wrapper class around it. Doing this
> in iterpreted Octave code would yield matrix operations which are too
> slow.
I thought that Octave would automatically generate matrix operations in
C++ that used the overloaded classes (as numpy does for python) ?

> I think the best bet would be to generate an Octave type at the c++
> level (similar to what swig-cppad currently does) which interfaces to
> cppad and already contains overloading of matrix operations.
>
> If the latter is done, I would not like to do it with swig, but rather
> try it directly with Octaves API. And the build system of swig-cppad
> would have to be changed anyway, IMHO, to make it distributable as an
> Octave package.ld
I defined a very lightweight (and I think efficient) interface from the
Cppad Swig library to any language. I figured that this would greatly
reduce the work involved in make taped AD available in not only
available in python, but other languages as well.
>
> This would be rather a long term project for me (but maybe others
> would be faster). I have difficulties in deducing from swig-cppad the
> interaction with the cppad interface. I'd have to learn how to
> interface with cppad from the cppad documentation.
>
> Problems are that cppad currently installs no library for me (see
> above), and that cppad doesn't seem to be available in current Debian
> (i.e. only the headers are available there).
See
    https://tracker.debian.org/pkg/cppad
>
> Olaf
>


12
Loading...