From 589f89c76803eb4b8934858c9f9f0edbe5778a34 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Mon, 23 Jun 2025 10:06:01 +0200 Subject: [PATCH] feat: start of lift, debugging, cleanup --- doc/img/ball.png | Bin 0 -> 28572 bytes doc/notes.typ | 10 +- doc/pres.typ | 16 +-- doc/rep.typ | 192 ++++++++++++++++++++++++++++++++ exm/ball.ml | 8 +- exm/dune | 6 +- exm/sin1x.ml | 72 ++++++++++++ exm/sin1x_der.ml | 79 +++++++++++++ exm/sincos.ml | 8 +- exm/zelus/ballz/ballz.ml | 131 ++++++++++++++++++++++ exm/zelus/ballz/ballz.ml.bak | 137 +++++++++++++++++++++++ exm/zelus/ballz/ballz.ml.bak2 | 135 ++++++++++++++++++++++ exm/zelus/ballz/ballz.zci | Bin 0 -> 6309 bytes exm/zelus/ballz/ballz.zls | 19 ++++ exm/zelus/sin_1_x.ml | 74 ++++++++++++ exm/zelus/sin_1_x.zci | Bin 0 -> 9262 bytes exm/zelus/sin_1_x.zls | 7 ++ exm/zelus/sincos/sincosz.ml | 45 ++++++++ exm/zelus/sincos/sincosz.zci | Bin 0 -> 1170 bytes exm/zelus/sincos/sincosz.zls | 4 + src/bin/lift.ml | 119 ++++++++++++++++++++ src/bin/main.ml | 42 +++++-- src/bin/output.ml | 38 ++++++- src/lib/common/debug.ml | 7 ++ src/lib/common/ztypes.ml | 151 +++++++++++++++++++++++++ src/lib/hsim/sim.ml | 8 +- src/lib/hsim/types.ml | 2 +- src/lib/solvers/illinois.ml | 21 ++-- src/lib/solvers/odexx.ml | 10 +- src/lib/solvers/statefulRK45.ml | 2 + src/lib/solvers/statefulZ.ml | 5 +- 31 files changed, 1297 insertions(+), 51 deletions(-) create mode 100644 doc/img/ball.png create mode 100644 doc/rep.typ create mode 100644 exm/sin1x.ml create mode 100644 exm/sin1x_der.ml create mode 100644 exm/zelus/ballz/ballz.ml create mode 100644 exm/zelus/ballz/ballz.ml.bak create mode 100644 exm/zelus/ballz/ballz.ml.bak2 create mode 100644 exm/zelus/ballz/ballz.zci create mode 100644 exm/zelus/ballz/ballz.zls create mode 100644 exm/zelus/sin_1_x.ml create mode 100644 exm/zelus/sin_1_x.zci create mode 100644 exm/zelus/sin_1_x.zls create mode 100644 exm/zelus/sincos/sincosz.ml create mode 100644 exm/zelus/sincos/sincosz.zci create mode 100644 exm/zelus/sincos/sincosz.zls create mode 100644 src/bin/lift.ml create mode 100644 src/lib/common/ztypes.ml diff --git a/doc/img/ball.png b/doc/img/ball.png new file mode 100644 index 0000000000000000000000000000000000000000..2d578f451610daa8fe0263d07157906b337dba8a GIT binary patch literal 28572 zcmeAS@N?(olHy`uVBq!ia0y~yU}|7sV0^&A#=yW}dhyN^1_lPs0*}aI1_r((Aj~*b zn@^g7L4m>3#WAE}&YQcHHK9-E+J88I_xhbXZ>kk^I8SJGaR&(seqY+XZE4w@YgxGq zZ$w^+Pha!)!pd7>(bkt&`stSCEZthVDj{~lu14;zhBw{}ii)q#&bxD`dh+gn{}L7~ zYSB5r^SjOSIp=?W{q)tU?a3VdAEq_McF!yOG^ctkb#QTUacN5o+s&Y)q@;98hpz!7 zz;>8%LPtl(i3AG$dYW$Y$rP)n^ZD(5Ffa&oAN}xRaevakKR+WjrEr?%-dgfxDaeA1mNj1U za(;h%Tk+u_yJ7jeJCo(5>d!60WWajo4kL`)>dLIukRqq>qn|-q@JTuI4{)O$sN- z{U$uYT2s9=3l|kw@bn*VmbETRxV9#8`JX5akwa57gFjpe_D}rx=O;tM?)Uqo>-Nik zzgN9pXu7(RlF}lRt*&Fo+(JS_ z=l=OK(>VRnQSo@0Pc{6uUj*j;y*GP9gpS7&p$;pj7scC?4qV~yKFYKIKL7UH%)9Th zF?6^k+dTWIdHi7k2g40@8@bPg{{H-rK3aVG{CVk$_T#s#-@SOza5by-{`={R%V5_!t$skPP*HnU|D$s& zJ(N1vy|G#nkqU~D@0;AGvWId|VVczSB~U6HkU`j~e;>&D0Fu^V(}T1*;@PiTjIGnmd^9eCmGJ)!}=>ww0l#*-R{Q$l~eg&+lBWt`oi5kUKL7up=O0(c|K&1HJ9B_tzNX>%ylTGt zKhKpv_*hYJJ+?gY>8YtZ-tYT;;mVbkK6(cJrBQCCad&?sHMFMeHFJi>HD_v%8 zKgkpG@iSHi-5ec(sOsqwbx(!svgDuSS_@tB-z%i*1Vz8{)us;n|by%9hJ#{PR>5oXA6$A zES=96JT$(k2`tr(*&(npgiG2iN8;-BrIQZK^6r1I=JscUr_bJe5pj2MahbxHJh#1V z@obJCER!X7e0o`U<8S=i7q(KNyJbMuzVxtWa9eEXvUsCm?b_WJnz-NTAA7ysR6#O) zhry*DhXuTA*D`KBt-zbiKT+eGK&Q~m_H!<`(iz^C%?@}LTewpE__fvzg_+!9c43PW zQkRIjy0&IbcaS;`@vW*AKdjtkG-ybV|k(G70+g+PA-O0tpMTIw$-9a_4a@Lsx zfBx$n-NwF4pxWmH+>aoy6qC& z(6ITvmc#XXj$3lMSI2UJ%F9&;IX` zVgG;OFsmY{jLx{k*YI-3)0dhX7re01-}Yz1_Q$;q2zNiyN@$(4T}fU)rc`9&Hcb;9b=UDsq{ih*j`hd*b^^}t#{eLjyE`B$GoCn=aUwaUy zCb6th*etr@SB}}Xbs2{bg)Yt9$+of|MSg?i{WH$@ComI+^LoQi{mQs(|gr9knOE z#Sag>-*r{5wWFg$BX@}`!_|%k#)Mnu#RWOC;mofYQ9{f4(=v;MjnNmj<_qjhw=}A& z4dO$QZ(e2h_TahyS>d~RKo#MZ31ti++DvQ>o#xM!Dl{4^THdo?GeUL4%{yv8R;{u4 zlOwC_Ng&QbA5*qNi{ANb2`HC`0P$l0rsf;1Snu)F94X0Ve z8L=(7@&Q}-wYqhp*fy!_bnT0~;;vdN%Tz%oo|;g` zFh!AxjUjv<>$Du>S@)l9Ziqrj8(|lJtd6nxvfOS;?3|@#nkHPxNh?7#H14sx$^@55 z(M&>uf)n{XtQk%TF)}w?&Trg(qUYFkN#<-$)2r}fsx@cOv8xXe;}%b*d>!o$EfUvis`G}G+gprAo29FkacIy39vAKv&i zaZ1dS(A^?ZY$(naVtrmIc*E*IIkT{!;KVdfYlc%2jLgUH^S;sAqOieqO_`<%4~o-V zy7UVcW;(MReZKRmUh6qjt0ta0z4^kvR&K|Rjt-C6OJo^Rdm9)Ps+)K>G8$LCezK;& zFZqiTs1X1vQzaO^mzB&DJ|4%Hbj-Zys$Of`0Yu)-)ML80zuvx2&UV(i%d6MDDbk;1 z-f~5+xuc^a+pkVXtoz|bclo1#etzDVa&l5|-P@9sz6Qnx zFLK1$s)Z%qJgqa>l!UO@M40h-lHr}=^R^i|IW|upuMS&#=+@Tk!vBB2FW$IOuv=XJ zSf8x5LFK0>m*=^EjdzNo`0v%5|sqyKThAq@$Y0 zrOoqP6fTr~`h3+w&*FUZw6%ZL?{h6Gv$e7N`oramlDDnRwR?pM7ux>*&R)6_)G)Ox zSnw@{MaS9JX4>b9x3{*=pHb{5BX51J;oZAD4d0E78}{#)FTYo*UiYy(e!-qSf2PF0 ze)*D-!D8;Y$@0a!@2*(8_Go~Qny^qOOVG+IlK(9odz0JG^Hi5{y|q?eb@jmGh5g-K zT!Aq$D!yQ=0|OO{ivC;zw>W?QtrHO!_YCIme`@2W);_Oz(n%H{HQ^~~UbS187(wms zYwP3t=l^-4ZZP{Ss2N=F_v`h?m*@X8$uW~&wR-i%D_2@BT?%42@apR7k8d`ge{?f_ zzU|weVtO$hAzH1x(q ze@dQho`39EuQbS%MSJ$heEs@Wo&7-V?{9NI?UWbqJ+|ON&aAM|`Ck){_wPL9?bhfZ z9klYw$@%*_-MS&oAuW)hf*+LoZHi9mo<6O=UuJdq`eVPozJ7>kR^HfGYaLhdkk#hj zkHk zKTq^px*_wjb=GvFHM{nIdGLZeQKx22$lIV7%BTspMtlR4m)^&P-%i>snRMK|L(Q^) zISkfRSh};EQM&VRx|>m^)AM_dTfG&Pl$4~VZDC&U<~`HD9hbcKWYj6ID3kmm*vN4I zpl8wzXK?e#+R=>@tRBz!l|8Yg{=m$i2Xt*7a z3bZ=kI&J$-ftHSrjwh;b8CD$SkYM0C_kGcGhWi{!%%}$K;MUn-v6jQ^*u+r1xAN;j zu0pg7OjXj`WOd)STm=_XFSaa^W!Ngk$n0>_MtxS?+wk2yhfzwuY>$6C5_RSyMwO#0fq^-S3BoL?GttJMX9WMSME90`%Cwe%J6ZaZ9;i}hLepg?xLb00d>G^L9 za=U=rIWW`UGdK}j=x8&lbSR@vQ3kWLJn`dfw-+c_i{idN0dxP6Fsp>Zop zhG9PXS=M5Q)rvC7Hk7<~t?|!IIi8XN=L}HOZTc4G1#-#>%kT4~@8kWeTe#KRkOd_i z{B9JleQ`ssZF`Hf({gE)@K*QI->~@QS5Oa4NlA%$av8&wWELKVf17PyvYeCwVMBuYGYx64D)!?M=DCYx&?QrTPiux(kq3GB%eI%Vie>8!QsQEfUfjVicmMds zsn1)M2i%(d!o|hKc~Tiek6ZGtz04WkXC<^r)*es#5PIv1p8IyDMzbYBnj+%jy+36h zIxzhXFJJihucJyb=Z(K&7PVnZgdCPEX5QE+64P4d^RG!fb$ZIdBb5Q+Q(h<6;qXndz4WzG_~dzbj_x;u?8w?9vK)LVD(l&4tFdvN#0km=gYJ&|T5?>ZN~s@kf( zu>W>ZJ^QSq`}XC3qg-5EK>c77Czd(3_u4{~IuFME)hi3UJaf_XQuC|H zr*e0$*cJ8jM!5!e_|(Oj$pY>F@66tqJxKx@NI&!v<`#+Fx=?M=bg1}Jz4)n&wy4_g z_iEp6_^q?zXHuGZj#+f0xBTh@tK%)(HnGO-;L?%YYu76mUN8yNeZ6!<>%m;>@<-F6 z^A28F8T@hT`Wja2@^>~rvo}WY^xOT?cy@k%zpD4Ngv-l(z2c%a7o96w^z`k9Td5Lm zk5@aEynS>0(f!aBMYZ=Uq|g1#Q;hrj+il6DRZ|>XT$bEnxN+@TTlV_BVxOO#JuK{R z)A;IDR!sHVtr@wwtPE?SwsK`(Uw8EH_xt@IHJ~Q+Lg)52oB4uD2@EIZCpw4g&S&;% z?|0X~m>~52Z#27J#fp%t<=^^t>u~S$G|Sn!-*>N4hKs~~@TkuuJEnWva-~n5I`yG9 zzlNE?p!8M9o$B{{Eoy!kY}&Nx!vSXg2P>D)JM`!0XGR9=vNs1V``h`^9~;*{%{6h#ndCqE+R%Vs z(%&GiAJ=1kzef0BfrNGWyCd`e|4E-y|L^AvyIQNL=;(_#Zgi|#wQBNw`{HLlZ|?1_ zKJmG)90M%k8$MqVseQ1rE}vgBjwbJyQkAXHsXV(%Tc&Mh?QHqP z-#a@9ZKXHC+ z6!}+lH8AGeZLS|$U-y2r6cp@@+x1NG=sO0xKOYW9Y|FWMvi^I84alu{^W%L}yzUElbB8UDWgiMO61p>&qP=V;4{ z-}7f(0ks@>TrQ?F{Qq&>{*izEFX3%DH=Eq$Dp}U=`=uqO8^zLpym|M#UEJE?>yCVT zd;8-dZv6)zkINtTHNR`}J}H{?N^jZ?h@wMlK425cSfJAsr{ah8-)a{hI&9i3BQaqnNCG-HpgAX$kN@oe?^eC^6dF1Ju zoU}D>M=i5P{FJ3-nulk7FW@w*+bk;DUArnT(0TbIEYjYSZl1e0e$UqXk+#fW#?`fDn#yGdn~&c3KU?mI^06JC z*Rnz@wx!u}lgb#j+|J!_yLYpeQT{gOzwW|^>y>!5CE5A@ZexG`e}AEdo?g)IIfBK{ zW=*uc+BtXsn~>cizwBLITW7seoWx(<;Q}$lWY(cbxg)=J2t{!A=CcOYYz&`#KDh4H z!X%UZ*ZM1FlrwgAXw)u|W#|=SynUzPYNlg-am?Ew6J0sa#Mk?`CZD_eb5>8BPMNTm z-nm29+HYOeYd!Fz@5sh>k11bPH%v`i5_&*AqBgDMJ@4V@Q>_ZgS9dIw7wz6t@Z+pT zlUexvv;U6A?D=!EN`)07z|>8jW+ zfz_8=Ho4xs@wcN)Y=eN{#BCnd3{tvGY<>DHoN|dQHetuFeOJn#CVlVHWAklWW%l~S zO=|rmUOgejdWF|%TbH5){)rlH2iPXqpY6>spZ`~?J4EQl$IE^(-%_5Q?9a5U`Zz6W zVtF4s`@A>DBIBR#Ysk7-#6M96)H@U0z-UpUW>K8WvAxiG#>K2ETaN!G?0Zr#{qr%- z+Pr^m|F87wiC(K>IRyn@_HJQbAgG+cu;F>sg<^$8Ws5&n_Evc7 zZg?JD_-5*YTT3OTU90Z60vdO}!*NT&Cbi6D|LInR#ez2k9WTkaUyS*#DS!Pi+nRLU zvq$&E7aDd%eYxPgMDLqG=OecZY=H*H`}Q3^w)?Z@qX)vaGya{m*?QWed(zIl11u$P zp7FD~xGa&$WOrEXc3|h_J03@W@;vZpW#%k-C*^Q-Y2^GLrh8rVh3%go-4`00H9gUz z>ZM{3C|MMWZgAXDbjD+*FUMB32b$Yf#P)s3`hTH$)%lczHM9J#oqn~u)BIV||MmSz z3$z*2^{SO7f4StmL=NOAt%St}jAreed)s#iq*$y}x$q-P;HvkRy}n&X`JT*uf9&4z zA36$xf||~^7%!|4G+_9VJ?DjW)|vaiGfeV1ly=|Z4FC7?eb-CPs+n7L_<7f}?^N|^ zd%3dRLl)F!4L!g%$^2^9G~dLnrzhN-=U!#-)Z2N#d8Gg9_nWlVdL>=e=Pp3?c?V!bgMM?A8G%W9LV#{i2a<%-me0wa@C+QaWjeT z!#(e6dm8=Kr?Skt&z1Ox*ZKJVhJTC1H;Zb-zB$^cUUfz9sKJ+O&QmHuuIIRsP_)eY zVzI@lGRdH<8v-p?vn7i@w$Jn3rzu^kZ_K8qq-1!BufeUafl*=kt)_3eHNOwPl=&Gf zzv}#kACFmV-Sqc;G_xpPlX3f8Zj+H?zpcwm&|pVT1Me~OldT>K>&hfANxtw}@I6E1 z@j8h*hAY`$jNWmXf3v%{vP`mV!F9PsXA&W)e{%3!hc)(nD(MG;W#TS3FXFzly!LZ@ zU`VYR`@>@_)y*9p$9!%vUU(t6;U?FO%{edZOjeX>265&{9GJ3S{^o4^FH(CgYSJpY z-mHkezPwCx*K(NI`?zm3$heBJ@bfF=uX27HwB$muMw>^v%)icWg0=hR-9H#vn{Ad= zf47%;+9>Aj@{>b^=@S!C(Vbg1qd{{ z-dmv11s?0^=V@k)6qFWSoLI(?!csHG_IA&` z!^ip>gFP2~|G;o8-h9{d|2+bGpUi)J=;898SvyJ$o%jQ*K(i|5O=cSLTZL=gdu%zj zzwbQ|p4zURU$VZ=_z>Hgso`l0GnwZz`-m%l5m+i--QiNGoY8!0W`^VI-^cbbSG+hW z=kudv+O1R%&GIhUJ{jvOyC>C)J3BI-ERki9WR!nC`_Z9<<%RPZeN-Fg{S6eE^IGof zbt%=VkJD~n{1CPOR_cpMfmI!s+HGAHNia%>e5;V2;Gt-@piK2uMx8#3`#Mp+`!Axe zYFs>hW&O{KzpU5~lo)vN2mad7?jfR_z+C(!@?!Ok8NvzY*E!DlTAEi}q;@mZjHk`_ z#oEWR#Y#$2l9}udZU^dedy<^~~bq!_9lMrZ2R1 zd!^{b4=Oa=F0h4JTyN{EEexr3d+TMiIyGA&tevy${f#<5!%J*yOwE69eCz(l+bt7W z%QMOM22bJ|de^rEY?`pT-u;mM&9(}G7q+K%KXY+0v4OZ^viZ^8i1{DyeY>z*z3F7U z%-8k?j?!ZIiYuXR1`j0To<~%p&V(z~` zfo{>Y=9QiDXY*UG=sAb?TytJx19GxfLbI`@?+e?MITAk}$X5lHG`hbvGEuFcCA=|! zsd=v_c);(Zr!|8Slik$aH~IoL7bZP;aFEgG$Cn$A;}29^>^{l##IO9}jrg7F3xZ(A zCb8rc94*e}s7dHqoAB#ZGRJqd388x1JnnW}^gb$@XR_}UV}PdH8^s_oNC~v#&dfZO z-@SrwWh)o*@4WmbYx+V@Rq?ci`!%~*;z9Eo=?ZTdPH_APSy!X&^wv;?cLLM)y`l$? zoKNGPf4t*j_ffuz(>WRUPwj5-GDzZ|xCxXMr5&Z`T+M8Umf<-#zr)=oQ-<(dj!2zvwH9od3Kkw(FVX)^-m|P)_Mxz^N&Gq1@yCdzQR=7ufx0 za;Ih1n@o(qxzw)Y>Z{%U9UUG$TbLD;ANXs(t!wec_1ucqw3Vz&}w&;}oQ@7*WHrz@T5S+MV$=})2 z1V5-(dh2)`Ynb8pKW=@sS;wW;$ga8PPh9l=r=u>!{2Dw=fbnn>;a{7i_ZUNjH!M8Wp?6u-1Ul{^Jb7Uz=3?( z^Ba@4)7j<^BCnONWk=lapL6rn{lCdym_+7*dSZGhL6f|LQot z_@Kz6TFxiu3$*IqIV|8^8FK1Fb-~oSQ;mnuw*(d~G+KYdc+$hF74IW>VyoE&c#lq$ zVpwB$Ao``tzev4^9g77U@3(uz?0Tl?%6{P3F}FW|p4%VKUca}k@bR&Seb(;|goekm zGIShGirAVZTKE6={l$wGF`b^Se|%5n=QbO0LoEy5sMhY`Tk>?T!kwjMGcFbA z&$@WJP?fvpM}hC=Td5hBURJHxuh}5(Q6Y>hY8u zdvI#|qF-t9TH?P$R?JB`_}*CsJn6xl5FanUZs#+pn8Kr?6;CF*Z%98c7gPK7YQ^id z+YRdf)gX@q96RQw6SqeK)Ioo8a`IHIg9e|L?GTiI5wM%b)nu<>L)Yobbgj8Dm10XK zseO^{J{W#-`|F<2CAN*<-|;T^ewZ;NZenlb{!{JSzkIpxBpmqQhW6eY|5wMo`IE_? z_9vO|Sp32_X6`O7N+P@8Z_2%GcIEo@$L#h$81?piaJp0X`|Xae*P=h3(q6A1!u7C` zUG6~b_q*+^++q)&&CWk|X{qVrwZO^ZcYAQy!O(S z#Zy6j=ilG2Cw;Fkf4J2D`oS={qhDXEC!DjbCN>o4^;l2*MFD4 z{LxQoeq}K~&yT+!6+Y|Zmsw%ElPjXFszT@a$6H&o?G2yLv9K#!!+7`ZsrUa^-(avP zeHF6j<1y(1tA6{_qMa@UPfiF?$~e1H~A(@d8exGbNi2cevx(iv*h~sA8Wn0H!$G!v{%WaPzYUKt^WhEYFTYCC(>stP@1F>p?GsBJDA zA3RLj>V0u>zmG<^Se{~{Vr&cV&1NoYJWaW-`{fndGq=5KZ-979y~aC zzI3OH-%jNpg?Bi9ohe!#d7)Op>W{kDnwMuej_G>ck85jWFi)^O_{{9!E7Qfr+v~RR zdo6G5V~P74#LOi&)9{_y(>gajtK+TO^W@rN`AZjlJi8$2xoYLd-vXVJ>Xcs;SoB2a z?PPsEw_MIF|6Wbe3+9GPmx88DpT2nCKDnh!mohZKm%sI|Za9C&`d(+_)?0C$GMSx= zE9!p9^*n#IJNT(z`NGm-iJzON+%IzXUsUV4r>t7}+^)O*tE8vw_nsb_cr^5U@8Q+) z*-lGz8&8`565Y>sviUxwgKM61qv7mq^T+3`-ydO@uPONS_~%dXSl7E=gQb zw_JYKtXUs^Jnnyd*8G0Rt50o{jZd~7jC*@e@vW6W>&>jD%ja#3g1huf9-R39Y2D?S zDvuBRn%%eXkCuN_+kWAgpJFk;?k|n_-^#Y*ldi`87VVZ!jm;v`4E!<{4d>4J#Z*3> zy1eGkzTfW-e|dR%W68^)8E4aUqPOu(ojNsQa~f|}R#uyhym^zI@oec0);cfp{SHKD zN&e$jx%?^5p6l7N@1bR&DEPZ#&r7{G|JLiNJhs?5q0#T_t?>SZFOvAr{cC4pc==Y* zzf-ikF1l^MUgm*E_pKkyKQVXl*N*$`9wED)E!?8R6v!mr$diBgQreXoO=;E7Tfa+f zy`GnSWJ>g=z%s-4H7S?k|L+qt-aobb(ayZB-V@K={(q_K-{QcSpT|%2#BG1K)Gp_L zp!&LKd$}<2``3M3mR~xeBydW{f1>DJ_O06mRHm{hHL7quD`0-%qdn`f+ws1RTd5w) zOTAg{?pE4SQ<~Op&md91m0xvw<%d_{(nfW;N=AB_><;UgZ=TKms5Kf7+pP#3bqwsHHut4%-u2mb2JunD#;IryDxq8CGT-P%8C z>zC-&$twk3TF(?F70$?iKK2d6qg%U~Z+!3)`1yRv|9A4?_iKavdYL;tggP}69xA6Z zH1NjfyPREXd1B_1A4k@wZfY0%yG-VMY@x{0=~X{$7W94i4V~+db&!g?S7LrB_m(t+ z^)~(mFJ$&j4&zPJo?v>jZ{6No_nbU5tL0v^oKZWo$N4Ys)whvpbN>481m(l5=}sw>9zHVNCf_RCTXNx^K@7`#-{J6Q@+PKWuOCRLFm$TZJ$TGYwlYUcnPt{6UIN@{j#qyMbALpmK z9G}N{qgE`!c!T)I&T_`fD^KZtRjBE)kJ-C(598JP>FnhXFZdr?anX6YO;}^{`R_`X z=AC#^s%iV97B;@LWMbitvoCI~o_oD|>mEs`Z&SX!`|qAU&HQR-K(I*R+bzA1`qsX7 z-g3R?<^2iob!E)9b57njrMdU?Y|wyWiO(&@34$L^{78?wSnOe4@3*&o(TY7f#>@(f zr(ZUXX3Ss8Y|GxeQ_5m0%eMb(q@Vxak#OUD=<4=m_op>S$Y1!9rt$c?S}LftNMi9> zJ~i~I!;VmC<^(NH7pp}YGbYdMT`jqnXYCIq@mXe4#fm5T+N~$rpKUxe^YDL#MQQx8 zJ2nYK@4MbNrA+SQ(#sFND++;{GTP~rq&F@wV>@uJw%K5&&z(9Krq<%SFYYa$o9;eU zNZ#Y&ja_zg{wx3d?{sedvzGrl*57O|E-lk^S>gdt+N#_I&=&q_JCvJM8%hGsW=T?Rq~F4f)o-bGA|X|NnkT zRn~WbOnK!5sY$xGo9;P(vOj-yqiIvuX`}CaQEvLqn{P1AKHs%2>hT(um$pLA@BeX} z`?q|LV14f`+b3sC7e7hT-}a|^c9!V`iLKQgTlxhJcuxi=xxW|ung8@iRZz|U^)9FM zzAkvG`&h^5#;5m;A$K3NADwqETtMsXo83i%H$)swPP8+2dRRidqqqB5)4Qmh|JUr| zbD3zQ{Pz9-zZ!1(>Cfh8IsXnzwqcjxRh!tOV)%FV;_AokJhmy{X9bolWWQZ^WnGr6 ziJ8$CflTT9tUk`4e2W%6T^qOcw8i3WYyW=mJ@Muh|GwtWXS$wMY_z}F8zN-j`FrsK zqq2YQ>S=$A|E;K>zv;^b%kR$~F3s*&l=3+=Z|Ud%5mw6H{1eMQ#vU+PGUue<3?EO2 zN0mEEgISh}r#|>9Z;}+IKEY?^g<6k)$Jt$2`j1}x;3dGCp8oD<@Y0WZ3WA!xki@pL z=xj&ZqzMFdJ*50{t;nVrANgKTP3x2r9ZXIZJ-+0;XdZSr?mg<<~zHv8NSElN+q!bb{$7Ayx zPm5WGm#~!`2`I@+yZ68PRH@`%p3^VQ53_LpivD6A=+f+2@`fYv*{6x`e}6CB^mW4W z8(jZc)cMxG4BIUtrFgTdV@tcB!Q8-eEY?bIi~gzJi?(h!#BKEX-~IliFGfb^inMRt zzka{E@%x`u+&?$Z*-(A_U#{OVxmiE?D_7Jma(`=-(US5;@sgO^ffBcSTI_79G7O6k ziY(g8`6TQ!>w!(Hb{uK_7F%@BR&I5OP{yqb{hQzG1?(0%`qXFJ;fLE7KK#1i)>46q zZQyRSoNWA3-jj_}O3t%w&E-EByXf49<1bQEip~`+4o}x^d_C`dqwQ^;^$#{L-}LRm zmFUeI-|v@X&Td>|>vHQL$Bc6)QrOP$3IFOl`PHHR`#cSmTJ6eax> zrKVIUpA5Y9-^D(yuJ{coPO|w!bxt{%)YP=vp7gyYr}Q>~#pm#;SsAV#n;18g21~FO z%VbS=JX|1p^TN((`P&Ce-pI0?;hb{j{Ga0GGeb)EitmC%Ne>Pq)3=*k=A^qrj4H<+c*_Tm2^%pS-{a9v^P(1sJo~m-zzB$c@_owstxGYu&4I@;1`6C=F@n+)l z^Ew~mSW?AzUfpfCc)H+bE5;dTUff#$_HDweYb^$kXaDS)>;Fye)cj`-r&1i-&u#m7 zQTS$c$CK&nDm`x9HZI%EyK>4m<&{30b~3hKZVFJ~^q=@^^Ir#r#yy*_N;7nKafR;Q zCie7f$)e?D0+pAy?95o6dtD*uB+I9Cjqv?qRT+EFLQj4N4#4rChFYnFTXLUD(!@%F6sXv?7wDT{a>7YU$SX)z}<)a7f-)P3AesjAM*E}WVegRi#x&x zcW6(I{9D~5R&`=(c4Ph=X8Cn*gLjKq@#=mNIN2*`(Ccxn!gv?|s(S4O-y*;dI#V<4YvcGGDC+CNhF{l+r zy;PheuACq>am`0Rm#8cExt0}|v$Q^6x#d51Z!c#-|MrE|^S>oMaQ-iS?8<)gcOQa{ zt~ZB>aIgzyO?PsUQUu3S&YsjVm#C|Mi@lFE7DS6!PQD_`Ai1|E&aapGt@YLaj7cU+ zZ#h-iCTq-mANuD`)T@PuKHFG)SuVHr^n(72LHrZv><~2R?bv>#DTS+J^$xxzuMa4! z$@^;7VAm+Up)c+pkGfU)%KaS~Cg1EXHt%_9HsjKVyFT8wCNH_|O!l37&XBrzshpB! zp>l?l;qv=!^PB~*uDie2zGu$Lva9R-lfw;8_eSIg1xq_^&*xnHY;|SU^o3UUE2Ph; zaa(;a@^?^@ng<>w+E^?WVQ^#JJ*NrR%Gm-e-2XoRFR{00-c`;Uo6fxW86a&~aQ&ag zq>U;?YD=H9^R4|FaQ3X-m4AKtY!f@LRdqZ$9D2Yc>0XibufC&!%PxIYKFj&!Oh7rS z!GWvg$y04l_Rae%_{3-0f4z76{=am(Y8shoz5iO(ijqLFtoqyCDKZA$zwh3S*mBKT zMSdOAGQrcUOcvHIw-Wr2kk9q8g6ZbtpMRExiOzlPaJww`NQW&uTVQ2xVfg>LImXAC zytlFH&Taxl$E3-iVPTuHeJA~v9C~%fHut`8ld*drC~4i;p)K>`miB+A3CFfO7o^R; z{A-(}?zamqukWn*D=Xf7!l~QVr6`fbNBNm7R&dN!Z$MSu%7o75M|9q)O%nn=r|JZnV`b>Ly-ERUX zRbxL-*ut#PJNGQ>jlE$lw;QQOg_lTowR|6};$OkX2DlRaWbEwG{&0Q8eg+l` zxxGt*G~XRGWqYS29d(b(+nT{_zkJ)_>_+YLldiqa*xJ4-@56(`_UydY49g}RSQlTm z_}$LD9bRu{lr&4TivOArTJ|cw($+nRV~^~idGpH?Z}VFPeioHqR-@^DF>1>#XO;ev zvygVzQH#%OILzG}EoVwIH1`%R$w;Ph9rjY_FxCkl0ZgEb;#C^|m#IVKKk#r=D0P@oT;1e-`OK!hcrmJlemmx%tMY zzh*L(Rvn+^b)4>nq{^jU%;KMDvx94f+fnm_P48}l?T4Gw)brS#)GSL-_s94^94~Zy}ueB@gX|uLZTKPl#LK%N(+(Y+z zXBp1iJy%lCEUNXiss9`J@$Ek^??=BXQtFWI7Pb-#6<^QXJ5?`==~C|AkJ6D~TpX<@d1 z#M_`tYS-o;>Arq_@h30&XaBMrA3fjy#eC_Ch&Kn1SxbxBn@_w29#ic$oXOaIyd@zy zP3^$E!+|q$)gn}*OisOU?Wyl#mbHFzOTcPFGDFyOQ{|8P3;xQjig6GySJm3@0xSCk4SHv-zR?yCoW^z_;VwhW z(;{81^e;+k%WriyemkS5B6aTK^6J1E&u8z9viG)o>3vkZG)JJzAa!e=r!3dG(noKs zi!ax$W!H=o}Ew@ zapAtsd4U@{%NatiZ#wUG>Y*_+--6uL`;s2`e&2d1$td|A^GxYAdX-ybVs{_8|9{FI?yBzSH!ez39k+yX=b_uC~14`8HFo zYDLmSarf(ASI4#`*~MIL{!3FbLP$L+%FfLAFuzo zFnztue_mJDTf#Sf7MP2z(F#x5`aWalC*2(%jvuXfoPE}Vr!OZt4rW1ASgY_uhhOg* z%H(BQ&L$h)k!J5X5dG9)y-M`KEqld}Em6IhHGQGwy&&yLfeTiL{r-J$PhmFmvU~H` z|GRwMcXLY21GkGd;`es!jod2LF$LUHj+mo;on`akmW)}K=g<3_SaJNM`Cr2-!RPCX zl*EGSSsj-sHy*raQMB6lc6!yzg_m08O!h8&>${NO<|eOqpJecVzrWw5<3BVjxi}YI zSu6d|R~tNAuACqx+lt!6K_LL+>xE1&U_GwxLEFUSN_RNfKg^uDIA zpb2%i114%=(0OeFrzk?8wLB0jfN-ZGrvh>%nj-zo6YFTQ1Z z2UGJxjUqM6OJ%(qlfEqKd85UBYuoDY|C5xsT&f>w?A7>uasIq2jp|k0yz8^y>m|H> zrMF~F+LfBfeb4Qm*RAKMEiE*e|M1Q&ap)m2NG*7uV1`yJ^F&+*LY8mxB4H>K82dn3`+LF44ba!}pDHXQCt5lxZ$c zHDWH_H@!Xg{(Fg{&AGN;1y(u-3w{8t75a2qfBEApCSLEjtc50ixy$(Go^0FlrsiGq zo`zTpW(B_IOL%;mc}9h>dA(TZA5Hm13t1+8dh^3tl*h{}Fytee6d|7i*ja zm(Pvgul2TiD@utTWj_$0(b6toCsF>T;TLqcFn;*Nl&9din z_=(Ez?-%P;t>9q~|30V8u>%xh3p+V}fJWJu`OIumpI`Gyc>lJnt6FPfcZ+?$UvIx| z@3$z>n14H;tkkENo`$K3#mviGxF(*Pm>Xpk)vO!%WBJLS1$K+YpH^tSd0{)@=k?A~ z`JO{b&TjYreGII6z3@_H-j)3;dbc<2`K<4|VxMLBs#baDk5eK73M<97+5R}`=$~)v zGIMz+$B)nR|J%%~erH+pu6Vxn)5i+~IQVToFudFS-VU^P+WdZv@ie{IuI2OVs!kN+ zS*llUXx`s$w=Fd~aQog$pZMx!kEbm&b8#;Gp>7fK-}aNm-s8(23Cy}DcjiawciYP? zYd-e%#D}jxK98IIV!Zt1W%CyHWyj|qbMO{Qzx!%Qvn6d|9sS2Q-71C7hdMsiC7|^YroBwQDnSa$G_~ zub%8bJI7LZclmp{PcbV)xYQ;~t_*qA^S|cH=gOL2$@#6a6OY}PvG1GD@q>#EZ`=us zc=PJNoPW6thfIH5%1Zas$NzbjxmZ7v+@JUS8q2Zx=O=%Z`pd90<-U0HM(XnQQx$Kt zOq6GuEZlFY=W6RL`~H8Vgz>(q+>;~ws@#flwlO{p|M}ryzm@;bitzC7CpjY8tSWV$ z7aCc9nUic5J;6iczs7XUuNS#QHkG^#@|j~H`2WxO|2zygw&hBnv;7`ZQc|Lzug}l$ z;f(S52krKM6u0Hyue((Y+8PbnJRGpg=J3t$Dds|dInVq^+qdZX-w*Y>FKkk6`)QxJ zXJPq9gAdpBuRNUf{CN7CS)Cdi{(W}e_=Ep(V)1Et=B%UD_j8@^l~}9QUtC|Nx%}!T z)uy%ga}SrkKE5*c)C$$d2jazQJ6u;ES6uCMXsw{fEw2dj${MW`(*V`ZF`z+bi!_%;I=~4}SefQehzbE_O z-rmk%{OpY6r|cbE5tjS1*oFR{+RvzPN67u0iDftgvxOw%WTQ293%)x%V{2c0zBRCB z-Y)K*FE^Z7H`{W^$xbP~NUh~Im%jK7=`xOP5H=X(+6KeJ3Y}vwx zUk$j<7p{&yFJF%KR*f1RhH^i83eMP3gJf&2nt@D@@cX&PdciK`Q!&z(UBmetc zC+9vobne-=o|x4izuG-HAM$@W_ug*-FOSOQ{Qe%c^|Zy&YwWiTUUq-l`@URfji2zb zdMozb|5gh>p6)c=)@A18b(J0Onans=U+(e94%ztmYt0Ooi#tEvv9*wO(?3=8!)Jr# zT8^z+|7~6#-@4Dfq09AEzQMZda;MgAR2K6}++bdt&neaW@0bN!RPu~Jm6qXZ{YSUl z{+m&`UA}T3Q|vc^)IK${Cbi0^=Qk}}@3(XR4(sjh6Kw8GviOtsMb7WN(}WZm z4%u$y{As(+|9?5_<&WF3QTL}l-db?G?owZu_Fse3vRXY|ceCB@>6r7`++19y>9Sn) zEyIb>1E)Ssxh24}{AG$Eo-|$8whUdE1bZM6o9Ksnt3u zYrfgt>z(t`@5bIs|4L2LB1OY?e4g8p6a2q2?2_Ku^-OA|mELbVw6=&R+;3<85N7z# zRO!En_@Vu;r2YxBmb{ZV@TU6Hr7B*gJFH?AU2k95iMXCiDg9UPr`mX|)KRBmjnCBA z>4zfY&mX)OUzU8R^jE}#A}5_GPOqocxiEQtF_-(lTh(%gQ`l5V*3GsjR&G4JWdF># z@FlVgJGmkv#iSpc*rs%6UdwKWH~&vjRw{ z0`eWo%ia2{)vDerW=~_iTdlN;w}E&4w}aENoAT$>U6%X3#o$*p+l5$ZSHIqG^#ysi z#SX5EFHsQkI24@Ln%^Y6GhOz-wW^rm!t{C>UVV0f%>Km>Wcbr+t@+p0WzT5LBb-+tb^>`EGQ7XQ(@ zyR5vaK=kI)%D=PkRj+tuEC=e@an28X{_aYn`fr2B(!ciIjt&Xw(Dc^1@!v)KnQg?H zGRx#^JN3K1@4L-2`DBl#_a2GzT0Ncg8R>@`S<39Mb(Wf6Tb!{@TKrBtOLuJF|Mgqd zVm}MyUy$lpeWdT8-$!=O-~C4(%vLxqD|N~yeP?Rp|6JZj0` z@Zm9I{(6Q7dYc~n`OQ`McbC?Um^7PuSKC?lrzBsslm(6Hg&sI{dY#eY=M!#lJFr&= zZ`k4ULF8n;cl>#6i>o0Px7()qzO3wx_|2Xua(;R`!9l&fI%b?=BDi zy<09s=iIfw#Yv40(pFWs9?VX0ej7Zqa&lkLqlqHd?iTPod}gt(dMDT9e5t?Lzt-pp zpZo2o_{=t>@8Fuzvh&m07rz76bXN889e6eI+sVe@7nQ8bXYCQ5 zrTFRG_R|6@`lPS_-^$hKo4DuaTCvtuCTw{uzK8CJ zn(>r$@HV`acbTyCG5fa6ihD-VY)zahNAuL`oqk`d`=R#nVv~<{+zXfI-wyV~eZH|v zE+zEVgEJHDJp;TBmaqqTa$A+9emh^K^Go`1#mmj^hF`zsoQYX=$9PuQbk1=5=+4G@ z5_T0Uf;WD?bi)I*AS`U&`b&)|F>lWKemKn1b$;@5zOyr5*zODNzM9#z{6gEi=2OSo zN?$8$*u}5A`qb^k!b3OZj<7xM;{TGq>HkRuwR*>YttnnVGWh3~?&fd*AO23l;?leW zzx1yqevAFNB68B~Y1uLRb_vN^J5~IBeKyl_|2+dHwL~%N1hbg}Ts^fMZ)}Yi&R68# zS!(E1vs`TT^{#(bsp3s6uh|1i7RIcvH+JH`c}K0HYt4(`AN!vNIO!I;>=xG8?ef3d zz9G1J#UCU2TiabH1l0e3d%4R!D5Czdq0iqaJGEM~z4>W7Yf4oY{pI}g|Jyqon{%7) zN&hJSW;^*zO7^TY%_Y9kwW7jQ5=%3fK; zw}}1TXC|-ywOvz!{LB}K>{=FmN$0H9WS_LNYo7PdS>`)Cx5(c3vrXl``LAD<+WJhM zIs41>^mjY=f8X>p-rD$=)nrEHgh?&uRyjX-al|c8D9~X!>!;<~$t)9$mM1LZo45B; ztrBmk+Jr8~b6RyTPP4jCm2{6geeT4Sf75GLM5N9B(rw?Cyd`_0`Pq)euYVWaoszhA z|4)rmdK;6!oq8j4di$I5%NI6kcI@3Bc<#=LlP=r#H}y0a9}ZdJ{^G`wXA3I6wlBe_1+hEB+uP?t;MUHFv%sZS8m>~ zPTy}I&D|E1pL*|GrIOY0zG&$S9c>}Omtjoi`!Dxz;(j9fpXt@gg^w!3j?{eYai1(} zt!%LC%j92A`h(|2IWOg2^Wwdy;*EJnJ*OvI>{`cvH`;l}m+em@>vO#q^3AW?BJ%$3 z&Ay6=nZH%iq(FAoEiF7WdtzCgPgUAs#mPwyz#<+<)1oimG7zC`#+K0yWUjvmQ3wF{oNPZR!)!n z`0aXB(Pf6IlXBT72o|J8EPEyR-S^Gp^c9v|*K;b@UXtmJdbB&n(5m!%7>Rv)ONa@{`Yb$~NOu%I|HH1rJ}HP`r!LBlBXl!uHdR<^NAH zXPC|?ul}W8lX9r_{M6+8-Ev>PKYNhlB$>?OGsz%=aV;Z9^`-8{=C_ih9 zr`YG}!q?6p_cA=S`tP61l^#98N7Z%EZJTAB%sgg1A%asskw$*IYl$fhtEtO||^IcUivwPtirqZMP z51cR!`Ms3s%iH(vTJFuUv)|_(o@u@-L3i5q7r*z^>p3W~nd~{K@x#LHfQn@CM7f9G zzI2@|S7G~JlbXX&J7?CT-7y7Wk_iTxEn$6bRg*Z%jy+sA_4%Ektyd?kz0mM%Kf6l( zy=7bX7RWT&PV_GOKS<|emtK)s#n-v_nmnK8g6RK#rC>i z+I^CJ)}j-$7p@U>Fg;qnGii2~!~Kb!?rn4BT>`&J{`Aucjpm(Ncd9Q%NAcwe?;X~S z9p-J9K3^&Pa(kuQ>V+@=2Qd79m=Sh9@!roZ?3Gikrp!^kf8UzH;{C*=E%{uEj~*H* zC?`znxww{H=y5#Dqz>JF6`#+ayQe=Fvb?kUS+^k|n)hW- z!?iVz*WEA8Jyr1IYn;V71I=eFiYt4=oR`_?_T(|&GW)&x(zFPJ&SUn!HvIUg8Chv7 zShbn$a(+qDT}trx0Qx@OpHtCpFTqh9K|*l3#||C(j3{$=mn`7DW=lfIvrAeFu1 zSC+NZJ&s3bL?6xFR#46qSTCeHe^U9K8J@{k1D;-GwfL}mvB$BA7bjaBisdZXCx7JX z>($2Jm4q&5H17N#)bRHH`vaT)&zbi6S&zmIx#Rn$@A*(C;;J;M;o`#9od+6L*f-Y2 z%Sf#i+omO)wzlMut@dC3Ju_eNpWMcrVIp;d@r*+7KE}pt95Ken-`&`n->1Up4)?WB#?2H>O4le`%?I+xNiv`#l4$^wZh9SL|DC znDaYPUFq##4a@tDo8N1@Je2vIDO(wz%(J$hKk4<~+zTpOW|WsRp1S;F`?-*MkE8pi zy-L%)u6)O8-wVAvQmw1QB0vAv@HS~UxYT`t<=<9?dwT_a_6F=#b}mXicin0A^@HjM z<@#l3AD&Uk{P_OH5dH0oey*Q+n(MD`VM(pPfj#d<(%zUFt$f$$l(X00u6tvU=YwYX z^2cwdo}E#(kTrbm?7r63T>&=#znqgd*pPbd&H5*gS6yFi@I9XE!^OvOf1b+E7I=B` zEtB}oGKQU87Mpv@{_fb`)Y&@e#C+xT{&9~NZWEoeGpXps#o5oE9Xx3vz<>O>Dbrl< zthI{iy@K=lkBT*Y`8|nat9xL`@^^m!Dza7-Roc~j4O}8LGxPB8{;~&WL|G+!XBe-{ znJ607EB5wwk^RDT`^`JIhMOJ;mp}T+K0NtbcueNsnB}og4qUBY=W?==Z%<_?bJVfI ze`eu~lb3v2Wy(?BU%Q>NdXmSHNd`|fwapl=Y;Q6S_d0KUpj>0`?BMV3I7-ykzABJmUF_vS2-s;)+gpZVfH<)#ACd# ztbgS$I;FpD;nUcc+wNPwYd*fq;H%Q?N7j#m&u>?}GSTTe-^+Tp#}|M33C!MiyHC~M z|Dfb`FXOcbBXf^jeZA_!9=(RM{Au=?6Cl;>)6f}$-+N^ltG~M?%w5!@%j~!G+7G?k zZ-RL{zNuZlb!J<2I!mzrgnIUXwWcZhmz|ei^r`3z(wX73z_M#n>l=o??QN^p7cR0#QI&-LiAeEH*3*R@sRFXnB?+*p&Vy6b*_ z&kfEL7N1F*M6)M$^w;)fzu#xkC9d)LzT?`DvD-4gHe6fa_~6Em{D8l*oFdn58r{&! zT(PrFCz$ueP3w&wpKi3--8g-<^4t5B$D*CT8QJQ%vVB=*IXB5f>*Zz?` z`Q`N|!Hlkj-wVF6wOeg#+wgK8U(AQkK{uAAq)CLXW}mh{vU6Q-;+u*?)wk{&GpU_? z*&|`h)>q<|xPIDy58kKKmR@1_zwqN(N$xF27jDf;Joo&q`#G6VHksQ-8Inp z$?IqNho@@)d4#3zfARbj)A{0n&6&#bU8%$s=TR}rI){1@d%XD>{8`G0B;bN%k6lIizjJ1(`} zj6Au{JiFkhx9E=_>9HT)*ne!CC1~JzF5$|->lU9l7OXd9>kZ`IYxKK*Rp0wfXJeMQ zKkxh4E4DSadtpkWkz4zf80UA9X)ZRMhRZwWtg3bAs9Mdo%;x64y(ytqrcYLxI<3BV z=(@tcwEf>*?Vr_V7oGWQJ@eig$Cs{UZ~kwe*R$_j^mdblUu>k$?XlbaJHkqf$Kw4% zp`1Nyc|U$Q@a%j~gJ%)zXC96%eSb7hl^ZZ^b~@Lyf0oi{CvyMQyH}*&F`G&jN<<{ zTifivF5UV;M|O8u^~Lwz5C0r5`us5eRw|#*-ZNd>_g|ZI@%lt?4tW;I#rwg$_s_?D z>GeE)nWyeDxt;wOyX~p(Wo0?9zRC}kSLQpdG3QL5YUK9)YMgUl7ITgJrP=H^_FDIB z;+|r7Xg2T5YXLV?!%A|*Zof)-zrHG<`t_5yEz#Yj1%BoW3uEU$IWKE7?ZCqORfR9} zPpj6yE77Uiz!DhGR9$!XG{=l3>DuNB59KRaHoOlIX-(Jse&6Bj$5_1?ci6tjN&M|y z{duy2j_WDOsR7lzzLw`(rfolFxAMV`d(MUmO9P+Wd^mGb9`97$EerRnw?z~rrSJbW z`%%ZYBY|IUo~&iLKgT-u#`ay|pPtvgzLs|8f1IJtyF1I{#sA(fdK3S+(K9)b#phDd z6lp=;o;!>^CC4Y4J1-Pw(0u)AnY3VzT(A<*oDX}qV$-c1&u$CkN&i3LeEI%g z4cqw2<8G%q!}udk@nqDNCmDWN7hhrcG>P40@}IkVpX6{)x+t>ex7xO#S_}DqJGieU z*)EIM?Y|g&Tf{bNe++|jZKm72_-(3(I4nMwxE;v2sC1$1SK9i{7Y8EOR*IgITx&V? z*5lndr&nJ)7;XP+iOlHdY{0#U$Zv#*GjLy6)p4Q$@IIEY{C~<_AT3L7xTeqrpkHssh07@ zic@_Qud4lhKQ+g?p{HT9^ZN^mH)?*VI&)J&Ra_mOof(Tz2H5|1w}FZ}e|e(9P`Q>6uSY->60 z)>|-bp8eGE@$~vO$%_w8xVL|m>2ONFcpziG%yOX>t_O3JZ*A+@EADfH*vH9f$Kfa_8LHvX0+yzF6_TS^9ye zRnj_?&3u0;e@OYwcb>z-nnzF9>F)CSs#lf*?gx2N>nDX**G$i?GHy6>Lc7t{nrVYy z@U<0o*&G_S*V*nb+~~e|*0olHuSrjP-ld&($v$Cwe9kqk>-&yLBrRlFrn7PR-zYD$ z&)TdOe_kC4XI?((K)e0Nq>KICDmhcP?BJSVq_jV%eP5!P>{t2TgQosB<)2F4Pn~{l zMdPsrt>xc8*dF*MTj&$HVvVco8byWu$65-G-;j1{yM_1RfP=cAHZ-e35z`jeBapvc_S`^+_#UY2C^-dJ~qccE?dB9Ba=O*gEP zHx-;UQk3*kx@7!mOV{Hd&jr613#8S%%68lPnZ6BT|F-9W@9*D92UVp0Uupj=clwN0 z^v=oP!lcaNlrBdp|KBGf{R|sDm?|$%v1K~(L+>|N+Utv(|9_c!V55M*cU95t_A=QS z4aIGe_iUB(w_p5IZXZ+p#JXvcxA0oUUa<|KvoCyoc;mN$+5WffhW_?1&;6BMe@ox; z;%4peSlX-HtsBHQQA%t7jh^K9y?dE{A2@9{ z=O^peWAC3a$@&~++pTqS(Xy>yI>Px-b`toUGuR8Z=d}Y{B!% z7tbg4{hG`n-@W*g*UWVZXD3DepVH2(-}^sh!dv@SN4A{Ky&!pYw{{~WW4gb-pj7ha z$l1G_jrqdYe#(E=^M2F$1lAx=UZ=Uv&vh%meOVmhDD5DBTuLfo{v!7+`mGzDfBBhl z-0XLvrdh%7>4F_sCLCSMToP6Bt#{9XhZp|8Zv5+`TypB}=6IoC{d3L_B6$8DK9gx5 zxtgwjFQF@a>pK%(wZ7_=Kduzzd@f2k>b~61h#k*L70>+Iw9EZ}L!J8=dLV#|DsdyKtM82pVm5tWu1i~5*MQD z7cooZ_$w^+Kltd;Tep)r!MquFORCs6y!Vhu71ydhA7H7y`v%+lcMshenAaA%G<#0b zdvac$A?@z8{%z%M_Rf}%d>4c}F0P-ca@*QS`N5mQkXaMo$V-^@w>DHLG5!{*T;CBeoO^$_m%02oTgV`FbnUxyY8R!n{%?^z zIPJ!U-v;mMZ?-Qgo?^aLDY)zTEWuYk3I#b&o*t05ed&6VHSB25KgOO0&E_>ntA6zJ ziG-g&Xy9}1Um+Zw|@n$IPgO}fp_j1)5}c}2`{8yZ!G`yU`5!^ zCz}mBu1z?{5#e3>;#T3WGK078Sq&b~k5M{1yJT71qXP@W+v0hI_?KVUsM^5#T`YZP zXp#hD@G*;1ycaUq-58j81!tt}Qx`fJ?byruaGkKML%p+zYx?HN3rgH}vsq0RtMyg? zmm7VLcX8ItPvSi%V{fLf$@q2p!jHYn?YMlo3eG-wf9>R|^(UQQAA7jH`Q}$|{ZG$l zC2ZOJn)$>I<`Wzi)my9_7GGAlG12Av@xmSHetkddd)^lvPGB|hf-X7u!|dm4+6xUwbI!<~R5|`$F{Mtgpg1`IO4d45nn{G1`lBG`v<~lsT$Fk_)$&E}B2{qHbzj^6a zy;^LsZKph|@l@+YXQZS~Rkm?9F8=Frtnu7;(gp%(7C)_9tt{<{W((&Ytw(?)l1sV_Z-3rmg+KapU#ZNB@4mcaM#geYfK=pNYA7 z`s<~cGCi$TUu51qmpq`-*(AlCRBQN!J*nV)P#*gPzvuJPolZyZHL&7cWnZT=Tjpui zl3&3Yb>>HYl*?r06)bumAo6#g3KMVOnH_;UWe@nA^kO%d>hB^yl;WE%JfOq<72%qUc8X` zQx(@88yzFnve(dJ)t5~f{UXx`K`MZ|gWMV}Yd>aR=V zwzB?H{yU?r(^u^6#NVb32RS0XSgu~bVbSZq0&~0Lsh)f<~qxvQ$#r=_S1apv6I#A<44db6SubV9t9XP6tq z^YiojU$5KE=WqK}By3GYBd@gCfee$RmzS0**ZG@BJ?2p01f97U;jGeS>axh_s{E-_ z91Fr$AANOowP2@9!RNE)ppBI4_kNSonKp5|-~+b{DqXXGOwza}(3!+LafQ#<_D$QI zr{8{i40Mv2`n(FIXnjx6vC)SoycJv9)8HwxSRCxVHHt!0SN@xjW8DzL>a%FauRfQJ z5T9sPi5dtB3JRWV&v9pvWSp!4IzrLZWf90lf(C)JB=$~#TBhN40A#WX=gmF_BPO-L zo;Yw&fXq=&SY!e^H4_?^Jq?|pbCZQO$ug#}_$)#-A`^6gG01dN(D4_-pu^`?QmmOy zaLjN)HR2NZ>`aI&4VlzHLF##;jAeqLK_{vaOTi~m@=gT1lLd5|w9BLkZ#fi{6NHX} zj=F>!;U$a+jiVehIyyQ!R1|LuI=CHB=>i=<2{&Sj8cHaD!pGBTn{oqcfF9+Tp#eU> zG`j;Fn#u_vH%@994ah+kkdR^K8@F$t=2PIjxuZ}yXr&0~nClqxRuNFfi*Q!a3jV8j zvy8!3PT#sgCFSj{|Brg7w;-~+%5A;{e#=wYC$!vnYj#b5GZEs)Qz7*otHH@VL`_-t zfcd=&=IZb7M63G17iDM$c6;B!BuW;fn$i!pB&>4D@-kXujA+qi$f z{QTc{o*NWD^O3PANH}qSX^`ffy5Da(TAeEX{d^u66x8&3{eHQL9xkgO7s(V}TB0&5 za)yuE#^UFGM;;gIL~Y>!-GCAEbxZMczbjX-9(?@qz|$f{x&GvLcXn3%dbxaK&dp7U z^Ve_LA|lt%{`Bcn50y#Hd+$HK{`2S04N8;n16&{lN8W73d^6US3{hyZhz)tHtM4 zy;|J@hjPzN*Z0TozN>eQbA9W&yLsYC8)xL)I-LD9_WASY&g-x9+f94@ z?L$=M^}OxUpe^A`%og$`Xiaqr4!-=LUR+f4VAj^YA7ADA*;`v%o#jrSKJD!8e!Sy; z@%z2sAKa^cKhxtgsJO@w5uU{?(d!l%8fy9@{ro)J#dWogj*NTa-cLCGr_TOP-glR{ zc_4)>J}a0rjApj+q|dK@w{vk_ZE2|~=rq22_u?2FlqMdSrW?)n{3__+>%ea*PgN%T z{QUg!uGi}h2m9N$ifV@?9B$(^km`-tRg!52G0<<))ZTAfqIAR7$Mt%sG|kT6C;9d3 z*NV4WuP1Db*!1xesAw(`I=QV8bkMa@>VrUOkcie{KJax&r#3KqxVX41F?s{OJZU8l z_!MOqtkS%X7x&v8`t|j7@m-!EgKGjOXRfhsczu2S@vYb6y1n)HisIbx%DUMu zT$`%D=Vh48QmoJ0`E**0dAEpy1ka;%``?n^-rO{P@-afEO+2oG(R-TCL0D_^~ET7Q4dt=gB@_x~&Hi}EbGdg)Tr?)UqoOG`_4e7lvsG4JlKgZn4HvYRq> zbAx))>uYPlXJDL;|0fi-TD1Ibskn|9x847r=b!J&Tq1D1Pu4j%S9f0ZyPXG*fBN*P z=!Q~I6KGNK>hSe$-rn5L&d;}h@;x$GZHTi?6S5UzB%| zNffWFZIuYGqo?ow4|2W8BH=%pk zB%j_8XNG|I`2KGFeKWqRN?mN1uj^Rs-v8)fyL{WtoMVrUc5f_w9rovKe!X*S?A(;^ ziZ|U3Ja@B7PzkK`l zY*(a8@z$+dBQ_*3c6E2R+uXZ%uW#Bkv5@d^emk@Ldp2)wZPmW*t1?r1WgA1n?{9A_ zKA$x&JZJg*L6`Qr1K#?3TPkekHHS@#IXlagTiPtA;bsmSL(Abr8OtJ;;AK7&zu2jq z=1b5JadLKMUK95EL9$pkE5i>qzsKA1@7wK&(YsUs|1apmfeodv!)kutzW*#HXx4ac~YR%vbTRX%EISHETd{`wCa3oLA&Ot!d~F~?6`sFP*Ot;UVh zKcBakKbRnJ-tMyg0!cdm1=QXc<^|SoAMk#b?bT)4H2$~-toVL z-jrE`u5_B6w~Mpt?}3WBQ@$%k_AxLW04=O7dg^s%rZKx%H>dGHrd~`VTm@+MWUL_moyjA6X znK`<+7hH=>_bijYv?OTd zl#pI<`|{`Hr!wGMrzWz4E}!#Ay}<~&A}C0bH0hu6W^9_wzmKP literal 0 HcmV?d00001 diff --git a/doc/notes.typ b/doc/notes.typ index fa0b780..ecc9b8f 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -335,7 +335,7 @@ composition breaks this: line((-0.5, 0.5), (0, 0.5), mark: (end: "straight")) line((1, 0.5), (2, 0.5), mark: (end: "straight")) line((3, 0.5), (3.5, 0.5), mark: (end: "straight")) -}), +}), $ M_i & : & a_1 & | & bot & | \ M_o & : & b_1 & | & b_2 & | \ N_i & : & b_1 & | & b_2 & | @@ -353,7 +353,7 @@ needed. This is done through the following composition: line((-0.5, 0.5), (0, 0.5), mark: (end: "straight")) line((1, 0.5), (2, 0.5), mark: (end: "straight")) line((3, 0.5), (3.5, 0.5), mark: (end: "straight")) -}), +}), $ I &: & a_1 & | & bot & | & bot & & & | & bot & | & bot & & & | \ M_i &: & a_1 & | & & | & & | & bot & | & & | & & | & bot & | \ M_o &: & b_1 & | & & | & & | & b_2 & | & & | & & | & bot & | @@ -620,11 +620,11 @@ simulation loop as follows: models has much in common with code generation for synchronous languages. ]) -In [Branicky, 1995b], dynamical systems are defined as follows: +In [Branicky, 1995b], dynamical systems are defined as follows: #block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [ A continuous (resp. discrete) dynamical system defined on the topological - space $X$ over the semigroup $bb(R)^+$ (resp. $bb(Z)^+$) is a function + space $X$ over the semigroup $bb(R)^+$ (resp. $bb(Z)^+$) is a function $f : X times bb(R)^+ -> X$ (resp. $f : X times bb(Z)^+ -> X$) with the following three properties: 1. initial condition: $f(p, 0) = p$, @@ -694,7 +694,7 @@ $f : S u p e r d e n s e(bb(V))$, $f(t, n) = f(t + n partial)$ in the non-standard semantics. Our representation instead uses -$S i g n a l(bb(V)) = S t r e a m((h : bb(R)) * ([0, h] -> bb(V)))$. Can we +$S i g n a l(bb(V)) = S t r e a m((h : bb(R)) * ([0, h] -> bb(V)))$. Can we convert between the two as we want? #breakl(```ocaml diff --git a/doc/pres.typ b/doc/pres.typ index 88c12a8..f4b52bd 100644 --- a/doc/pres.typ +++ b/doc/pres.typ @@ -145,7 +145,7 @@ Solution: simulate assertions with their own solver. Nodes take an additional reset parameter `'p` for their reset function: #align(top)[ - #grid(columns: (3fr, 2fr), align: (left, left), + #grid(columns: (3fr, 2fr), align: (left, left), ```ocaml type ('p, 'a, 'b) dnode = DNode : { @@ -200,7 +200,7 @@ Hybrid nodes are a bit more complex: #slide(repeat: 3, self => [ #let (uncover, only) = utils.methods(self) - ODE solvers are discrete nodes producing streams of functions defined on + ODE solvers are discrete nodes producing streams of functions defined on successive intervals: $(h: bb(R)_+) -->^D (h': bb(R)_+) times (u: [0,h'] -> bb(V))$. @@ -249,7 +249,7 @@ Hybrid nodes are a bit more complex: #slide(repeat: 3, self => [ #let (uncover, only) = utils.methods(self) - Zero-crossing solvers are discrete nodes too, looking for zero-crossings on + Zero-crossing solvers are discrete nodes too, looking for zero-crossings on functions: $(h: bb(R)_+) times (u: [0, h] -> bb(V)) -->^D (h': bb(R)_+) times (z: bb(Z)?)$. @@ -456,7 +456,7 @@ When receiving a new value, store it box((10, 1.25), (3, 1.5), double: true, content: `state`) // (h, u) -> state - content((-2.25, 3), $(h, u)$) + content((-2.25, 3), $(h, u)$) arrow((-1, 3), (r, 1.5), (d, 2.5), (r, 5), (u, 1.5), (r, 4.5)) // sim -> bot @@ -583,7 +583,7 @@ When receiving a new value, store it // solver.step box((6.5, 1), (3, 1.5), content: `step`) - + // solver.state box((6.5, 2.75), (3, 1.5), double: true, content: `state`) @@ -757,7 +757,7 @@ If no zero-crossing is found, there is no discontinuity. content((-3, 1.5), $"Signal"(alpha)$) content((5.5, 2), $"Signal"(beta)$) content((14, 1.5), $"Signal"(gamma)$) - + content((-4, -1), $alpha:$) content((-4, -2.5), $beta:$) content((-4, -4), $gamma:$) @@ -795,7 +795,7 @@ If no zero-crossing is found, there is no discontinuity. The output signal stops at least as often as the input signal does. #linebreak() This calls for a "lazy" composition: we request rather than provide data. - + #align(center, cetz-canvas({ import cetz.draw: * box((0, 0), (3, 3), content: $M_1$) @@ -809,7 +809,7 @@ If no zero-crossing is found, there is no discontinuity. content((-3, 1.5), $"Signal"(alpha)$) content((5.5, 2), $"Signal"(beta)$) content((14, 1.5), $"Signal"(gamma)$) - + content((-4, -1), $alpha:$) content((-4, -2.5), $beta:$) content((-4, -4), $gamma:$) diff --git a/doc/rep.typ b/doc/rep.typ new file mode 100644 index 0000000..fe4061f --- /dev/null +++ b/doc/rep.typ @@ -0,0 +1,192 @@ +#import "@preview/cetz:0.4.0" +#import "@preview/finite:0.4.1" + +#set page(paper: "a4") +#set par(justify: true) +#set raw(syntaxes: "zelus.sublime-syntax") +#show raw: set text(font: "CaskaydiaCove NF") + +#let lustre = smallcaps[Lustre] +#let esterel = smallcaps[Esterel] +#let simulink = smallcaps[Simulink] +#let sundials = smallcaps[Sundials CVODE] +#let zelus = smallcaps[Zélus] +#let TODO(..what) = { + let msg = if what.pos().len() == 0 { "" } else { ": " + what.pos().join("") } + block(fill: red, width: 100%, inset: 5pt, align(center, raw("TODO" + msg))) +} +#let adot(s) = $accent(#s, dot)$ +#let addot(s) = $accent(#s, dot.double)$ + +#align(center)[ + #text(16pt)[ + *Internship Report --- MPRI M2 \ + Hybrid System Models with Transparent Assertions* + ] + + Henri Saudubray, supervised by Marc Pouzet, Inria PARKAS +]\ + +#heading(outlined: false)[General Context] +What is the report about? Where does the problem come from? What is the state of +the art in that area? + +Hybrid systems modelers such as #simulink are essential tools in the development +of embedded systems evolving in physical environments. They allow for precise +descriptions of hybrid systems through both continuous-time behaviour defined +through Ordinary Differential Equations (ODEs) and discrete-time reactive +behaviour similar to what is found in the synchronous languages such as #lustre. +The #zelus language aims to reconcile synchronous languages with hybrid systems, +by taking a synchronous language kernel and extending it with continuous-time +constructs similar to what is found in tools like #simulink. Continuous-time +behaviour is computed through the use of an off-the-shelf ODE solver such as +#sundials. + +#heading(outlined: false)[Research Problem] +What specific question(s) have you studied? What motivates the question? What +are the applications/consequences? Is it a new problem? Why did you choose this +problem? + +The simulation of hybrid system models, as done in #simulink and #zelus, uses a +single ODE solver instance to simulate the entire model. This has various +advantages: it is less computationally intensive, and simplifies the work of the +compiler. Unfortunately, it also raises a difficult problem: sub-systems which +seemingly should not interfere with each other end up affecting each other's +results. This is due to the chosen integration method: an adaptive solver like +#sundials will vary its step length throughout the integration process, and the +addition of ODEs in the overall system can influence these step lengths, +affecting the results obtained for pre-existing ODEs. This is particularly +problematic in the case of runtime assertions, which are typically expected to +be transparent: they should not affect the final result of the computation. + +#heading(outlined: false)[Proposed Contributions] +What is your answer to the research problem? Please remain at a high level; +technical details should be given in the main body of the report. Pay special +attention to the description of the _scientific_ approach. + +To solve this, we propose a new runtime for the #zelus language, which simulates +assertions with their own solvers in order to maintain the separation between +assertions and the model they operate on. + +#TODO() + +#heading(outlined: false)[Arguments Supporting Their Validity] +What is the evidence that your solution is a good solution? (Experiments? +Proofs?) Comment on the robustness of your solution: does it rely on +assumptions, and if so, to what extent? + +#TODO("Justify") + +#heading(outlined: false)[Summary and Future Work] +What did you contribute to the area? What comes next? What is a good _next_ step +or question? + +#TODO("Next steps") + +#pagebreak() +#outline() +#pagebreak() + += Background + +== Hybrid systems + +Hybrid systems are dynamical systems whose behaviour evolves through both +continuous and discrete phases. They are used to model software interacting with +physical systems. Continuous phases are described using ordinary differential +equations (ODEs), while discrete phases can be represented as a reactive +program in a synchronous language such as #lustre or #esterel. + +As a first example, say we wish to model a bouncing ball. We could start by +describing its distance from the ground $y$ with a second-order differential +equation +$ addot(y) = -9.81 $ +where $addot(y)$ denotes the second order derivative of $y$ with +respect to time (the acceleration of the ball), and $9.81$ is the gravitational +constant $g$: the acceleration of the ball is its negation. We now give the +initial position and speed of the ball: +$ y(0) = 50 space space space adot(y)(0) = 0 $ +We have just described an initial value problem: given an ODE and an initial +value for its dependent variable, its solution is a function $y(t)$ returning +the value of the variable at a given time $t$. We can find an approximation of +this function using an ODE solver, such as #sundials. + +As of right now, our ball will fall until the end of time; we have not said +anything about how it bounces when it hits the floor. To do so, we need a notion +of discrete _events_. These are modelled by zero-crossings: we monitor a certain +value and stop when it goes from strictly negative to positive or null. For our +purposes, we choose $-y$ as the monitored value, and call the zero-crossing +event $z$. When $z$ occurs (i.e., when the ball touches the ground), we set the +speed $adot(y)$ to $-k dot "last"(adot(y))$, where $"last"(y)$ denotes the left +limit of $y$ (we cannot specify $adot(y)$ in terms of itself), and $k$ is a +factor modelling the loss of inertia due to the collision (say, $0.8$). We can +then resume the approximation of the solution. + +@lst:ball.zls shows how such a model might be expressed in the concrete syntax +of #zelus. + +#figure(placement: top, gap: 2em, + ```zelus + let hybrid ball () = y where + rec der y = y' init 50.0 + and der y' = -9.81 init 0.0 reset z -> -0.8 *. (last y') + and z = up (-. y) + ```, + caption: [The bouncing ball in #zelus] +) + +More formally, a hybrid system can be described as an automaton + +== Executing models + +Executing such a model is quite simple. There are two modes of execution: +discrete and continuous. In continuous mode, we call the ODE solver in order to +obtain an approximation of the variables defined through ODEs, and monitor for +zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in which +we perform computation steps as needed, until no other zero-crossing occurs, in +which case we go back to the continuous mode, and repeat, as seen in @automaton. + +#figure(finite.automaton( + (D: (D: "cascade", C: "no cascade"), + C: (C: "no zero-crossing", D: "zero-crossing")), + initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3) +), caption: [High-level overview of the runtime], placement: top) + += Runtime +To solve this issue, we need to redefine what the runtime of our hybrid system +models is. The runtime is the piece of software that handles the pairing of a +hybrid model with a solver and the different execution phases that need to take +place. + +To allow for independent simulation of the assertions, we need to simulate them +with their own ODE solvers. Since assertions can themselves contain ODEs, they +can be viewed as hybrid systems themselves. Assertions can also contain their +own sub-assertions, and so recursively: our runtime needs to handle this +accordingly. + +== Hybrid models with assertions +A hybrid model with assertions can be defined using a recursive data type: +```ocaml +type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a = + HNodeA : { + body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec; + assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list; + } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a +``` + +Notice that the list of assertions takes as input the inner state of the parent +model: assertions are checked on the model state, and any variable which is +required by the assertion becomes a state variable. + +== Solvers as synchronous nodes +== Simulations as synchronous nodes +#TODO("talk about the new runtime") + += Assertions +#TODO("talk about how assertions are done") + +#pagebreak() += Annex +#set heading(level: 2) +#bibliography("sources.bib", full: true) +#set heading(level: 1) diff --git a/exm/ball.ml b/exm/ball.ml index dcdd63c..72e804e 100644 --- a/exm/ball.ml +++ b/exm/ball.ml @@ -40,7 +40,7 @@ let step ({ zin; lx; _ } as s) zfalse = of_array [| -. 0.8 *. lx.{0}; 0.0; lx.{2}; lx.{3} |] else lx in of_array [| s.lx.{1} |], { zin=zfalse; lx; i=false } -let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode = +let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec = let yd = cmake csize in let zout = cmake zsize in let zfalse = zmake 1 in @@ -49,10 +49,10 @@ let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode = let step s _ = step s zfalse in let state = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in let reset _ _ = state in - HNode { state; fder; fzer; fout; step; reset; horizon; - jump; cset; cget; zset; csize; zsize } + { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset; + csize; zsize } let errmsg = "Too many arguments for the model (needed: 0)" let init = function - | [] -> bouncing_ball () + | [] -> HNode (bouncing_ball ()) | _ -> raise (Invalid_argument errmsg) diff --git a/exm/dune b/exm/dune index 9eb098f..4eecb1e 100644 --- a/exm/dune +++ b/exm/dune @@ -1,3 +1,5 @@ (library - (name examples) - (libraries hsim solvers)) + (name examples) + (libraries hsim solvers)) + +(include_subdirs unqualified) diff --git a/exm/sin1x.ml b/exm/sin1x.ml new file mode 100644 index 0000000..a6f2447 --- /dev/null +++ b/exm/sin1x.ml @@ -0,0 +1,72 @@ + +(* Example due to JB. Jeannin - zero-crossing detection *) +(* the signal x(t) is expected to have a zero-crossing at [t = 0.3] *) + +(* to draw: ./main.exe -s sin_1_x -stop 2 -sample 10000 | feedgnuplot --stream --domain --lines *) + +open Hsim.Types +open Solvers.Zls + +(* let hybrid sin_1_x() = + let der t = 1.0 init 0.0 in + let t = t -. 0.3 in + let v = sin(1/t) in + let o = t *. (v *. v) in [that is: t *. (sin(1/t)^2] + present up(o) -> 1.0 else 0.0, o *) + +let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a + +type ('a, 'b) state = { s_x: 'a; zin: 'b } + +let csize = 1 +let zsize = 1 + +let fder _ _ yd = yd.{0} <- 1.0 + +let fzer _ y zout = + let t = y.{0} in + let t = t -. 0.3 in + let v = sin(1.0 /. t) in + let o = t *. (v *. v) in + zout.{0} <- o + +let fout s y = + let z = if s.zin.{0} = 1l then 1.0 else 0.0 in + let t = y.{0} in + let t = t -. 0.3 in + let v = sin(1.0 /. t) in + let o = t *. (v *. v) in + of_array [| z; o |] + +let step ({ s_x; zin } as s) zfalse = + let z = if zin.{0} = 1l then 1.0 else 0.0 in + let t = s_x.{0} in + let t = t -. 0.3 in + let v = sin(1.0 /. t) in + let o = t *. (v *. v) in + of_array [| z; o |], { s with zin = zfalse } + +let reset _ s = s + +let cget { s_x; _ } = s_x +let cset s l_x = { s with s_x = l_x } +let zset s zin = { s with zin } +let jump _ = true +let horizon _ = max_float + +let sin_1_x () = + let zfalse = zmake 1 in + let yd = cmake 1 in + let zout = cmake 1 in + let fder _ _ y = fder 0.0 y yd; yd in + let fzer _ _ y = fzer 0.0 y zout; zout in + let fout s _ y = fout s y in + let step s _ = step s zfalse in + let state = { s_x = of_array [| 0.0 |]; zin = zfalse } in + { state; fder; fzer; fout; step; reset; horizon; + cset; cget; zset; csize; zsize; jump } + +let errmsg = "Too many arguments for the model (needed: 0)" +let init = function + | [] -> HNode (sin_1_x ()) + | _ -> raise (Invalid_argument errmsg) diff --git a/exm/sin1x_der.ml b/exm/sin1x_der.ml new file mode 100644 index 0000000..e4bf116 --- /dev/null +++ b/exm/sin1x_der.ml @@ -0,0 +1,79 @@ + +(* Example due to JB. Jeannin - zero-crossing detection. *) + +(* The same as [sin1x.ml] but with an integrator. + + d(t . (sin(1/t))^2))/dt = sin(1/t)(sin(1/t) - (2 . cos(1/t)) / t) + + d(sin(f(x)))/dx = cos(f(x)) f'(x) *) + +open Hsim.Types +open Solvers.Zls + +(* let hybrid sin_1_x() = + let der t0 = 1.0 init 0.0 in + let d = 0.6 in + let t = t0 -. d in + let der o = sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t) + init (-. d *. (sin(1.0 /. (-. d))) ** 2.0) in + present up(o) -> 1.0 else 0.0, o *) + +let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a + +type ('a, 'b) state = { s_x: 'a; zin: 'b } + +let d = 0.66 +let y0 = -. d *. let v = (sin(1.0 /. (-. d))) in v *. v +let csize = 2 +let zsize = 1 + +let fder d y yd = + yd.{0} <- 1.0; + let t = y.{0} -. d in + yd.{1} <- sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t) + +let fzer z y zout = + let o = y.{1} in + zout.{0} <- if z then o else 1.0 + +let fout s y = + let z = if s.zin.{0} = 1l then 1.0 else 0.0 in + let o = y.{1} in + of_array [| z; o |] + +let step ({ s_x; zin } as s) zfalse = + let z = if zin.{0} = 1l then 1.0 else -1.0 in + let o = s_x.{1} in + of_array [| z; o |], { s with zin = zfalse } + +let reset _ s = s + +let cget { s_x; _ } = s_x +let cset s l_x = { s with s_x = l_x } +let zset s zin = { s with zin } +let horizon _ = max_float +let jump _ = true + +let sin_1_x z d = + let zfalse = zmake 1 in + let yd = cmake 2 in + let zout = cmake 1 in + let fder _ _ y = fder d y yd; yd in + let fzer _ _ y = fzer z y zout; zout in + let fout s _ y = fout s y in + let step s _ = step s zfalse in + let state = { s_x = of_array [| 0.0; y0 |]; zin = zfalse } in + { state; fder; fzer; fout; step; horizon; + cset; cget; zset; csize; zsize; reset; jump } + +let errmsg_many = "Too many arguments to model (needed: [bool, float])" +let errmsg_few = "Too few arguments to model (needed: [bool, float])" +let errmsg_invalid = "Invalid arguments to model (needed: [bool, float])" + +let init = function + | [zer; d] -> + let zer, d = try bool_of_string zer, float_of_string d + with _ -> raise (Invalid_argument errmsg_invalid) in + HNode (sin_1_x zer d) + | [] | [_] -> raise (Invalid_argument errmsg_few) + | _ -> raise (Invalid_argument errmsg_many) diff --git a/exm/sincos.ml b/exm/sincos.ml index 17131c6..ab62149 100644 --- a/exm/sincos.ml +++ b/exm/sincos.ml @@ -34,11 +34,11 @@ let sinus_cosinus theta0 omega = HNode { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset; csize; zsize } -let errmsg_invalid = "Invalid arguments to model (needed: 2 floats)" -let errmsg_few = "Too few arguments to model (needed: 2 floats)" -let errmsg_many = "Too many arguments to model (needed: 2 floats)" +let errmsg_invalid = "Invalid arguments to model (needed: [float, float])" +let errmsg_few = "Too few arguments to model (needed: [float, float])" +let errmsg_many = "Too many arguments to model (needed: [float, float])" let init = function - | [t0; om] -> + | [om; t0] -> let t0, om = try float_of_string t0, float_of_string om with Failure _ -> raise (Invalid_argument errmsg_invalid) in sinus_cosinus t0 om diff --git a/exm/zelus/ballz/ballz.ml b/exm/zelus/ballz/ballz.ml new file mode 100644 index 0000000..08df6f6 --- /dev/null +++ b/exm/zelus/ballz/ballz.ml @@ -0,0 +1,131 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +let (+=) r v = r := !r + v + +let g = 9.81 + +let y0 = 50. + +let y'0 = 0. + +type ball = + { mutable major : bool; + mutable h : float; + mutable init : bool; + mutable z : (bool, float) zerocrossing; + mutable y' : float continuous; + mutable y : float continuous } + +let ball cstate = + + let ball_alloc _ = + cstate.cmax <- cstate.cmax + 2; + cstate.zmax <- cstate.zmax + 1; + { major = false; + h = 42.; + init = false; + z = { zin = false; zout = 1. }; + y' = { pos = 42.; der = 0. }; + y = { pos = 42.; der = 0. } + } in + + let ball_step self (_, ()) = + let cidx = cstate.cindex in let cpos = ref cidx in + let zidx = cstate.zindex in let zpos = ref zidx in + cstate.cindex <- cstate.cindex + 2; + cstate.zindex <- cstate.zindex + 1; + self.major <- cstate.major; + if cstate.major then + for i = cidx to 1 do Zls.set cstate.dvec i 0. done + else begin + self.y'.pos <- Zls.get cstate.cvec !cpos; cpos += 1; + self.y.pos <- Zls.get cstate.cvec !cpos; cpos += 1 + end; + let res0, res1 = + let encore = ref false in + if self.init then self.y'.pos <- y'0; + let last_y' = self.y'.pos in + if self.z.zin then begin + encore := true; self.y'.pos <- -0.8 *. last_y' + end; + self.h <- if !encore then 0. else infinity; + if self.init then self.y.pos <- y0; + self.init <- false; + self.z.zout <- -. self.y.pos; + self.y'.der <- -. g; + self.y.der <- self.y'.pos; + self.y.pos, self.y.der + in + cstate.horizon <- min cstate.horizon self.h; + cpos := cidx; + if cstate.major then begin + Zls.set cstate.cvec !cpos self.y'.pos; cpos += 1; + Zls.set cstate.cvec !cpos self.y.pos; cpos += 1; + self.z.zin <- false + end else begin + self.z.zin <- Zls.get_zin cstate.zinvec !zpos; zpos += 1; + zpos := zidx; + Zls.set cstate.zoutvec !zpos self.z.zout; zpos += 1; + Zls.set cstate.dvec !cpos self.y'.der; cpos += 1; + Zls.set cstate.dvec !cpos self.y.der; cpos += 1 + end; + Bigarray.(Array1.of_array Float64 c_layout [| res0; res1 |]) in + + let ball_reset self = self.init <- true in + + Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } + +type ('f , 'e , 'd , 'c , 'b , 'a) _main = + { mutable i_73 : 'f ; + mutable major_62 : 'e ; + mutable h_72 : 'd ; + mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a } + +let main (cstate_80:Ztypes.cstate) = + let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball + cstate_80 in + let main_alloc _ = + cstate_80.cmax <- (+) cstate_80.cmax 1; + { major_62 = false ; + h_72 = 42. ; + i_70 = (false:bool) ; + h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. }; + i_73 = i_73_alloc () (* continuous *) } in + let main_step self ((time_61:float) , ()) = + ((let (cindex_81:int) = cstate_80.cindex in + let cpos_83 = ref (cindex_81:int) in + cstate_80.cindex <- (+) cstate_80.cindex 1 ; + self.major_62 <- cstate_80.major ; + (if cstate_80.major then + for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done + else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ; + cpos_83 := (+) !cpos_83 1))) ; + (let (result_85) = + let h_71 = ref (infinity:float) in + (if self.i_70 then self.h_68 <- (+.) time_61 0.) ; + (let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in + self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ; + h_71 := min !h_71 self.h_68 ; + self.h_72 <- !h_71 ; + self.i_70 <- false ; + self.t_63.der <- 1. ; + (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in + (begin match z_69 with + | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64 + | _ -> () end) ; + Bigarray.(Array1.create Float64 c_layout 0))) in + cstate_80.horizon <- min cstate_80.horizon self.h_72 ; + cpos_83 := cindex_81 ; + (if cstate_80.major then + (((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ; + cpos_83 := (+) !cpos_83 1))) + else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ; + cpos_83 := (+) !cpos_83 1)))) ; result_85))) in + let main_reset self = + ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ): + unit) in + Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.ml.bak b/exm/zelus/ballz/ballz.ml.bak new file mode 100644 index 0000000..d248fc0 --- /dev/null +++ b/exm/zelus/ballz/ballz.ml.bak @@ -0,0 +1,137 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +let g = 9.81 + +let y0 = 50. + +let y'0 = 0. + +type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball = + { mutable major_50 : 'g ; + mutable h_60 : 'f ; + mutable h_58 : 'e ; + mutable i_56 : 'd ; + mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a } + +let ball (cstate_74:Ztypes.cstate) = + + let ball_alloc _ = + cstate_74.cmax <- (+) cstate_74.cmax 2 ; + cstate_74.zmax <- (+) cstate_74.zmax 1; + { major_50 = false ; + h_60 = 42. ; + h_58 = (42.:float) ; + i_56 = (false:bool) ; + x_55 = { zin = false; zout = 1. } ; + y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in + let ball_step self ((_time_49:float) , ((y0_48:float) , (y'0_47:float))) = + Printf.printf "STEP (%d)\n" cstate_74.cindex; + ((let (cindex_75:int) = cstate_74.cindex in + let cpos_77 = ref (cindex_75:int) in + let (zindex_76:int) = cstate_74.zindex in + let zpos_78 = ref (zindex_76:int) in + cstate_74.cindex <- (+) cstate_74.cindex 2 ; + cstate_74.zindex <- (+) cstate_74.zindex 1 ; + self.major_50 <- cstate_74.major ; + (if cstate_74.major then + for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done + else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1) ; + (self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1))) ; + (let (result_79:float) = + let h_59 = ref (infinity:float) in + let encore_57 = ref (false:bool) in + (if self.i_56 then self.y'_52.pos <- y'0_47) ; + (let (l_54:float) = self.y'_52.pos in + (begin match self.x_55.zin with + | true -> + encore_57 := true ; + self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end) + ; + self.h_58 <- (if !encore_57 then 0. else infinity) ; + h_59 := min !h_59 self.h_58 ; + self.h_60 <- !h_59 ; + (if self.i_56 then self.y_51.pos <- y0_48) ; + self.i_56 <- false ; + self.x_55.zout <- (~-.) self.y_51.pos ; + self.y'_52.der <- (~-.) g ; + self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in + cstate_74.horizon <- min cstate_74.horizon self.h_60 ; + cpos_77 := cindex_75 ; + (if cstate_74.major then + (((Printf.printf "idx: %d\n" !cpos_77; + Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ; + cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false))) + else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ; + zpos_78 := (+) !zpos_78 1)) ; + zpos_78 := zindex_76 ; + ((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ; + zpos_78 := (+) !zpos_78 1)) ; + ((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.dvec !cpos_77 self.y_51.der ; + cpos_77 := (+) !cpos_77 1)))) ; result_79)):float) in + let ball_reset self = + (self.i_56 <- true:unit) in + Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } +type ('f , 'e , 'd , 'c , 'b , 'a) _main = + { mutable main_ball : 'f ; + mutable main_major : 'e ; + mutable h_72 : 'd; mutable i_70 : 'c; mutable h : 'b; + mutable t_63 : 'a } + +let main cs = + let Node + { alloc = ball_alloc; + step = ball_step; + reset = ball_reset } = ball cs in + + let main_alloc _ = + cs.cmax <- cs.cmax + 1; + { main_major = false; + h_72 = 42.0; i_70 = false; h = 42.0; + t_63 = { pos = 42.; der = 0. }; + main_ball = ball_alloc () } in + + let main_step self (time, ()) = + let cindex = cs.cindex in + let cpos = ref cindex in + Printf.printf "main:cindex: %d\n" cs.cindex; + cs.cindex <- cs.cindex + 1; + self.main_major <- cs.major; + if cs.major then for i = cindex to 0 do Zls.set cs.dvec i 0. done + else begin self.t_63.pos <- Zls.get cs.cvec !cpos; cpos := !cpos + 1 end; + let result = + if self.i_70 then self.h <- time; + let z = self.main_major && (time >= self.h) in + if z then self.h <- self.h +. 0.01; + self.h_72 <- min infinity self.h; + self.i_70 <- false; + self.t_63.der <- 1.; + let y_64 = ball_step self.main_ball (time, (y0, y'0)) in + if z then begin + print_float self.t_63.pos; + print_string "\t"; + print_float y_64; + print_newline () + end; + Bigarray.(Array1.create Float64 c_layout 0) in + cs.horizon <- min cs.horizon self.h_72; + cpos := cindex; + if cs.major then Zls.set cs.cvec !cpos self.t_63.pos + else Zls.set cs.dvec !cpos self.t_63.der; + result in + + let main_reset self = + self.i_70 <- true; + self.t_63.pos <- 0.; + ball_reset self.main_ball in + + Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.ml.bak2 b/exm/zelus/ballz/ballz.ml.bak2 new file mode 100644 index 0000000..05842bb --- /dev/null +++ b/exm/zelus/ballz/ballz.ml.bak2 @@ -0,0 +1,135 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +let g = 9.81 + +let y0 = 50. + +let y'0 = 0. + +type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball = + { mutable major_50 : 'g ; + mutable h_60 : 'f ; + mutable h_58 : 'e ; + mutable i_56 : 'd ; + mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a } + +let ball (cstate_74:Ztypes.cstate) = + + let ball_alloc _ = + cstate_74.cmax <- (+) cstate_74.cmax 2 ; + cstate_74.zmax <- (+) cstate_74.zmax 1; + { major_50 = false ; + h_60 = 42. ; + h_58 = (42.:float) ; + i_56 = (false:bool) ; + x_55 = { zin = false; zout = 1. } ; + y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in + + let ball_step self ((_time_49:float) , ()) = + ((let (cindex_75:int) = cstate_74.cindex in + let cpos_77 = ref (cindex_75:int) in + let (zindex_76:int) = cstate_74.zindex in + let zpos_78 = ref (zindex_76:int) in + cstate_74.cindex <- (+) cstate_74.cindex 2 ; + cstate_74.zindex <- (+) cstate_74.zindex 1 ; + self.major_50 <- cstate_74.major ; + (if cstate_74.major then + for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done + else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1) ; + (self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1))) ; + (let (result_79:float) = + let h_59 = ref (infinity:float) in + let encore_57 = ref (false:bool) in + (if self.i_56 then self.y'_52.pos <- y'0) ; + (let (l_54:float) = self.y'_52.pos in + begin match self.x_55.zin with + | true -> + encore_57 := true ; + self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end; + self.h_58 <- (if !encore_57 then 0. else infinity) ; + h_59 := min !h_59 self.h_58 ; + self.h_60 <- !h_59 ; + (if self.i_56 then self.y_51.pos <- y0) ; + self.i_56 <- false ; + self.x_55.zout <- (~-.) self.y_51.pos ; + self.y'_52.der <- (~-.) g ; + self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in + cstate_74.horizon <- min cstate_74.horizon self.h_60 ; + cpos_77 := cindex_75 ; + (if cstate_74.major then + (((Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ; + cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false))) + else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ; + zpos_78 := (+) !zpos_78 1)) ; + zpos_78 := zindex_76 ; + ((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ; + zpos_78 := (+) !zpos_78 1)) ; + ((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.dvec !cpos_77 self.y_51.der ; + cpos_77 := (+) !cpos_77 1)))) ; + Bigarray.(Array1.of_array Float64 c_layout [| result_79 |])))) in + + let ball_reset self = + (self.i_56 <- true:unit) in + + Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } + +type ('f , 'e , 'd , 'c , 'b , 'a) _main = + { mutable i_73 : 'f ; + mutable major_62 : 'e ; + mutable h_72 : 'd ; + mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a } + +let main (cstate_80:Ztypes.cstate) = + let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball + cstate_80 in + let main_alloc _ = + cstate_80.cmax <- (+) cstate_80.cmax 1; + { major_62 = false ; + h_72 = 42. ; + i_70 = (false:bool) ; + h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. }; + i_73 = i_73_alloc () (* continuous *) } in + let main_step self ((time_61:float) , ()) = + ((let (cindex_81:int) = cstate_80.cindex in + let cpos_83 = ref (cindex_81:int) in + cstate_80.cindex <- (+) cstate_80.cindex 1 ; + self.major_62 <- cstate_80.major ; + (if cstate_80.major then + for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done + else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ; + cpos_83 := (+) !cpos_83 1))) ; + (let (result_85) = + let h_71 = ref (infinity:float) in + (if self.i_70 then self.h_68 <- (+.) time_61 0.) ; + (let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in + self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ; + h_71 := min !h_71 self.h_68 ; + self.h_72 <- !h_71 ; + self.i_70 <- false ; + self.t_63.der <- 1. ; + (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in + (begin match z_69 with + | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64 + | _ -> () end) ; + Bigarray.(Array1.create Float64 c_layout 0))) in + cstate_80.horizon <- min cstate_80.horizon self.h_72 ; + cpos_83 := cindex_81 ; + (if cstate_80.major then + (((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ; + cpos_83 := (+) !cpos_83 1))) + else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ; + cpos_83 := (+) !cpos_83 1)))) ; result_85))) in + let main_reset self = + ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ): + unit) in + Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.zci b/exm/zelus/ballz/ballz.zci new file mode 100644 index 0000000000000000000000000000000000000000..08fa125c8a8c6a515d7c97cca6557fe41356e53b GIT binary patch literal 6309 zcmZpfx@;cI57NoU}WTVT(DsZ3nwEZpQ8gKBfkS9Bah>TDJ*ZaA5e)km2Cq;N-yYA0d+N;Naj4GRa}V#08UBSQkuU;dGe9 z!ZBfk16TyCb%GPae`g1fN{0UnCb6(MOo9kIF#LA`34=^r_Q1(u!G*~SE=*bAqf}|I zpt^R$6c%F-M~G4f$0`S?rH+*gSXkg32A2gatO%YHOc8^#14x;J1ITO^4!DXMxB^EQ zs}9Csh=#L#;VeHm%OB1PfV1G{GQ`5g;v5#Rh=E)Mw%lO>3kS%L4F4S#uy8v#I82zZ zz)i6V%#~ZfA`F)kgR>+dESQ1yFa|>mj8y?+RD*r#2#OPgQz93zusbYZQG+RBh+4qH z2IoY;Ss`#%AeU2FIs^BSSD;Q|JN~6^8{ZtP5BW z9>eEZ77kEgEMSp`o9zi@fh}IZq5}y#C*?|Y190g1E?_Z%Dp|ndzko#qq+sF#SEWjI zhW`s#j25t1IV@m_U%(OtQMG`@$YB9X9E9VrfF%*6egX?O*swH+@B)_H1uO{*SfUrO zApD;K^HeR2!Qi@prF;PkDD>GN8XXp}aKdHN7O-$REMTc{aC8J`4N#oBC{}_kLYR@f zfJGD$X^Ef|12-W7g_nfFOGV+ugTt5;6nIdVK*Ng_nG0brV1c<7ss9v|$R%Q%}c9EHwyqYaAA^ zoCd3NSWv~_3@(RQ_AOxHb?|aK{|*;Ge1`=W4zfsI*tWpI;ljd+u;NjNrP6h7(xg@(e=y zIQFGaf+D zjVxCuIIvuT#1+d$2bS}Uj544|Vmar)a)yyn8Wb_I4F6e9vz%IRVFD{7;{ry8EsP9X zH%wt=V(tXS*v2xI#ehX0;`wb1{~cHvSsgAourkUmn0R4=gM)*k z1B(KPudra^1r}CN)y2xFc)@{X$%2U&96$;jSe7|BIyo<3aARcH2J*Y!0tPn*cNX1A ztc;~B+AORrdMu1A+8{H0Kzcy-vmCtOz{w*V4R z6C4@-JAr)R0E#{)hXs%}>I6`To?v7+#Uck$=itbY>408KUI*Fb3=0?;o-i^zXJzVu z8L~p4SE>;X@mBCq+a8?Cs&_b*Ol}_;T-vK3PAn?bE%u!|S4 z2rXdYhnrXqXH~;lRoHETCFBxlX$WrFD#P6XZR0vDV9|s`f(Nus-T*V0kzv{b7DJ>+ zn1UM2lTo=-QMofniv))SEK?B{H!NV8ie&M8u+c1YppFN18(A_yMIxwVU%)aA)Ww6y zEnryzV}M%bpbj}GIKe>)D-IYLZaXYsWtM^{1eITx85y2~Dgh={hYJoY#}-Vyz{;oq zssuJnVP&?Dc64xDz{qf(k+J*2`U4K&igQ6Ai`fE3hL4O4pI8}vCb1+S41a*+@H?QS zhVa&Vus@&?!^rR#mHQi&`v=5@x)JPtC>zp(KoWt~)~t-W6C7BX99fxjSeZ;%nT!`) zm;q^nGgm298h}gc-;4}@S()=Du~Z=(`~=Ct4-pOqjZ-i(+ymJP4?+-!;Xl;SZcyhF zZ0dhTMkZF~`bjMC@|%(28Iq}w&`rG$G6>C7P(=J#{wfWIQp7F(H8)!Zdn;i9UL4$(Kj~+ z5`Aw`qA!T${{luvUPeYfRwmU+tc+Zc;6SctUPImJ$jI;l>KA0`R|x6nh#UqIWcYu9 zgMoovaDfA;TB;UV=l-y&#eM-JqXZ+P3@h^c425te69o#|>793a*%*x2(;IM#^QJImk-N9kO zg=LWDJ&WN2Mn-K$MjaNhNvwBIdv% zk_c*-FltR=VaM>+cO8hgic5+z^U^{7`nKUPi!P``%gXdm4I=C?!GV!c3*s?mW`xHy zK_1fvdF(4I)0YJoCQVV~T)@Z}$;cSBVdDZ;7A=H+hW`xz5#F4D@W({tA{w6m`#|Lu zLQyZsvj}cKDt8jNJO$Uj=otrG7J$ML%s@Ec7g9L?1j!=||Bh<-4-gk6guzBZMHVnJ zdV+n+qVK?>o5Jva5{ocKFtZd*aA0M&+;Etsz+uBimV5`6Tt-G+P;j#3IIv_fGU|ZB zQIFw2E3*+RvmrQv$Sh!F%wc5AWjQs8l`&uvE29`F87D9@=5E-?QtH4`oDPmt&q=I| zLWoq#9HI+}&b-v}oXosbkiEeh4zqNDld&xCm@9ML9Q?axgv;_IdH*+ zNgx{=85x^dnan4#GM<5i9#TeU1P>j7dnJsF&;Ui2W&s5fM4ACq4nV>fSp_p(1;am3 ze;-o?6GFv*kO~kNR=uIH7cepwGcq=B*vO*fz@m@^4xRUtShyB2GA1xGHnXs@Ftad0 zqN|yOjfI7Uc>yD%0V88O3;QHi#&0aFEKHyxhIx_$3zGvY^CSmH2M3m$Frk|l99V8I zm;f5>0gaHevaC@n&P>ls%mEFGizYZafQCpI{>M8yII=P?xG;g`70XK&o=L2XJ(!^e z9$Nzk3&TrjsDbkLXK<2dWmyF62!oRu%V!ob7A_V}NMv*}GIp~vl|qaJ=~>_a>VF~S zEJjcR+R@PwQsi~BGF5>*&dSIFXKBM(olq92b~(Yw*bVA!vD{*1YH@G`RbQa~7Rx#0 w{#H4ts{+np@ctIdW>9~N -0.8 *. (last y') + and z = up(-. y) + +let hybrid main () = + let der t = 1.0 init 0.0 in + let y = ball (y0, y'0) in + let z = period(0.01) in + present z -> ( + print_float t; + print_string "\t"; + print_float y; + print_newline () + ); () diff --git a/exm/zelus/sin_1_x.ml b/exm/zelus/sin_1_x.ml new file mode 100644 index 0000000..d185ae5 --- /dev/null +++ b/exm/zelus/sin_1_x.ml @@ -0,0 +1,74 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers +(* open Zls *) + +type ('f , 'e , 'd , 'c , 'b , 'a) _sin_1_x = + { mutable major_22 : 'f ; + mutable i_29 : 'e ; + mutable x_28 : 'd ; + mutable result_27 : 'c ; mutable o_26 : 'b ; mutable t0_23 : 'a } + +let sin_1_x (cstate_30:Ztypes.cstate) = + + let sin_1_x_alloc _ = + cstate_30.cmax <- (+) cstate_30.cmax 2 ; + cstate_30.zmax <- (+) cstate_30.zmax 1; + { major_22 = false ; + i_29 = (false:bool) ; + x_28 = { zin = false; zout = 1. } ; + result_27 = (42.:float) ; + o_26 = { pos = 42.; der = 0. } ; t0_23 = { pos = 42.; der = 0. } } in + let sin_1_x_step self ((_time_21:float) , ()) = + ((let (cindex_31:int) = cstate_30.cindex in + let cpos_33 = ref (cindex_31:int) in + let (zindex_32:int) = cstate_30.zindex in + let zpos_34 = ref (zindex_32:int) in + cstate_30.cindex <- (+) cstate_30.cindex 2 ; + cstate_30.zindex <- (+) cstate_30.zindex 1 ; + self.major_22 <- cstate_30.major ; + (if cstate_30.major then + for i_1 = cindex_31 to 1 do Zls.set cstate_30.dvec i_1 0. done + else ((self.o_26.pos <- Zls.get cstate_30.cvec !cpos_33 ; + cpos_33 := (+) !cpos_33 1) ; + (self.t0_23.pos <- Zls.get cstate_30.cvec !cpos_33 ; + cpos_33 := (+) !cpos_33 1))) ; + (let (result_35:(float * float)) = + (if self.i_29 then + self.o_26.pos <- ( *. ) ((~-.) 0.6) + (( ** ) (sin ((/.) 1. ((~-.) 0.6))) 2.)) + ; + self.i_29 <- false ; + self.t0_23.der <- 1. ; + (let (t_25:float) = (-.) self.t0_23.pos 0.6 in + self.o_26.der <- ( *. ) (sin ((/.) 1. t_25)) + ((-.) (sin ((/.) 1. t_25)) + ((/.) (( *. ) 2. + (cos ((/.) 1. t_25))) + t_25)) ; + self.x_28.zout <- self.o_26.pos ; + (begin match self.x_28.zin with + | true -> self.result_27 <- 1. + | _ -> self.result_27 <- 0. end) ; + (self.result_27 , self.o_26.pos)) in + cpos_33 := cindex_31 ; + (if cstate_30.major then + (((Zls.set cstate_30.cvec !cpos_33 self.o_26.pos ; + cpos_33 := (+) !cpos_33 1) ; + (Zls.set cstate_30.cvec !cpos_33 self.t0_23.pos ; + cpos_33 := (+) !cpos_33 1)) ; ((self.x_28.zin <- false))) + else (((self.x_28.zin <- Zls.get_zin cstate_30.zinvec !zpos_34 ; + zpos_34 := (+) !zpos_34 1)) ; + zpos_34 := zindex_32 ; + ((Zls.set cstate_30.zoutvec !zpos_34 self.x_28.zout ; + zpos_34 := (+) !zpos_34 1)) ; + ((Zls.set cstate_30.dvec !cpos_33 self.o_26.der ; + cpos_33 := (+) !cpos_33 1) ; + (Zls.set cstate_30.dvec !cpos_33 self.t0_23.der ; + cpos_33 := (+) !cpos_33 1)))) ; result_35)):float * float) in + + let sin_1_x_reset self = + ((self.i_29 <- true ; self.t0_23.pos <- 0.):unit) in + Node { alloc = sin_1_x_alloc; step = sin_1_x_step ; reset = sin_1_x_reset } diff --git a/exm/zelus/sin_1_x.zci b/exm/zelus/sin_1_x.zci new file mode 100644 index 0000000000000000000000000000000000000000..42b4e831f36d4800c963051a35613c59ba1a365c GIT binary patch literal 9262 zcmZpfx@;ccN?L@rLmg7aY`!p{xZC8xA`0 z#tROP6E-a1U|>*ZaA5e)$ne3z!NJji;Xjhldj|&xCy>bw3nsBJPh2pGg=4`a7A}Vg z8yp-KOag12;K=ab$pHiy{!fI7fy`s91_^@9gDP}j`0oN0a#*nJfs@063lkg`1c$O% zKomGQGK4ueID+hh2{1Bz1Ubsd(c!}S0}cx=99-b8RAO*p(nKBxF!*4;{Xi=M}}B9D+Zg{%3i=?31KZ@v2_5s32YWj zND;z-k_%WO;4CjF3#@zri$5wi8JP=VFJMW5FrefDmJSGO0ZWg=0+t#G2SP4jX@oE! z8ea_MaPbN_s|L=h z#h#lO%HS&M;jAhGg@rUI-NIFupg6f0#Z9HiJVu5FWG+Jys1N{^gNTBHArHBnV8{ij zfcY;UZet;g#mKN3!E8Y=8xhPV1hWmnoPuC>BAD$6W(R`Vg#%^uY5|Kg zM95(Qiyef8aB(w&*@|FpLNI3_m{SqVX$a;l1QXVPb7W+ghY+2OV9r4>7a*7mNeK73 z2yOEb%=HN75(IM*g1H#MT#8_>Loin&n9C5% zT=|WEg*6s+8q|KWWZV3a8}*|mO?nE2*O&xk`K->*$xX>Di*NpgNQmT zVA%<0?SZrQFJL(Y=NyKx7O)&xz*4t>Wsk!GmXinwZ#+ySvc z?E+X{zX}onS3R&6E?8y(3l~%Z%6)*$g|I>O3s@Gy1%(QzoepX~urhKwEMWNwl>jSR z!1523yM6&nC@8~&TlOxBC1}L_|y8FkS3Ul*`4$}M2oy?`Zg0gL4V zR>n85P+(+uh+sZKFdrkBPaGDoGRcCf5U@)fklV5t<>+&Sz847QOVk_&bsea&4>l5^ z@D+mj8o_*nL%S@@AXX-0gz5(iSQOzkJ0rtekjodaGFc(Fr`*(vQj1G-;Gw#JMP&gi zWBURYt_7@&O(0{z@x%meCoN!Q@&$zgR02dX{D*bu7#Y4ot()M)@E=?@BPF4@Vut?? z7eG2394?$6G^%gBH!tA&uok_f;}7l`Yiy?*3C1(!h(S?GWSvaBw` zWN}P;!A&8EUNKBreS}^KP`?D^Z*V+eDm6zamBf^_K*&mC%7S}V5JyO1%32}x%3#V` zBV-jYWx?GWh{^JpvbG4likPzC79T{fGN!COLa!30ELakvS7iaq3`EAoRO$#;3o0WR zRWW5@V{(p+jA{;`LIxcCD9vGb+03F0G74OWf)g8<1&%^63tXCkS>P~(jSwQI6H zVl;w=b&bK`vHb!DM+O&G<~yv6+N_LP3u+l$82&r3GTMXcc}D#O6E8S`Dtb^o&-}z4 zRL@tU)bpOKj1w0y_%Zk+hTGwz>Yfb$C$TaeoW#ngeqq}J2Zsv_CoW*&QPS0$;IP2K zVZ&huhYcGSu!z9bGcvX_{D%okBMG-M{NJ#VMZtkZmXWcA;Xfl|8;Iy&_&?EMB4~`9 zmHFWUh9HIzR>tBq)H4n0j9GVR?r4C8URD?rYk%T8h4DUh`nglX@3W%7-@Si1@C1(LB z1%c-m7(Uo9V8~}EVr6NW#L8?miIs8Yf(sK^<}6?+V<=}yoy5xQG>MgQo5KQzJcf!5 z8(BmgScJks>4MpO5-a0ehXoATU~w}C7LyRLxbh@c#tCpQuq>G1z%mcwe3rQmEVCII zXTjXfGRuKw1|#DPkb_yaJFskFWSkCi<65|B%mfru3h|5+Ze+(#I6c!C4VPpCl$9a#1=GOmFe zw9kQM4WSavcpt^}E}8bqvR_|Ni*nb#D?YlM^U+L8T-eSQ*9r z!BWzbSQ)2GfE6{YjHaM)gO(7ijK&VEjE0PiJ7FQm%4p!g%BaW4xE&OBtc zjN3o~w*y4%0)<`*D`PS|CKj+Vb}QF`ig{MXdPsaMU}c68)VpC zkYQDZLong9nv9I_X-!5(_#6QvBYcXJkue9G3s4$4ypZ+*!uSjXGZRI( zFhr%p0v0hiOA^kKfwSb{EG0M#+^GULhrpd6Fbmv?2D8AOdr)J?0@~LGMM1;@77K?3 ztjyr90ophbxM7Z*yxuEmX@OGhV|eQNhLn;H!qV;nR>r?b;`{El-$FkgkWho=$8Iaj5OB`4hF*2S88G06E zs5C2+6exQ!H!NUeILF9vo|VaO5-aN(RwiLk_BzGLaACtnmNgD6s~7@61qSQFNuVi7 zMutO-3>QJg%svNJ#yt$)V7V!iSXx24_AxSC0?9paU}e0=-~pCvn8Y#x$@So>;2I;t zO;$FyNvy2ZlUU|Kb>7^tkp)x@GdP3wRZU`936j~$$Z%`JMph;c2UaGwYOqwuB$gV7 z1q{iI47WgSH|DLNW-61~f{7PcnYVxzXKa|l$`1a7M;My=a3s?ETXJz z(^*7Vgjqzv0rb;>Mbv?n`KN=U!vwf9Z5A*x{AOhM$I3Q;5-W@1BozN9OeQGIS2F}VIykbj@Lrg};?3g4 z!ZL}KMP(9;8jBAL8w(Q)BRGg#7BGSU6HCk_R@Q|qPM`pOz{tq7VFGxRJjsC}#9`tB z7QO{6oFH0Y0gJB#OBxFo3mXeFE7J@Z<2IDRA_-wQEC^;~WOA6WfR*hUE87(ZP`4id Dwh!nG literal 0 HcmV?d00001 diff --git a/exm/zelus/sin_1_x.zls b/exm/zelus/sin_1_x.zls new file mode 100644 index 0000000..48cc71e --- /dev/null +++ b/exm/zelus/sin_1_x.zls @@ -0,0 +1,7 @@ +let hybrid sin_1_x () = + let der t0 = 1.0 init 0.0 in + let d = 0.6 in + let t = t0 -. d in + let der o = sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t) + init (-. d *. (sin(1.0 /. (-. d))) ** 2.0) in + (present up(o) -> 1.0 else 0.0), o diff --git a/exm/zelus/sincos/sincosz.ml b/exm/zelus/sincos/sincosz.ml new file mode 100644 index 0000000..ebbd09d --- /dev/null +++ b/exm/zelus/sincos/sincosz.ml @@ -0,0 +1,45 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +type ('c , 'b , 'a) _f = + { mutable major_11 : 'c ; mutable sin_13 : 'b ; mutable cos_12 : 'a } + +let f (cstate_14:Ztypes.cstate) = + + let f_alloc _ = + cstate_14.cmax <- (+) cstate_14.cmax 2; + { major_11 = false ; + sin_13 = { pos = 42.; der = 0. } ; cos_12 = { pos = 42.; der = 0. } } in + let f_step self ((_time_10:float) , ()) = + ((let (cindex_15:int) = cstate_14.cindex in + let cpos_17 = ref (cindex_15:int) in + cstate_14.cindex <- (+) cstate_14.cindex 2 ; + self.major_11 <- cstate_14.major ; + (if cstate_14.major then + for i_1 = cindex_15 to 1 do Zls.set cstate_14.dvec i_1 0. done + else ((self.sin_13.pos <- Zls.get cstate_14.cvec !cpos_17 ; + cpos_17 := (+) !cpos_17 1) ; + (self.cos_12.pos <- Zls.get cstate_14.cvec !cpos_17 ; + cpos_17 := (+) !cpos_17 1))) ; + (let (result_19) = + self.cos_12.der <- (~-.) self.sin_13.pos ; + self.sin_13.der <- self.cos_12.pos ; + Bigarray.(Array1.of_array Float64 c_layout + [| self.sin_13.pos; self.cos_12.pos |]) in + cpos_17 := cindex_15 ; + (if cstate_14.major then + (((Zls.set cstate_14.cvec !cpos_17 self.sin_13.pos ; + cpos_17 := (+) !cpos_17 1) ; + (Zls.set cstate_14.cvec !cpos_17 self.cos_12.pos ; + cpos_17 := (+) !cpos_17 1))) + else (((Zls.set cstate_14.dvec !cpos_17 self.sin_13.der ; + cpos_17 := (+) !cpos_17 1) ; + (Zls.set cstate_14.dvec !cpos_17 self.cos_12.der ; + cpos_17 := (+) !cpos_17 1)))) ; result_19))) in + + let f_reset self = + ((self.sin_13.pos <- 0. ; self.cos_12.pos <- 1.):unit) in + Node { alloc = f_alloc; step = f_step ; reset = f_reset } diff --git a/exm/zelus/sincos/sincosz.zci b/exm/zelus/sincos/sincosz.zci new file mode 100644 index 0000000000000000000000000000000000000000..e531551d2cd8905615a83b010c30eb8d1cc752bc GIT binary patch literal 1170 zcmZpfx@;c<14|tP1EVkl1B*KY1B=}S_2A6BW~zW^39h(-rTUk3+Akm)df3dnjV zM~4gR4>&HkaBzXEa&czfg$WB>l#}y|7nGE5n8IS?<_J?A2z8~SKa3F%W5mE1u?`DZ zSU?V8`0oG{^IO2e3u6bt7{SO2SvX+Ap$ib~5QhaUtY~HhEMVb)8yp3aU%(;?(;T^g zg$*tl4weLmn1f>&GA9C=6TN_ibpZ=M$Qp34JD`yW?>K;@U{)n8V3BfIz@p{g2#P)j z2Zsq07C5swGW-V_1S1!)$ShznS-`>x)0Mb@#l~R)iwRs8*hm%|h$fh-IEMu+_7EY5 z1uWpG28AYz>jD-oh)qr`eqa%(s!1$a6D~}e$in~y3u+l0So$WhSi@p29U3YVoEZK) zgG@$vazPd-9Xeb9i!QiukR|EDwgnCj7Zy&0C2-|BUA+ko3mhCa9CmQnuyFwkCp>MJ zG5iM!@WBO28UB|u{GaGB5tQ~=niepGFobWI!m`8NaS}`PB$hc23m5_!!Z&PW5pZDP z^I-Tti6sQ$#);sF0_CD?uuB;JJ7z+{10jPV zM+%}3_6rzd84_5oPGX6G*ze$&=dggq9jajggF8b4!+!@BSBDD@EK3$lyf6WjQyf^9 zUT|PpyI|r4P`Y$ exit 2 +let args = List.rev !modelargs +let () = ignore lift + let m = try match !model with | None -> Format.eprintf "Missing model\n"; exit 2 - | Some "ball" -> Ball.init !modelargs - | Some "vdp" -> Vdp.init !modelargs - | Some "sincos" -> Sincos.init !modelargs - | Some "sqrt" -> Sqrt.init !modelargs + | Some "ball" -> Ball.init args + | Some "vdp" -> Vdp.init args + | Some "sincos" -> Sincos.init args + | Some "sqrt" -> Sqrt.init args + | Some "sin1x" -> Sin1x.init args + | Some "sin1xd" -> Sin1x_der.init args + | Some "ballz" -> lift Ballz.ball + | Some "sincosz" -> lift Sincosz.f | Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2 with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2 let st = if !inplace then (module State.InPlaceSimState : State.SimState) else (module State.FunctionalSimState : State.SimState) +let output = + if !no_print then Output.ignore + else if !speed then Output.print_h + else Output.print (* Output.ignore *) + let sim = if !sundials then let open StatefulSundials in @@ -57,7 +83,7 @@ let sim = let z = if !inplace then InPlace.zsolve else Functional.zsolve in let s = Solver.solver c (d_of_dc z) in let open Sim.Sim(val st) in - run_until_n (Output.print !sample (run m s)) + run_until_n (output !sample (run m s)) else let open StatefulRK45 in let c = if !inplace then InPlace.csolve else Functional.csolve in @@ -66,7 +92,7 @@ let sim = let s = Solver.solver_c c z in let open Sim.Sim(val st) in let n = if !accel then accelerate m s else run m (d_of_dc s) in - run_until_n (Output.print !sample n) + run_until_n (output !sample n) let () = sim !stop !steps ignore diff --git a/src/bin/output.ml b/src/bin/output.ml index c397aa9..e60ea97 100644 --- a/src/bin/output.ml +++ b/src/bin/output.ml @@ -1,7 +1,6 @@ open Hsim.Types open Hsim.Utils -open Common let print_entry y t = let n = Bigarray.Array1.dim y in @@ -13,6 +12,16 @@ let print_entry y t = Format.printf "\n"; flush stdout +let print_entry_h y t h = + let n = Bigarray.Array1.dim y in + let rec loop i = + if i = n then () + else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in + Format.printf "% .10e\t% .10e" t h; + loop 0; + Format.printf "\n"; + flush stdout + let print_sample samples ({ h; u; _ }, now) = let step = h /. (float_of_int samples) in let rec loop i = @@ -20,8 +29,16 @@ let print_sample samples ({ h; u; _ }, now) = else if i = samples then print_entry (u h) (now +. h) else let t = float_of_int i *. step in (print_entry (u t) (now +. t); loop (i+1)) in - if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end - else begin Debug.print "C: "; loop 0 end + if h <= 0.0 then print_entry (u 0.0) now else loop 0 + +let print_sample_h samples ({ h; u; _ }, now) = + let step = h /. (float_of_int samples) in + let rec loop i = + if i > samples then () + else if i = samples then print_entry_h (u h) (now +. h) h + else let t = float_of_int i *. step in + (print_entry_h (u t) (now +. t) h; loop (i+1)) in + if h <= 0.0 then print_entry_h (u 0.0) now h else loop 0 let print_limits { h; _ } = if h <= 0.0 then Format.printf "D: % .10e\n" 0.0 @@ -30,3 +47,18 @@ let print_limits { h; _ } = let print samples n = let DNode m = compose n (compose track (map (print_sample samples))) in DNode { m with reset=fun p -> m.reset (p, ((), ())) } + +let print_h samples n = + let DNode m = compose n (compose track (map (print_sample_h samples))) in + DNode { m with reset=fun p -> m.reset (p, ((), ())) } + +let ignore _ n = + let state = () in + let step () = function + | None -> None, () + | Some _ -> Some (), () in + let reset () () = () in + let i = DNode { state; step; reset } in + let DNode n = compose n i in + DNode { n with reset=fun p -> n.reset (p, ()) } + diff --git a/src/lib/common/debug.ml b/src/lib/common/debug.ml index e34b3cf..8db8349 100644 --- a/src/lib/common/debug.ml +++ b/src/lib/common/debug.ml @@ -2,3 +2,10 @@ let debug = ref false let print s = if !debug then Format.printf "%s\n" s else () + +let print_entry y = + let n = Bigarray.Array1.dim y in + let rec loop i = + if i = n then () + else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in + if !debug then (loop 0; Format.printf "\n"; flush stdout) diff --git a/src/lib/common/ztypes.ml b/src/lib/common/ztypes.ml new file mode 100644 index 0000000..82ee092 --- /dev/null +++ b/src/lib/common/ztypes.ml @@ -0,0 +1,151 @@ +(**************************************************************************) +(* *) +(* Zelus *) +(* A synchronous language for hybrid systems *) +(* http://zelus.di.ens.fr *) +(* *) +(* Marc Pouzet and Timothy Bourke *) +(* *) +(* Copyright 2012 - 2019. All rights reserved. *) +(* *) +(* This file is distributed under the terms of the CeCILL-C licence *) +(* *) +(* Zelus is developed in the INRIA PARKAS team. *) +(* *) +(**************************************************************************) + +(* Type declarations and values that must be linked with *) +(* the generated code *) +type 'a continuous = { mutable pos: 'a; mutable der: 'a } + +type ('a, 'b) zerocrossing = { mutable zin: 'a; mutable zout: 'b } + +type 'a signal = 'a * bool +type zero = bool + +(* a synchronous stream function with type 'a -D-> 'b *) +(* is represented by an OCaml value of type ('a, 'b) node *) +type ('a, 'b) node = + Node: + { alloc : unit -> 's; (* allocate the state *) + step : 's -> 'a -> 'b; (* compute a step *) + reset : 's -> unit; (* reset/inialize the state *) + } -> ('a, 'b) node + +(* the same with a method copy *) +type ('a, 'b) cnode = + Cnode: + { alloc : unit -> 's; (* allocate the state *) + copy : 's -> 's -> unit; (* copy the source into the destination *) + step : 's -> 'a -> 'b; (* compute a step *) + reset : 's -> unit; (* reset/inialize the state *) + } -> ('a, 'b) cnode + +open Bigarray + +type time = float +type cvec = (float, float64_elt, c_layout) Array1.t +type dvec = (float, float64_elt, c_layout) Array1.t +type zinvec = (int32, int32_elt, c_layout) Array1.t +type zoutvec = (float, float64_elt, c_layout) Array1.t + +(* The interface with the ODE solver *) +type cstate = + { mutable dvec : dvec; (* the vector of derivatives *) + mutable cvec : cvec; (* the vector of positions *) + mutable zinvec : zinvec; (* the vector of boolean; true when the + solver has detected a zero-crossing *) + mutable zoutvec : zoutvec; (* the corresponding vector that define + zero-crossings *) + mutable cindex : int; (* the position in the vector of positions *) + mutable zindex : int; (* the position in the vector of zero-crossings *) + mutable cend : int; (* the end of the vector of positions *) + mutable zend : int; (* the end of the zero-crossing vector *) + mutable cmax : int; (* the maximum size of the vector of positions *) + mutable zmax : int; (* the maximum number of zero-crossings *) + mutable horizon : float; (* the next horizon *) + mutable major : bool; (* integration iff [major = false] *) + } + +(* A hybrid node is a node that is parameterised by a continuous state *) +(* all instances points to this global parameter and read/write on it *) +type ('a, 'b) hnode = cstate -> (time * 'a, 'b) node + +type 'b hsimu = + Hsim: + { alloc : unit -> 's; + (* allocate the initial state *) + maxsize : 's -> int * int; + (* returns the max length of the *) + (* cvector and zvector *) + csize : 's -> int; + (* returns the current length of the continuous state vector *) + zsize : 's -> int; + (* returns the current length of the zero-crossing vector *) + step : 's -> cvec -> dvec -> zinvec -> time -> 'b; + (* computes a step *) + derivative : 's -> cvec -> dvec -> zinvec -> zoutvec -> time -> unit; + (* computes the derivative *) + crossings : 's -> cvec -> zinvec -> zoutvec -> time -> unit; + (* computes the zero-crossings *) + reset : 's -> unit; + (* resets the state *) + horizon : 's -> time; + (* gives the next time horizon *) + } -> 'b hsimu + +(* a function with type 'a -C-> 'b, when given to a solver, is *) +(* represented by an OCaml value of type ('a, 'b) hsnode *) +type ('a, 'b) hsnode = + Hnode: + { state : 's; + (* the discrete state *) + zsize : int; + (* the maximum size of the zero-crossing vector *) + csize : int; + (* the maximum size of the continuous state vector (positions) *) + derivative : 's -> 'a -> time -> cvec -> dvec -> unit; + (* computes the derivative *) + crossing : 's -> 'a -> time -> cvec -> zoutvec -> unit; + (* computes the derivative *) + output : 's -> 'a -> cvec -> 'b; + (* computes the zero-crossings *) + setroots : 's -> 'a -> cvec -> zinvec -> unit; + (* returns the zero-crossings *) + majorstep : 's -> time -> cvec -> 'a -> 'b; + (* computes a step *) + reset : 's -> unit; + (* resets the state *) + horizon : 's -> time; + (* gives the next time horizon *) + } -> ('a, 'b) hsnode + + (* An idea suggested by Adrien Guatto, 26/04/2021 *) + (* provide a means to the type for input/outputs of nodes *) + (* express them with GADT to ensure type safety *) + (* type ('a, 'b) node = + | Fun : { step : 'a -> 'b; + typ_arg: 'a typ; + typ_return: 'b typ + } -> ('a, 'b) node + | Node : + { state : 's; step : 's -> 'a -> 'b * 's; + typ_arg: 'a typ; + typ_state : 's typ; + typ_return: 'b typ } -> ('a, 'b) node + +and 'a typ = + | Tunit : unit typ + | Tarrow : 'a typ * 'b typ -> ('a * 'b) typ + | Tint : int -> int typ + | Ttuple : 'a typlist -> 'a typ + | Tnode : 'a typ * 'b typ -> ('a,'b) node typ + +and 'a typlist = + | Tnil : unit typlist + | Tpair : 'a typ * 'b typlist -> ('a * 'b) typlist + +Q1: do it for records? sum types ? How? +Q2: provide a "type_of" function for every introduced type? +*) + diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 87406bc..ac8ed4f 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -8,6 +8,7 @@ module Sim (S : SimState) = include S let step_discrete s step hor fder fzer cget csize zsize jump reset = + Common.Debug.print "SIMU :: DISCRETE :: start"; let ms, ss = get_mstate s, get_sstate s in let i, now, stop = get_input s, get_now s, get_stop s in let o, ms = step ms (i.u now) in @@ -26,9 +27,11 @@ module Sim (S : SimState) = let mode, stop, now = Continuous, i.h, 0.0 in update ms ss (set_running ~mode ~input:i ~stop ~now s) end else set_running ~mode:Continuous s in + Common.Debug.print "SIMU :: DISCRETE :: end"; Utils.dot o, s let step_continuous s step cset fout zset = + Common.Debug.print "SIMU :: CONTINUOUS :: start"; let ms, ss = get_mstate s, get_sstate s in let i, now, stop = get_input s, get_now s, get_stop s in let (h, f, z), ss = step ss stop in @@ -46,6 +49,7 @@ module Sim (S : SimState) = let s = set_running ~mode:Discrete ~now:h s in update (zset ms z) ss s, Discontinuous in let h = h -. now in + Common.Debug.print "SIMU :: CONTINUOUS :: end"; { h; u=fout; c }, s, { h; c; u=fms } (** Simulation of a model with any solver. *) @@ -55,7 +59,7 @@ module Sim (S : SimState) = : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim = let state = get_init m.state s.state in let step_discrete st = - let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget + let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize m.jump s.reset in Some o, s in let step_continuous st = @@ -97,7 +101,7 @@ module Sim (S : SimState) = let _, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in DNode { a with state }) al in Some o, (st, al) in - + let step_continuous (st, al) = let ({ h; _ } as o), st, u = step_continuous st s.step m.body.cset m.body.fout m.body.zset in diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml index 705f20d..23e7b52 100644 --- a/src/lib/hsim/types.ml +++ b/src/lib/hsim/types.ml @@ -69,7 +69,7 @@ type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode (** Consider a node with state copying as a node without state copying. *) let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset } -(** Consider a model without assertions as a model with an empty list of +(** Consider a model without assertions as a model with an empty list of assertions. *) let a_of_h (HNode body) = HNodeA { body; assertions=[] } diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml index 434be45..3f5cdb6 100644 --- a/src/lib/solvers/illinois.ml +++ b/src/lib/solvers/illinois.ml @@ -3,7 +3,9 @@ (* part of the Zelus standard library. *) (* It is implemented with in-place modification of arrays. *) -let debug = ref false +let debug () = + (* false *) + !Common.Debug.debug let printf x = Format.printf x @@ -121,7 +123,7 @@ type t = { let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c = s.t1 <- t; g t1 c f1; (* fill f1, because it is immediately copied into f0 by next_mesh *) - if !debug then (printf "z|---------- init(%.24e, ... ----------@." t; + if debug () then (printf "z|---------- init(%.24e, ... ----------@." t; log_limit s.f1); s.bothf_valid <- false @@ -152,6 +154,7 @@ let num_roots { f0; _ } = Zls.length f0 (* f0/t0 take the previous values of f1/t1, f1/t1 are refreshed by g *) let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c = + Common.Debug.print "ZSOL :: Calling [step]"; (* swap f0 and f1; f0 takes the previous value of f1 *) s.f0 <- f1; s.t0 <- t1; @@ -162,7 +165,7 @@ let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c = g t c s.f1; s.bothf_valid <- true; - if !debug then + if debug () then (printf "z|---------- step(%.24e, %.24e)----------@." s.t0 s.t1; log_limits s.f0 s.f1) @@ -212,7 +215,7 @@ let find ({ g = g; bothf_valid = bothf_valid; dky t_right 0; (* c = dky_0(t_right); update state *) ignore (update_roots calc_zc f_left (get_f_right f_right') roots); - if !debug then + if debug () then (printf "z|---------- stall(%.24e, %.24e) {interval < %.24e !}--@." t_left t_right ttol; @@ -280,20 +283,20 @@ let find ({ g = g; bothf_valid = bothf_valid; match check_interval calc_zc f_left f_mid with | SearchLeft -> - if !debug then printf "z| (%.24e -- %.24e] %.24e@." + if debug () then printf "z| (%.24e -- %.24e] %.24e@." t_left t_mid t_right; let alpha = if i >= 1 then alpha *. 0.5 else alpha in let n_mid = f_mid_from_f_right f_right' in seek (t_left, f_left, n_mid, t_mid, Some f_mid, alpha, i + 1) | SearchRight -> - if !debug then printf "z| %.24e (%.24e -- %.24e]@." + if debug () then printf "z| %.24e (%.24e -- %.24e]@." t_left t_mid t_right; let alpha = if i >= 1 then alpha *. 2.0 else alpha in seek (t_mid, f_mid, f_left, t_right, f_right', alpha, i + 1) | FoundMid -> - if !debug then printf "z| %.24e [%.24e] %.24e@." + if debug () then printf "z| %.24e [%.24e] %.24e@." t_left t_mid t_right; ignore (update_roots calc_zc f_left f_mid roots); let f_tmp = f_mid_from_f_right f_right' in @@ -303,7 +306,7 @@ let find ({ g = g; bothf_valid = bothf_valid; if not bothf_valid then (clear_roots roots; assert false) else begin - if !debug then + if debug () then printf "z|\nz|---------- find(%.24e, %.24e)----------@." t0 t1; match check_interval calc_zc f0 f1 with @@ -314,7 +317,7 @@ let find ({ g = g; bothf_valid = bothf_valid; end | FoundMid -> begin - if !debug then printf "z| zero-crossing at limit (%.24e)@." t1; + if debug () then printf "z| zero-crossing at limit (%.24e)@." t1; ignore (update_roots calc_zc f0 f1 roots); s.bothf_valid <- false; t1 diff --git a/src/lib/solvers/odexx.ml b/src/lib/solvers/odexx.ml index 31637aa..30252be 100644 --- a/src/lib/solvers/odexx.ml +++ b/src/lib/solvers/odexx.ml @@ -51,7 +51,9 @@ module GenericODE (Butcher : BUTCHER_TABLEAU) : STATE_ODE_SOLVER = struct (* {{{1 *) open Bigarray - let debug = ref false (* !Debug.debug *) + let debug () = + false + (* !Common.Debug.debug *) let pow = 1.0 /. float(Butcher.order) @@ -274,7 +276,7 @@ struct (* {{{1 *) "odexx: step size < min step size (\n now=%.24e\n h=%.24e\n< min_step=%.24e)" t h s.min_step); - if !debug then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t; + if debug () then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t; let rec onestep (alreadyfailed: bool) h = @@ -288,11 +290,11 @@ struct (* {{{1 *) let tnew = if finished then max_t else t +. h *. (mA maxK) in mapinto ynew (make_newval y k maxK); f tnew ynew k.(maxK); - if !debug then log_step t y k.(0) tnew ynew k.(maxK); + if debug () then log_step t y k.(0) tnew ynew k.(maxK); let err = h *. calculate_error (abs_tol /. rel_tol) k y ynew in if err > rel_tol then begin - if !debug then Printf.printf "s| error exceeds tolerance\n"; + if debug () then Printf.printf "s| error exceeds tolerance\n"; if h <= hmin then failwith (Printf.sprintf "Error (%e) > relative tolerance (%e) at t=%e" diff --git a/src/lib/solvers/statefulRK45.ml b/src/lib/solvers/statefulRK45.ml index d85658e..997d0ee 100644 --- a/src/lib/solvers/statefulRK45.ml +++ b/src/lib/solvers/statefulRK45.ml @@ -22,6 +22,8 @@ module Functional = { state; vec = init } in let step ({ state ; vec=v } as s) h = + Common.Debug.print "SOLVER STEP"; + Common.Debug.print_entry v; let y_nv = vec v in let h = step state h y_nv in let state = copy state in diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml index cc8c255..f6dc60b 100644 --- a/src/lib/solvers/statefulZ.ml +++ b/src/lib/solvers/statefulZ.ml @@ -15,7 +15,10 @@ module Functional = vec = zmake 0 } in let reset { fzer; init; size } { vec; _ } = - let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in + let fzer t cvec zout = + let zout' = fzer t cvec in blit zout' zout in + Common.Debug.print "ZSolver Reset"; + Common.Debug.print_entry init; { state = initialize size fzer init; vec = if length vec = size then vec else zmake size } in