From ad9e84e93eb68f107ae972d33458fc98e312097e Mon Sep 17 00:00:00 2001 From: Niclas Rieger Date: Mon, 2 Sep 2024 17:02:34 +0200 Subject: [PATCH] feat: add CCA as cross-set analysis --- .../_autosummary/xeofs.cross.CCA.rst | 46 +++ .../_autosummary/xeofs.cross.ComplexCCA.rst | 50 +++ .../_autosummary/xeofs.cross.HilbertCCA.rst | 50 +++ docs/api_reference/cross_set_analysis.rst | 3 + docs/auto_examples/auto_examples_jupyter.zip | Bin 104784 -> 105206 bytes docs/auto_examples/auto_examples_python.zip | Bin 62435 -> 62435 bytes tests/models/cross/test_cca.py | 292 ++++++++++++++ tests/models/cross/test_cpcca.py | 67 ++++ xeofs/cross/__init__.py | 8 +- xeofs/cross/cca.py | 356 ++++++++++++++++++ 10 files changed, 870 insertions(+), 2 deletions(-) create mode 100644 docs/api_reference/_autosummary/xeofs.cross.CCA.rst create mode 100644 docs/api_reference/_autosummary/xeofs.cross.ComplexCCA.rst create mode 100644 docs/api_reference/_autosummary/xeofs.cross.HilbertCCA.rst create mode 100644 tests/models/cross/test_cca.py create mode 100644 xeofs/cross/cca.py diff --git a/docs/api_reference/_autosummary/xeofs.cross.CCA.rst b/docs/api_reference/_autosummary/xeofs.cross.CCA.rst new file mode 100644 index 0000000..1149b43 --- /dev/null +++ b/docs/api_reference/_autosummary/xeofs.cross.CCA.rst @@ -0,0 +1,46 @@ +CCA +=== + +.. currentmodule:: xeofs.cross + +.. autoclass:: CCA + :members: + :inherited-members: + + + .. automethod:: __init__ + + + .. rubric:: Methods + + .. autosummary:: + + ~CCA.__init__ + ~CCA.components + ~CCA.compute + ~CCA.correlation_coefficients_X + ~CCA.correlation_coefficients_Y + ~CCA.cross_correlation_coefficients + ~CCA.deserialize + ~CCA.fit + ~CCA.fraction_variance_X_explained_by_X + ~CCA.fraction_variance_Y_explained_by_X + ~CCA.fraction_variance_Y_explained_by_Y + ~CCA.get_params + ~CCA.get_serialization_attrs + ~CCA.heterogeneous_patterns + ~CCA.homogeneous_patterns + ~CCA.inverse_transform + ~CCA.load + ~CCA.predict + ~CCA.save + ~CCA.scores + ~CCA.serialize + ~CCA.squared_covariance_fraction + ~CCA.transform + + + + + + \ No newline at end of file diff --git a/docs/api_reference/_autosummary/xeofs.cross.ComplexCCA.rst b/docs/api_reference/_autosummary/xeofs.cross.ComplexCCA.rst new file mode 100644 index 0000000..0c09990 --- /dev/null +++ b/docs/api_reference/_autosummary/xeofs.cross.ComplexCCA.rst @@ -0,0 +1,50 @@ +ComplexCCA +========== + +.. currentmodule:: xeofs.cross + +.. autoclass:: ComplexCCA + :members: + :inherited-members: + + + .. automethod:: __init__ + + + .. rubric:: Methods + + .. autosummary:: + + ~ComplexCCA.__init__ + ~ComplexCCA.components + ~ComplexCCA.components_amplitude + ~ComplexCCA.components_phase + ~ComplexCCA.compute + ~ComplexCCA.correlation_coefficients_X + ~ComplexCCA.correlation_coefficients_Y + ~ComplexCCA.cross_correlation_coefficients + ~ComplexCCA.deserialize + ~ComplexCCA.fit + ~ComplexCCA.fraction_variance_X_explained_by_X + ~ComplexCCA.fraction_variance_Y_explained_by_X + ~ComplexCCA.fraction_variance_Y_explained_by_Y + ~ComplexCCA.get_params + ~ComplexCCA.get_serialization_attrs + ~ComplexCCA.heterogeneous_patterns + ~ComplexCCA.homogeneous_patterns + ~ComplexCCA.inverse_transform + ~ComplexCCA.load + ~ComplexCCA.predict + ~ComplexCCA.save + ~ComplexCCA.scores + ~ComplexCCA.scores_amplitude + ~ComplexCCA.scores_phase + ~ComplexCCA.serialize + ~ComplexCCA.squared_covariance_fraction + ~ComplexCCA.transform + + + + + + \ No newline at end of file diff --git a/docs/api_reference/_autosummary/xeofs.cross.HilbertCCA.rst b/docs/api_reference/_autosummary/xeofs.cross.HilbertCCA.rst new file mode 100644 index 0000000..6f696a7 --- /dev/null +++ b/docs/api_reference/_autosummary/xeofs.cross.HilbertCCA.rst @@ -0,0 +1,50 @@ +HilbertCCA +========== + +.. currentmodule:: xeofs.cross + +.. autoclass:: HilbertCCA + :members: + :inherited-members: + + + .. automethod:: __init__ + + + .. rubric:: Methods + + .. autosummary:: + + ~HilbertCCA.__init__ + ~HilbertCCA.components + ~HilbertCCA.components_amplitude + ~HilbertCCA.components_phase + ~HilbertCCA.compute + ~HilbertCCA.correlation_coefficients_X + ~HilbertCCA.correlation_coefficients_Y + ~HilbertCCA.cross_correlation_coefficients + ~HilbertCCA.deserialize + ~HilbertCCA.fit + ~HilbertCCA.fraction_variance_X_explained_by_X + ~HilbertCCA.fraction_variance_Y_explained_by_X + ~HilbertCCA.fraction_variance_Y_explained_by_Y + ~HilbertCCA.get_params + ~HilbertCCA.get_serialization_attrs + ~HilbertCCA.heterogeneous_patterns + ~HilbertCCA.homogeneous_patterns + ~HilbertCCA.inverse_transform + ~HilbertCCA.load + ~HilbertCCA.predict + ~HilbertCCA.save + ~HilbertCCA.scores + ~HilbertCCA.scores_amplitude + ~HilbertCCA.scores_phase + ~HilbertCCA.serialize + ~HilbertCCA.squared_covariance_fraction + ~HilbertCCA.transform + + + + + + \ No newline at end of file diff --git a/docs/api_reference/cross_set_analysis.rst b/docs/api_reference/cross_set_analysis.rst index 970b6a8..f84fb0c 100644 --- a/docs/api_reference/cross_set_analysis.rst +++ b/docs/api_reference/cross_set_analysis.rst @@ -8,10 +8,13 @@ Methods that investigate relationships or patterns between variables **across tw :template: custom-class-template.rst :recursive: + ~xeofs.cross.CCA ~xeofs.cross.MCA ~xeofs.cross.CPCCA + ~xeofs.cross.ComplexCCA ~xeofs.cross.ComplexMCA ~xeofs.cross.ComplexCPCCA + ~xeofs.cross.HilbertCCA ~xeofs.cross.HilbertMCA ~xeofs.cross.HilbertCPCCA diff --git a/docs/auto_examples/auto_examples_jupyter.zip b/docs/auto_examples/auto_examples_jupyter.zip index 6ee59c6a06c061714fe8f2741a7924dccb185221..186ccd5aee0f3a2528da5c57cb29ec78c2089b0a 100644 GIT binary patch literal 105206 zcmeIb?~i0xdfzu2*#b=1unbGE0893D(<|0c(>48LhQpz@7aq=%)U3GVgflDdYKG~) zU3I&=cB<>v)~%ZEUV=h4Z}L`wY(sz{Brg&m5S$ni_Jt7ml?+2b5W_G0)^TkA1vU`K z0b)S;e4lg9J?GYs?&(?f^bAwnZL+HFx##zJ&hvbqU-y6GXa2&^U(r8*^S4*N|N0+4 z_)mZQb62kD?{D$jdes~3_lxVpemUOml!wE9alBoWd+pxvWU%}BYuQSt==ZCYcd|#X zWf|Wy_sjP9WLWU)%3(fw(k+h$E7$DVVKL6T`8em%=U=#=s&X>wsH|$>^~b(!w(@#* z*R`0v_woB#KFIqgRj-QAt(eGzRdH)GdM&)GxAUi0J zvT>OWi_u;=I?VF9wP)oZRy#f@vLBVSOXJ<_uo(5aJ^Omp#`}=c949dy@4iJQ_2b^ry%MMW?8$d~}ix^Ujlezo@KG3E&jGsa=b`S%mR1u%WPsgKnRgAaGp@vemy1i=H&ri0C z;~~qiosY+(YU};HUlmu+oSgG-zxT_L+e^i9uI(PS(xfVE<=*?`&g#9vn7wdVZ?{Zq zMVGUkonn;t2IK9?(;w~ZX!UGfZBt+7I?|nNi#2J71JW*s#bDcayY+Y_cKh+l)!1#` zyMQTvsG}ODtWyrgOc`f49JXu8NxNpp;%ULj>2M6E&Rv!#TaaU36j`^ow^xiflsr2u ztKk6))yewB)1qIsvyaNLRxux2Yw7~$y???{rU|-mRAz_t=b*Y4=b$PfCfPy$w8(xi zsd$q@-Ggo~e_ z>gVZ+WY3oL(Np-Wj3rk z4`N4)ymOH4bDL)39yCfF{_2O5d{l5wxL?XW`>9jHe-z(^M?#u%<;BPO#v+-6e8#n!EaAX|-X z_FL<>u0iV#w;r$9LPNDLV9Vcw{uFsv^KYb&6JV?BnV2TENcbTgldO{u*f=cyJ{KeF z8<$@M(-zVW&+ve&(=Ra}RH^5xvMCm-UzSg_ZyoubSD)WDd%Y2aT66XCP^@am?#=t# z8vMA&F0LS9ZaL$S?Ws*oD_QjLs?F5g*X&99N}QbVa%>mos^3z}TkBWb;^tbNvd>1c zOR{nO)6;eZ+uXX4%{9)npjFoIJ#nm6T()f3>+obUgt9Pia9e|HR8+;`Za)yNe%Wz1 zqg(T2z?0qasH6h;g&(l>xI&8@b{`=Cru|R;aC7ZX{)S_hB!)AKEt1R`@Y)@a&s@EOAbz>N?{H zT+2F#`H&+qKI-+kxyEOniq|qpP%Ei;%f++8rWn9@a%W56)m5|-q?ti}Xae&6lkovJ z_NEu%`uSjgf;>q1m7zWI4?=nO;1C)!5o=U%Ym|fqDGfKl<_CyZsASuITTVr!kKL4s9C>GU&b_$m0{({-+{x0}eUN zPYzj8^A$Y%_#k^P6&vFIj-O1p3SwT?>-P~lkRgn%{TTs({Z~a^$-YSUjLvgc@XWIe3%!6-Icyu>{AlH524vG5h%T==BW7Lq&K2J z4OC>W7Cpaup+Qj|>;Hy4G*sT9f#*v9_^p zmm?HOZK@6WQwf!eXP|S!ZJBk;!lZDM?GJKLNGl%pw=t=HFRPl+aHrpcew6)kzlXRT zE4}1BSHaPq0pU>O=}oSxHXmZY-HZIYy$=gC+&uyA5q{AQfDq-#ahQ0&!7jfof5YPDdGUEn(*=Xc8Qg9%MMJM}x-c4;R9#vqZha!=l^c4iK9R&{m^pk!yu~55=W(dD{~0_sd;@cKL|Q z7u?%!(YP_cUod_YwK@a%eY(Y67sMhjJ6PY zh*9k8dmlacIIGl8kvi}_;bn%?NRUbzzho-0*V_m1Ic8Uw>^9kz>`=A97&fjUU+h;v zR$Hz0wQJeUtI3vhJ8YIK&X2e9<40?sB}L;-TU31gvx^WHKPqe3bBx+VJV}?G=CClf zj_odpIG2J%B~V?vCmfbG=o+>6)PCwhYlPI*_iod3m+?~9%T%h)6a=7Q$ZHebf$n^k@=ov^h#F#BvpDb{aWUs%< z^iH$Bh=?tqc~h^X8)wafX^v)I3=tG%fu=ZJq==p_p~jl>OU=_#U`j}Hr&j6+)Z~a@ z+`2Cv9|c?F1a4!0baHw^41mayV+WPQ`&K+P>G~zNlq%^KMi`-0TVik>o9nB-fW!4k zVKg?#skK66`hL;u>*U4m`6YOHkT)jXh+n@khx1w6o5`%tc`tCfv8FS4-dIqCbG9R= za>Hl7oD$&+*aN2%SYslE;$V0IYru6ZjyP&OzA9w3-r4Elbee`!8OpeNi(YI{)g|Wo z>SC__Kl{i3=3Ves{r!!3nClBdTMYmmb48?7IKoqnp7GS6XE70!L)#+6-!gahU`HZvAw);v#}jfyg|}_YoTg`>|4!d;6esunkoQNZE

KwQ9}@(rfd?vOe|bn!x9T`nh5Gt zT^ttvVF1Fzakh>B@?75IjYCPkiLsxc`^4BU?~S$-;9iI^Kb^UY&M$A?AN6ES33&nQ zm|zwp6p&-zfFQ1Mys%}%99D1+r63&)FRv;Vxfu_k&8BlR!BconD-!D z-C_@mE#J6dt+n103oQ&*cXq%5AkYFBRWEj&0pa1-i1EmN%ZL`reo@8V9y$WO5s4+p zU@^jeF54}}M@8Yl>RQ@T1HRMFh7dHRs6c4E>7a%rbB9=8x7U#skfWD_3EhGUGkS2R z)q40tb;ES}z*ZX3YNsqmz)RBEK^if~iZ#|Ef0&W#1T&vQ)}VWWcDUC`VVw6_z*>oD zVA4dUOm4=uZumJ%=411WrE13`Of0ESh7z z`?K>9cYA6C%PW_8$w1ODM9Og~II9=uCsNvd06|U);yb{ISiW6D(W6l@=dp``)@I>g^G)e3v#~Y7SReETFeu1trr4N| z-U)H#%u%w$8sA&yOKus&p+QPU`631)+|hAZ@<|vJ)Y^vul2nqg5;9!KN6z$x&3#aI zWtDIOd5~QfeVG#Ig6TB^*h9MMPIF4d5x4wVkbSnlG;$a({%^Wg-4Kyenw4WPOlFP} zV#T;G{XDemer|^gbldvb@`!x0e=pRc^$sdQCtpjz^i*3mu1q>qNL$Apguk66CdGnH38CC)4f6eWp*NNwxd&vy75npu0LmJ#Iqm%X4 z8VuS#?nqL=jK&8q&$`ybsv^Y7#lT>;-Uqgr`uoLc z*!6H^%C;Bc*!!Cp?jW$Oj<|pYzor=H0v0wNy@iFnwMI>bGv5q8Ma~pL2k|oRtJgCF zj~dpk;@}h2-%M0UdF|k-i8ExqN`&`IUSR2&o7mfoqb)_uN)z8j zKH{qwC(}gy@j5l|q}|i)i9oAy6EvHuaockV0g%XAbW(LRxZHt2&d+W_q?zPBXFyUo_Q;Ped)t zQEk1tky=3ccU(&}ne<>3xeKZc1>=vYU^j9cW>PLli5EaW?ki-0afIvfuwC~45rYbx z6gvc9@dJKINOU@3PcWn1mr2a$@ePUjk5}ZY^cc^U$19(7-^Cp|L1|opr^s|Fu!@J* z%{fo=V?CYw&=60=#zU1*2WMW-yIuTm`{igCJtT$PK8=0+B7TgYuta4f62|xhp5LlX zQ*X67D{rNAi$`salNqvZ=Uq9XwvG(b-b(kS%}uzeT3EZn*1F2K_zygOuRA_)a$&Xg zNaW7_KecT&!`Vk_P^0pVw)U$>oWY60D$5BqIKhHa=i0n zj0hR!oa$)$owVIJF|UKX#c%$QL{?pUI5jh)h2;EMixeJjANTWJE|1464Wsh-9MPp; zOh1n8zAhFlxm0bK*Roa)XaK}qUDc2alo5vSw|>+cw(?_;`FPZ?rhbIlg*o<*c?Ekh zC@T4zmLuHj+VWP2&8nQ$zi`0UvNzv!2Vv^9+zlH~&L>st5li(bonO}hP0Y!eCFW0E zQp4P-jnBi>KCl}g`~00R7^hB9K5q3Cs36_jW?XGuKMB>^%%`)W{Zmbdh2oa%JE9wn zT5VG#^I;eV?lZi|5*t^)8m`aw4~F~9d8DyRY3|jP=KlWg|G_VQ8v`2s%>-+*>$m8> zz+g?rHXzHTwk7d!$4`3Rdhc-9BQgt|-NyhV<-S>^ypOYnsB3gWlqEbQ0aqog*U>g$ zi!s2yUmjWl(|TBhfVY^c4#h}QD=6+USgSX!H+05mgQ~420>l__xQnyzSP_}b;T#n8 z*UKXU$KZ9{ouG>;;}BGJZy-Nw0?i2FaRboXPP?ZdKZZ`h?G>X@|(`BT7$S24K?8k#s`3C#Q<1coylRlP}k?3&csx~^-9&FpF8uh zyYbVzGV=)?Z^fQC0~|Msti>etu<^NCi~Npvu~2v!_1qmU7tuzyr&%GSNA;bo zwVvFl!DhGFM;M|FMN&$(W#|G{cbu%&rhhy&@+KVYjln)O3>t~gg^QEPPQAJ22a5~9 z(PaW!=S8dV)!i)#=E0_}1__55V5GV688@0Wm1qJ4xy1z=k5Drfo7O`rrA%v{wOKiGlAi>fn71Z-^~+!9CCMzw|UgGJ%#H$20Vo#S0Vkc99hz_I*#;L1&mSIS)KSo9 z7~ODGKKtlyBqLB4c=;ToQ|UPFRw~yumo8LYkM06ab%4kXPLa>F#f*!W-CKr#RfbsT zy4{4Yv~YqUXkoY&bu!|b!<%s7!i&J{VACWfDh7kI8tg(fxYJCleKRl`YC8a!vBH|U zERGL)yBOJYd&3fj$_$i7_IcG<)X0T@mk1@Ju)(dM+YPZ&CfN8%V@9wBs_jHgcp!+U zd97pQp|iv&stH!+viK0Mhklo?nhlZT`(wQ@LU4cjmx4qQ&k4@Mi2^pd$b-W8H!t_(bB4Dn{%Ez5b{?+ zA!na{$r$5|`+E~R1b|El1E^UwmJi{V2n;$JJg-8a!RJC6jqJ^<5G3^66~+@#BY)fo zIep8yYM=qr*q>&IhsU~H?7!^683<66g=sVpEVf|C+1fLrgZcWi7TkEof{}m|=3pW2 z$mkaUkIDP%i{EGF&VGMm@%t>|ocH4-8aTLJA7t5#w*pbzx=rt&=D@8Y1maF556S9? zpLU7=)?ioNZr5B$a%!P6pzmvKJmnv+1P)6z)F+<+5n6qFL*~(Uj%DBZ=9_7kl55uO zk@@v|9%iO`Hq|+*XY0+c-rpM2oe$-oq4ROVg|T9zdobR{*dBMzR8!N$Uu_vY4rkyi z?BTGtMfAWW3j37^ANTM6>finFdw==L75%+FV^fp($F06N^O}VgGI0mTPYVk#tz|CW zS|%<-{6yY^XE%63;FQ3Vczb>2El%Q0&6qlfO*3I+%(dk{aCf>7i4n#H2w<-Qta7@+ z$vN#!&TnoKc_o%fug4$K`_qZZzDSUnFJi@(Rv}BPkXO(u#Dhz?fmiTcpm@-J;geS45F+th!Z&y)LUmA8J4onFvTsws%v|4Hzh7& zKUQvDSsd%~{&WsPZ<=;wyMXEB$sqFznkXUPGZ9RjwwOv2TMXBt`c}V>#Nkq7U?eA3 za}U)}2JN>L%C2$@_*Ieyu?Qh7W9+4XbfO)Y|2GDSbY!|LMrhM@H^FPKkzVYoaGB`# zrD~i%Hd7xQ(}>b}Ew$8*^Z&AGou4Bb2%ZVV1G$(iO0R;L6j>IDlh3Aoem4++nS&p< zhGfUvR!oV2B;ec)KEODD-9Ve!D&5{nJS;H zc#y}EqD}#wOh7`nAc|*Bh^X9iPlz=0NnhL)Q{!{Zhu({Gb&5teHZy-dpwr2GOcUVQ z<(tvS*S;51J7UHVVbO*dJ=X+qS)6AARDAvfaPcG)Fuft3VM4&wpLtizOv>C%L1$Iz z^UjOCdydT!nycsM0uYP<6E0HvYn#>@{~J!@A2%PLF7qNs?;aSPNt-mI*UOC4zrq*H!2zfUB^;)94PJCNT-6i*&M{$Ck;#(Zy-8zhRHK!M{-$_c; z{Mo#sP?o!qqnlTPqnDJJUs-ae;Uz$?=egoYO}@U7nzc0bykepJ%`guKY;!4zN>{C?OYTW6rJJl>y@vMsHNFg5e5`AvSw5Ir{id2$-JylK)moFh!|CYdd;ZUt<+ngbBd|79mp4ks&`RQP=5-JxJS zp^<*EhsWz_`e(J-+rvoU8`);!6|JB3PwvNeiW`rL$uJYMUcyHdcqVP-%p6kEp=F4? zmJu+LjF|RiYP96)N>c`^D*f|J)lJ{TarZsDve*fPOKmZNEd=?i2+I8KpZ&)_|0jRx z$`$?n(ljUo30fsAUKRGIt9bv1Yq$S0}AsM1{LP# zIbUR_{vL^K`|n{wcyM``!Ej6D$WK~yB8OSt%?UPvJ>@d@u#^;Q&)G8Ixb)h&k{OZ5_N9hI0`KPhEdKq z+az|@y`H6AEH@+~K^*_-Joek^9%-`Jm@L4$G+o#EfiU>7^88q_66cwI8qg4Xwj`xYp*)l0KUbe zl3=*vIu&ktzaB158~85c4T_64UU_(oS64WBiyq zqkK;>-n8<=3I5sDfh&IZU_=G)jwXYuizf5Y2KC1+`HD2saw)$=M0V{GN`g!Z@rz7$d7)e6Kipdoi5ltWMR$m-DYoy;!9rAFmX-|t4`OB zGxQ@wRR=0+#8D3cB3+A`b`&{M(57w2!L{3(&g_$dNCpEo)!lZ+F!O#pOt{tsyG1;} zPu^i+$F_8+w&luW)mo6t{0BIGtBV)|;Wy5-@|awGVLHrTYh)pAlZ9puJEyNT$ueAQ z1+}7ucxXL%Z}7B7PEE};5wUN94xeJ$-|rQp&st#| zZ013Gv?}nL32&F9{p&}??)4DS?|S`8{L~(H_tKeB)ekwU9g?zee#|{xN13c8{nF)V zyzi32`&y@ZH5#65#183O8^5lD$!@z>UhkHj>ug@U+pZss5Bt}lXT@<_AMr&u%7q)7 z;#{aqvto;IfRCH33d}g7o)%9aJh610Yk%YV@SESb&We09Ruey+cd71UUBR|;nOX3- zdqbg9IRHati8rP*b)}>|c&3dCOLM=>&h7L*6*;9c!Xqki|9CH_cx!%&*qYpkcQ8}KZvyQfZ;ky6t++(d^9{$vl zJzzwzD_vT5zx%B_>Vab#vJdWmcU7k!(#0w4`rBLc6(<~qg&(*aI>i|qDbD(b%3~3N zc8abIW&%fO16-yk{<^KY`WE=a*nRgb@P@Vt3!q$T+k{;dZ!l!iu|>lhiTs zZC|gQbPeJ7A`)qQy`k1K=jZ(Ru%lp#ePcUEmxb7WT{Yq6 z<%)+)wTT*mrV{-AN{=8oLM4w$=$$ zgXO46Nj1A1tM{y9A|E*c(Mw22Bu&c5>M(@m^ztsEP%Sc81L5vgKBFly;WD=#VJ+A^ zjk{cdB>7+CH_vlSZHi_Y&1`EIbi5oUyGz&fx5Nj0V>NRRA6;92XFZLLt0Kt?xq|KZ zOS75j?j>y*Ba1II>fpVcdB%v%I5?P+U*)FubQZzH|9Vv)&GmU2nABLdyGklaa>u^%iQ#Fb;B$Y+YgeWSf!K0LACYy=e z^!@sW z6Ez)=K>EHe4c8yvn=Lx#E8w+M)sVi{CGV$#H&m3q`w_Bqqk2DpS}x)`Iz$5ZD8p3Y zGnessuLVR+3&ZC)c_W?(5xtQjE|IiGx@U}L?P9Qp5@{^)6MI%Sh*CQ7e=iAW^vP}U z&_p3WXS`X9k4YFlg`fABQT-SthVB_JTZ;*C_#faRQj9JE^0En}>lNLkId_sMPTUzQ zIw%C*l-Yhh*?0Mfe9~5Ky6z9abP(u#?|`0&XLm>(WS^GgByI)yGf+&&6#SI{o9*ml z_1f)NCw*~<`x6DdOVBK4N-sdtus~LBE||DjtK%2!NR#PrLe$JIp%R&;(8j&Jj3@{$ zU)ofyArVSTQGH^S^8q$|92#vY0UK0yz|HR4vL%JSC38vSJVofy?e17Xu?D-n3a$G{ zFFq-d!}<==$hlPpT~nwD+9JyOTK42nNti4(7*I895_=g7PuO6zK?*-%NBw?zB0jUQg&$_cSJ4hDTpc>zvvh)#n$Sl30<3)jv*RDfm6JR}xe+W<4MuFiz^rR%Lpq2P-0P{B2Z#E;roAj2sdBvDka%)qP3 ztkZ?dq}QfFfi6WmdOPz|d`rlIHXZdmnHy&l7AI7@;+zeZe{s7kny>;0k8z}Sd8bW-9FS{D$-vBjt8H?(ie z%gWQ6h3Bd!pO9=~=5i9X-CQfsK^FxQC}n5ok*|v6bq<`mX4Vs*iQC}7bJ{}IGZ979 z2ycw3$246!S8-iYQA|3bu{1Qlc{|A5pGn+bm%}So;x#>1iKTdi2HSkV+t-S*$z_z zf9#cqRm%71-7|^E^KuITM{-jaY}0B5qu+8fQ=iN%7q~Wy-+nj<{gC!b%fM4cfd+>aTfRcaoDnMl%lV ztYNWrlJ^7%q`NInQ>`5nU*7iL+7@v*f&T{ENx6P8Xoar{={9q0VJ?GV8?0J6AmVGy za@z{ZMaj>L{yW^gN;e0w?3+oq$e5BSgciJ=cXf5HHBlHm!li2;@sf&@pP%6z@=lw) zArfz_XK82oo%&?-Z1scx5J7RQd=X^JTwjwkGT=y)!=c)W!-W&Bwv~%y+ORrKeEf8} zMbvsDagk-^nD!}HUr+I>C=s(M%vK1$B-{LE{!MiAW3KAoTvyjm^l3w%YzTa7Z|a)} zO<`Fn^sQ>(VOrw+xUMsfL)^mgRxCPa>3p1g$uzvEjMMFbrTDFiJ zWNOb`oAUwP{g8|`(he=YAseD&hw4n8cER^Q63?S8ZIqGqj$X-U79L(l1pmB*2t4hB2(Swi4 zW7N++>ix4nELR_(J7ffd>|@dr^#yNQ9wRrA%3MmC5H0(Z-5}{OxM$bD9xxDR>KEw0 zl!-~7s_~?Ig5W@&YA1jm6~z-M&w^F(W?L1b?Sr$pRncgeXp8^=H>ihB4aT}B4gN4Q zgZDmp@ZjOon{8}<2UTHV-32QO{thv44)YW2vq-9#k2?oYKB%FJ`0PdW=%AeR;V6pi zhcXQc@#E*N%ImLV2Z;%$aycdPmp74o4-or+ydI<*g(s=aB*WIK&1~-mh#}f1Kkelf zHr_VIeP^*Jja~NO-ktAc-)?7jz~$T7cXBdxXK&+JaQoJ+RlvFpSu>LKQoX%QF5%L% z1xJsW7vi{f=478A3*y9Uk1seN$(~=oeS1wLX=802yVe97Hx_u7GBkJcW$u!nG5{>{ znac>Qc!FdOLy)s{!;e=qCE{QnuV_-FJG__$@}_R-^^%}id4G72@0bE_XBQ&B17#Mp zBb7&>%;@$JVC?R&%R-p4tYAbX5o)x0gv%n~LkE+%an>@6Z>{!2BMs}qUsT9w{#DB# zC1MwD2Z3X3p>ek3JHwmn#VVm$prYRJJGPW_4}gq7UHW`L2C=IBuOQS$H{DyZn0JY; zBu=^Eq02T^$UR+-25(9RrQg+d=r& zeDhL;6{|O1o=vTrvYDU95q#K~Iy4R}yRFBT)*@lm?cJtZ=A4*FaO>Q%?0r|fX5MSr zE5#07;aXf<^vbNbmasCRrTA0Yf^&PTyK{m5n-6swu$Fc*f8aFWfV@Nl?E~^XOg3M_ zy@WSG!%t;g!pTI3^D5C&dF_nBIp3R^AK87NgR)K2HnXu741&Xqlf4%z_{75VJuTd8 zt;YU&Q)*|-gyexhaF|fZ1rXQwVC=Ds02F1`1{Uy=uRc2KGnVVY-iy{ZBdTY{e`=6E0*3YlbJ!&kG%kKD)oKlV=H*lrOU}XI~wM-~aG0Pk!MRad!1Lla2YH*Wbl&=mlnD z{&m-4BLi{Vt??5QQiA(rs5rzJ7so|MscQ3VUOfhjn7?MM=h0^E9SU2i~d zNM$NG1XN`zkPI#9o)JsbFOFf358+A@y#Q7OF_E3#jm(pJDGLv?_Av1PS4u%Qs3c&5 z|*U**)LOtN&S^k_|Kuek(p5(k8a!1C83h8S-e2Q1N33+58*D|eNv z*B~G7*W@8J5Ru}fl&{vw(A;nR@Ay5N<^~wXrloVjJxp6NqN8)MW_bF+6JbM!a%l#S z7_-zlHgR}G&*!>R1F0LH7|A#lzrD+tpYOi zEYdSBpezWm$-0ChbI3DLQL;yD*2EUmL#e-RpP=n|Qh;g|h8*f$41l^_kS5eXU5p`J zPGlbNf4>H?#?AqIq(JdLcO1dZHm@+3G!W@M@SgUy?9K?sK|l?pG``)%B=1`GX?yki z?d(^wZ)23yc~T%Y0TF>Xg$a#~>RK?y*_<;px$`hkb0XI`zl-f4^5FyXxzsS+-ZPQ0 zS*fOdg7%&xX*=B+o6u#C(S6D!v1s6BPkY4?eB|wG+0D0a;e&hg_S?jszWp}*WUNxd z#|QUaIiig*AbHr&2SX5uAhtyw{GH7A$@8chelVA3bV_RgO%G2<{2};RCwPx4i z*jF?d5Zx$y2y%eau3N|E0pe89|5VrNxp7}q)%^R1h2F|J#; zHf>y+)8pFwjP&e;2(VSU>j}vHQ{X`BkB+$pDbCm&!&&2I7w-Dx@*W!mXCTOFwy}w` z^ILDN$=mrY-=!OAmtrN(hXt)DlK}YOym0~#KI@TZ-AS`Oi}J|Q<}xuQ?=J?ri%T*d z(7Wy7i9Tub!sx`o`hxT0Va*nkr>W{VPG~|bb-K?J<{a3YERM2w4jima+1M;!&YbM$ zLLOysuG6)%<1bNeFd1VRw$UYMfmYed%z9m~9U%=v;bJ=+2P^^ez!Qyem7txKsCK9#l>woMBV6g7wmhODXNmvpv_^N3V*4)=wutM(U*$bU0)ft;!D-eMM6o&{vqUJsqBO zQEEVza`ghm3d5^*=&C{l)sxwztFb4B*1r9A`|a!{GJ`r57rlTLy*IGM!?=noRTqht zL_RW2xm00=I#29YCWn`~^`XV*NIzuaNoQSf7Kt8Cm z#GAORh7j9`MkIFL@fhB#NwRNR8!Ua-*?G8dOT?aWjW>qQmtC}SL~{3WUJ1}Ozgy1;f30#-r?q!Qz;Uw zc5h>kn1)A-ArTkX!g56P0rt7R_}f1cqy*Z*^~H066|yc3xbp)3Br(4tTp3x((!FTU zH>tHgGuq9m4EF2g5unbY2Erh5D<0=2gsZ5ba6~eO#^_ zN5aP|aH!nIgsoGAoFnLvJ4VoIQ6*W;x?_`RQ!pp`$-`O!34txICW$mfA3Y1$q-oC6 zNGA(e%bZuVlA&m=-~hpoofHd?_p`HPvg@UqIOg2UglsE-e2K%SpNhdnz#5efXEy0} ztTm1f-&C~r%!CMJ&C{D2n~w4ur(IA4-~n=ygA;87&aKxRK%WfCBbu+WH$C_JrfKu* zAV)RC7Gm*K6YQ7rARItm7#@Zbkr(#W*Jkk(}=;F?fnU_9Z;enfDbNaQ;PCJvv zWre&MMCKwY3LF9O$t^13-eh2wrB1Bd*|7)P7}$ZLq>aaZvP|Pqi9?w>YQ+&IxoG^} zfd{a-jcO%S4)m7G2~>v^Phvv+2@K(lSgR^a<{5+Z!M^OIJiNitGN`z2G7%KXou*JZ zs=?cfy*!jz2e?ldm1c68xl$7oC(@fPz!n<8`8uPz-ZO;;NBVsi;%n|_!bjA3J)Ip6 z&D}@|eO#ef(W`flplEdGY5F{Sv=@PcBcJ+g8ZEGWOAw;0=dgi7%P^QQi5Y*9`A?QZ7EfHWZ&xtzjHm>=A2uEgZt6R4Urnia) z>a-fpSyPJ&m)o=OqtM*-x7Rq%`|t(C4ux7v!8J%Qt6c>Y#Qw`i1AtBLPYS}7vJn;g zwN|_Mj6|q`fW87=uAbaOPJpm~EVWGRtiOACW`@VyHnC%x3oH*quQXM#+ZGB!^i%pU z`WX43WQ&>zz(uJN8k2?PG;>Amm~1mlbCdZ;hjapmozz1|FB1l{jq5kYoXu`eW;VJj zKX0p>`x|UD*RT{#4G*({`miu4RZ~jmuY4^T4Jb`Vd|kz!1r+9kre-M?FbvK*9$ zCHloeQUAq!9mppj{1O)am)1wk{ylYgDeeCM!6w5t@l5h>(*+)s8hu z@l=tsU9ev%wPUS2JKqC-Q58Eo-zRdR1u_l_0c&LF&}d+kUDXN`U&CaKa8;p5;p|)5#kqQztQMlNw?kvlzS+@B6nv%LXph!DyMQg>Z#TR<^9m&ls*|pC{Chn@EXiOk ztLIspWw*dW&pOXroPz+#EM-b=RucXfgtmGflWeRyl;?Q27+=H+ z+?S<~VdC|PK;^?4j6Ec6^5w?MBSOLzV67soWGeeEVR;UV>OdS^tDbJQZJRf|-$g<= z?Ex90-4dmt`(&FozpnK~b?L#}x{}rQt_)M1PMO*0+W@eR_iwQ%Bl0@_s)=?*9e%EQ z-fC2n^wn9bh#4a`sf277oAlk$yB`Q8eZ2DJcu}}!#K|?9`K0^qWU3kX5jE;OHS9~R z#dtP14coaNB_{ZRsPq za%UVkHk$d_G&l3|voQtFSK;P@3gPvhufiJ(D#V8Jtl`DA(8Y*H2_e=YH!2!0w(Hp%OLv4$X1*Yjl%9E7M#Ye$ zsS`c*oP^l1>w2bt!tFDsXa4GWb1J}-doPp5xm5hCCmgiddY!A~Z1^a>T_3!y({iXM z+UiA%yE!w@o=KfE8lgNE-77BY3!P_wQ_tzy%QMe$85hm?jp-S8BF3Ci)$^UKqZXgd z%v7A^9@XrLooX)B5zY+#)+0_eC;b?LlOJz>Hr*phMv81JV-m%yvrXE**+6z(M^Gtk zduXR6Oeg7SY+oQS=A#GWZ4xRVU8LBu0mSqP@{bp>BEy#;=Zt)_?c^{?7mO{l9VL zivE6i3NH4cXkH*}{0S-5WF~V89vnAq{PYa3En(zrB>=H5l(1Y55?Ou;BSSA&mN0Us zXf0vn?9pc?^sJypA9V`c-X94)!87-jhLPpJbz)TMY_tT z%LI++AW~C`*&CrRM>TX$P~6yCiW@A?8eJ3yZ&Swo6y&Sy|1Vg)=JpIMUPYc2i&xnM z-A(ZHxSaE`czZKdY9#NTxmnMzR$OB7ODuls3Rq(C(m^~Mj=#j>jcZt9@vBQL{@fvx zmRP(vZqTm<(D@RJuN44GEZ&DbTw?LMaP&{mn$JP4mstFYS7}CRBYu)r;MKt5k%feo z8a}(kAVdHc1rc2m(2^(Y2U0aeaO}Niv6TXCv z2La|oM?wuVUc$!C%1e{6G~q_JC;9iF=*!Vx7>SAdU82cgB)&Z$-hF_xc#`kiqoF%y#@1f@;j~AUay-H`!a7a4$``Pt*9k+A#9cguB)?oS!rZgZDS$(u6x& zJk(e`bF(zzUWQYF@n>_0jlL{Z-x7@?mhwl@zXt`0UQuQv<~7) z)<_~s?w6MIg@lnlv~QlBUU$~IOVsXIuc8cC?jeEV3f)aY$B;+ny-yxIc=+@tIUNS$ z1L6}bO0C?Y*i}Q-W4#u&Wc_xXX?cS*MkIXcJtaV<8;KeuUK%kDeDWjmOI7dirD$X( zNTPe`QAl85Sqlik1e~v*?~Kt_jY$M-zmuznEJgbzNn~7Zo-jL%J+l!PRg#9tCroID zi1g~$_j{v?WW+8Deyj$>DC1 zgifOiQytez`@G{YG~qVbEfD^p1Y|xW>5*rD4mA~d$K2U07&e)YzHYg1wSIA2ev)7` zDz1_~@n6fy{sEcMAu(j!)soVQij8GizQY#e=iZv9kqs#AqH%sVWc;HTYGD4t#*R8|SqWIecr2AgY+608P_UDAPX1Yiestpq7BRyX3 zbj=--Pw>Uc_4bUP@EbT;m2N4PgB)SjC7sH;TVu|C!NISyO2+oBms(A_MnX>|{cLvN zGWiryn?AQ6QN~O}YeN@XvF}13LKgasL`Kl>(<$kjZv`Sc^-`FyntYNz@%B`iCJCYd zmuRd!x$&&?J{5_-@uqT-#Mffv6{8C!ZT{8^CYhCs>>b%kS6}3&MX`h$SQKCCT^> zTPg9>w)1?sji(yubvsl79(=8+4-O=d8R~_R39f~aVYt9qPlFf zZ>D@*GnhP4k#AfhdC5N1dyD@-ZF}AE!8Yk#du6rti1&F;skQdTXUSQ~kGZ`W6tA=u z6Pe`m7QK1DdITpqC#buNK(A=#GtuWa)9CY^{Fr*Dt99MHk#xXC@uw&$567+Asbgzq zmSJ-(H{~jID3-snT>A`u3i-c1Qqcd-Za&)P8LcO6#IvHj4RzL3o#Z6qJ2Rm;>e-ub zx<1XXfMgA)RRCf=uY#UZ?IhjUIK7sQ`L$>^Evhx6Z1oNqt{dvLtR&wv32Mv{gj5Gu z;tJnt#J!e=g*9O}YRwZ1Or4(mIE|d0jNCs|+-|Z;^b2Iji;*Aiz>W76yZc0p!W%GS zor9wDL`i+@T8%W&?g~=lWGzR8UuybiN70Ug1q^qmt7`FcrKDOd9qdMq2sfrv_o?JP zlG&gzx|b#qT}hVB*K*osjkVlm>x>Yx?mJPiv1@!qr|s8R(MOnmAJk}IqaHR`zaRxy zx0M7TT_)6gbTy zzl8E}DOp1KODKN{<$J~M63WLWdI{w(q5LJ3Uo4^gY-I`MdpS@r;u6ZQ5x*srzl8Fa zP=4R($(K<663Smj%Bh3)5H{xpN6Ptc{`>#=fBXl3{mK>nWkk=}BojgpOB(E-Ty*u_ zvK&{ah=%PK6D;RJZ+{@$$y^qKkOwmg!<6W~kKf01H1D5O7(zYmJq~+Y_fyV^v$(MNmJ8cvwd8Z;iKkY;_OQ>qQrO?gY2re+b?iCATxSK#5ft8@F18m zcMi%DOIBHOr$)if>atSj>ha)L)!uJow!XHO^zV|l=U`j?V|{O~HI9rIXm@2bg>_`S zDz$yz4HqRd%T|sf+;U%%9lTq)tgJG$Z2U}MKEG=zGP`BdJ=p1`h5hID4<~ZQ!c^zc zT6=ASgM=yNt(#aBJzg1bJhJTd?DoyCpK;U0G5?S{7;D{@nLmX#u)^eV(*GaQ)M=%8 zoRp5O^s4Q2vu}sa^+SNgCZ`}I)@rri4!u_tUe{nVYW|pkXAbs5#snd-qlhtM; z?IU(((Q5uE)98nz33g$f$ze4)Y*kNs!$Ix|P#e253=oN**dOGrG|O%>sNo3C0XwPB z-}xf&7Z^^bqr;du?+*`h%^t2kZP*I#KQ>HK$c;h!Yf@BK6Qh|gjjlX|eAGMCXbOgo zouVuIZ=)u#U@%9Ujjn42s)pTrTti1Bob~$co7wy=76_L(sOxY~<_mPFAm~_n#k#Q0 zuemln9{Pm^j|?w6W(~t`)$;Pt49-ZrEZBY=;P{<=6Hd*7lNm=(MfLqcCE;>;_6lNE zZhKq5M6w3MwkEdK=qg5#qq57E4ILDj;CP+MADsU(doVfVC9ZV4sn8(fj!>CVqBmZB zc~%qKUHlC3@$&O_SX2ix(eM<*O{v)%yF-e(tP!#BARmX*8&^Z7#H_-{*vXaq2vi5X zPOCN!M&cSzy2X~9!|ES{ZS>aa+bO-N&1RvBvVs*=e7qt^W$|YnrCm^~cfVh{KAJf1b%VWKouR{UyF%Cyl_n98r zE^RhtKeKHU?+E>*e;PZq(GRzWvp3lwH|M+`VtdmCx>L|HaS$jLT^EH_2%C`X4{|Pk;P#SFY&qZz;{4|K}$^)`xXx zMdRQ-@ogm9+b+sIXR~+zH=1?*+dum8-@8o+Uj6-YsEGgmp9MVtQd+|L(8;-4DO_m#_7hbKlw{nuITTV8kOZ=OjSmQ zR>vZ@^Jo4~R{r7WpI^OlMSsni|Mq`7Rnz_v;k)M7@{f8yd+(nzF#R?A^WGm$)uR4X zi~I8*{^iLp++u$8*IdJ2`)5-%9rXG{92_s+<$vdoe&1EG5YO-}s9^tM>mo{~P|rD_4g9;A?#Q{{UTEFN6R9 literal 104784 zcmeI5U5sVfb=S)TkP-+1J2n>B;iyJtZcVCg*RAT0>0xS6PkUz4u|3n7ZhM$?HC3nX zJymtw_nvco&%O2G(FiO-BA$3aLP&rTDaeZQf{?68JOxB75u^kqRunKVi2%uO_yqwH zLj3+~?|sg>A5~p7J+A4RzT+uR)%n)YZzT1rsp?CkxtZ025Bd?KMzy>*6{Nz$%aK!-~Mtke7=Z^VLTW{ zJI7JV&B2h_;=4(dCGDgr;=ys$kK2#py`-=vyLsFRk7)ESZDxypMjqOK{Z!1J#dh9F zsu@gQ&n&|4$pjjlt zt-P-h6^%|>^t|W0Yj<%AcSY#haYg8n*!mqu3 zYvoQhghqNL`9)e&h7oOVCxbZ6hFhb9Z*Fh1%9nVZ8A7zd+BL%_Yv%nV+wwzgY%W$q z-CSJS%y|BSmgv5=ewgxho(-9DctP02%e5W6Y}@Z3VTZQa-V^&diXm6nMUq6Fbayuy zu*tdB%ZvU#3*3&n$wAUBn$es2P^%vg?JN2;e4~2|YoKX9eV9i*`mtXuZ)QafGm7@( zgCzR;sNfB54%1!&%Sie~)F_Osilt`s&VJfevy3J!vge0Uk{RQIUBRse`3|pdW_8{d zh_OG&+ey-i^kgUD)5UNQYX?RB(V(9fiS_w#nE{Ad!NGRY%;miw~opjapxiUhr2AxN<>M z{bq7B4C|X_>5yYB{bphdw3|A@+ursEA~dnR4as$qu8ufPPuma8=>9%jb$k0rOqM77 zAq2C2#7Sh7mc+%e&aJrH%~|EWXvB>?D*8z~-90w?b`;9$LbX{Ts)8(xz99{E6o)R{ zUmJ=o;E&gEy7>-HoY--DkQZ=$j%c^tW^J-tJfL$Gsv2~n{vdAijgtRz254iYd;6PN z61VrGJx=LJLWoAA&0leA`qsKX42`iq1EwBXvt{lO?< zWX(&>t{pMJu#={f-I^7xG_oz-E;d@rEce*Bm!km_p7*wx_+hef14eEwHKDo2FD}Cu zdmEdJw)k-Q3)%{I;ActP(YhH6WoH;Qcp+fx$=guZS=U)X3ZH6=YZkWe8sJ8(^$peM}&5x#L@ zL^q7I;6>C;A30I4THdIiwz)Iv!=acfq`xd0Bt_EO=?2Ev&D-t-cC^8uM?1qo&Wp$$ zUxySqqLUak1my&!Vso^*y3<-;wg1LZU`@+3FIi>OjR!DS-pL`rRO(K`=V7)8EVSw* z6i$a>w6?k$?ca@!wht2+J@RXQC>Cl=aWvSC+X=!(ZD*LcZd;?i8KP>u93*-4C`tO( zhCSZyCh;JgPbU0mb=-aJsLkvUnaNFDwc!wL;vIw}$qoo%2x`2+%ue{cx7`jyttKSyLtVaQgnOiV2}^&ysF(og$+%DK&U6>x40!C95IFG*QfN_ zU1V(3$)S^E)K~P?^5w-IGEP7%hCFuVXl0`P^21ZtchEhgT1Z<314ssHee$XP@g7pt zu-F1YaZmUG)Mm2k|KfVH)oNb#^Tl9n&K9#BK^U=={i{;W%L`k*G|Tnf)t`9(Zi(pJ z_{`$F|LoQeKl4YQSy<4|=O+=J1o?0aew1}C7|?keW#vFDGlPX7n#C|Hb)!Jz816@} z)s-Ec9H#081vE;#T~s189TR3@Dt!+2itWCpRqmh`sg;J>l716VBZSGW7XtBhWYr1t!Gp_SUU??L4Pji=h zfjN}4z@-AqKr7946-|N%OfO~8f@Z4qG|8uE=g5`3uZP?pQa&pE#owwfT$pCTlBT_0 z(n&dcB-8^pmRMpXb90RJrBfg)*plvb^BqB!@qjZP>w=xuo7uYwV>bqo?#W_-gk9vl zeAW-h9y-SzB#gUyeA_?C|qoC~AT1956Yc$~W)56&322SSF~ZP(wqVr0QL6 z;@nO5fWD5PBGXJqJJO1Z28hbqGMdv~0a&-uXss?s*OsIo8#zwP@bS?`eDq-T;TUP2 z^jPcRMJc~;<|SA>#CRoHsRJcCL=X})>f=!1aFQYj_HE~tLlm=)oj#shK`KzfEG7l6 z=uSr`?CG3Mhw%wXTVO~bPuVWOcmFdS!L%cIuXsyqG= zD)&>6Li+DNVE8&3bPh_?w`0AFe(ko2lS6G-=tArQq&qQ-5x^MAg$+>xrvn?u252p6 zPnw<~3AE9T_m+6wwiv54c~75N`ipTaE7CzU0TVbc%8I?Zxpr8@~qc zibujB%f;Scy_TbFtKW`2N_owA{b|g;V&YQm#M$+9b=43m+ky4_;Z!(>6#gWZ#QP9L zjlE|`Y-eEcFmXt>+-aad&miB_s-46XCwS+?qjkbs_0Xc@xfjPw#j)U3&SEU)sY7s# z7PwN%>O3%^ggWjZqYlmusGOYg9}D>ovCV-FCZ=h}rovA!NTAbag&>ZG+oKbUY6mJW zcM@aL@Xrki?oNF7O)f9y;3L*jyH=JbQWgl3^%c;F!7AYMq)uk2FKBKt3Q z%~=nD$G$Ml{ee;Eh`52~&j2F+7w~cXzyI~`t^V-S3k&-B%d-&i3&z3?t{!rj)z6L+ z^!Dy&HTNDAYJ24bgj-M$o*(9QhwDcGC7qK_a~xVnla+;B@BrV(8Sg7248IN+NUqo) z?fp2zEEL`2sFO8vGrNa{XgC=4qFWwQr?TgFn-)VUz-mNU9eAB1cffhrD~B&r7^Qe4wM| z1%)Oeo9f;kqnlR3417s4`c>$KgqC)4tPkEgXIO`z-+aJmOuh8IYex4*I|XPzs&V(& z^&)6?>Q5)J(*-w6jr%xFFmkHTc&#KzI&GX)aR94MAm@JN?#iQzo0yWiSV?IisE*Dr zK^~%g=w*P;A%-uFvpr5E2-;U|Y3^PjSG-xLP4&^7S*Ohr@@s zX}XJrM{kHP(nWkUqn%`Um?X|bQ0kIu*mu(TrGiN%1^AYCN7WGI;vxsYw!Fd|f6Yd8 z18&VI!Q5BMVfOV|^XmZI(+F1Dc|HJ)lcg1wiVL578$7lhMqC`@;?-krJI5II)3(5` z*IAI##A}e-NUm&dMmDbcIjPOZ#?{sBZ3E`uuQJpn(bZM4VR|7B>f8d%0H>K88?E^e z=DWLc+XgScTJzDxCobV+a0ivyf_w2%+T-La_S0Q+F0@68jw%!jN^f1KH_hk`0KTkr z$)x$`IT^06Al@@Nzc9=;i_ANGR;sU^LL#m>TD$5VwNEiAcSL~s&Ow^XQPa5~Lq>;C z3QHjyqY+@3*K%3H@!`)Wyw_NEe3HuJi^X!Mm0rccRFK6E$VCNtntvw;l3XVaa3RwA z1h4qbY;Ss;6?MuH`Vt>Dur8d>VuDz_Kgk)QG@gVRi8_J2ocBVoA0Or;xgd% zXEIHtVg?Swa`+{C*F4T~;Z=ib2XsqxbNBYG=*pPH=Ln6*@J08LH)AX{SS-seR%^o> zO`@N}kqiQIF+GCOwztReVad0G(7Z!gz}5^#zXf=T2jenkF^z&bB{re&Ob#}WJhfq7 zbAa?r`EaQ|6Yt0kh2HZ~Jn=ACzzDIh$>>I7!S9}r6LNs7GKYdqnYtK+bL8LcZSO{F zTwLzkrHl5Z-WTQ*%IwtupBl)&nJI9q^fd<67uJ8Z4#{4fg8HXCk~cGOBm0mcMQp_r z@FDu)9t*hsgtgb*tfcL6bMOW2>|3CKvT?C%@C8XT1P$~wGylxS4uu!9$9OU)lEX)J zNS}e9!XgIoNdiGQgR226ZY4`-t&0F7Z3WRl?j!NQIVM5e`*}ycA;)Zgx!v*Y2}Lfr z;iG`vr(15>~J+)E_pw$0!0F+x~S3FgC()HYRg^>QzaKS@e2 z^Q#}>QryX+!vR(UotLA&Eyh7wj9hq+7Q!Hg(t>TLWx0tZ$caF38UnWYB9{mV@U zjP@)7 zM*;5qV8QApp0x?Y$@|C6Hd2vq=B^Zm+aT--l9WA!TaWy}6a}w%M6uZL@EiOs$Xo)> zbnN>%Ci_Vg%iF{pp`h%V`cZ~BYjVs^PykI_;)RzFz&OG<6{5}~<`GPEFTWDaYk^tl zb1mNFN;8#=fEUDpV8P&{G3|rwdStQ zey+7ndy>_cF3R8Uu2;<^eh_DiFoJ>+A|);{MbvFahBPzb4-EAhvKH*1HH zP(=^_GajyKwqxGM&_&SKqY2Gl7vJJilU-m= z%v*6sVb+a9L+dw0Z#I^2f;HGN2+gjs4gQ0+k#>grPRlJe9*E7kf92O~PSOwb0rk>< zq_I~#U=O7aiEw2g*g_F2H}os;o1g@#B0%ALN#p7g+c7@EA(Twwdt)p(d1JmDJw7>!Tl3hnB>L%f`S_jGFfVA!{M~xVm18}@p zQdbkqA!c!Cd^7Dg;v;MZ)z({5G^o1%h3xj zIKuLe$FBRt2R|wr517gang{nG4ea>DJL5yHy)iZw)%7qH_w0C!-h1hN`YlR~hmBMz zGDJh>$ZlOa4$at%X>?k*VnkLSH>F>Z3Om+{EfNnJGD27sFEI5cVXuUvsQLB&9&7K+&uK$qgB0w!Z-LyR!aA}N_I z#1V{pmT33UHceU!Brq~Sa)$u%p)!(L;sTJ|pUV$PY(zAFXN0{cR}Lg0&WhhB0g_N3 z*I~WoER0HOrLNV-YKU^xx61RT+r^;kfoPAtxqXv_Z;zzVpD&U33 z)-}uMDj1*0sjyfMMv%MW5!E@(tCkvbEwrOpPr6$5P{sy^4het&`^PKez3Vp?0h&?z zE7q6D6qZYr5Sdk+wuSaXnzp%Qut=C(NSAC_ieDSZ=L#|-NmXLbVPTn^MZq1JzL(_) zBs4hNnx(?X_V@#EbE1XPyTWV;YNHRa4k}~Z-*D8r4^*v87d&9es2nl6X>*q~LrpD+ zpYkvk-w)`dr+h>c-?zr1M%qlsOktlULx~&$*e?pq%fr%onViXKMt3}ls=35QBQVxP z1*IcW58WnknU8etnMZ>5L^$YPiW)7>!r%ow-X}zlHD%=ag9CTL9hY0JPoJ1}Va*-kD1(0Se zZ{ZoBpk>$?>D7qVhVT^{{((UuQ7c_Bk;>B&GEpBt1cZF(gS#I9@i@mqJ_*$HafK8r zu^x5gf$$@_pi0-%AJ1v$M)yM2GHvSFyV%>bdF_l4gR+v%L>OQoT{OF417es06{_Ag zF*)JpIXNIZ(rfU=HMG9KN6VM2 z`{IOMJ~rKD*qfR@z0={L@ovKk9Tz;M*Zo0$_FiA5hx!yRouSuysv1wj5JDrDuo|JM zvey^1I|C4GfPkXIE%#q_I9)fCrZ8Px$KpEjnx&BpNh-rBrd5Mi7V@q`7C^|r-$+7B zmE6zXw1dxZHglKu&2!-En*mCL5QDIm9POt&ILdU=evSa6a}Bvfzpt3cS4pp5#lde7 zh|?Z)D_)(*CJf-gQ~G9_kWDgIi@(PaLv&BCixD2;^2(4*<>J}R@5jl>azNq$N6n`@ zB;_KH77~CHnQau6Y=@w1hwV8Q^XWP>pZI*D9q2Do`7j%)vg7D@n1yPfGV$j2j1t1!6JFs-=t zxChoA^FYNwNtAFRI>x9b;IZ0!t&`tlCQg2D?d12E{MvhHA%eEH$_*i3=?2{06vU02 z^zI-A?9)b)^d_A^-Zfm>A@^NC9+`WMQCdR-K$yO*Hi?nnTnrM0-Yj<#sUx&n_J-Vs zZym|c@r4)aU2#;+TF!3wYtEcJ)ogg-^kyNA+0A+mJlWas{3)$E1Q6fT{oxjF+XO)h zPaj7GZCX9XV~`+XPe0ut^WhvD6lP08xMwlV-S2}9{-+=O;eY?VKew=;pXa9hHENM$ zD;Mv&aiTYdG`SORB1l|XG2$|MVrVg_H1T*xi4UTCBjB0^^ZEzWjw?xYb60_R*6UEm?7yH69RI`b|M})$Tv*W0Uz=11 zoCcnBJ_1$X9aH{iI19|Yf&DRy?yU%XS;Cq`>@i-dE-tLhOM$|kD))w2a0{|E4g_PC zl~MY+EcRwJB$ajP`Rk^@RRh46$r~s^PA%Zs!@uKZ{%8jvm|7fEfx*%t^mZzvhOg6D4nP$hd{^;mK!k zMVf5WfY=iq!p|YqJY5Lwg;+J;GS?4+l?W(auKPIt=q;NcAX@}J-Q3g$GkA?nj(8cnIv zOg|R18y;~+&f36gya|SP;y}g0KJY+g)#cl+-?VhDaq9VM4rW*vt7gZ?36HIfbCQXM z?!aEpIeZfW=lJ6KMSitjwUuC8`{^v5sA5 z+BUlKSlcKC`0`o|KxxsD@BZnw(N*lCc=el!$DB@_8E$h?(%9`X!E{k;UA`N&JEe?N ze*}(jyOO5rci2+3JL50o@PeN}>=yeEL%iUk^`= z1VGU{<;B)nV+v{}(NoN!`k1qWwu?$pf1FXX+Q|>nnhh~jH&Y7YX~U!2(a7PNgtIb?GHZ+%2uiL zo?g{IXL>!HRA(_fI%Q^VdBoReUSO+MVC#2)8ovAOmEG_D_`-sIJ~s(92%4gfFhwse zwqSbQMK5aHcMyV(<*^%}3+~ngxPbe#Ww{I3!ZWT)k5@21$rJcB%FcFQ!$EBS!`GPS z-o6)I-QF~)leHku9ecv0{nG!IgBA@y2O+=U2EHgEh^#>DBbstAoSY=6?g3jqUda<< z20GVP&;`yRYyT)!0ayNiT5e2IgeIAtR%aUq&o=L~E z6|%{H`!~M#?|!C1Kq@~|iGubHv+vx%GU1alNwkh6B$ zCrtG=<{+zsUZ!`C8gIPA%WP|ls`Rp92wB4~UfBbkR! zaN@V8g1LuGw|&bXT=zvP6TcIKImk&oDw$L%aQDh#;*?KpC((C=tP&wb#4Hh!l-N{_ zF4>=yn^XocxRSy-t?VpO52~E0EQnmJ+MI`n7>!CEb3rD)5jQn$1f(v=X>k?<1Rd1sj%SXo>h$-e`@> ztY!NYO4iKU*9%MaXHkVkuH@2AR;8LB+OyWg74p-Iug2sxZb$d1a{8e4aI5K)Qj@EZ z@CdCmfow!BgMVu(%PIn85?SIu+k<#lnW(hV{Sk4j#lCy|_WpndULB0GqJtIh!5ZyP zmPG9_ak6E$R());B2gRdFdV{lty;1RmPH@Qis;7m7b8l?(LYU~<%V^3PIr-ov%p`Q ztv&`JMU^6s8DM9fEQ1)Bt&o`-h-WM(bx4lqzNBf!U?F3qTeN0V(-pm)5Sg7pB)6Lp zd86ZQGfb%Fop=%=V~<~A;fA&}m}eD6Va-?~zxeAw?3y={4iQhwv~ip2aADXMM>DW= zp|yo%uG*)sC4~uGthl~%pJl#pJ$NlUNU7?m86|u04Mc`(*EUC#aa?avc(#H6G&%4# zGYxL-VYELS_KTOUTp6&ah2Z!rZ7t}IU(YKyZ(hGa0u=YZmpkbJVdfw&Xth<>q_}b^;@&w+0{sOV%=v~wW(PF zsavG3b{)kYy16Utch}W9_pSL}4Eo1wRR{D}ldzcm(M~hXuXOVE73hsXqAUBuUiS(t zCOK;ASHg0QnGiXo5PUu38nj?-!ii+;HXD68Y`=Fs-*TO4zI>(s%F9=n@mDnY=biey z826z7SFnO$rstK}5nAo5Xt_ytp=?O=i*AKTuiTB>kF*P5X3|{JHlbNQVKm$giV_~= zXY`iNhr@jE$nlwFGSoWN+1k4`7{tfTVZLQ?wKVM$0a&&JHZ6g!imkTQ-`=izp#7|h zCc0v6kX-$Y6VXu~84SZlh3~`zMa>67d{u#85F1%|{Nd%8tuNwDEU?R^PcO-dUQdSZ z!WaV;L8}Z?U@=3&IPVyj?RBwHpS^$`uJkdt?|xMw^vog&9pI~?cZKhx->U{Dz^;*| z9<``+MYfZYo0734o!7Qe25OGF*LlFA{W+3cR`#%|RO?cKYD!w&?NYae?I9x98oD$A zwK4`i)XHibP{M_HIOPoBi_|e-8@lwcB}qmP-gw9K&T0odZVUPsuK53sZ9nA|@VkrT z<^wY0=sdXn)i2)Cm)*G=y>a)eE83;76n0I=ABdWtYGOI_xd11x2nt{)U zJd5rns#TPqw_T|&9bHuZ{x^Nr2c0FBN|*nRoo3O8d61FM-lsH4q6eCsYGc||o$YAy z^T;;1#2bO?GE(Qo6cE)Kg!c9cBBOkyXvpFjT&h0we{#Yg857H>!cM~hn4MiQlbZe3 zbY+HbzkTaXmkk{+ybjqv{HsmBuv>uuX`vJ~XyqKt~}9`Oc0N|r`c}o04wKxs7}J_7j^c(yb`&K50+amwUn$kuPG#ooUv|z8q`-jfKx1RhsKsB@s zU`}#3CW^9<4wN8<8XPj};&*8Wg58wLV)bjFWt|+VP{Y4o^4AP-O58UU$IlaQlxkls z1DNtkdl;s?4dSo!#p}KjgB)ZA9*35qLx3_+uQ#Ej?#*c#5NDjkF{hiPs<&*_*(LE*6xiS;^UVPiHz+jn*&g4SJ`ris=<9fgJx zm{e64Y{pAW3K#gXZ>ZiAoHgU5_TR9jg2ypeB=b@-5}npg1k@3~lNOku2YTUAf!Bvh!Sl2MO(stIEii(I-B67`i2jUW`-rm~nYnQk{f6dKqkk)LgeMpTMsyZ8*p<>RYJ>xL+nR~5P5)QZ5}tBUW?QxU-6U&- zKY7=@kcFLO{U-RblC)QkyZc+pf5a2-CEb@ehgFFR#@RIwX0bM9Ne?Z0x8`EW++o5X zP!`>H(C%occKHnZjdz+P+?U#6JyECln|98Nfd zaveY9i>~d0io5<-n4eG(@WhXv|y3LvRbOwTF{{2?^q**uv2j3~5rHZgc#Q4~B69lp3-{ zIvvc?Gu9d@+tn|q1ODXd9-FWQRUBL2>KQ0)fIz*B5O)(F3XM?YbAUUSkrJB3)(_K? zAw%{A<%t~a?Q_0dUHyU-%^Ry<;QYA6=KyfGPzqO%gPrBXl+!#;lY@ajIL6+6%n9<8 zvOzzXamMNxmS-hmuZHr7yTbcmuKaS+1nIu9v(+ZCOUB4nn>UPgO>he!^>M!-h^lP0 zq{|P=Gv*-Gy+&Bc9I6(N9v~s1Z5EiJn_(FIJa<^`!IUwiZ3TNDE6MsKE@qt(^jR=#o%D;mAd zqPM8+(Up^$6$5eFmx}a}<6lcW;haXEdQR5X7XrHA{L8|zC-YjZUj>1x#~iv8f_19f zVUj$8FP$*+-r=ZV#9a^!w-Oo=lZy!*_<(lk6v6?G7{Fjg`PbgQcklkewI+_ES&>+d zX90BrghMisUVNGBURYYjDtO{UjfFbL#6O8Gh21Y& zKVhTPXuoyy=Bn5co|HIM)=*~?ey6#FQS>-0knj+h(BI-r<1QHh;Qc6>^X#7WAN3*+C=J zS3V`a?7dc>JVs=dZ{u_}j2h;iENZ!LtfwV3L=V~dv1G+##5BW6BiITO97fx|GrZYK zR>+_NH}+o2@O2a~H4<~4fM3A9y8MDxQ@;v^Z2Z@~C3kLD-9`$T8y?m{Og@$jUuTfB z36=a?D4_VaYVFtof;bjVVpA36Y-p2XM-T(0FoyMC+eqE>-Uj|_{_`?MB`cS_z`RcB zE(Bf~VB8^u!dq7EM0v0^+|c>~4%j>6k%P75Hs}qU=rrjbE7b^Tc9!W&D4u($RP;`I z)IFdVYNAnSxw@nWTtrs+*$WOWyszDYfANGO4{GW3DtIekZLdb0H7O~)30{4wfKQ*IHW?Q?o;$`jylTkxcn+Aegt9vQn54Jyyud08 z^u8-?_l(A}eu`Muzk2=ipZ!hJrRry@!0LY5-61~bBPg)?bMDJyMMs6xqo-A|lSzcR zMdEdGl(bdFGLB}wXF!dGByMI8*2}<0q*pnESU63Gb7DYsM1e9GQzkN|@CBVB+N|B= z2+_KS%z|A44}hwcmfoqhvw9`Z4fEnKrGn@|%n&iA5JYDM^cbAhjxl{%jaK9}VxClN zk`Ck^Bg*Y@SCdN0gD@Br1RZ(Kl?pMPj{?!?u4N+G%<>(T=`?`J#|?SVjx zFOfz)I*kmTMi6Po32A4s#&@}zR z2n+totAsmc_#W2bF8@3eM}C~5FPO=Wb(-(cAjzk*K%h z0=GkhqFw_vYA#2&1_aCjd!UQ)*KJ(kmZNu?EAKX=zaITO&Oz-*3Ca;b6S!5F7{|y> zz_$PfP0lQI)Z&#jj<4bff%bRLVhPnvM=Pcfy2dn(BP`}@5u-%s#6cH5g7*%CXVO8} z4$|ZhiSFidbnV3(E-dug%@@h?d-FwvyD}{FRmErx@y)p3)Zp%@o>67>P^W0R$mJu# zblxS@kvO;LK07Ii-cgXSl#R}9Qp>bjI24)eP97XSLLe5n*;&gdC9;`$u>if zRv!69g3_mY?-~n3XmL8;xOUwQ;l_=18^U@R!umt%9TRJAbAQ`!-@B(kh1LaaI*mtC zr$rvKa_+0QeRXD!7viW$R?28?o#@)@*H?Xb?RDRatLjAvN5lNd5YgOS2ucuTr$f{1 z%k@gF+VZ%noI}!nx*#MsU;)7jKI^7C&DmjS;~AG`1$1dPPY?RRuL6gTD$CnB;HxIj zq7fIvb4>IT;!*Ov)*jm4_NNv{h`gE%Ug!{rFU(wgP7e*@L*xE350DL@REFS|#q5XZBG`odU z&Bb}?BPwpaCo@eSe`eBS!Yb-4T!(;MhIh#7;A5-_diD+3o;=4IKV6$?yyRT)^#yFI z7Re_}#Onb$g_3d1)5Z<-5l0Di)T*MohMyb^;~6|xXhE5Z=NxnOB>Umb{7~y=;RvF4 z*IiHiW_l7!_USlmM@(Z!mwRn(#1{$?UeMCIutk%GV;rdRM}ZMcuVAopdWc7+UrkGK zxqUS!khcq(_S+5-a~@Zf+pAY&Tsd6tukFR*E@ofs$9`G3=t`Iod0{*@?kN+AmMFl| zh=jY^wDm4gHBAk z;Qg}FwdvdX0D2-Be58ChS;>br$=0P-P&_)$CB>y6BFT&U*w_&HV+SK3DLnKIT4jJe zd@9ZJ(%UTL3V|AWjC?|D8B&+o*#rMw{t`jyukPUnO5QQ5*G5?F> z$7qQ~c4Q(QofyM>K&k=gUw?QhrebtpN8mSJ`Ylk_Y1r)_04?DU1u4F$Z&jyNyT49r zP7PPYJ;ruG+9)X~06jz$ZZ0B;bNZ9Q zO{r}T3F33x3LY;mB)4C;jUo{MW7j2*j|C_Ob)y&o?a(JLwvO{I2Zw2=fsa|#idHB* zTG}$A24|X75?^{=^@rGy#mFfe=G|mV3xp<(g?=fi6`5wtug&k>S~qduIC3VNNx9gUF2+Or7HM_GOdvnrw&yk!1@ISZ6ohe~WNB+V&C zxLL(P$${Un&p>Vh&Y)xvrbmsbaCB*V+k9a+7rEi)(4us$bHo> z6kEo<-BD&{Tc>tzZ`+M6eEq<3>L0b~9}aRN$u!aS7|JaNg{J{z0Lw+F)p4_+`KbD|j7;*WjngUA=Eak6J3*K$*a*LVTtX=Y- zj2AbVspT=rVy*Q>+Q9HPZfYvqJ5ecdJl=n@dOSE1R5mT={zLOu#kj*6pcgOqIe?AP z?W@Us?BQe-7`A#nc++@*(Vb%^_5D?kd4@m8=I!LCAb{wAk&J?~ z;1UCc8r&A+Qe>0SQnMUxq3dqE!Lh}jbjgFq#RccRU;6H^uF3tzYRbb1E-W_ZK1`+T ze&m)XaJWhiAd8~&TF{AXBoyF0FQGbVmk2%57Mb&~PTeAS5-W42Go!euO%VZ7Z)zzz z;nl4dSJ~ow$OGirguBZ^mVwFbD8nD_XFeSO;&OMCkjfK6RUXvRuH&oSU{?at3+TRj za0k7BgWxUMeN2?$@If$gzSEpSRYx>t9413AG$n}Jmg7NuSB5nDRLMDzu!;N}n7X6B zTMA3_*b!Ak>tPzlg?@}u$4J^AgDv2RWyU73JLLv&&&}Eb7N3VW%~38Zn1#rlxC7>q zCeKi$$v^qb;=6z4=re;q`1Hbpe)xA`%_8QmNEO@aU3lv1uc*W+nJ%lMRi)p)PP0>P zM~$y2Lp8p-L^%L!fG!wLEc-ZY;2hY>(}nSW2p;=mDIL85Zrr|Rnob*GB<~+L+oW^z zPt{AH&GO?9uY`QC{^jQJuaR!cKWw8Ce$4;5kT+WF2IWR^IncMI8i=-iIV%u_7tuls9S1;?2;RN z!c1)WnUF@qGL;H}<%p>OQf>v68Ppi7EwiB60}U>=_}jl;n>szA=tygbDI9-;wQwIC zo7~zv=Qtam#HRa>a)Rkx-c6G(pr3$tkO!zEA&HFSe>52-8S~H&uKY6inO`0!5-;GG zhUy`l$Kg7z7!J6*UUqtvX+#b}g_q;H1H4yvkih*XIZ#4xm#h|`# zBx1Du$-&HFjrScE56HBm)=*hWaF7wLD1$0m%!*RKSDD6YD+Dm$-S!xq|J>k`=bR}C z^1w$mM>erCiFR|!_w1j_@ zkgEsIQ;=*)y`iEY#5Ngq5@WAJu^I?`@wl~qv=le61soz#YIjI{w(_;?ZD zsk(4q1&(0aIIX^MivVcfp%3)_b$gw@INv1+nWSzNMB!4mzB+jI4Pnuni%%Ak;+mI) zI@bK#omWQ_%_Ar2V=mmq{wVE3vtwP1LN}r{k?6^JP*h{Z;O#xpLC|#?1lPC{$R`D> zAn2*8u9k|E5&UgKzfx&;D(Yp;)oarjXVq(CDCcWv{e*^)AxoXG7msVM~ z^khlM8?!2lGl-MsWxbpiFHFwpoe&~RT2@W=hdM9*1e@{YAN{*ayTs${}d2-=at|7Klxu? zSkTYsCqQ7IY3CyVgWslfnOtU0L48H&XLxK54xbokF$agmL!(8^sDpMw#Pb78VD6vA zIXE1$a&g#@m6ei;Bm|vG2U~<9HB>~rf6$- zi{uiLs>$fzNh*SJDLOdzi6#Chq^UVVpKPoHHG7{<3Yd@;-EjYvIYK{2=mnPxh<6n} z%+PQ$Bsxdv=Lo%g#AJinaE6#ULQk=l4*{V^ClCT^*y|jrf7X%uKcO(dl)JXYzs=pX zJ+MY459NOQ=mg8TySA`%uOkK|bX*g75+A}jB0opupO&ru@j>M0?%H*$Ux3AHV6Z$T zP5^^->v3SPo;e>FY;VlLU_wym?%E%c0kQ6PdO`MtIT$<#gTrY(2ZK2QsC{;V%lF(} zd$3o`-L=o*_dR#lp1W%+QxZPdb9e2N9NFjY+N5mnD(b%Gn6}sL%7bGAX5L~NyntZB%yr*L(cl9UUzx$;h`_n&-WA?(r4*-Fu+_h^UaNT5ZvE^akp1W&5 z#)oqb4$r~i5*%lzpV|3$p3B2{|EXq$r^~?z0B^zrxODx!WVt-%?%GMQF?ZK~n!NXB zAoTK{J^`WEt;Zqsdggoxy}dC<=;sLihZo{AN9ZLUF4mpRy~@Yc4>6RAjeefwD0y$5 zP~{``s+qfM&)v1>303wi(L*uXu@n<~#e;c56*B=~RjVaJ3EX`~YK|YhTk=Ub^cD{<6`I%HzIX5b!8NKoP_~7P04G&?Zt&QWH9!o~fc&Sen`Oc0HHsL~czQr0NY`^q z`t6G>S?MDxOBF9!63=>_0uBr%4Gm=stabrumq-Ds%eRKuv4#{cw(C^#p;pmeP6--U zxhKp5V~lKsNkK6pUoN4k%muHfgMvcIuJWi)Zmou^3cnb^Mk;)IiPIeVafKw{${SAV z?l%7|6MEOuf1jn!uPtxwou@Ko_?ws#t!|;}G8?%V&53q|K&kS)mKv&so+u9At?E|| zk{&s2t;mo{-Z{E&TN!EcdhPKbX`Wq$--*&6ojfVB=ZDIjPPrIWl#8G~qqTE5^F3V{ z`E)&I&FGeA6(lIb^ji!M6h~FoV+K^`T61rZkNP{uJfmt)O`A#{i&jp%OAgLB^RMvm zIp*z2LKZ?P+YVXRd;uTF5<__v@(rjyB1NZ$#()a4Wg7|gtmug)J-*xvtzs|~h;d6d zo=$$9gqW??HBy0d8u=9S3QS&EU$vrszKHy+S6g~+?IsCDXTKW07)z**hQCE-xb5!wxTsqWG z)WK?%apGv7#D6Pfes`V`s-2d^LrSJm>^p8JzO*t&6e$7yiGG)}j1)YsIF;jS)$1w; z*a)TTs9Q-%J?dwwK#;PgkdVG3ZbbDBo#Iw!PD!=9R{>V0?^EHSvjwyyHwFrX}&|w*+Zr8rmK&Vq;HLtrOG1>Lf$3fusrPx z$Wn6X3n8g^_^hQFjVV;bK;YJr27;KPA-Kvxw=!R=LSL+viYrLrts?4MjL^JU3JviN ziH)tV%4BSZ{iwH-kY1cmOms{u=BuzH_=HnqCqYSK=Tk4SjcCSd(n?Q5joX$c!4f`G zKj{lrX;OYRk_X~w7q&i(8mlD5 zIEr50h+4`ZOmw+xxver))pV&l)odYe?NK!?6J_a{KrErtcJd-H4AP_L3yt}aDn@BL z?8X!Zv}ti==~tUKu8q|+%v3^)J>X;dU+Dk*eEsnCYM8W^1aGfZB|(kP;M z)889s7u?|R+RQjXP!hI`d|EGgvz>NWyW&HS1Ns zVGUlfZ*&6_;C_Uvw-^JyFmss&~4_NR&aC-(Z~${!^$H zO`}D^;&ZBTmlQ z+02p5W+Wl^0|n}x%8tIS?C9 z>o2A%4O9E<+UakvjlZo~vh31~qLKC(jvLBylv9D3A}SIGNINdCSYb$LP(F$Ek*G#% zN5NcKaN{EyDl=B|l)s2vwMBJ-Zh{VWar)OSq}8r^bF72XiV(2pITvMp#^!V-ZyLsBumk*#Kd8~>?vE4_&VD-gt;dihTgg@5}r;;r7k4QUNf zQ0_X_%BadWa2owq4ZMka_I`=j+4~S+H=(RlXR5#}{sa%r$$HAkVDS+Jl;*DTc)5RU zT;+e2%;u8-{Q4Y#p9AnkOkptO=s5sCjAj%9odfXFyypP?9Dp~$b`HQ}2%UwY&jI*3 z0B9i(woc2HLE`|{E>c(DFDugSKL7y1ZO zmOM1@dy<#j&M=lmJ;!?9N%SG}^kfR~+)_l$PLA!0FS^I?-oDiS;7@gWn~MOopDy?6 zCrKJ}L7qc*z4K{(&uN#Pxxe!k zyJV<+c}K!K#ZTSL;p-|!iMz3`Rk&hNa@(O}ug%R2pXlJcd~LPj(qf?oy{yRbQR5#u z-xZuCakj+~x=44r384a1MUTi2CubIJ1h?Y$exBozDi3No@47x%Rlu^}#78yvnz)Lu zuJW-F8ENP5i6qqbG)1_v6)GJ24M)4@^XR&7(P5jSr93=T+ zH;WspBg?=20g?rNRqws@evm(KVy^cdN8q^I-;Xs5gw!+Fyy9!ThV=86;W51A1UekB_4ei7`|OeGaDY77sMX)lSJ*Wt>_lce-Exjv z?q&D;SzZhq%f2s`ll3_*_p-)VmzJ7PK?jl!9ntn@d!0izcoIzY-l)gx9Mg6Z5u8Oh zLE-GdXz`Hd#X8i_~GSy zD02I0yHPq3qv=^)$rGwgz9{}xy5Yi27F5R-Z%+x`U;#m~}s z>zs|>yisr6G0}Rz$D6D5E*$2NYVgFpd%NsrycS$wyLo#xo*Q#YEHp3h&uT@a}4*8O*`EIe15OC;S6(j_|a>Q0jW0 zxWoc^uDmwGy|kV$!oOJ;9hU0{IH#a8e$2tUF0$BIZiF};Y!Nz!PHq9IF#ZRdi(K4X zyr`J;0Oz0zm4La9nR5u}bMVf^!p*@uZv&cxcV|xuKL_s&W1NKe=HT5NypvUz^dEEZ zZVujoEziNbaYzQ#Zw}tg!8@T@{38r_4&Is3abphNy%fFYbMDWp#3@mL{<%a|;y&@G ze)`8;ow%R+XL;UwU-}(buA8RY2k*(o?{0q}f{RO)9fB);>Tm2!~ zaP>o8-149AeCxk8W;8Y3$oKFoy!A`J{tJKmCl?m<^VzZq_2`*RBm`gv}wNB??zqAg28;M0w) z%xv(#(BSbu{QjSB{>6m_{rt7D2LF9`qQS#NxPH=^-PrHY*mu9Zvisd1Us%x3=f)aa ze}AGeIW6>Jf4=jrAO6Y0!diQf7lqp}6!@AR*gvI-$>NNmjqiNxfBvuWCIUsxZ0b*# zqmz5{D?jy#g|%8Li&$Km-7*<()N?gL`TA`WY-Z;j_%-#Lv1Qi3y*B0o1c~@&Et; delta 316 zcmWm8u}i~X5Cw3#RHaEs@%z%a1VKR%!AaDq#SU!>Di|G|92AQk911SFYDe`PluRLl zE)H^Acd?3t#X+Q#2JB>$wEh8JZ@Ks0^b0#m*irQL6(wa34bBh``feKLEvltAG|&V! z&2t@ElC~x@qwp@NI&)&co20+RU<5uTUFTj@c#^bJ4pJ~KNL@R!;aO6>@}ojV&_Qj& zgbNE@Us6y99j+|8y}Z+*Ey?w>I($jmyXStICDCXb;8_$k`MZ){KPnpd7A=MQ8W;j^ zx(?Iuw_9QvE=0==!^5b+E^I^&vrvgwS%q#~W)ymHkxlp=7np?8Ns&cZO`HKJ*R!#e oOtA-#Iga&anZ#?JMZA<4gpb@j3y1r8=HOZ1#_j$bi#Trn1KLh@t^fc4 diff --git a/tests/models/cross/test_cca.py b/tests/models/cross/test_cca.py new file mode 100644 index 0000000..b603e9c --- /dev/null +++ b/tests/models/cross/test_cca.py @@ -0,0 +1,292 @@ +import dask.array as da +import numpy as np +import pytest +import xarray as xr + +from xeofs.cross import CCA + + +def generate_random_data(shape, lazy=False, seed=142): + rng = np.random.default_rng(seed) + if lazy: + return xr.DataArray( + da.random.random(shape, chunks=(5, 5)), + dims=["sample", "feature"], + coords={"sample": np.arange(shape[0]), "feature": np.arange(shape[1])}, + ) + else: + return xr.DataArray( + rng.random(shape), + dims=["sample", "feature"], + coords={"sample": np.arange(shape[0]), "feature": np.arange(shape[1])}, + ) + + +def generate_well_conditioned_data(lazy=False): + rng = np.random.default_rng(142) + t = np.linspace(0, 50, 200) + std = 0.1 + x1 = np.sin(t)[:, None] + rng.normal(0, std, size=(200, 2)) + x2 = np.sin(t)[:, None] + rng.normal(0, std, size=(200, 3)) + x1[:, 1] = x1[:, 1] ** 2 + x2[:, 1] = x2[:, 1] ** 3 + x2[:, 2] = abs(x2[:, 2]) ** (0.5) + coords_time = np.arange(len(t)) + coords_fx = [1, 2] + coords_fy = [1, 2, 3] + X = xr.DataArray( + x1, + dims=["sample", "feature"], + coords={"sample": coords_time, "feature": coords_fx}, + ) + Y = xr.DataArray( + x2, + dims=["sample", "feature"], + coords={"sample": coords_time, "feature": coords_fy}, + ) + if lazy: + X = X.chunk({"sample": 5, "feature": -1}) + Y = Y.chunk({"sample": 5, "feature": -1}) + return X, Y + else: + return X, Y + + +@pytest.fixture +def cca(): + return CCA(n_modes=1) + + +def test_initialization(): + model = CCA() + assert model is not None + + +def test_fit(cca): + X, Y = generate_well_conditioned_data() + cca.fit(X, Y, dim="sample") + assert hasattr(cca, "preprocessor1") + assert hasattr(cca, "preprocessor2") + assert hasattr(cca, "data") + + +def test_fit_empty_data(cca): + with pytest.raises(ValueError): + cca.fit(xr.DataArray(), xr.DataArray(), "time") + + +def test_fit_invalid_dims(cca): + X, Y = generate_well_conditioned_data() + with pytest.raises(ValueError): + cca.fit(X, Y, dim=("invalid_dim1", "invalid_dim2")) + + +def test_transform(cca): + X, Y = generate_well_conditioned_data() + cca.fit(X, Y, dim="sample") + result = cca.transform(X, Y) + assert isinstance(result, list) + assert isinstance(result[0], xr.DataArray) + + +def test_transform_unseen_data(cca): + X, Y = generate_well_conditioned_data() + x = X.isel(sample=slice(151, 200)) + y = Y.isel(sample=slice(151, 200)) + X = X.isel(sample=slice(None, 150)) + Y = Y.isel(sample=slice(None, 150)) + + cca.fit(X, Y, "sample") + result = cca.transform(x, y) + assert isinstance(result, list) + assert isinstance(result[0], xr.DataArray) + # Check that unseen data can be transformed + assert result[0].notnull().all() + assert result[1].notnull().all() + + +def test_inverse_transform(cca): + X, Y = generate_well_conditioned_data() + cca.fit(X, Y, "sample") + # Assuming mode as 1 for simplicity + scores1 = cca.data["scores1"].sel(mode=1) + scores2 = cca.data["scores2"].sel(mode=1) + Xrec1, Xrec2 = cca.inverse_transform(scores1, scores2) + assert isinstance(Xrec1, xr.DataArray) + assert isinstance(Xrec2, xr.DataArray) + + +@pytest.mark.parametrize("use_pca", [False, True]) +def test_squared_covariance_fraction(use_pca): + X, Y = generate_well_conditioned_data() + cca = CCA(n_modes=2, use_pca=use_pca, n_pca_modes="all") + cca.fit(X, Y, "sample") + scf = cca.squared_covariance_fraction() + assert isinstance(scf, xr.DataArray) + assert all(scf <= 1), "Squared covariance fraction is greater than 1" + + +@pytest.mark.parametrize("use_pca", [False, True]) +def test_total_squared_covariance(use_pca): + X, Y = generate_well_conditioned_data() + + # Compute total squared covariance + X_ = X.rename({"feature": "x"}) + Y_ = Y.rename({"feature": "y"}) + cov_mat = xr.cov(X_, Y_, dim="sample") + tsc = (cov_mat**2).sum() + + cca = CCA(n_modes=2, use_pca=use_pca, n_pca_modes="all") + cca.fit(X, Y, "sample") + tsc_model = cca.data["total_squared_covariance"] + xr.testing.assert_allclose(tsc, tsc_model) + + +def test_fit_different_coordinates(): + """Like a lagged CCA scenario""" + X, Y = generate_well_conditioned_data() + X = X.isel(sample=slice(0, 99)) + Y = Y.isel(sample=slice(100, 199)) + cca = CCA(n_modes=2, use_pca=False) + cca.fit(X, Y, "sample") + r = cca.cross_correlation_coefficients() + # Correlation coefficents are not zero + assert np.all(r > np.finfo(r.dtype).eps) + + +@pytest.mark.parametrize( + "dim", + [(("time",)), (("lat", "lon")), (("lon", "lat"))], +) +def test_components(mock_data_array, dim): + cca = CCA(n_modes=2, use_pca=False) + cca.fit(mock_data_array, mock_data_array, dim) + components1, components2 = cca.components() + feature_dims = tuple(set(mock_data_array.dims) - set(dim)) + assert isinstance(components1, xr.DataArray) + assert isinstance(components2, xr.DataArray) + assert set(components1.dims) == set( + ("mode",) + feature_dims + ), "Components1 does not have the right feature dimensions" + assert set(components2.dims) == set( + ("mode",) + feature_dims + ), "Components2 does not have the right feature dimensions" + + +@pytest.mark.parametrize("shapeX", [(30, 10)]) +@pytest.mark.parametrize("shapeY", [(30, 10), (30, 5), (30, 15)]) +@pytest.mark.parametrize("use_pca", [False, True]) +def test_components_coordinates(shapeX, shapeY, use_pca): + # Test that the components have the right coordinates + X = generate_random_data(shapeX) + Y = generate_random_data(shapeY) + + cca = CCA(n_modes=2, use_pca=use_pca, n_pca_modes="all") + cca.fit(X, Y, "sample") + components1, components2 = cca.components() + xr.testing.assert_equal(components1.coords["feature"], X.coords["feature"]) + xr.testing.assert_equal(components2.coords["feature"], Y.coords["feature"]) + + +@pytest.mark.parametrize("correction", [(None), ("fdr_bh")]) +def test_homogeneous_patterns(correction): + X = generate_random_data((200, 10), seed=123) + Y = generate_random_data((200, 20), seed=321) + + cca = CCA(n_modes=10, use_pca=False) + cca.fit(X, Y, "sample") + + _ = cca.homogeneous_patterns(correction=correction) + + +@pytest.mark.parametrize( + "correction", + [(None), ("fdr_bh")], +) +def test_heterogeneous_patterns(correction): + X = generate_random_data((200, 10), seed=123) + Y = generate_random_data((200, 20), seed=321) + + cca = CCA(n_modes=10, use_pca=False) + cca.fit(X, Y, "sample") + + _ = cca.heterogeneous_patterns(correction=correction) + + +def test_predict(): + X = generate_random_data((200, 10), seed=123) + Y = generate_random_data((200, 20), seed=321) + + cca = CCA(n_modes=10, use_pca=False) + cca.fit(X, Y, "sample") + + Xnew = generate_random_data((200, 10), seed=123) + + Ry_pred = cca.predict(Xnew) + _ = cca.inverse_transform(Y=Ry_pred) + + +@pytest.mark.parametrize("engine", ["netcdf4", "zarr"]) +def test_save_load(tmp_path, engine): + """Test save/load methods in MCA class, ensuring that we can + roundtrip the model and get the same results when transforming + data.""" + X = generate_random_data((200, 10), seed=123) + Y = generate_random_data((200, 20), seed=321) + + original = CCA() + original.fit(X, Y, "sample") + + # Save the CCA model + original.save(tmp_path / "cca", engine=engine) + + # Check that the CCA model has been saved + assert (tmp_path / "cca").exists() + + # Recreate the model from saved file + loaded = CCA.load(tmp_path / "cca", engine=engine) + + # Check that the params and DataContainer objects match + assert original.get_params() == loaded.get_params() + assert all([key in loaded.data for key in original.data]) + for key in original.data: + if original.data._allow_compute[key]: + assert loaded.data[key].equals(original.data[key]) + else: + # but ensure that input data is not saved by default + assert loaded.data[key].size <= 1 + assert loaded.data[key].attrs["placeholder"] is True + + # Test that the recreated model can be used to transform new data + assert np.allclose( + original.transform(X, Y), + loaded.transform(X, Y), + ) + + # The loaded model should also be able to inverse_transform new data + XYr_o = original.inverse_transform(*original.scores()) + XYr_l = loaded.inverse_transform(*loaded.scores()) + assert np.allclose(XYr_o[0], XYr_l[0]) + assert np.allclose(XYr_o[1], XYr_l[1]) + + +def test_serialize_deserialize_dataarray(mock_data_array): + """Test roundtrip serialization when the model is fit on a DataArray.""" + model = CCA() + model.fit(mock_data_array, mock_data_array, "time") + dt = model.serialize() + rebuilt_model = CCA.deserialize(dt) + assert np.allclose( + model.transform(mock_data_array), rebuilt_model.transform(mock_data_array) + ) + + +def test_serialize_deserialize_dataset(mock_dataset): + """Test roundtrip serialization when the model is fit on a Dataset.""" + model = CCA() + model.fit(mock_dataset, mock_dataset, "time") + dt = model.serialize() + rebuilt_model = CCA.deserialize(dt) + assert np.allclose( + model.transform(mock_dataset), rebuilt_model.transform(mock_dataset) + ) diff --git a/tests/models/cross/test_cpcca.py b/tests/models/cross/test_cpcca.py index 20035d2..258ff14 100644 --- a/tests/models/cross/test_cpcca.py +++ b/tests/models/cross/test_cpcca.py @@ -272,3 +272,70 @@ def test_predict(): Ry_pred = cpcca.predict(Xnew) _ = cpcca.inverse_transform(Y=Ry_pred) + + +@pytest.mark.parametrize("engine", ["netcdf4", "zarr"]) +@pytest.mark.parametrize("alpha", [0.0, 0.5, 1.0]) +def test_save_load(tmp_path, engine, alpha): + """Test save/load methods in MCA class, ensuring that we can + roundtrip the model and get the same results when transforming + data.""" + X = generate_random_data((200, 10), seed=123) + Y = generate_random_data((200, 20), seed=321) + + original = CPCCA(alpha=alpha) + original.fit(X, Y, "sample") + + # Save the CPCCA model + original.save(tmp_path / "cpcca", engine=engine) + + # Check that the CPCCA model has been saved + assert (tmp_path / "cpcca").exists() + + # Recreate the model from saved file + loaded = CPCCA.load(tmp_path / "cpcca", engine=engine) + + # Check that the params and DataContainer objects match + assert original.get_params() == loaded.get_params() + assert all([key in loaded.data for key in original.data]) + for key in original.data: + if original.data._allow_compute[key]: + assert loaded.data[key].equals(original.data[key]) + else: + # but ensure that input data is not saved by default + assert loaded.data[key].size <= 1 + assert loaded.data[key].attrs["placeholder"] is True + + # Test that the recreated model can be used to transform new data + assert np.allclose( + original.transform(X, Y), + loaded.transform(X, Y), + ) + + # The loaded model should also be able to inverse_transform new data + XYr_o = original.inverse_transform(*original.scores()) + XYr_l = loaded.inverse_transform(*loaded.scores()) + assert np.allclose(XYr_o[0], XYr_l[0]) + assert np.allclose(XYr_o[1], XYr_l[1]) + + +def test_serialize_deserialize_dataarray(mock_data_array): + """Test roundtrip serialization when the model is fit on a DataArray.""" + model = CPCCA() + model.fit(mock_data_array, mock_data_array, "time") + dt = model.serialize() + rebuilt_model = CPCCA.deserialize(dt) + assert np.allclose( + model.transform(mock_data_array), rebuilt_model.transform(mock_data_array) + ) + + +def test_serialize_deserialize_dataset(mock_dataset): + """Test roundtrip serialization when the model is fit on a Dataset.""" + model = CPCCA() + model.fit(mock_dataset, mock_dataset, "time") + dt = model.serialize() + rebuilt_model = CPCCA.deserialize(dt) + assert np.allclose( + model.transform(mock_dataset), rebuilt_model.transform(mock_dataset) + ) diff --git a/xeofs/cross/__init__.py b/xeofs/cross/__init__.py index 830df25..456927c 100644 --- a/xeofs/cross/__init__.py +++ b/xeofs/cross/__init__.py @@ -1,17 +1,21 @@ import warnings +from .cca import CCA, ComplexCCA, HilbertCCA from .cpcca import CPCCA, ComplexCPCCA, HilbertCPCCA from .cpcca_rotator import ComplexCPCCARotator, CPCCARotator, HilbertCPCCARotator from .mca import MCA, ComplexMCA, HilbertMCA from .mca_rotator import ComplexMCARotator, HilbertMCARotator, MCARotator __all__ = [ + "CCA", + "MCA", + "CPCCA", + "ComplexCCA", "ComplexMCA", "ComplexCPCCA", + "HilbertCCA", "HilbertMCA", "HilbertCPCCA", - "MCA", - "CPCCA", "MCARotator", "CPCCARotator", "ComplexMCARotator", diff --git a/xeofs/cross/cca.py b/xeofs/cross/cca.py new file mode 100644 index 0000000..524777f --- /dev/null +++ b/xeofs/cross/cca.py @@ -0,0 +1,356 @@ +from typing import Sequence + +import numpy as np + +from .cpcca import CPCCA, ComplexCPCCA, HilbertCPCCA + + +class CCA(CPCCA): + """Canonical Correlation Analysis (CCA). + + CCA seeks to find paris of coupled patterns that maximize the correlation [1]_ [2]_ . + + This method solves the following optimization problem: + + :math:`\\max_{q_x, q_y} \\left( q_x^T X^T Y q_y \\right)` + + subject to the constraints: + + :math:`q_x^T (X^TX) q_x = 1, \\quad q_y^T (Y^TY) q_y = 1` + + where :math:`X` and :math:`Y` are the input data matrices and :math:`q_x` + and :math:`q_y` are the corresponding pattern vectors. + + Parameters + ---------- + n_modes : int, default=2 + Number of modes to calculate. + standardize : Squence[bool] | bool, default=False + Whether to standardize the input data. Generally not recommended as + standardization can be managed by the degree of whitening. + use_coslat : Sequence[bool] | bool, default=False + For data on a longitude-latitude grid, whether to correct for varying + grid cell areas towards the poles by scaling each grid point with the + square root of the cosine of its latitude. + use_pca : Sequence[bool] | bool, default=False + Whether to preprocess each field individually by reducing dimensionality + through PCA. The cross-covariance matrix is then computed in the reduced + principal component space. + n_pca_modes : Sequence[int | float | str] | int | float | str, default=0.999 + Number of modes to retain during PCA preprocessing step. If int, + specifies the exact number of modes; if float, specifies the fraction of + variance to retain; if "all", all modes are retained. + pca_init_rank_reduction : Sequence[float] | float, default=0.3 + Relevant when `use_pca=True` and `n_pca_modes` is a float. Specifies the + initial fraction of rank reduction for faster PCA computation via + randomized SVD. + check_nans : Sequence[bool] | bool, default=True + Whether to check for NaNs in the input data. Set to False for lazy model + evaluation. + compute : bool, default=True + Whether to compute the model elements eagerly. If True, the following + are computed sequentially: preprocessor scaler, optional NaN checks, SVD + decomposition, scores, and components. + random_state : numpy.random.Generator | int | None, default=None + Seed for the random number generator. + sample_name : str, default="sample" + Name for the new sample dimension. + feature_name : Sequence[str] | str, default="feature" + Name for the new feature dimension. + solver : {"auto", "full", "randomized"} + Solver to use for the SVD computation. + solver_kwargs : dict, default={} + Additional keyword arguments passed to the SVD solver function. + + + References + ---------- + .. [1] Bretherton, C., Smith, C., Wallace, J., 1992. An intercomparison of + methods for finding coupled patterns in climate data. Journal of climate + 5, 541–560. + .. [2] Wilks, D. S. Statistical Methods in the Atmospheric Sciences. + (Academic Press, 2019). + doi:https://doi.org/10.1016/B978-0-12-815823-4.00011-0. + + Examples + -------- + + Perform CCA on two datasets on a regular longitude-latitude grid: + + >>> model = CCA(n_modes=5, use_coslat=True) + >>> model.fit(X, Y, dim="time") + + """ + + def __init__( + self, + n_modes: int = 2, + standardize: Sequence[bool] | bool = False, + use_coslat: Sequence[bool] | bool = False, + check_nans: Sequence[bool] | bool = True, + use_pca: Sequence[bool] | bool = True, + n_pca_modes: Sequence[float | int | str] | float | int | str = 0.999, + pca_init_rank_reduction: Sequence[float] | float = 0.3, + compute: bool = True, + sample_name: str = "sample", + feature_name: Sequence[str] | str = "feature", + solver: str = "auto", + random_state: np.random.Generator | int | None = None, + solver_kwargs: dict = {}, + ): + CPCCA.__init__( + self, + n_modes=n_modes, + alpha=[0.0, 0.0], + standardize=standardize, + use_coslat=use_coslat, + check_nans=check_nans, + use_pca=use_pca, + n_pca_modes=n_pca_modes, + pca_init_rank_reduction=pca_init_rank_reduction, + compute=compute, + sample_name=sample_name, + feature_name=feature_name, + solver=solver, + random_state=random_state, + solver_kwargs=solver_kwargs, + ) + self.attrs.update({"model": "Canonical Correlation Analysis"}) + # Renove alpha from the inherited CPCCA serialization params because it is hard-coded for CCA + self._params.pop("alpha") + + +class ComplexCCA(ComplexCPCCA, CCA): + """Complex CCA. + + CCA applied to a complex-valued field obtained from a pair of variables such + as the zonal and meridional components, :math:`U` and :math:`V`, of the wind + field. Complex CCA analysis then maximizes the correlation between + two datasets of the form + + .. math:: + Z_x = U_x + iV_x + + and + + .. math:: + Z_y = U_y + iV_y + + into a set of complex-valued components and PC scores. + + + Parameters + ---------- + n_modes : int, default=2 + Number of modes to calculate. + standardize : Squence[bool] | bool, default=False + Whether to standardize the input data. Generally not recommended as + standardization can be managed by the degree of whitening. + use_coslat : Sequence[bool] | bool, default=False + For data on a longitude-latitude grid, whether to correct for varying + grid cell areas towards the poles by scaling each grid point with the + square root of the cosine of its latitude. + use_pca : Sequence[bool] | bool, default=False + Whether to preprocess each field individually by reducing dimensionality + through PCA. The cross-covariance matrix is then computed in the reduced + principal component space. + n_pca_modes : Sequence[int | float | str] | int | float | str, default=0.999 + Number of modes to retain during PCA preprocessing step. If int, + specifies the exact number of modes; if float, specifies the fraction of + variance to retain; if "all", all modes are retained. + pca_init_rank_reduction : Sequence[float] | float, default=0.3 + Relevant when `use_pca=True` and `n_pca_modes` is a float. Specifies the + initial fraction of rank reduction for faster PCA computation via + randomized SVD. + check_nans : Sequence[bool] | bool, default=True + Whether to check for NaNs in the input data. Set to False for lazy model + evaluation. + compute : bool, default=True + Whether to compute the model elements eagerly. If True, the following + are computed sequentially: preprocessor scaler, optional NaN checks, SVD + decomposition, scores, and components. + random_state : numpy.random.Generator | int | None, default=None + Seed for the random number generator. + sample_name : str, default="sample" + Name for the new sample dimension. + feature_name : Sequence[str] | str, default="feature" + Name for the new feature dimension. + solver : {"auto", "full", "randomized"} + Solver to use for the SVD computation. + solver_kwargs : dict, default={} + Additional keyword arguments passed to the SVD solver function. + + Examples + -------- + + With two DataArrays `u_i` and `v_i` representing the zonal and meridional + components of the wind field for two different regions :math:`x` and + :math:`y`, construct + + >>> X = u_x + 1j * v_x + >>> Y = u_y + 1j * v_y + + and fit the Complex CCA model: + + >>> model = ComplexCCA(n_modes=5) + >>> model.fit(X, Y, "time") + + + """ + + def __init__( + self, + n_modes: int = 2, + standardize: Sequence[bool] | bool = False, + use_coslat: Sequence[bool] | bool = False, + check_nans: Sequence[bool] | bool = True, + use_pca: Sequence[bool] | bool = True, + n_pca_modes: Sequence[float | int | str] | float | int | str = 0.999, + pca_init_rank_reduction: Sequence[float] | float = 0.3, + compute: bool = True, + sample_name: str = "sample", + feature_name: Sequence[str] | str = "feature", + solver: str = "auto", + random_state: np.random.Generator | int | None = None, + solver_kwargs: dict = {}, + ): + ComplexCPCCA.__init__( + self, + n_modes=n_modes, + alpha=[0.0, 0.0], + standardize=standardize, + use_coslat=use_coslat, + check_nans=check_nans, + use_pca=use_pca, + n_pca_modes=n_pca_modes, + pca_init_rank_reduction=pca_init_rank_reduction, + compute=compute, + sample_name=sample_name, + feature_name=feature_name, + solver=solver, + random_state=random_state, + solver_kwargs=solver_kwargs, + ) + self.attrs.update({"model": "Complex CCA"}) + # Renove alpha from the inherited CPCCA serialization params because it is hard-coded for CCA + self._params.pop("alpha") + + +class HilbertCCA(HilbertCPCCA, ComplexCCA): + """Hilbert CCA. + + Hilbert CCA extends CCA by examining amplitude-phase relationships. It + augments the input data with its Hilbert transform, creating a + complex-valued field. + + This method solves the following optimization problem: + + :math:`\\max_{q_x, q_y} \\left( q_x^T X^T Y q_y \\right)` + + subject to the constraints: + + :math:`q_x^T (X^TX) q_x = 1, \\quad q_y^T (Y^TY) q_y = 1` + + where :math:`H` denotes the conjugate transpose and :math:`X` and :math:`Y` + are the augmented data matrices. + + An optional padding with exponentially decaying values can be applied prior + to the Hilbert transform in order to mitigate the impact of spectral + leakage. + + + Parameters + ---------- + n_modes : int, default=2 + Number of modes to calculate. + padding : Sequence[str] | str | None, default="exp" + Padding method for the Hilbert transform. Available options are: - None: + no padding - "exp": exponential decay + decay_factor : Sequence[float] | float, default=0.2 + Decay factor for the exponential padding. + standardize : Squence[bool] | bool, default=False + Whether to standardize the input data. Generally not recommended as + standardization can be managed by the degree of whitening. + use_coslat : Sequence[bool] | bool, default=False + For data on a longitude-latitude grid, whether to correct for varying + grid cell areas towards the poles by scaling each grid point with the + square root of the cosine of its latitude. + use_pca : Sequence[bool] | bool, default=False + Whether to preprocess each field individually by reducing dimensionality + through PCA. The cross-covariance matrix is computed in the reduced + principal component space. + n_pca_modes : Sequence[int | float | str] | int | float | str, default=0.999 + Number of modes to retain during PCA preprocessing step. If int, + specifies the exact number of modes; if float, specifies the fraction of + variance to retain; if "all", all modes are retained. + pca_init_rank_reduction : Sequence[float] | float, default=0.3 + Relevant when `use_pca=True` and `n_pca_modes` is a float. Specifies the + initial fraction of rank reduction for faster PCA computation via + randomized SVD. + check_nans : Sequence[bool] | bool, default=True + Whether to check for NaNs in the input data. Set to False for lazy model + evaluation. + compute : bool, default=True + Whether to compute the model elements eagerly. If True, the following + are computed sequentially: preprocessor scaler, optional NaN checks, SVD + decomposition, scores, and components. + random_state : numpy.random.Generator | int | None, default=None + Seed for the random number generator. + sample_name : str, default="sample" + Name for the new sample dimension. + feature_name : Sequence[str] | str, default="feature" + Name for the new feature dimension. + solver : {"auto", "full", "randomized"} + Solver to use for the SVD computation. + solver_kwargs : dict, default={} + Additional keyword arguments passed to the SVD solver function. + + + + Examples + -------- + >>> model = HilbertCCA(n_modes=5) + >>> model.fit(X, Y, "time") + + """ + + def __init__( + self, + n_modes: int = 2, + padding: Sequence[str] | str | None = "exp", + decay_factor: Sequence[float] | float = 0.2, + standardize: Sequence[bool] | bool = False, + use_coslat: Sequence[bool] | bool = False, + check_nans: Sequence[bool] | bool = True, + use_pca: Sequence[bool] | bool = True, + n_pca_modes: Sequence[float | int | str] | float | int | str = 0.999, + pca_init_rank_reduction: Sequence[float] | float = 0.3, + compute: bool = True, + sample_name: str = "sample", + feature_name: Sequence[str] | str = "feature", + solver: str = "auto", + random_state: np.random.Generator | int | None = None, + solver_kwargs: dict = {}, + ): + HilbertCPCCA.__init__( + self, + n_modes=n_modes, + alpha=[1.0, 1.0], + standardize=standardize, + use_coslat=use_coslat, + check_nans=check_nans, + use_pca=use_pca, + n_pca_modes=n_pca_modes, + pca_init_rank_reduction=pca_init_rank_reduction, + compute=compute, + sample_name=sample_name, + feature_name=feature_name, + solver=solver, + random_state=random_state, + solver_kwargs=solver_kwargs, + padding=padding, + decay_factor=decay_factor, + ) + self.attrs.update({"model": "Hilbert CCA"}) + # Renove alpha from the inherited CPCCA serialization params because it is hard-coded for CCA + self._params.pop("alpha")