matt-weinstein / adigator

Matlab Algorithmic Differentiation Toolbox
22 stars 5 forks source link

Bad Matlab syntax generated in Jacobian file (Matlab 2020a) #12

Open rschiemann1 opened 3 years ago

rschiemann1 commented 3 years ago

Hi,

first, thank you for the powerful tool! I am currently trying to make efficient use of it with quite a large system model, for which the Jacobian is required, so that implicit ODE solvers can be used efficiently.

However, I came accross some seemingly wrong differentiated code being created. I have a right hand side function according to the following definition:

y' = f(t, y, p)
where: 
   y' is the time derivative of the solution vector y
   t is the time
   y is the solution vector
   p is a structure holding numeric parameters.

In my code, the function f is called 'sysfun', the parameter structure p is called 'Params'.

Creating the Jacobian file with the following lines of code finishes successfully without any errors or warnings:

ay = adigatorCreateDerivInput(size(y0),'y');
adiOpts = adigatorOptions();
adiOpts.overwrite = 1;
output = adigatorGenJacFile('sysfun',{1,ay,Params},adiOpts);

In above code, y0 is a vector of initial values of the system.

However, when I then try to run the generated Jacobian file (e.g. with the command sysfun_Jac(1,y0,Params)), I get the following error message:

Error: File: sysfun_ADiGatorJac.m Line: 5898 Column: 1
Using identifier 'cada1f' as both a variable and a command is not supported. For more information, see "How MATLAB Recognizes Command Syntax".

The relevant lines of code in the Jacobian file are as follows:

cada1f1 = 0*eps_rad.f;
cada1f2 = length(TemperatureSolidInd);
cada1f3 = TemperatureSolidInd(cada1f2);
cada1f4dy = y.dy(54);
cada1f4 = y.f(cada1f3,:);
cada1f5dy = cada1f4dy;
cada1f5 = 273.15 + cada1f4;
cada1f6dy = 4.*cada1f5.^(4-1).*cada1f5dy;
cada1f6 = cada1f5^4;
cada1f7dy = sig_rad.f.*cada1f6dy;
cada1f7 = sig_rad.f*cada1f6;
cada1f8dy = -cada1f7dy;
cada1f8 = G_amb.f - cada1f7;
q_radiation.f = cada1f1*cada1f8;
%User Line: q_radiation= 0* eps_rad*(G_amb-sig_rad*(273.15 + y(TemperatureSolidInd(end),:)).^4  );
%User Line: % q_radiation=0; 
cada1f  = TemperatureSolidInd(bedInd);
cada1f dy = y.dy(Gator1Data.Index1097);
cada1f  = y.f(cada1f ,:);

In above code snippet, the line before the last line is problematic: cada1f dy = y.dy(Gator1Data.Index1097) is not valid Matlab syntax if cada1f was used as a variable before (which it was just one line above).

I hope the issues I provided can help find a solution to this problem. I greatly appreciate support. Thanks a lot!

Best regards, Robert

matt-weinstein commented 3 years ago

What is the user line of code on which it is failing? Is it by chance using a function handle?

This looks similar to the underlying issue here: https://github.com/matt-weinstein/adigator/issues/6#issuecomment-652001616 where they were using [a,b] = f(x) and f is a function handle.

rschiemann1 commented 3 years ago

The relevant user code is as follows:

q_radiation= 0 * eps_rad*(G_amb-sig_rad*(273.15 + y(TemperatureSolidInd(end),:)).^4  );  
% q_radiation=0; 

if ~isreal(y(TemperatureSolidInd(bedInd),:))
    y(TemperatureSolidInd(bedInd),:) = real(y(TemperatureSolidInd(bedInd),:));
    y(TemperatureGasInd(bedInd),:) = real(y(TemperatureGasInd(bedInd),:));
end

where: y a vector of doubles (the solution vector) TemperatureSolidInd, TemperatureGasInd and bedInd are vectors holding indices. These vectors are constant. There are no function handles present in the code.

if I get things right, these lines of code are differentiated to:

cada1f1 = 0*eps_rad.f;
cada1f2 = length(TemperatureSolidInd);
cada1f3 = TemperatureSolidInd(cada1f2);
cada1f4dy = y.dy(222);
cada1f4 = y.f(cada1f3,:);
cada1f5dy = cada1f4dy;
cada1f5 = 273.15 + cada1f4;
cada1f6dy = 4.*cada1f5.^(4-1).*cada1f5dy;
cada1f6 = cada1f5^4;
cada1f7dy = sig_rad.f.*cada1f6dy;
cada1f7 = sig_rad.f*cada1f6;
cada1f8dy = -cada1f7dy;
cada1f8 = G_amb.f - cada1f7;
q_radiation.f = cada1f1*cada1f8;
%User Line: q_radiation= 0* eps_rad*(G_amb-sig_rad*(273.15 + y(TemperatureSolidInd(end),:)).^4  );
%User Line: % q_radiation=0; %% To do:  RADIATON!!
cada1f  = TemperatureSolidInd(bedInd);
cada1f dy = y.dy(Gator1Data.Index1097); % <---------- **bad code line was generated**
cada1f  = y.f(cada1f ,:);
cadaconditional1 = ~isreal(y(TemperatureSolidInd(bedInd),:));
%User Line: cadaconditional1 = ~isreal(y(TemperatureSolidInd(bedInd),:));
    cada1f1 = TemperatureSolidInd(bedInd);
    cada1f2dy = y.dy(Gator1Data.Index1098);
    cada1f2 = y.f(cada1f1,:);
    cada1f3dy = real(cada1f2dy);
    cada1f3 = real(cada1f2);
    cada1f4 = TemperatureSolidInd(bedInd);
    cada1td1 = zeros(1010,1);
    cada1td1(Gator1Data.Index1099) = cada1f3dy;
    cada1td1(Gator1Data.Index1100) = y.dy(Gator1Data.Index1101);
    y.dy = cada1td1;
    y.f(cada1f4,:) = cada1f3;
    %User Line: y(TemperatureSolidInd(bedInd),:) = real(y(TemperatureSolidInd(bedInd),:));
    cada1f1 = TemperatureGasInd(bedInd);
    cada1f2dy = y.dy(Gator1Data.Index1102);
    cada1f2 = y.f(cada1f1,:);
    cada1f3dy = real(cada1f2dy);
    cada1f3 = real(cada1f2);
    cada1f4 = TemperatureGasInd(bedInd);
    cada1td1 = zeros(1010,1);
    cada1td1(Gator1Data.Index1103) = cada1f3dy;
    cada1td1(Gator1Data.Index1104) = y.dy(Gator1Data.Index1105);
    y.dy = cada1td1;
    y.f(cada1f4,:) = cada1f3;
    %User Line: y(TemperatureGasInd(bedInd),:) = real(y(TemperatureGasInd(bedInd),:));

The code line marked with a comment is the bad one. When I comment the isreal checks in the user code, the issue disappears (as would be expected).