ADiGator through GPOPS-II: Undefined function 'eq' for input arguments of type 'struct'.

[Friedl De Groote]

I am calling the function 'ForceEquilibriumFTstatevectorized' (see below for a copy of the code) both in my continuous and endpoint functions. Calling the function from the continuous function only worked fine. But when I call it from my endpoint function, I get the following error when ADiGator generates the derivative code:

Undefined function 'eq' for input arguments of type 'struct'.

Error in cadabinaryarraymath (line 189)
if xtemp == 0; spdyind = []; else spdyind = dyind; end

Error in cada/times (line 3)
z = cadabinaryarraymath(x,y,1,1,'times');
Error in adigatortempfunc2 (line 179)
cada1f2dv = -0.5.*cada1f1dv;

Exploring a bit further, I find that xtemp is a structure with three fields:
xtemp =

   id: []
 func: [1x1 struct]
deriv: [1x1 struct]

and these fields are:

ans =

 []

ans =

    name: '-0.5'
    size: [1 1]
zerolocs: []
   value: -0.5000

ans =

  name: []
nzlocs: []

If I replace
FMtilda3 = b13exp(-0.5num3.^2./den3.^2);
by FMtilda3 = b13*exp(-num3.^2./den3.^2/2);

I get a differt error:
Undefined function or variable "xtemp".

Error in cadabinaryarraymath (line 189)
if xtemp == 0; spdyind = []; else spdyind = dyind; end

Error in cada/times (line 3)
z = cadabinaryarraymath(x,y,1,1,'times');
Error in adigatortempfunc2 (line 197)
FMtilda3.dv = b13.f.*cada1f6dv;

What really confuses me is that I am able to use the same expressions elsewhere without getting the same errors. I tried also to copy the code in my endpoint function instead of calling it, but that did not solve the error.

Please let me know if you need more info to explore this error.

Many thanks in advance.

Friedl

CODE

function output = musdynEndpoint_Fstate(input)

q = input.phase.integral;
output.objective = q;

% input
AnkleMomentArm = input.auxdata.AnkleMomentArm_data(end-5,:);
KneeMomentArm = input.auxdata.KneeMomentArm_data(end-5,:);
HipMomentArm = input.auxdata.HipMomentArm_data(end-5,:);

Ankle_ID_data = input.auxdata.TA_ID_data(end-5,:);
Knee_ID_data = input.auxdata.TK_ID_data(end-5,:);
Hip_ID_data = input.auxdata.TH_ID_data(end-5,:);

LMT = input.auxdata.LMT_data(end-5,:);
VMT = input.auxdata.vMT_data(end-5,:);

% parameters
NMuscles = input.auxdata.NMuscles;
params = input.auxdata.params;
FMo = params(1,:);

dFend = 10*input.auxdata.controldF_end(NMuscles+input.auxdata.Ndof+1:end);

% states
Act = input.phase.finalstate(1:NMuscles);
FTtilda = input.phase.finalstate(NMuscles+1:end);

[diff, FT] = ForceEquilibriumFTstatevectorizedEndpoint(Act,FTtilda,dFend,LMT,VMT,params,input.auxdata.Ftparam,input.auxdata.Fvparam,input.auxdata.Fpparam,input.auxdata.Faparam);

A_init = input.phase.initialstate(1:NMuscles);
FTtilda_init = input.phase.initialstate(NMuscles+1:end);

% constraints - mild periodicity
perA = Act - A_init;
perl = FTtilda - FTtilda_init;

torqueScale = 1;

TA = (Ankle_ID_data - sum(FT.AnkleMomentArm,2))torqueScale; % why should this be +?
TK = (Knee_ID_data - sum(FT.KneeMomentArm,2))torqueScale;
TH = (Hip_ID_data - sum(FT.HipMomentArm,2))torqueScale;

output.eventgroup.event = [diff TA TK TH perA perl];

function [err, FT] = ForceEquilibriumFTstatevectorizedEndpoint(a,fse,dfse,lMT,vMT,params, Ftparam, Fvparam, Fpparam, Faparam)
% Try to vectorize it
FMo = ones(size(a,1),1)params(1,:);
lMo = ones(size(a,1),1)params(2,:);
lTs = ones(size(a,1),1)params(3,:);
alphao = ones(size(a,1),1)params(4,:);
vMmax = 10*lMo;

FT = fse . FMo;
lTtilda = log(5(fse + 0.25))/35 + 0.995;

lM = sqrt((lMo.sin(alphao)).^2+(lMT-lTs.lTtilda).^2);
lMtilda = lM./lMo;

% Gaussians
b11 = Faparam(1);
b21 = Faparam(2);
b31 = Faparam(3);
b41 = Faparam(4);
b12 = Faparam(5);
b22 = Faparam(6);
b32 = Faparam(7);
b42 = Faparam(8);

b13 = 0.1;
b23 = 1;
b33 = 0.5sqrt(0.5);
b43 = 0;
num3 = lMtilda-b23;
den3 = b33+b43lMtilda;
FMtilda3 = b13exp(-0.5num3.^2./den3.^2);

num1 = lMtilda-b21;
den1 = b31+b41lMtilda;
FMtilda1 = b11exp(-0.5*num1.^2./den1.^2);

num2 = lMtilda-b22;
den2 = b32+b42lMtilda;
FMtilda2 = b12exp(-0.5*num2.^2./den2.^2);

FMltilda = FMtilda1+FMtilda2+FMtilda3;

vT = lTs.dfse./(7exp(35(lTtilda - 0.995)));
cos_alpha = (lMT-lTs.lTtilda)./lM;
vM = (vMT-vT).*cos_alpha;
vMtilda = vM./vMmax;

% load Fvparam
e1 = Fvparam(1);
e2 = Fvparam(2);
e3 = Fvparam(3);
e4 = Fvparam(4);

FMvtilda = e1log((e2vMtilda+e3)+sqrt((e2vMtilda+e3).^2+1))+e4;
Fce = a.FMltilda.*FMvtilda;

% Calculate passive force length
e0 = 0.6;
kpe = 4;
t5 = exp(kpe * (lMtilda - 0.10e1) / e0);
Fpe = ((t5 - 0.10e1) - Fpparam(1)) / Fpparam(2);

FM = Fce+Fpe;%+dvMtilda;
err = FMo.FM.*cos_alpha-FT;

return

[Matthew J. Weinstein]

Hi Friedl,

That is quite odd. First a possible explaination:

I use xtemp for two different cases in that function (bad coding on my part) one if the x input is a numeric (for example, times(0.5,y)) and I use it as a holder to turn numeric x into a cada type structure (lines 55-76). I then later use it to define the function level sparsity of x, if it is known (lines 147-177) - this should redefine xtemp given that ~isempty(x.func.value) && ~isempty(y.func.value) of line 143 returns false. It seems 143 line is returning true (in which case y should not have any derivatives), but if y has no derivatives it should never get to the line that it is erroring at (dyind should be empty). This leads me to believe that y.func.value and y.deriv.nzlocs are somehow both populated at the same time, in which case there is an error in whatever function is defining y. Why this would happen, I do not know - possibly do to that b43 =0 cancelling derivatives and a function not recording it properly. I have attempted to replicate the error using

function y = myfun(x)

x2 = sqrt(x);

b13 = 0.1;
b23 = 1;
b33 = 0.5sqrt(0.5);
b43 = 0;
num3 = x2-b23;
den3 = b33+b43x2;
y = b13exp(-0.5num3.^2./den3.^2);

but it differentiates twice just fine.

Can you show me y.func, y.deriv, x.func, x.deriv at the time of the error? Also, can you show me the lines of code in the produced endpoint Grd functionn starting from b13.f = ... and leading up to the FMtilda3.dv computation? Also, are you using the latest version of the adigator code?

-Matt

[Friedl De Groote]

Hi Matt

Thanks for looking at my problem!

I am not so surprised that you were unable to replicate the error because the same expressions in my continuous funtion differentiate twice just fine as well ... I get exactly the same error if I remove the zero constant.

Below you find the additional information. I am using adigatorV0.4.3.

Friedl

At the time of the error:
y.func
name: 'cada1f1dv'
size: [9 1]
zerolocs: []
value: [9x1 double]

y.deriv
name: 'cada1f1dvdv'
nzlocs: [9x2 double]

x.func
name: '-0.5'
size: [1 1]
zerolocs: []
value: -0.5000

x.deriv
name: []
nzlocs: []

From the endpoint Grd function:
b13.f = 0.1;
%User Line: b13 = 0.1;
b23.f = 1;
%User Line: b23 = 1;
b33.f = 0.5sqrt(0.5);
%User Line: b33 = 0.5sqrt(0.5);
b43.f = 0;
%User Line: b43 = 0;
num3.dv = lMtilda.dv;
num3.f = lMtilda.f - b23.f;
%User Line: num3 = lMtilda-b23;
cada1f1 = b43.flMtilda.f;
den3.f = b33.f + cada1f1;
%User Line: den3 = b33+b43lMtilda;
cada1f1dv = 2.num3.f(:).^(2-1).num3.dv;
cada1f1 = num3.f.^2;
cada1f2dv = -0.5.cada1f1dv;
cada1f2 = -0.5cada1f1;
cada1f3 = den3.f.^2;
cada1f4dv = cada1f2dv./cada1f3(:);
cada1f4 = cada1f2./cada1f3;
cada1f5dv = exp(cada1f4(:)).cada1f4dv;
cada1f5 = exp(cada1f4);
FMtilda3.dv = b13.f.cada1f5dv;
FMtilda3.f = b13.fcada1f5;
%User Line: FMtilda3 = b13exp(-0.5*num3.^2./den3.^2);

[Friedl De Groote]

Hi Matt,

I switched to the latest version of ADiGator and it works just fine.

Thanks for your help!

Friedl

[Matthew J. Weinstein]

Excellent, no problem.