Message Boards Message Boards

0
|
12676 Views
|
22 Replies
|
9 Total Likes
View groups...
Share
Share this post:

Problem with using NMinimize with a compiled function

Posted 9 years ago

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.
*)
POSTED BY: Tomáš Hruš
22 Replies
Posted 9 years ago

Ilian, thank you for the clear ansver!

POSTED BY: Tomáš Hruš

The arguments of NMinimize are evaluated first, which means

krit[kk, tt, f, X, Fexp]

is evaluated with symbolic kk and tt, resulting in the error messages seen. Instead, krit could be defined as

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

Another potential problem is that d can easily be zero, leading to results like

 In[5]:= krit[1, 0.05, f, X, Fexp]

 During evaluation of In[53]:= Power::infy: Infinite expression 1/0 encountered. >>

 During evaluation of In[53]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

 Out[5]= Indeterminate
POSTED BY: Ilian Gachevski

You could put in a Print statement in the function and see what values NMinimize is calling. Perhaps it's passing the function values outside the range you expect.

POSTED BY: Frank Kampas
Posted 9 years ago

Problem with using NMinimize with a compiled function

POSTED BY: Tomáš Hruš

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

POSTED BY: Frank Kampas
Posted 9 years ago

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....

Attachments:
POSTED BY: Tomáš Hruš

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

POSTED BY: Frank Kampas
Posted 9 years ago

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.

Attachments:
POSTED BY: Tomáš Hruš
Posted 9 years ago

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. >>

POSTED BY: Tomáš Hruš

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

POSTED BY: Daniel Lichtblau

agreed.

POSTED BY: Frank Kampas

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

POSTED BY: Frank Kampas

Frank, ListConvolve in Compile is not the issue.

POSTED BY: Daniel Lichtblau

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.

POSTED BY: Daniel Lichtblau
Posted 9 years ago

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..

POSTED BY: Tomáš Hruš

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.

POSTED BY: Frank Kampas
Posted 9 years ago

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}}
POSTED BY: Tomáš Hruš
Posted 9 years ago

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?

POSTED BY: Tomáš Hruš

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

POSTED BY: Frank Kampas

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").

POSTED BY: Daniel Lichtblau
Posted 9 years ago

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?

POSTED BY: Tomáš Hruš

NMinimize is expecting a function it can manipulate symbolically.

POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract