From b49ea2458736dcdeb877513992ea6a3933c8aa96 Mon Sep 17 00:00:00 2001 From: ricardoV94 Date: Sun, 6 Jul 2025 22:27:47 +0200 Subject: [PATCH 1/3] Vectorize CumOp and simplify infer_shape for vector case --- pytensor/tensor/extra_ops.py | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/pytensor/tensor/extra_ops.py b/pytensor/tensor/extra_ops.py index dc92238010..c27602b7b1 100644 --- a/pytensor/tensor/extra_ops.py +++ b/pytensor/tensor/extra_ops.py @@ -13,6 +13,7 @@ ) from pytensor.graph.basic import Apply, Constant, Variable from pytensor.graph.op import Op +from pytensor.graph.replace import _vectorize_node from pytensor.link.c.op import COp from pytensor.link.c.params_type import ParamsType from pytensor.link.c.type import EnumList, Generic @@ -360,7 +361,7 @@ def grad(self, inputs, output_gradients): ) def infer_shape(self, fgraph, node, shapes): - if self.axis is None: + if self.axis is None and len(shapes[0]) > 1: return [(prod(shapes[0]),)] # Flatten return shapes @@ -473,6 +474,25 @@ def cumprod(x, axis=None): return CumOp(axis=axis, mode="mul")(x) +@_vectorize_node.register(CumOp) +def vectorize_cum_op(op: CumOp, node: Apply, batch_x): + """Vectorize the CumOp to work on a batch of inputs.""" + [original_x] = node.inputs + batch_ndim = batch_x.ndim - original_x.ndim + axis = op.axis + if axis is None and original_x.ndim == 1: + axis = 0 + elif axis is not None: + axis = normalize_axis_index(op.axis, original_x.ndim) + + if axis is None: + # Ravel all unbatched dimensions and perform CumOp on the last axis + batch_x_raveled = [batch_x.flatten(ndim=batch_ndim + 1) for x in batch_x] + return type(op)(axis=-1, mode=op.mode).make_node(batch_x_raveled) + else: + return type(op)(axis=axis + batch_ndim, mode=op.mode).make_node(batch_x) + + def diff(x, n=1, axis=-1): """Calculate the `n`-th order discrete difference along the given `axis`. From 70b1fabc478281882aa7a1fddede47232332edd4 Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Thu, 24 Jul 2025 15:03:35 +0200 Subject: [PATCH 2/3] Simplify `Subtensor.infer_shape` for reversed slices --- pytensor/tensor/subtensor.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/pytensor/tensor/subtensor.py b/pytensor/tensor/subtensor.py index 6c881a0312..c9961e0564 100644 --- a/pytensor/tensor/subtensor.py +++ b/pytensor/tensor/subtensor.py @@ -948,6 +948,9 @@ def perform(self, node, inputs, out_): out[0] = np.asarray(x.__getitem__(cdata)) def infer_shape(self, fgraph, node, shapes): + def _is_constant(const, x): + return isinstance(const, Constant) and const.data.item() == x + xshp = shapes[0] assert len(xshp) == node.inputs[0].ndim outshp = [] @@ -961,10 +964,17 @@ def infer_shape(self, fgraph, node, shapes): # If it is the default (None, None, None) slice, or a variant, # the shape will be xl if ( - (idx.start in [None, 0]) - and (idx.stop in [None, sys.maxsize]) - and (idx.step is None or idx.step == 1) + (idx.start is None or _is_constant(idx.start, 0)) + and (idx.stop is None or _is_constant(idx.stop, sys.maxsize)) + and (idx.step is None or _is_constant(idx.step, 1)) + ): + outshp.append(xl) + elif ( + (idx.start is None) + and (idx.stop is None) + and _is_constant(idx.step, -1) ): + # Reverse slice outshp.append(xl) else: cnf = get_canonical_form_slice(idx, xl)[0] From d049dadf9f1b73cb2dafc9ffcde446fac206cb8e Mon Sep 17 00:00:00 2001 From: ricardoV94 Date: Sun, 6 Jul 2025 22:27:56 +0200 Subject: [PATCH 3/3] Add gallery notebook on VJP --- .../autodiff/vector_jacobian_product.png | Bin 0 -> 36611 bytes .../autodiff/vector_jacobian_product.ipynb | 1071 +++++++++++++++++ scripts/generate_gallery.py | 1 + 3 files changed, 1072 insertions(+) create mode 100644 doc/_thumbnails/autodiff/vector_jacobian_product.png create mode 100644 doc/gallery/autodiff/vector_jacobian_product.ipynb diff --git a/doc/_thumbnails/autodiff/vector_jacobian_product.png b/doc/_thumbnails/autodiff/vector_jacobian_product.png new file mode 100644 index 0000000000000000000000000000000000000000..8b530c77cf9e6e6093959d683250243e93fd66b1 GIT binary patch literal 36611 zcmeFZ_dl0!A2+TNk|HUEtdOiEAw*V12%%wT%SsYil~p7mNkRxA2_ab_J4uodk|ZH} zi|_O3^S$qX;QsaUaCP<3+xvZ<$9TP->o{*|X&&1|%R);*L9t2oxUvof#o7-1?;s5o zezLCdr#1eM@{*#eE)BlCXe@5xzjwN-oN+znWaWCr%*B$z+R@3uQuva&i>0OGB^xK# zv32Egc##Nsk&=t0nX9dnBd4ydgC<WlPTe2RYR)+jEMEiXG$>Js>T5Ncw;%rwv-BWO4yGLK3Wa% zc+X)PetbfxQ?7ickrwB5C+1@7gemW_IUg3)2p)zl_zJgs*>#wL8(&GKfn2dS@pasi zu0xr837FUMk$=(IsU-35wGBIRZ;@B%*4)(~FL)vrXir|i9kt(LD|tow?f<>S|Hk5f z>*9ZVC>5?%^hbG@<)`E);-f7zU13$FA;Kt^V+Rj8w`tF)^FOh zN$Tn}|Kj4J?oP>;yt_Ob@DbL{@hVbmH3m+l>(g2{_*J~@3LD86zXVyx^F-MmYhGP-ena?8A7 z?`b9#i@dWuk*Ysj4&6+;X)1bCoiDlb%dPmRqsiCnEfi0j*!c73Pv_Z@voG`WE3EPo zOG>E)LhvddSb%O*md9w#Q66$h!Sym8sW-8rF?&cw&ZH~wB`>~VVfhCL_a zmCv4y5nyG%TPeADzpU)->S~qg+1c*_ z%*IDBYZ(7=&TE5H9&h_Es;jGWvWGd%Ij0|&ique+xo_wzy?_7p6VW1E+e7-^TwnKV zbaa5~1@~6UH49$s$;rtN5)zK_Gp(gKBeY#tPcN~wbXQwj8*jTeZ$k5K##^Q?E{q#D zZoFB2{=<-~vlQF-)l$5H3NF+x`(&1W$nW2FC5yi^UtR2{3JeVV=2A3XI_6Dpl^4@} zH29eoKPMAy-|x1}U%!4GId!V_RFs>LurODgW^!>c>*;J`6^U3q$J?&byJknbg8nYe zDP_0IPs-e`?Ps9ui)pyQqpBzpT{S;1b)$mw4%?x%K~)o9J@tPrjCbrw+I1p==k59Z zUyCkL7drJGJ$;(_s;8&+sZ)8~N?om#3P*&tS5;Q7m0KD)Kn8>*R57y+i#HjI;=+Xs zWcjnOu;}jX>ew{7Ad59fMNPe5RP z(Rt-LC$ciyyRu0-Zu8&2f8B?(3coo^iAqWaO$|2pTgxx)H1_>_*=w%2w|b(#w$QWA zfO`ks_l79J`9a?mYpncl2EB{hHRGh$haUA~m2~_Tc~Y*Tce$zeTlUKDZ0Fyf_P%`i z^6BE#$ppF3loZ~tB_3PHO6SYxGUqW2`l$zmbw7PZJ4*HY{@x8hnf}8JGT1uwK~dadUtefzk-v~ zzRu20l+i+)wgQ?DF0hC*7a!`z3Eqf9Ed@7($BM=y8H zIJ>qb`B+PDZ&+|JEybHerJ%b_Ru}tgHhC=$rHRHiJxx#7PI;oVGUwpnfDQ1bu8u}- z;m1~tGRo6lbI!{e>guM>&J5R&Q20(C^YyQ(sY&ahJ}=I)C9HP;;lrV4gaR`%1jHpI z1Wf3u{LY*`TY+)CJfM}8?LB&72RnOBckz|*#({X-XOkD6Q(EZ%resSuj2H9z7LReL#k~(aNXXjOA%8U|Cucv+jThZj`vl5+!@|S2U=))u z^`)hyUkaT#Py_H?V!ZF_)vGUjSAE{Sdq;-m#EBCbj#{oohF{8k<*~u{A2{G=Q0Qo8 zWpy-v?4&Oz&05q@5yMy4(0ath#UGcKbD-2WPkZA-e+~|E1#`=?(`RI4P*9+eV#U;3 zFrv6&Zk&q@rQF9krDbFihEwZ@XOghZG4W(nXNKD`0Y%zwv7#1RN4bZGhhO=w%AuI| zRRv8XbNr?8U5NI*RX!h5PJTlAcH@Ba&v!cqR#$)Qk;e+m);U?sD(ymZ;~KWe@}Dc^ zNSs3HFi;bMl#r#mzUD#2k#Q=VZL~#zwr-v9Jn# zJ%WdF?HZaJ6op-zLaLK~)7~}CKBl3uscf!vtYU1ge{!JlL45pCk?8&l7OW%7Ck?oD z{1v(0H8=OO_Rdood(Y^7YHwG{Y~zk?I`8eh3+*#2D@)B&CQg$&t@YOH*V5t#527b8 z7STI(&|)6s66ojtbo*h?dV0#B{ug#vbxvWFvFWxmGbgs$C4AwnG8KI@nVYDWZj+W%ac#aEl^^W>?>b((SGeK z#kQGE)}5iLsfmVit@~W}l`SFFJ8}!xC`q)4U|C_=CAIKNOH1E-zm_GGDl;>ahRW}Q zT|V|X0c8Ly_%_ZC(o8#Lr+h`Jt_({EJtkb8@8bA31z;BkRpx-uy!I5SNMV|_xZ$Yw zh&DTKdm2pXi*s4{^LUC5-V)!fL!XuqUn+M-x# zi;c3SB`+!3u}y*KV8D0*d+uden3=7Sbo}~epfS3oqvMuR&}Q!;Z%)TTR%&{7WsEuB zo;?#Lt1~9^)%$$bC@_3LRYt?94V)V+ovjn*U0PbwHZTYQ?)#qCEa%rexb(3%vWzXz zKFLUtN##M46_Acxu7CN97mr`w_&qY>U)A?z+G*EjYE=EGj@mdixhG~eHZ)hpdIIV# zDp4pU9d_s#Bx-4Ox^AMQGxhRf3s6)@uPF7H;OHeU8po%1DVgO-73M`I5bci)NJ zb5bCSX&)-W%*@Q;@`qaL5j?;f=BOpi`;3)=kFjT4T3Qm-&&Ub}RbBt{$MKPzHw%ge zYCgURP`PkH&|bnoTjFA4|6#9rQIs8A6oH(mjewV#D^tJf-TtHrZ*kxb=KlTr_jrD` zx5n8U_k=M=OtcmD67mI}$B!MOb}9OLXKaPD`Eyuks2O0?wr$&f4GrzJGZ|6m^FuW* zS^87FgM~#=oCP?j`r9{7jDz;+(?O@wH2d04n{-bVPxGoOa))bKspOvX^<~ zOuTkRsK8n&o<|i+1VFmM+?Or%D0TxidNy{ru0V>v{lThx9rgqjVKy$#{V_r@U=RSeKcR;S8+t zvary6b>**6X_q3Tv|@k5G26l}j=c*6J#vR%YaF<0PkjfK)<}EtZoP$bU*&p~hU90@ z=m2Z|TZe>xTWN$;-?)E&`|?o2%3j7>l&Dy_XEn0h1usXXeB={OXjWqI-}k7rYhEV& zWNcT+E;{zPFY#_xm-@@~mJU8^v%{7Rz~Uwa2dzjnMjiG0afW=FE|`MXbjz_FtgK2N z9*2D9t2pE?NruNgEh^ecaQ(zo!YUh9p$eu|X?FYhx>Te7K67(=cG*`B2P6T3fd^k* zoe{XhA-_4h-C+8&%cRb9%{AfRB$sb%C{WYn-3k{SrM77@sQ_J*4vMBU?k?*-loEY! zaq>LM_KU93;-o9?IjRv)Uix_5*Wy>z(Qe>&Tf1;*A@|4NDUHdc#LUc%C{DP>DzY70 zTCU%@b0* zm-arIwVywK9*^x^rthuLnU?eVb#U>OF>c0N{IP-Y zKWWT$!q~o*{mNe+Z=aiNX`+`e*mGDiteA>v%EQ*yR^ig60~CZE95`?Q8%Jfy!&8qt z{LRRS5DH6_&LhxE;n8}IkZOWHou%mVI_N-ZT3WQ{|E>&nezF6WM#Wi2y&0I8ybS>0 zv!)*INA;;jyQ3TXPM0rkLs@M3`n4pBI`V-!pRTU%J702BaZ8P(c!H$g(s=t1z5Au} z#aG?eJ!N~6@bh0kH5FwV^+-=oS2%TQ>sLp~wbv<9)6(`!O9#Jx{W@v6q_kmhs&Awt z@1`=h`$RSK#pz!y!W>F?Ew(>raMEPY$065-r3nvksFKMB!82#hyaAn=sG#O}GwS9v zr0$(g)*Tn<6y{S&Ny*06wxV%gIbCG^ONTsbG!Re>r?LgU9*@CY1YQENqEvKz{CHcL z`+9x7DrohqZbu-Tkb7QqwrhQWj?rVP>Xt2sVV`%smut=Mq>BFZ_{#d%g96??#Hv0kkNaZLf*nZ{D z*NUcC$%~8Aw(h0I&zZJEUmtZxJDr-En*TfIn{6ocy12NobX#1vV=a&};Za3JhaSmX z*^FB^&c4!kRlBzT;_a|7OnycG#N+4B9Zv_MP_9d8M!nX?CMWnIKmT%Zke&Y<_|`h3=`Ii3`=2F z*|I%kJ=puz(z2(KYxMN=2p1$(5Dej6QxkfHd`=UmLjeXE!(i^?!yXzM+8{90dhXo0 zO7N}_>6)N}5)w(}<*VtdWnKq9fA~NFe1mfSwfGAC<;#~#-v7JzGn4AtoD6S1d|)7$ ze|f%Ez9xM4F}Jtfh5DJKAp)mO^mvZ-wVsq$0cinSl4U>HV%@>3s-sw-nFd^Ry=zqD zceBuaXdu9m^M8G8#90Qa{&@Vwd+A$xrGU|rk;lC>HZ~sTiwrxwy#$0i{P+ zE`-PF>Z90LB@|~vk<;3~#``r+o~#<2QUJz1pYyn&px|j~&d|j4)x&`5p!a?i6$;AS z;kEJKA2bhcxupy&*Ee&urKe|eUPp=GY3+D@^t?a4FgzpW`ilYlFJWkj zAVK-nh2{{6_q87LkNnfxP`c1I|Mi`nJ1c(t;4$`^s7|`#EC7Mg%67BMNKi*6)WHn`G5Z(wXpXdUPQ11hJS z+#M))&QiCDdczW$MB1t5Bln;5A*15Zbg+Sb&0IxQc+)>I5?tWt8H%;lKCFVuEFd88 z;^j+!Eq?6%(!DRL&+jMTB}P37CDhy7`=snHm+W7E{^-vJkPCf#YisPHNh9?F9X za?*NH&eGCSrRX!WA1(u0T|+|yUO<=*iWo-G`D8Dr(;m6Zm%9h1JU}1F06vI`x!yYT zRxrC{@~w=c0>i!Bh!Ey|1auPV2kcYwZz=8r6H1ydw5SZWb;8k-kEAZH1=TEc?5@ne z(9RH6tLLY-?pj1sHt|;^WDa^q{8&-eJl)V39-_ z8XXVw7MWUD_-99QcXpjt>CFE)B#wUgucg?C zCxPzcvrjHBb=YsOs;c@Yok&UDG!zLU1dphXz%$R6OToIAA?@t>^R>uNeko7Biq)xP zLOs?n{w-=OR8Am{SC@yk6D^e4#7jGB^G?=-sI%$ms`+Du-9K19qbxql&aMJCO-x*~ zImo;?NdMy<_QSNQ5zebC-i9UaTR`@qC<6nXL4_ne0P3^XOr{A@nC3oRm{fgdbpUt; zEP4Ej@7$GFUJGLAuw^q@PHQMYIGk6OJ%sk{>-X$kVkZpQ((@Q*h6sYWxnk-X8fJEO zYDFDCf7aHX;EM!geXUatGPM%HpPl~n=^OBhoPvTt6f$*9O)(4R@IOwgvremS#ncMZ zFZ%k-N>*2vh$v?4vtU^=+j-%nEUg3+t(m<&9f#NWn$7I8H_dDBn1QcC7J1QOZ*$>- zUz?o@R7do`iAc%FuFduOa?FkO9y@z>XOGX41EDAwV}pV%x=bp-0$-Gq88`e5eG&cvWL(e;AZ9TYkbp^sbS({LxpzCl&*7NTX3smG9cqs&i0KMx|TN@X9 z*zOp0?ilq93CGN~U}$zkIU@5*ezEoH>&MUjPZl8UBv#OUbs7(SSPkUHWZ9`Hc7z^< z__6&mpA2x3t?dSkD*zpky^WJo?c$>J_X(-Vg_7K>-)F9dNZg52Pix(Kz!*|tNJt3G z155$rC=fTynUc5n-#fpi^YlO$Gcq$TLZ`&_RSZr+@>Yo80hb3T`A>Iquua?MS8=VT z=4i^h3@v{C`~t{61`>C=LEX&GpI~pmnls zroOXFCJ-y#Vb0l&&aZ+ZyS+Hs@eLF=FgB>uo0D5;6c}(#*GxrWF?5A5`^ZC$guG@i z!AewIfOE#N%Xk%hnt%Gx4R_gp=1EpoC|5Df?22HpTO2?snhXS`?5S1Y;s>S4NK`T$`jTe zN@i8NTDUXz2;Blx5sC>)GUO%t-9C?J)vTEfEMc|E#+ z_TAP=RQ&Mf!HbhWR002CZ~+Em`xC&Gm9>dte6n#LfoNf_?=a`+(hWm_Dbp3|GSSOE z7mK@5v;vz4Pw{a=h=_N}o)^N!fY{uRj^+YoIrMmTYPTxm{rmUHB7)YeE3h6ogp@;Y z6vzV&fo%loXT`PK_8`d^<(9(h1e}q(Wft^@1153E?P*1?1#Mo(bFGFB}5b+ zMs10_`dAAC0|VfZwGc$Grb|kmmB^caA zxP%5dhU|9YXtnv?X&$Vao)(8QL}6-axgJd&`x+QPXO^x$_LO8GBwRqsyP9dJ2D(mR z>1}q}CMNV4KXhPIQ5ujIXQ`6phR3sO1MNGvz{rA2w=gxBlbgF1AJJbM28I5x-raDX z4N3gxfiO_bMOi|@`8<9cxj!>=`PuX51l&V165T;pXCX!u3bSvPTkLF4i4OSB|9|0?6%iR2q^AN`{MW_}|0G8rYuBATKT~OXpEmBkYq~bJ*d^h=U z3*{QD-3A374&UYdXos&dH!mfA)Y8_@N$0f%Six?=m%GyZo3CFVX`J2 z2N3tKGGmi*qXV>2v(6-1%luOLswegy)XMXiSP1U06&NTX(1ES;Fm6Dzg!GwSw#5+K z<@s~Xv%K3?v^%J)?`DaRVfs;JL2TqBoR8lv%Z*Y=K_(A9x`l@F?s1-bp+fW~kFa{T zZ{Pmm^XF7HB}Z|V2@p#TC^%>c@F$O_wKmVVL+=MX35MN81Y_tvdW%QoZSySCZ&jZM zFl-2d%&KB(IdH!qoEA_qs4AgLOe?#+|5r;gX(WJcax3%w!mJ!EdB6j)#;+Auu9yuVn-(2b-Ty+s4~(9tcCx9h5j%$78YQ5lNFaLBMKD2J~gC8>~_+GhG)2= z8_lqWh*6H_3~v&^7*i6=m|JR`?J~t1%nhqTvfJyf^p0HHW5nk{|L~IIK;;JgfC94o zsnNpac(kMfG#t6ycYOPnn%}?w9xBhvIVSAP&3B3z=#0R248&h||Wy)ju+t(f@uCI40Tqpt9+$hKA3G%zrr;?Cm>w^Yl> zC=~4mw3f(fV4;9tbRG1v&kdiPInXgV@(~CX!o~+E;twCPV4HC#Gz%JOaEIS?Id80| zx9x$FW{JhAnKluX z?X_M%b0#On7HqsWjN=$e(F|SCp3oNDR%$9$_;E_eQK+d5l-Dgo1(o7L0bcnM3Y2NLAAN2j44neysdm0@KAzi2~;*mL93U{I*@9v(Xh=))eyQUVow+&&D;AP9XL-u4St zx>0HU%*&^Wj>9YO?0l#bq1T|Stt(-$hj9a-nPv4QK_5pcwm?oAh~i)NQ%e;Yidg20 zKr-yUKYvx5-a6OT;8xPj<2-9o?SMca8GQb%PXyD}xt8IXhQxg^!$E~AmX~D#IMkPK6~~S2m^2)o0Ri4!Obs+f4NG7RiVK8Uy9K_wUaG%0l=8H`UP4XfCNb z8(;4O%@FOv8A5~Cd|%Lx9Xkk}!Df-h%_)LM)H~d1#KnbK8NfpyNQ$tuWO^E8FV%vXF@Rv9S{_hFEkOdzG|}v z8?#|pDz*tl2w<2no$VnhPdqsrxnyFS*0O5;(e>0ngaYO0xEThCKaVQ5C|a?9@7QUH zv%g)Xi4hC01{ZQtC|q#%a*bo5)49WlXh36zvN}NJ9VKefTEIxT20gOjDAS=Ivco#1 zhURX|^A~};p%6l`H^;kUi()+A+}gSa^}E2bTM+IJB$%is;YkHakSSh{+<|MtPQS7F z5NNntwFH%rq#=%#_4M|xzUz%th;`!zatHT-aDaY-oDlvpK#`&o)ntV&EfHY9;D7qy zHrn-+ehQHC6qr;nJ?5yq6sBS^sB(}8MP+0{&=!eCd(yB&X!c+YM5brY)QjY5`Y!^X zGFhTiK&o}%<)G}JS}@H&2ipOGsD%`3@i@&r=(I#qfJtgk^Nq^S>_k^pJ%Pt){b;Qd z6Bn`eOdD%oL>`ZruB%~>Z7Sg9m9;Z_lNzViGcVk%yq$JE$n@@nW~aFi1$r=}&{TCq zxSK69y}S3bghCZo0ERkx^r-V_=RwyiSE>eF%&?kP$LUueKk6TE5X^p@o9ohL6`G|# zoZ8$4&<15N(0&_Kw}QP)w7~x0;EFtq3Mwk=A%_ErmATe^w##p*sU76v$bU{%*cFPZ z_pM^DDbWnmTA`Ehp547yl@W?Ws__r~OuaUVb9Mquv_N#|_-tq&FqVQgGK&!H2zEWH zsY!otOxZyUyOLM=si4QxC~U5Egn1Ub~+LneQ;p+ z#{^gYI(_Y>N_}i`_{iTXZbeQ`nwgvT7rcAsB)RbIF&wYu@@~;c614$|;#e6E!p`Ph zC{B(W>x;23n)6H#8Zi_}^S3WGrjKn(y|+C>^`8U|1!OQh+hN(>L+Z@N8`lt$2nbl8 za!)ye2P_P>v#%NFsu=E6tOCYuRA|^&{!GNjIwe37L1Ir=-UpUAaia15l^!sBINrb_ z=wOKR;4k=3xWb_52u+z>{n^n!!QZXJqgshT3g9S~j7fxo1~N?j6aTuhcWof9E(|hY zt&a`ksJKhr{nMu@75EoOS^&^u3mp}hK$Mba!E27`!M{QaR=+^Di6bx062M9G%3~^3 z=X2t0aLhk>`9|xIGC~IeQ8#+Wj`2la=MUyxSa>2qt{v9}HV0CvcI7TulY<9+3rs^2 zQp1l4d4&Owa3O#gP%Jep)xusM*4jI}mUNgcq8n4*c|mQ590$n&s>EylFl$MyBdBT( z2~Ac9)cINq?%cB3ch}l>0KxsW-zVHvqbUe$(b>x@%MyB2itWy>gQqgvhe^FA3J6KY zA%exDs-q?yd~gGo(>)^t10{)A6;wyj9m?p)22B~)xMB-5(4{~HaA_!k*mP#d_>`3C z!Lfk11{^|63z$Y|Lo}FW82iaBH?aryc`Hc7HmSAWgNzPL91su?aX3IjE^eytu2DlA z{moBFtwVqf5&Hu~yr7JO)NmckN~0vP#?13J&k0~FU_eq+=F zxeV5WNC^MnefMAH5pBo4w5K-y(;558#QoLW307aeekI-wXfo;zgv5cDRXW^-j@k$?8^DEF3uW|{pKPRs~i?zL$2w35)g@ge&3xr9G2%`tjl ze(3Xf7f>e{G%n3lnlMw*ZGpAX3@!7l_rA9}5#g`Smp7!;YbyFei8`DX`|F+6RO(%$ z-v>J;Q{`gnEpBu>TG#bxiR9+6hL#+=72n&VF8;DVpvnwyKF|zg=ve649B2p02|$9VVPL3R-M_s7r*?vWqZ> z&mFnvQup3MES67PD)R5@kn3xWxD(3S%ChX%b^05WBpmR2H^BdvJCIiya2f50vG=%i zOV&Ve&R)1JX2f2iIK68P(xsmgF{5HJ_b#M=Eqi`X$xS4>>yy|gBuVV@#Vq9AOV4EK z^I9|t>z@@#x>FisV^iz-%FrUa*@Bgh%3U^IgUi#*EYZR`cVj?hD!un-i>>KOf4oaB z@)hgO`PPZWyi2_+eC1AOW`bn%_dXV%_ZFE*@D`|mej-E3F8-^qYZOhmH}60Py^~1& z;1uuS$fmm*x(5cT(MLZ_r4Pb?!ii&tlaj@CM1Y@jO0*!=49PKHOmj-`N~= zojP{ym8zW1W~H>&ivPkn@X1*@w%+ppUHfW`>ISF%;)nJjSI=3h5F0^T@rN32q8x zj`~3_w~#}cQ&tvN`0;x`RP6&4`@-)w?e@ANp&D^-)n@RcU4EFq{X%MiuJ4r&``K`p zE$lus)aSC?PD!$OB-`)x5!nVkz+rEZw}C(uQkx{7(PVY%T-WGjaC9VEkdj4F3_osh z>O#8CgmY1yaal#@3AIB{NBQWpX}8e^M68s;HbGw>Xm%;0$v(7N9_VXM)&p;;w!VhY zjVHn&N=V>>0`sWeK$WZN<4?aU%%LYZ_Fa>Q5gi9>N;vPLl}2it-C2>`s=7oVMv`Df z4u<&uCMFyNYZZ5wbt^rob$`WYz!%J|tSwZ(Em<+~Mf$Rtpj9rfY9&Y$RB3P%SXAI% zw;{ah3he9KKHqlo(~- zkh|0B1}7r0drw@MOs&<~(1AG6h7Nb<k=$!v3jLba@DH^oIQoZB}7bCa5f%XavC|g-i>o*4|5}|~sCy}A+u2aBi zMMMcHsfn)MmF*-~tf$vj5NPjL@#d}(PcfebGo~70DBqPqU+1OS(X^E-Sf21j$q(-I z^cV+@I22suW2N_zPtg_m?koiy6eXY5I;4?*BLzuo@P2rP1FP_G0!MaP0aFt70ERiC z>;Dw2*R~yTZtM%*Od-`0aVlC34$G7ZUzk`fOB=N;O@@VqnS+{OEJ^b3ef;4NXbH>p z3Cp)Wn(h6Ed#wx%isVWrEAD1JH`BXHSJ_Fh++f}QK?oJvsZ(G znMnfKePIB)7B}bvUOY~mcho`4=FO1z;2}lP|JF&Uzdzi1B6OK+Pjl|j%+>oLBcez{ zS-*dz05=QMPE3%*#Kckd?T1rB$M%OFJ@wH+Y`wylm*FQj->+t(4S3vQ^8qOolhws` z^n1|gH-5uLSXRI@_*UoJ3Tcg)i5MJwACK@~T)ax8dkC>AvHtFck}L`gYOe!LG>U)j zibiM<1`4D+qW=-%45<|timYD%a0o;w;6L(D*f{VYfGLqv zW<;lG&DB?WX3}vYMR&8{e&?T)jAFHymGyHYw*Nn<7cmL()c061Wv_YcVcAGDwGj5) zuFg)*U~c#<$ZD90#eC@O%vfFqeS%=IH?8+QNQ{LA<>AAJv4{aNYTzlc0?bd(%_T#) z;P3Y`+_Z``g`t_3T=Dl^r~e6I$j2e6?6-A~=3~qc0xRg0fP`4SP~C`sV)qeXkl-_1 z)))JNh|-W=)gOPjUA!L|zCR@(dBiHuqt#rv+69Yf9l|2T-OrZ`bNKGwjBoDmF=w?F zJQ2zt9QMwFxj0$-zmzWnGL(o5QoOyPvHw@^Y zAE64sV}VYo3DY+7{!3g=5d-$KQ(fB<>>yM@N3gNC_wO3*(LJ$^jw%QGyf#zR-O_OM z{wAw_>iNwR&dHv*w8(nh>BU!Z-}`F8$<#z}NYEi=l%Ki6&wCSWW3)D~B7nk>$$W5P z&=A5XVN_r=K+#~+10X{SfyIH)z0PJw!N5J;&UFg0@y*X4bVhZRcn=7gTyG!lo<0w4 z3R=+hZpR>4`7dSpR2iXdk(6IhRgOk-@5tvY17h7_+eMxFBk@_WucF@pu@XYCIy#D) zcI&rS9YC#l;+8(*?ki;xeZHtpb2rUeAa86CRF~)WbuYU{ans81HG!1Wta>TD%l_$? zH(|CSHgAfZ4dVz3kByDZ5oF!bVd|f~AEg5Ry?*4#-kqvxtp;M53F>@(b4!TqN#r|e zr0DQ7(L%~Z$$xQ1NXmOr90{6l6YfYXYVk9{)F$a(skrw^cTOQI3j2s{vy)wS{oTLI zZYp&WG3uz}!0u{?GYubc&`@HE$%+XM4nAh*X(Yl7Q?tw#{tF|eS?iW;;A^PI_)1WbaE zA5ZK9j?8?r7$aYATq1HyC2P+iY8rtUffoPMX=H|O0sT6qD zO@H%_S4+9q?wF|cC8Et?!QZ@SUm+44JgOv;hD8ar2(68@^CDe|xjW z)zy^~*?Nvlw(s9z_o*gHWbO<0Sik?#+g0NduCuX&)%ul7(WBfA z_X>>vJ1gQdo0pJ?c$~>nXFEwQe;IWH&HwgyuD7;$tSkiUmYg1XvvFUEmRO2LLK>4_ z0gozNtbZ9|5<9EKA{(z^b@~RQ=f=_G3YA2egefAMj<4g7M^LIp5PFUz4bU_pzr!=zy`Q;I zOSmO#w?Fr$Qqg`3`w|`;<_^aFmFhkkm&cGlQDKnU$C3oFNmHnARPoHYV@Pb$A&mj< z6sBAbpJ&Jow_Z{(xPTNWQDGtItjve4?$}p4P4pDlV8m`FqC8Sw2*8j%VPv!mvTmX0 zm#z$^brS4`#UvQh`}fZkpu2y<3yhMzTsoz7XF=>Ic3+c=D^ow zBkuUH$h5A9eX{2oVgQI%C8wqyJ$paOBm}LV90K!MmU8SVxdmqnZ7ABS?&vM5QCwb% z)hYrZHgM*y!1Ihsbn~F(jF523-K_Ru$)?93M*u6FCi2M6qH$dO zQMU!UpuCIUFwooQ&*ZWY8YR)J-i7y9u5HS{opRLpqsT#ICLTO|xRwG=J)CLBtDEXd z(`962E<$%9<^{T+k&zMaCGX^>AMI~*%xlZW$Dj}qbq`Y3u3ftb3Y1HzMcNUUhK8Ch zU9W**!Rf0BjgcD9LeOT?f!miR`N2d5FuSd^H(iP<#GIU>upjD$_ zB3OQ%lLlGfO04kzbTJNt3l8Eez@T6=B*_DwiQ`{;P85~Z)tDM(aCt%`<5~+Ej2Co)vKG(ZAfMks|*W;WJR&fGP`zLLJV^p;_SuRh(M$PTA9SxHHHWPR`AFWM`dk9deCzzfVVHSWGuaAc%Zy-Qg}#mnEiG?!a=BaXpzQFhN!&K-(XZuVt_tHeA z{yT7&<5T1%=UaQ}ZsS01+$0A1Ki^lJMMJN~NEKWkoeyt}m_6cLen?J2J;5NS=xm77 zEF4R=f0w`WUk?AvpIb6<`)&)ai*8f_spDrlpY7!0u$3K7u?mi}^mt7*QTMR2v7uHP z>D8#Xy9w~+?)3dZPMIRk4DbEfai9KBNp$`*;YG+(l2A)IP8*Q4vI5j8MEv?Et|jq9 zu0@0f`8LJ0)|`bX3=wU>`SvtJ(ao%mg|Yd6Q0Fo?oNPnz9r%MMHzEHGS!ZmU%|XaW z?iAFhxY!4k#C=(0pOR76GTXQQ z22E5ov;4!fb_(fONsS%^S(#Dbz1ypllve zX|!VPkpSM=FiSqgH(%b18gb247N`kCRl*Fy*%??T=(?`b-RrsqIH-?4_C^j65~h{* zp6Xi5wAcW+)J;}ZsJm(`^@v_0u#WG?wk%94*OQ)FGe)Bz5J#Nkw-CdG(+5iWL?@lQ z_{jAmB-=x`c{9mS+h^cXU=5RcgN!EN9@OVgb_YAVx-J4p@iPTSNKvBSSct_W&E}J7 z=J*)`Ow|sKjx{)lLfXyK zHalWo^Dwdv8iq0*t$~ueOhc7&! z(?iF#H}Lpb@!fVXdhfN%MTYa0%w-g?NuWByYJh%+Bu!qoqK7OyxoI3?Af0zxv)6CvlZniejE&V|T z5`0jKNv@WOc0Je(Y9s}PerCQ+qcH*juos}4Jd4vvsgVM&BS&b#m&ib~u(EP7(5Mw0 z7b|*9#B07t&4kdPQ9E#45lBDT%>+SW7Fr^EcHwG(}ca$c^JvjSER`RWd1t{ zikgb`uXXzL=^tM3Mx5uz&Ae8ZdV@gm&Sf1%>+~Ee4Z_LPl8G8tI3uJI2W(^>@8?)0 zOcqO8B~&dQ8@9G!zPt0bb>Y{WkR|{;{$)#I)TtmPR3OH>7E$Evq*3+t@_*#d$le~c zUCSMAbMp00hk??$o~kIp(*ka?eC3|m!S<<}B&W1yUkz@&lr=UpWho;NMSCoay&9p$ zOq3riIluDNrSj}qQJUg=!M8^8ted~%)F0RtiQ{6cA*J+RAnLIefAEU@*2B1Jxa`|@ zn&W0dStUiWW{^e+PWo3F?Ioh0N!{7f@qc+Jfmv5V#;b5@TT3A7nGsjNVTgS+$@8K1+Wiw&zDyJ03=)?~lSrlS{x6KJ!bTF08RCj(?!UcX zKQFqY@kZfC(e_j#>(+efkCckqJBQq`+rzlZ>?DAPyLQ|jG24KDr9~c%$ zf@5>~bYV87`a$sfXP$M?DWB*Jty}p`fsscb=?-K#fI*BSJc&oL9ucFj>tEz^KF3sG z6`^JtdyXE2SAoj>r??xZ7D7lUF#nvrqa&Ux@#l-j0W1NQP^2-GQgm>^3R(JN3&$n; z>TT`pDgm>Ho%3546t$U9!})`6Ag6}_N)x3P4P~XD?tmuAKZS;Ff^Gk0be8&~0#rai z0-WMgDYQz_*^AOQvB!6E&&j_#0iVdT3h)pK5_^N(m-H*Gw!Z?Ki2_X%ivsmB+G07M z4vphzGE8oJmWU2}pEsUg@gxbDlMXK`Zcdc143|?CY%}F!+yF%YzRc6M-w%iCAos&A z4tW0jxgC9qJ$@Ja0w_OG%uWP1xKDPYz{W~{ed(Ml~L%}J*_kszo zA|L=mgPN+&cPhNj7-SuTM)5c~xj&}RA%^YEm0oXV4z$@*9eYRpwfL2xlo?+g;-%ZP zDeA|8#tkak5MEGM5w$|>iT)-Mf|%6JIN5Tg#{)%wAu%Dk-r@z)m2jMqRL2{yqw-5! z+v~5j0e6XJjsPrN7c@B`SF|1a{ksy8OJ{OC8RyB5 z-C#Tx%?y2!&6?IYa#5Uci z?a%q_2JI2A2IGHLJK-3VYFm)AIjq*N!@dR)3scb@cxnNtGkK~A@>W_*0Z7<<<{m$5 z`g0mG08a2L*O#3vGF5gHXIWn$ZoN?Q-ZO!^|6C=LSZb>+{2F}L=avD2w(jKUq<`IC z`IhNO4<;9>DJaNNzAIj;vbj4SOkIby#F@}6@5X#caR_HmsDi30-o2wEnmHaX@LAp) zv=Pl6yH4uLA2z^w_NN?NpEl!GkxPeIh89fhX`DmrxY6;SprMqg_3d5E0jTFk>e=62 zVFkKFq8`sD$(q%26NXzm1bX?MRN|HQ^X1|!&jfBfx;B{b+z4JAa&WHyyn(aQ4 zrQZ*9Fh3=~ddaZNYbQBoH6gz`C69uQ(@f;y0KIrN2u}Ntsh)IJ9Kw$w!px6X(IZ!e z^a2pkIF5pEDPL^Mevy;2Cj7WR_8gnMkJm}pOg`B|dPSKB7WAX{nbhn_@VbEsz#E|8 zlBYhDm-{w5_8D}eR57ekS`@TOht%)swBw{V1$;ag28>nrEJwM1A1%Ec#<2c}?pp&`q>_Y{l4 zsp7ZIKC=gaOH6U9Ogiq4#vtaIe!fV4uegR15979jb$j(1-6|uio`k~mEt&rwLWq~p z_H`Pq&qB?i@QJfDm94VE00$@rdR+$&MsjySyM3rduBU) zl?jPNtS%547Ev}rJx{JqpCm#X;0WXkq&RT+cMTkN zWCO5%z+5oHh%uq)VPdv!DLadcvm9_eCz80Tkfu*(NW7nU$hHg zQAoh3w;#TKr3iaY-svvc(d(PUc zpzBBkl4sK3+R+_Qcgr8QW^4*YlCM{#!dS~t(ye}5=af} z^^aRZhCuqB9%s|q*c+Z(7C(Ms``Uv26K%lvTQ2KkbfZ8U|BLB=BO3q78ujU&safj!7k_dPKcKato#ZIr9sP34n*g&rwut9v;PcU})2rX@pcT6!7-3 zv&e^|en@*i{zZwGB z62IoPJW*?JcZ7$r0@;@`?C5A^3#ZlR@VbN^`e)?jp8*iFIelr07gf8+vaZdAb=z5ou!_EabmhH@);1 zY!ugbO_`gVg$7?jS{Khj`Y-940}hK_KsvswE6{6h@o~6&KYVg(Ik{uCEQJ;%Yr?N z#sk#@q7?YVh&&a@Ffk?b@*=Z5OmB!`)3SsHj?j%N+^pl%Wvb!u&U_(c`MNGL0?@c- zive>oBXT?4d-BHtu~6{+^4#dU_A{Vt zfNTa29CFSaKi9|Z|FwI=wt$P5OqC__IQ0Z0VDFMb4``jsUE8%{r23l1Op=j~ij)+b z$~1q-q|u5x2|BOL$w2lJ(lXWO%|iH3bUg7t@!y9`Mqc=18E0jVDb6I?vS4N4 zNY;Z#kKO|O{6}5hOPoRb16;>XmBGkRLs3Z1Wm&O4NA3IX;0^2x^qIT5TPmM}ks(#` zLVi!AD)hMR!@57&_b}07e~9kiPnsrTS$MDI$>?&*uv$F&2S6*Z$;v!Ha_jm8JS+pb z{HbPn`r>XJd?!am^3NF++zYo5J1c%neUDXc>!Ll7|5N@h6$v+n_@$UQWr;1j{p}(e z2S@-J&)h>N!(HRTaH0x-Vl@?L%|NKtAl}|?j>sRzi8(^l&UFUrDK$0=Q%! z|Lec}G?ANod&6*^*<|@U^DCDDS~#LuSRcB&2RH-H;((ckem z+SDqV$_DuXw_CFCGZwPgGY8>I>W~lG3#1*8_9Z!g91KPqhLQ&sWAFDA9W9uJQn%~e zmGUFlgh}?y9Zz!}jl^+xU;$@r0PrLa?)vnJ5?W4tq_~DTMg%1g8xq73M5v|k1RJgWSB6rvr1L2DhmV}w){(!X zUIfBkS}P`K`Y;Z#i9FDy$Pgm}#0zcEjI|Yk zckbLw%A%TJ#4$kfLB}a5t{-`Gu2%n}=c0BkvOGXv1O;LGOil|#f$TwLwA|VO7ZWL{ zx*C0crq4f$bnM7g-HTIpla9l4z_MmP2vp9Foy+8w*&YJfIX8R%p$H4@URF?Fm_W#{XOP>aXiQXT8J@KX9nymU0rm8ba%;E5GjxpQ!A?V+~1h;6T<0 z&5VKF)kMJDmPdqdcR@t)h^bFGZ?|BJt>H$}|Bcy46+2Z6Jd-3zRy9CT zJCYmEJc7cLvv)n)=AfG+hy77&q5T6NWFBU!lz5;^1Jxx!v7$^ZrmklhB-}xjTQeMytKHklztb&RK8XB)4p7 zMLmEbbPp^t{CO~1B==whpd*C55c^2GfKl#PaXC%~wW=!a$f%bYT@EyqH5NrN`7(nLlp=RHL%0$5zfE9OaXMWWJk9YQtkh1 z?_At+&ek`cBxG_3F&RaYi}!*{Lq9Pa17@0DaWAkq#1ZODWg7PD8*ez28*MWSVqrf$ezZ@|{OgP1VQRxmql0`|@L{ex=6{jR^bgg^3FB+CXZE zciQUE29U%!C&@=rR+ELPT+}n(&M`F%UaU$+wz;tM=q8ei!SHA}di@SIaO7-{+QZ~3 zc3#=7h=u%G^+}&qI!)@$1`RQ_$;#_@?BP|Th=rzJo4-9`$=R?Mfw|3#wXUR%lFo_? z9~>d5iz!)}r?8Uk%QLpVy;1YnNOiHL7L2TXGS4g3{mXofh=eRC+3LqFkGG!q7FVMt z+V3u@+jHCD7EBV<>A?)-!i_%<6}$=}wMRpR9=-H&&zb%9Bxn(E4DG{D>LjtR`zP;a zEC?AC2F!u{AECm-!OX7orV0dM%mNqQM}#x@{=sVkANYOQXN`!3BYu(Bj3_?nc3=BS zq**D+ik!sw?@8xjI!bDwVP@oMVXC*y*bEYiqpz*|co#c~r%{!0%}Y;r4*IsD^KRZk|DBL60G%JNY&@5HT&vqAHLk?Wj0{4#SovL8`W{lTW;!fAb%atP!J9n^i$DZ&69W0S7B)l zLr+>J(dQD!Su0dUyfUr(L%R`iMXy3l4FM+16E8gQTQcqIB6$jxl5!s=1$@4$t?f4l zOL;itGSF%A9k`Ifz68yd>Jul8)iZi`_J}!q`_*>u!-dsj?axH^QM>2Z+`V4CdDr!e z>N5r@1n0j5yZp-b+!W&+*40C{8>yBfEwDm*W<)Te|C%rX!$4Mw)8gOgx)z_LCDQvj z_D(zCz3Hoz7N$J^FBd@R!NdU~CGTIL37|U}6Aj8)a6?$6Jg4?$e~zZu`k`-rqZhU$ zlr9NLNlDJ`zK05{uN^#Rb0MbwqeVp}vm}aJHart@!CGUzgxq9N+${UXm%#~I6fcGa zkp;D74c(3Qz;828To(kdnBH^4VCgxlIOkzMQ$S%|h^xl2fV{<8Kb6!a(#5lxDQZ7R zO9J-`4Fujk!s|m|zX>$xS^g7x&M?dHoIX6nDXLe*v4!VVSbi`%_8Dy=gI=KK&P6b( zM^ZvQnPkaG7ICF>fGZQkK?-KpqfVQur@#H&w71IwX@hVIE{f8rMF=M*n;?Mnjwwd^ zF9LtTuwrnC8+Br-`YJ{eq?1%z27Ajj?C&z{MkM!YJ1WyV4!Xbcu>y>dUb5HFhcig0 zsG&9GlOvv_Q7<4pGrO_wSz?oEFebKn>GG5^18KHpO0&9vo1e zRBiWU_pQ`EuA?R1FpGg6D&h=O6sgSTY52|8h^^mm0TqEF6Rd;T0!v4S8K$6^qBtre z`u#Y~0!gIVKuC0Ssn`^19$VK(t)*jTtkv4SRHZZNpezrXKLD}Bl)Z<+kD_)UEcm`sy zNrQu07pSOA)y|^1!=Gu{#vo)+o3?JfR6h@-rGOl4YPik{L1Z(%cBjlwGrXHb#8WwJ zi+)+^COO|cm`}YXEp0Q+&f^B46Z#Q`%dzO>Krn*JNUUO$A4boqb^M`HP^P2AJOdF9 z4V%*>=_Rd*WO-U4%cKmg>2sJ1`C)ribas00eRu4ayeqtHWWSd#?ObV@|HxQcMcM#} zLzP;Z{d>{ZN2aa~o;||Bd&A}!=RQtp=t+r@@M*#2grOF(IzbCU*?Nu|afGkTrr$6l zPjappfQ%yaUB@ST&qZ%Hh_tTgye_@6@v8S*RUt+OV&*YQ<~*tW8QN&`u6NRnaqQPy zOjj5_gA-7QGGnGo=Pz_A0s;g;U@QFQyhBPayse`v(AFtiBEcfy=GjEc1=o5Vi-KfU z(X_D)+6So))gJ*Pv!4$U4~i7KAE417s?@SnpV3j#D?sseiTEh?2dOlSJ$$u)g9-7R zs17zu%RhshVNkBfm6=D+MoY& z*Na&UC}C6cEJUw{8F=BVe=|%o(O!t3`{xXMJ5xF@d0y|4lP%6`Gphl>g^D-Dz4|Vo zKjG$&`uFQ~{MlE6*76^~hj~=3h&kIZ7DikvgzBUdA`EGU$HmS=SIiEdi=`(R6m|Wo z#z_Nu#xE{>y4I-pSODrTmJ1k5hFutsU0jIqM&sL6AIe+|=uK8tWo8w7BqeCkJ3_q& z0~NkFq8mi~#mnzA2}TR_-FM%;z*m5oufZs-Xdk>KgmV=#HWqNjjH`lhp{V#0)(R7n zQH$hq9CJbZQK;%NaKjfjo zGTJNnc!>F@uwkC+`OhGhEIlWucUpO$(1eN!)2p6ea|PYMu;kpfj(<0_wfEU4%Ub`} z;~qxxu?MqpAK98+`>CUaQ)+en9qEPDq=Ygw&LIC8#x{n`+L)KHSc`2zyh9v_E}O5Q zt>Q6XLrH)O*?C=CWM;!tUvzpV_WGe%ItYgy97?_flcuh}|6@KB!sLpIWD7_!a8rcM zF317^Y4@OaJ%eO`0hyzbO+OU`ed;Mh5>ab>d*i-braI|m;LgbO{4<+Xi136Tw=Q>! zAySnGZ5{p_QhU&IQ$l|w-*}Hb+r93|&S_^mYZ(~|i))i4n;*dv9dudgJzE@?xKNQ2 zlfv$qY|ohAYk+ibvIw48c0$pLF)rgjTHXRBDcKDAG`wYz-yS#t(?O2dL&}ANXPG41 z7Xxf;-E2-3wL4wEoDB5aB=h;wGTdM$KiREhb>mRC6Cy`Rdc39c9eV5nt8 zZpa1cwHG|0=Ku`A4B>fF<51~9Ak$X1h-QQWHY~u4$M?h|ZRu<@-pC$~m>cWeT z>YNJr#DcS-hf-q%1?-Qc`NDN`$jx^>P9nu%BGo^{lo1P3RYB28U6uW-Z3;*&^mP}> z3GZy860LKnHPkdHfRWgLv@A?y$T7k-O)yrdiN<97 zO|4536yKF3PZVU(=&JqdMGKn&JoTb38xvS+L+joRglKZA2*{8-dIPYNK!|(bzMu&_ zK)BaJSkD(0T2S~cMYAL>cg_sMtETpKajeWw{5*?vJcf+b`|IJT1_kmP^u-lN8+T5}RBTrJ>m=J0kEW_=pEz84 z{2;SmYz}$?uc4+Lp6i4T9sJ1Fsc!o;C{Y3*!)icYqOZeFz}XVGbeH-7-iBI)Du%KP zp}4U7HS8Jw&h;oP+?Cx#u4?Zbr$U_>VrBR@4aH)n*W z#_)}(TZH(WcuEohE9&>*eR7$Lz}BdN7<2+<4(`ZpboScwyR*9VJMH74hn)rn9>%B$ zQb(UY-g^q05X%dc3PQn~^0eA#-p}umpaEx(Fdv)g_0M0&ph;fN$qDhF)eqhZBgyCm z=8vM3Oxomn<;o%tJ$+gg^*yj1Ln-0L(zEXmyJmlw8H(E+Z3X>SB1T5RpX$Df2l`}6 zSO59zORzb0*XX^S>eKM{g3Us$sd>!kv9cZDG;@1qDbRE(WPjUU*`sMRsE6N;>d=D* zp$WQ->S3I1%c(kWO-Hv!40R`2iRsYXS#i&^E8h+>6P9aq9gO+nPsv7K+D5hm?Gry? z7`t!9`0A*Ze|;L{gAaGs?*K#XiZ;zqyX%sL275gwA7n(vanTFT{xjncV{8%JXnZOR zP<&t!{>!How2=>-h->}%!`}&79P!b{wy89{>z0EHtQ<))h!aNl7zXYpS&u)uuvDKa zgB^x1B%=(^s2fj7eL6z)MhQl^(24xaX9t}1yx`~I6&DuCO#T%(VOWswi9pAdEMn)}Y;% zNl&99Hvv8}A|bd^VH`X*hYW)`NwE7<=$uO>)?A%gj}BhD~!dvnAeRjyCkf zeK#0i3J&KtW-l?O%yj|~vUGb|_MB4v>ucPtOV!mgindz{Pb z*P=?R{tb_rakX)|pn4E;Pw!g8(|_;Vb#n-aU|*^msas@LoYP5d&NDqJF*-%5te>3h z#X{v?azr>5c)Ygtqa`um5R1+B@UyoWz(iH5Hb9IN7f8ghIzv(syz3LQ+MW0qA4m z;}r=5l!j88HmDFSNYEj|ECrP)^|wfCcBvD!*Hrbs-=ZJA04A3hn2y-Z5Dgrgj-}f7 zO={26PrCw*ML8>elLj4y=7)UaPjkfAZ+zN1l|V;ctEs->2h#%>2Zq9}fsf3$j8G0E zwT6GyqTlg)Gkyw>`m4GxX2Q*97X9_!>Mb60Pr0mV?upf&XFFfl758ZYVgL~~D;dF_R1Rqy9*ww-fW1vF!zux`TEL}3XS*vcyRIKFe?eX(w2D@;%th_* zYv*k)h&xOephF~tlqO6g`QG$rR9E4xMw~b*L4Xxf*5j(G1u>p8<;TMj$nttZeI&+x@|U$@TC@^ECUh($nEU!ujs)I~V?jLMz~VT8&W zt_Q&Xuv;x(0o{*eSwUa0h*j(^`(ZUwtLKoB+=Q;)F1;hpn^kzCW}I4zHxrx`abK*=_;UC9Kt@UMhNYY84!ukAa@o%uWjB1u}REy z1#dvD>+iXob7GQ^JmS3y*)&X@TT6uIFfO3&dn!2BQ3pG%ZJl#H0|HL9Dh3V(KcLnZ z%m?-2n0ki#GkGVkHFpZ|>#yusGdXbHks`)`ghr>38)7O|@e1HJM2N@L}J00Ew}(Y)Gw8-%ZK+wk38|#{E8+ zuyYBCA}7^0(}@4nc{F~Leki9lv~BK|8|I1D-;8@9nf(BY4`EMZpGI8wmTma$?T6hc ziZQr+)F~bQH{Bc~(+~Y(ndCM(^i*gYusqc$BD#yR^7XOsH1?b94Md#V|9(JS))Mua5P zAsoO?L)|+oqN=ZX_MxBx8DW`s+g#A7vy2{*X<2ccP=cehn6olS6|V1eT<>1Tx<@X( zN`E0ZA4xVsIZCjQ5gj))G1^45u{cG}_1%z0c}*x}WrYqEGx+)dwE<&C!=9kFXpnfr zW1Sc^F1FD5`XL6lJVu)fe;(;LvEZ}%k`ZHRJoQ!Dub=QAdZ;zl<^odez`w7P4j=Vl za|8F?cBp7kl?w5Fb8A}Aji>p$)Ol7sF^HXzM;||Z8X6oNoT1$JZO@@&Qdl`sjUbuw zMfDfHb^XdHrFrC)z9Aa@gYC@=e^c_B@QqX3~qh#5=Ky#eb@I1IQ1l7Pg z@4%h}Tl6`BGF4WiEHsbpYRah-_hCVQGg|s?L*_~W)pL5fckeFjV8JCUC*8CoY;{A+ zCg=#?LAZcbr~_6wQ=jAAoMeno=yAaQvM)EY_eqQJcnl4%>AGx0uFxb3j5 zz~{kNdrwlG5P&kYH&vHr*EmfX7-A|KN|GH-0IL8uCgYiF?-?}lLf?WMO}qjSk+_8g zvq)O2tOLuHd0;{~zjJ!A{D~j>7TjDIRH+E*pdyR@PwTj)NHw5(Z*O)^#1esPlMc~G zMf@PNBv1_L4Bh6p0i)v=dUW1hg#NNqH~!S)pK$7Bmju^^BdgH4N75NmIM-KXJ-yM_yyHMX4*atmK}d3kR!q@z_ZuM zC;$g&QKU_PFsAYw*ACg+q zNGUv+rhG$MnXwbIDbL(qj^xd_FohFqrW(owTaWvNWH=m--X0G;SIrU`-6QyvGZEvZE^<|-dQ9;*aE3Fpr5sea=ME!z#c(=?;!c2k=n@J2*RL2>+ zX1Q9fAOa#ZfEvcP0*vIbf`ie9kPsmnmLD{*G&nm z72?2y6}v2H=|vX}3c_}#;F0`&dG*ZpsvW||7#g34xKRK~D@jv-A$r{!;WN#YL>wb# zgJL1JjL!`bh7XO|&-Zh^f{0s`5>N&*r~D)OXXYvpo$~P{=ZjfB>eRlaAOAi2d@3$T z_`%(KCL3Ke60%0h3x?PNOT_Cq59hf!Kj_KgbVHiT>&(pytBVT~@mth~!aD;R38MlW z`WeE=wl}s4^Xs@ze~tOjmYsK4{6f9}X|TWBs4|~^8c!~sI(RR=2L*4Ov~2&oj-TgI zSm&FTv&BS*Pk$h8qhbZE%HmT7fhW*(iOV2&qf;U~A%Ac8yty?8n?)>Rxo1jeW!uy}AWC98fk@zTvVn>_F_|%T%A72Cz#Yvf(Uq7cQ_YeE?PDm};q$ z!?inTA}F_661uPS$4+P68Mz5x4(9P-B>>@IEW#pQ@P;AQ(4T;R??sV-7b2|FCC}Cp z()k(DydZW%W{(=Au~N3_nH=4$Q83R6*uikbgea%mL*M!saML$$hh;g_2`k6 zCHr|psumVfxcD^h+Qz$KG#nBJo<5m%_l?+HE=x>SY%t{=1F zaHRTk)Q>u#I3^4XGF!ILsgIm2t!e)f<1hXza|b?A;FN=Q2`>aT4hWxE5MheAVPQT* zHY`d6&koReRpj~xGbp5gx_B-unNKfFLn!Vup;F^tWo`)vN}dEz7osZRNsR<<-Gtdi zilrVqXp50)QiO_%7*m6gZN3Uv019_rbmav1_*b+T)X7j!tcH6!UsDMJ+o5t10|ZXbnk7H(aWT1M ze&;j?ErCr7v-N!P+H|H6>M z27JX{K#9T{$(+4{0XqOXurVVZyL6Ao`iChf(8Ms~WwtHMIxZ~v!DTr89{nx>N!+go zhXDL7eMLjGaqQX60tZpJf&{3bDU{+X$0{poA58c#Y=8BTjfZaZQGgzxeoD4we!KRn zFx^FdK?kRPbRGj>0Bv~G0iLVeg6yDlkQYp#=$EHw^ys6k4un%F_9djbmq{TKaqG)3Ui7om)S6@ z-kOaop*q7H#?&;E-G8|+Bhz=Bh3St=6SxDH1H{&3j!sJ9-H0uk5HSX#-cGQW!_fZPq3}`}XqMh) zR6yfN{#+Tf{_MevL@}Zz-xIptoav|bQ*INT5fD+F((`~0oUWe<9glq~5)QQn--HM{ zjaqDbz|8IhT^bBw-^A`|W|sL9Nyp-pGApMdOY1@Wzq7jDFX5AmE`sgOF1;9@L$koK zgIx#=Cmr5PGiz^x77vk_i5pISZ8J`FblFBnE12jA$^qjch?1D0j7sC6)Hcc>h9i%>jL8;Qi$1OoOkvzX%WSSP$n$D`R`!z zS2m6R=+L;kJf~^o{>oYkvCZK2bSLq*;Vq$X04Uw!Q&+ zv(m^Yt=FeytGK_J4GR+BRc6&VK5z}@GUqS%0vgRnR&m=)K zY%M%y8b4;t7-5*j1A-5@R=A`{8_(;0TT#~Qr!tp&IXmj!Js1-wE{H%7`&HZXdoKga znW`{pXCq|x^XoiPKc9nyy9}reJZ@Q(3oi0~_h3qrIc#0w)<=#(-z^$J)x%M5>l zEpMxzJu7aCWNA5xRJ!Ssthk;FOF`_!1$z`z9$4&%F%L-frpR95eFjM7WY?8{gEX7( z^wDuZ#4iNx+}IX7$waP4{7Kz41z_b2-B8QWh+&8Xm?n7ImPr7Tc*!G86}L;;rS3DyTn z?6ep{*CMCO`&eOBa_#d5UH32RmrslW`yYMDK37;Ya7X9&2=9;FzJ1%+AcqiA^Cm6t zNcq|zBXxBOv>K|}oz1^Z)Eh3>V8N-j9CZ8%#Oc~ zLG)I`|HI?+rw&o@)N&y@)&o+c{4e4{8wg9(N6j^1I`{c?H5gb&=*XuGcpI^wkJcuA zrYd(!XUlU5ZL0Ed%h9ZuxS<;9n-b&ips9&+LEwVs+_yCZH+rQky9rDKp>9ovl!WEE z&zNf%r~d$bJj5~m>cAf>k3K0a9Z2Upqo8s4^qpFSA$K2hKfmxsYDlhs)}s)fiRFKK z(3ilpJ(Uif-c}?ZHUu@bC0Nc~&{+84FLZjoX>0kMS5RJFK9?Wh9yJ}ZAR$DH=YSQ7 ze$&<3BgWTl0{kLBF?4PQH8JPs2>jD6lOhVM%_LzmeztrLMyL1h-Zdqc73?T68iE~& zo~~|x9&fJE%+HBPN0`g$;7KO_!-uVH%6(g%ljBb5UAIG0@;lw#A%YTyg*Ge~+m-H0 z8PuT}cW5w50Tj<lpF*IO2?lG&O`grHd3e^g$ zj|3G069>p7)Fj?v|9*D8C!2Brt)`xVL2G7sSM|-O*nSD!7xOm*v6tBokDomo4mD2u zxaI0xuL+h(e@r;c824UHjiV^f2MoYKK$VwpzWik*MH+acGPN%;c`$QDzwKd+bGO?V+w3_+i~!dujKE2l8##Sr&THLmN-~XcB#E6T@ePS{Zk-j}5*2HeIegpcZ(r zqod;kW-$q?C%&yxS($2^n36J%UrWfz>4-etB23Wu3E?Oy;125bb@H8ROSxGNI`X5F z(*`ElY%X^n_IX0(syCFgBC*pSILbQqyF6xCMTv+Yl?31fKVX2z^MLky-TTe#j}j{AA4Rvy$*)Os!m;Kqhhk=1#ded1NRtlh!b;jJFSe%d7^G^Zg*J$@gJU=N_(}RJ5spk zPu_MbT$0^ds^~1^x0>A7Z`OwqWd54?1FE51;MpJ@z4}9nOO_qjCJz2P*_N~qB5*WG zqnExYH#K~hmO&2XEVJCujq`SDs5E#Q-#F<+U&!4qfjeRqE>8(@laZO{{H(W(OvKl} zYG!nk84+D5FY}X=`mq1~qs+YjJ&^xxhX3yep{w8YA|*RnN3CiZ{%2vf$n=tl!|wkB DbdNxo literal 0 HcmV?d00001 diff --git a/doc/gallery/autodiff/vector_jacobian_product.ipynb b/doc/gallery/autodiff/vector_jacobian_product.ipynb new file mode 100644 index 0000000000..2ca6d2c7f2 --- /dev/null +++ b/doc/gallery/autodiff/vector_jacobian_product.ipynb @@ -0,0 +1,1071 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4dd1aa8a-f595-42c1-ad75-716f497c726b", + "metadata": {}, + "source": [ + "# Vector Jacobian Product" + ] + }, + { + "cell_type": "markdown", + "id": "ff301ea1be9ef0a6", + "metadata": {}, + "source": [ + "At the core of autodiff is the vector-Jacobian product (VJP), or in PyTensor's archaic terminology, the L-Operator (because the vector is on the left). The Jacobian is the matrix (or tensor) of all first-order partial derivatives. Let us completely ignore what the vector means, and think how do we go about computing the product of a general vector with the Jacobian matrix?" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "35d314e1728148d5", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:19.996660Z", + "start_time": "2025-07-06T20:20:19.594782Z" + } + }, + "outputs": [], + "source": [ + "import pytensor.tensor as pt\n", + "from pytensor.gradient import Lop, jacobian as jacobian_raw\n", + "from pytensor.graph import rewrite_graph\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "e805b85ceb6e8667", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:20.067703Z", + "start_time": "2025-07-06T20:20:20.065308Z" + } + }, + "outputs": [], + "source": [ + "def jacobian(*args, vectorize=True, **kwargs):\n", + " return jacobian_raw(*args, vectorize=vectorize, **kwargs)\n", + " \n", + "\n", + "def simplify_print(graph, **print_options):\n", + " rewrite_graph(graph, include=(\"fast_run\",), exclude=(\"inplace\", \"BlasOpt\")).dprint(**print_options)" + ] + }, + { + "cell_type": "markdown", + "id": "2d3afc0b692bd81", + "metadata": {}, + "source": [ + "## Elemtwise operations" + ] + }, + { + "cell_type": "markdown", + "id": "3b08aefde765313d", + "metadata": {}, + "source": [ + "The naive way is to create the full Jacobian matrix and then right-multiply it by the vector. Let's look at a concrete example for the Elemtwise operation log(x)." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "153da4ec5f8b5a71", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:20.128996Z", + "start_time": "2025-07-06T20:20:20.124242Z" + } + }, + "outputs": [], + "source": [ + "x = pt.vector(\"x\")\n", + "log_x = pt.log(x)" + ] + }, + { + "cell_type": "markdown", + "id": "91012e92f0eafee9", + "metadata": {}, + "source": [ + "If x has length 3, the Jacobian of y with respect to x is a 3x3 matrix, since there are 3 outputs and 3 inputs.\n", + "\n", + "Each entry contains the partial derivative of a one of the outputs (rows) with respect to one of the inputs (columns).\n", + "\n", + "$$\n", + "J = \\begin{pmatrix}\n", + "\\frac{\\partial y_1}{\\partial x_1} & \\frac{\\partial y_1}{\\partial x_2} & \\frac{\\partial y_1}{\\partial x_3} \\\\\n", + "\\frac{\\partial y_2}{\\partial x_1} & \\frac{\\partial y_2}{\\partial x_2} & \\frac{\\partial y_2}{\\partial x_3} \\\\\n", + "\\frac{\\partial y_3}{\\partial x_1} & \\frac{\\partial y_3}{\\partial x_2} & \\frac{\\partial y_3}{\\partial x_3}\n", + "\\end{pmatrix}\n", + "$$\n", + "\n", + "For the elementwise operation log(x), the Jacobian is a diagonal matrix, as each input affects only the corresponding output. For the log operation the partial derivatives are given by $\\frac{1}{x_i}$, so the Jacobian is:\n", + "\n", + "$$\n", + "J = \\begin{pmatrix}\n", + "\\frac{1}{x_1} & 0 & 0 \\\\\n", + "0 & \\frac{1}{x_2} & 0 \\\\\n", + "0 & 0 & \\frac{1}{x_3}\n", + "\\end{pmatrix}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "7c359f386b4cb918", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:20.195333Z", + "start_time": "2025-07-06T20:20:20.172499Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True_div [id A]\n", + " ├─ Eye{dtype='float64'} [id B]\n", + " │ ├─ Shape_i{0} [id C]\n", + " │ │ └─ x [id D]\n", + " │ ├─ Shape_i{0} [id C]\n", + " │ │ └─ ···\n", + " │ └─ 0 [id E]\n", + " └─ ExpandDims{axis=0} [id F]\n", + " └─ x [id D]\n" + ] + } + ], + "source": [ + "J_log = jacobian(log_x, x)\n", + "simplify_print(J_log)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "a3441a71c6ddb3ad", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:20.517814Z", + "start_time": "2025-07-06T20:20:20.237690Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 0. , 0. ],\n", + " [0. , 0.5 , 0. ],\n", + " [0. , 0. , 0.33333333]])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "J_log.eval({\"x\": [1.0, 2.0, 3.0]})" + ] + }, + { + "cell_type": "markdown", + "id": "f57772bc5c85c82c", + "metadata": {}, + "source": [ + "To get the vector-Jacobian product, we will left-multiply the Jacobian by a vector v. In this case, it simplifies to an elementwise division of the vector v by the input vector x:\n", + "\n", + "$$\n", + "v^T \\cdot J = \\begin{pmatrix}\n", + "v_1 \\\\ v_2 \\\\ v_3\n", + "\\end{pmatrix}^T \\cdot \\begin{pmatrix}\n", + "\\frac{1}{x_1} & 0 & 0 \\\\\n", + "0 & \\frac{1}{x_2} & 0 \\\\\n", + "0 & 0 & \\frac{1}{x_3}\n", + "\\end{pmatrix} = \\begin{pmatrix}\n", + "\\frac{v_1}{x_1} \\\\ \\frac{v_2}{x_2} \\\\ \\frac{v_3}{x_3}\n", + "\\end{pmatrix}\n", + "$$\n", + "\n", + "It is unnecessary to compute the whole Jacobian matrix, and then perform a vector-matrix multiplication. The `Lop` returns the smart computations directly:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "776b81a3c3039098", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:20.968461Z", + "start_time": "2025-07-06T20:20:20.953510Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True_div [id A]\n", + " ├─ v [id B]\n", + " └─ x [id C]\n" + ] + } + ], + "source": [ + "v = pt.vector(\"v\")\n", + "vjp_log = Lop(log_x, wrt=x, eval_points=v)\n", + "simplify_print(vjp_log)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b42553ba90bb69e9", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:21.211452Z", + "start_time": "2025-07-06T20:20:21.043061Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([4. , 2.5, 2. ])" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vjp_log.eval({\"x\": [1.0, 2.0, 3.0], \"v\": [4.0, 5.0, 6.0]})" + ] + }, + { + "cell_type": "markdown", + "id": "cf099c75e01f75f8", + "metadata": {}, + "source": [ + "## Cumsum operation" + ] + }, + { + "cell_type": "markdown", + "id": "a662a8e6b30e079c", + "metadata": {}, + "source": [ + "A pattern that will become obvious in this notebook is that we can often exploit some property of the Jacobian matrix (and that we want to multiply it by a vector) to compute the VJP cheaply. Let's take a look at the cumulative sum operation." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "f1c068fed4abb1a7", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:21.240090Z", + "start_time": "2025-07-06T20:20:21.227054Z" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1., 3., 6.])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = pt.vector(\"x\")\n", + "cumsum_x = pt.cumsum(x)\n", + "cumsum_x.eval({\"x\": [1.0, 2.0, 3.0]})" + ] + }, + { + "cell_type": "markdown", + "id": "a9b41bc35d7c1716", + "metadata": {}, + "source": [ + "The jacobian of the cumulative sum operation is a lower triangular matrix of ones, since the first input affects all outputs additively, the second input affects all outputs but the first, and so on, until the last input which only affects the last output. If x has length 3:\n", + "\n", + "$$\n", + "J = \\begin{pmatrix}\n", + "1 & 0 & 0 \\\\\n", + "1 & 1 & 0 \\\\\n", + "1 & 1 & 1 \\\\\n", + "\\end{pmatrix}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "87a87483-e92a-42aa-b42f-509739063888", + "metadata": {}, + "source": [ + "PyTensor autodiff builds this jacobian in a funny way. Starting from a diagonal matrix, it flips the columns, performs a cumsum across the them and then flips them back. A more direct way would do cumsum along the row of the diagonal matrix, but since a flip is just a view (no copy needed) it doesn't actually cost us much." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "1a5e9c70f2681436", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:20:21.366061Z", + "start_time": "2025-07-06T20:20:21.290442Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Subtensor{:, ::step} [id A]\n", + " ├─ CumOp{1, add} [id B]\n", + " │ └─ Subtensor{:, ::step} [id C]\n", + " │ ├─ Eye{dtype='float64'} [id D]\n", + " │ │ ├─ Shape_i{0} [id E]\n", + " │ │ │ └─ x [id F]\n", + " │ │ ├─ Shape_i{0} [id E]\n", + " │ │ │ └─ ···\n", + " │ │ └─ 0 [id G]\n", + " │ └─ -1 [id H]\n", + " └─ -1 [id H]\n" + ] + } + ], + "source": [ + "J_cumsum = jacobian(cumsum_x, x)\n", + "simplify_print(J_cumsum)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "2c932d61-40af-4c65-a45b-d397d5b0e3a0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1, 0, 0],\n", + " [1, 1, 0],\n", + " [1, 1, 1]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "J_cumsum.eval({\"x\": [1.0, 2.0, 3.0]}).astype(int)" + ] + }, + { + "cell_type": "markdown", + "id": "789fe3a70f42f6c5", + "metadata": {}, + "source": [ + "The left-multiplication of the Jacobian by a vector v has a special structure as well. Let's write it down:\n", + "\n", + "$$\n", + "v^T \\cdot J = \\begin{pmatrix}\n", + "v_1 \\\\ v_2 \\\\ v_3\n", + "\\end{pmatrix}^T \\cdot \\begin{pmatrix}\n", + "1 & 0 & 0 \\\\\n", + "1 & 1 & 0 \\\\\n", + "1 & 1 & 1 \\\\\n", + "\\end{pmatrix} = \\begin{pmatrix}\n", + "v_1 + v2 + v 3 \\\\ v_2 + v_3 \\\\ v_3\n", + "\\end{pmatrix}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "530901210f4afebe", + "metadata": {}, + "source": [ + "The final result is a cumulative sum of the vector v, but in reverse order." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "e4611d5f792ecd5e", + "metadata": { + "ExecuteTime": { + "end_time": "2025-07-06T20:24:45.832808Z", + "start_time": "2025-07-06T20:24:45.809493Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Subtensor{::step} [id A]\n", + " ├─ CumOp{None, add} [id B]\n", + " │ └─ Subtensor{::step} [id C]\n", + " │ ├─ v [id D]\n", + " │ └─ -1 [id E]\n", + " └─ -1 [id E]\n" + ] + } + ], + "source": [ + "v = pt.vector(\"v\")\n", + "vjp_cumsum = Lop(cumsum_x, x, v)\n", + "simplify_print(vjp_cumsum)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "e7695008-1622-45cb-ae36-5b9187d2abcf", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([3., 2., 1.])" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vjp_cumsum.eval({\"x\": [1.0, 2.0, 3.0], \"v\": [1, 1, 1]})" + ] + }, + { + "cell_type": "markdown", + "id": "d5a03eb5545625b5", + "metadata": {}, + "source": [ + "## Convolution operation" + ] + }, + { + "cell_type": "markdown", + "id": "7b27b04c-1abe-4f78-ad4a-35ff35614093", + "metadata": {}, + "source": [ + "Next, we shall look at an operation with two inputs - the discrete convolution." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "43f42d4e-46cf-4a19-97c3-5c9a6888d824", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0., 1., 1., -2.])" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = pt.vector(\"x\")\n", + "y = pt.vector(\"y\", shape=(2,))\n", + "convolution_xy = pt.signal.convolve1d(x, y, mode=\"full\")\n", + "convolution_xy.eval({\"x\": [0, 1, 2], \"y\": [1, -1]})" + ] + }, + { + "cell_type": "markdown", + "id": "26224b6d-0164-4421-989d-b277e1ec38cc", + "metadata": {}, + "source": [ + "If you're not familiar with convolution, we get those four numbers by padding `x` with zeros and then performing an inner product with the flipped `y`, one pair of values at a time" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "56075ceb-dc1d-4e17-8e95-e43dbc28783a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0, 1, 1, -2])" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x_padded = np.array([0, 0, 1, 2, 0])\n", + "res = np.array([\n", + " x_padded[0:2] @ [-1, 1],\n", + " x_padded[1:3] @ [-1, 1],\n", + " x_padded[2:4] @ [-1, 1],\n", + " x_padded[3:5] @ [-1, 1],\n", + "])\n", + "res" + ] + }, + { + "cell_type": "markdown", + "id": "bdb3a78a-ca4e-4adc-98f6-65663ffd0ec6", + "metadata": {}, + "source": [ + "Let's focus on the Jacobian wrt to y, as that's smaller. If you look at the expression above you'll see that it implies the following jacobian:\n", + "\n", + "$$\n", + "J = \\begin{pmatrix}\n", + "x_1 & 0 \\\\\n", + "x_2 & x_1 \\\\\n", + "x_3 & x_2 \\\\\n", + "0 & x_3 \\\\\n", + "\\end{pmatrix}\n", + "$$\n", + "\n", + "The constant zeros come from the padding. Curious how PyTensor builds this sort of jacobian?" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "aee5b3ed-1533-479a-abba-6375435c5268", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Blockwise{Convolve1d, (n),(k),()->(o)} [id A]\n", + " ├─ Eye{dtype='float64'} [id B]\n", + " │ ├─ Add [id C]\n", + " │ │ ├─ 1 [id D]\n", + " │ │ └─ Shape_i{0} [id E]\n", + " │ │ └─ x [id F]\n", + " │ ├─ Add [id C]\n", + " │ │ └─ ···\n", + " │ └─ 0 [id G]\n", + " ├─ ExpandDims{axis=0} [id H]\n", + " │ └─ Subtensor{::step} [id I]\n", + " │ ├─ x [id F]\n", + " │ └─ -1 [id J]\n", + " └─ [False] [id K]\n" + ] + } + ], + "source": [ + "J_convolution = jacobian(convolution_xy, y)\n", + "simplify_print(J_convolution)" + ] + }, + { + "cell_type": "markdown", + "id": "046b91b7-47cb-492e-9a7b-5c92de0cab93", + "metadata": {}, + "source": [ + "It performs a batched \"valid\" convolution between eye(4) and the flipped x vector. In a valid convolution, there is no padding, and we only multiply the sub-sequences that match in length." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "1ab05d9a-33fd-46a9-abfd-4de25cd6e520", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0., 0.],\n", + " [1., 0.],\n", + " [2., 1.],\n", + " [0., 2.]])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "J_convolution.eval({\"x\": [0, 1, 2]})" + ] + }, + { + "cell_type": "markdown", + "id": "05577d83-091a-48fb-a383-854f69c4cab5", + "metadata": {}, + "source": [ + "Following the theme, is there any special structure in this Jacobian that can be exploited to compute VJP efficiently?" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "89caf646-0213-4aa6-9aca-b7cb0230c8bc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Convolve1d [id A]\n", + " ├─ v [id B]\n", + " ├─ Subtensor{::step} [id C]\n", + " │ ├─ x [id D]\n", + " │ └─ -1 [id E]\n", + " └─ ScalarFromTensor [id F]\n", + " └─ False [id G]\n" + ] + } + ], + "source": [ + "v = pt.vector(\"v\", shape=(4,))\n", + "vjp_convolution = Lop(convolution_xy, y, v)\n", + "simplify_print(vjp_convolution)" + ] + }, + { + "cell_type": "markdown", + "id": "6b369467-9ceb-4f4f-86f7-486801ff3275", + "metadata": {}, + "source": [ + "It's just the \"valid\" convolution between v and x flipped. Our Jacobian has a [toeplitz structure](https://en.wikipedia.org/wiki/Toeplitz_matrix), and the dot product between such a matrix and a vector is equivalent to a discrete convolution!" + ] + }, + { + "cell_type": "markdown", + "id": "4d568e5b-310d-4a4e-a89e-198a337a379b", + "metadata": {}, + "source": [ + "$$\n", + "v^T \\cdot J = \\begin{pmatrix}\n", + "v_1 \\\\ v_2 \\\\ v_3 \\\\ v4\n", + "\\end{pmatrix}^T \\cdot \\begin{pmatrix}\n", + "x_1 & 0 \\\\\n", + "x_2 & x_1 \\\\\n", + "x_3 & x_2 \\\\\n", + "0 & x_3 \\\\\n", + "\\end{pmatrix}\n", + "= v \\ast x_{[::-1]}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "b4811469-cb72-4303-98a9-cc75140b323d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 8., 11.])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vjp_convolution.eval({\"v\": [1, 2, 3, 4], \"x\": [0, 1, 2]})" + ] + }, + { + "cell_type": "markdown", + "id": "d5623e88-9210-4d2c-bcbd-6d37b991a01d", + "metadata": {}, + "source": [ + "## Transpose operation\n", + "\n", + "For a final example let's look at matrix tranposition. This is a simple operation, but is no longer a vector function." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "cd3d779b-6bc3-4462-9904-82f6177ad839", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(None, 2)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = pt.matrix(\"A\", shape=(2, None))\n", + "transpose_A = A.T\n", + "transpose_A.type.shape" + ] + }, + { + "cell_type": "markdown", + "id": "60a08160-10b8-494d-8031-621a8e7d0521", + "metadata": {}, + "source": [ + "To be able to think about the Jacobian (and then the VJP) we need to look at this operation in terms of raveled input and outputs." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "ab4fc0f0-882e-46ff-8129-08731e26a025", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0., 3., 1., 4., 2., 5.])" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "transpose_A.ravel().eval({\"A\": np.arange(6).reshape(2, 3)})" + ] + }, + { + "cell_type": "markdown", + "id": "e3043a1b-179e-49c1-9e3e-587740b8ac7d", + "metadata": {}, + "source": [ + "The Jacobian is then a (6 x 6) permutation matrix like this:\n", + "\n", + "$$\n", + "J = \\begin{pmatrix}\n", + "1 & 0 & 0 & 0 & 0 & 0\\\\\n", + "0 & 0 & 0 & 1 & 0 & 0\\\\\n", + "0 & 1 & 0 & 0 & 0 & 0\\\\\n", + "0 & 0 & 0 & 0 & 1 & 0\\\\\n", + "0 & 0 & 1 & 0 & 0 & 0\\\\\n", + "0 & 0 & 0 & 0 & 0 & 1\\\\\n", + "\\end{pmatrix}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "4765b01a-8078-48ac-ad22-1ce1dabf9ea7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 1., 0., 0.],\n", + " [0., 1., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 1., 0.],\n", + " [0., 0., 1., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 1.]])" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "J_transpose = jacobian(transpose_A.ravel(), A).reshape((6, 6))\n", + "J_transpose.eval({\"A\": np.zeros((2, 3))})" + ] + }, + { + "cell_type": "markdown", + "id": "db3f0b31-9514-4724-ac77-755bbb0c2be7", + "metadata": {}, + "source": [ + "PyTensor builds this Jacobian with two reshapes and a tranpose of an `eye`." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "c65f6740-f809-4fc0-831f-61bd2e09f510", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reshape{2} [id A]\n", + " ├─ Transpose{axes=[0, 2, 1]} [id B]\n", + " │ └─ Reshape{3} [id C]\n", + " │ ├─ Eye{dtype='float64'} [id D]\n", + " │ │ ├─ Mul [id E]\n", + " │ │ │ ├─ 2 [id F]\n", + " │ │ │ └─ Shape_i{1} [id G]\n", + " │ │ │ └─ A [id H]\n", + " │ │ ├─ Mul [id E]\n", + " │ │ │ └─ ···\n", + " │ │ └─ 0 [id I]\n", + " │ └─ MakeVector{dtype='int64'} [id J]\n", + " │ ├─ Mul [id E]\n", + " │ │ └─ ···\n", + " │ ├─ Shape_i{1} [id G]\n", + " │ │ └─ ···\n", + " │ └─ 2 [id F]\n", + " └─ [6 6] [id K]\n" + ] + } + ], + "source": [ + "simplify_print(J_transpose)" + ] + }, + { + "cell_type": "markdown", + "id": "0e39d268-12c3-463f-9ec7-e8b48bbd95d4", + "metadata": {}, + "source": [ + "To recreate the outcome of `Lop`, we ravel the `V` matrix, multiply it with the Jacobian defined on the raveled function, and reshape the result to the original input shape." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "2d09cb7b-60f1-4d22-adbe-fb88d487ba86", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0., 2., 4.],\n", + " [1., 3., 5.]])" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "V = pt.matrix(\"V\", shape=(3, 2))\n", + "naive_vjp_transpose = (V.ravel() @ J_transpose).reshape(V.shape[::-1])\n", + "\n", + "vjp_eval_dict = {\"V\": np.arange(6).reshape((3, 2)), \"A\": np.zeros((2, 3))}\n", + "naive_vjp_transpose.eval(vjp_eval_dict)" + ] + }, + { + "cell_type": "markdown", + "id": "9164ca89-a49c-4768-a5c4-f16e16ffc0b4", + "metadata": {}, + "source": [ + "Because J is a permutation matrix, the multiplication with it simply rearranges the entries of `V`. \n", + "\n", + "What's more, after the reshape, we end up with a simple transposition of the original `V` matrix!\n", + "\n", + "Unsurprisingly, `Lop` takes the direct shortcut:" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "e29d3283-77f6-4539-8d41-851848267cd6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Transpose{axes=[1, 0]} [id A]\n", + " └─ V [id B]\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Lop(transpose_A, A, V).dprint()" + ] + }, + { + "cell_type": "markdown", + "id": "6de26bd3-4e43-444a-abf9-eee682d201b0", + "metadata": {}, + "source": [ + "## VJP and auto-diff" + ] + }, + { + "cell_type": "markdown", + "id": "4f92c6a4-27d2-4854-8b66-038eca62ebb7", + "metadata": {}, + "source": [ + "It is time to reveal the meaning of the mysterious vector (or reshaped tensor) `v`. In the context ouf auto-diff, it is the vector that accumulates the partial derivatives of intermediate computations. If you chain the VJP for each operation in your graph you obtain reverse-mode autodiff. \n", + "\n", + "Let's look at a simple example with the operations we discussed already:" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "1a8a50dd-d62b-4229-b885-e96f4e43f76e", + "metadata": {}, + "outputs": [], + "source": [ + "x = pt.vector(\"x\")\n", + "log_x = pt.log(x)\n", + "cumsum_log_x = pt.cumsum(log_x)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "8a32e30e-5885-49e1-8973-e734d66f05b2", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True_div [id A]\n", + " ├─ Subtensor{::step} [id B]\n", + " │ ├─ CumOp{None, add} [id C]\n", + " │ │ └─ Alloc [id D]\n", + " │ │ ├─ [1.] [id E]\n", + " │ │ └─ Shape_i{0} [id F]\n", + " │ │ └─ x [id G]\n", + " │ └─ -1 [id H]\n", + " └─ x [id G]\n" + ] + } + ], + "source": [ + "grad_out_wrt_x = pt.grad(cumsum_log_x.sum(), x)\n", + "simplify_print(grad_out_wrt_x)" + ] + }, + { + "cell_type": "markdown", + "id": "5c788690-c1b9-4462-b78f-b4c16bfdc0b0", + "metadata": {}, + "source": [ + "You may recognize the gradient components from the examples above. The gradient simplifies to `cumsum(ones_like(x))[::-1] / x`\n", + "\n", + "We can build the same graph manually, by chaining two `Lop` calls and setting the initial `grad_vec` to `ones` with the right shape." + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "b7b81d74-fc84-4de9-97df-0cb47b745825", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True_div [id A]\n", + " ├─ Subtensor{::step} [id B]\n", + " │ ├─ CumOp{None, add} [id C]\n", + " │ │ └─ Alloc [id D]\n", + " │ │ ├─ [1.] [id E]\n", + " │ │ └─ Shape_i{0} [id F]\n", + " │ │ └─ x [id G]\n", + " │ └─ -1 [id H]\n", + " └─ x [id G]\n" + ] + } + ], + "source": [ + "grad_vec = pt.ones_like(cumsum_log_x)\n", + "grad_out_wrt_x = Lop(log_x, x, Lop(cumsum_log_x, log_x, grad_vec))\n", + "simplify_print(grad_out_wrt_x)" + ] + }, + { + "cell_type": "markdown", + "id": "a2e9e121-7363-41f9-9c1d-d6792db36bcc", + "metadata": {}, + "source": [ + "Similarly, forward-mode autodiff makes use of the R-Operator (Rop) or Jacobian-vector product (JVP) to accumulate the partial derivations from inputs to outputs." + ] + }, + { + "cell_type": "markdown", + "id": "eff94009-64f2-4c5d-b4d5-5ff36a86b7fa", + "metadata": {}, + "source": [ + "## Conclusion" + ] + }, + { + "cell_type": "markdown", + "id": "123a10f2-ef4d-4cd9-86b7-b2b7975a1c31", + "metadata": {}, + "source": [ + "We hope this sheds some light on how PyTensor (and most auto-diff frameworks) implement vector Jacobian products efficiently, in a way that avoids both having to build the full jacobian and having to multiply with it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0614f4ce-6a30-4fdc-8d74-bea8203ff3c4", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/scripts/generate_gallery.py b/scripts/generate_gallery.py index 15e94ca7f4..3409cb1b3d 100644 --- a/scripts/generate_gallery.py +++ b/scripts/generate_gallery.py @@ -57,6 +57,7 @@ folder_title_map = { "introduction": "Introduction", "rewrites": "Graph Rewriting", + "autodiff": "Automatic Differentiation", "scan": "Looping in Pytensor", "optimize": "Optimization in Pytensor", }