## nonlinear optimization functions

Librairie C++ de calcul formel/ C++ symbolic computation library

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

parisse a écrit :
jeu. oct. 11, 2018 11:49 am
You forgot optimization.cc/.h
By the way, it would be easier to update for me (perhaps for you too) if you could put these files in your graph-theory repository (I know it's a little bit dirty but I'm not a computer scientist )
I put all Giac source files that I wrote in a single repository here (formerly the graphtheory repository):
https://github.com/marohnicluka/giac
Now everything is fetchable by the single link above (you just bookmark it somewhere .

I have just one more contribution to statistics tools before the 1.5.0 release: distribution fitting by maximum likelihood. The fitdistr command which does it is implemented in optimization.cc. The corresponding file fitdistr.tex (alongside with fitdistr.png) is in the doc directory. Short help:

Code : Tout sélectionner

``````# fitdistr
0 Lst(L),Fnc(D)
2 Returns the distribution of type D which fits most closely to the i.i.d. samples in the list L.
-1 normald
-2 poisson
-3 exponentiald
-4 geometric
-7 cauchyd
-8 weibulld
-9 sample
-10 randvector
-11 randvar
fitdistr(randvector(1000,weibulld,1/2,1),weibull)
X:=randvar(normal,stddev=9.5):; Y:=randvar(normal,stddev=1.5):; S:=sample(eval(X/Y,0),1000):; Z:=fitdistr(S,cauchy)
X:=randvar(normal,mean=5,variance=2):; S:=sample(exp(X),1000):; fitdistr(log(S),normal)``````
I also needed to fix a bug in find_randvars (vecteur.cc). The correct code is:

Code : Tout sélectionner

``````  gen find_randvars(const gen &g,gen_map &rv,GIAC_CONTEXT) {
stringstream ss;
if (g.type==_IDNT) {
if (rv.find(g)!=rv.end())
return rv[g];
ss << " var" << randvar_count;
identificateur v(ss.str().c_str());
rv[g]=v;
++randvar_count;
return v;
}
if (g.is_symb_of_sommet(at_discreted) || is_distribution(g)>0) {
ss << " tmp" << randvar_count;
identificateur t(ss.str().c_str());
ss.str("");
ss << " var" << randvar_count;
identificateur v(ss.str().c_str());
_eval(symbolic(at_sto,makesequence(g,t)),contextptr);
rv[t]=v;
++randvar_count;
return v;
}
if (g.type==_SYMB) {
gen &f=g._SYMBptr->feuille;
if (f.type==_VECT) {
vecteur F;
F.reserve(f._VECTptr->size());
for (iterateur it=f._VECTptr->begin();it!=f._VECTptr->end();++it) {
F.push_back(find_randvars(*it,rv,contextptr));
}
return symbolic(g._SYMBptr->sommet,change_subtype(F,_SEQ__VECT));
}
return symbolic(g._SYMBptr->sommet,find_randvars(f,rv,contextptr));
}
return g;
}``````
The corresponding f.type==_SYMB clause at the end of vranm (void) routine in vecteur.cc now goes as below:

Code : Tout sélectionner

``````    if (f.type==_SYMB) {
gen_map rv;
randvar_count=0;
gen e=find_randvars(f,rv,contextptr);
if (!rv.empty()) {
int nv=rv.size();
vecteur vars;
matrice R;
vars.reserve(nv);
R.reserve(nv);
for (gen_map::const_iterator it=rv.begin();it!=rv.end();++it) {
vars.push_back(it->second);
R.push_back(vranm(n,_eval(it->first,contextptr),contextptr));
}
R=mtran(R);
for (const_iterateur it=R.begin();it!=R.end();++it) {
res.push_back(_subs(makesequence(e,vars,*it),contextptr));
}
return;
}
}
``````
Furthermore, I propose excluding case 8 in distrib_nargs routine in moyal.cc. Namely, weibulld(k,lambda,theta) is evaluated to a real number, the value of Weibull PDF with k=a and lambda=b at point x=theta (not as a three-parameter Weibull). In contrast, weibulld(k,lambda) is a function and should be detected by is_distribution, but fails on the mentioned case 8. Perhaps distrib_nargs should return 2 on weibull distribution:

Code : Tout sélectionner

``````  int distrib_nargs(int nd){
switch (nd){
case 4: case 5: case 11: case 12: case 14:
return 1;
#if 0
case 8:
return 3;
#endif
default:
return 2;
}
}
``````
Finally, I improved the code for the sample command in prog.cc, now it can accept symbolically defined random variables too:

Code : Tout sélectionner

``````  gen _sample(const gen & args,GIAC_CONTEXT){
if (args.is_symb_of_sommet(at_discreted) || is_distribution(args)>0)
return _rand(args,contextptr);
if (args.type==_SYMB)
return _randvector(makesequence(1,args),contextptr)._VECTptr->front();
if (args.type!=_VECT || args._VECTptr->size()<2)
return gensizeerr(contextptr);
vecteur &argv=*args._VECTptr;
gen a=argv.front(),b=argv.back();
if (a==at_multinomial) {
if (argv.size()==3) {
vecteur v=argv;
v.insert(v.begin(),1);
return _randvector(v,contextptr)._VECTptr->at(0);
} if (argv.size()==4 && b.is_integer()) {
return _randvector(makesequence(b,at_multinomial,argv[1],argv[2]),contextptr);
}
return gensizeerr(contextptr);
}
if (args._VECTptr->size()!=2 || !is_integral(b) || b.type==_ZINT || b.val<0)
return gensizeerr(contextptr);
if (a.is_symb_of_sommet(at_discreted) || is_distribution(a)>0 || a.type==_SYMB)
return _randvector(makesequence(b,a),contextptr);
if (a.type!=_VECT)
return gensizeerr(contextptr);
return _rand(makesequence(b,a),contextptr);
}
``````
I do not plan to propose any further changes in moyal.cc, vecteur.cc and prog.cc, unless a bug emerges (which I think is unlikely).

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

I hope so, because I'm afraid I will miss some changes
For the documentation, it would also be easier for me if you keep a copy of cascmd_en.tex and aide_cas in your repository, then I can run diff and just cp if there are additions only.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

I have done so, now the cascmd_en.tex from the most recent unstable distribution is in the doc directory of my repository, with documentation about kernel density estimation and distribution fitting included.

Edit: in doc directory, I also included a modified copy of xcasmenu (en) in which the Modification submenu of Cmds/Graph Theory is added. I also suggest changing the "micro" prefix in Phys menu to "μ" (line 484 in xcasmenu), so that when user clicks on it the symbol μ is inserted at the cursor position in the cell. Could that work?

Edit 2: I had to crop some images in doc, they were at size a4. Could PNGs be automatically cropped on export from Xcas?

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

I think I tried the micro prefix, but it did not work. You can check by copying the modified xcasmenu at the right place (/usr/share/giac/doc/en in Linux).
png export is done according to the size in pixels choosen when exporting. This is done by FLTK, not me, no idea how to change that, but it's certainly better to leave that to image processing software.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

Indeed, no luck with micro prefix...

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

I added new command called bvpsolve in optimization.cc/.h, it approximates the solution of a nonlinear second-order boundary value problem.
The changed files are also aide_cas and cascmd_en.tex in doc.
I also autocropped three images with prefix random_ (also in the doc directory).

How can I pass the first derivative of the unknown function y(x) as y' to the command? dsolve supports it but I haven't figured out how. It is much easier to type than diff(y(x),x). So the first argument of bvpsolve could be entered as, e.g. 32+2x^3-y*y', not as 32+2x^3-y*diff(y(x),x).

Edit: I implemented an alternative finite-difference method which is more stable but somewhat slower. It is tried when the default method fails to converge. The changed files are optimization.cc and cascmd_en.tex.

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

I will look at that end of next week, I'm going to a conference from Saturday to Tuesday...
If you want to be able to enter y' instead of diff(y(x),x) you must declare your command to prevent evaluation, this is done with the _QUOTE_ARGUMENTS that you can see in the declaration of dsolve. Then you must take care of evaluation yourself, after replacing y' by diff(y(x),x). It's unfortunately not straightforward, you can inspire yourself from desolve_f code.

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Unstable source should be updated now.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

Hi,
here's an update to the optimization.cc:
* added commands for calculus of variations (euler_lagrange, jacobi_equation, conjugate_points and convex). The command convex is generally useful as it gives the conditions for convexity of a general twice-differentiable multivariate function.
* added kovacicsols command for solving second order ODE with rational coefficients,
* added icomp command for listing all compositions of an integer n into k nonnegative parts.
* all docs are updated accordingly.

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

I have made a small modification in optimization.cc in euler_lagrange to support t as a default time variable if x is not present in the argument, since it is very common :

Code : Tout sélectionner

``````gen _euler_lagrange(const gen &g,GIAC_CONTEXT) {
if (g.type==_STRNG && g.subtype==-1) return g;
gen L,t=identificateur("x");
vecteur u=makevecteur(identificateur("y"));
if (g.type!=_VECT) {
L=g;
if (!contains(lidnt(g),t))
t=t__IDNT_e;
} else {
...
``````

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

Thanks, I modified my source accordingly.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

Hi,
Xcas crashes with the following commandline:

Code : Tout sélectionner

``````assume(a>0);
extrema(a*x*y*z,x^2+y^2+z^2=1,[x,y,z])``````
One sometimes needs to run the second command several times until the crash happens.
Here is the backtrace from gdb:

Code : Tout sélectionner

``````#0  __GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:51
#1  0x00007ffff109342a in __GI_abort () at abort.c:89
#2  0x00007ffff10cfc00 in __libc_message (do_abort=do_abort@entry=2,
fmt=fmt@entry=0x7ffff11c4d98 "*** Error in `%s': %s: 0x%s ***\n")
at ../sysdeps/posix/libc_fatal.c:175
#3  0x00007ffff10d5fc6 in malloc_printerr (action=3,
str=0x7ffff11c4e60 "double free or corruption (out)", ptr=<optimized out>,
ar_ptr=<optimized out>) at malloc.c:5049
#4  0x00007ffff10d680e in _int_free (av=0x7ffff13f8b00 <main_arena>, p=0x7fffd4001690,
have_lock=0) at malloc.c:3905
#5  0x00007ffff72ac31e in __gnu_cxx::new_allocator<giac::paire>::deallocate (
this=0x7fffe1b8a280, __p=<optimized out>)
at /usr/include/c++/6/ext/new_allocator.h:110
#6  std::allocator_traits<std::allocator<giac::paire> >::deallocate (__a=...,
__n=<optimized out>, __p=<optimized out>)
at /usr/include/c++/6/bits/alloc_traits.h:462
#7  std::_Vector_base<giac::paire, std::allocator<giac::paire> >::_M_deallocate (
__n=<optimized out>, __p=<optimized out>, this=0x7fffe1b8a280)
at /usr/include/c++/6/bits/stl_vector.h:178
#8  std::vector<giac::paire, std::allocator<giac::paire> >::_M_emplace_back_aux<giac::paire>(giac::paire&&) (this=0x7fffe1b8a280,
__args#0=<unknown type in ./libgiac.so.0, CU 0x231e6c5, DIE 0x23f4a72>)
at /usr/include/c++/6/bits/vector.tcc:438
#9  0x00007ffff73882c9 in giac::in_zgbasis<giac::tdeg_t15> (resmod=...,
ressize=<optimized out>, G=std::vector of length 8, capacity 8 = {...},
---Type <return> to continue, or q <return> to quit---
env=env@entry=536085839, totdeg=<optimized out>, totdeg@entry=true,
pairs_reducing_to_zero=pairs_reducing_to_zero@entry=0x7fffe1b8a280,
f4buchberger_info=..., recomputeR=<optimized out>, eliminate_flag=<optimized out>,
multimodular=<optimized out>, parallel=<optimized out>) at cocoa.cc:12883
#10 0x00007ffff73893f6 in giac::zgbasis<giac::tdeg_t15> (res8=..., resmod=...,
G=std::vector of length 8, capacity 8 = {...}, env=536085839,
totdeg=totdeg@entry=true, pairs_reducing_to_zero=0x7fffe1b8a280,
pairs_reducing_to_zero@entry=0x18,
f4buchberger_info=std::vector of length 0, capacity 1024, recomputeR=false,
convertpoly8=false, eliminate_flag=false, multimodular=true, parallel=1)
at cocoa.cc:13073
#11 0x00007ffff738a47e in giac::mod_gbasis<giac::tdeg_t15> (res=...,
modularcheck=modularcheck@entry=false, zdata=false, zdata@entry=true,
rur=@0x7fffe1b8a220: -603741888, contextptr=0x0, contextptr@entry=0x555555b3a1a0,
eliminate_flag=eliminate_flag@entry=false) at cocoa.cc:13753
#12 0x00007ffff727d2ac in giac::gbasis8 (v=..., order=..., newres=...,
env=<optimized out>, modularalgo=<optimized out>, modularcheck=<optimized out>,
rur=<optimized out>, contextptr=<optimized out>, eliminate_flag=<optimized out>)
at cocoa.cc:16192
#13 0x00007ffff6d5d872 in giac::giac_gbasis (res=..., order_=...,
env=env@entry=0x7fffe1b8bb40, modularcheck=modularcheck@entry=1,
rur=@0x7fffe1b8ba5c: 0, contextptr=contextptr@entry=0x555555b3a1a0,
eliminate_flag=<optimized out>) at solve.cc:5922
#14 0x00007ffff6d5e15f in giac::gbasis (v=..., order=...,
with_cocoa=with_cocoa@entry=false, modular=1, env=env@entry=0x7fffe1b8bb40,
---Type <return> to continue, or q <return> to quit---
rur=@0x7fffe1b8ba5c: 0, contextptr=0x555555b3a1a0, eliminate_flag=false)
at solve.cc:6060
#15 0x00007ffff6d5efe8 in giac::_gbasis (args=...,
contextptr=contextptr@entry=0x555555b3a1a0) at solve.cc:6982
#16 0x00007ffff6d64fbb in giac::gsolve (eq_orig=..., var_orig=...,
complexmode=complexmode@entry=false, evalf_after=evalf_after@entry=0,
contextptr=contextptr@entry=0x555555b3a1a0) at solve.cc:6389
#17 0x00007ffff6d990a1 in giac::solve (e=..., x=...,
isolate_mode=isolate_mode@entry=0, contextptr=contextptr@entry=0x555555b3a1a0)
at solve.cc:2376
#18 0x00007ffff6d9c0b8 in giac::_solve_uncompressed (args=...,
contextptr=contextptr@entry=0x555555b3a1a0) at solve.cc:2739
#19 0x00007ffff6da0bd7 in giac::_solve (args=...,
contextptr=contextptr@entry=0x555555b3a1a0) at solve.cc:2775
#20 0x00007ffff7528731 in giac::solve2 (e_orig=..., vars_orig=...,
contextptr=contextptr@entry=0x555555b3a1a0) at optimization.cc:120
#21 0x00007ffff753c260 in giac::find_local_extrema (cpts=std::map with 0 elements,
f=..., g=..., vars=..., arr=std::vector of length 3, capacity 3 = {...}, ineq=...,
initial=..., order_size=5, contextptr=0x555555b3a1a0) at optimization.cc:1320
#22 0x00007ffff754098f in giac::_extrema (g=..., contextptr=0x555555b3a1a0)
at optimization.cc:1668
#23 0x00007ffff73a9010 in giac::unary_function_eval::operator() (
this=this@entry=0x7ffff7dd3600 <giac::__extrema>, arg=...,
context_ptr=context_ptr@entry=0x555555b3a1a0) at unary.h:200
#24 0x00007ffff6ca168b in giac::symbolic::eval (this=this@entry=0x555555b31f18, ``````

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

I can not reproduce, I checked on linux 64 bits

Code : Tout sélectionner

``````assume(a>0);
for j from 1 to 100 do extrema(a*x*y*z,x^2+y^2+z^2=1,[x,y,z]); od;
``````

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

I reinstalled giac from giac_1.5.0-37_amd64.deb, but the above command still crashes Xcas for me in Debian 9...

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Really strange. I checked again, this time with Xcas and on the VM where I compiled the debian 9 package, no crash.