From 6b86cdc1e4d4e9793b05dbea925050eb61b65f36 Mon Sep 17 00:00:00 2001 From: Pariterre Date: Mon, 16 Dec 2024 12:28:14 -0500 Subject: [PATCH] Starting an example with torch --- .../_l4c_generated/adj1_l4casadi_f.pt | Bin 0 -> 6326 bytes .../_l4c_generated/jac_adj1_l4casadi_f.pt | Bin 0 -> 11699 bytes .../_l4c_generated/jac_l4casadi_f.pt | Bin 0 -> 8642 bytes .../_l4c_generated/l4casadi_f.cpp | 111 ++ .../_l4c_generated/l4casadi_f.pt | Bin 0 -> 6342 bytes .../custom_non_casadi_dynamics.py | 513 ++------ bioptim/models/torch/torch_model.py | 1064 +---------------- 7 files changed, 289 insertions(+), 1399 deletions(-) create mode 100644 bioptim/examples/getting_started/_l4c_generated/adj1_l4casadi_f.pt create mode 100644 bioptim/examples/getting_started/_l4c_generated/jac_adj1_l4casadi_f.pt create mode 100644 bioptim/examples/getting_started/_l4c_generated/jac_l4casadi_f.pt create mode 100644 bioptim/examples/getting_started/_l4c_generated/l4casadi_f.cpp create mode 100644 bioptim/examples/getting_started/_l4c_generated/l4casadi_f.pt diff --git a/bioptim/examples/getting_started/_l4c_generated/adj1_l4casadi_f.pt b/bioptim/examples/getting_started/_l4c_generated/adj1_l4casadi_f.pt new file mode 100644 index 0000000000000000000000000000000000000000..63d2c1b8d6c99eb002913b5b5c0d473771b3ee40 GIT binary patch literal 6326 zcmb_g2|SeD_a8e$S<9~Ml$bHrn31vWF${$)$xc~HwiZid$&#&-d7c!J zrI01bk|gqKLEg0a&!o4Pdh7rDFLyq7?#w*re$RK#z31F>k13Xk83JKvhy2upAW#V2 z%?pkrs=49=@NNX0hl(4Xf>-wQCK_N^A@*NeB%>=1MjeT}!?@fIg33@lBkrJSh z`k9pxg!v~@t|T{imA`f1RKAi>@i?loo&@drHcOalhGw)9TP)F@t3PYhE#k8XMr`D2uLKY5LNIA;@sXX;$kw!L;HEE5 zFPu1PM>V6_2@03FnBSb`S&B}l#7i_5z%beZB_%65&$WH?m4Z}_`gKEPWyHqtAFCAP$!+DZ=sj_ayQ}A~ZO>l#+Rl~TDeW+_))XA> z&M;2Gb7+Y2rEE80%rfpfdpdY`LEWlOBzN`kC76UHqZ)`&Emjr2Y5Xdxp4kA=Du# z=c>IX_h5Q^&o=EkJydA@mm*$~B%iDLd^!Uy{BqvL7jhJHU9B9Obr=TYISDN2=YAIx zL5>WC~U`tX_MR`R4$uknAc=Nshy~3BZb)gSV#RT{N2=^+KA8Yhm*lXm)yZ8 zc9&ba`)`=GH}*e|?msCR_&l}o?S7}mA}?u}WRPn@AdYcFC1GH5D$!eA@`!i4lg^Yc z)!!yD`ozr9KrDX^=eexy3$QRo`c7-^;S_((0D&m7{cfjOKzRQDcAB!AyG!6<+*i-x z`sO)n{*z&QfGFCwrC7ndhVXgh?CSi}b{ZHkWh!g^N1iH>$x zR<2F;r9WYv&u@KX7_pg<+n~K9R5v4wTpAQ%2vO|071?WHD#H+>Ad;t_m>0*Ike2P_ ze?K*Z849_#<)afnq|3-3icy2<9FoVXjwYe!af|z>-!vDY?(KVbvysGl)X}F+*h(o}(iB7Y~oE z3R(t~pbuUR$=t*RP*R@FEgbVk{<-SkIsD1GOXE2+@bu4P*_a1o4|=OI2oYv7QpCQ| z_bhh8^**7u6rCKB4_s?gd24jBMPB-@vt!W5JzH&$$X8!)Bs7LLC+;a|KC3-@sb(kq zo^xX}2eoiWZdCqt{M5rhJ@YZ%%H_{m&4JNMxUma{eHXMQY;w*VB$)F)+djl^&y7q= zU()L^3vYj+*IpFzl0$ypS9#E-s{vg%*0dv}IcpyY?Ot5tRBPQVEZJ*l0@OuY8 z%Am#Jr#Qm$-}MM89>iTjhNUem>pxB_eg96s@iNeJ9q@CFFc=e@y*!sx*QDY&-}d%c z+ezm^&f;;iS5akwu{qS9mWj9CNTu$I)!lyDIt=2!?8zeSibS?T1SWn!)Nyy%@|4L8TqF8v$ z$kUqCD~EI`%O`yHR=S_I5O0mSFIY32J^8S~iNF7TrDKDDy!CFJGIz>6fi>cfS4JVY z*=`#f4ws{~rz=#&Q)QAil}oIe6&2WcovaWGdvCq)>gLJt(eAzq>xu!QbozU(zsTUh z)h7)X!?CvVB@IEjltJxf<6&gG{>nJ><* z`gk3q0CnSg-dWY&Grq;zKy0bAR~BB~*vQK>=C)Ip>Vk+zLlwTGA|IRCI-WMr^V)vi zqSe@D3-Vw9_TmHnzqp(;Y~UX$RsyfvFft8=!MgL7s=F#KS?xRgg%Q_EQL*Z(*-pLIMD;Qaa?AC| z@w{h;-|4LQ`xd%YzVMxMc;Ra2USU?Q@b;a(e#*m6_0LIGh52nei@XxTW0O|g9IM6% zcjtfdy}Zgj#}$cU$DpIPb1}i#bq>kv9TMbfh9%!pNN(0L zzUBLfDHPAQ+#Ie2VSCCW7Ax&1aOybYgIi=R2KlFUS}YcMx7E(YzSs^$@SM=n zwuLKbzN;5u0n+iWECN#VjvO6CCq82G2p7{kp%yL99Co|*i%kp4y?8K%NBsz&8vj)B z-6%$l&aTkUGrPPD*gxz>ASCMqq^gthLO*qUlv>`6(8SD0sK=EoZC7z>k$^4Y*~5&x z@ljv8YrpJkLE)z{Liw=cyisx=wFjSlW(-Qt2DuC4`o&~to`;{>lJsb2*D)~zHdgjz z-X)9hmk;#*JmBFG$0H#HGwHzvGog;y<6y9+$wNi>I0Q{vqjwLMJcWVGiQ2 z6PmPk-g>oMH|Ri1Wt`Si!xm?kY0*0~K3=c#i=JH`6TCWL657~E9n`%ELWuP)@q9`U1VZqGT)U}AfPVb;UN>pUbg-ogCt%z0CQ)8E>&2lOkOdfU<4D%XY2+K1~Zq zhtLZki!lQ=Vs-}g;jbDLeklgvXa^h} z)E)}Z^-7%Li&CKJr1#(H<^GLcRRav?`gJ|U#|&a?@1g31lfg0VDp1lg6O49i05PO$ z5ML{WcDVlm6*XT8_LEe=#x14*18^4A91`E;-K$Dn9SK3TKe0lI>?@*L)ouZ*kzuIY z@NHljVJDEDn~priwW3=o7zL#7vZqch*i)sRU8YV{-@~;JB z-;1C#Bfg-TIx@gWtT{MIj7Oyz=773xX<+Q76p)qqIw&HX0xEO~p!r%GP#K6aFrYXP zWa{fc4o92^*Nn=uM;=Z{C6AmkmYj_q$KcF=hu|wpk$Q>r9qR z(@3RVQ0V=4>-;Bv_dv##i2v*UQ~P26t!@=VS8h)PC92E8305a`=@)PjLK-CZdJl5YU?pud~Dzq^|r!Ic6fvnrZX%*kw?G<}{7o@;n; zGP@@;ZL%4egO2kRg_7X8CJNmks!HahiMr85xxR{iQwjYpzG;Ivg1q?~GWXBOJR3yS z$hX@F7r?^k@|K7hlZ^etzb5c-_HYV`hV0>TiDp9r`F zQ`|}PMKr*u{!f+Ds1OMh3elZTAqM$On3blU-hWL&4kGsxg*6)lxcd|EL_!FjzOirw zHv^2+f4~P0*Va;3hihtTX&^MzRT1iHFm*UwQ&meHfk4nkZMd4M7EE26&d~I+;n{!U z2LmIG==aLw`tqL;{lSLzt2(49mX7jIsc?O1$7lmW2JPFADF5ri;g4$9CtMaA)UML# z{lBUGr;4;j;2-O)&q!!V-?}!z`ELZiXL&yYtxqvF0`1xKJ0SRvK9b3#9rZ z(E3AXBTydq?-(I|1X_P8Yy{Hc`5lnjk3j2#jrclKLm(5s1EL=(KSHe!8ylhe_1rEha3n09?f;rjOJ@9_ye kTGHB4fwhq~kA6x*EWZT+`f+HAWuXT&2p8=|mtT|pFCa!@m;e9( literal 0 HcmV?d00001 diff --git a/bioptim/examples/getting_started/_l4c_generated/jac_adj1_l4casadi_f.pt b/bioptim/examples/getting_started/_l4c_generated/jac_adj1_l4casadi_f.pt new file mode 100644 index 0000000000000000000000000000000000000000..96fc6b38209a3137061caf31bcbc7b5695187abb GIT binary patch literal 11699 zcmb_?1z6O}*Z2mFsh1waMZLoFasOM50reT z`CY&f>SW_!1rY*59@)D9T>)%?{;Zznt_a*3Dro2qX;qDKHzaPW#tkW=8&dz1l^Q_u zFG*QATUzn_t_Q;Ng9wkc7thZ`ZrcCD^miiMF5dFe;y0uC|1a4F@0?ae@jA}*(f!Q< zJ$Q|VIPnOmNCS@i;X_>L=sg95gNjW?DoxeSiHUtHpUM2u&50kRCJ*f7YA;Rdx+W=M9~Cwmj!TK0kIqI&hZu>OniZ$ z&qRFB?7sTk%F{+V20WvzV{ceHXw7H(xW7byM2t+`Dbgt+fm8Z~q@ANv*X3*P;gNhL z^(N&-(N{Y&-*g!z4D)AQOZ2$G*L+yz2Ze>rk4qV`gG;VO)m9HRMyGt>FFj5Y8X6={ zXqyAl1|_RyZ+p*z@kjM=ge&r-qJVpvhikbd1yt$*AN>hXd7KpCYdDNe9}3-k#-s_k!c3ty)h$ zqpT1q!F!IzWiw5R`Z&SxZyi^z-`D5A#tnk_@1s4!9xC-UumG5npj?iabC!8m89wN7 z;~Xr)O4ogUg>jTxCy*&BfQi^w&J&4GFNC{y9aPmSfMrtQ61gf~vPewJ|Cn!pdVh@_ z7rO8s>5hq7h7`U1MOdc;A@`K(Ko>kfRo=T%!=yC&0BLqWCkW{ud> zdJc>WPhxB^FBeiOb=Vc^0h_gHNY$lQvNj~H#QBrU%J*A8l_(O7OQWqI_%A*lKdWer z%ussu30EE93GWvG)oBJ~Ir7U{FN!5`NO#E6~`3w@u_R=C>DPEOn0!cXCEkr!C(amx4 zUduXvmgP|m3WO|tyPm_*cniWha%j=W^3ojO5~M19Ky}oF)H~deGd#p__kjNH5yRap zE~cZ3?56p1sj%8=4kNqnn5Liwl+GIXb8I2pth2Ov^1G=rk=8VaPUQ34_s;%~1aPWS=1kTOq)13}fw{EM>cBFWHk=&tyLJ6vjdDGVo zdg(CZ9qVwP&7GMo!_itJ+tE*im3WpbcLxklM_R5D>jAUPAoeJs7{DmA{(T)jT>5m*%2G; zj*pb=Fne3L>|=-bwVSe-T?w!hR}90H{jow1?nbS&Hr&I9>~B?%4QqVLY_0wF4K6@# zb0xHDwv$ISJY6+RKYcKSMQ87-ZtZV}M;`SxIcX?)HKSc>ey6DSsg&4k^SEuht4#hu zdH)QH7UKqQ9nMp_<>3tkF|Xu1VOnN|3Ypa0NIcu+LK;&OjNTQ!r2G@P6X%GZP2_yU za~?FMJzNdOLMs^e`^sCmWh^=?_bUByC239R>I7|mxHJ4FW}~6A!$P_-Pyhfg;(vz8xL^eu4EDNS%{O(tM@UK;z7Lzjdtf z$-=ge?w85!y_wMK54qp=7a@~p{<e%ny=`0)7A`s24p)Co%c=EModu#C-nmpdtm>naL*xzRX2vMWkDsG3jrdhb~2U7MCp#~xO@ zyrWv8HI#9zU#&5(Uwzg(UoiS+c~-1VA@tO%_IY2L>w8Y<#oBmwzMg6+)*54jiq&8& z3cJv>)4OH(%PwAiEn_Q+nFOY9eOF15Ub&kZC>yrqnr5*z`CdfP@4eW`d0K3nEIqEv zN;uc&Y^(Yw5)^Hs6Spph7xY%bq=QNdTR1Y!LdCS_56m?VCQDb^Uin7toJsTL?kIZ1 zuS8N?(!?*8Hrn`>BrK5?R=KDn8a_T6e{GL&^1!xYG@{mMdg7IkrYZ@co$ISNh*jU3 zNJ^Co5A9}^BUO~W)23b0Yy0CW^B-lA5c`j*0hzNNuvPkG3)eR=S*~7q8PzdNKRNrP zK1Gu08KwRu`hhQs#QMDy%F|x=4z0cuCW1D?ml)A$zJnh!FlH%mWp2}HYdNzQyV*1cxB%L<6@{*`owm>Iak;s^#*6$u9`tGI3bU2^x`JrY{Zr&BgQqSi7ponn*P*k8Tp@0P*{!L}w?pgU-@Y=xuhWA1qdWXG$lOR*q z2a+u}2bFgJdlJNn%@L7RteSr9HYN~^d{WPTXp{<$ws@kvAOZWzRoEMIt13iZA!MUJOaAbM*jow@n@T+- zZdd*f8dh&9#JpoB?bq3x&NpQlHNxu4EyeZwMJtRF>c`3-HmIo&ZNYq9YxE{#8;xdM zs#4=?^_umRydL(k_}EyyoHMT-;;Js=_EtOj^uc{jw=AfwoM6wsQP2JEb!tyyj4i6@ z>3bmuq>UTn8`J&9#i(F8E70-$EmdPr6K8PG)32Q}m z8Xetr=%=9PxC&1{YzgTMyL@fdy+&4PWsNy(R!B_d@ZyabaYAPpS#YCRX?8JkpfGz1 zxkECaL9wh^p{irJq*=Evo`D=o4L(|+E*2ssw4Bflu|kiz1{(d^Na7d^*I_g`U|+91E|H+zKsvu?&6%$sz-BI9 z$gr9@-gGwys#I-R8E%kOH}Ll1h(sAnfNpc zMp2ulsvaLvOXm!z2bynk*HuaK7l~2O_S@{)wV~%)<*{t9JoM(hfe+(%uVMPe8o+v@ z3n>DdV;hclvY`{^YBi>_UyECZ9z5!K#2jr}D=nTK;mN$xljFxAKVHm)qEo&y*|9(Q z@m}J~@>RuP*4B82|Yzfy&JgXmL zG+{BOoh@pD!G;kZ{N8$PaCR@-+P#Qj_E)>NSIIsH0yOMb424QA3N z1;QGW@%$($nz^*{ya8;9ie38B%toU@%7=!PI0Y<29OdO3E8w|+%YUsO)m`0!1((S_ z$ap2Ur^2&-LgQGWqR*bvmS}pk3F;RcO67qK;n?h01$}?`P?%8PuI(b3DVxtn%gsyj zl;9F`$-w3@h`RYXB7JAirw=)MIYDd*jiX|oZ5do0O|gbWE1hARu4T}H0GR_Jhj5K# zPkqfr&*o$t{}?gfnbfEqNJhoo>xhwV+?mPDJfah~p`vW?K zGLIM5i{xoIh}Dxc+|Ak;-%;Xbjo+V}H9UduEL*kgC?!ohp(QkVKuWpXww~#@KA$x= zjUmgxu^lY7-AWL*W&t0DGvMc%#bdoYnf0iE9WT?5N`5<<96lBaQw$9@8kvPCHj_Ls zO`=gNQSP3kBpWY)bGT_2fiQaBnz=U^>L$N{oS;qQ5%8S>RHsCyx5djj4Bk4{%f(^)XM;IoV- zozsW(c&-Vj51-wcyY|%X%to%>ft)bnbhw6P7(1YacWQgPcyUA%xi3zcEH7!mgFb!f zsN#mjQujqDaT>5irl@`_QH@?ocD$rucAgS+%3S0GBh`6ZgeXsY+tUaJzoi(xH}`C_ z&|kbn|NJ3r*0Yf`sNpUr*SZ|9)qnakDz$F{h?Aib%#b*MQ#0HpU7*l_p!O$xR5PF&(7*K@}N2g3H}p4UpAUe>1A74;y#G5Ehfy-hg`XYkC^P? zbY@85+E%n$>-apU=3cAL^NPf?p6P#|g3CYHHR@k;ct=|D6m-w5I5B*hG~-f6K7g1M zLuxgpX*Inm0e=udumL>4P)yL)e6;X%kpVWR2=u81vY>Uy85r?+FCeq^d>T*^eQS=^V6FM}_2bQj59CFs7y%prOu`6|fp z+cvWb$dK`nX!v+-vwvcAAVH7AQQXlH$G1r_E?jMv0{vW7gghMgIa5gVXPar$Fm4TY z7kxRhaqGY~c7}vHlP9Si86%D@<9-s~HwMkQIZy-FP?ztXil+rkrjKerF9|zLnbEn5fT!lm^P0m&i^VTvRt4n>Q-TH& zZi(R!nsBCDA~yR;vp%aPpENBma;(u-$zF9yb(h8j8Eu~EQljY2v)T6L<$a+3^iUW` zA(ldP?0LnwBbkSYf}KtzIKZdTe^Kb*HYQnq4!rl;-;JZSCZyxj_X&lbth&@_DvXEn-Q5QO(d?b6gXmIz-;^w-Ym#!$CVhvF)9~u+7uP* z0t43QZ8kTyl<{TumHI|GEh-&qRWlkkU7lQ}kQhG1zW2p|$L~e77SQ5+$k<*W4Up*^ z_Wp6zy4u_T`}aGkPBLfXJ78a3nG@fgMha)c&kfT&^5X-(ECi&E14If|>NPgPp$_bZ zg3gmSDqWjBpMcvoF!d9P*U8B`#5w&O0WW@D=h=dT4OI)<%_QV0sx z>WJ&@J@CL&-r-ogw6?2iyBElUC2XU~Zd9*w!EjOLI5`3AYfL z*O0H&pv>kY5sY9M=VK<;U^mGC$KnH>`lN9VHJ&mYzs8xy=%~Sb)NQpd)I#<0DF+nS zkNRa$#WKY-<`tm3v4-|#Fp`S|RS4x0<-S@AzEe=gsbbZNN6p`n3^ge+f?t?J> zC0#G<2jPYL{-UVhry@(D?HJCTANO(EF~4+ToYS}C(n^BGTxg}ahQ4wEC+T#8^TA=c z9){iL!=~+@wABw==`N61E_faUp&mEO+PVZt2|w$4Ns-K4gIX9_aI6ZfekpZEAPOj4 zBcdsj`hsT=tr^($@T&0Zn@h=VypP(YXznM8HBk;enBuuihCjG@C51%${r zowk0poC_z!X7HCi&b$q3_zs{!+`Z{EeJ^GI{&}G#9x^ zgdk4UWg8jIDzt9fm_(!y&>XX}Hqp>wUNkyY@TezJNGc6pF&f+4np88sC`qA+>|n3l zrZEw0TPhdNi$XqNMPd)a(TA$xwZ}_D6$Yh}*S2phc%oUNF?|nfUe%0RKJ8GhO?E~Y z2%!uEG3|TAI2F`5i5}wi^3oLQ+lQvY*p(rrF5NCVOfHD>@K1x+M#*1OyH>kL7YRl2 z3tt7AjJl%&>yB8dz`zT6s3=k)vEQg?fofrmFVU>RQaYKb%5nxpL#pgIL$Y1cVi9bj zS+{!Dj%bmiB39pwu3?4zCxv{fT^s(uVX9aiL>7x<1z=j7yeV>yT96Z7sUkU3pAHZ= zR^FMwAWjoE)c0*(Td$cYTA>0e1`dc_8CTQG9P2XqwL9Nsnv}G`*9z9JB0!cnIW-Hu zEQOL2v!PM}RLr`1);UoUV)U$_QFfsJ>PVI=*$~USkm^>2WK+52Z>}3p5KGNHLU=CA zS$(s+?!V3dPUk`)PkSU6!1n%$-13}QLq1M$T8l^shnQys-H}#HHO{=08It@%)@wuR zOM$@gNMD(dy>_GIy<(TQo3`&Z$*j$eshmHfXr2;;KF)cDL%7ksQR+7uKx&n%ecppC z=*JM2&J>Hg7WDy_*m7N~#hhV_sV3q0N`wB9e)O1{OlVU>*O^*pP>%| zBgW`LbJY|Am^C_C8A-#|d@X7bcD$cn@ndV@!guef|qV3Y+}p3f(!3nBK! zOMDpt`R_#)_OcupS#ov2&~1uQc|>F0OReN5a!?z_tvfEl>0{k?M8XB#$U~j1vdi_c z8=`^kd`9fbX`1J*jF;4bN4&G@%OT8y9T?qyi85>sRdSi=^ahe{FIXicmk*g6(qD&? zQJ3tOFdx$eMk_)z&nuWGhU1nISC~S$dHkQA{E;D5p^5I zoW_!@T#0G=NQ4Hm8p$q+@k7{lVhSE0PFd9X0S~_77995yfr#aL(I*s%w9dVW{NSz_UktL2#XuS5RagfE3nh@768TE^JFIVlQzozBbGDSsP^awB{MZ=k z(A+`k?`liMA)|4{g;POeTx8}<@Hk>I|3=oE1zpdzA%iHAv_Qu3W*#7hqXO`+b8n%qZq$ZSGnL)j5@qf$dDZ4I&|Z;tpub& zCfQS~1xmwC@ci|L9&AIRUMc*z#PdyAxpHZ{S@s)YPIPqtzaW-m)Zs4loZ0`bLbA~cZtG@tm|cK zH>2NTaJeTBUypu5|NfieGBM^T4Wf3)Ve}Xen|K2rH6yQDpG!8Lu=W=6rVcOwj)x5i zp_d6Ne^yUHXyvM|H}I9(>6-TT5(B+?JqCvi07x?W?God!Z^vl?cmC!Q!^zzP>g3_h z!z(XM|L>O(P3m7IvLsZ%d9EQHMxB_jzGY4r6*&nA;HwI{*z|_?y9`Nu_Wmp}?fwx& zG%_MyP2~i8IXnyx1}1kTQ?kG|i{xRbrdW5$E+k<`mp$+=5&5v7qac`wbV56_03m#! zq94Ymu?TM`Yk&<`hj&cM*|e9q^1>4q2HFEXrr?X5k70%pcCbeg2CySW5|}Q05e9c@ zhx49r!^caGmFYAig)5@XUk(33>n`sQ(Ku_)h-?%(4AY!q_(lj-rat z5wN)dx;(1`L61XW1egY}1eSO3vL*=_I4TDAB+M9w--iajweaahL(BfnHUDkBz3Trv z3(NeAh56*A!T)aIf4pyYC{1=V_^3rI( zb9_dC3|Fat2WlXXfiK*xhSTe$!9y&Y;L^@DaA-XP*m&p?6tw#iKHYLLAQ79C5=++a%Y9#zYU3|3gayt%KVrHVmGpq~om!v_kL@uL_ zPlrYw46CIbODhK*Z_w++lZhSRhYu;hr)-1pPLoDB66Z2J{{1K%S!V=JhNlHy`Ai8` zbgG2U%jLtd4J$$73?^{C&3ZUScO=}O-v{3ARu2z5rvRq~d;_(=Pk}#G(SolyM1T?$ zGvMNu3GlF@I5;YD3!H*H4$eA21SaTd0;LF)!`(|f;YefeC1wJ$;o&T)@FdDNaL$Yf z_~=#>e76H1d{FZmq_v(4S3$Og3qm$P0Xai(z0CsnJ<=gK4gVxu8C(D-+(rVsUl)N8 zOe)~Wt$A<-rFoEEQUkmLJrQnY+7GvQ5Cli3dP(+HMilpl=zkaj|0|An zBX0XcEAU6#&sy<+ff95j_>We?|EZPe7f~Ts;u}$-8?7X_qKH3`fj0yEeLSRpXuYx6 zALAsueg6~FpZk#i(dWj(H@$x9M{(QsV-A2n+E9PeMd|S)cTFjUdJF*(0JyDxzhLPB zynpF%+_xEZJddp0+;3)7UYhft9qvyHnU5EE>+%MXzH7H$b6fwO0rlTxUYL7(SUKN1 zlDzc)Hh;pw&cnmO>Xt&N#40)Jjr!aA_Y~Lx4F6yRcPlqLsDqsk^wzl{c9!zeH2-1^ zK3-lCp}Ty1qN0Mlg8aN9q5?vEy!?Wqydna;e4-+vqC$egBK&-Uw+t(I?TY@0U$>L> z$I0X8^h2uKhF{7*yvuK05jjM^_KDL23V3hDnldGXJxe?GeE|DxLPW>)@B zRsWH5^BajjEcx^4j`kN4G&ev0n#6z4%zp;^`B?TV*i*WH1LONM*w5#uU%`mz{|$`) z&tN}q(!YYW{41EipTU0KS$_q~WBfN}H?IAwEC0Ox{R$?^^lxB7e+K(GwB3yBr?T{K zV7L3{pW%LvQ@_G#vivo-U!vom;eL)e)IV_p0LK3Y?!Te#&rm-H2lii}^loCtU%BmX z%>Ep*G*nPgIsen228rvxNZ9}R_WRTD_MgU0IlqlaH|3_^?GXW>`JXNbKi;9-MpFRx L&F5DBH`)IOAx`Ohe zL|GzBA^Xzu%~a3xptrvF{g&JBGQ;)1|L0uioa@}@zK)>|1tkbXOAGqxItZc%ok!Ut z-DK=h-YB#S(oqtP!l5L*T-^{VG$5<**F|!BBogO|wRc7$CBA)1IQmOCVNqVrNOw=P zkDG%8@-Kip%EQUc0VxGWx}n_d&?qnlBoY=Z>|=)^lQ7hwp{4j^ZBmn{2Wc1*^$3vl>DZKq@%y&PBn=0KlJ=d4GFIRgbM7hh5T60fBJF7T+cHP za-d<|d}CVMjgtC$wXvDio$3c4jtp5xh|}oCzPuowS~F9DTy}7?zpkCgoZ$4jqdw}H zfa7f!kfQ+n808lU{Ego2XNn+5xuzbr&1J-vVgD${>dP}31?eI=BAV@6nSau`3G47g zvPjNtuubTFZP0?ty-PQ#P*@m~_s6zqVH{cWI#b^;d4aBNy|-Ins_)4anG0I{u)VQ- z_l-x~88ZO|iPAV>RLiY_(E}lGwzEXK6*?oc$1%7(rlQS$<1H;c@jeC(TAEr4be2xH zQq;pLijXe>fgzh8SC;ByG(uhXA6JE4!V4tA#^r6Oc(z}>u=l(}7k13MOwIX=U^xTq z%*(cG9%Mcmu9nW3)-OeyMf=2yo@mbTSO%VUT@hWS=Svb7e1c$pGjt_15`FB2c86?w z@D<1Wu4BDx&%<5S#O~8+_pp^vT#0OaZcwQ$`@9G-XELXy;tC&=ElJFrUDPp&lCn@t zX(%;*`4ZarEJ8d^#KW%X3GU`{u^EMfLZipmfKx5&el4pHwjl!_ni?J*Mv3G>k>lqb z+5F#4=QA9lt_&wk7+9rTX-;G>;y*bmT&tEy!5mUNJytAUM8Q#)be!iPT`-HwQ}tvv zg~9+!FJt21>qbxwj>*KeEJyxPzk z_F_1tT-v8sQUU#Fsq%MhIcpCM*Ia|X9?N-%WjCgJkUu1u5B^|0V}9?QRkcc~Ritt0 z9X-+R&*5g8y(PImudx#ima?2yF&PWfUHmm;;WZ~59)xSiDqVMzV{2Ck>E00YxMwiL zQ88O{c7Vc{TVsgl%9or3g9 z2>X7Q8;n6*|7W>D0_|Yu}EidT4wIii1$w_*z>5j}tRnmU%w#Zbp3EQ1Fmouf`b#&Kage`Zzuq zW>Mt#sKU@m7ta*3l%3XUo$^`nvsrm+>U{U4&Sy+}$J(px`yUM(D-0+CwkU5WgxUHS^#2nHNHPNy>qY*@mcnmZn^19PF+MmUJWI$CN6U(|! z-EtuPR>-;^GAD-qLCK_5dJuVXDWvO@{9(UEfh4eWLBNxhj!P}B`Zvh{cYs+|`jy># zzYa&Kz5=O3dNu5e)(v@Jx((og7)>r~`?l~8@*;COQ?XOc_}25; zZ}S8-#^p}S*5btKrNiRlUnQ5K%!l63osLi-lg{k0RYXSBb!a^ueqHQzoaJR8e{(4C*k%Sp~-mSR9%y)Rq1&f^FM35@I5FlGbht*jpRZn|Cm=$X?%M; ze?7qzTq>%B@D+@6=LR!%I2U3J?^|r$OzvWLPKrZKc`RvP7~ubiLgutp|h=2%CgBrW|l02&gGu=z_eZupOL78X6q=kb+PzshU+C|Cgvfj9o94V zK}-o$(Hr2>h_?QBMw?9PlN$b+D#I?yKl7+~Tj^@=;$kuP!$- z4@J+<_KWo)T?+1>YfT6gfX|AR;Y*O&8mk&5C2%X#&xy0M3D~rdirIv&(k-8yEwQ^5 z^7BUdQV%O<6HeB!HKk@`T>0Q2?>MAuaSq27wDn@?swquqSGT|98|X^w0wvJt4183gBFaEIQ*9(^q(K;ViO=!$W##@6~pJ-vDU zlb6`}7f#(-@3RZ<_HSsL28$e2u|)ZLcghExWQq_`tm)JHqx^|()OmW(lI@;ti>h1Yw7L$_mOw5TP{pw1lkFn zNjDm%8t~JrTu6!$k4l@9ziK8~1e*(&+MY~)Uu9rxLtn)ylTI7X-?}(j`C=9Eg|l%m?PBOxENPd;>giI7kBprG#Z#*-Ab#Y%tLXWCH6K$EasIWbEg%i1`X<=YOn z%c=ydT19J^+bYKEOzIZ2(x~&zINi}!Z|wGs=Lo23Zx-v;d+XP%otfF%nZ$))${!6W z>=*I#>N-Wvq76gU{Mk`I=g zxUk3Q#1=Hdg92smEp$j|-EHFzkpDUo$5i-HL>lqT>@d4!pfq^-C3v5%g8A*(KINyG zC9X}QT>X^2ba#2~?R$iS(gbY=x78;l1Q9C1x2Fi9ua?r*N;2x*lIrU9f>#f1ksthO zG#N+fU6t^NNv$aTrF0GN!ocO)jL-G%9}aJzl`vU^o25KQ74FpOm_4_c^f1)V@2R2L zV7wv5q`Djg6Q;)w^>b5rV-1}rq90ySExRJtwEyFZd7xnZMY$XJ3fil6b!XpCv~UG; z6$Eo7ZV$W^sKVOk!tZ6>i)c34Fkq@<+Gp1 zPrlc7a(H+0X#!Wd34>%=p$nq)k(K5&CNE{($7HbI*L?Y-5UVn`W=ww?^VCC?0wyMGnEp>^3u*j9_^7$EU?@YPz zt10!C9UT_uJiYklW?itQM|dh3-{*+u>?bIX!Ns)3~WB`Jb zmS&RJ9;MzZu+RUKt?KP+Y!PR{pvR09Wc8}t)VO{~sLvOl7ZM)%+-nPu&S|ICVZQVj z7w1pAt#&p%S4XGhV)LWJQ&%6wdpry5t#*vH?2V~_ppHCmeC@rWBQzOWWv0a4anL(OLb1Z&^Je=i)tsDb@Z?pb_|xS8 zY0KJIqRlczT9x8Abu8an!jcuT<%03j#W<-OZB{cTtxIh4y@E0?PVSdAyhM-LX59Mv z;ziO$tD;MR&KD8R;f~H>W{al}i--GKFzDEd;2Nc-BQkE$-%yn=yO>-FHdAFR3%g&r z|5N1=kIVZindBZj4v=%MG?@)hlTBX=p=|<1oQi^|9R**4Q2qfSe}!W8ZPXv2deC^M z^6*;Yx!Qeek3$qdG-C&b_k*Iy9y2L`{TiVP`-N;7^qWK!_EFyqQ$%sihjVg1r<56O z5>O)Vwk>+HxV9otPqkoM^v9ye=JuP>1&&$hCyk9SZ)UV`$m0*SFU#9j->;Eka4$`+ zg@Hq@<6iSsGUesnf1woMY%?V&{Jun$qWLrn!lF75=ukd2mNMJpN9*=J&j8RZqr7XQOtzKjyi^ygC@^e9H5{z}4RE+&)J1yee2;Ha7>s5C@$tI#^a- zoKow<3VRuS?x+_X;@KQksStWzAaX%IcT12#2O6l!wRJXa8y1<7(L)=?t#&ffQ!dAq zQi4B?%_)S>`)jV3gY?i)d=5+AJoV|V%hWNi)KEExlr=~7)Z8!G`shyS8*RLET%ij) z%Mon(W1$xGqSdw*S9{xY`^?ev{j9Lp;UQYK@`(N?IoYkl`@lTz>IOA-7GXwCc}H!} z<~2tO56`GNUEnaY!y&lNx7idJEv6>d)wi&r%hgc1Y-f)iF;C=68r&ASFc|e(&`TBp z9H6S_^R9^1KZmdDSD!bKS{ID$f7N^M2}dC0yc@GkwXBv+VT(`uP2^j-Hk<60NVG>v z%ZKE;VAd>TJS5B&@=$H8@5!_!dB5EXKluzS4#kZycbLb3D~=<$bs!tiE&9l8$36JOZc!-rhoSg; z`MWNdPrYWT3vC$|>TEjs<;+yM(S{LCiREP!%OH()L{1oc%v>Sl0Rz#9j`8EiRrrJ) zB(%mJwfF2Je_*@f9CviSq(fPMtVw}W%c#V`{KbsA&@-Y zLyQ&?Fn(nE9Jwk-1ONVDgDKdp#UQaByclE}?d{5P*~0^`aHT1<;*^e=9W8z-AM%Jf zIHQpkVeo*CE&wr=mB08R*|}Gg!N|BX7oUKhw+Hu%;^f0FaU8c4s2SsnvDO(&FwaS$ zE?T^AE`8CQB4d+&jV0A1$lRq!8ZA)0|V$K75$%sQ?v*T;V{`EAO{9ZJOV&9yw?pB-hre}dY z-jwsr(^0v7*)78cR9&J1USGeGkzT;=iZQiJGJrrSV&7lD{`H{|Cy3?e7cdWR9LfXd zEeS!WFz$Vks?qE^Uk)e8t0dM_O0xhR zrR@Z$;TxcyqnglN5!KMA?o?lZfdGlbQh*+$>BX5Nzign8MislhiO6v^fLFNbFGN1(Czr7a;Is?VAMm`Xy!Spe%q&m#iUl zW*WA&eh-FPi6k)7TM!b3T7lvkWdb}lo)B``n!waT1CUI)wpm@fNIc^r4bnN+uS`k% zYg0lIDwN+$i5_ABj;b)jMLDYotJgXR@aY_a=Dt#b$@IsDlpqGUe;X9;+Y&?|$vtwF zIb4`1lQe!O_rR~^N+DDjcDe8>GbNx?)e5UV6c4;UQ33FpqyQ1<8bHOf5PZfro5{1Tw-6*yF(Hde>g5hOe&zVJ*X^Fiws92Gc4wLZxyrtO?2kl(>iy zu4i0VzRIwvRDLj=aQ(Pd!`mgR2EK{HhSBkl4b8Ool#_?k#3@Yv=)Kp{Q` zTx;tEC{20*jswQ<@u|abO^-5QL_G)CXITbQ;H^GWaskDIUBD5k zK0q6u3$V;lz`Z}+hLPEn0+e;P01d4Xm>IqrXrR3gIGk+(UZtpPf-A`5D8TD0htv z1R{;!`C0haewIe42<-Cnug*n-M~mV8X}kb4r5M03jD$ndgW8$caz5N{08~8-{cS~;{TT4{_-dPNA)GD zZ}MnpQ4{X5|4aGu|Nrv2KoY-P|M%vE=zm~=i{!%g-`f^QADoMuHyDW}9{ahtp^-QT z4{uK_(pQSurVtjkpdlK`n}|d?Ibg`X-vUOsd7^MznQs+aOZCdl7Cw|c*h(!X>NyEcFr9k&3*F{nCKsp zum0iW!F_8@t-1$2#F(lN3V;73u#V?El*SR&FjhoSOqlg-GRzL)1j^r15JNj)8>ss_=b@ z-VRt7l$%Q+iWJvK7c@eJ|KHIE3YAlklRY6PPaG%Ypb%N86a*qCB@307mXSSiLPn1G zOA6c88mkhq!RzyL-Idv zI_{Oav$8hcCAFGZ8vb8N{imJrUH(6&+gXls?vgLd{4eDHUCr<1*;!uh=DEo78y=!- zce`?D`Ldg5<99q#dwF)28@qW%*?yys^j@Bw<-l$pfA-(-$n53W`FP*WbCejdKl%AD zJIn6n+4-2>&6D~&9=U&?hlArc`pEC)*_j);cEwNXZ+J+L>b*=mv*2zf>BB#1v@5gj zW!jm{h)V3Z@IN#CokRDs>`XYvc5r||>_qkUcvCYVjYp9-%n?Hwoj9Z-0Fj o5w|+&s{OHU;(erw17yof`V9h+N()0BDpGO=F%o}B`0v>N0^$7#`2YX_ literal 0 HcmV?d00001 diff --git a/bioptim/examples/getting_started/_l4c_generated/l4casadi_f.cpp b/bioptim/examples/getting_started/_l4c_generated/l4casadi_f.cpp new file mode 100644 index 000000000..730d5df7f --- /dev/null +++ b/bioptim/examples/getting_started/_l4c_generated/l4casadi_f.cpp @@ -0,0 +1,111 @@ +#include + +L4CasADi l4casadi("/home/pariterre/Programmation/Technopole/bioptim/bioptim/examples/getting_started/_l4c_generated", "l4casadi_f", 1, 6, 1, 2, "cpu", true, true, true, false, true, false); + +#ifdef __cplusplus +extern "C" { +#endif + +#include + +#ifndef casadi_real +#define casadi_real double +#endif + +#ifndef casadi_int +#define casadi_int long long int +#endif + +/* Symbol visibility in DLLs */ +#ifndef CASADI_SYMBOL_EXPORT + #if defined(_WIN32) || defined(__WIN32__) || defined(__CYGWIN__) + #if defined(STATIC_LINKED) + #define CASADI_SYMBOL_EXPORT + #else + #define CASADI_SYMBOL_EXPORT __declspec(dllexport) + #endif + #elif defined(__GNUC__) && defined(GCC_HASCLASSVISIBILITY) + #define CASADI_SYMBOL_EXPORT __attribute__ ((visibility ("default"))) + #else + #define CASADI_SYMBOL_EXPORT + #endif +#endif + +// Function l4casadi_f + +static const casadi_int l4casadi_f_s_in0[3] = { 1, 6, 1}; +static const casadi_int l4casadi_f_s_out0[3] = { 1, 2, 1}; + +// Only single input, single output is supported at the moment +CASADI_SYMBOL_EXPORT casadi_int l4casadi_f_n_in(void) { return 1;} +CASADI_SYMBOL_EXPORT casadi_int l4casadi_f_n_out(void) { return 1;} + +CASADI_SYMBOL_EXPORT const casadi_int* l4casadi_f_sparsity_in(casadi_int i) { + switch (i) { + case 0: return l4casadi_f_s_in0; + default: return 0; + } +} + +CASADI_SYMBOL_EXPORT const casadi_int* l4casadi_f_sparsity_out(casadi_int i) { + switch (i) { + case 0: return l4casadi_f_s_out0; + default: return 0; + } +} + +CASADI_SYMBOL_EXPORT int l4casadi_f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, int mem){ + l4casadi.forward(arg[0], res[0]); + return 0; +} + +// Jacobian l4casadi_f + +CASADI_SYMBOL_EXPORT casadi_int jac_l4casadi_f_n_in(void) { return 2;} +CASADI_SYMBOL_EXPORT casadi_int jac_l4casadi_f_n_out(void) { return 1;} + +CASADI_SYMBOL_EXPORT int jac_l4casadi_f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, int mem){ + l4casadi.jac(arg[0], res[0]); + return 0; +} + + + +// adj1 l4casadi_f + +CASADI_SYMBOL_EXPORT casadi_int adj1_l4casadi_f_n_in(void) { return 3;} +CASADI_SYMBOL_EXPORT casadi_int adj1_l4casadi_f_n_out(void) { return 1;} + +CASADI_SYMBOL_EXPORT int adj1_l4casadi_f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, int mem){ + // adj1 [i0, out_o0, adj_o0] -> [out_adj_i0] + l4casadi.adj1(arg[0], arg[2], res[0]); + return 0; +} + + +// jac_adj1 l4casadi_f + +CASADI_SYMBOL_EXPORT casadi_int jac_adj1_l4casadi_f_n_in(void) { return 4;} +CASADI_SYMBOL_EXPORT casadi_int jac_adj1_l4casadi_f_n_out(void) { return 3;} + +CASADI_SYMBOL_EXPORT int jac_adj1_l4casadi_f(const casadi_real** arg, casadi_real** res, casadi_int* iw, casadi_real* w, int mem){ + // jac_adj1 [i0, out_o0, adj_o0, out_adj_i0] -> [jac_adj_i0_i0, jac_adj_i0_out_o0, jac_adj_i0_adj_o0] + if (res[1] != NULL) { + l4casadi.invalid_argument("jac_adj_i0_out_o0 is not provided by L4CasADi. If you need this feature, please contact the L4CasADi developer."); + } + if (res[2] != NULL) { + l4casadi.invalid_argument("jac_adj_i0_adj_o0 is not provided by L4CasADi. If you need this feature, please contact the L4CasADi developer."); + } + if (res[0] == NULL) { + l4casadi.invalid_argument("L4CasADi can only provide jac_adj_i0_i0 for jac_adj1_l4casadi_f function. If you need this feature, please contact the L4CasADi developer."); + } + l4casadi.jac_adj1(arg[0], arg[2], res[0]); + return 0; +} + + + + +#ifdef __cplusplus +} /* extern "C" */ +#endif \ No newline at end of file diff --git a/bioptim/examples/getting_started/_l4c_generated/l4casadi_f.pt b/bioptim/examples/getting_started/_l4c_generated/l4casadi_f.pt new file mode 100644 index 0000000000000000000000000000000000000000..8c503c1ff319972d9f106b1ac61a79636e677e91 GIT binary patch literal 6342 zcmcgx2{_bS8y`cHeUBlDQbsWbV~d$7H&Yt7L`Ak~ni*LJ%^0pFw@|W`>=lYc*@cuc z|08rGvXn~FqD9N41to6FH$!U1&He8EzVG?|&vTw*`Tc*t_nh}V?|J7PJCd*n3??cH z`_rQUlYud?R7wDa=0*0X%ulLOeZ276}X)#gtjBi75%lrP-rdJ4&Jci z+NH>~>(7R*h9UUvqWaV5`qK(z{Yew_J%aW5P2m3jtBunppm&B4t>D~+97rYmL}3&! zKN!;)_Pi8(z}bQlb>>;1%#e$Y=FMU~tzK*TZu}?1K$A1Ek|OUyUzv}6)+k80WxOb1 zm`cpPz05~NH`!vZXsz8UbF7B9Rp;W~g%@6)xfW2Gy6WKYx?Ofj=EfD4Nhb>AU&P9m zUkn>REYl7b+NF~~g)#@caGhFNy zIG((CRBYX~7k8DjYxWHE!7^4@pBa;h+W6@{d{<*PNzWo`U;M}W$s-$EYV+dVMu%(e zygFuju{1+Elp390Up|0sB_7xnR5H-`YpCVc-tbTycF zUS-IO9b2`)ThXc!pF7X^mL1i?8>eOXO5s}WZyw!Uz-|q-DhO>+Xn5;t_FB{-GlVGib#WybK7VO5wFqx7 zmJG8zcIBXM4kx}vVZV~h(nLM9+EKzG*S?Fd-~8IRSZD8hr0-Sa_*$RVN32}7yIy{N zWxAVmoc9Y`EtTx7d6MPHn#)%Y=RJ)Zi2HEJQ)ykH?P)vR{hJ(G_ZMa4r*!3@^ygfz zuC?zTv~Il<>HXW2oWq+P#pm6C+)YTK&ZLy5z2-RU+Z#^CrN!D@dZY7)&(me&H$w5* zw)5hotyS!cLb4Tn%35BITe&MpfqOV3MLzBGJudK-Ir z>sg%lihHgpuBjjVRF|CIk9Z(`Q!fhl{`m*yrq?wWkBb@D`}7QZiIr_!Z=TR@b5Oyj z#N8nXLr>dUtX8r3P$@S4MzK6N$J9-5F28Hgg1(|9K|DTG&wg2^qf+^kA?@xe!0H)S zd!kgCRN8atbbCXx>nYvd!kw(1)Of9sV@GafW!YjzRPCf!EU4VovA1WgOBJm}>cRoS zfZ^xE+e#cO-=&(Jy?t1r*DcTXhMuPBjwK#<`*PYWj&CeUE9ogMpE$dEf7XcnMNrq2 zHDW)8FHze*{0Lh<(1cAqRc87qh{nuV953-3Dg{})uc>K)yu!H!3AMX2v>xj%`3)}} zrk*}XG%;o^FO0p>-7I~hT-D~55!7OC{5I?qvFua<`rMr%QF-*5L$V1>v@uJfFxE7G zX&2k2DbpdP!Z*&Z(55<7^@z>SoH4!jLYqFaw|skHVc3JBwzBlDO6RNWFsaii_U}rQ zp4)p}>hS=y^r|P}oxb3jO+xAtdBugPPn{UVgH7r^*)bh$T0v_(inQy^D?j?3yYk4j zc^Sn~!o;?GS?dH_rB#}~VZ=!dJaU^`v%~V*YXwA-D<>y-H??TRlL6!0 zC*#CE-<==5jk0Q;{R7_ZMya|~#jSfUVf_2g&>y28Cf*~p({^OO(!Ub5v#znI-#ak4 z@}6B|>kTx$dLZ}Yku-aVk{)T0U$;(au&cfB&XbUPZCV}P%VqY{n6;alWBYA2Pwi`I zjH>I(Sm$hZ>RI5htWEKmigbdbv~r*^8Fe%iopp1={`ln{1&0#%onW1}HL+E`*}e5$ z!WL=AKcw`@qCU!T&NdxhAr9A~_GthE!L=jf5wb>#HvXO~BooOSBHPVB@M zEO~RX#0yuN`6Q|ERJ(_9Kij&aASLvE+e^9w#_D?V$s-aEg=9<`Zpuvj;TJ96iEw`t zu@K%Vy`WRfO=G3$XuYbF1gCTErwVDrHrZQ^2xKG$8L<%=b)Z7bq|q#LBr9`l zx6LKKt|(eK+!rYRW%Pvmbt`S+xve+JA+JoucZ5GTi*Q5QxO==fnen&)RX6aHzZ_~C zJGDXKSK5RDn>q@=IAVP)EVkkHr{})#u$Vx7pX_I44#8IvC*a6i+eQ_NCyYAo71&AT z7)JOQyc{zq%G{-9y1xIM5PN^UdxFxggqH@#EY|nG4;o!-5sgnMt0a}b?OU$s&RF9p zHmH3|XIHgw_RIQ;gSP4S#w|KNtdwL-qAYhsu6*A?Da&yKH`LWVaL-G5k{sG~M$)sL3^3HhlJ{VC(w=e&wFD<8w7 zDG3O8Ld09x6>FVc(KlUt6cgGvoBMtr+I`y0SP_<$q?+y2L z`C%U%ph)J)I0>5uzY3^nR<;~l(i(m_fNXTMK<;wxm&to?5kdBnGXe&qDNNmaKd#i2 zX5M>#0c?sNJ3t>pwEB1VSIzoQxGWrrki*=|alRx0HTUa6YAQ;2SePw-{6z?8-Pw*C z2pPcj1l+|dc6M5ptNB4?9UUM_Je89=UmJRHk_aK(B+ZfIE1{8z%iwiPE)+Qu37J|Y z)rlA>fVSdR$iVIas6$pm9p(EuH?4Qn6)-U%>3&;X1p794p!*B7HO31HigAKQ7ArxU z!2<~FtOMaec~I7=3<%oxiX#`P1U4Oxz#BAWg3P2goF)v8Z+(p+ENLBweCvjBE@8*O z99tpI?ibJT6JwQl%IH2wUcw1V(ry8TH8_Y6mjL}5?E=X)&jq|{;{AB_`8V7elQ(#E z=4`tvDOgQR@2UaOO7vS*KKMPrp{x4VINh(d!1|G zlIvL@tlJL`UF`xQn>qnf)}HWaU_N23UkSKpoeQM4mf$THxdMY1RY0sU7K9syf;v_e zi2k5T$k_cEUw1VP>>=5MM~oPJ(%PfIl9mLbPbLC{$a$cuk_b?3iUj$~HTX26A`nmz z2!yX+#ogJR4fboNgOvGAK=)`2=z3WLhB$JBw-xnx`)9`iNyHOi$FHQgp z<#wQMcoVE6oB#^H3ljo9p2Q2e76Xym9I(db9)3$oHQ*Vu>WbsoN-x7HED6J}+u6=JYAXu8yyAf4`I)Ji zWO9gGt@6fS`C0eJel{Xnq5kCOX@Ae0!ygtBZXeME@*>N@VicKxIT}t7O^zZ&nIsZO zLq`c3v_L{8R&_G&p|q0S|3&;^e-wXTl+*}i8gW3qq8L0xxDg6Jvk1W{41!}Kn~)l&Xm6r`90Ut|cpD9JWVV1IH#|@~k6Ovwzv{i( zH|=<5uy3)h*mvMSBJH11@H4FZ~T?tus`;j3DFAkxBWKlM-2;& zHR{1Uy}6=w-985Yg}&+kM;{3@{Nehingi&6xC49ODC$((fE>v7Vg!hjS=O8aJ52pioC$a0o;MhVDeNT7ovAwm znIa@h_D1F$5)4H`ki^+2k`M{Oj-^ul84Nm=?d9(mAijY`qqFF=)m~J#I1`~`&$egI z_2lY#3VBZP88AgXMYu0-WQy^4CRy0Ao>Ne90jMEUf(xZ_p^}r(uae^5;8Fr`Bc}9M zW|==T%L+iTOgS!8k;^PU#SEWZ%y(u50k|R0G|9#m+O7furT#;nGkY~!eF+_&uDx`H=mKSoH{a{e$v^ZHBrL3CCC*DXXV z?Y|pF1{iVP(Q%3md&=`VZ_KW#2(#eV+#bTF`}6!yv@-p>dei52dS-aB*$g_*gnhWj zD%?_y+nxF=6EuIc>;O8;i^A{_z)0>o+zw!I&ChOifLUxRZe?*3ihn2xDe~ z!I~Ih%&`~~tdXgqnE}tj70(r~f6xCyaIT@>I)D7VH*T>xOTo_~ z0h5@8gX=~?xAB)6f*irKaTw0z;4cdVIkrj7Vu#U84*n}!kmJK_9N3v0{MWA_M~?I? zb{Nm(;J-ozIgm25aG3mkj>EHYn9k(j&w)sR_(013*bd(HIg^Dys|m7@)nnthnHZv-%W78o#5i9u>S%!?y5%s literal 0 HcmV?d00001 diff --git a/bioptim/examples/getting_started/custom_non_casadi_dynamics.py b/bioptim/examples/getting_started/custom_non_casadi_dynamics.py index 818b0cdbd..a32513ede 100644 --- a/bioptim/examples/getting_started/custom_non_casadi_dynamics.py +++ b/bioptim/examples/getting_started/custom_non_casadi_dynamics.py @@ -1,442 +1,191 @@ """ TODO: Explain what is this example about +TODO: All the documentation This example is similar to the getting_started/pendulum.py example, but the dynamics are computed using a non-casadi based model. This is useful when the dynamics are computed using a different library (e.g. TensorFlow, PyTorch, etc.) """ import biorbd -from bioptim import ( - OptimalControlProgram, - DynamicsFcn, - Objective, - ObjectiveFcn, - BoundsList, - OdeSolver, - OdeSolverBase, - PhaseDynamics, - ControlType, - InitialGuessList, - Dynamics, - CasadiFunctionInterface, - BiorbdModel, - PenaltyController, - Node, -) -from casadi import Function, jacobian, MX, DM - +from bioptim import OptimalControlProgram, DynamicsFcn, BoundsList, Dynamics +from bioptim.models.torch.torch_model import TorchModel import numpy as np - - -class CasadiFunctionInterfaceTest(CasadiFunctionInterface): - """ - This example implements a somewhat simple 5x1 function, with x and y inputs (x => 3x1; y => 4x1) of the form - f(x, y) = np.array( - [ - x[0] * y[1] + y[0] * y[0], - x[1] * x[1] + 2 * y[1], - x[0] * x[1] * x[2], - x[2] * x[1] + 2 * y[3] * y[2], - y[0] * y[1] * y[2] * y[3], - ] +import torch + + +class NeuralNetworkModel(torch.nn.Module): + def __init__( + self, + layer_node_count: tuple[int], + dropout_probability: float, + use_batch_norm: bool, + ): + super(NeuralNetworkModel, self).__init__() + activations = torch.nn.GELU() + + # Initialize the layers of the neural network + self._size_in = layer_node_count[0] + self._size_out = layer_node_count[-1] + first_and_hidden_layers_node_count = layer_node_count[:-1] + layers = torch.nn.ModuleList() + for i in range(len(first_and_hidden_layers_node_count) - 1): + layers.append( + torch.nn.Linear(first_and_hidden_layers_node_count[i], first_and_hidden_layers_node_count[i + 1]) + ) + if use_batch_norm: + torch.nn.BatchNorm1d(first_and_hidden_layers_node_count[i + 1]) + layers.append(activations) + layers.append(torch.nn.Dropout(dropout_probability)) + layers.append(torch.nn.Linear(first_and_hidden_layers_node_count[-1], layer_node_count[-1])) + + self._forward_model = torch.nn.Sequential(*layers) + self._forward_model.to(self.get_torch_device()) + + self._optimizer = torch.optim.Adam(self.parameters(), lr=1e-3) + self._loss_function = torch.nn.HuberLoss() + + # Put the model in evaluation mode + self.eval() + + @property + def size_in(self) -> int: + return self._size_in + + @property + def size_out(self) -> int: + return self._size_out + + def forward(self, x: torch.Tensor) -> torch.Tensor: + output = torch.Tensor(x.shape[0], self._forward_model[-1].out_features) + for i, data in enumerate(x): + output[i, :] = self._forward_model(data) + return output.to(self.get_torch_device()) + + @staticmethod + def get_torch_device() -> torch.device: + return torch.device("cuda" if torch.cuda.is_available() else "cpu") + + def train_me(self, training_data: list[torch.Tensor], validation_data: list[torch.Tensor]): + # More details about scheduler in documentation + scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( + self._optimizer, mode="min", factor=0.1, patience=20, min_lr=1e-8 ) - It implements the equation (5x1) and the jacobians for the inputs x (5x3) and y (5x4). - """ - - def __init__(self, model, opts={}): - super(CasadiFunctionInterfaceTest, self).__init__("CasadiFunctionInterfaceTest", opts) - - def inputs_len(self) -> list[int]: - return [2, 2] + max_epochs = 10 + for _ in range(max_epochs): + self._perform_epoch_training(targets=training_data) + validation_loss = self._perform_epoch_training(targets=validation_data, only_compute=True) + print(f"Validation loss: {validation_loss}") + scheduler.step(validation_loss) # Adjust/reduce learning rate - def outputs_len(self) -> list[int]: - return [2] + def _perform_epoch_training( + self, + targets: list[torch.Tensor], + only_compute: bool = False, + ) -> tuple[float, float]: - def function(self, *args): - x, y = args - x = np.array(x)[:, 0] - y = np.array(y)[:, 0] - return np.array([x[0] * y[1] + x[0] * y[0] * y[0], x[1] * x[1] + 2 * y[1]]) + # Perform the predictions + if only_compute: + with torch.no_grad(): + all_predictions = self(targets[0]) + all_targets = targets[1] - def jacobians(self, *args): - x, y = args - x = np.array(x)[:, 0] - y = np.array(y)[:, 0] - jacobian_x = np.array([[y[1] + y[0] * y[0], 0, 0], [0, 2 * x[1], 0]]) - jacobian_y = np.array([[x[0] * 2 * y[0], x[0], 0, 0], [0, 2, 0, 0]]) - return [jacobian_x, jacobian_y] + else: + # Put the model in training mode + self.train() + # If it is training, we are updating the model with each prediction, we therefore need to do it in a loop + all_predictions = torch.tensor([]).to(self.get_torch_device()) + all_targets = torch.tensor([]).to(self.get_torch_device()) + for input, target in zip(*targets): + self._optimizer.zero_grad() -class ForwardDynamicsInterface(CasadiFunctionInterface): - def __init__(self, model: BiorbdModel, opts={}): - self.non_casadi_model = biorbd.Model(model.path) - super(ForwardDynamicsInterface, self).__init__("ForwardDynamicsInterface", opts) + # Get the predictions and targets + output = self(input[None, :]) - def inputs_len(self) -> list[int]: - return [1] + # Do some machine learning shenanigans + current_loss = self._loss_function.forward(output, target[None, :]) + current_loss.backward() # Backpropagation + self._optimizer.step() # Updating weights - def outputs_len(self) -> list[int]: - return [1] + # Populate the return values + all_predictions = torch.cat((all_predictions, output)) + all_targets = torch.cat((all_targets, target[None, :])) - def function(self, *args): - return [args[0]] # self.non_casadi_model.ForwardDynamics(*self.mx_in()[:3]) + # Put back the model in evaluation mode + self.eval() - def jacobians(self, *args): - return [0, 0, DM(1), 0, 0] - - -def custom_func_track_markers(controller: PenaltyController) -> MX: - return controller.model.custom_interface(controller.states["q"].cx, controller.controls["tau"].cx) + # Calculation of mean distance and error % + epoch_accuracy = (all_predictions - all_targets).abs().mean().item() + return epoch_accuracy def prepare_ocp( - biorbd_model_path: str, + model: torch.nn.Module, final_time: float, n_shooting: int, - ode_solver: OdeSolverBase = OdeSolver.RK4(), - use_sx: bool = False, - n_threads: int = 1, - phase_dynamics: PhaseDynamics = PhaseDynamics.SHARED_DURING_THE_PHASE, - expand_dynamics: bool = True, - control_type: ControlType = ControlType.CONSTANT, ) -> OptimalControlProgram: - """ - The initialization of an ocp - - Parameters - ---------- - biorbd_model_path: str - The path to the biorbd model - final_time: float - The time in second required to perform the task - n_shooting: int - The number of shooting points to define int the direct multiple shooting program - ode_solver: OdeSolverBase = OdeSolver.RK4() - Which type of OdeSolver to use - use_sx: bool - If the SX variable should be used instead of MX (can be extensive on RAM) - n_threads: int - The number of threads to use in the paralleling (1 = no parallel computing) - phase_dynamics: PhaseDynamics - If the dynamics equation within a phase is unique or changes at each node. - PhaseDynamics.SHARED_DURING_THE_PHASE is much faster, but lacks the capability to have changing dynamics within - a phase. A good example of when PhaseDynamics.ONE_PER_NODE should be used is when different external forces - are applied at each node - expand_dynamics: bool - If the dynamics function should be expanded. Please note, this will solve the problem faster, but will slow down - the declaration of the OCP, so it is a trade-off. Also depending on the solver, it may or may not work - (for instance IRK is not compatible with expanded dynamics) - control_type: ControlType - The type of the controls - - Returns - ------- - The OptimalControlProgram ready to be solved - """ - - bio_model = BiorbdModel(biorbd_model_path) - bio_model.custom_interface = CasadiFunctionInterfaceTest(bio_model) - - # Add objective functions - objective_functions = Objective(custom_func_track_markers, custom_type=ObjectiveFcn.Mayer, node=Node.START) # Dynamics - dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN, expand_dynamics=expand_dynamics, phase_dynamics=phase_dynamics) + dynamics = Dynamics(DynamicsFcn.TORQUE_DRIVEN) + torch_model = TorchModel(torch_model=model) # Path bounds x_bounds = BoundsList() - x_bounds["q"] = bio_model.bounds_from_ranges("q") + x_bounds["q"] = [-3.14 * 1.5] * torch_model.nb_q, [3.14 * 1.5] * torch_model.nb_q x_bounds["q"][:, [0, -1]] = 0 # Start and end at 0... x_bounds["q"][1, -1] = 3.14 # ...but end with pendulum 180 degrees rotated - x_bounds["qdot"] = bio_model.bounds_from_ranges("qdot") + x_bounds["qdot"] = [-3.14 * 10.0] * torch_model.nb_qdot, [3.14 * 10.0] * torch_model.nb_qdot x_bounds["qdot"][:, [0, -1]] = 0 # Start and end without any velocity - # Initial guess (optional since it is 0, we show how to initialize anyway) - x_init = InitialGuessList() - x_init["q"] = [0] * bio_model.nb_q - x_init["qdot"] = [0] * bio_model.nb_qdot - # Define control path bounds - n_tau = bio_model.nb_tau u_bounds = BoundsList() - u_bounds["tau"] = [-100] * n_tau, [100] * n_tau # Limit the strength of the pendulum to (-100 to 100)... + u_bounds["tau"] = [-100] * torch_model.nb_tau, [100] * torch_model.nb_tau u_bounds["tau"][1, :] = 0 # ...but remove the capability to actively rotate - # Initial guess (optional since it is 0, we show how to initialize anyway) - u_init = InitialGuessList() - u_init["tau"] = [0] * n_tau - return OptimalControlProgram( - bio_model, + torch_model, dynamics, n_shooting, final_time, - x_init=x_init, - u_init=u_init, x_bounds=x_bounds, u_bounds=u_bounds, - objective_functions=objective_functions, - ode_solver=ode_solver, - use_sx=use_sx, - n_threads=n_threads, - control_type=control_type, + use_sx=True, ) +def generate_dataset(biorbd_model: biorbd.Model, data_point_count: int) -> list[torch.Tensor]: + q = torch.rand(data_point_count, biorbd_model.nbQ()) + qdot = torch.rand(data_point_count, biorbd_model.nbQdot()) + tau = torch.rand(data_point_count, biorbd_model.nbGeneralizedTorque()) + + qddot = torch.zeros(data_point_count, biorbd_model.nbQddot()) + for i in range(data_point_count): + qddot[i, :] = torch.tensor( + biorbd_model.ForwardDynamics(np.array(q[i, :]), np.array(qdot[i, :]), np.array(tau[i, :])).to_array() + ) + + return [torch.cat((q, qdot, tau), dim=1), qddot] + + def main(): - """ - If pendulum is run as a script, it will perform the optimization and animates it - """ + # --- Prepare a predictive model --- # + biorbd_model = biorbd.Model("models/pendulum.bioMod") + training_data = generate_dataset(biorbd_model, data_point_count=1000) + validation_data = generate_dataset(biorbd_model, data_point_count=100) + + model = NeuralNetworkModel(layer_node_count=(6, 10, 10, 2), dropout_probability=0.2, use_batch_norm=True) + model.train_me(training_data, validation_data) - # --- Prepare the ocp --- # - ocp = prepare_ocp(biorbd_model_path="models/pendulum.bioMod", final_time=1, n_shooting=400, n_threads=2) + ocp = prepare_ocp(model=model, final_time=1, n_shooting=40) # --- Solve the ocp --- # sol = ocp.solve() # --- Show the results graph --- # # sol.print_cost() - sol.graphs(show_bounds=True, save_name="results.png") + sol.graphs(show_bounds=True) if __name__ == "__main__": main() - - -######## OCP FAST ######## -# from casadi import * - -# T = 10.0 # Time horizon -# N = 20 # number of control intervals - -# # Declare model variables -# x1 = MX.sym("x1") -# x2 = MX.sym("x2") -# x = vertcat(x1, x2) -# u = MX.sym("u") - -# # Model equations -# xdot = vertcat((1 - x2**2) * x1 - x2 + u, x1) - - -# # Formulate discrete time dynamics -# if False: -# # CVODES from the SUNDIALS suite -# dae = {"x": x, "p": u, "ode": xdot} -# opts = {"tf": T / N} -# F = integrator("F", "cvodes", dae, opts) -# else: -# # Fixed step Runge-Kutta 4 integrator -# M = 4 # RK4 steps per interval -# DT = T / N / M -# f = Function("f", [x, u], [xdot]) -# X0 = MX.sym("X0", 2) -# U = MX.sym("U") -# X = X0 -# Q = 0 -# for j in range(M): -# k1 = f(X, U) -# k2 = f(X + DT / 2 * k1, U) -# k3 = f(X + DT / 2 * k2, U) -# k4 = f(X + DT * k3, U) -# X = X + DT / 6 * (k1 + 2 * k2 + 2 * k3 + k4) -# F = Function("F", [X0, U], [X], ["x0", "p"], ["xf"]) - -# # Start with an empty NLP -# w = [] -# w0 = [] -# lbw = [] -# ubw = [] -# g = [] -# lbg = [] -# ubg = [] - -# # "Lift" initial conditions -# Xk = MX.sym("X0", 2) -# w += [Xk] -# lbw += [0, 1] -# ubw += [0, 1] -# w0 += [0, 1] - -# # Formulate the NLP -# for k in range(N): -# # New NLP variable for the control -# Uk = MX.sym("U_" + str(k)) -# w += [Uk] -# lbw += [-1] -# ubw += [1] -# w0 += [0] - -# # Integrate till the end of the interval -# Fk = F(x0=Xk, p=Uk) -# Xk_end = Fk["xf"] - -# # New NLP variable for state at end of interval -# Xk = MX.sym("X_" + str(k + 1), 2) -# w += [Xk] -# lbw += [-0.25, -inf] -# ubw += [inf, inf] -# w0 += [0, 0] - -# # Add equality constraint -# g += [Xk_end - Xk] -# lbg += [0, 0] -# ubg += [0, 0] - -# nd = N + 1 - -# import gpflow -# import time - -# from tensorflow_casadi import TensorFlowEvaluator - - -# class GPR(TensorFlowEvaluator): -# def __init__(self, session, opts={}): -# X = tf.compat.v1.placeholder(shape=(1, nd), dtype=np.float64) -# mean = tf.reshape(tf.reduce_mean(X), (1, 1)) -# TensorFlowEvaluator.__init__(self, [X], [mean], session, opts) -# self.counter = 0 -# self.time = 0 - -# def eval(self, arg): -# self.counter += 1 -# t0 = time.time() -# ret = TensorFlowEvaluator.eval(self, arg) -# self.time += time.time() - t0 -# return [ret] - - -# import tensorflow as tf - -# with tf.compat.v1.Session() as session: -# GPR = GPR(session) - -# w = vertcat(*w) - -# # Create an NLP solver -# prob = {"f": sum1(GPR(w[0::3])), "x": w, "g": vertcat(*g)} -# options = {"ipopt": {"hessian_approximation": "limited-memory"}} -# solver = nlpsol("solver", "ipopt", prob, options) - -# # Solve the NLP -# sol = solver(x0=w0, lbx=lbw, ubx=ubw, lbg=lbg, ubg=ubg) - -# print("Ncalls", GPR.counter) -# print("Total time [s]", GPR.time) -# w_opt = sol["x"].full().flatten() - -# # Plot the solution -# x1_opt = w_opt[0::3] -# x2_opt = w_opt[1::3] -# u_opt = w_opt[2::3] - -# tgrid = [T / N * k for k in range(N + 1)] -# import matplotlib.pyplot as plt - -# plt.figure(1) -# plt.clf() -# plt.plot(tgrid, x1_opt, "--") -# plt.plot(tgrid, x2_opt, "-") -# plt.step(tgrid, vertcat(DM.nan(1), u_opt), "-.") -# plt.xlabel("t") -# plt.legend(["x1", "x2", "u"]) -# plt.grid() -# plt.show() - - -# -# -# -######### TENSORFLOW CASADI ######### -# import casadi -# import tensorflow as tf - - -# class TensorFlowEvaluator(casadi.Callback): -# def __init__(self, t_in, t_out, session, opts={}): -# """ -# t_in: list of inputs (tensorflow placeholders) -# t_out: list of outputs (tensors dependeant on those placeholders) -# session: a tensorflow session -# """ -# casadi.Callback.__init__(self) -# assert isinstance(t_in, list) -# self.t_in = t_in -# assert isinstance(t_out, list) -# self.t_out = t_out -# self.construct("TensorFlowEvaluator", opts) -# self.session = session -# self.refs = [] - -# def get_n_in(self): -# return len(self.t_in) - -# def get_n_out(self): -# return len(self.t_out) - -# def get_sparsity_in(self, i): -# return casadi.Sparsity.dense(*self.t_in[i].get_shape().as_list()) - -# def get_sparsity_out(self, i): -# return casadi.Sparsity.dense(*self.t_out[i].get_shape().as_list()) - -# def eval(self, arg): -# # Associate each tensorflow input with the numerical argument passed by CasADi -# d = dict((v, arg[i].toarray()) for i, v in enumerate(self.t_in)) -# # Evaluate the tensorflow expressions -# ret = self.session.run(self.t_out, feed_dict=d) -# return ret - -# # Vanilla tensorflow offers just the reverse mode AD -# def has_reverse(self, nadj): -# return nadj == 1 - -# def get_reverse(self, nadj, name, inames, onames, opts): -# # Construct tensorflow placeholders for the reverse seeds -# adj_seed = [ -# tf.compat.v1.placeholder(shape=self.sparsity_out(i).shape, dtype=tf.float64) for i in range(self.n_out()) -# ] -# # Construct the reverse tensorflow graph through 'gradients' -# grad = tf.gradients(self.t_out, self.t_in, grad_ys=adj_seed) -# # Create another TensorFlowEvaluator object -# callback = TensorFlowEvaluator(self.t_in + adj_seed, grad, self.session) -# # Make sure you keep a reference to it -# self.refs.append(callback) - -# # Package it in the nominal_in+nominal_out+adj_seed form that CasADi expects -# nominal_in = self.mx_in() -# nominal_out = self.mx_out() -# adj_seed = self.mx_out() -# return casadi.Function( -# name, nominal_in + nominal_out + adj_seed, callback.call(nominal_in + adj_seed), inames, onames -# ) - - -# if __name__ == "__main__": -# from casadi import * - -# a = tf.compat.v1.placeholder(shape=(2, 2), dtype=tf.float64) -# b = tf.compat.v1.placeholder(shape=(2, 1), dtype=tf.float64) - -# y = tf.matmul(tf.sin(a), b) - -# with tf.compat.v1.Session() as session: -# f_tf = TensorFlowEvaluator([a, b], [y], session) - -# a = MX.sym("a", 2, 2) -# b = MX.sym("a", 2, 1) -# y = f_tf(a, b) -# yref = mtimes(sin(a), b) - -# f = Function("f", [a, b], [y]) -# fref = Function("f", [a, b], [yref]) - -# print(f(DM([[1, 2], [3, 4]]), DM([[1], [3]]))) -# print(fref(DM([[1, 2], [3, 4]]), DM([[1], [3]]))) - -# f = Function("f", [a, b], [jacobian(y, a)]) -# fref = Function("f", [a, b], [jacobian(yref, a)]) -# print(f(DM([[1, 2], [3, 4]]), DM([[1], [3]]))) -# print(fref(DM([[1, 2], [3, 4]]), DM([[1], [3]]))) diff --git a/bioptim/models/torch/torch_model.py b/bioptim/models/torch/torch_model.py index 0180bd35a..abdfcf28b 100644 --- a/bioptim/models/torch/torch_model.py +++ b/bioptim/models/torch/torch_model.py @@ -1,20 +1,24 @@ from typing import Callable from casadi import SX, MX, vertcat, horzcat, norm_fro, Function +import l4casadi as l4c import numpy as np import torch -""" -INSTALLATION: -First, make sure pytorch is installed - - pip install torch>=2.0 --index-url https://download.pytorch.org/whl/cpu/torch_stable.html +# """ +# INSTALLATION: +# First, make sure pytorch is installed -Then, install l4casadi as the interface between CasADi and PyTorch +# pip install torch>=2.0 --index-url https://download.pytorch.org/whl/cpu/torch_stable.html +# setuptools>=68.1 +# scikit-build>=0.17 +# cmake>=3.27 +# ninja>=1.11 +# Then, install l4casadi as the interface between CasADi and PyTorch - pip install l4casadi --no-build-isolation +# pip install l4casadi --no-build-isolation -""" +# """ class TorchModel: @@ -22,1050 +26,76 @@ class TorchModel: This class wraps a pytorch model and allows the user to call some useful functions on it. """ - def __init__(self, model: str | torch.nn.Module): - if not isinstance(bio_model, str) and not isinstance(bio_model, biorbd.Model): - raise ValueError("The model should be of type 'str' or 'biorbd.Model'") + def __init__(self, torch_model: torch.nn.Module): + self._dynamic_model = l4c.L4CasADi(torch_model, device="cpu") # device='cuda' for GPU - self.model = biorbd.Model(bio_model) if isinstance(bio_model, str) else bio_model - if parameters is not None: - for param_key in parameters: - parameters[param_key].apply_parameter(self) - self._friction_coefficients = friction_coefficients - - self.external_force_set = ( - self._set_external_force_set(external_force_set) if external_force_set is not None else None - ) + self._nb_dof = torch_model.size_in // 3 self._symbolic_variables() - self.biorbd_external_forces_set = self._dispatch_forces() if external_force_set else None - - # TODO: remove mx (the MX parameters should be created inside the BiorbdModel) - self.parameters = parameters.mx if parameters else MX() def _symbolic_variables(self): """Declaration of MX variables of the right shape for the creation of CasADi Functions""" - self.q = MX.sym("q_mx", self.nb_q, 1) - self.qdot = MX.sym("qdot_mx", self.nb_qdot, 1) - self.qddot = MX.sym("qddot_mx", self.nb_qddot, 1) - self.qddot_joints = MX.sym("qddot_joints_mx", self.nb_qddot - self.nb_root, 1) - self.tau = MX.sym("tau_mx", self.nb_tau, 1) - self.muscle = MX.sym("muscle_mx", self.nb_muscles, 1) - self.activations = MX.sym("activations_mx", self.nb_muscles, 1) - self.external_forces = MX.sym( - "external_forces_mx", - self.external_force_set.nb_external_forces_components if self.external_force_set else 0, - 1, - ) - - def _set_external_force_set(self, external_force_set: ExternalForceSetTimeSeries): - """ - It checks the external forces and binds them to the model. - """ - external_force_set._check_segment_names(tuple([s.name().to_string() for s in self.model.segments()])) - external_force_set._check_all_string_points_of_application(self.marker_names) - external_force_set._bind() - - return external_force_set + self.q = MX.sym("q_mx", self.nb_dof, 1) + self.qdot = MX.sym("qdot_mx", self.nb_dof, 1) + self.tau = MX.sym("tau_mx", self.nb_dof, 1) + self.external_forces = MX.sym("external_forces_mx", 0, 1) + self.parameters = MX.sym("parameters_mx", 0, 1) @property def name(self) -> str: # parse the path and split to get the .bioMod name - return self.model.path().absolutePath().to_string().split("/")[-1] - - @property - def path(self) -> str: - return self.model.path().relativePath().to_string() - - def copy(self): - return BiorbdModel(self.path) - - def serialize(self) -> tuple[Callable, dict]: - return BiorbdModel, dict(bio_model=self.path) - - @property - def friction_coefficients(self) -> MX | SX | np.ndarray: - return self._friction_coefficients - - def set_friction_coefficients(self, new_friction_coefficients) -> None: - if np.any(new_friction_coefficients < 0): - raise ValueError("Friction coefficients must be positive") - return self._friction_coefficients - - @property - def gravity(self) -> Function: - """ - Returns the gravity of the model. - Since the gravity is self-defined in the model, you need to provide the type of the output when calling the function like this: - model.gravity()(MX() / SX()) - """ - biorbd_return = self.model.getGravity().to_mx() - casadi_fun = Function( - "gravity", - [self.parameters], - [biorbd_return], - ["parameters"], - ["gravity"], - ) - return casadi_fun - - def set_gravity(self, new_gravity) -> None: - self.model.setGravity(new_gravity) - return - - @property - def nb_tau(self) -> int: - return self.model.nbGeneralizedTorque() - - @property - def nb_segments(self) -> int: - return self.model.nbSegment() - - def segment_index(self, name) -> int: - return biorbd.segment_index(self.model, name) + return "forward_dynamics_torch_model" @property - def nb_quaternions(self) -> int: - return self.model.nbQuat() + def name_dof(self) -> list[str]: + return [f"q_{i}" for i in range(self.nb_dof)] @property def nb_dof(self) -> int: - return self.model.nbDof() + return self._nb_dof @property def nb_q(self) -> int: - return self.model.nbQ() + return self.nb_dof @property def nb_qdot(self) -> int: - return self.model.nbQdot() - - @property - def nb_qddot(self) -> int: - return self.model.nbQddot() - - @property - def nb_root(self) -> int: - return self.model.nbRoot() - - @property - def segments(self) -> tuple[biorbd.Segment]: - return self.model.segments() - - def rotation_matrix_to_euler_angles(self, sequence: str) -> Function: - """ - Returns the rotation matrix to euler angles function. - """ - r = MX.sym("r_mx", 3, 3) - r_matrix = biorbd.Rotation(r[0, 0], r[0, 1], r[0, 2], r[1, 0], r[1, 1], r[1, 2], r[2, 0], r[2, 1], r[2, 2]) - biorbd_return = biorbd.Rotation.toEulerAngles(r_matrix, sequence).to_mx() - casadi_fun = Function( - "rotation_matrix_to_euler_angles", - [r], - [biorbd_return], - ["Rotation matrix"], - ["Euler angles"], - ) - return casadi_fun - - def homogeneous_matrices_in_global(self, segment_index: int, inverse=False) -> Function: - """ - Returns the roto-translation matrix of the segment in the global reference frame. - """ - q_biorbd = GeneralizedCoordinates(self.q) - jcs = self.model.globalJCS(q_biorbd, segment_index) - biorbd_return = jcs.transpose().to_mx() if inverse else jcs.to_mx() - casadi_fun = Function( - "homogeneous_matrices_in_global", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["Joint coordinate system RT matrix in global"], - ) - return casadi_fun - - def homogeneous_matrices_in_child(self, segment_id) -> Function: - """ - Returns the roto-translation matrix of the segment in the child reference frame. - Since the homogeneous matrix is self-defined in the model, you need to provide the type of the output when calling the function like this: - model.homogeneous_matrices_in_child(segment_id)(MX() / SX()) - """ - biorbd_return = self.model.localJCS(segment_id).to_mx() - casadi_fun = Function( - "homogeneous_matrices_in_child", - [self.parameters], - [biorbd_return], - ["parameters"], - ["Joint coordinate system RT matrix in local"], - ) - return casadi_fun - - @property - def mass(self) -> Function: - """ - Returns the mass of the model. - Since the mass is self-defined in the model, you need to provide the type of the output when calling the function like this: - model.mass()(MX() / SX()) - """ - biorbd_return = self.model.mass().to_mx() - casadi_fun = Function( - "mass", - [self.parameters], - [biorbd_return], - ["parameters"], - ["mass"], - ) - return casadi_fun - - def rt(self, rt_index) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - biorbd_return = self.model.RT(q_biorbd, rt_index).to_mx() - casadi_fun = Function( - "rt", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["RT matrix"], - ) - return casadi_fun - - def center_of_mass(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - biorbd_return = self.model.CoM(q_biorbd, True).to_mx() - casadi_fun = Function( - "center_of_mass", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["Center of mass"], - ) - return casadi_fun - - def center_of_mass_velocity(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.CoMdot(q_biorbd, qdot_biorbd, True).to_mx() - casadi_fun = Function( - "center_of_mass_velocity", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["Center of mass velocity"], - ) - return casadi_fun - - def center_of_mass_acceleration(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - qddot_biorbd = GeneralizedAcceleration(self.qddot) - biorbd_return = self.model.CoMddot(q_biorbd, qdot_biorbd, qddot_biorbd, True).to_mx() - casadi_fun = Function( - "center_of_mass_acceleration", - [self.q, self.qdot, self.qddot, self.parameters], - [biorbd_return], - ["q", "qdot", "qddot", "parameters"], - ["Center of mass acceleration"], - ) - return casadi_fun - - def body_rotation_rate(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.bodyAngularVelocity(q_biorbd, qdot_biorbd, True).to_mx() - casadi_fun = Function( - "body_rotation_rate", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["Body rotation rate"], - ) - return casadi_fun - - def mass_matrix(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - biorbd_return = self.model.massMatrix(q_biorbd).to_mx() - casadi_fun = Function( - "mass_matrix", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["Mass matrix"], - ) - return casadi_fun - - def non_linear_effects(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.NonLinearEffect(q_biorbd, qdot_biorbd).to_mx() - casadi_fun = Function( - "non_linear_effects", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["Non linear effects"], - ) - return casadi_fun - - def angular_momentum(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.angularMomentum(q_biorbd, qdot_biorbd, True).to_mx() - casadi_fun = Function( - "angular_momentum", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["Angular momentum"], - ) - return casadi_fun - - def reshape_qdot(self, k_stab=1) -> Function: - biorbd_return = self.model.computeQdot( - GeneralizedCoordinates(self.q), - GeneralizedCoordinates(self.qdot), # mistake in biorbd - k_stab, - ).to_mx() - casadi_fun = Function( - "reshape_qdot", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["Reshaped qdot"], - ) - return casadi_fun - - def segment_angular_velocity(self, idx) -> Function: - """ - Returns the angular velocity of the segment in the global reference frame. - """ - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.segmentAngularVelocity(q_biorbd, qdot_biorbd, idx, True).to_mx() - casadi_fun = Function( - "segment_angular_velocity", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["Segment angular velocity"], - ) - return casadi_fun - - def segment_orientation(self, idx: int, sequence: str = "xyz") -> Function: - """ - Returns the angular position of the segment in the global reference frame. - """ - q_biorbd = GeneralizedCoordinates(self.q) - rotation_matrix = self.homogeneous_matrices_in_global(idx)(q_biorbd, self.parameters)[:3, :3] - biorbd_return = self.rotation_matrix_to_euler_angles(sequence=sequence)(rotation_matrix) - casadi_fun = Function( - "segment_orientation", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["Segment orientation"], - ) - return casadi_fun - - @property - def name_dof(self) -> tuple[str, ...]: - return tuple(s.to_string() for s in self.model.nameDof()) - - @property - def contact_names(self) -> tuple[str, ...]: - return tuple(s.to_string() for s in self.model.contactNames()) - - @property - def nb_soft_contacts(self) -> int: - return self.model.nbSoftContacts() - - @property - def soft_contact_names(self) -> tuple[str, ...]: - return self.model.softContactNames() - - def soft_contact(self, soft_contact_index, *args): - return self.model.softContact(soft_contact_index, *args) + return self.nb_dof @property - def muscle_names(self) -> tuple[str, ...]: - return tuple(s.to_string() for s in self.model.muscleNames()) - - @property - def nb_muscles(self) -> int: - return self.model.nbMuscles() - - def torque(self) -> Function: - """ - Returns the torque from the torque_activations. - Note that tau_activation should be between 0 and 1. - """ - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - tau_activations_biorbd = self.tau - biorbd_return = self.model.torque(tau_activations_biorbd, q_biorbd, qdot_biorbd).to_mx() - casadi_fun = Function( - "torque_activation", - [self.tau, self.q, self.qdot, self.parameters], - [biorbd_return], - ["tau", "q", "qdot", "parameters"], - ["Torque from tau activations"], - ) - return casadi_fun - - def forward_dynamics_free_floating_base(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - qddot_joints_biorbd = GeneralizedAcceleration(self.qddot_joints) - biorbd_return = self.model.ForwardDynamicsFreeFloatingBase(q_biorbd, qdot_biorbd, qddot_joints_biorbd).to_mx() - casadi_fun = Function( - "forward_dynamics_free_floating_base", - [self.q, self.qdot, self.qddot_joints, self.parameters], - [biorbd_return], - ["q", "qdot", "qddot_joints", "parameters"], - ["qddot_root and qddot_joints"], - ) - return casadi_fun - - @staticmethod - def reorder_qddot_root_joints(qddot_root, qddot_joints) -> MX | SX: - return vertcat(qddot_root, qddot_joints) - - def _dispatch_forces(self) -> biorbd.ExternalForceSet: - """Dispatch the symbolic MX into the biorbd external forces object""" - biorbd_external_forces = self.model.externalForceSet() - - # "type of external force": (function to call, number of force components) - force_mapping = { - "in_global": (_add_global_force, 6), - "torque_in_global": (_add_torque_global, 3), - "translational_in_global": (_add_translational_global, 3), - "in_local": (_add_local_force, 6), - "torque_in_local": (_add_torque_local, 3), - } - - symbolic_counter = 0 - for force_type, val in force_mapping.items(): - add_force_func, num_force_components = val - symbolic_counter = self._dispatch_forces_of_type( - force_type, add_force_func, num_force_components, symbolic_counter, biorbd_external_forces - ) - - return biorbd_external_forces - - def _dispatch_forces_of_type( - self, - force_type: str, - add_force_func: "Callable", - num_force_components: int, - symbolic_counter: int, - biorbd_external_forces: "biorbd.ExternalForces", - ) -> int: - """ - Helper method to dispatch forces of a specific external forces. - - Parameters - ---------- - force_type: str - The type of external force to dispatch among in_global, torque_in_global, translational_in_global, in_local, torque_in_local. - add_force_func: Callable - The function to call to add the force to the biorbd external forces object. - num_force_components: int - The number of force components for the given type - symbolic_counter: int - The current symbolic counter to slice the whole external_forces mx. - biorbd_external_forces: biorbd.ExternalForces - The biorbd external forces object to add the forces to. - - Returns - ------- - int - The updated symbolic counter. - """ - for segment, forces_on_segment in getattr(self.external_force_set, force_type).items(): - for force in forces_on_segment: - force_slicer = slice(symbolic_counter, symbolic_counter + num_force_components) - - point_of_application_mx = self._get_point_of_application(force, force_slicer.stop) - - add_force_func( - biorbd_external_forces, segment, self.external_forces[force_slicer], point_of_application_mx - ) - symbolic_counter = force_slicer.stop + ( - 3 if isinstance(force["point_of_application"], np.ndarray) else 0 - ) - - return symbolic_counter - - def _get_point_of_application(self, force, stop_index) -> biorbd.NodeSegment | np.ndarray | None: - """ - Determine the point of application mx slice based on its type. Only sliced if an array is stored - - Parameters - ---------- - force : dict - The force dictionary with details on the point of application. - stop_index : int - Index position in MX where the point of application components start. - - Returns - ------- - biorbd.NodeSegment | np.ndarray | None - Returns a slice of MX, a marker node, or None if no point of application is defined. - """ - if isinstance(force["point_of_application"], np.ndarray): - return self.external_forces[slice(stop_index, stop_index + 3)] - elif isinstance(force["point_of_application"], str): - return self.model.marker(self.marker_index(force["point_of_application"])) - return None + def nb_tau(self) -> int: + return self.nb_dof def forward_dynamics(self, with_contact: bool = False) -> Function: - - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - tau_biorbd = GeneralizedTorque(self.tau) - - if with_contact: - if self.external_force_set is None: - biorbd_return = self.model.ForwardDynamicsConstraintsDirect(q_biorbd, qdot_biorbd, tau_biorbd).to_mx() - else: - biorbd_return = self.model.ForwardDynamicsConstraintsDirect( - q_biorbd, qdot_biorbd, tau_biorbd, self.biorbd_external_forces_set - ).to_mx() - casadi_fun = Function( - "constrained_forward_dynamics", - [self.q, self.qdot, self.tau, self.external_forces, self.parameters], - [biorbd_return], - ["q", "qdot", "tau", "external_forces", "parameters"], - ["qddot"], - ) - else: - if self.external_force_set is None: - biorbd_return = self.model.ForwardDynamics(q_biorbd, qdot_biorbd, tau_biorbd).to_mx() - else: - biorbd_return = self.model.ForwardDynamics( - q_biorbd, qdot_biorbd, tau_biorbd, self.biorbd_external_forces_set - ).to_mx() - casadi_fun = Function( - "forward_dynamics", - [self.q, self.qdot, self.tau, self.external_forces, self.parameters], - [biorbd_return], - ["q", "qdot", "tau", "external_forces", "parameters"], - ["qddot"], - ) - return casadi_fun - - def inverse_dynamics(self, with_contact: bool = False) -> Function: - - if with_contact: - raise NotImplementedError("Inverse dynamics with contact is not implemented yet") - - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - qddot_biorbd = GeneralizedAcceleration(self.qddot) - if self.external_force_set is None: - biorbd_return = self.model.InverseDynamics(q_biorbd, qdot_biorbd, qddot_biorbd).to_mx() - else: - biorbd_return = self.model.InverseDynamics( - q_biorbd, qdot_biorbd, qddot_biorbd, self.biorbd_external_forces_set - ).to_mx() - casadi_fun = Function( - "inverse_dynamics", - [self.q, self.qdot, self.qddot, self.external_forces, self.parameters], - [biorbd_return], - ["q", "qdot", "qddot", "external_forces", "parameters"], - ["tau"], - ) - return casadi_fun - - def contact_forces_from_constrained_forward_dynamics(self) -> Function: - - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - tau_biorbd = GeneralizedTorque(self.tau) - if self.external_force_set is None: - biorbd_return = self.model.ContactForcesFromForwardDynamicsConstraintsDirect( - q_biorbd, qdot_biorbd, tau_biorbd - ).to_mx() - else: - biorbd_return = self.model.ContactForcesFromForwardDynamicsConstraintsDirect( - q_biorbd, qdot_biorbd, tau_biorbd, self.biorbd_external_forces_set - ).to_mx() - casadi_fun = Function( - "contact_forces_from_constrained_forward_dynamics", + return Function( + "forward_dynamics", [self.q, self.qdot, self.tau, self.external_forces, self.parameters], - [biorbd_return], + [self._dynamic_model(vertcat(self.q, self.qdot, self.tau).T).T], ["q", "qdot", "tau", "external_forces", "parameters"], - ["contact_forces"], - ) - return casadi_fun - - def qdot_from_impact(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_pre_impact_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.ComputeConstraintImpulsesDirect(q_biorbd, qdot_pre_impact_biorbd).to_mx() - casadi_fun = Function( - "qdot_from_impact", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["qdot post impact"], - ) - return casadi_fun - - def muscle_activation_dot(self) -> Function: - muscle_states = self.model.stateSet() - for k in range(self.model.nbMuscles()): - muscle_states[k].setActivation(self.activations[k]) - muscle_states[k].setExcitation(self.muscle[k]) - biorbd_return = self.model.activationDot(muscle_states).to_mx() - casadi_fun = Function( - "muscle_activation_dot", - [self.muscle, self.activations, self.parameters], - [biorbd_return], - ["muscle_excitation", "muscle_activation", "parameters"], - ["muscle_activation_dot"], - ) - return casadi_fun - - def muscle_length_jacobian(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - biorbd_return = self.model.musclesLengthJacobian(q_biorbd).to_mx() - casadi_fun = Function( - "muscle_length_jacobian", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["muscle_length_jacobian"], - ) - return casadi_fun - - def muscle_velocity(self) -> Function: - J = self.muscle_length_jacobian()(self.q, self.parameters) - biorbd_return = J @ self.qdot - casadi_fun = Function( - "muscle_velocity", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["muscle_velocity"], - ) - return casadi_fun - - def muscle_joint_torque(self) -> Function: - muscles_states = self.model.stateSet() - muscles_activations = self.muscle - for k in range(self.model.nbMuscles()): - muscles_states[k].setActivation(muscles_activations[k]) - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.muscularJointTorque(muscles_states, q_biorbd, qdot_biorbd).to_mx() - casadi_fun = Function( - "muscle_joint_torque", - [self.muscle, self.q, self.qdot, self.parameters], - [biorbd_return], - ["muscle_activation", "q", "qdot", "parameters"], - ["muscle_joint_torque"], - ) - return casadi_fun - - def markers(self) -> list[MX]: - biorbd_return = horzcat(*[m.to_mx() for m in self.model.markers(GeneralizedCoordinates(self.q))]) - casadi_fun = Function( - "markers", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["markers"], - ) - return casadi_fun - - @property - def nb_markers(self) -> int: - return self.model.nbMarkers() - - def marker_index(self, name): - return biorbd.marker_index(self.model, name) - - def marker(self, index: int, reference_segment_index: int = None) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - marker = self.model.marker(q_biorbd, index) - if reference_segment_index is not None: - global_homogeneous_matrix = self.model.globalJCS(q_biorbd, reference_segment_index) - marker_rotated = global_homogeneous_matrix.transpose().to_mx() @ vertcat(marker.to_mx(), 1) - biorbd_return = marker_rotated[:3] - else: - biorbd_return = marker.to_mx() - casadi_fun = Function( - "marker", - [self.q, self.parameters], - [biorbd_return], - ["q", "parameters"], - ["marker"], - ) - return casadi_fun - - @property - def nb_rigid_contacts(self) -> int: - """ - Returns the number of rigid contacts. - Example: - First contact with axis YZ - Second contact with axis Z - nb_rigid_contacts = 2 - """ - return self.model.nbRigidContacts() + ["qddot"], + ).expand() @property def nb_contacts(self) -> int: - """ - Returns the number of contact index. - Example: - First contact with axis YZ - Second contact with axis Z - nb_contacts = 3 - """ - return self.model.nbContacts() - - def rigid_contact_index(self, contact_index) -> tuple: - """ - Returns the axis index of this specific rigid contact. - Example: - First contact with axis YZ - Second contact with axis Z - rigid_contact_index(0) = (1, 2) - """ - return self.model.rigidContacts()[contact_index].availableAxesIndices() - - def markers_velocities(self, reference_index=None) -> list[MX]: - if reference_index is None: - biorbd_return = [ - m.to_mx() - for m in self.model.markersVelocity( - GeneralizedCoordinates(self.q), - GeneralizedVelocity(self.qdot), - True, - ) - ] - - else: - biorbd_return = [] - homogeneous_matrix_transposed = self.homogeneous_matrices_in_global( - segment_index=reference_index, inverse=True - )( - GeneralizedCoordinates(self.q), - ) - for m in self.model.markersVelocity(GeneralizedCoordinates(self.q), GeneralizedVelocity(self.qdot)): - if m.applyRT(homogeneous_matrix_transposed) is None: - biorbd_return.append(m.to_mx()) - else: - biorbd_return.append(m.applyRT(homogeneous_matrix_transposed).to_mx()) - - casadi_fun = Function( - "markers_velocities", - [self.q, self.qdot, self.parameters], - [horzcat(*biorbd_return)], - ["q", "qdot", "parameters"], - ["markers_velocities"], - ) - return casadi_fun - - def marker_velocity(self, marker_index: int) -> list[MX]: - biorbd_return = self.model.markersVelocity( - GeneralizedCoordinates(self.q), - GeneralizedVelocity(self.qdot), - True, - )[marker_index].to_mx() - casadi_fun = Function( - "marker_velocity", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["marker_velocity"], - ) - return casadi_fun - - def markers_accelerations(self, reference_index=None) -> list[MX]: - if reference_index is None: - biorbd_return = [ - m.to_mx() - for m in self.model.markerAcceleration( - GeneralizedCoordinates(self.q), - GeneralizedVelocity(self.qdot), - GeneralizedAcceleration(self.qddot), - True, - ) - ] - - else: - biorbd_return = [] - homogeneous_matrix_transposed = self.homogeneous_matrices_in_global( - segment_index=reference_index, - inverse=True, - )( - GeneralizedCoordinates(self.q), - ) - for m in self.model.markersAcceleration( - GeneralizedCoordinates(self.q), - GeneralizedVelocity(self.qdot), - GeneralizedAcceleration(self.qddot), - ): - if m.applyRT(homogeneous_matrix_transposed) is None: - biorbd_return.append(m.to_mx()) - else: - biorbd_return.append(m.applyRT(homogeneous_matrix_transposed).to_mx()) + return 0 - casadi_fun = Function( - "markers_accelerations", - [self.q, self.qdot, self.qddot, self.parameters], - [horzcat(*biorbd_return)], - ["q", "qdot", "qddot", "parameters"], - ["markers_accelerations"], - ) - return casadi_fun - - def marker_acceleration(self, marker_index: int) -> list[MX]: - biorbd_return = self.model.markerAcceleration( - GeneralizedCoordinates(self.q), - GeneralizedVelocity(self.qdot), - GeneralizedAcceleration(self.qddot), - True, - )[marker_index].to_mx() - casadi_fun = Function( - "marker_acceleration", - [self.q, self.qdot, self.qddot, self.parameters], - [biorbd_return], - ["q", "qdot", "qddot", "parameters"], - ["marker_acceleration"], - ) - return casadi_fun + @property + def nb_soft_contacts(self) -> int: + return 0 - def tau_max(self) -> tuple[MX, MX]: - self.model.closeActuator() - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - torque_max, torque_min = self.model.torqueMax(q_biorbd, qdot_biorbd) - casadi_fun = Function( - "tau_max", + def reshape_qdot(self, k_stab=1) -> Function: + return Function( + "reshape_qdot", [self.q, self.qdot, self.parameters], - [torque_max.to_mx(), torque_min.to_mx()], + [self.qdot], ["q", "qdot", "parameters"], - ["tau_max", "tau_min"], - ) - return casadi_fun - - def rigid_contact_acceleration(self, contact_index, contact_axis) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - qddot_biorbd = GeneralizedAcceleration(self.qddot) - biorbd_return = self.model.rigidContactAcceleration( - q_biorbd, qdot_biorbd, qddot_biorbd, contact_index, True - ).to_mx()[contact_axis] - casadi_fun = Function( - "rigid_contact_acceleration", - [self.q, self.qdot, self.qddot, self.parameters], - [biorbd_return], - ["q", "qdot", "qddot", "parameters"], - ["rigid_contact_acceleration"], - ) - return casadi_fun - - def markers_jacobian(self) -> list[MX]: - biorbd_return = [m.to_mx() for m in self.model.markersJacobian(GeneralizedCoordinates(self.q))] - casadi_fun = Function( - "markers_jacobian", - [self.q, self.parameters], - biorbd_return, - ["q", "parameters"], - ["markers_jacobian"], - ) - return casadi_fun - - @property - def marker_names(self) -> tuple[str, ...]: - return tuple([s.to_string() for s in self.model.markerNames()]) + ["Reshaped qdot"], + ).expand() def soft_contact_forces(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - - biorbd_return = MX.zeros(self.nb_soft_contacts * 6, 1) - for i_sc in range(self.nb_soft_contacts): - soft_contact = self.soft_contact(i_sc) - biorbd_return[i_sc * 6 : (i_sc + 1) * 6, :] = ( - biorbd.SoftContactSphere(soft_contact).computeForceAtOrigin(self.model, q_biorbd, qdot_biorbd).to_mx() - ) - - casadi_fun = Function( + return Function( "soft_contact_forces", [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["soft_contact_forces"], - ) - return casadi_fun - - def normalize_state_quaternions(self) -> Function: - - quat_idx = self.get_quaternion_idx() - biorbd_return = MX.zeros(self.nb_q) - biorbd_return[:] = self.q - - # Normalize quaternion, if needed - for j in range(self.nb_quaternions): - quaternion = vertcat( - self.q[quat_idx[j][3]], - self.q[quat_idx[j][0]], - self.q[quat_idx[j][1]], - self.q[quat_idx[j][2]], - ) - quaternion /= norm_fro(quaternion) - biorbd_return[quat_idx[j][0] : quat_idx[j][2] + 1] = quaternion[1:4] - biorbd_return[quat_idx[j][3]] = quaternion[0] - - casadi_fun = Function( - "normalize_state_quaternions", - [self.q], - [biorbd_return], - ["q"], - ["q_normalized"], - ) - return casadi_fun - - def get_quaternion_idx(self) -> list[list[int]]: - n_dof = 0 - quat_idx = [] - quat_number = 0 - for j in range(self.nb_segments): - if self.segments[j].isRotationAQuaternion(): - quat_idx.append([n_dof, n_dof + 1, n_dof + 2, self.nb_dof + quat_number]) - quat_number += 1 - n_dof += self.segments[j].nbDof() - return quat_idx - - def contact_forces(self) -> Function: - force = self.contact_forces_from_constrained_forward_dynamics()( - self.q, self.qdot, self.tau, self.external_forces, self.parameters - ) - casadi_fun = Function( - "contact_forces", - [self.q, self.qdot, self.tau, self.external_forces, self.parameters], - [force], - ["q", "qdot", "tau", "external_forces", "parameters"], - ["contact_forces"], - ) - return casadi_fun - - def passive_joint_torque(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.passiveJointTorque(q_biorbd, qdot_biorbd).to_mx() - casadi_fun = Function( - "passive_joint_torque", - [self.q, self.qdot, self.parameters], - [biorbd_return], - ["q", "qdot", "parameters"], - ["passive_joint_torque"], - ) - return casadi_fun - - def ligament_joint_torque(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.ligamentsJointTorque(q_biorbd, qdot_biorbd).to_mx() - casadi_fun = Function( - "ligament_joint_torque", - [self.q, self.qdot, self.parameters], - [biorbd_return], + [MX(0)], ["q", "qdot", "parameters"], - ["ligament_joint_torque"], - ) - return casadi_fun - - def ranges_from_model(self, variable: str): - ranges = [] - for segment in self.segments: - if "_joints" in variable: - if segment.parent().to_string().lower() != "root": - variable = variable.replace("_joints", "") - ranges += self._add_range(variable, segment) - elif "_roots" in variable: - if segment.parent().to_string().lower() == "root": - variable = variable.replace("_roots", "") - ranges += self._add_range(variable, segment) - else: - ranges += self._add_range(variable, segment) - - return ranges - - @staticmethod - def _add_range(variable: str, segment: biorbd.Segment) -> list[biorbd.Range]: - """ - Get the range of a variable for a given segment - - Parameters - ---------- - variable: str - The variable to get the range for such as: - "q", "qdot", "qddot", "q_joint", "qdot_joint", "qddot_joint", "q_root", "qdot_root", "qddot_root" - segment: biorbd.Segment - The segment to get the range from - - Returns - ------- - list[biorbd.Range] - range min and max for the given variable for a given segment - """ - ranges_map = { - "q": [q_range for q_range in segment.QRanges()], - "qdot": [qdot_range for qdot_range in segment.QdotRanges()], - "qddot": [qddot_range for qddot_range in segment.QddotRanges()], - } - - segment_variable_range = ranges_map.get(variable, None) - if segment_variable_range is None: - RuntimeError("Wrong variable name") - - return segment_variable_range - - def _var_mapping( - self, - key: str, - range_for_mapping: int | list | tuple | range, - mapping: BiMapping = None, - ) -> dict: - return _var_mapping(key, range_for_mapping, mapping) - - def bounds_from_ranges(self, variables: str | list[str], mapping: BiMapping | BiMappingList = None) -> Bounds: - return bounds_from_ranges(self, variables, mapping) - - def lagrangian(self) -> Function: - q_biorbd = GeneralizedCoordinates(self.q) - qdot_biorbd = GeneralizedVelocity(self.qdot) - biorbd_return = self.model.Lagrangian(q_biorbd, qdot_biorbd).to_mx() - casadi_fun = Function( - "lagrangian", - [self.q, self.qdot], - [biorbd_return], - ["q", "qdot"], - ["lagrangian"], - ) - return casadi_fun - - def partitioned_forward_dynamics(self): - raise NotImplementedError("partitioned_forward_dynamics is not implemented for BiorbdModel") - - @staticmethod - def animate( - ocp, - solution, - show_now: bool = True, - show_tracked_markers: bool = False, - viewer: str = "pyorerun", - n_frames: int = 0, - **kwargs, - ): - if viewer == "bioviz": - from .viewer_bioviz import animate_with_bioviz_for_loop - - return animate_with_bioviz_for_loop(ocp, solution, show_now, show_tracked_markers, n_frames, **kwargs) - if viewer == "pyorerun": - from .viewer_pyorerun import animate_with_pyorerun - - return animate_with_pyorerun(ocp, solution, show_now, show_tracked_markers, **kwargs) + ["Soft contact forces"], + ).expand()