Problem with using NMinimize with a compiled function

I tried to make faster my minimization using compiled function, but my code doesnt work. See below for the complete code.

Please, tell me, where to look for the problem!

In[1]:= ClearAll["Global`*"]
SetDirectory[NotebookDirectory[]];

In[3]:= data = 
  ReadList["test66-2-0000000949-MOD2.out", Real, RecordLists -> True];
X = data[[All, 2]];
V = data[[All, 3]];
Fexp = data[[All, 4]];
data =.;
{Length[X], Length[V], Length[Fexp]}

Out[8]= {4781, 4781, 4781}

In[9]:= (* X, V and Fexp are List of the same lengt *)

In[10]:= (* Now the function crit, which will be used for NMinimize *)


crit = Compile[{{k1, _Real}, {k2, _Real}, 
   			{n1, _Real}, {n2, _Real},
   			{d1, _Integer}, {d2, _Integer},
   			{e1, _Real}, {e2, _Real},
   			{F, _Real, 1}, {X, _Real, 1}, {V, _Real, 1}},
  Module[{ kernX, kernV, Xu, Vu, xModel, critRes, Fmean, i},
   
   kernX = k1*Table[(1. - i/d1)^n1, {i, 0, d1, 1}];
   kernV = k2*Table[(1. - i/d2)^n2, {i, 0, d2, 1}];
   
   Xu = Sign[X]*Abs[X]^e1;
   Vu = Sign[V]*Abs[V]^e2;
   
   xModel = 
    ListConvolve[kernX, Xu, 1, 0] + ListConvolve[kernV, Vu, 1, 0] ;
   
   Fmean = Mean[Abs[F]];
   
   critRes = 100.*Mean[Abs[xModel - F]]/Fmean;
   critRes
   ]
  ];

Print["Test: ", 
 crit[121., 87., 1.2, 2.1, 163, 87, 0.26, 0.98, Fexp, X, V]]
        
During evaluation of In[10]:= Test: 6185.17

In[12]:= (* OK, function crit is working *)


(* but when I use it in NMinimize,  I'll get a long list of error \
messages 
in oposite with using almost the same function without compiling *)

In[13]:= res = 
 NMinimize[{crit[k1, k2, n1, n2, d1, d2, e1, e2, Fexp, X, V], 
   k1 > 50. && k2 > 10.  && 
    d1 > 10 && d2 > 10  &&
     n1 > .1 && n2 > .1 && 
    e1 > .1 && e2 > .1}, 
  {k1, k2, n1, n2, d1 \[Element] Integers, d2 \[Element] Integers , 
   e1, e2}, 
  
  MaxIterations -> 5,
  StepMonitor :> (NotebookDelete[cell]; 
    cell = PrintTemporary[{k1, k2, n1, n2, d1, d2, e1, e2}]) ]

During evaluation of In[13]:= CompiledFunction::cfsa: Argument k1 at position 1 should be a machine-size real number. >>

During evaluation of In[13]:= Table::iterb: Iterator {i$892,0,d1,1} does not have appropriate bounds. >>

During evaluation of In[13]:= Table::iterb: Iterator {i$892,0,d2,1} does not have appropriate bounds. >>

During evaluation of In[13]:= ListConvolve::kldims: The kernel k1 Table[(1. -i$892/d1)^n1,{i$892,0,d1,1}] and list {1.00136*10^-6^e1,8.00371*10^-6^e1,0.0000222445^e1,0.0000370002^e1,0.0000517631^e1,0.0000670016^e1,0.000083009^e1,0.000100025^e1,0.000116025^e1,0.000134003^e1,<<31>>,0.000730019^e1,0.00077828^e1,0.000842307^e1,0.000945801^e1,0.000983019^e1,0.000984017^e1,0.00098402^e1,0.000979975^e1,0.000978019^e1,<<4731>>} are not both non-empty lists with the same tensor rank. >>

During evaluation of In[13]:= ListConvolve::kldims: The kernel k2 Table[(1. -i$892/d2)^n2,{i$892,0,d2,1}] and list {0.000175059^e2,0.000265539^e2,0.000362456^e2,0.000368983^e2,0.000375018^e2,0.000390574^e2,0.000412792^e2,0.000412703^e2,0.000424728^e2,0.000462279^e2,<<31>>,0.00108777^e2,0.0014036^e2,0.002094^e2,0.0017589^e2,0.000477701^e2,0.000012517^e2,-0.0000505149^e2,-0.0000750124^e2,-0.000287488^e2,<<4731>>} are not both non-empty lists with the same tensor rank. >>

Out[13]= $Aborted

(* It seems, that the kernX and kernV remains unevaluated, because 
"Iterator {i,0,d1,1} does not have appropriate 
bounds." and "Argument k1 at position 1 should be a machine-size real number" -- what does this means?
 Values d1 and d2 should be integers  and sould be positive (larger then 10), 
so their values are correct. Value of k1 shoud be a real number > 50....?

In tutorial/ConstrainedOptimizationGlobalNumerical is written, 
"Constraints are generally enforced by adding penalties". May be, d1 
and d2 are not substituted as integers? But this modification of the 
code of the function crit

     kernX=k1*(Table[(1.-i/Ceiling[d1])^n1,{i,0,Ceiling[d1],1}]);
     kernV=k2*(Table[(1.-i/Ceiling[d2])^n2,{i,0,Ceiling[d2],1}]);

doesnt help, so this is not the correct way.
*)

NMinimize is expecting a function it can manipulate symbolically.

Thank you, Frank, but it is still not clear for me…

I tried a very simple compiled function - NMinimze works with it. The non-compiled version of my function - also works.

It means, that some compiled functions cannot be manipulated symbolically (and cannot by used with NMinimize – for example my function above), even his non-compiled version can be manipulated symbolically (and can by used with NMinimize) – whereas other functions can by used with NMinimize even compiled?

It would be much easier to address this if a self contained example were provided (one that does not require “test66-2-0000000949-MOD2.out”).

Does the code run if you don’t compile the module?

Frank: It was my mistake, the code without compiled function also doesnt work, you are right.


Daniel: I can say, that we are looking for constants of material with given parameters (experimental data).

The problem is to find constants of a model for given experimetal data.
In the imported file are experimental data (specimen is deformed and the force is measured): displacement X[t], displacement speed V[t] and measurd force Fexp[t].

Our hypotetsis is, that the force depends on curent values of X and V and on their histroy (implemented by convolution with kern as a decreesing function)

I am looking for constatnts e1,e2…etc, which are leads to the best fit between the force given by model – xModel (the symbol name is wrong choosen, it should be better Fmod, becuase it is force given by model) and the measured force Fexp.


So the result for me is now – should I write my own minimalization procerdure?

Do I have a bad day? This my older code is working…


In[1]:= ClearAll["Global`*"]
SetDirectory[NotebookDirectory[]];

In[3]:= data = 
  ReadList["test66-2-0000000949-MOD2.out", Real, RecordLists -> True];

In[4]:= X = data[[All, 2]];
V = data[[All, 3]];
Fexp = data[[All, 4]];
data =.;

In[8]:= odezva[k1o_, k2o_, n1o_, n2o_, d1o_Integer, d2o_Integer, e1o_,
   e2o_, Xo_, Vo_] :=
 Module[{k1 = k1o, k2 = k2o, n1 = n1o, n2 = n2o, d1 = d1o, d2 = d2o, 
   e1 = e1o, e2 = e2o, X = Xo, V = Vo, Xu, Vu, jadroX, jadroV, 
   odezva, Fmean},
  jadroX = k1*(Table[(1 - i/d1)^n1, {i, 0, d1}]);
  jadroV = k2*(Table[(1 - i/d2)^n2, {i, 0, d2}]);
  Xu = Sign[X]*Abs[X]^e1;
  Vu = Sign[V]*Abs[V]^e2;
  
  odezva = 
   ListConvolve[jadroX, Xu, 1, 0] + ListConvolve[jadroV, Vu, 1, 0] ;
  odezva
  ]

In[9]:= krit[k1o_, k2o_, n1o_, n2o_, d1o_Integer, d2o_Integer, e1o_, 
  e2o_, Fo_, Xo_, Vo_] :=
 Module[{k1 = k1o, k2 = k2o, n1 = n1o, n2 = n2o, d1 = d1o, d2 = d2o, 
   e1 = e1o, e2 = e2o, F = Fo, X = Xo, V = Vo, kriterium, o, Fmean},
  
  Fmean = Mean[Abs[F]];
  
  o = odezva[k1, k2, n1, n2, d1, d2, e1, e2, X, V];
  kriterium = 100.*Mean[Abs[o - F]]/Fmean;
  kriterium
  ]

In[10]:= krit[10000., 3000., 0.5, 0.01, 3000, 2000, 3.5, 1.8, Fexp, \
X, V]

Out[10]= 58.8654

In[11]:= res = 
 NMinimize[{krit[k1, k2, n1, n2, d1, d2, e1, e2, Fexp, X, V], 
   k1 > 150 && k2 > 10  && 
    d1 > 10 && d2 > 10  &&
     n1 > .01 && n2 > .01 && 
    e1 > .01 && e2 > .01}, 
  {k1, k2, n1, n2, d1 \[Element] Integers, d2 \[Element] Integers , 
   e1, e2}
  , MaxIterations -> 5
   ]

Speak["Calculation finished!"]

 

During evaluation of In[11]:= NMinimize::cvmit: Failed to converge to the requested accuracy or precision within 5 iterations. >>

Out[11]= {20.3068, {k1 -> 150.043, k2 -> 10.9215, n1 -> 1.94636, 
  n2 -> 0.712322, d1 -> 12, d2 -> 13, e1 -> 0.679406, e2 -> 0.95152}}

With compiled function odezva it also works

odezva = Compile[{{k1o, _Real}, {k2o, _Real}, {n1o, _Real}, {n2o, 
_Real}, {d1o, _Integer}, {d2o, _Integer}, {e1o, _Real}, {e2o, _Real}, 
{Xo, _Real, 1}, {Vo, _Real, 1}},
  Module[{k1 = k1o, k2 = k2o, n1 = n1o, n2 = n2o, d1 = d1o, d2 = d2o, 
    e1 = e1o, e2 = e2o, X = Xo, V = Vo, Xu, Vu, jadroX, jadroV, 
    odezva, Fmean},
   jadroX = k1*(Table[(1 - i/d1)^n1, {i, 0, d1}]);
   jadroV = k2*(Table[(1 - i/d2)^n2, {i, 0, d2}]);
   Xu = Sign[X]*Abs[X]^e1;
   Vu = Sign[V]*Abs[V]^e2;
   
   odezva = 
    ListConvolve[jadroX, Xu, 1, 0] + ListConvolve[jadroV, Vu, 1, 0] ;
   odezva
   ]
  ]

with a message

CompiledFunction::cfsa: Argument k1 at position 1 should be a machine-size real number. >>

but cca 25% faster..

To repeat what I earlier stated, if you provide a standalone example that is failing to work, the chances of getting a viable response go up substantially.

Given the CompiledFunction error message, I think that the compiled function is defaulting to a symbolic evaluation. I would only compile the Table functions and leave the ListConvolve outside the compiled function.

Frank, ListConvolve in Compile is not the issue.

It would be easier to figure out what the problem is if the compiled function only includes generating lists.

It would be easier to indicate what is the problem if there were a stand-alone example..

agreed.

I tried to prepare another example, hopefully better to understand. And without compiling.

Please, what is the problem on this way for looking of unknown constants?

X = Table[
  UnitTriangle[i/10], {i, -20, 20, 
   1}]; (* signal == deformation applied to a material  (here only an speculative example) *)
Fexp = Table[(1 + Abs[Sin[i/10]])*UnitTriangle[i/10], {i, -20, 20, 
   1}]; (* mesured force - reply of the material (here only an speculative example) *)
f = 15; (* sample rate *)




(* Hypotesis - the material behaviour is given by konvolution integral of the signal X. 
Konvolution kern is a trialgle with height k and leght t.

So k and t are parameters of the model. 
I am looking for their values so, that the model and experiment should be so close, as possible.

I want to find vaues k and t, which minimizes  Sum[Abs[Fexp-model[k,t]]]       *)
model[k_, t_, f_, X_] := Block[{d, kern},
  d = Floor[t*f];(* length of ListConvolve kern *)
  kern = k*Table[1 - i/d, {i, 0, d}];
  ListConvolve[kern, X, 1, 0]
  ]

krit[k_, t_, f_, X_, FF_] := N[Total[Abs[model[k, t, f, X] - FF]]]

(* 
Manipulate[
ListPlot[{model[k,t,f,X], Fexp}, \
PlotLabel\[Rule]"krit="<>ToString[krit[k,t,f,X,Fexp]]],
{{k,0.328},0,1},{{t,0.4},1/f,20/f}] 
*)(* Working as expected *)

NMinimize[
 {krit[kk, tt, f, X, Fexp], 0. < kk < 5., 0. < tt < 5.},
 {kk, tt}]

…a long list of warnigns and not a result.

First messages are here:

Table::iterb: Iterator {i,0,Floor[15 tt]} does not have appropriate bounds. >>

Table::iterb: Iterator {i,0,Floor[15 tt]} does not have appropriate bounds. >>

ListConvolve::kldims: The kernel kk Table[1-i/d,{i,0,Floor[15 tt]}] and list {0,0,0,0,0,0,0,0,0,0,0,1/10,1/5,3/10,2/5,1/2,3/5,7/10,4/5,9/10,1,9/10,4/5,7/10,3/5,1/2,2/5,3/10,1/5,1/10,0,0,0,0,0,0,0,0,0,0,0} are not both non-empty lists with the same tensor rank. >>

Did you test krit to see if it works as expected?

Yes, using construction

Manipulate[
 ListPlot[{model[k, t, f, X], Fexp}, 
  PlotLabel -> "krit=" <> ToString[krit[k, t, f, X, Fexp]]],
 {{k, 0.328}, 0, 1}, {{t, 0.4}, 1/f, 20/f}] 

can be seen, that it works OK. When model is close to Fexp, then krit is small, if modell is far away from Fexp, then krit is greater.

For example see attachement.

Looks like the iterator error message is from the function model.

Yes, it looks so.
I can use it in the manipulate command - its OK. I can sipmly write

In[6]:= krit[0.5, 0.3, f, X, Fexp]

Out[6]= 2.39701

…it works.

I can draw a Plot3d od krit as a function of k and t – it also works…

But using krit in NMinimize doesn’t works…

Problem with using NMinimize with a compiled function