From e9b0322096af209b2c7dcccbb360ab1e8485d0b7 Mon Sep 17 00:00:00 2001 From: jpacold Date: Tue, 5 Nov 2024 13:37:32 -0700 Subject: [PATCH] Add dynamical phase transitions tutorial (#2546) * Add dynamical phase transitions tutorial * Fix docs build, minor comment edits * escape character in matplotlib axis label * Expand introduction * More text revisions, fix heading levels * Rename tutorial * Clean up some notation * Variables for value of depolarizing noise level * Add brief summary section * Add summary of other topics covered by paper * Tags and inline documentation link * Fix indents, reword some commentary * Fix documentation link, revise introduction * Link to paper via reference, fix tex typo * Add note on choice of Ising model parameters * Style fixes * Notation (specify that $\alpha$ is the noise scale factor) * shorter import statement Co-authored-by: nate stemen --------- Co-authored-by: nate stemen --- .../_thumbnails/loschmidt_echo_qiskit.png | Bin 0 -> 47420 bytes docs/source/conf.py | 1 + docs/source/examples/examples.md | 1 + .../examples/loschmidt_echo_revival_zne.md | 425 ++++++++++++++++++ docs/source/refs.bib | 9 + 5 files changed, 436 insertions(+) create mode 100644 docs/source/_thumbnails/loschmidt_echo_qiskit.png create mode 100644 docs/source/examples/loschmidt_echo_revival_zne.md diff --git a/docs/source/_thumbnails/loschmidt_echo_qiskit.png b/docs/source/_thumbnails/loschmidt_echo_qiskit.png new file mode 100644 index 0000000000000000000000000000000000000000..206f8a38cf51fa964be6666820844e9a8e68b03f GIT binary patch literal 47420 zcmce;bx>7p^gnt4MFdooZc$OBC8SX)MOrv?O1DTiNT?tnBAudyhwcth8VTt}x;v!r z+VAf>_x@(?{O+IkopD}9j_2&XpXXWYQ|r7}lz(~!p9&v^LS2!5Ch-D=!d6G2us-79 z!r#2+rHg_82{=fqJ1EZ|p>KAHdQ{jiytZ{nv~5vlA?h)Zcn zk*B$rwYr`^qh_*xK9Vlid|HvK{bfOTAz|SgL(R^IV`Xxo;kTc* z)m-H^>3FR4`E^~cfJ_wY6D5*{h6XFI&ce*#9F;Wp`=Os|>#|l>7%?1|^Cb$Hf?kg{ zf1u>{p{CCa4Gq<aNAp!Y>#4-d$bVAcQ90Bsu1qj-q~rgyV&!qv2kmwnSSTz zwMUA1x*}N7CYe*tw%a)vtKN%!`}U3GIA*Cg zBa3vYPIY31ii*l2awl?@_d5Ei^yK7#7kyMljr#${=jg{J`k2zRjEt-+)DKmvoYpK# z4BJyzORYz<&nlQ$Z|nG{_Q@v+w-$AY$~=39NAI!u#P{_4!2A5IZl3A>>aelu+snjv zyI}=TQekuIp`uU>Hmld&-uUp4JaFzASR zP+nD~{8j2olKZ*<7Z(?PLhq;JO0Fk)SnOt)T(Mq?NtxDe z|1?;V|3I>}wN?21Waq-gi*PAEd3hYCM|(5fU&IeXOJ-{n`$SF-CcTCKP5HrFDHfT? zO?d3csi5b!=h_LW9}>c?TMXv@eoODAK1I{$+NFMU#=dk)G6?_Pe6m*L(djO+^=K(_ z!6_*>MbG{jX6db%u&CuqWc6?VB36hg6%Z8MUi=zX=zd_su2*XnY1odjN}8QDdj0x! z&^h+I8{zC?O3dcw=CE0Db5z?4job0p*uU%5iOw{Knf`4EsNNs3gbh#b+WG63#7jfN z9cyKh0C+etBJ-{HlE29(3UNNVHLRF6zF&`q%sXn0{+(%j5j? zAf2nRCQtbIV8e7US7*$5((}TF3q)jOP0@N@(r_jGj>{ySot+m$E5pJ*;Hy;XJj6r)Mm;Y3TQkzdzqn^1r(2<>l3pAV6HB z+cj!iC-*y!*Vg@PT084o9!tHlf=DBr&sbqsE;h~LDR`sRp#neHw!OJJLMi_?0vbNh z3PsKTj^}pYw;IURX|LL-4K6F=)AK&%bsl#lx1aC$+f460I~T2&ZrB_;=DLvBSj#5^ z2d$i5hW3$HLSLWCU4^8M+17~f=q$2HkDV(Vwf}C(#qk(yPW?D+Gc`4hh>0=lsi>)` z=`q8%wdeQv#|gjvF+MIbxerT)^=h`|IyMfDjHi`Gir2xo^Tqw~u{r0z4HwZFz2Cb;kF(9|D%cwu=&oM9 zT9y-Bkq+ndZPxz&{&D~Tx$NkzTenpB`1tlZg*JGs+gFY#_-qm+4@5owz3T7jQMz*N z8lSb0k`XUx+5<_XZ2&ePxqpSXSO@fTW*KV3OW6H zcikd{ig(+Zn#<_j^6F|rUe@7Qg*`?-qbRez?T=;Tt|JlgH1|Ng$5)Y8?EnV5{4P%y8;i}>_O zMN&xxPU_z7Qg8a&g9i^V5fRE?&TV%@iCJAq(0Nl+`g^%m>&*G5uyXK3_xjaT1xi`l z*2gN8X$4+?J6bE7#MGUiRUN^;%i^`@Q>3M(?H7N0H8!V=No4633FCgR!lMGSo?CQo ziwsR+^m<0axL2q)yTz{+=>NdJb^CTO>+#`2a*^xKyhqsH%3yS6s9Q}ERI&%~(5ZOL z^heeg7Mh{7cR_D>_WXI9a|i5I;~rN&o3R#Hyzg3Nx#sRW^S>!M^&i(3Vldk)d39Vy zEkr74#)l6Fb}2aZQSaZszgYS4u3X&bjfMNl>8O00@yZ*MYIOJR-RrF@Dk`e*Iu(Rx z-O$@x^jKtDeLHS!6;Ay^Fh{LE+~@B4_;Kl+t;jj8^Vko#Lb-#K!Ij}cV}hS=uXPP# zk7O#{hC5Qpt|WEL`BR;Xg3`J=UKRcB-8->^Xeg60@%QXtqrwLD!C;U*bGZDc*j$-i zw;GRRL&R|DcvSPi@I91^R19-o%{~>JOK$5tjiS51WwPN}Toje*5P@Q9y+5+F zv=mUj9^0+=Oxhin`&upJ2BGxaW@8IUmy0JG8mxA; zC&a~x=Fq##&p%XLlTY~}Adyc*BoQhh`OqNTg!0YLY7ZG0d@mA^65WX^YrR$DofK2O#bIFL*2${8jcUJ zv$MOO|GT|CxYqMH>S#n>Uf#&Wgls5@k9~WdgDp*}Zwc6XD;-yEbs7Brm8P67Yp^xlh_qLSwGq;5*RGLB{;fS;=kDKA zjjFv1E#U)=z$e(Vo4@%c)wS>0({EZ_H=KHoVa?+;>7YQ$o@$O}Bdinj-)k}4rdR7u z_o@AGuN1bqW8v4m-S)V$qrT$BYB&Asa&=A3Wjwr=GLhT2{ggB0V>MlEqaq_EwY3u@ z7ht2M4-a;Cb>($yd%5OpNrtenbru+-OF4ageT|Hb!{g#K^=YE>uOdja=488_M${ul zLG;l7@Qd` z{MD*GOTe4s1oU`N=^ga{D_45S~&?_M0rUwXM_U~KzPQIV4aJl~e z7lM~fLVxH_`Cg*C2|JBVtCU1gz4jaQHh@jDCx<&1_eo20zGZKDZSCyL#F{4kg+B~; z?*=-pjciT0Z#*#TPPtJ-o}I>zY=SD`jfH<(GXc1dVKd%^4*)n$HdJOg$R1poi6tg3 zE_YsM(_d<>1^7?u#fulU`bHZQHR-4GBnv}@%vv$SajevAUB_$c(6U`o5nk?!W7edWP#w7H5_L==V*Ve9ZuX%h%G#2 zkN3?Bu~qWYdTKSgdMUyEgM;MW-rmY;^3oFo11}QXi*`?Sk`C-)!{q8zMR6)6Ns~SlN&Fh!R4rpo$pQBh5RGq!)fDcwaL($L1#V$auX zjiNvxKJkeH_n#|h=KD87t#{q+dD>%90=?#KnoKmAW-g#!*!}&huDaf*j?kViLBrBo zRF8r(3y2~Xd$h`NCCc_yZrX(BfzCtQ30Y{SF7vUb-AfJ#T81;r<-9)jGMS^RBbIwO z-vGyX%2eMm6NwqgUQ__NfZj*}VautUHWTQ?<{O>^B!6a zDW1;9+xA84x=F(Yux^fwVfY7`TUX*@W53sKx3O5jiX#B9zp+u$x;H)K2)YtUW-#MM z!1UNpe4utDJJ`KX59Q1D>_t zo^2ItCU_rYJznWJX5S@31GNUaMJyX1Q~_PDL*p09%EoI?mCTCwykzkKeE;eQ%*)WspegNqK?;Wd1R98oX=Q7TP)L zeoagy!J76Dpotk3hepx@1nj0SAlPMQ1kXU{hW=#+<`6hn+fcoTNl63z{BU%toQwyu zH3IJbB0>9MFbTr0f2JJ>#8gP2qW&*^#)vQqk#PJYrF|F?<;K963wPx%bcuKG@W#f* zvMcxFk>4|uqtVjx@=cSIbrO^tjP;e$5)u;c2hxb)cM|;$hud@5!5_55u(RXKD=M&> z$dy&e$>CR+JYZwMCFBa8*h%|1MpWB80SNrRq#&hf263(yo&S?e)C9K^B$qcBk(g+` zV)5=peo>KlR{uyPm6(bKa;vQA7-p#YckbL7s$_1!7Zwruu2WSydrF*2_zZ3(UtRar z+{#KYfX)+x+;4zU0c*v?P@qmuPOP+VGhpbDt03!H*-M|oNQF?9=7usc0!&4{eED*7 zZ!dCPi~-cF?AIS@uRoFZ@bFy&9d7mH)LQX2LlB2;MYrq4Xk6UCzrW8jmaqGURq6Y)GW~vrIvIvG`2e}zP3Rh&0Gm%vPnAdl z`m@v~_hUdBLg*h24UJh}=5^3)?(euw$iwG;`9bOv7l}4&|8JF;%kZi7)6&!bq{~r@ zo}btUx4jJvSi42TV@8DFnDP_YRg$1M90JtCGMlIt5OP|bnhASk0u9p`Eyp0ixCuPg z0I1`N1uSYf{^%RDc9ru+y6t4`Vd)hLP9G@MqnT1tPj8C{k9o6lWqUB+lHi?RTMLPg zkI%_s@XwUO$oi0yl1k=imTX&R0Wp=9mPP=lo2-e+UHJKsZlH8|C{J?zd~fJ^mYyx=L=V#pB)1!FOtm*Pu0CCL$H6;a4`7Du39#aLqelX zMqOj#dBm(kga2jN+4*oTS65eE_J>X7F7xetpcA2jYZ%DW3w!(a5`s58k9Pk88T$!( zz~|2->$S%=0Dweb{r&+vqvSHegEp83lx? zvN&iSz|B*a7XH^M$1lnQ6IeVr7Cn}_PdI4@i-7s`thTW?qRm< z$}I^;=c8S72SHlc-F()T1KE^7*wT{$SUiXSe699?VQH7Y!s5hrK}ze|NQs1$)Tf9D zaUe0#>^jWR$3$W(9;6TJpG!;s1^AhsEG{mN5SfXIiM|}IX!G85yS-&}?fIEID#>m6 zzUc8Px31?Nn$C6h`q|M6hibO^4|wX+D`BbsptTXI^ann z2NlT72P)nmR&Id}f1-YcfGhb|Y4{%lj0LgSaZco9^CzoDp?^zDi^n9{3T$Aw8rQ!8 zS}SXJ=n+PCZd3{|mCQ^s{H2k>-QBv``fn9R%OJ zd4n&~1#}b6-K;%dm62o{ zwYTi_H{j=^#|@0nLZ`p_N3iZYm@3aJEaFbICJJ38bu-f~MM^kR4NxMkS5uFZJhr84 zPqt-XrJ=NW)Y`=X`bh(rK}au}0Nw7-H)&#!d(^9(*x))+01EQhjKu)g3tl$4oqvg|w`3LV${Gar(1=zgfH zdnDq%KMRL$4z_ar;NW1fo)XDGj#hJ`^F-^f_r&IZzFh`UiXJ5MbQMl|&)w^~y1EW4 z1ITA`>er)SDNM7pZiuNo4$@TNR^pBL{MjGQHegFoCZpFHa}ya~8z`)2t%Wr6Ua=>O zZ=FG1IIX;qF*Q|~6Gl`P9v=S0%!~=pb8}Bm&p`im)thz73>_@ivNze(OJ5Y-_o~q~ zb&_Qb`c;l}$jDn*QEZN`f$Y5E6p`9zMgQ3)N=|p`3mHSmjQ^2;;4%u&VXG#U$?Hws z*=Ak{x(wz?Di8B2YlFpCDe<$FLqd)|-9Ss5(u4VELfIkAr&YW>5yB0Vdp%cTC7Gac ztC0P9gT|^(PH1(!7SKJ?)ZDaBbtCZy8_u%s+w(P+a*CY)+Bo~=d}g0-#oZo{R!Z=N z-h||6h=i?3Kd&+iw)pQO{qNLcP;cySDB28{)^P)1OdTxA0uwUBCZ5yW&Qx|F1iOCPvt=zE$7{1u<|1bP$kRX-gyTrT!c{J19j; zoZ^UbZx=U7;o0AnOp59~x7^AV{uvW(v3t|a$^X5S5t<%-AR*-&*b_mi`EGj_T58gG z1~@nL(IgDl>SuAfPc>GTi!<~rs8HU;yz_NIdNufHTdx{12Ih$FA(eS6{-2+duW_6H zq!cq$w$@?Hl9rQ8gW~7`Z(N!iMAQ%nxBbk5$G?_)>c`|_a@@uGo~x~Q7Weqm(jqF| z%dA)K`IBbzzls&5w=w$j!U*Am3uwLOK~ZfT7u#d22UzW0{(tyQV(2o_PT%X{k; z#$Rj*j$vt@ZvNgPp?kbXflIoJ%R9N#rx#`?$kmg&lgX%l(-Yj9-PPfUAASfUdl1KG z+X)TfcY=U@jDWpKZ6OHS4gLKkEzLj-v-i0088A0Bm^0oK+H%L8jqY!MeM^1pK4|bO zu07LzKSe=aYGvMJW%)AJW4F>IGF$87whBUK!Xu-%lEw-(A}G`2X-|(DP`aQ+#;NJL zUi7h>X?kEY#vgcvTI$)eXSIch>SWmT;r_PA*;%1SBvFGfY(dOjoTajL0k2)3V7WXy zQI-g){i}DZhkOPXdrET?i6TKDR~xEv=L5Bo+q!qjArqsSvWH(kI!bM8yRn^}N|YKD zj><}X{VGjPQO_b{_4e~i@s?G8nAqbj+tmR6F%#YsV#}DI)3-YedwAYFO>3}7ur%Gv zj8E)-Ve4hxT{4|Ae)(k4E02MB%Z>sSmH`G}){jB^Zo}hw36NPbX{YbIGPId8^*uOU z^st-sD+j((If;_dZfqr6aNu{8E^ugl=H`bBsnv0^O zk(Y!9S)2xw8E+610~cosdf;lStxW^>^c8kjKvXol7=wg)>&=6pY^6s>#uRl+_xgiP z)J$+u@9WNO+BGjLz1;KEj^QrkrMg&5#W=Z>5T0{++wOkE3h-+(AiAJZ3)s&Sf~I$y ziRm5aP1ot@bW^IqECBg{dHYK*DIH$=w{M7^A4(&D{TK??(@>i5CtG4F+;BfePNBsX zJRf~ltyJW$j#3?>h|WiuRs4Q5MG~hPNjc-q=0&@0{Ytz-F*cHXZ{4b1k7Mtp!`phU zDwD5LAuh*jWrs)~yL|b_tjsJyQ-?HoGe2J8}8(6_s_P%BHdC020zNk?RqV~ix}dD z=4i$REmc?NN3atp1%eHScjZc&*Xa?zpx{WQU<3ZFGqb;ITF7N1LRHDZ6!O!;Vhq|zC77D?`_KcE#@Jgj8&j?hF|Gd7lE=Y+?*vR? zt;5D3ArOJex&b``Z&pum`;~ZrjObLtMK7rrZ^B{dPk2!*^SUM;EmQTP)7b zj%`6q{gou*K2+_>4aV8-#$+;>Qd(+qcuJvAN+qDbFy$+Trfxp9VHwalPJLix9O+); z&lL4j=5Nu}{2lLCuOqh(i>|)<*UVM7`~Oiqt%mY_!4f8-p^*hSLA}uEGT_GPsULo- zRe_uwLPGHX5tXW(tUWwDGJ$fOot<64#YG?syLP7FkLt`J#E>VFLq$G4Mb;l0@P` z;>^{nlb4kx1W2J>@j3uB$s0FrGz5{dSm~gZt}tOV-zZ?1neX4f49YcB5<6$-GDk=J z5iURU(>epGD++|Dpf{-)ruy@<=Z((ExA#M+{^avL%YaJ~H<_|zz@$yJwmQ1TI9tW$??w z+=7KC*LCp)FiPVC7ID7SFc(O(YAHz6Am)J~E|mCb@UBci>jjyeg9n7c>MS-4a$X~- zt{5zu_Gxj)X55HyX=@F42nqCArQupXLn(AgOxNm2jFkmLre*GD=kQRR(Rb@2{jTUH zn;%;DH5a#ecF!KQ=V}pVC?w;ifb^1Xjj%wd%f%M`n8$|z*dG83y4B?h;@|Q>j-}H% zSbc|})_^VN4TOM6)iYquq$9=vKNlWb8_+Vp#6+Dh#*TutATF{afEOYZv!JT)t_+3& z>z!s~F+j35S%sHfQOizWesd6aXKrbi^!C{3xto;k7;pYF_~m(yV}9xCN3UHbdFzp( z%FfxwTZwuqjR?#vO$UD0;?M3b+Ywzz3$vAIJBJrgSaN0$&b3e zNk{yHI7?e0CY9=wuVAE7NMN~ByjVoG6>6I55+Z-CTJ$#LvIGSPR_Kkwk8MND99s*+M5^l z0HP-=JUaD{kEaG2Bq1pY!maMz167^0w6q6aCoW*Jfn|`=_ye5;44}wi=>=E@kIqlH zL0Z@^QkUGjg{5jt+U@)CgDgsWg>tCUS_m}}CVSx~FyKP`!=DWPzXOC(5!W6cm@>Hd zZ@F$`&Fpji0b_E)>)2ZKaF$q1OlUkY9!tsE(3c;gCk9rM3 zEmN(LcRV;+$=d;g`tQa@*ZJ8&-HpVr;M2$^i_-d-CbSA{WC%Z=KFc_udvtr7$%1XOiKi5}Tytv_Vgb1s z*#tAI@7_>cWX7?-_b!3^NB_#+eg+<_F?h-YndwBNeDz=khq=IOVj<3e#4fV(Y8!y% zO9tc9YZZehRq4FZiTDtZ7XbswbgbO2%z6*gVR6antGv@y47N?q$sv@CmiJ+AS69W_ zvotiCPN6wVBR-vfK!6g&zeGxRc12|zeu^CQt7gF3pDf@5QN#I1xF?H?~RyHH5&-@qYqc; zTT5>HXI4Uq`|;)?0ZP5v<$;ZjjjGTA6az}|Hb5AcO%&|6fDcoH4}%O2D4Ts4G|S{Z z#G9^@lA;$MfGhUMw(j(yh}&*`XJ;rN^gkd1+^WLk=H`|LzyBapl_ytZan>!o zeaXqBN-=;aweT2hPZWi~>pL|$2`dUpBv2gQ%4ZmY#j-k9(TiCi2jzn>|M^Rh;vwb( zUjDx~!11?We-t>aX->E<&=xf#;u|2o7&`De=hO{Z@}0jV%GjieXJql?mPDmn14RaN zlTkSrQ)pp;QOip*_fRvSJPTk7#i!Ii{P!=V-vUe-h{%Db2M3qjtpmi9R`F{hzYN1f zq_R|VC?34LfYfzZ#Aba?61NVJww}q#@`xaTrakxjmRD9Be#qI@r}WHF(bV_cehNen zP+KKviZKZZG$5eLfj!W^D)?XE3~mo#`szqYR!8n_Of6i)@aSl7t*3jb&A7#g2`9x+ zr6UVec(&ab=y5z{tgo{BEwE;amm2u_`KQ*`lQxuw&9H@tKM7rJf1ud)DC1nbpp2t9 z`=xD|d=$Vs6)v#FA<9R@R5g$#{N%}#^v_g(f*7BF6~k3WA#oi7a`sv91!mxY2%qlu zeoIe}i?eMGqwj*-lAU`?NZ9~crB(=H0U$PpORYYfj|Zjn7l|2K+Lu=D^3$~h@!{yy zuLF53myaA7b3KnyQ?8ThZ@hG*0S*^%7gfK7t58YyJ`?!<{X4uV;vhyJ?rm!pJhUD$ za(Wp_RHke&c#8(icW9F{Ib{>Kn3;p1ffgG4#OQGk!ly_Qz&GLI2U`iU zQgT4V8{9W)KwDv=*WJ7P<2uTHL6K&;bV1*RANGD@fjRe!Qmc!k32*$YwUrqBnQGwy zDyQ!wzD9A9kp|{}Qe{17v>A-q{70&u&&uE;{macLeeY; zqzg#{{Bf>r#V80ui0YnKBYx^cT&QK9?gg*Cw4AK_%oQI#5Q5!|#ML0v2|9oNo3>4% zQwV3o@R;9#XP*Gs%T=n615R7hl5srdZy;@Wus#mkF-<#{g*o(3(*@9v71cW#@`-;` zEK$*%aQ|~Ze}HnRtwk9%OT>$)g871!c;Jc;LsnqZuz=QA=2*g<1{(_w%cXi5jv5u4 zAEXo!VWISo=sp&xry2qB#K#tlZ89)VY~UwLe6O3=%w8tMdr$MV;rk#);rj`Pmlk%G zC-b<@r#u=WFr5vvSXfvHUf*ePHenV0AACJ%AYkt9vSFq;pXlRfIxe)RK$ z4mW2e?k+}ja%VTOSAg1iAa&=n(c)$(qUu{TAR@Yj@`J;vi@bm7;pz~1sJcM}{0&-;d%7r6BG z`{=*2GFKg()$os!=1C$U93u$iR)ZCa`asPWipVBbL!1B~wj({4zh_Wzo5+SQ>~Kt< z_k#7k1e*%lJMv5|`k6J0&6Ki_RGr}rUQ7Wi+k{wll2z!|8v>k5UrdZihcbpnqxt2K zOr;qlfFXix2q7xr{UO6f(YKjN-!5ih%FSt_It4!bR8tN)@1$L1pUIbSbed};b=sKt zjfB4b2VQq$b%53m=5}VwQze{S;8<oCbIic zvU&DbGdYzu7Iz)<*4x*c)aU)2zF zT$X^+kN-It78i;?^L7GUPb%0h0qaF&Wf_nvLKWLiCP7EvfgEZJy&1vUSVcs4w<0)1 zl%gKxIa!F++yRZf6ObA|u$$V0aH4ro0MIV|_aq*i_4M%TG!{&po>bYBi&)O3TfB=x z$*3o%GBqw^sCyj)|+ zA3X+qf`T~0XE|I5^g(WQ49X>hn?67+64=uxHN~+chwG}OL1K^BuJ(1bVm10Cnqsdw z(LY^?AAG~;Hx{tgZ=cuPVZ+_SL{%L0aB>HJ{Jx^=60B%m@_*EYbU7s@%q#IStE z?TV|bQ$RB@1x$t@`Ac*zDYEN1nAVmS{l*|tpSf|rIflG=McgU^8EGGM0F$Fo&qKby znF9hJaH1n@9j-qoh|0Vf$Vb&l|MOg5Z1VzcA9cO|)a8>0C)eO_bxj7Xz#KuKMzY8D zWd;TY$Oadg_uhqA{Lgw{oE*0tFk&I0r~oJ#(BU;ChzlYXkeV*c;CBp!US+`KgL;Pq zyr*Xh^R(|*J-0K=cYpnw3s473ujPBplS;g-Pd@U+qvsDm^ym?MYo{|v=a+HmCRS9Y z`f}yq+~Iy#^a?M;Uy-o%Z;#Qdbi!tPD;^+QsXj7w{4^-?uO8>xtd?ylQ(&LUOVN`S z53C%r#Avne%aMFTr;%#922kC&;uXiDVB0b}j$Ho~|4LpIH%X~?77{8_4byjGVjizQZK-n{&Vfo zV1&y-EdIxjA7r$|P_XXgO`-&8v68)lEDzTCFJ77;Lz3!N7`7}Bz)VQ_4$vb% z{3N{fwFxT+vht?K*u5LnC9~4zjSNS^${A(4@7?)a8R{B4?tQd$4+Mtio}x;QuuU3I zBIqQZHHG-DqYC7k?9_R=KB;?wTY|X@@0$3?I|<Q(MPh(BQ$W7i`zs{+p+Az&ERMeOA#~@#1iei5F|2p~<$f2Z zi7XqHH*H)`CiZR{u@~xNw0BiqR|o?`XpB8RP|8TH{}C6s9==FS94uUBnJmg*YUUdk z)#S+1vE5;5a18t&&^7bkp+x^3_yq2BL!DaZtG_-T7j^FOycbnM{afsMP%56X#te@? zp#1e5{-tYl5=tb1q+aNIhp!Ntm`*w+K9N+~balceQ$#+(1XePWsR!8Z}q);$~Tew$e>eXQ$V4KQHp^4_d1W0oF2u;2;WApDUAN7ngH> zmzXNhj227%Zn^37x8X&ij1q}5OOyjKzCASLm#T?c0rOFMmfAKx=-ZTDLI?3zEFd8= zKUabcpe(!l6hk^-r$Q1C*G0ASl95$Q-bqE}duKGcKBy`0Ez_*76Qd*GdSoS@Unol_ z4zcGhHF0IklK3X^**ncY)2P&TGBlOx&X|d&H*hW>lpN85k-s3}Z#G(*larePiA7$G`f%t5!tGl-cl?*oMVLO0Ev-MK5bZZ2b$MlEpnY2=Je7W)t1PXQuX$3q>pAn19_D_@Tra$R3!nCL_Re-g=j0?VFLTSp zL~XYURv8AwAc4Mgi(}v`TfIv$Lq+!hZ(ygG7~7qWN`A8%#`ZCMh8GJEM6g>2sKbDu z-@Y$aQ&&?0O|I~%gRedBSpLNH!C#Q#T<0Ry{sFhE@Hpi1_2qj1Yz!<&v~0~PjY*U4 zQyiuW%WIg&Z?f#XbDnAJY7EKyF75M)uoKd#tdd=jLGmFG)ZPSvVPxC0&bw0YQi%0KVcfWj_>QPck z+N17QKizt%{P-RP5fk~^N!>CTFh&;@+6dj=GJ#Hl$gMthvGfc!2nx3UcPTLA+fH-j zvvqlCoO%uPRtcY%kLgClll2BS6?EgrD3h(*LR^1N9FDVD3vb_M6I75UV?PICzPYt! z0Er6OSgrs_8^-l~ecB7+$5mQE3D`@>Jk~QAnRaJw;Pf$%yhS1B0g;9ql$5QI-qM~r z0Z9~W_LYsg$X`H`|&+~B_k>0%Kf*<{{JR!I_5p+0C zvl6-`j)1*W(=+T#b=fwks*`_q14U^H8moT=wu3dVN~`P_)rB*vcYEJiv-|vusE2fSNq!45oJg)1P`eL zD6I>OK0s0uOiI#3`REPJ|IZYQ{@KYP67@&^psFDmXf<4zCRX+m9YFH}gu$0>pYMYT z^RCFW3mF9fXIO;@jPib$fr}pu^Ac1dn>3oDvZ}S!1vEk*Vq%C&`$4?!tt-)Fg9(ci zbe3wmb%{kk_RYI@(^7Ml$Ib z_Alf{6k&J_aHd=s;7xYza#TG)oh7jBYG)uyqnYY{AJ9VFwjSC~=p1ZZ}O3)iGj0*U(2#4^>x4 z4zjTe300jeD1P=B-8Pw|BE)TLTOY4-g(NpvzIP#2T>tZ@A0%nF0DF?TYCnID4=G1L zel_*A4JLf5xNon~&%-R4|7>Spb_eWoQAj#<0e&-~O1|e0MO-c1(E@)85kebBPr}3IU z_rV9s>wv^FGF4bxpaId0vKm!zVBZA=iDmR&Cm~T;+=~ah8V(0CWnexEgwDp@wiLkH zP>!V<-W8UzmF4Pc#nm+edVvkbJ}yDcF-duqobYnL$^odKyh$etr=PC-GzMcf7)L|F%d ziVQsPoArnSfG6^0|EnIgTV@lTWd$hH;`r_ zW*rHq909yAYdQpC6`XQGpv$PKiiIRY5H?tKt7AY93qzDqF21!VlAp=-9qD=S90I%t zIf4)T6YIQvDc>i}taIyj{rK^DBAp+MqowgGzQ?|k8JV@voiUI`INNwpr%J6&h?`qQ zT_v(x_{}f>OT$@tq(h~rE4T~P{~7Pip=nSRIueD-9#g{P9@roNr@p}47(xfrqnMk= zqN=`#zuj&i;NW!KHX13lZifYO+4;@4lW02I5>80N{}Gl$8I&EbvqMu7NdqM43Q26} z58+_*y_)nq0N2^Zf{hxKA|%O^S9`<^lI=ZV;e-s^)*Tfk?MjFD;8mpdm@PV=pO)|O zg&J&Ah-Pa~5SMj6%RW6?;V-obiHm#w>-?17;&SDCC6ZOuzw|GJ0Px@YUXigR-x=Eh`N{XdRZy3Z35YaqYCvi$&G^unBsmRbf% z>r|ecRlPn!K~x&5^iYNUy@gAJH-%1a`09BFSSpO(5$~*~c)nP3kZNM2xDvv*>9y(m zl9ry5QaVx4Q4t*pIyEs7k)jkY>0Daz179Eqj=qWJv9Z$WHDWk!0$GHEr2@nHEn-cRmeQgf5&ha5h7rS#Tt$}Cnm@QBKZ6VXAMn*=jKp}`?)ewVp zft5}UzA*&d?l3a`1q%kGU+MNl^QEzZuOXJWK~jE`?i*k1>~GHuJ;>$IR?nV!`hfzA zg1fJ|sVPd5RFV-~7Tzi9p~w*dfpnv?z+s@BPth``eOIpIa=%uQI?Mq96=p4wrzz>e(NdtJp$4%oZ$m zHcgBVZ?NkY+XYA%-`#nMxvF)-zSE?KWK|PZ23y;HOKkK_xXG|bMkycfr}X&fD$KIv z&5wd~LM~8+OtZ|jgx}8AED3@M;cvni7G@^%|K7ZpnLC1C&Cn}4;RgG~TbQ!y?1g|i zpBs_N6C2wDtOv$03K!*fDJ1qBNUZ`q`a0S;p6`kR-Mg1B>)%ED0d$r3s2(=lX)-`? zG6eVvGm|p!K;akMw0;sybLD*2>3XG!iQU==FO0mJkGj2oE}7Loy|iSmnbiem`sc^W zD(o1{r>H0?@Qfi)X;w1wx?&yLsuDU9@^3R68w!1?jVONSOp0;8f$mgE0?;%YkccQu zL{SJiW#?qwv68^fUWNq2cLmEE_8)pfEr&-7ywAOm7%R1!tja@I0?KHICq_yI;ONe4 zB~)gQkBQ)YjEWDKsEAm*RIU}A%7myEh#U#ljU0Bi1w^GS8g6G&uWv1&u83-~)vb-Q zU_~)K1lF>%^K)=8dUe%GOYQ0BH(+R1lxX`F^jiV>_jyqN=ic2c}~ zWQ$!=N0f8Wjv0#Oc9BF@01tOKLS@@qiFPmEuM^ozy zd`pKpW0PKHNZ5v9&*U*HzUk45)%sk$pM6l@8CtAU^)B|nAbbHCazcTpbt&)$D;*ni z9sa0fIiL+)k1fgCqh%Ik_yhv6Fo=l)hq-ZZu75c)t>%;ddz(ZcVDpUHn&%GP*c zL+|Ebr1C3Hs!Nw{>AJ<1HwHt<=~~_~KK-u=n7WRL#-)aTSEr+B9$(D||$H|&Ss7D9%VHmb$j zIuiq;K+;B;L0jYRalfF{kEMCjPS-RW7M3#X z=mjtAy5$NV(@4WKiI)65TWMO#x9A#l3wy)RkUC@^05i4$w#m_i)LnYFAxxmQKx(`Z z=bCUmoI-w>!$U@=;5dWZThnbg0aMB*oeAc9M`|iTYDDW-n_??d;ufGlypKSidmgdu`Opn3z$-Hiy~o9g^cm6_jP12hBn37Z$z|zInjBFEjj# z1&dXhSgC%b{||jg$VS8HXg63`O@PS}#S@mtcs~a0+i(~Mchj!*sGza5!U=lwRSSjE zrDTo`l!Y&2FrR_sC#DdTX1nz<0eF<)_BVY@l0SM+{J4?U+xT)wjYus=jhRzjTx39z zt^#BavCJDWq0B};b2Br3s#%PJf`XI#ka8)~czb$sVg_O0+A(1mTDb>y60DP07@W6y zWcK_krnlevl8-Y!Kdt}}6QWdE7^)%lfnnt7zki<`h2Q=d;n&%9txC|*N~_|T?FTe^}KJaRk5a zBQoR#Hq*{-uR=szTpI+9u@KuEY8iJV=wTV}gM()P<;~QcpBcghAbu$LA;ZgL%y;TO z>dp>P5`TF3eWC4a>>B7X`47=(IV@F}L_{$=yy+<xB5p#nz)yx1WYWM(2MlIFK~Vndza#L3C%h&ifHqQ1kHh9&^Ygz-(+$Kby+4 zlIrnRI}rKT*IM2a-+Lw}S7MjrPlWFIuABkD8_5&^;(OWBap!bVDP{ZWD~`*ta_S=J z#U(=+1~7pEr+{_v{M~)6kzt~+>;knrl1eR3PCp+!xO_l!B`CfD!0uAd*S9%Ykl%i} zt(k{pKWJ!arAO=Nd07ao<>*)a-dJe~u(W8sc*r~(4Z38%XQoidya9~lg8vbR)NB3=cOktrCs+<@X2P0v!X5#bZ zTDT-6W&!CkB=E1Ls`2^u<$TPSRlyLfO?lNVi(5G30Yn4Q5fO&V=H(C0FsnOy0Zr*# z5cN$z@rV(i?j)kpDiyy8gFSZmSDQQ^d!q{3b${LG=bz*f0^b6+N`}v1j{!Q-EL9&7 zXS&6v`CSRB`z6`S98=PSZ;n9OLxvHMiL{sX}G*M5sue%Bkk zBD6WtVI4=#Ur)sf%Th8EB4#!o=+cS$vu^%f+>-7YLQGp^dGQUs!58wbSgH(c>OLQSVoaDF;&7E7 zI*5u=kX+2c8dH8ks~I8pAl~K9()Hzv!6L^!f1=1SNV~MQwcWbI`LRCR>Q5N$ZR=!! z3~ufa2nn|(DrZ!sDethXqZ3ymQ9(;8k<^|YV>>_m-ugaE`d7M9oQNvg;Gl8EN|aP? za=w-+NX*QNndPeX{BsXMEM55zOGM8+>7!qaJL~#H+AWg=X0!5?W#4)YQv)DV^SkXb z0YUA{x|dvrjs=u29wh&?GM-ysw*T=Jp-=I5xgGiyRg?9sdn;k=t_}{~bVl!1Vla!r zSE=g{O;|0Ya8cJe>KY)a(h~YE|DT6sNZyFoK2^c}SQmj&C6#7FtLJ24-rhc6lGlvp z<06c?C??oBiE_&LESw%vg6=n@2l@e&udDT7Jt|==-v=oG#$p0zn90=8zQQ>OA7X!K zPJR8O+$HVjl}0mE_zcWN(9x#8_opfoY!oUN6+&F4?7vUOrl8P{_slymUe6pbBk8vI z?9s&?h`ZyE7sx+qDd5E$-dNSpX2e;Y%EaeLOvhII!GTznUiDf&2%?Cr2q1%zoxO7; zgIqO#d5n4bHM#GBd}p145*l=2twXivcuI`$Wk?#HQDc$+98|b}#7SWZCE^O@1oqHk zAFib73fr0D;Fj$#yceS09F(bv>1KQP9VdQkBdf-odeMAeaDakvqc8o6kGXV;*4#yP zP>1e4cu=4+c!M+^-J}n5c}P+V4BI9MgSND_aqWd)x^SVp&f6Q|Rf4KtOEQ=X-<6-& zljzT#cfUaU+0_w=eXfUo0H7Gaj`-gPN(C2P1Y*ZCn4Z_3J`rV?r5(!dF67mw#oKk7 zABwzyPkk*9!>l&lD}x>!z2aaC$|a&*z=XZkSb6vZ!=}&6ZmHAEjsyen6*Cf`SR)B3 zAp6KTIaK6Lkl+g-nh5j05c0o@M2M#9eUV@R6?`$oR5FS9o2;$LiTj5K_bZMeQ*&qg z54QvIH4qEcS+TZXu&z`1r54a)vx%ztL`Mv-OZnknHrJS4lu5$t`;NA0=*i3!BMxl{ zQDud0R}v{$KLWKLzR%zqprjf2)&e7#(}Rp*5Ja}1tlvWhpkPQwRz-yby#EEbwi2lt z$@?F=M-dF_IzA?uiOH<5aXW6fN-Xwnv`x#hOA{*!gYnoObnWdb3(bpz%=weXm{Y!V z@f5h4|BI}@46ABuo<-Jl>SASlv}ba!_nAc!K}AR#T?DGQ{L2I-a% z=~AC_?fbfq_kE7%lVA4Mwbq<-%rXA)JN^BJx*Uh{MW&p7dxpo))$4XB)e#w#Ky@7g z4w2g&ED5jNi2j0vOcv*4Z@>6CKb>WFEtqud(6&SDr(#UU zwbhkwC}B%Kd1AlbPe~%3wW8}rba{Mjtho5v8i|*evz{JV+F`y=!r5i3)zmgGn|Fbj zm*mRx&!kA#8*FO!7%vyPLxk>W^Bl*21PzV5xJ1G)tsa?FV7#YYj)f>O+}Y}v$K^6J zwvt)M>*}w!(dK@U#`tHS=!a8m9nB;+fdD7ubXi}s)q)!U5^59 z1M)kz1EV-rxE^oPx5ZcP4vs0la^RtOHes3MDHfge&En=kOog!_ATaV z-r;YD|2dGq<*ONPwa{`VYIF6W>25WE{a26*+h4rw3v2DT)z7l-O^Q`rV}p&c8us53f2?aTc&>=p)8WiY3;HMyIt(gsK?TCObX5q1fII7NRzZ&K5ok=O4$9N7iDxRPH=q zv}07alLaXLjyll8Bhydv(4&5ecyaA5^10d#c~h$Of&2aPf?A19wjC%Fs)s{f$wQN& z4O{~fP+Vn{NDK{e2x1f!6i_>uiMo&Ffl#>F6PbC7S3--3TCPGA*Pd?Rc4M2gFc;p; zw(P|%N99^ZYx;*DB6hYXIVbALuAJTnE&Z+&-b@r|z^8H|dO^u-$(Y5!gTMQ8LV9~G zBsooC{j+jN{MJIN1ki<`P#6Fe^^b)u*@hKz{;JTlhCB0xiYltn1I7FAW%?4R79T1o z3hBfEpbGJEFccef1@tjf;eUyyLwxO`&gbnG%F8YkVmsd`M!BK|!4M$A#>JJ&9|gsq z5TOzjP+X4sgjky1rK9gEzzC3OUgI8KT&may>~U&>KV_d!shzA zDb5pyE%R~XIK+0&J8?r0#$mG#7Ds1C2=3N=v&aS%=t0hB({1<&l{u}IR~g5I_&gx&2OxLj-3o;=OT%MVfv^LgYiz45Jo z#oF=?>Yn%PBf#shTo2rl#4uHas&;Ek@%AW=?OqQWtO|Gz|&8W^kM0&4d z!%>CguC%n$PqJe91lwY5Gq;p~j}RaVv+puUKUnh-Y1}hT-QD__F5Qf0^m1+2JG6>e z+gei$?%bTh!ld8o^=K0||7})Km8*^IR!)0yWa1ye2VyC^W;O{&f(sQQC0F(i52ICG zsz$Y(n;5X4qF=kVIh=(K2L1ODmVz-fHMjyEQ2mEpbL%uWq;pZbdw9d3phod0AF@Bf z&#~r0G|}_yt9a+xfNWWF&wJVw6Sc^j9fpwaW)v2tdiwNfgY%Z*Ou7lxBAxNc6x`bX zb+Y2Tr)N~Har8JOm=TutWn1=aru47K2ue^*1=w6uUfljWJC`v`HMb&~XGf z4%d^cja(Yb`-Ni!t|*xTyx5!B+{4?)V0zH@+%}-(=I)8+9el>Y(uT9|*zv>vy8NrV zAM=0fZljf7^r+q%t+dk?x6AZW!hT9B&E-PeYu!=URhaa{4yCt30C#ijGqH_hid@L9 z$iXK-vNGTGmwD`TtV%<4gwGxX$|AYY>>?F{>ZvkTo66yc-poee>u~em3 zmI7+txL{#b7f038XT;|{VdkX$=`m7VLrpo2sYci2IC~I(pC*KQE*MaD`D}l6LO#m+ zA48*IAB>3=Xe7h22OArE307V9Ncr$vu+!soMpHw0i3xfLqx=W}C-v2vLRpEi%&Q#s zZ=22nubL8(`WP6bNtUm-x2SnAs;Mxd99KLBdnUT{`s>*qM05n{zAe#^KMgi|NFqkr z1GlPyPIoPAD!z7i_bv3Ye6_@Wmz|xB>e+y(bC;L*BZNjP;noDhB;IZ~r{Jt&w|Wsf zjDR*FEp-h;WAE$?sn=A{0W_;yI{1rkp{nwcU`FjjFnby91>lZ!gNak@VEGz-R;yiFe^|ak ziT(-EL3}sLHKH1bl3*$`&PJ6zP~sX~X7lL^2EcTp5fKqh z)V(dgm6kvr(d?6TiKi7p>Sli48ZONEMCI*uKgl-{4>&>jx{4F$CzF!HVa9QF*(MX@ z5Q5*pl>)u|fg7KscDPpA!D_CxM$NxlnnIy>`-H0pwi$UP2ypKgd8rz~MhB>_qCgw` zA1=`e&S{6BIHIEWsM_IxeE6&K!&B77hN5;wTjx^qhLSEm1+T$=!Z&QmHyyyuK&@7p z2L5pZe|t+6X?a7QJ6`y@35IV?^A>M)h?THz2oAQeTPpR<1x$Tp<4@&XIL3-+puwoy zNcV!2Giooby`iq6XcCWA1iBfxY6+mV?t^V@V3*hRWGWOkRIri4pu^vzlG(;H^2977 zbFtQur+*?{ZPZNwhpErnwH&Z~_1FGFo)rbo1wBX-V!Jv;P8Ib|d14skk3TGp{o=A0 zHuqXOo{xb>%;URthY|wnc&KyWGgdWKgKA#Bfd!%0Gy#C!slCMLAbSdmJ^Gh@H(-18y6_xX>2|#c%A)Rl0`aAG|(ODQW zY!n}2Wwjc~!Wfa&`nd9LJceIyxb94~gU`ugb*koQLT3=(oB!GB-krDG$CVT`WXRt@ z6x7@2`~2@Jas!=k(c{;RZ$S)C_!0#E7Q!J4<0OfCF+Dk{v%}9(yX=C*X=i*Xx7-vb z6H)lvV}l8-=Nlev2i#o`+^DO6&{`~Q41OsJzj9o}i7eGeN{t%m+4%eeOf)<#_eK>GIXL1znsLLzDp->s7R>5O={^ywD^_S1m zrKw694Jn|2&(zwzgX2p|`3ImSlPqeMVc}Os;U;-H+akfe|59d^$`wn?fzBCHum=Z? z(85AsY$?NO!~1d!QYobrlUt2KZF0>PJL9YKQNtFIzc)7a384;$lSQ>yM28HsH095* z^uHb28JJ%FQ}$tNlk;9RTjg?_39sp`K$eJ(19XOza512tApn0id5e~l1ryM7k9}>- z<1FyXE%}^`A^AUrMG12HzPrP?ky0A0;oB49EGGo1ZyqkWl3n?{AASD~@#Ma7`Op3FW!>RQ4pZ%Nk|+m%&mO5v z^MQ{elC=-^zguh`C7wGX-alQ;A=_UuKPyJxxw@UIW8(d{7@g+thzf68-?yC){1MqV z7uSrZ|IB5Bgj;1j{Q*w9aPCshcDgZ(x?+Vv;BJno5Ms0o$m@IqAUVPCe7mAAmjvtlV_8GccRVOt0**z zRdMnW=}kHWCAPios7<#;)YeT6W88@;f&j7~I!m8Y$`AhgXJ5~W*z;5LOi|{P&>!+MfzH>X_uKFM|%bbFq ze>Ef;0MEki`}s*1y#>{itSYoL^LoxmAI0dpvC(p$A^Hof@q3^mR~TnD-W;D8I4{GA zO)x(C5J_0rj`3Z-buKIHArRm@X zWD6^+R`*<~+dF+2xtO5Ls5)18)=0pYB9*ZBs>@j7mZQ1XW3quOyI{N(-QyPpTBvpJDGH373#b zA+j5sGZm(iH>gqNIp{f}T*zvx$!nL#*w~v+X3h8%!=H~cHWp9QOe($;ejB1 zz(oAleg5p%32~OKS<#OmXE#USy#BBT~v;yrlv*)Dcb%?yb_3=%+{7!V4QX|{W+P1vlNPSJ%2)# zl#~SFg*C4?)JN&NzlR3VFq{?D>z6HQlo#WY`rn#dMhM_OyP>p28kde{tG6n2!+#ZT zMKEZ+dhY+9hfu{(qAs%DdF-^fY_-$DYWIJEBj&C{QDX$+3U*$bMd&Hn%;kTdkUZsJ zTAhV}Z0_xPAKdu=W1;(l&?*Cnqll;|6kAdw>@%2#P+KMS><8oP>}1wXv(7$CdeR*N z-@o_wABtT=mlw~iVYy%^i^kdXl|EydTpE(jkOeyQ? zd}K2eH>q$2Txb!R;;MamF`3681@Hn$pJy)BGl1_Kity6@z44=-< z)=XyT8u=^LvqOwwv&T8UH@(X&CU*n#RX5bQoR43|r|NVd_oGl9rmUi{ii@kCh6{FR zv`${qily9x)=g357ko&)B9t|7<;5{krMt2mC>6_y_;g#`arQq(te>ABH0^WY&HT9> zEgnLRS|#zD+IO!uH#e_Y>7$`ECTQoPm@|s`pb4Q0^eAT;07&EEmqw747OM*gZy$DaTDl1&nO@I9rA*hPsw*A;kf{1c=?Ys*xtX4iIm zXs4NMcOf~|n(5--et5@?aT&hd{wI6Wo%Ddme3|OU=jr*sk_U4?y~<1YJW6L-i|Ml- z6>utjW8ZN-s}z|u>BOJfZm-;nlV4MhKF3ZGmJ7Li9#cGnRCvk1gKP7@HluQegnzq% z0eFrfgLFPjH4rMF4%6&Vw2E{ zHdi0A^9S25cB0=-POC1++@GyqIIv-LRnN|T)ZP8|pLqEVJX5mxmT#ZtqwHOV9k}m? zS9jl2m{nLQXG@g5k?3gf2Ow*(Z@L9nhpY{vu00^4okMO>5SBw!GwFZ%Wq(Ko^4KJE zp{3>I+yh;s{}f{10Jw+(oyX9QF}Cz!8NBdLb}ir6J_;E;#+$tM_O_&r^ZWRyTqs9+ z0vQJM%VHB=Y`bT@C~kk^vxlI^>G*$DrU z^zpW>uzbxt#YdqHObsx9vQ16qPDDW%$_Ca6bUzPJ2)~-90tyyqe?He84M_WRYuhzT z!;}UT4-Gujhxp#D&Re&jPE{)u$^04am3W3fwwBrIByW0k$J@n!T5hHM0NJnuVhADa z2AP0|Gxb7$j{akNd}SyWsS^U{!M@ty=UD{*9=(?pBmDR1iCpmXB6geX#uN4pEz_Wf z{yk&Q=6{?sy8Uq^x7+nBGS&@0;TUEV$hkl}+zPZkF@JY)RyRZJDi-xLVaKd@DD3>q zS=ik@)70cWw6<)9vX~=2lWxNoso8AkIDs#MfeIXUbap~)5=H)X-t)dk-V2tEW`%M1 zv{3?__+ee|Yx$u1@wc|oU&88wADlVq7WeBDWJ*3;!5|=fa^aCQub4U#tb_wZ!2*pQ5>U`Xofl&`}e93idnz- z6TG9y7%C0=S|pu-p+n<*WNu8;crjG3F+mq;(B7K07>bfZ4aF ziPUxZ#ZuA8b(bM~^WuU($W!d?{bwGQ_WZdw8<|Z zcxXP7jRDB)@AfZc0fy9n@FL>Dkr8tlX_21$(s>V<|J|T515~sdLK8rl7`3i%;^_F} zDS!ZvPD)wjG_$Z{eG5q$w2tfW#lg<>SBA}uqt!ldNSNg0hs4?y9UX|N zRBIi{jr0j=XF{r@{qD+gr^K-&$Qpo$nzlc-!w~URD~5HK%MC z{#%plCH+}9u@u`Vp7q`==WY;_Jjq0yyHn&uewYNq>`Uh&5@9cM0oW9wFOX@%x#aUO zSutK3`Vt-T&e#VNCj?561rrcqh<*7w;S_Zu5YCz$@}dqAX?jI~2sA?S~HM8SdffhHKf^`gE0A z_L9Ff#l$?%R_iOXEloWhd5P5%@#c#QSRE$ja14b-aCv(hPK+f3@qN^Or0d41>}GNt zp6;Kyx=C>|=M&irT>P@+S|u}-(aTW7yH_I z=!Od%Ivy>b+YD8TTZ_4BiOg^Pjh_Ck<&G&C8-Sb+#v*HPJH=&J`@bgRq!Baorc)5R zbO8UKhR(pbE}U{eqI9CR=#cX115z2FTguNglv+UVC1JOdJ}*_?OtO9crncBs0deHu zuo)ZspY2B3Amr63WAt}OTF>9*TP`EMWyA1~p%Q~@BPnM0sVmx7WDM+`mWS|PQwx$f z-SBV6#~ke{lc)Tlx#Wp?Vhne&IVrmpNAXg%pB}65tXCKahEIKi*$ycup8gNk^<)V1 ztbsC2O?0MN+%38ew_FqyM5xfASlxw%6^30A%L|{F1#TIH&lZi`yl7R}$n>3}R=K4! z8eYj`iX6gW?UFl|e1o&=+PAX$O?&mc&`Ke+>O=r?Z!wQZa5#dT4iQHZ9DkMeu%j-= zx4XT>=Fts&raZ~gG>pY(y*eA)SZ{FOdfnx2B(NDziCjD&`#)Cu+l>o;%*DlAcSr&a zdwOc6+Vj4cH#h}2`{z%!4|)TxYjhY6DqxL1xQFiljMIp_xY!<6hJpwXSdm#$$oNn2 z8YI=^-wotyYgpDglyMvhmQNhA%U$2g?=i7wle^Ny)ViqR9ha5n%h=V;L2Io+`eiKO zF{x8dDqo?_^k0UZPAyDz!dH9Nu9xPUE{!p|Et%m?GF&$j$VyA2g`>hMM7bCmWtG+e z?h|h{7A7(pg*L%=dxb5fSRPV-Bd_x z(!$#xd}4pu6tD*s5EMNc;SHiD*=?p?X!3P94$dDPSYtl=^gN<8q+F03zIHs z2euSi9n@+D;&G!C)>1NR&~L823gUl^)C?|Q*-7-vS!Bn4LU}=HvFF=#F;qnJM%>`! z%A;@bqky0Kzq{no{1`rDu0ivKFP!udf5qJJ6%NTzhXx{9*Qzz1N-TV;Cs6zXWsaXA zh?F!EMn1Sx5oVfs3csI65@s}A34nQE1l1~lZo^C%LKwEBWS&Hd#WA|kxtls;rW1P< zVWoiNoayGo;-DVoNj)d#b+um@I=}UX`y4?O<#o6p8a92i1}!by$UG=3q?F$MkXvx= zbAk0Q%yIjBvyEoD5!UV8qY5+M7tr+uuwGRj?L-cK&HCc7=;hFWy$c9nz|5n8;}PD~ z9ex1s0A{*{^`nCp1(-%GK2s<0K53cQ8}d{jlIA({!D2 zAL(Q=h0590(1Wpf21%z0|XCHXlDbHQEM<%==`OTad$^Zn!b zC7RUa#jU7_WhIgsUJ8Hte$QpyHda#Z&Q1r8yVziWIfA@)fca{{63p(iJSLzFORm79 zs*#6`X2ecNT&cLJ;7T>NJ;<`~$h4R)A}}x_qOdS8K5Kn`;Cw;%UrtS=1H$KZxP*;` zd(KkTp-s2lhVQ+*4A=bB{jPMcJj{NHTzRjX!INlg-kAle?TY_NlW`$4auLU1oOpkj zrC_eG?fjCrsB7qwqw?4pkF&-`s^eg-D_Qo&Fx0MqeSn%p1nP_tq(pOD>|Tq(N@Y08 zPy6g$`NEqc*RQOk9j?215J8hv!G>8W>rk{Z=cNMmU&A~pxV+E1*@FKhlsc&B_L&_n z;Z-e$bj}CSuW@z0cHApv`Tc`8y=%K00r#EF!nWBM?+{fM?b|imxR94+CKV_}=ii|F z7JY`8gR&7LWe(l11)*>2>U{wvcl!ebVpngLUFZi|SBKu+ma;9`^GaC}}_wb-F6_G2>6Z7XDYy}@+;;l22J4YfbJS?Sy}L9KoV4OqhP z=tL=EQZ-1U;4yyUtxeEAKAvbj#>W%V-~XI?5;k>OFtW)SLcPqwZ3m~$YzLygdWzC{ z!sV*ItX2-SZf3q^=Ym_QQ`FT7M6K@(zI|iBNEmal4+(Ixr755~V7RyIXefOAkvUpV zwr#kH;pG(q_db)FG)gEYx?B%3-!#W7ymm|l<48v%voL~Sf*-0yS)((S zJGqvVuhAw32;L7a0ZLw1ewOpWe{ z3e@t#mEh?k>&)kCpDI>4R})7h`4OH$vM9)C_9Ja;&jY3+WTf+32dC>QQ=>$1+nS?( zmE1fpDQ3UR&|&RwLS^QfV*E+Zg$iwdGFN4Xn(|+o@6od_+Xp~kxrUU*S*i%u47mf^ zZFqexU588UcuaPeYDz!6(F$NXO-T_$FmnL-)sBGGPEAYO8ZSzib`)v(8{WW`lvyyg zMfEPNNPGdDktqzm^=%`4^oPvmH+dQcm^+haqg#twzvPz1V6uz~R%X&8+jm1=efHQe z8WX0Cpl6BrN>%N!Eb*b%x(l2)1Vp^1&KG=}wEjncfSQvD;1*EanA|#^*ZA& z7NBv`E~oMotnlbF{qF{P1sXSS>N`ii>_6U_o-A8w-sC}e2P)f_r3A#fxweF*MW!D{ zc~%JhIP=(d-|p_)z_iXRU~qZ9(r?Mj%Xn`sHNMIrijq1cj%wRnbg1vV(f zY+H1qhsfDF(@!_om%^TAvw-3KC{dpFq3GIN&J&XaJ(s1!z@ZW%vTM zqN3j7H5ZKe%Kb2tiC&5|u|!_xICc0CZVn+qwHr#mp<@YPIPo^YdnFIQ&DLbi)?ka_ zn8~~_ez_x$q_x>Y_lb8XR_c)o{>f{7<->t^QR7je^V3dVfb+8x?CHY{39^)(oeSb5 zoEjfcK4|+XYj6Bu=lmMfsThZein6Uu} zo6@mUCOJi|@+Y4KUHE+|CZn;!7;*An$bVO&XGwvyA~EcqU>f1ykhiBJOiaXFKA4DB zzMtO`eOSl_pkO=Cb$7@jWUnYQG= zO?%aIg{GU%EuNu;9jJ6Sbx{!-5bdD$lmLUM5{}}7 zN#P_cmn%CNji@vAZwx~(+wokJVKjn48z=)80ho-Cjl77YQ~JGFWA)`f z>T>>c+ABvMudVhKH=FBN<_lNp&B-Z|lq))3Z;bs8kK!zPLIdz0W3qsge==c0DQH!hV-`Gd1%#(R<_K=!TG%<;_NHhk@Kf5^9f&SDpy0hB;^+7MMf!c~)q+{QR<~dC zK;rxehFJxqdq#I8|2t59?1g>wn-gwVqOs^93EkLhXYQ3o-co`UueCd)`VmHfO?9jEbwfT+5CaUCF2HoF}U;DoXhVW`_ z8t=$di~J+b-|+7Hs5itTboyK87rA@Q>0xRTh!vt50?H&^V-d%|Y04~I=l&L1aUWLa zxHlhfWH&|qJBb`>eW8kve2#?E(gPTlHINc@4r4?beENEBCd9_X611@2$;`3va<;76 z>KltKDO#8+>Ue7*9EwSY|E%ZvVo{DhsUGv{<8rs4smkIDH=YeK^tzbKv!NV<4bf}v zzke}z7%HO~eOKxNUW77?64`8OdM@z}=wUE~14RJ@N4VJ6Z}|ikrrhN+7slSvvtKg~ zwN(b6F|Voz*YqFI2TP+aovSN95MuyNdsVY0*5OEg%T@x7Bn0%q_1mdKIVp4uvoQ*| zzsQ`6=f3nfsgS!;iRkT~u^{>e*;Wjfy8=UK=Awiu{GU39aDFdl#_1!};ZGl3$OxsL zAU@ibTUkjU&0BLp<(wzKuWfC81F{VQH6^XsMLU}U_0hQRG9>k+XqO)(XR8|du6QZG zvR37z|7ugLQ-7Ld8O^4c0tj} zomiIx*M-v3+x#Ut2Z*ScF|@D#^Nbs|%-|jZmR!_3wUPp!U+BY@vVTtX;WhKkZ3n1p z!7@1zZ^Q?RU4P*zqwR}pY)UN0ua-S@G7$dqjCHw`=xi8?MAR?lJ_^pND&BalW@va9 z)j>6Y(s2CkDMP`fY1#fMG4vqIwoIz$jVda6|H9 zR_%SZ|8A`7UJ`q2yq0^f3vthbva$(DS3kj6=DU1B&z5O{7S@I zBX2xc*N}L1UKU|2^r{ThY4atscNb>Tj}Q7WPj9AKxnUY8pE2@>+OWm+BC9t@5Hc43 z+S9H|{<4&F@EiqOGJAG`KOKRK^w0N> zMXi4Qg#H>=0Zaj~CBOUXZiLRN23*mb=_G@DZ|3qcz@aSi#4Btco7tGmlfbr(a!})@ z%)g|Lx~3HsnC%~=*!jK>`|Xt_dBVzfoaRvE2vKX29eT|EA@_FzTJBv)Sc}03md9n2 zs6bZ2?m&b$>TYus&0#EC6AjR{^&e-RmE!#h-rm-z`eQq=quSgH#@ z#OfdjK{HNcX?p)d$w)5rM5m=T&Gp(!a`V0fl4Y>Z7J? zm65$TeJO;{E7u?36vx%1s!I2}UD8)*)}TB)_x_OU!P_%K^6t|6H>py3Q}N9b4-Q|$ zZ~-Z3op9(1uW(u3N4&Rd3SY`|rY>u-q_J)V&C$7SS+ zrK6Ce23ml;mR2%bhB9Ei;aL}Mi!Kuw$Q$bFMAXTcJzViZcOu)a-+1SARQaw3m0GUt z-KGQbi*w6F97KwQWirqV^lgl0LLIb>2y*d52i-5nFz?9w<+wb_F&{?J&B&s%Uz^q@ zRd{wPfSMmkQ?~x=r(J4V{Yh8KNdeLYM-{y|>G^JI(ba*9x-h24`x1drE6<>a$IZV{ zttg$r38l-ZNwhbAF04LG*xdp4Ew^Ht;e{085}0=dS{HofFi+|i znrw5&S&i+^t5U7*2Q$mT4Fg31ESx+uin;#CMGd2tJfG_1mFJk9RU&8_dA|lAoPoId zHT(~;xwim%4l(8bOAvP#nv-37nEm@kCA$XZxPkSq$Q61bcE_WDPKm{znK`<;oQJuN zhbBTbmKap{FwEsRgd9x$i1PP2<_!k)+tj>I$ZropoWU$go4A(ViYS6x^`0*|#Szj3 z$n^ZR!0i4$I%+Xf^#EG%1GPgN{Y2cV=p&_iqO0{Gw#l(|ER;I@KmdF5<@$M@J?!#k z%dDRknJ3LyupEp1oL!e;3HRJPPMo%39CcH8`uMT5w5qC#iYSffzUxtFu)G1|4Rk&? zi6{3>Gnoax&t0nP+9Y48gVkCBQ7@cqKMN!JQveBDD( zGoN~@jlbPOKzKmIVBcjduLQw++}a=8AYqYHk9g$2E}Y}?4+gOw5)NV+fSzzmN7Ees zKr!SPufH&O)Ql!Nb@kX%X!0&SPr&I(3s$Lu6-~@uPic2$ zXSwtD%XB2)IJ!W+6*6oA00qkmqwB4G!AJ4z&#d5l@0C*%x1kF=f45Cw_Hcy?k!od@MZ#_Dm=;OA;rPM6(rslV{$0F@YEcxEmc--PWEh>2&r`=U0K524{E8u2eX zHnM>R0Xp4(oNY)2W~{Ex-o7W72TfbB%z$_u10A#E8Vg=3)JRVV;zbG`gaV@$(wo0e zPB@(XO#e5pEZf~yZ%%Bn@&{#*Fm2y$9~cp|?Kmg(x;FOS%wfu%lO`Kh)+c-?Fx}h% zxEW1NO;GALZ}EW_n!T7cjpEYJ49&$Uwne8GZU{Py^BZ&obS8QWKOg0f{QgV(YsGS- z_r29q4R;sTqcpS}GRh1ox?uVY3fzApgalORt9Ndb=Mh*V^H(G_F3ipLtY5|U)1VLi z^`>Y}37h6@{%_U0{a-e8Tq|Y6{|TFDI!8kpc=^|_ z49;hb(;ywel(VocA`17*zg9{w=YQ^%y3K_^yVTKfbDSl{yu`t$Qs0T@t&!k~x9v(Y zj*<_%UK0gs@XOL>qTlaDyGfc%hYl0Ok1L$$zVoiRxFT0!OJVolj}EbVTPnt%`D-J} z+_=_DA&Hal{L)S9FGkj^*mXCmN#fNZ^2MuoZ+M;I9Lm(@g+_{*Wq&x|H2<(AZ}%TX zamfSW%xCr8;Eym=c1O0HoBq@|SfqiOO+-$vZhNYoh8k>D1th8g$XWl#=OK(mm9y`0 zag7|I2ns@z6_!f8R5Z!FXql>}c!NU7Cm~GDI>YEcmu}PZO@sn~-c#<`GXr4U`#YR0 z4sdkm64h^UB6fU8oE7VpubjHCLV2uSqo!lruuclFn^0Ykw#0ACPRrVPKCuL}MgT1EA>RRNKoU}o zt4|i5U5R=gnFEwl3cAIxcvTNXsXIV8=-~J7$3^3}XwxK!-BJA+;FTXt!I$QOBsLJZ zP@VCCcxc}JSEHBu7%g>Zir5_b(>F`$I>!HlH@9;9k&d?1!C@8Cc4UkrjU*2iZw{5rcfXdo zXn*xoyq@!Vi3|jw1Y3G_Z+g0Z;(_8D#Vuml^UGhVNMvI56eEiMa|Iq&5Vln7iy?YSiL0O$IT223* z-TW{Mx6Oa<%%;SMQX6W=2L;4Du0wDRYb>zkH#)`RV9;eby(}b%(7LIms=5+5Y%uac zINDVU50mdy<6%Ke35zxE4YWIsYlYrV=K%kQ8k1O1h(donOw=@eN}?a4A88MGgXig@ z^}{|;IfjYPFd%PQ{TYeu%FbJ#xd!8quvmw6HbPhPZqC-pM@D{WKGBLt1I#{@`cwgMz|9bW9mb6%Lt+>uoMOBxG=s zI3e`WQy+(50%xE0=Xe=jy7uMRNZh&n?B$hJqpS8n1{bJmy_Lzy6j&wv`WMG~X!CIP z+x&0N3t$V=3B0oL21lwGMm__H<6zEhfi&?3yJ7_3k8@E*8y-K$Np zf9JGRDP*VD50gN>g?|4I!Bj^OL`WX`BSM|+0xm8|NXGp!#k5o^l&g#xW&^U9Z_ya+ z6mr*+^bP2$b`&0v1y7(h;S)!xX7d3aS1wzeH#}|}%KxOevr1q3HJ|)K1n+QvCA;yq zmJb0uAM3Llzkb=ti365d1*XIlyI74UgOB<1`AWHB23n?>fg8eABnh|Nmd%tUxV#14 z7w;=JypC$G<`lwt|7JKBCP^QRtew?gB+=r2%KWiecuUP~J8&4{CFQ15RC^Zogr2HH z@_wyfTepb7rw8FOb}~D+I;em3E_v|0-lSjxszGj+#w;?C+D*zL+Z zJa3C)aA^?34M=;Ofg!QbU1us>SrRjgF8N65nS^^xF{}E)E{uAWszc z|Eev7bhk#D5sy(}X#YREUxO>zg_Bs(euKq3J=u4z1zxuNs;SG%qqH9%pyJ9PQHG4nLYjf*m#?VzRm6rc=QK7 z2+;eTTKW!~?TMT`g;l-19lBnXzOJowny*KTKsr8MR~{$I!ckf`>FOd9+tO{cff>MW zRQht}x%|3QKgUqTrWs&-;#r94eL?PbjjdVLjX>3c1V486IpdU+XR>^rNY#WOcDF!ivvob#%rYAWLys=apZZ!R^ zo5`b{$|gqjtI9X*xOrlOvQUe>Cc516b(c3^2U)YsP=B}Ccm7Ov3r&l+3tZ_KOyv9W z?yg;@QJfVg@3R({aYwp0eSYiKy!H6_`kiXl(?1bsUK*i+LtUR1rw=UEr|Br_CAk-}Kzl3a zrit#)EGZq&mfJ~72_xs{S(`TGXi=WqHiM}O#u@cw5r&!#86+kiqfv#{rqOnwMOa-> z4NbW7FT9L9(h*93Y-R9#l&d}?Mb*hx2*1+O^!rG?j;|-bd6fMHlcZ~FyZhX$b*l9* zdFzalAB)%|ZnhsIOO#Smd+OZZYcnaXil8Q-0D+Ir;Ip&`yaPToe8q~suaT`GXkwS- z#Og(-RvCdmyd0R_EN%bKuY6{jJE7CtoAeo)6uotY+5(H~U+)eqyssp;-zopO2rku< z9oOf-Rbk1HuJ0?`TEKc#@BBV^BrA-$V~Kaw zy=$k@v}p32@qEw|UH;NdcTY$t2~enp+)2eX#3!w56SB3^mn=c^2}18<#-V=iOUIkEamdV7GJ_UV_c~}C;}&=oa7bNVUk5sq$QeMi zA~yo*>_340?XVFqR9NW9#UFx86o(j2m<>~2l$+l#l&n*Hy526W6QST>TqpTa?&PO7 znhJ#&PE@7S-Udw8^-k8U_vugMqG_P%4%X-3YF6+ohzUV{dT1QmTX36Q>=u&E{hu}u zdLG@Kt77ZeS2LuA%a=>i(`(cRuMi|5Ei)z08`lFSW~oO9Qa+Z?-RWn}^JjoKfz z2_RMrB@XdHeZJQPYWU8={o&!py47$0uqJCGm>aG;({&5}bcsq2APskEZF|v3ASh@^ z-+M=4USbJ~1j{hN1qQKwZ%}rFnXEB*u}QUogkLs@W^-xq7^l3gqXsslw?H4{Y!uVL z-oh?$SejyyDOAd=V2t;9FZ;FDmhe#r;klLJ8U>W4fhx`f9xn`<562Vv*3$C=o-;OH zOmyV%ON4Ej@IX=u`AEXT4gYNa?1IdKL?!pyCIOE3^0bttH%;`QH4tm&BU9gv`O;}P z?C9?mpl)-$`MP2bULWt$f*Bwn9BvpgNr$hIVt%s zqEsb&h|l$Rzb^0_PPmzSHbp+|n@{cAUo?%#tzA0Cxm0>(^m|-*_K&Cj#sLFSUfqic zmX@+#)#;}*En#XS6&$?=~pT(zf&X%@TF%8)m zhWX0c#imn%f%NWQ|3x(h>Q6Tb{f)?T&e_`u*GlQ7Xf;*D|L`#nXIZ-0zV;FO;S;GfgS^yk0j{ z3}r{W-rpMgu~!co;x7@sC5?BO;=Pyen(WS>C_OauORV^pwz@sAhF^B@>KV=WTTOG- z%JG!~N=oCdt~C#9m^xNvJC}TR<_{W|cJTTPqUFr&N>nRxDg*#GvMbC^WP!&1iGf zuY%|m=CQU-%8R+FqQFV{0L2Q&<*(S7(0GA!ueB{3`+fM$o^7lI)iiPDZGTazOGd8p z$0k*I-HI^jBzV+r%^z?Lx6dQH_5E`*>LL}#y#n|j2&tXqxPb~mrLHd8C=s6BcWLGY zBO5+-=Po)=QtM(3ytcHwEBh*Kn(Sr5fNbl*2CJg?oF#wDV2j7h`hBbHok+Qwx!E!n zs%)L+I|sLH>;{%Tq^%^zvWclADJj)AYwESlXjY=oh9Z;2OOQV^2DGR&`4zehJhH{9 zmJ@f^E)>55TVdyX!{SvIcF|BNBa2W3q83poSiMnwEbnljJ!Dr}>MX;8iR+L*=5vW= zZ~OhezN5?EE0315Wn9t5Cb@K7AEoo!LGA_1_n+^XRc`2sjZ{cH9-!lVUl(FDRTLM3 zgI3+Kc=2!Xr$gIGM`)i^*z$Zm@wXJ;KIOM#=o(}cPa3A5pThj;;Y>8px>&zkxxeh8 zfzRiN8=YRE;@s->t}@+)`a5`?I;Emf?bbfx{x=CL$I#=c95@D&*Sz*%Y{6swTp`te>fIQ1&0xw%b&BhLVi&yx&%$)YUAB5p0oPm&#u>cN z64zhkm|@Fs(zL6U6w~|@yV7s@*3=sxF_dFo%w(uZg)bq&gELuUvzuO|uss0*T<(IW zgwZ9MtdX&Mccp@}HP%iH(u9I6`8Tv>2o|@p%^ZgsuFiZgv~JlRYG`-VKKvxzZ5XNg z<~!p`2gY*&yKIUm!O5b8OLtgwgF^ATYB$_@6>WJ5IJyN}(KlPX;vSK_WS83AF0)%J zbQ8MLx!TlkY|Y`2`dAKIv-3Qoqr!>EI5YK%{Z&bE;e8?nDdGLDSgu5#y^ECMKto=b zG#VBCC(LlFx^LYBnKhjKSV{{9!TgO9I2hRSoFZpJZ35!x?_Y-%LP3Gy@W2F*o?dO% zrD$Vl(pRq_cs5z=QT+hc(zGG2>YGzx#+=s#Q-rS)GL&*+pjYYo#7)hHhC^(>jTau@RW{k#3$|5MpnhD8~*378Z? z1*AbjYKD{&kVaBMx;v!18>FSXBqgM~yFt1;B?Y8Jy7tT$`)9xF+Fk#Q&d9vH=RJL% z=f1g;D+4d-Jy@B^W5`ytLnim0{usY(CW2 zO{0PpEG;)l-+m$^q4feG1Q-@ZN#SV(o6sSxS^9Uc88=bP)Ncy5lj=ujZLwxox|_hj z`p)94{NZX7gI~RwBnqz;7$q9Tu5{G~R(PmA5BDUpRitXerQLPf6%1IyCar8+J6o1+ zwBL3RVQ^-km;4Zt#qhRs(*Kr^!iPv-HM$ZwyGn~pDs_3RYIi`DgI4ze6r8qo=fqV0 z?>dWGDCuVsC)|gklAqb_U7|jLGgj;+W=3vj4L7q(xg>lVC8uceL=)hItQL$nj_$i(>?hEs|y}c#H{)CX!oQCJNFxSD`(FT&Nv~4l#!EifH2mX3 zRhc$oYlemIjC$oMj)2vze=c74I9Idc(pUoO3PZ0aE&lT&$p)EzA{SMTrl#m`@&z9B z>HU+FjAPI2f(h2lO*f*VZ4T?QvkK^nDeL9L@p~#pHw$9yQbxb}U(woWDYnjkvK}M^ z{JnLR6F1`#D%%Uo8dB0;dy>2qFsGj9nyG6+oj>Rd)vZK`SO|}keR9USsh+Nvf^92P zzFRWU!H*9$uRZTEEdt%bQ~(C4G_Wer$kgv9e5e#TeLaxklYpg=ew3+Z%`a_i<48&* zNL?a3YO< z*+3s1^T1p%#GP+4nED839Z&DtO`(+P1JMik_fc7W)>#=l3K#z7N9{#uQ%+|oapR;` z_fHZ7rB$4N{$NfGVK;h-i7R>)#-su1PW_X4{QiE4(V;@5a?0Ix>tSP=+Y7RGCIQGI zD$>u1R4L1R_uMkh*is0hqmepT>s_VYt=h`R^sW3A$OntkkPdQXaf`?3L%RmOcJ+GU zW(h`Si(Tk}3`75XVWyX`h{pdn9+mRQw9Q+%Jbm3_Mx5TsAtT}Dv_0jt6ng-N5(i>vzDk331K$2xl6M)N@lhmxP)^(KMK(iZ+r0# zr4PE1YD6)ZHrOZiDvVEsoS5-!k)4%Ex37dI?WD-EHlnDghWGxZ7N2Bp96cnpIj6w+ zG9McQYJn`lK6n_vUYc%GDZsZAsWa8`77dAZPThkl^u2(nZ6D=lu9q5pBtv>K0g3aH zGs(@4wB@b$9Gf{gDIWLk651MtBeW>C?%D8Pu_9ZN#y=gv|Ho-w^A7Him7qp-GGA!tV=wvRfxVubvZs z850QBlacW*fGCBci0{Koxdky!sFvc{l@w|H^g2u;<{)xqWDYpsh^4Rit9M?u+T0+7 zYm`>&mz(K)6K}WFS4nN~%y3W>Pi9LxSoJe<2B~yptAtwh(OtbaW00p|BRTV~%v9B# z$z74*Ta4510IWOx;7v=5oZ#r%0XkaA07V%;ry~%S014HJ;-U@=sI`DUasL^(6ya4Q zl<(popj=cV+-#3Yw%ju8*n4qw3 zCQ4bbhHTbF7*1_}?9jGAR!)#;TYam}7C_%QH9_|#`6y$_HxFvD;QJKp_KRMlZPPY5 z?)i&?X&l90Utd!#E$001vEHW?GVJEj72b$DTc7?W_}f`cEneWy9lVISHLJc zo?x#OzIU>EoIRHwY&-Rd9(Q4r%We@K?P<5xQ>kv)fRG#{52Rhlz8kRApcOk2-rbI= z&X*F7GJTc%yN#YT(XT)Hmx8U0a7S9V%SLUIEgv>8LQ(EVjL$RpsYIXgB-4hMNg@;F|?S9_aoob8$-UAyg*%Y~c)dd?WfWZ8xO zZHh|CJIwdT5MQO0E4Cyi9iyqtiJs*Txl>h9r;<_zib4oBU`M+^K8JwQns{|}^?Fnl zZ-7on-oUz9#@`X{(3_B(oa5%T&C)cwI_hz2GwJ5__OY>`g1OV* zY(%#!yelpv~hi+}W%Yh}yfUHj5TYhKd-Pq~ESoEL6vz959JO?>U zxcouYZ*@b3cbX~7JIlLDS8Fb}NoaXWBQmkMEPbihtO12^Sc%-h5Oc-Y%!PbrB%^ne;fz z{!~8R`tX|QEDbW}Mav&a z1}Y@5LQa>c$AN%1W*~P)Ktyx{k-S`h-kxD-A|@K+{O40!l!u!Dt#zommF+os0!Ycv z&xZo+m@99nh3jG_JG05FjB?R0*|wj~==N>Mq;qwWGz|+n-4ebFaK-%8-0@@@P|1raIHLXCW_4xKa z_w{8e_75r-d5Wryjgkk}j@Hus<`lQ79g#@Wz=+~rCZQn%2I(Fa0g$m^c`=ZtIlhd@&DX&=k zshXmgZECjHD$g_6L07 z0GWHyepLtU#23f?VCUwS_URM;k+ePaRp{80ky;7{8$KW8F3eVeZjVxNA4@JTAEI7! zPR?vc4TLJ`zcMnAtAR%i^FGv{f(faQsRPfFSO@xt_CE^d(8*T4cIUB3J1?l)F)ySG z(>y5jS|B+-+qLs3hK?NMZ>cyQ5M)o<{8Xw-lhHrj=xs=fv2ywJF=`wAjbXnV!;W~m9)beU{2jdim|xRhb9JDYUDl#CTxvibb0mUrmPtqFF5 zb2{AapTE}pU)S}9nAg6nN=e+5QLdTj46RX@(j+Sak`A^__fjicG}Q>!zev9!b!;);IS%Cu~}8J9T;^8@zsj^kjBm zDoj-_HO)DXWg`5YOWF^{GZre**q@{y^n&$~l#4UVS5`OMtIPrJjBTYy>fmE0UXH2q zBkE8b=V-O@kn3vcR{!CGaz)2f#DM1NYX<5-;ZS@l=}{F!=Hw1GT%zg25B8okx3O(? z40ftnquWcUiG-%V&qXRd$pUPVJRE9Q8()eymXTHc>Q@&m8Zngikr8$@;(Oqn8W=Zl z8mf8D*kE1?HKQ(7f~vtX`9RI`Zr*^H?gA*k47=nP8 zzrfK1^5=jX1?2k%13y>>VwvF9I0}-S)8-iUMc*sTajb2vfSFw&=G1j0YEysG-2Y=5 zAtvUNnAoCK(n)s|%@eP0DiX$GqIRy2W0qXkd?DYFcxy5o3znk#Mw+*!6wszdty$r0 zU0IS!{?-NXaxJlPQrv2pG&#PHOA$@w*=XxEAN}PteYK7M8DVKlOIP3H7Okl*`%dRj zlMQ7JTXBB#Xiyo-h^vV4^5XvTcQI!OI4ayBJ(yw&I+0plT^6{jq++$|CTKb~@2Y@` z+gcCdte&3^vaziEp(t)qTumhdC#g}K<9+*rUqnBb5jk*=?ig62X}5*;ueYOlL6nBx ztU0i%35z+BJuGxzMC4)`Man97aeIEi)^w=dzcA?^ebJ$x8eFiUnndlNSbZ@>9!^+? zU;p;cc@Q;EH_h&1uQykqd=a#<9+Wj8gs|DH^1^{bR1j_wPRtbxl76==*6M9^f$vNW zXjP63f$T9GD3W212RTC)(kSq*3%h%-aff!eh|A~xwP<&!HN84pxd2}hbcwfbiOKv# zIuKoa?}NN4xwE;CjRXCAQ0|&mCLd%_`>c;PJ-4?~%16dEQU4M&zVIxvw_CGWb~pRf=CI8?>WadvCE;V#GQuIV zHk@aF&yh>aD+2Td97};E8X==Y#T^HTXy<%1!i1U@lK zz0zBGzmfLhE~lU1yTD3}Nhck1`3jBdC&AehD!#zF}d{01}dlfq~)j6dv&jCK=fn&}2s+ zTlo3=69PjTNUj_8vIJ8`uu01yTgj?P)#+YQ!dIUt(Y!kF(5qPO5*R#Sm5-^y7{2)| zNdLh!c!YDBo#}gZj^E2UwbB*N;hb+-?82=lw_z^24zkt>^ugp@fhr!B7`oOzWoRjH z(4yXO+pt|2GP!PAhf+GGm{*UKGyRgE^eVKy;~sfzY7bq92x@+DmMG)4<}#!4>2|jv z>zL@YFynC;S1OHak67y#Ljrpp5EfLn>7O=7b~joVu*fy zQIs?_*~3fN80i2PbEHMmU^AjMn3+5bw%5w`Yn6uIPSDc^9=M!Ku0^3ZgANoCl}+oVLFr{C2JY8nr%PQRSS<#Ss6pKc8dz)7mt z;SZzEsb{$FC+L3J6e2C8_~^2%9u>m5x`Z@#Y;`k|Oon61zckmXS~&_{Ri`oj)!cQF z8D^vpSsv_ceEuRQai0IwMZ0qw?=N^fp8ZAY4_3C|3mr)ass&XAA8$W1(s3kmZxP*C zbXvh>G^`6DW_mk*73em$f0jZZBs}mXJYe`vz0iPpUBHd?lwLK~`6O$^A>q6Ev`{>b zmmnOw_jTZpAMNw@ZbKwtX29&F#_4&W%^9T9Gd~KQfo9A0Wf6f>^SCJX}tLgMd{%Fd8*xLLm?=9I7)?Qo4v+MXm;1!`~%U2a9{`eWrc$X>wex z-Qw@a)MSEYFKjp0v&7YX&Jk>TjcK6Vss>%D?JW2HGDuz)%_~PT`~=o1rgdzqM|rsM zo_EVR9ERVSb{GFX*zP5dIF~wol*5HvveIbe_k9$eNi@1DkQT9zhw2MuCm4q zRyGk#N(RdouabEh>19jmR)t?+?iTka7xq&BOdR2mXe+x?E)^Hj;Jq=-x2m(K;vsb& ztGMRs(W?hKe9`h|YQOiz{O~T@!xk)tPgVh5F3gH*10B zLn3EVKSH+OH7k(q`^x0uQCDftHRIGe?cl6mj?;T5OtI}g$Pu(YV{z)m3i+Gua^E^Q z;0`tKWG5A#(%oJ;hs6PX>3DSGz}NisP}xSvSe-I>h*jT*Z=qE0gZsg zlLSpi2HSdg;ac-gY1LDuPmsTn-e<$8EUf8mv;iS=^<&yZmN7l(KGmtsQSNIF+(s*yvlX7`RZoau6lrtcwGMmBtc zt~S*&^y2)kHIrjO`aEd0SpE2@({(#f?CoMlp_1RrIV5sBt+{M@Lg%zpgF>T}nn;qd zJEiKgnbg2HtIFx@I$HIo9CIx=_D1KrS;M{ISY@w%N3h8?1 zhf7HyLFH_e=NOZO;F{z1g{R zx&uozTl8M09-Gbrjlg;u%(TGxO}z5>tUdkiO}(*Ve$&0#e2bKIl(_y9?sr=Y{jg>Z zvcHg+U3Z1cNYtD@bf&30@ET}_cygT7o}H5ZsitC1b#kq^Yudc}A`MeWP0@Vl?K{7C zGfX21eY)>Svqk&Vr3GymqP1pjGSYc}kos3hj11idy&pgqlcJ$?x)`GK_{}krB zZ~}V^gG5CDW_ENaH*GAE+82~0Dm9^}b{JH&SVA~Kt=U6DPeBIEhD|h*6znhAU%<9Y z0H)AW`5thP4f~A7AQrgWz*{&hurC5{g(kzk0=y+Ey|1@IZZekza9$|l6tEA)2M3TJ zjw1bH!3VUosGk9y_n)79nS2fYEh{_b@xZZOBXC{kN-dJKr(v56t>eu5FsWWOmXDfVr`Vv;%pT6&)50^Di_$j zexZQy{reuV82!-@zP+4@&j z($R3~V3TlQrX?XsY!C;kGjVwwu_)5RUveVkT8;LB_`04iL4*sS0~G4E^MfD_TSv7n z^tgc@z|TQ_@uGKOA*D#U_zw_h-XAg=lj6kcC42_vVHQBG1LvqcV0LdhUqKFN;^YuW z-^4dm9uQ94ck;7)q1lBA06Zz!*y5lWu72_HIN~6E+zg;!JV9;@P>6R-eET8<6oNS* z)d>jSVK!P=P9)Igg(M{<BtbICU$~?=ncc? zVP@?IPu?IgCZ-eU{PqC0|HzOiA_!>uDow^nQN$oe%0fa2aA1cPJyHY*RuaGl;*qh& zC*O1r!oc2I%!Zp}dL@GAG7cIx$Rd{#52i*N2dx_P=)*w14hJ(;@wgo20)_&54bx+TtRf!R zT>wvY2e@1xOe+tf)qxyC0Os?3;_nM4v+uZh^}9ia?+0n=T#z^dFyzt@Y(VAb2S4rN zXo^k^YFZ#ONOp(iCxhfx5|U6FK@kxon7RSAmgRv0}eHngN~vX^M)0VFJcR3@1NohliU9>V@A~UnEiXI_cnH z@;@MVf>i@NMh#|=@Bx-@1aD2FVHFHE16-W>Y~J{8yg=&^{8R!Lj%}cF!2sAS=NhTK|gur^pWh0;o0mVE#q5qV8 zka-~tSV*gkK7Cr>-L9oq4m622=F_#}` z|G@qgi7xe2{>eYaqXi($1bAkc_W{_a#e*pWXu3b4!>;xdH1t?4RPv+Le~fFvz##@s z{{k9pU^m`F0;-J+Xg0(oL73U$aj-G&NRK1+U$kc`Dynzf#NZc)H8`PVeGsDslh;Tl z(6N9%>K|9h$%S3zAH6>5{~QPYQ ZNE with Qiskit: Layerwise folding ZNE with Qiskit: Quantum simulation of quantum many body scars +ZNE with Qiskit: Simulation of Loschmidt echo revival ZNE with Cirq: Solving MaxCut with QAOA ZNE with Cirq: Hamiltonian simulation with Pauli gates ZNE with Cirq: Energy of molecular Hydrogen diff --git a/docs/source/examples/loschmidt_echo_revival_zne.md b/docs/source/examples/loschmidt_echo_revival_zne.md new file mode 100644 index 0000000000..1e00b5ad5b --- /dev/null +++ b/docs/source/examples/loschmidt_echo_revival_zne.md @@ -0,0 +1,425 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.16.4 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +```{tags} zne, qiskit, intermediate +``` + +# ZNE with Qiskit: Simulation of Loschmidt Echo Revival + +This tutorial replicates some of the results from {cite}`Javanmard_2022_arxiv`. We build a circuit that simulates the time-evolution of a transverse-field Ising model, then run ideal, noisy, and noise-mitigated simulations of the circuit. + +In particular, let $\Lambda(t)$ be the _Loschmidt echo_, the probability that the system returns to its initial state at time $t$. $\Lambda(t)$ has a quasi-periodic series of peaks that are flattened as the noise level increases. Here, we demonstrate how to simulate the Loschmidt echo and use zero-noise extrapolation to mitigate the effects of noise. + +The paper considers some additional effects of noise, which are outside the scope of this tutorial: + +* Let $\lambda(t) = \lim_{N\to\infty} -\log(\Lambda(t))/N$, where $N$ is the number of sites in the Ising model. Dynamical quantum phase transitions (DQPTs) occur at values of $t$ where $\lambda(t)$ is not analytic. DQPTs are observed at different times in the ideal and noisy simulations, and occur more frequently in the noisy system. + +* Noise weakens the correlations between adjacent sites. + ++++ + +## Model definition + +The Ising model that we will simulate has the Hamiltonian + +$$H = H_{zz} + H_{xx} + H_{x}$$ + +where $H_{zz}$ and $H_{xx}$ are the interactions between neighboring sites and $H_x$ is the interaction with the external magnetic field. Specifically, for $N$ sites, + +$$H_{zz} = -\frac{1}{2} \left[ \sum_{i=0}^{N-2}J_z Z_i Z_{i+1} \right], \hspace{0.4cm} H_{xx} = -\frac{1}{2} \left[ \sum_{i=0}^{N-2}J_x X_{i} X_{i+1} \right], \hspace{0.4cm} H_x = -\frac{1}{2} \left[ \sum_{i=0}^{N-1} h_x X_i \right]$$ + +where $X_i$ and $Z_i$ are the Pauli operators acting on site $i$, $J_z$ and $J_x$ are the $z$- and $x$-components of the spin-spin coupling, and $h_x$ is the strength of the external field. Here we will set $J_z > 0$ and $J_x > 0$, so that the spins at adjacent sites are correlated, and set $h_x > 0$ so that each spin prefers to have a positive $x$-component. (Strictly speaking, since $J_x \neq 0$ this is a Heisenberg model rather than an Ising model.) + +Assuming the system is in state $\ket{\psi_0}$ at $t = 0$, we want to compute the Loschmidt echo, + +$$\Lambda(t) = \left|\bra{\psi_0}U(t)\ket{\psi_0}\right|^2,$$ + +where $U(t) = \exp(-iHt)$ is the time-evolution operator. + ++++ + +## Reformulation as a quantum circuit + +To simulate how the model behaves over $0 \leq t \leq t_{\textrm{max}}$, we divide the interval into $M$ steps. Letting $\delta t = t_{\textrm{max}}/M$, we have + +$$U(k\delta t) = [\exp(-iH\delta t)]^k \hspace{0.25cm} (k = 0, \ldots, M)$$ + +Next we decompose $\exp(-iH\delta t)$. Up to an $\mathcal{O}(\delta t^2)$ error (since the terms in $H$ do not commute), + +$$\exp(-i \left[H_{zz} + H_{xx} + H_{x}\right] \delta t) \approx \exp(-i H_{zz} \delta t)\exp(-i H_{xx} \delta t)\exp(-i H_{x} \delta t)$$ + +Now we observe that each term in the decomposition corresponds to a series of gates in an $N$-qubit circuit. For example, + +$$\exp(-i H_{zz} \delta t) = \prod_{i=0}^{N-2} \exp\left( -i\frac{J_z\delta t }{2} Z_i Z_{i+1} \right)$$ + +Using the fact that $Z_i Z_{i+1} = I \otimes \cdots Z \otimes Z \cdots \otimes I$, we can rewrite this as a product of [$R_{ZZ}$ gates](https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.RZZGate), + +$$\prod_{i=0}^{N-2} R_{ZZ}^{(i, i+1)}(J_z \delta t )$$ + +Similarly, the terms $\exp(-i H_{xx} \delta t)$ and $\exp(-i H_{x} \delta t)$ can be rewritten in terms of [$R_{XX}$](https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.RXXGate) and [$R_X$](https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.RXGate) gates, yielding + +$$\exp(-iH\delta t) \approx \prod_{i=0}^{N-2} R_{ZZ}^{(i, i+1)}(J_z \delta t) \prod_{i=0}^{N-2} R_{XX}^{(i, i+1)}(J_x \delta t) \prod_{i=0}^{N-1} R_{X}^{i}(h_x \delta t)$$ + +To compute $\Lambda(k\delta t)$, we repeat this sequence of gates $k$ times. The circuit is implemented by the function in the following cell. Note that: +* We use $\ket{\psi_0} = H^{\otimes N}\ket{0^{\otimes N}}$, i.e. the spin at every site starts out parallel to the external magnetic field. +* The $R_{ZZ}$ gates commute with each other, so we can organize them into two "layers", one layer for the pairs of adjacent qubits $(i, i+1)$ with $i$ even and another for the pairs with $i$ odd. The $R_{XX}$ gates are organized in a similar way. + +```{code-cell} ipython3 +from qiskit import QuantumCircuit + +def create_circuit(delta_t : float, + k : int, + n_qubits: int = 6, + measure : bool = True, + J_z : float = 1.0, + J_x : float = 0.1, + h_x : float = 0.1) -> QuantumCircuit: + theta = J_z * delta_t + phi = J_x * delta_t + gamma = h_x * delta_t + + circuit = QuantumCircuit(n_qubits) + + circuit.h(circuit.qubits) + + for _ in range(k): + for ii in range(0, n_qubits-1, 2): + circuit.rzz(theta, ii, ii+1) + for ii in range(1, n_qubits-1, 2): + circuit.rzz(theta, ii, ii+1) + + for ii in range(0, n_qubits-1, 2): + circuit.rxx(phi, ii, ii+1) + for ii in range(1, n_qubits-1, 2): + circuit.rxx(phi, ii, ii+1) + + circuit.rx(gamma, circuit.qubits) + + circuit.h(circuit.qubits) + + if measure: + circuit.measure_all() + + return circuit +``` + +```{code-cell} ipython3 +print(create_circuit(0.01, 1)) +``` + +## Ideal simulation + +To get a sense of the ideal behavior of the circuit, we will run a noiseless state-vector simulation. First we define a dataclass to hold all the simulation parameters, and functions to visualize the results. + +```{code-cell} ipython3 +import numpy as np +import dataclasses + +@dataclasses.dataclass +class SimulationParams: + '''Simulation parameters. Note that by default we use + much coarser time steps than in the paper, so that + the noisy simulations later in this demo run in + a reasonably short time.''' + t_max : float = 8.5 + M : int = 25 + n_qubits : int = 6 + dt : float = dataclasses.field(init=False) + t_vals : np.ndarray[float] = dataclasses.field(init=False) + + # Only used in noisy simulations + n_shots : int = 2048 + + def __post_init__(self): + self.dt = self.t_max / self.M + self.t_vals = np.linspace(0., self.t_max, self.M + 1, endpoint=True) +``` + +```{code-cell} ipython3 +from matplotlib import pyplot as plt + +def setup_plot(title : str = None): + plt.figure(figsize=(6.0, 4.0)) + plt.xlabel("$t$") + plt.ylabel("$\\Lambda(t)$") + if title is not None: + plt.title(title) + +def add_to_plot(x : np.ndarray[float], + y : np.ndarray[float], + label : str, + legend : list[str]): + if label == "ideal": + plt.plot(x, y, color="black") + elif label == "mitigated": + plt.plot(x, y, marker='s', markersize=5) + else: + plt.plot(x, y, marker='.', markersize=10) + legend.append(label) +``` + +To get the ideal result for $\Lambda(k\delta t)$, we omit the measurements from the circuit, and read the probability amplitude $\bra{0^{\otimes N}} H^{\otimes N} U(k \delta t) H^{\otimes N} \ket{0^{\otimes N}}$ directly from the final state vector. + +```{code-cell} ipython3 +from qiskit_aer import QasmSimulator + +def simulate_ideal(circuit: QuantumCircuit) -> float: + simulator = QasmSimulator(method="statevector", noise_model=None) + + circuit.save_statevector() + job = simulator.run(circuit, shots=1) + + psi = job.result().data()["statevector"] + + # Get the probability of returning to |00...0> + amp_0 = np.asarray(psi)[0] + return np.abs(amp_0.real**2 + amp_0.imag**2) + +def run_ideal(params : SimulationParams = SimulationParams()) -> tuple[np.ndarray[float]]: + echo = np.array([ + simulate_ideal( + create_circuit(params.dt, k, n_qubits=params.n_qubits, measure=False) + ) + for k in range(0, params.M + 1) + ]) + + return params.t_vals, echo +``` + +```{code-cell} ipython3 +# Run the state-vector simulation +result_ideal = run_ideal(SimulationParams(t_max=25.0, M=200)) +``` + +```{code-cell} ipython3 +setup_plot("State-vector simulation") +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +plt.legend(legend) +plt.show() +``` + +The result shows a series of Loschmidt echo peaks. Intuitively, for a noisy system we expect that as $t$ increases $\Lambda(t)$ will approach $1/2^N$, where $N$ is the number of qubits. This means that as the noise level increases, the peaks will be suppressed, starting at larger values of $t$. At higher levels of noise, we may not be able to detect the Loschmidt echo at all. + +In the rest of this notebook, we will run simulations with two different noise models, and try to reconstruct the peak at $t \approx 6.5$ with zero-noise extrapolation. The next cell re-plots the ideal simulation result over the values of $t$ that we will consider. + +```{code-cell} ipython3 +result_ideal = run_ideal(SimulationParams(M=100)) +setup_plot("State-vector simulation") +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +plt.legend(legend) +plt.show() +``` + +## Simulation with depolarizing noise + +The next few cells run a simulation with depolarizing noise. We transpile the circuit using 1-qubit rotations and CNOT as the basis gate set, and optionally use gate folding to scale the noise. + +```{code-cell} ipython3 +from qiskit.compiler import transpile +from qiskit_aer.noise import NoiseModel +from qiskit_aer.noise.errors.standard_errors import depolarizing_error +from mitiq.zne.scaling import fold_gates_at_random + +def simulate_noisy(circuit: QuantumCircuit, + noise_model: NoiseModel, + n_shots : int, + scale_factor: float = None) -> float: + # Transpile the circuit + backend = QasmSimulator(noise_model=noise_model) + exec_circuit = transpile( + circuit, + backend=backend, + basis_gates=["u1", "u2", "u3", "cx"], + optimization_level=0 + ) + + # Apply gate folding + folded_circuit = exec_circuit if scale_factor is None \ + else fold_gates_at_random(exec_circuit, scale_factor) + + job = backend.run(folded_circuit, shots=n_shots) + + # Get the probability of returning to |00...0> + counts = job.result().get_counts() + ket = "0" * circuit.num_qubits + if ket in counts: + return counts[ket]/n_shots + return 0.0 + +def run_depolarizing_noise(params : SimulationParams = SimulationParams(), + noise_level : float = 0.001, + scale_factor : float = None) -> tuple[np.ndarray[float]]: + basis_gates = ["u1", "u2", "u3", "cx"] + noise_model = NoiseModel(basis_gates) + # Add depolarizing noise to the 1-qubit gates in the basis set + noise_model.add_all_qubit_quantum_error( + depolarizing_error(noise_level, 1), basis_gates[:3] + ) + + echo = np.array([ + simulate_noisy( + create_circuit(params.dt, k, n_qubits=params.n_qubits, measure=True), + noise_model, + params.n_shots, + scale_factor + ) + for k in range(0, params.M + 1) + ]) + + return params.t_vals, echo +``` + +```{code-cell} ipython3 +low_noise = 0.0005 +result_depolarizing = run_depolarizing_noise(noise_level=low_noise) +``` + +```{code-cell} ipython3 +setup_plot("Ideal and noisy simulations") +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +add_to_plot( + *result_depolarizing, + "depolarizing noise, $p = {}$".format(low_noise), + legend) +plt.legend(legend) +plt.show() +``` + +As expected, the Loschmidt echo revival is weaker than in the ideal cases, and we get a lower peak. Applying gate folding suppresses the peak further; in the next cell we do this for scale factors $\alpha = 1, 2, 3$. + +```{code-cell} ipython3 +scale_factors = [1, 2, 3] +result_depolarizing_scaled = [ + run_depolarizing_noise( + noise_level=low_noise, + scale_factor=alpha) + for alpha in scale_factors +] +``` + +```{code-cell} ipython3 +setup_plot("Scaled noisy simulations") +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +for alpha, result in zip(scale_factors, result_depolarizing_scaled): + add_to_plot(*result, "$\\alpha = {}$".format(alpha), legend) +plt.legend(legend) +plt.show() +``` + +## Error mitigation with zero-noise extrapolation +At this level of noise, we can use ZNE to mostly recover the ideal result. Running the noisy simulation is expensive (especially as the scale factor $\alpha$ increases), so rather than using the high-level functions for ZNE, we apply the static method `RichardsonFactory.extrapolate` to the results from the previous cell. + +```{code-cell} ipython3 +from mitiq.zne.inference import RichardsonFactory + +result_zne = RichardsonFactory.extrapolate( + scale_factors, + [r[1] for r in result_depolarizing_scaled] +) +``` + +```{code-cell} ipython3 +setup_plot("ZNE for depolarizing noise model, $p = {}$".format(low_noise)) +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +for alpha, result in zip(scale_factors, result_depolarizing_scaled): + add_to_plot(*result, "$\\alpha = {}$".format(alpha), legend) +add_to_plot(result_depolarizing_scaled[0][0], result_zne, "mitigated", legend) +plt.legend(legend) +plt.show() +``` + +Increasing the baseline noise level makes it much harder to reconstruct the peak with ZNE. + +```{code-cell} ipython3 +high_noise = 0.005 +scale_factors = [1, 2, 3] +result_depolarizing_scaled = [ + run_depolarizing_noise( + noise_level=high_noise, + scale_factor=alpha) + for alpha in scale_factors +] +result_zne = RichardsonFactory.extrapolate( + scale_factors, + [r[1] for r in result_depolarizing_scaled] +) +``` + +```{code-cell} ipython3 +setup_plot("ZNE for depolarizing noise model, $p = {}$".format(high_noise)) +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +for alpha, result in zip(scale_factors, result_depolarizing_scaled): + add_to_plot(*result, "$\\alpha = {}$".format(alpha), legend) +add_to_plot(result_depolarizing_scaled[0][0], result_zne, "mitigated", legend) +plt.legend(legend) +plt.show() +``` + +## Simulation with realistic device noise + +Next, we run simulations using the noise model of the IBM Nairobi device. Again, the relatively high noise level makes it difficult to recover something close to the ideal signal. + +```{code-cell} ipython3 +from qiskit_ibm_runtime.fake_provider.backends import FakeNairobiV2 + +def run_ibm_nairobi_noise(params : SimulationParams = SimulationParams(), + scale_factor : float = None) -> tuple[np.ndarray[float]]: + noise_model = NoiseModel.from_backend(FakeNairobiV2()) + + echo = np.array([ + simulate_noisy( + create_circuit(params.dt, k, n_qubits=params.n_qubits, measure=True), + noise_model, + params.n_shots, + scale_factor + ) + for k in range(0, params.M + 1) + ]) + + return params.t_vals, echo +``` + +```{code-cell} ipython3 +scale_factors = [1, 2, 3] +result_ibm_nairobi_scaled = [ + run_ibm_nairobi_noise(scale_factor=alpha) for alpha in scale_factors +] +``` + +```{code-cell} ipython3 +result_zne = RichardsonFactory.extrapolate( + scale_factors, + [r[1] for r in result_ibm_nairobi_scaled] +) +``` + +```{code-cell} ipython3 +setup_plot("ZNE for IBM Nairobi noise model") +legend = [] +add_to_plot(*result_ideal, "ideal", legend) +for alpha, result in zip(scale_factors, result_ibm_nairobi_scaled): + add_to_plot(*result, "$\\alpha = {}$".format(alpha), legend) +add_to_plot(result_depolarizing_scaled[0][0], result_zne, "mitigated", legend) +plt.legend(legend) +plt.show() +``` + +## Summary + +For the system considered in this tutorial, the effectiveness of ZNE depends strongly on the level of noise. Under low-noise conditions, we can accurately recover the ideal Loschmidt echo signal. At higher noise levels, increasing the scale factor $\alpha$ almost completely flattens the first revival peak. As a result, the extrapolation procedure can qualitatively recover some of the signal, but is not quantitatively reliable. diff --git a/docs/source/refs.bib b/docs/source/refs.bib index f7b6762353..c96b3fb3d6 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -395,6 +395,15 @@ @book{James_2021_statlearning url = {https://www.statlearning.com/} } +@misc{Javanmard_2022_arxiv, + author = {Javanmard, Younes and Liaubaite, Ugne and {Osborne}, {Tobias J.} and Santos, {Luis}}, + title = {Quantum simulation of dynamical phase transitions in noisy quantum devices}, + year = {2022}, + month = {nov}, + archiveprefix = {arXiv}, + eprint = {2211.08318}, + primaryclass = {quant-ph}, +} @misc{Jurcevic_2021_arxiv, author = {{Jurcevic}, Petar and {Javadi-Abhari}, Ali and {Bishop}, Lev S. and {Lauer}, Isaac and {Bogorin}, Daniela F. and {Brink}, Markus and {Capelluto}, Lauren and {G{\"u}nl{\"u}k}, Oktay and {Itoko}, Toshinari and {Kanazawa}, Naoki and {Kandala}, Abhinav and {Keefe}, George A. and {Krsulich}, Kevin and {Landers}, William and {Lewandowski}, Eric P. and {McClure}, Douglas T. and {Nannicini}, Giacomo and {Narasgond}, Adinath and {Nayfeh}, Hasan M. and {Pritchett}, Emily and {Rothwell}, Mary Beth and {Srinivasan}, Srikanth and {Sundaresan}, Neereja and {Wang}, Cindy and {Wei}, Ken X. and {Wood}, Christopher J. and {Yau}, Jeng-Bang and {Zhang}, Eric J. and {Dial}, Oliver E. and {Chow}, Jerry M. and {Gambetta}, Jay M.},