Partial Derivatives, Abstract Functions, and Code Generation

I'm trying to develop a large and complicated Jacobian with abstract functions, and generate C code from it. This simple example illustrates and reproduces the error I am encountering in my much larger effort:

% Generic variables
x = sym( 'x', 'real' );
z = sym( 'z', 'real' );
beta = sym( 'beta', 'real' );

% Abstract function of alpha and beta
syms Q( alpha, beta );

% Define alpha
alpha = z / x;

% Create my function.  Result is Q(z/x, beta)
functionToDifferentiate = Q( alpha, beta )

% Differentiate my function.  Result is D([1], Q)(z/x, beta)/x
partial = diff( functionToDifferentiate, z )

% Unsuccessful experiments that have been tried, at this location in this script:
%
% Experiment 1:  Throws an error
% replacement = sym( 'D__Q__Arg1', 'real' );
% partial = subs( partial, D([1], Q)(z/x, beta), replacement )     % Error on this line
%
% Experiment 2:  Throws an error
% expressionToReplace = diff( Q(alpha,beta), alpha )               % Error on this line
% replacement = sym( 'D__Q__Arg1', 'real' );
% partial = subs( partial, expressionToReplace, replacement )
%
% Experiment 3:  Doesn't throw an error, but doesn't work
% alpha__ = sym( 'alpha__', 'real' );
% expressionToReplace = diff( Q(alpha__,beta), alpha__ );
% replacement = sym( 'D__Q__Arg1', 'real' );
% partial = subs( partial, expressionToReplace, replacement );     % Doesn't replace anything
%
% Experiment 4:  Doesn't throw an error, but doesn't work
% alpha__ = sym( 'alpha__', 'real' );
% derivative = diff( Q(alpha__,beta), alpha__ );
% expressionToReplace = derivative( z/x, beta );
% replacement = sym( 'D__Q__Arg1', 'real' );
% partial = subs( partial, expressionToReplace, replacement );     % Doesn't replace anything

% Try to generate code.  This throws an error, complaining about the term
% D([1], Q)(z/x, beta).  How can I jury-rig around the error?
ccode( partial )

The problem is the hideous syntax D([1], Q)(z/x, beta), which means "the derivative of Q with respect to its first argument, evaluated at z/x and beta". I cannot generate code from it, and I cannot subs it.

Any ideas?

Thanks in advance!

[Paul]

Hi Chris,

Ideally, what should the generated code produce?

Suppose we have an even simpler case

syms alpha real
syms f(alpha)
partial = diff(f,alpha)

partial(alpha) = 

char(partial)

ans = 'diff(f(alpha), alpha)'

In this case, should the code generator (ccode(partial)) produce a line of code that calls to an unknown function, diff, that takes two arguments, where the first argument is call to another function, 'f', that takes one argument?

[Chris Rygaard]

Yes. I'd like to replace these functions (and their derivatives) with calls to undefined C functions.

Some of these 'abstract functions' are table lookups, and I can compute their derivatives from the table data. Others are functions which are computed iteratively and do not have an analytic form, but I have an expression for their derivatives.

Ideally, I'd like to replace instances of these functions (and instances of their derivatives) with C function calls. I don't need Matlab to provide the functions; I will provide the functions elsewhere.

I'm willing to replace these functions with uniquely identifiable strings. If I do that, I would load the results in to a regular-expression facility to recover the needed functions. (That is what I was attempting for Experiment #3 and #4 in the comments.)

Answer 1

[Paul]

I'm not sure the following is what you're looking for.

Not sure if it can be automated for general usage.

And, AFAICT, it might use undocumented functionality, so buyer beware.

syms x z beta real
% Abstract function of alpha and beta
syms Q( alpha, beta );
syms gamma(z)
partial = diff(Q(gamma(z),beta),z)

partial = 

c = children(partial)

c = 1×2 cell array
    {[diff(gamma(z), z)]}    {[D([1], Q)(gamma(z), beta)]}

symFunType(c{2}) % just to see how it's defined under the hood so we can use this general format

ans = "D([1], Q)"

syms dQdalpha(alpha,beta)
partial = mapSymType(partial,'D([1], Q)',@(u) dQdalpha(alpha,beta))

partial = 

partial = subs(partial,alpha,gamma)

partial = 

partial = subs(partial,gamma,z/x)

partial = 

●

syms dQ
ccode(dQ == partial)

Warning: Function 'dQdalpha' not verified to be a valid C function.

ans = '  dQ = dQdalpha(z/x,beta)/x;'

[Chris Rygaard]

Too cool. This is the solution I needed. Thanks for pointing out the symFunType, children, and mapSymType functions; they open many different paths I can take for my project.

Just to help me understand: Which parts do you think might be undocumented? Everything seems to be 'legal' Matlab.

[Walter Roberson]

@Chris Rygaard

The D([1], Q)(gamma(z), beta) part is undocumented.

[Paul]

At the time of posting I wasn't sure that the argument to diff here is documented.

partial = diff(Q(gamma(z),beta),z)

See this thread for further discussion.