From 43682dd4ca7f0dbd7ef11dd6af2560c4ec3f325b Mon Sep 17 00:00:00 2001 From: Dario Coscia <93731561+dario-coscia@users.noreply.github.com> Date: Tue, 6 Dec 2022 11:40:42 +0100 Subject: [PATCH] adding new tutorial (#42) * adding new tutorial * Update tutorial.rst --- docs/source/_rst/tutorial3/tutorial.rst | 205 ++++++++++++ .../tutorial_files/tutorial_12_0.png | Bin 0 -> 20608 bytes .../tutorial_files/tutorial_14_0.png | Bin 0 -> 17833 bytes docs/source/index.rst | 1 + tutorials/README.md | 1 + tutorials/tutorial3/tutorial.ipynb | 302 ++++++++++++++++++ tutorials/tutorial3/tutorial.py | 153 +++++++++ 7 files changed, 662 insertions(+) create mode 100644 docs/source/_rst/tutorial3/tutorial.rst create mode 100644 docs/source/_rst/tutorial3/tutorial_files/tutorial_12_0.png create mode 100644 docs/source/_rst/tutorial3/tutorial_files/tutorial_14_0.png create mode 100644 tutorials/tutorial3/tutorial.ipynb create mode 100644 tutorials/tutorial3/tutorial.py diff --git a/docs/source/_rst/tutorial3/tutorial.rst b/docs/source/_rst/tutorial3/tutorial.rst new file mode 100644 index 0000000..f76346a --- /dev/null +++ b/docs/source/_rst/tutorial3/tutorial.rst @@ -0,0 +1,205 @@ +Tutorial 3: resolution of wave equation with custom Network +=========================================================== + +The problem solution +~~~~~~~~~~~~~~~~~~~~ + +In this tutorial we present how to solve the wave equation using the +``SpatialProblem`` and ``TimeDependentProblem`` class, and the +``Network`` class for building custom **torch** networks. + +The problem is written in the following form: + +:raw-latex:`\begin{equation} +\begin{cases} +\Delta u(x,y,t) = \frac{\partial^2}{\partial t^2} u(x,y,t) \quad \text{in } D, \\\\ +u(x, y, t=0) = \sin(\pi x)\sin(\pi y)\cos(\sqrt{2}\pi), \\\\ +u(x, y, t) = 0 \quad \text{on } \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4, +\end{cases} +\end{equation}` + +where :math:`D` is a square domain :math:`[0,1]^2`, and +:math:`\Gamma_i`, with :math:`i=1,...,4`, are the boundaries of the +square, and the velocity in the standard wave equation is fixed to one. + +First of all, some useful imports. + +.. code:: ipython3 + + import torch + + from pina.problem import SpatialProblem, TimeDependentProblem + from pina.operators import nabla, grad + from pina.model import Network + from pina import Condition, Span, PINN, Plotter + +Now, the wave problem is written in PINA code as a class, inheriting +from ``SpatialProblem`` and ``TimeDependentProblem`` since we deal with +spatial, and time dependent variables. The equations are written as +``conditions`` that should be satisfied in the corresponding domains. +``truth_solution`` is the exact solution which will be compared with the +predicted one. + +.. code:: ipython3 + + class Wave(TimeDependentProblem, SpatialProblem): + output_variables = ['u'] + spatial_domain = Span({'x': [0, 1], 'y': [0, 1]}) + temporal_domain = Span({'t': [0, 1]}) + + def wave_equation(input_, output_): + u_t = grad(output_, input_, components=['u'], d=['t']) + u_tt = grad(u_t, input_, components=['dudt'], d=['t']) + nabla_u = nabla(output_, input_, components=['u'], d=['x', 'y']) + return nabla_u - u_tt + + def nil_dirichlet(input_, output_): + value = 0.0 + return output_.extract(['u']) - value + + def initial_condition(input_, output_): + u_expected = (torch.sin(torch.pi*input_.extract(['x'])) * + torch.sin(torch.pi*input_.extract(['y']))) + return output_.extract(['u']) - u_expected + + conditions = { + 'gamma1': Condition(Span({'x': [0, 1], 'y': 1, 't': [0, 1]}), nil_dirichlet), + 'gamma2': Condition(Span({'x': [0, 1], 'y': 0, 't': [0, 1]}), nil_dirichlet), + 'gamma3': Condition(Span({'x': 1, 'y': [0, 1], 't': [0, 1]}), nil_dirichlet), + 'gamma4': Condition(Span({'x': 0, 'y': [0, 1], 't': [0, 1]}), nil_dirichlet), + 't0': Condition(Span({'x': [0, 1], 'y': [0, 1], 't': 0}), initial_condition), + 'D': Condition(Span({'x': [0, 1], 'y': [0, 1], 't': [0, 1]}), wave_equation), + } + + def wave_sol(self, pts): + return (torch.sin(torch.pi*pts.extract(['x'])) * + torch.sin(torch.pi*pts.extract(['y'])) * + torch.cos(torch.sqrt(torch.tensor(2.))*torch.pi*pts.extract(['t']))) + + truth_solution = wave_sol + + problem = Wave() + +After the problem, a **torch** model is needed to solve the PINN. With +the ``Network`` class the users can convert any **torch** model in a +**PINA** model which uses label tensors with a single line of code. We +will write a simple residual network using linear layers. Here we +implement a simple residual network composed by linear torch layers. + +This neural network takes as input the coordinates (in this case +:math:`x`, :math:`y` and :math:`t`) and provides the unkwown field of +the Wave problem. The residual of the equations are evaluated at several +sampling points (which the user can manipulate using the method +``span_pts``) and the loss minimized by the neural network is the sum of +the residuals. + +.. code:: ipython3 + + class TorchNet(torch.nn.Module): + + def __init__(self): + super().__init__() + + self.residual = torch.nn.Sequential(torch.nn.Linear(3, 24), + torch.nn.Tanh(), + torch.nn.Linear(24, 3)) + + self.mlp = torch.nn.Sequential(torch.nn.Linear(3, 64), + torch.nn.Tanh(), + torch.nn.Linear(64, 1)) + def forward(self, x): + residual_x = self.residual(x) + return self.mlp(x + residual_x) + + # model definition + model = Network(model = TorchNet(), + input_variables=problem.input_variables, + output_variables=problem.output_variables, + extra_features=None) + +In this tutorial, the neural network is trained for 2000 epochs with a +learning rate of 0.001. These parameters can be modified as desired. We +highlight that the generation of the sampling points and the train is +here encapsulated within the function ``generate_samples_and_train``, +but only for saving some lines of code in the next cells; that function +is not mandatory in the **PINA** framework. The training takes +approximately one minute. + +.. code:: ipython3 + + def generate_samples_and_train(model, problem): + # generate pinn object + pinn = PINN(problem, model, lr=0.001) + + pinn.span_pts(1000, 'random', locations=['D','t0', 'gamma1', 'gamma2', 'gamma3', 'gamma4']) + pinn.train(1500, 150) + return pinn + + + pinn = generate_samples_and_train(model, problem) + + +.. parsed-literal:: + + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00000] 4.567502e-01 2.847714e-02 1.962997e-02 9.094939e-03 1.247287e-02 3.838658e-01 3.209481e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00001] 4.184132e-01 1.914901e-02 2.436301e-02 8.384322e-03 1.077990e-02 3.530422e-01 2.694697e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00150] 1.694410e-01 9.840883e-03 1.117415e-02 1.140828e-02 1.003646e-02 1.260622e-01 9.190784e-04 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00300] 1.666860e-01 9.847926e-03 1.122043e-02 1.142906e-02 9.706282e-03 1.237589e-01 7.233715e-04 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00450] 1.564735e-01 8.579318e-03 1.203290e-02 1.264551e-02 8.249855e-03 1.136869e-01 1.279038e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00600] 1.281068e-01 5.976059e-03 1.463099e-02 1.191054e-02 7.087692e-03 8.658079e-02 1.920737e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00750] 7.482838e-02 5.880896e-03 1.912235e-02 5.754319e-03 4.252454e-03 3.697925e-02 2.839110e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 00900] 3.109156e-02 2.877797e-03 5.560369e-03 3.611543e-03 3.818088e-03 1.117986e-02 4.043903e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 01050] 1.969596e-02 2.598281e-03 3.658714e-03 3.426491e-03 3.696677e-03 4.037755e-03 2.278043e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 01200] 1.625224e-02 2.496960e-03 3.069649e-03 3.198287e-03 3.420298e-03 2.728654e-03 1.338392e-03 + sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati + [epoch 01350] 1.430180e-02 2.350929e-03 2.700139e-03 2.961276e-03 3.141905e-03 2.189825e-03 9.577314e-04 + [epoch 01500] 1.293717e-02 2.182199e-03 2.440975e-03 2.706538e-03 2.904802e-03 1.891113e-03 8.115429e-04 + + +After the training is completed one can now plot some results using the +``Plotter`` class of **PINA**. + +.. code:: ipython3 + + plotter = Plotter() + + # plotting at fixed time t = 0.6 + plotter.plot(pinn, fixed_variables={'t': 0.6}) + + + + +.. image:: tutorial_files/tutorial_12_0.png + + +We can also plot the pinn loss during the training to see the decrease. + +.. code:: ipython3 + + import matplotlib.pyplot as plt + + plt.figure(figsize=(16, 6)) + plotter.plot_loss(pinn, label='Loss') + + plt.grid() + plt.legend() + plt.show() + + + +.. image:: tutorial_files/tutorial_14_0.png + + +You can now trying improving the training by changing network, optimizer +and its parameters, changin the sampling points,or adding extra +features! diff --git a/docs/source/_rst/tutorial3/tutorial_files/tutorial_12_0.png b/docs/source/_rst/tutorial3/tutorial_files/tutorial_12_0.png new file mode 100644 index 0000000000000000000000000000000000000000..e59cd75bac8ecfe98a0a6b9274622f4efa1f380f GIT binary patch literal 20608 zcmZsCc|26_-~OQNYstQbEM>}C$P!t`zElcRmYNub!pwxMSt|Pw*+RyaJ!WJkR1#&0 zF=jA^in0u&#ZolzcmYU#cO|Jh(yB9*v4uMb>U*E)5RhpAMli*I4@aYBHR9dptf_( z2ltQH)4l&W)i3p0Z2xtZUV`TSQ6cbv=l;=*@{Q9|hd43F)N5VU2|lf}(_FBs+uBq# znVA%FEUx(RWC09c@r+O0v4f>Lv?4#aqNNi>e~ZX0nsgZdNRgx2H*hp($pk_-rpZR2 zbiEoRSK;K>@zrbkt#i*9nL4d0fjFiRd@G`*fbKjLB0M#2FV-2MqP)tuGo9c@OL47;!&^&YE9dsCSemHf5%IaFvy( z1(uxefr#t9^8!FI%yl}{m61S zZ=FZJ;~c4VaHR$>z>(G9ro0-yWi%b#gdNa*cGxL&D}=B^RWTN+h#zxifKI3lc&NIp zflgcxZm0uab7e(WV#h8?$2Tm$*P*jG(P&6%*7_7)K;U-TDklU@1+P78bRmmgwdnh3s}wAj|G!M%$TUMKBi?S#dx($ z@1#l4oQ_gEcRvjs3@xXq)Ldn(EICE}7|S#{8k)Mj@w+0>_jLIqJuVYDB9Eo#H&)9Q z6IwY=h*FkiFt>^(Q&m_=1)U17<{Xb4y~4~3w_p7^I*#0_y)+^QiG0mR{iC|~vmfM? ziSZ(Y!dH$BhR{7(Crs+yvyrkatZYbM_lr(!l8W0;NB)6u`uilnKD0n^dwjU7Rq zodm2n#J0S({X#?3FrVNrYV6NP^makD2`;Cu1P^MHRa3)N#*^D`x;GbCOzE8f=y@S} zg#XI1@Nb{eqkyMY8+hr?{F~nX)r0YFU`4wXhY@@f9;-L>+XWK3)qT2^refUfvAeiG z*TSz+b=daDlQ%F|j|OBCXAgh6?eaL+BJgm|s*~2@u`gB-)yI=o=?}>Um|cq;qr%BZ z$SiuxdT50kVyd`=HVjNqA3%j96Ae8fY3zcOkK=OGu2#0y7x)F>gJN zXnJ+X(cm0Ombn{O5dCrRs>ad)CI&)^4WW&6yrvAAs4RWOWSnY+VrEp90xeEY6`h@KI_n@N<~(Unw^dcwkn7dXG>c*~vpmat8ZbjPP#FtmK>3#PlCp~cQQjENhoP5bl&VW%T37!B#B;1w0$ za=TY@6)JL4ZlqD&mex|Lrdb;4hE^MlnJYUxR;xw_Z`@^VU}D;SaXMzZHE?7l@5hP1 zWTIrU9L(MA+Udk;p8#iuT=x1eH|AC~W;gtuUvN&iYd;DWy;&t{?p~uGe&oLEWowH! zYG4rPkKzb8dDIK5HsI~)d-?beJpp><@#;gtY#p(;?S7yjobQRQFS3DYGGI0p*!fO@62D*Pd z$yr+CE%@hz`oXISU;-w*S%v)ceA)~(-o$A|FE?I#<<`ydJnfb%B!ylM`<7 zfz0k4m_oVA4apX5m>&Bs=pOnyZ`6Zd;IkMg?$$OxuJYj272s^i`A@Qq*aXdTF$>EV z2SH_>UFNIK?5ny4Y2O@m);iC9AG+rHT1-Ow4G`G}l*sS0$NZV_!k4k!-2ryTCbT_( z1Q2&)T?Ee4_UqB1(3@B0qpi4^*R+9zP=PmdqfX=0QHQGx=O$!Zia=!?iBr9gxHK>sB*8P_97;e111!tVkdxMwBrx+WD@5J|pO9J<& zE6hJV4goUDjxSkEMKba}^lcl`a^hgsabt4hY)kuwUwx2=c$o@aasH?>kcbr25`r%0 zFmlTH!jE675fHgvWCD8_pP!YAf6C=qrN}X+1y3y&hV59493)6LtQwI~BC!038FB(x zQGDROTo4;5;^<&i`EIPjyf1gC5U%CT3fxtYFt^9YCjLt}C_WMRJXshvroX!8RKpfC~tU%>LR9q;wmEJq>LQ@q&0%TGByGnx*5^COFC0jK-a zz^6Q^<@yU|Mhv^>w$#5k zm9@+#(Pr{_;~CGJ>BedjX{s)t6pId_)a$P%xvh(+<2Zc zw+@-%ZJsu3R{Y!b>sZYczd#;bU8?cJ{iQA(%oXP?%_jT)d$Mq7K}K3ujm6{hFM5T; zwLI)<5a4fW!Lt?ZKnSf7oOw~5zwpMU%{yvuj;k+nd1f?zMcj=MH|~;iitn@7K;$=q zweM*kusOoRH~outq$k)a$jH1;%Y+@=Yzja8O>z3Din^1Cp;KEl=ZR^~$tZ!e6 zVq(lseilpJUz6to9cn%rlchJ1}t?Uu*B@RN?d%n(aCd zY+BZ=I+e2J-}2;4SjB{{sR;tRUgDul;6=xC_y_N)Ii4db(?@SYqisF&n>!BpURicTjW zfNk|tOHZ;A(h)4;Up=5c&%WrsXhvI)a1|gx9i5I%e@o>M`exG+cK`eD*c2bw&B94R z9a;cvA9tDQ5ML5l12B}=D!$#Fx^}A7ePimoIrg_n%;K;`>EG9PYa8oCXFauC7^I6I z{8!J>&|6|~_o&^+h{r0AT-FsY<0%D`$A?zdhE~2UZykqk{W0HoeEOoPDH&gKfF{}4 z5qMbod;FKZnT&#l>~#@HQ(pK}Y}nAs+XubithO)Ak2&j#khLae6!5Ncrx~dL6VnuV z_PWM+TwbauRbDSN(Sur&JM*jMjoU*X?cOXRFuRVs*2Uo~>IirVc1M|2Pu#1){!)P& z-tN8@zLir~8b!lD+M5y$YI_}3_UM1$Fy#{CmcdC3xOTV8l@tJjjc??4>{Kvgs<9>G zbLHh1W{B#X&x>13a?Ur{U0{|gSAIHU2kgmPwA#wR5@PH)z560GHWhnsaOG8_)CVp6 zVb>+mqg}2}xuL+aW{@7ad76AUV*Pj$KXqguC~a)*Be!6h%96@qq3`qB%-zQY|6ntg zz~L>2U!c+hm9E{~-gL<(+C<-{Z_e_yNdM#S{1F=GZi~$6_nF9{Yp(I--f18<;RioL zwz<%;Z78~Y!Q>of_W~(unofJkn}Bf8C^5$>r3`y2Z*}Z^n@BhF;BI3@06!g1JFASc z5!8=V?sN9cH-9A{8%%p*qmkr6cdo__4_+-BU7y#1yYD?EuIaa0=UQcphJJhHRaFd9v#Mm7d4Pq`4TzhHz+txr6*_RGa>m zv!{5}Lv7VZyp0YGIR&zsOCMEO7~-r+*FEA|V~nZVQV#E}Gn}t4C3@6nZ1uef&J9l} zKKQ(WQ|@#mEj(LN^;9WwD}jboU=J%tS5~9^cn;M4_T3fYyb4SZ z_fsq|l+LGQ{&tx_;@t+F)|-&+;}CJe&TNLbek%@!I=`_srp!vit|FvAT$qnvteDGr ze<8J!|M~~0@OaCSknC_ARu2P4ZFvr^9EH%PqwBDl5CVA+WgBT$akLiRaa@_Mv7`m* z{>IF-j!=|*a0JYaTeA5&K~pTJRH{-(r%1hpPtTJ5f-M5w7|NZNq+iONV}mnSpy*_5 zY<4*E2%KuOU5_;lqETLR(MX&6eJlRMkxW-Zk3`xx0K`Dt9lW-Nr)D#4w*N?O0rkNa zYRat5yv)HpZLTBPOS%^6BhZ4+HS{oYNQ;x!n8q7Z^w9$pLv(B=c5Cp?gm_9&h3HoVfxni->{K+`aWY={d<8lJghlW%OTJteg-1dDvSkdi=TT^61kjn%$Wv4+Gx!Tk;~41Nwi+ zv_+#A+;R#qUR~!nX}Oio6Mh4Vp9KG0SC`0}hP(FH`qmj=ax8Cex@8AFf5Kiw(?Y&) z1~1TVVfx_0GN*Wun{2~S!>@@HG6%TiOfEQuAX;Z(v>M^YTsalWryvf^oY4;TY&>_6 zCTf8YlV1{=S|(@{~R9STdQxRedl9k;HkZWTz!0PQtNc{ zt`!(+poco_Tyw$1Po1;p26ofX_0Jx=-%Hh&LeU;g64~D_@PEbB%Q`zR3!=vdyMkd& zJhSA#m1n*@uFIv|ur+cRaRlQU2e!O;OY6u9-fFcOUwrR`)*_=7KVdUlZ7Lpgc)II& zDlI8tu4*owhb$CMJ~*7px`t`IwIvQQv?i@31^Y8y7f#2n_y6;^Z&HY6#~ zH%jaDbM^ESc@H;I&7Li)jSZdKE z7~&4ub~@*gZ(u0EdeGS@S*~U-hHI?b@<{tQlEs0LWw06BJ)`u_#3-4T_h`C}4v4_* zA{()p#GO&hRi!0Y9okyU(24^LKQxGsj{Q1n=F~~J9DmQ8ZfeyxdhdE5RJz4nbMuJO zCs!RUi)>tQ;pOtG+@|LblCnMR$FB1cJV{!(0c&^aI8+XUYp;OZPCC{_|O8D}2v?sM))`)nSa43-6tF{vw0KY0Y{UAZ~Wkh6-Mmo zwvyxh*1Q{hKG+MMH_OG?2b8AY<&b_*{lTmTUteD|uN&V(hf2}6etSakhbq5M$~|qQ z7$FbYa8$Gp_>fVL9~vJI4X|C-4TU**Yn5Ad5_peykgF}XUNr@IP+jam-de?a2C;YT zthXNW;h~OE4_mm5@&)IfzH2&u=~d+mkPp_Q!rWe|Xj4Kx8Hy6~Rda1P9W9wTedtz^ z0f751RnQse6@=HX+MWOjRMKToIbK?;S{2PQ{W_tY3&u&|0HmWDyL~4kH{195-Cjq6 zNMgTq3N7xWO}-o1lSb_-9$6Q-)VDtW`%dONcXDi5YUi8iAhO#Zsp8!A4wJLO zJ$TlB$grO{1B z{xH#*Cd1&pe6>H_zqGU$V~&8A-wP$d}LxPA4Fhkps+6c&k;KR4l)(R!N}vh>z4TJm zB@4{VHTo8iU~D4eu{tWOCj|`wzaR4{A=nma{{qP8NrBax^GqrvvbkDcC4jqvWwEM8 zMs^BLv`7iyz)BVZ?WAVwkkNj142rNL%j|L|-H;k@Cq#A<|1Y-Iix^C@q)%1e`B??|gHG$GUFCE3D8@LXdEuuWI%zX`Tnb!9ZCJ#t`hdX~zW1aS0~Zhk5^5*mLk zy>%e!r3oAg=Em!l8zrbGP2EUcZ`-L9MWsw!v3ncM?24eIGh?>vX@e_9FjOva)(>V4 zA)vEW8`I;?!eOz>ONi8I`61-Rp5&>3pr-@3QY)Pqoc@0$Wqcc?GN%O=i$K z8@b}QFkmvHy7UN!D#DJ{4QRIx!8S*Kn8P+zAhAzJdAL#Ijv{+~FjOnHNJUn0$+#0A zqY(Pr>(L1mw4`##v73!@4h~g($A&KPMNqmF7Wn-dDdFrTQ8lC5vX>8h?kGb=A zU+l0jjzys?eMm?>G-fhwr~4_9Qf}A65FJs4)7?VP+`iyew{hFH{&uB2nm8Cff6SJV zPuRJLX)JufEWvfl1AL+mbLT2#4e|3aLL^8h-&aKQ>+gFbTw5H`klPy34{;huD49eu z^fa?Ayn8jQKB+Q=X8DNpQB&l#esk9*xaz7z-zwwRJ*N}_ZlzuqO!-regSjeHky%U+ z!o32&DBmewX+UTx%;lTK7EkeFvDi70T8>pLem9f`WACVuA1L}FsQE_T)w)@rCetN@ zY>+_q4)4dsa4@8GerY~izR)6>Q=Qa-b z&=SLra8BLKx9!)|ng$yV6&4mBduF8o75juSIUD^_Og~xTc`oW*?gb)9L8Zka_`zp= z-R%n2KumhOv7lz9<0oAe6u;eF2hh*3kxlx}Uo80gDz@ifu2PbVm%_o~BL;DehTJM# zuNptD7krV#H2oW_tiV7zN%c2s$D5A5MpCrX>_G#N=B|F#V1 zyVN3iiCF`I*L(Qq)ye)rwn~IrvcvMxtVR>-TPZRJ+-m#o@V>@Zf~^8m*3T{d#A%NU z-d92iFf{Ie(?xVVkkEUMZE7(`PsF60d(k_>PYA`n=rcg2oJv;gf|cS#S%!(%znc{rjWo=?W?tGb_!3yM<40}(Og|WxGdwr&(cP1wo-Z#b!w`+_e0WJdh z{he`3kxQEP@qU>W>lrnPp_axtpW}*T6iccm#|E-uaw$kiZtrc3DPpaA(XZ`DsI%gN z3+)iz)S~D9>MMtxk#Am!m%@N}zTR~kWJ|qKqfzU28(4b1luO?I+sp3c(&|hXdX`c< zfCkHJA~V88g2aWNX1$$mooTDF2a6MRE|GHQ0hSc2lR+0DN`j8s7)*nSbD06*5>F?L z-tX{s*xLm+V!PnOoCZfC`~|1Fh~rP5|6onhILl6wi+V@?93Y7Md3-+$ZCAf=*$!EK z^_s7|k+Ep~SGf)PkHw>D?|5I-hEFBOtO$P~p)86%ENH79K;_mC>1&T?#MaMwi9fQC z4m&?IY*%>kK*PF-Nf@iFn(+C+_gpUqU*YGv?W}uV_k6`yVPvrm39pl}lbsJb(t_3} z6|YBXXC-4`hxp5@!Z)NO&vCu!*BWuO)thFAN_jPFUL4Gtr2ZlWa5z5{@(=6N$qfRw)UDsXCqrOTf01JNIY@DLR-3DoovJ%~OYG$Y& z5B|S!LL>7J=1B-FHV}3iketWl`UhvO!)Sv+v=1D2Yw{jh+*)~_o@P^~R6-`KHJzPu z37h+@BjG5xaaS%wXuwew$NZ79TUqff*|J7DfG9VWF{Ldl3FX#1)^&-*gSiS`@`A+r z!#);z_0s`l0YlHe+-=@XnVbB*S(hkC(|wuAiTLcXmA1jPNqfWp$m4iinx+?nBJ;Ot z8}a!NH*V*lY;YJU+^9EW_C+CffR4$ig_;lujlucY2IejYy4#Ip3cn!SN`;x6y(!=z z-e=(b@Eq0ftp(KjBT$uchrEEVkrn6vBp%9UPq5&%uGp&=n>iY+7RQ9Xpg|yw2dtB}k^-wgbAmdmve;;a9$2 zjjWOXPVnXj#fVx8f%Mgj7AGUnK&cdl3js2~ILqJwuRkof`H?cZ+m>|kFtxLFv)wK- zaDzi2>byw*5nQkd-1GDV+kF#(N=h$34$-p$&m6ZLiT?guHBZj?+GJbs5~8O^710ml zQbig@4;p+h1Izy&&;2eXoLnFagHa175#wtcZg~-n1+04o^y@I-af<@mJ0}K%^;=6| zn?5?NG1ty{-PbOp~enHvmq2(@E}2h&&wb|P4l z0LtS%Apb6z>+Dd+5#3Fr=o3$zS@tX`UTmi9RW%tL7}Le!A@754T7aVHI0H(a%gD-5 z0}_UMj70}rSa8f?L!1az54l^9SwMN`osB|&TT^Eunsf|`o&uBU5oY5x_ zJ?#5oKm>+wKt-eyLjvy1D$!X_me)e;3`3{PVW>20H)j28x}(=mO&wbODAhS_74x?7 zs;J`p4^IF^%l@6Ttwk>Jb>8a$WXnQ;%HYwiL{82jXflDRUvIcF^lb=4R77{-rD*rf zkG9((xVqTk!5G(yCYc>mdgrdhj^f0z6F8A!;WOX6B=!mKE2j%5Z>kvI{NS)wbL|%T zv`*=!#Kj+xK6;vKCkIzP=nGhR*w%-r9m0VRM<-2*avKYV<)!3_$KT4I}?p1^fz42J$??p4Igks)5Z5q z0_j1I!echJ{$NauJK_pt5iqWcu}%5w-F8zh4#+=G3Od^iEf!BY5KRd#p*IgMf^dyo ziBprIK|0TM+bK)6YwxaJ8z<;IV_b)=Wrx>b_q@6-j$rzRG3)6O0*Vvti?ki(PSPMC zjYlG{!0?lU?{2=_5rF^}n1aDgiuEx?-qN_f>w)` z?-z6U((or+=i&;&@9DHtI(9fPGnz~9mKTnqV_ywsQ2Y&lJ;1!X7OlVQQF^VZi+wQ= z69Xklf74TbP}{Hmu*1ix-0$AFari>{m5@c>XH|K{`d_yZn3>D0=mxC#@QXcth$?c@ z3gn;aLG5g`4z!ybPBE;#Ks__8@8kf7BLDb8#GOvG+L6d>i>NIikn9{q__UvRT6T-A zr*}c#44}n&+aiq9rob5%=r0NroT6VG>{En5Aa8y1)ll1LXvgPUnapNbcBj6PC}0Gz z?Y+;|{DRao`@E^CDJT5vt#_XrIm<&IgP#~iZro2E4FFAC%fMv|N{ZoDGQM`PUwXgI zyJk_Ek+v%mbQt4z6EKIB3E7{^NG;*Dk-cfj8OP?cZ5$`_NCN)x45}}N3Js%HU43Gd zucN||sxZM0itHH_1=TJoksbM39P$|&E8xSy@@ID+B%!hzgFM>5|Jc2|DwJp0T;f|G zBH&1BDZYna8~$=0@L?%rfUN_Gf;a~b#8}^A;Cb&I`y;g8(rX=ze?c% zrc}AQY4d?Byw5eD?*t&s%5&rIXLt*NcAkQLYq)U(fwU)*aJF_PFTI;F>k*Ai*ik#$O8iLJ7$#Y_vv#W zgo9LX92W8 zci+)=sAt_sp$lWas@;^7W21QbAD5eGiZ`Xyq$@?@amqh#IC8G8hfy9Q(3g|1@F?&L zgUZwpR)JvxanMCya?xR%^5cPir}y1)5`_V?4zd4)l?BN@=7=5|vAk{&2gvTAxQs2lqF{=PX+*m?1k$JFBL)0FnJcS)T4v~Y^1 zIgdok!8o?hL{7F`jYp$XBjxzizHfsc)kJ_Ph)bM511+k3_{t`~9RG3QvY3&Rh+E$4 zys=ykP_tpz@k9Q{$cr~V#J6gqt}1d};oP_Wd7nUiUp6A=$vlm}Ap#g7K&NaBh=<6N#%HC-5mZyx94r9Bo6|S+(4P1vrpUZ`w4L?OTTuB-A5UpW&v=@5gGQ9Z|&Y(;6e=0zx-erh3uILS zrf&P0f|qWLcD0Rx|mP~oJip5di?!*X#mE{vaVZ&&__{EZbb1gP|h zA&FuTF1Q_=PX)_sF7b`~eSC%_s>mt@mMzM!`>@*d!J^ivr3s|Y!nr3qUs-t1eJEumorH`SzuaN@W8fT9!rMormy})L=$L^(K z#nCi#kbkBJRn{tSpyhc$2~ed=6D9r21AsfHszA@S3=VoMF!%Lm%PEroTcGHx?>d97 zE*P>OfNh4szQOc=eewg@T;#$M)PK7qcEoy}VRI+DKPrqTC=ZIL z1K0MAdWbDy#-tF}9)%ZJU3+Y_M90jWT*{7Hl|EpiaC;=a=fHkBS4Lq|bItKKTH{0B zfSDd!l0fg(jsMT;?Q0@$RkSAf!sE;qW67ReLoBphn;=2e(jqkg9&&xT8SZ7VaZ-lNA-x)p{HY6fm(3;8gv;X3jT5v0)7J;8^2ChV(g7WsKa8PVQZXw2%@#8q zok^7c%5sz%VPmZE>9pWcKCAAgHnD{fgCpq=Kcw5^;g zoYJ+&#(Plrsa}t6zmOH!_2Z8`E?YxoGN$2B+NejBqLLn&7JQ4dIuScD%9hoIajY?y z^^$A74B|gX?-PX-(hR$KUQJTDF4xkZ}%OALAm9O$l(?* zH35!8HayWk?+OPU&fDRtP2C)e=XmqMhZjDKJGFr4&=Q7^ZH|!yg6IDye9iN|ofy|7q$KT<`KdA6K08%m-@9CP7ddkr&Xo@{LFK9R% zI64hvcHZ5`Dbq7fZx;!aZ*6q8*Qlyo{q>IM8Uj3I#n{_ElO4uj{3@_89xP8c^A8lf z`*NH$`3_JFKS)$6P8?M;3^$94&eWa9XFGHZuTScuHWl>-x>WC+Ov+|RX8Tc}UvZtG zcpD}lxCp=FzJ33^{5;*hox?fb*TKHU!E9CxFqV3>S+4drX+Y56^I%+-L+PVd+bZ;z5Kb5aGxyu z^7hBUPlPFiS&;fBZj}Eu>Ih=pL*u`*h+z`4S zAmht0LdUE02&S39Xvi;v7biT#>L+2AV5X=+RoPw~BKD6EB!d)FNkG3q#QW z>N5;sZ2&@jA~2m!+N!^8kYs)9l@*SxR2i6|tb5E0Ws^uQm)#s&Xu(05r&fKXwL&4HO3!3CccT6&QbOLeVW5jN){bc6kQQXmUVaXhRUiN|K8|@bJV( zaN`A%mibS?TTChMn;%d5z)%eZ^h>O(rhXlMLn~4`JIMufIIQm?HAJ`dj^{rZ( zk-P8s#c6*Oi73&%S$qi7Ox&Fp%{qciNHcbf1N-X}4}Lcgu*SXLqyLoiMxg zDC$69dn~}U+}TnrBpbc3Us(CKeGpJbRmD%W9#&~UZUw5mrvpQ+6%{}tld#9#iKe8y z?VG0h6eq13B)7)-8rct%S#CK~Wp>_SfhE}gh3#6q;tGtfk%S=Iq~`*XUd@QZm(07= zmtKTde1GHe@V~j0585B|*P%K)I{Qr>(XR$6`=|^}XSKjPpqMlHd!mY{!SFv3PdOd? zBGq|ujiB1Jtaos(KI*VlXN+r(NBcPv9XBurpE8h2ABkKglqd<_6^@ zFeKk@zmcXr==34&d8soF`N-S%jaW6v9+z$1E$4Fmtr(rE0W3xDb_VtrJby1A>jZc% zVVlRL(mA_*aa{fS6YqY0`r&lIIL77u zJB~*%nlkD9%F5wGj%8gaC6()JuGdbD@DHjB5Yd4P>{wqWb%VuP|m zil77c@gj=AH=OT9PnD8utmv+)O#A5?mulwZ_*Zuu0$pkZ2!Ld=qkre$z_HdHynCyTk!Edre8PMd7OMUchWeOvdXsZQbZgdVxy4r|JN4w#Mxx6nKnV1W1 zcK^{z0j`@5KzS)v)eo(GY#ZWR*IlH!M#^mg^OwBJa1i8c;zhH=d}KnZp5dVHl;7X) zxr`a`(rRkk4GH2UZHy;(@y$iq_Cu#%NK|#1+^igb6%bk+Lr0~^6`vnFqQ zalUeIhuOGgbd_t z!w@y7dIcU?oAxR1syxkCRet28RbUUKva3Qe0imE}6afY7R;K;6sHwPeTg7Qx*4kpx zT}KehFXeqW!dz1|oj1iKQ0ATxUeQ|1b#Z=(5j1WTZGapkoDfcU`Zw!2mCr(ad9c5{T4~{EAEBKat_)Hi_D9^8EmDu! zv0kF3yH>Uqbqa@{db1oEjnX^%T2!5%2YFI%*a&3y<;7IiRPvPNvNlftL_))x?%Kuq zq}P!L-Z0=w9~N-22K8a6pXGcnn>NSw$%40!34O?oZg#-E_&W00=88Ncv-7&yWpd~4 ztr<2DccyxO#DK{}GO0JY*l{NGa)?>i@4xJ3vA!9Us_-wbgv}E@@bR}FZDiOBhkWl{ z$8||$n+2N>J1%XroahWOhPXxEvVEt-nMQpXEBJs}^D0|O)*m^5kn^P+mWIX=s7WTv z%0{fXvx(d|8Be1_8f36XX>I#wk)S}r-Fu_&iyGFD7K&`!e`K7u-n_pJeV!w%B4}0EVXHnt(zYV2; z3RhPUO&@I3XkG&k}ahBFrhqpS!r4jS0VljH+H)R9$JHhhA(gP zw_7ws4fsZv^BDj&`^oDy2qTZ8U%K;v^(7(21CBH-e5vD<-8M#FyS3Gc0q2CHqS?@U zFNv|;ggf`J6+QZ`_7>X;Z<*vn@X>uf+TvN%z%ubF%0)iK%P z%Icq}$@xe2<3wC{FCkto&Ci8f@fV{(BH^wKirfdsmsz25^#RAvORRpmJO#0sK~+A3 z5MKV5=VssTOW9_JZ61Z}xmA2YhvcYs?ftkXjtCSD-~1-su9v^bt^?dD^pMUXe@k-i zc4W$#=XnJ{39LgkKJ>EOV%)N32*LaFd3jr=8BNUS2vG;E@tkN1gZ=+3@Ip3*mr13M z3zX?!nKzLW%-1#OyRL=L^CMI0*PfoU2FH>K*+3Y3FvX^^+MEWdB|92k>KVV$z*uf> zE)Vi^ldoJQf=_B|cC4K$veoSCl9KgfFaf9o|3W}hc1?_L5%WkF&!y8xrHlRwZ(Q*X z-3E>M^gK_#%5Zh{c;!%=?5Z?+qe0p{df7amaU#iN<@s_vQYO?dv%l)X$Cjot)l7{m z8lUmOQ(Scj)l&;Xf%U^v&YQ}9XI!75Yll4oj`l`QAxp4w=0~T&51ncDD`wK{Mq)(E zW~iBMN#_;U6`PHCbGF1szMTV#^MW`|jvVtO8U~#=FS_g2&=FnnuDam`I)HjG3-sto zIm_u^j_<5E6T>L=Sp2eTL+a$HO{J^@_V3<-khHwrb7w7$cvK*MJpF$|f5#BIGl_I&h@pd07TirDt za*UJ~GoZshx}wV}zEUY;#Z`2a3KjWE_xk+Pel+}A)bGU;Z?iu{o&<`7pQMGTQ;$B% ziuLRX)Jpu+|5ZrT-;jcb(M1Y0P+8w9`n_qYfGgYJ&YT?Ow6I8N_sjJ^ds3nP=fYB7 zI@#WF+K}~EF-Q29n`_cO#)wJ>NxoOJgEL+v3V=V1Y(;6|J|5ldb2P&Sann9Z#N&@! zJV_q3LByPSr=N9OHNp6XPvwMx^cyS5%MQbR(|KU4A?z_TJv@?d zx0Qjshu+c>vR$Xt*v-PqyH19RkO}L~r?ZO=Ii7fAORMese9E<`>FI*_+o`!e-R7^A zUCdB=YVjc*qrt&Oy;j<@Sw{7VUv4mz+vl<1boY7LqXLA&|nzj8AGTZ zyjCA<8hq9BVQk;K6S?8P3oGuORr?<$ggkRM0&O};cgkDmh0NxMJ0P-WWoDVo(;Fk- z=3N-c52gxq0g&a&*UH_L>THy?eYTx)(L`bCr8*U0;!>tr5tvKLP-iG$czCqDEp&09 z;n{C*Y3Y)T$#;xn`5QT7GxKjtG?y^TYj(<=sNtVYIqQ<>y_o!Tb0dT3JUfx7Xt%=) zPbr!;CVB>_xO%S13B(lvx2zs-EdG%dZkP^S@69sQ(Md-wLk+&}Zu~HPo{k?Kbsx`DiMGeuXRqt2P;-Z0!u z8J7hw5t+M}krDY=f$>?uZJ`PIr^l%_uRDZtPvGc870h3@d)mlHC21-S=qev-4c2Y_ zXd5}4TQuX+{RwDa?)oq%qq1bJ0)-N{SEY^&BQbcJZB2-{UWnoLtM`>}>}vA_YXBI3 z2onfeW#7=_w}O|F?HU+U-U$d^CCMCHxn74_4$8YY!;o_wJ4>y>I{&7UoyH4@=tn8r z2O+UcrfH4w=AEWbK-FQ5nGwT|CfU?d1Q4oM-CgUd$-ULij5mbK_EdZu&$e zPL+w(B6n_~^VcKNhqX>c>}}6=?e{I3X);*NCPKwm&xH8bQG6Yx`R20@N30qKufJ;@#KGCf{-^;&H@UpYn3tOjZ? zk87Qu3oMA7eo%#D+F^FnvG<08^09o=EO@$L>kGi%p|n(oF}Y;S{8z)aigDR!1otJq z7v4oA-Oef)GMsQCxE9?=IN#Uc$&Bg7j#YVeJCLLU^Sf)S9Tl9fF^@p;O`^3_0!FC! zvX!7-{?rQzIQ_0VRV}v~2$&X{LWp>et>=q>m>_d9M$8sP-I~(!eb0HZ^LVHMHpNN?ksV0Nml6t}9-WZA$Oa9SI&~9c7WbhoE0&vNu(iK`=X{kH(Z%e& zaUP9K4;Q7{g{36#PbDLWdXFdSO*Hl-GQ4g`VZ=P(ZFuCyOxW98Cr9^QxO>Sc2kofr+;@gEW+j|qA)Lr%afx}iOG?QOVu(?YUo z5MCsd$oS$dYuQa!)StlH_U{7Qh%GYnSHGf%40Q)q?&l9T=_e(6H|J%>!d*_59J{56hJ-~ElM0{nRKGgHu z*pSwf+go(~o_iBRi-H9Myy|hyLkamd`33+i#jeuKJ4O6cD&D*_VG7Yqh)jC8QM$2}gI=41l z@$NWz(0Dd%DnuYa;p3N7#`ez&p!997ei4NB!lm2fYg`RFqpY_ZyOniOpvh+vd4Sje zF|JDK^)(5uG{U9B_FZ-X>aL!PgO{|9NQ>L9C`uoL{(SU>bTX z5Kwd7bC5DG9O5+Hw@3sa#O|gyi?d;F*;DQ0L`jQ?B&!|qqvF_>mnQ`vFXZXjhofqR zdq>yL3xqu7g_slK-U3B3LX0rJh}pXD?a$3pM%RUttFi*RbPo1rX9C6dda(P`Zk->y zE(jb{H|A(6N?>$)q#(6GLwXFc1r~b0Fe0a5<_J&DJ}zId5Dukh(tHo zCy+D)V{N@;El>4!fp2k3wc;ztZ%=*J?3_|H?)nss>HK2}kgN)gLQ~RuIG|GB1Z7_z z0FQu9*m|(zP`&C^&Rzj&6Jw$<#Pjm&gqn_c%e!VdPqveN?m+QZL6ianaUStpiTXgp zW%CGJ2$_}K=L*~k6Vf5uo?L%uTdQ3Mc!N{GNESv?2V;WmL+>sPpnlC{# zXJMJp6TaMSOkSqp$V zNNR}t)?)v7Vb+4Vz5~0>UaWoo#J-IHj zk|FXZcBZHZe-sdAJg!q}ea07{J$hG#zna1UE8VH6k3xCmV9RoOsijU9T(O3gHL@Uf z$&|E7T&h8|^{9S)@)Q4CLEOumM4z5R7mV_n(YtZqxmB;@@`s%ip{*E`JfV>d)fiCD zb=e?(hKk8FxAU8IBYVMez9hXdxH<>f#2NZaomPVC_~zEYqjGX%ZAx0B4IkQ;vqRwH zB;5oW@j%d4IvBV3Tf_+1A>$O6nG!fDzngkrnROaYv2EdkIxvvXro;-G#do7K``B+v zs59pXWA@0Hh0l|-jIqgi4c68P+pYjV8z<~ok67}RyGfB^__;01?!3X4wG|1q!e!Rq z&b56Wnj4kKN^1$2Q%vfTj6le8Aa!OO(d*+Akq8yQ9IXF;> zoN!?W8?c%5<`92B*eO4n)gD*&ey)B7sy3#g`Ii{sK6_RycQ{aymsBe07#*+Nj%zIS zLFZ&vSLc(d>gybJhSka1Otf&&ZAGn-FSGeFr{nDFwNx4V223T;9nAML57Z>zhyD*O z^4jN_6D_)GCFVTe59)N^7Zq4_J3V00--WKr*eAjfNpn3z+L44yHq(1~v}u6!nhN4x zhAW|QZUup!zm^?POVL2tZh1dy6bidG`5(4peV>`?OGfBYxGNx4rZaRe1VtQ7o7dL; zlW9Q3WaL8Kc|2pdMi%y*gOfpxZ9(ppH6>d%BoWLKO`-5WCYm~vL2y@MU5Jc(fa8Tz zsg)G#z9WsZjn_!Q!L_8N-EC$-DsPi#ZC}rmIbK&+`N}|<6{0nw!XoW2*tr8z28tUe z;M*iv)u?w*-eU7b6Z}-u(@3VJJQ)s5AehCkzJVl^N-cqfw1{ z!-MZ;bhV;`Z*U%|E%W7csR@ygbBq8K#fJ33ty0Tqsk(iQw*d(^aK|ai#xfcw>!ROb zrU#DO{qpfvvU>#+u zvr?2uAEwf}tnWF0aZ=(LEiUU2)aRa@Q)B$JBiy&;@#O4+{2E5(32jGMo8kF(m~!)R zlSKzOyQw$VtLlvNWFn^Mpp z24+2{bl{uhbgL<8$50Pi2aSA+Gp~v054W4s{lSXrSo$yAIGkK+SxKp(Pn~rc2yIb* z7gXL7H^!uRAfFkP%*C?Z?kuXn=du^@@XUqbr+LbpU%qg;P~dwBTTY_PAB|V$VBzl{ zts2oTJ^Io1x2;={hgTyhraR zEu-3Uc|f=FL>0s}E$lwW_&Gvwzl3U0L@Kci7eCY|sd87%F4LQ9i4k75Z1L^1$2lOj`R zYj^J$>(AL$BdwKCcuw9h+ppD^4@R znoq|&m9(CUv|td+ML*%&JC0MzVy%ydnB@ea`+B4UurThn*v3n@aaO8p!mo7R!0s$w zDdP`#X0Lzqt3S4p=d}o9HawoferTGP!~iQaFZOCe;XvEry+07+92xL`Q@sZ6wI&6Q zTcM$b#bcKA=6KFM=liaMZ3S?8^YqV($l}2Ov9<#M%m!#Q5Fh5smSP&RN@6W$X+5~ z%@a@J5j1@QzNqg&-wwrA7g~NG`@2Ga;8UG}wU{YGdNY1R=do%-b*x@R#IM9C4O_M# z=+Xs^T;b1j2R&jW^p|tPcTng+AYJDNJNX_z2&7kKW4-t%`wp0Zz$^aIy8jIVaZ%W4 i_+-Pu_)mBqV@*upZ6Wn;whNtq1KiBYwA$G7>c0Wp0_sTs literal 0 HcmV?d00001 diff --git a/docs/source/_rst/tutorial3/tutorial_files/tutorial_14_0.png b/docs/source/_rst/tutorial3/tutorial_files/tutorial_14_0.png new file mode 100644 index 0000000000000000000000000000000000000000..f9b6870e71d10e2404931db58f220d6ae615f591 GIT binary patch literal 17833 zcma*P2|U#8_cuPJa+7qo$Pl_KvW13($&!>cVQe8w2_aid);fB|>4sh~w!eFoi zI@(uIFc=#F2HWNE*FNygH)|d<@JGq#>J1-b4@V#0TlXDc`nP;M?|Jy#bH4r01Bd(G z&K~Y^(ify<&;E1Q$H&tfDI??d-xo-G+;@@@n}Od3i*R^qTX@4@`;S1cU9glFa2V|B zF`X;wCVnsGM*V|sqU*L6S65#Qm0e4y;Bc!p7AucCHoW`P#oZC1Z(fOip-BF$38*u# zGgq&3r|wN;EWh;QIiyec^%jHw1SqBJ0!Qhy4rj#FUu9F=ZLRO zkJb{3Rg#NcJD(prsgkj9yNJ9S5CqG-shBGby@7?k{R?aw{C$1da4+;Gj^~jC^f$BU zZyMvq5$AaJjQ(eUN`)yJDt@UcyX=T=!WtAN&-! z{Jqj~3q{-0`=?#+d<}{^`+^h`=55B`l%&sI9hKDpJ^%kgVcfESDfAP3?%U#HRWh(YU&IKkZp>qe(~t_AfdQezq7M5MHRz( zexkX@7FtLwOI!FQi^VeHMTxY+v%$Of^T&>C^C9ztroFb{^Y(_~vyTagsZ;hH#c642 zxzgKA)2wT2YeOae*yq4!QcPa8a3$>zd;jORa~*B61>0QLWvoqSP+LVaiu{+RXnlGz zl=;yb5gmo`pfa4w%8kxc1q1V&UQdN-i_Isfzt@)JLpFn?ZoIlQIy&mFqCnh)7P*xC zs}CbX`|f7FcO#|wdeI=Sbc*+<=$fe%S7zM;#=z-?*E7|vwQh!!;QZ`hyIU(Naa-zpsY64B zI2?DHR}&vT$CXrPW}v+IMA&(()<0k&aGZE?{*8u^LG`LXy)d7CB;(f;_xd$Qw2#l7 z=;;WsSnwi~!J62AXrsSI*O~*+9ol;;v1{71W@~eO`NI)|nfmqle)9&FX|YY6rNU>= z$I`VwF01@#WJ4|w7%%r|9CTip9a2(b7o-SOHHkUsS#PuXsEwk4nVVa2deBcDoj#e6 zH9d$t-CWf0Ylp}#B;zb*;`q3slT$v}S9$P8fS~VZFBUUjt@vu>G3U+{d4Cl(u%HsR zIuD?zai)h#&Ckorn_V^t`gVl1)*Vv9AT`)vFqjmC{j*?Je}+WzH9BS?Favq+$$o6iZ@7#f81QAHEL|AZ|{z zc*&~EJA6m$-QVP=}6t;^=a&2_Rq;e2Argn%`hlss1Rq#6LFt`$S^cbLUK2vGOb-X)SFeP@}BsqKA zPwc&F*m+2B$LzQtx!_)}7aZQPpv=DrO8hegc_j)G?dP8B4)`5U?z?f~w-P1!M1Xx$yZWhekKC3)c{euo%CL+JB$mX($kBxD$8 zWM)dXN~HLq9TJ{{Megr=w)sCPdL6y(^0s#+l(^$=47tEUN*zYqq`my7VCrAY&Jv1X ziMD;F|G}6mJCBNZ40K>+=Z0_Kh8C|+>CjOEi+BhrheilwSeatL8}Qpp%Qk%;9vZTs zqqy(d1*YuT)*Kio0r(d3imK}z`m5)gfn0%}#@o8$JO22niN!ujOym2N^U#upNB>Xi z=A5DLGa>3;?g-+ex?yKb*dMCJ!1w$&mJ|P6_6qv=Nh|(8KDN0V#7cObgjPA%22`Wm z2_b^t;mcul`GdF^{Mo#1u|WcUeEQjrJIXp_vh+ofYhy=(K#3pXcJ2|{z6bS(fscS} z-^!!#f2aYP;D4@3O#P4ULq8%=5j*D$v^?ATj@o?%t2gfWcroQiEk3u^{HiGXAHV+R{%@WZ`J?tg&Hu;d{s%WSb^oY4 z#AnYx+upmqDJ=d!C35SJ!-4quM}`m^-Xdsv?LT%*4G!2Ys}wKg0_`v}6f6h_CYrh9 z2Drg}Z0hpy5^jRf&*xJ92HKG1gUbEGX$`|pjf!mDKVk-M;+7pkWM#+PA+K?=Ch71W z)qe)qYPO?B+mbNdF_$0Ufh}D|e-_*e7G(P)73decKz``C+q!|(5j8~;ppS-%fH||Y zf`kfg=p$SoBKp@9g2@YgbSYAC|LU#8i(JPF|6Dzx?^*?Y4Ue5!Yh(!28+3 zn@mOuinlLOIz8LsfN8Gue&{~#_UA}KQ_@Hjq6UkQnMeQLZQns`NL3awl$%ZG1MvEI%|qUwEc~e6eP^Yipxsv%a1@-~g#a<1IdV&v=GT zwcU?h_E<^UF&&S1pQ0SeC`2>AJx;Faw=*f4yhFcOpoK_8ks!Xx(Q(=TyIWFX7Kt~V`i4^ zS!(k_kh16Zb}AVq9D@>0H}Ux*n-1P~2t{&NOC2^4 zPx<6s<*s3JaKWZ-x->lQTkF8CjBbH1Zy1U{6*YyZEfer6+-Z{QB=U@RQ1=|$>T=)< zW*F79kL|~btGej-S~X!Ec$fVtG3kikmjgQ zaPJ)QO4xc$EMf|&L)AfPXgp3hYk}7a|Hu<{MVL>@w*ImpW3-PFQ(nE%s@i#C!DX7* z?L18kbYIUOGK+o|#2Tq^->}okN11W@=d?$!(3oco9z>?y*0msRJ{dS(ITW4#W`MVn z-g73 zEv}Q^T4pnIG@3XgW0&9Qok03eWIXxxMSR9BSE7ZhM}p7N-}Ao294is$ukW6&K!Nmq zhpNEE5-s7YKEEIBwLD*3atHmJo~%y~ZzwoER;n%9Suj4mQYM!e_@p!IY{*XCu^ z59D2;9AbU@=skHX)1FsfczL)>sCaC+(6lpccwms2FkfI9_>w*o&9FNi{n*A_{v_$0 zI?2VgHs$aSi7%u57F=BP%iQ%H4}|JUb|>TRtn@J?p4WB8G^5M!v>raNx+Usbd{~@# z<5c|ouQMw?53?6ln)g-Lvo)XcnBfwo9M>z5NvXn$P>$;th^7p5ZK_bS!`Xbp*ncVG z-l==6WFx<{2HTi_HOJ%{Oh{n(&1vQs2Ub9Qh~HSxrwr5!xk=prnQ(j;&B&}a?VD7m zrph)qVej9#CLaHq)oOfPo+MuMvHze@0-KGwXmUy6fxZmwh5gfIwDYJ5|@O8!LMITOH&_Fk4pB6Q7nfj%hc0T($`eX*-u?^)( z5P>6rhugm_@oVPltMtGeL{kfc-``q`XY-q?Ij|2~<-+>AD3#tP0WUewE6tRl6Xw*+ z_7ZL7z*@3s-m+BP5MzB}Ns-_Dt1%|m7Hrami{0mnV!`QmyY!1V2nd0m51xkV<1;GJ zkP+Cd?V&h{bjmheGe>6HoTL#RnaE$hGhL8yu0ZuyAGh0G8fBGad2X^Kkpag{xi#Mn z%C&`@#@u`)?f61Uek(eQO*QZ0VQRMP1KLB%2zio{aBtQ9wqj=XaiUKb;)qYjAxaAM z=wMv6VJp(pKDilfGKVB^`*{ZhDJS>Mk0n_K-gcP_Clo6X2Is+tY?pR{5sDy%UYrKL zo?H5r^4vU%vP+(_Po7+Vxu33)+|CL1(K7McO%8HSj@s%J_>!xRlP*LBVS zi#QR_d@gL>H#az(7-8njC$W56w=j+(6dk#sOwESc`-C2P#)O^Q ze_UI+FwOdq8{^-*TILNE=qu%QWwsZQr-%E-Ncn?giNoV^M~5ch`X&3FT`<~OIMRSo zlqs@0OSg<)LTmD!4cLHKNmdCf@cS;D+CQ|1kIc?L;$m{eiVU;9SGHHq1vMM4c*hol$*SL&Sg!sC*tIy%lr4%mgEhZ-*v%! z4ECrTbTRieAyuxp7+Yo7crYjox4Pj^!}YM) zWpIwbH6>3U%oSjb)Ycf-ZY}XkrQ77|@m_0VxC$QPre|~~n`X9!GCzNQMi1pDu&tl{ z5r56X&3&Ptn5@)A0y*EX;Ap@s#d>IMY?((@kb2Zi@sRQDch?2gsx&8v5z4L5q3HlZ z+c1<=-hf~@vLT`4_{?5|(aFc}_)jH8dL%A}xvRWSL^6sC~ zdByfNG~ewlM@prbvwjQ;tD79|CP{CMWj4+$0D|PTV_kuo}R#ZQqmVVDW zXELYXkTg3CW|t$;p6PJM7hRhL;A6uI?IxLPrn8a(=rDQhR}BFn%jVAy<b0h8 zjRd^s*09)T>N?H9|B`3jJTsXd zDN^pEk(y9D{4&gTtPs4!sR*QVKw4Kuf5Fn!`I&z`Ukx~bc=gNQAQQnIeP2gkcq}vD z^F}x0da8Wzjr4d?=lUE6(=@RD`twlk9&{F@xb-42xRd6C|t(PYy<0aR*VPB zGjDk-8DqQ^G_}niFYA34^pG0`i`92FuK>Hvexy284kF{Y{UGk#_wBX|{MlxGb+*-L z%#-WM79lei**-PL9FbV=Epwusb`mT!ukd9=cY9DyWQ^x&QD;#t)AYF|T+jE&yqPEB z@?6D0PGShK)N759u#*KPp9Z9LptEtNdoR|p2`+fC)Di+w@@0~Woq)vC<2zaXggMsS zPe(zDVWH$P{odUDX9^Rt-U300fPCtGdRD{7kdnko-@*%%lpOwJk9Ki>^uhwtNXCwX zj3wUAUX0AWevkCgmoi&4hf(IVe6^i>cF9fUtKiVdZU=lihXs5(_x=ri0h}?wq9fE2 z9-b;>2Ks#eCJs5Ty4Q7jzN$4>ebuN%fwrIWhN4UZwo>zG__f%i?fbc2pm@*Wm6Chz zMhci1ixbHgtASz>3cA_;nM*oNSpFC}eRs&Io1!qDHpWeVR%6bHOBL$X53$L16Wqv^ z22rA&Np^Nh*uyy6S<{fp2TYXlwaR@N&nQkB9dF zkT>*Zyn`WmaLlQ=PFjFFu{!QCe+iDBu%fD8OcbMg!&;)cCuUwhja{ZjXJnnIU${bf zR?C`Ik)aoP$n=sD-OOHaL(c3LAP9-#>d@5V@yF5gR&4YO5ArN4!6uu%u+m8*1!_T) zCL6OSjyd0qw^bOcvvk^`9^P0PlG-NW{Uwk(%5kk1Xnxa~n~iYM6;$Vea{$X(n&Du* zSF6*hYzWD6d8CS=uM&OceTAsjk}1jWA#K_NaviSOY|D;ql6Cq_kc7s6Z)1Ch%8vUi zXc5vpLDNy+nT4qpRa2Ge{G)P7wr|F#94$5fy!y?qiE^UkOILT>y!K#bdAiNG8T}`Cn(iaOnr6h+v?Xnj>~`MciKU2sg7{9Z&^FHhuJZUzJvwgaaMHv{6Ry1sbK zcz8(A=xogGgQj*|i-9+-z1z;KEEbkAPk5-Ax%xg}jH(4$Qa{&U1EIJW+8;df{q`dZ z?rMavSLEFSx6N|+Sz~^5_Ah_P@9I0?kycFU6bxvV3h6ghvmrEZ3tzM9H5W)pBe96e z2M-=3RL5!Dzqa^nllz{&h87RJ2DtGOHT%5u-^?La$A}{o$d=S$lE+wT26W)*B1Db1 zy+P&|a1dC!E|d-0@N$RgS+||!@urM-x?%1MW{1naeHny-oz?;>#JHYOD5{38s9&D< z*Uc+Z4W8bQoJftPRfsASl&RU!4}-LpGwb;2Di`C6`gaK zlp9B$Q@@|zrM=`YQ0MgZjtY~M&k&_bn;5%k+g5Gyr2Qsq`Am!wnLjH^P>FxRV9yl{ z3ircl0$OcOgij0R_b5wNNapVZ?mBPv*moyARw|;3>r%{YjY1`MVK3n9LHv3EIi0kx z?hV z7DtBz$%Q?Il%>^UW`UI;WL@z}hkTNhKGc(-QlBmPvi=52P83v9xb*c~@*xh2Y(2;e zW^G=zJeO-$sZwPG+`vTjY?uLKRWgNgLNB5E8Zb?b+cO^NpOwuv=4K;_#R_@J6O!aW z%bV>s9^GRVS-Ix;^WFjgPZSVpMj%~R_N2)NK6-$@%3PnWpm#YbFk%}5)op86MXb+< zyezqK-lAEbvmI_szYheo0uHfQ20ZuvY=Qe-{e#V^h{wVakFf5A-NMd^{}zh;`7XDFLOFBvz+Kw7tZuZm%vk!5>1bNWh7sy96o9-^=J8t_CK!7#9sykhO_kEAD+w^=-E#@>Zf<9@v{cB<>AW zNRdMnfuPa6J7`J(vt+$=lv&MsS0_WuA*sGJ7$@rzQ!naSp904Jc&1>gwbdTyqUUR} z93NbDz1p?JiW<8tL>1P*6OC?PuLXqnf7*!5zH~Fz)=<#4%88r) z%^CgLj5U_c5Y`dBE>f@?9^5Sdba^??zk1QFdFCMC0WLR{NgF2(n$L$^T|j2(ET}j~ zILc$k4xK182{^&vNqykjiE{_c>}6T243F`|kfA{8Lm@duwdC?rip%+|Hqys_1_%yv zr3GQ3Mj#y!746PynAj(kP;RGNe6x9mM4BDwSma6h7|{UJEkJ&u$md_{r-&L*e(Kkt zoBQ-x4NeyrP9Y)Ik@-{_27JHQoB#U~c6R=@aQXmMgAtuIlwRJ%IPP3v ziS!)4WeAA_!RcTnlm0c)m!5oJj`guWRahrJI66p3%R}}43%VpGV7%2>zVM8;-(J{1 zTEGP>;o{o7DC46TFa6`?1};6*^o*_-g(uEhkU~^e#g3tN8y>?VKp=QMGLJ%A@`$GU zN~SCmoi5Pq8*B=_;#DF4WKH(dEo>KfU7~$c(q_EafmF#Dikh>ed@j7I?|x{wjvS1fK0UaA9ko$r!C|2QHZHhg+IRJxC z&?PsjH?hgp^W}U2-^j;COMS|MnZE)yek?AhSS=8W2^P~rBd@xG#x?Vp;MCx{dy%Re>l8>5pLHmZIb`z_&bwaDJwtu^ngL2tp(<^%2P|9{VtvNn zJ>9QnLbv1=p=LD{#+jPm8PAjmf6kdX3t+`$ErTBDHHZryt-4n>|7$$}^V3(RWg3Ze z8&phF%QwpJD56eQ8}>57MU4RfmA@>dVp&PlVps)~I|f#Q5~bOxCAsP|lN7JmcdSHj zN9VA-IPPN>poRHR!GOv>rffZE2?ejBwyu<>K^q9qVmDi88KX{no%asP1=-xd zBXKg)%ca;+!L0Pg5tv=a=|kUJ(b#+3k>Y5rBd^UI`q3aP)F*Jzk4y^|uZ_78l9g;% zNy8;-q_GmM!CAa8KseBS%m6{}EAe0v*QQ^bz&EvDd~z<~ngxXOj2$=Dvb_cpG`9up zT`ec!z^dYHtYPSqgbn~mI$a6|&|2L+U>lG$@`k^>GAKOPGz0V<>8y!axnV9`g02ZD zcXx`m8KD$YLDED!;Kq8tG$O2{P9jc-hM~U=_;^|?7UNXhcERic5O07(=uZ7@Qp4)s zk2V1{YwTeh0)*fQ)`oZI&^rHViuo^7L2a4Kz_4XC zBbtmOM2sWO;ORKwfuR%B&EOBLn_5v{L;2tDl6qlD43}3IRLVjDvvdL5HI4UrMdbD? z8_(Z0Xz9kl-u;+28sOgVUi_)tN-ZepFj%-n{=GJn+u~KSMo17(*};{i=}9br)alY z5?P#Ds?n!27zGkvOx>OufMXT_LLg5zvHjYA=xZqWknQqT$sT)#vIqErg2?C$F&&T& zLm`7oeZ^HsejyF6@v7@#yvqsp{62H3dCtA%!$!+0AfGq~#k8U-<4jL9D4t*u6A~Nz zzB#t}30p{ihm(^_Gx~DPH0bqgzoYB-1n$xG-xo&9?i1&2>#V)cm@(M$uTw>qL}l}N zVCU#RQl$(^Vn=#xK|vxkOE&oy#kSTDmoVmstA9ZAiV#7uCKqgaM%;r@(ImI+DB>un zc0v;T{yGUId{T$N*)y;Y#$^p07a=xpM&^W?wOY_P2LQUNS5_RL7eFU#2EOlt8U78q zyQazW%u3M|Zo<8*Y_OK62uR9_-~?%8@Jr~_@>yKX=wl_MqZ0SQE_Fe|^EElR z<)u>kBuL0fd;PE7Fi$ZwkLt$i_a*d*{C+a=B0;ZyO93DSI8@?Z*rn^*I{x?5n9rDD zM1>T)PcpE~FI_4X>Z=<|sq*$gPgQ@n^RNarT<^0~z{Q|OyA;Q`t*a9AMB=ymA&K2E zEIQWZ5o<2Yp!fTJQAR@ZAR~Kg^~I^hJB97LVRkK0T>nDgvCH{F&J#eFE_v;#Eva9e z_ut-gH^4N=l`_WW!Z!5R;A}p;owLAM`(a{MTYke_nAM+amZi>>W4bt5|#NefSzt2omru zDDs5TC6~D_A#3k-*{g@ozS;v5vxOw1skD&Yc3T$~`qgJ}mV@=yRCF6B+-xhj9WJn9 z&IYrCgQ8g+qi1~c2I}N7GAIS3N~aNdNAN7>g1Pv|L3voHtVsgEPQ7CgQ-fVy zU8uzS+Epd2 ziKBpE1Hxbf90JV?^>)MliOh);ZO)N7=!7!}IV-N?14^bM-_e4YEy7-qhBNC&(v9?d z%Lhbu!Tt%)Oi`fb#L~YO{RVVDmVgR*xfimKnYO_{gwf3A<(RF&>zL9j(7K_SDKet{ zD3urfgpK4(ekN4V%=PuLk5~2mihuyRO67`@f$soDFzieOjS699KV_f_`}e60kCeQ% zSXZBRd&V+2#PqDWVRl$h)O5C+S@?RkAZ~8E9l_xpJli(Nv()9X4Q7?(O5q7ltQLmA zRKrf=Gx<~JnF@`dCFl}gA-CUhly7JMaN20~&%s;}h{bPT2UNiL!OX&~J(Rzo$2!mf zY*+oAaD>@cA)U6K)=VZ!qUd-S%>HknR6^4~%r~Isj?5Qn3i(D*E)U89b*mLlU7G=Z{BdiATz zdnIXvSU0z%`k#ftzPG`%4QjuO=@SpqG#YM{@T3mw0Y$gt0zRS@-8cU3JxGL6J&ZOu zsK00aUQ;OvUV9sRlsjO#O071>edWep*ik;PPj_yOko5_jDNx$hjB7=o9A$$1`hipB z*xKU;r~Fb7&V!~3gY=IFg)Ayt=4uMi-PXo~lt z((iK+HcEx}gBrSKjyKM^Cs8Ij8^RVxLj1|lmNmtIp|Cwj|-?z9H4*FI?jlKc( z37CGQQ~kYk8K93lveE?)5!87BK;mVsRMLRF0jBR>8z>(sxeIn5 zSWB#vnzzZpl}M!&ZJF_fe2B46yHfiLqB8>H<}IpB4;GCX(eDHB=cu!|0hu!iMYPOh@UeF9#bXg%6HyEvFW z6TEODN)9vuoMQoyk+~tsd$)Z)*cvk;rz;GQ5Wu2G&l7$E|6ob9wvL%yoWgC$baq3g zbi)hanV9L<6HpYrxlv+zwqy?sE(@6cSv^vGPGG=jv~R3(eEr~A>dy6^i;G8=W6B%> z%wo0N5KWSSRO`^A@;FGfuH^yDiSS%zRNM`e%vLuc_S%OdG&B6#djr}_53r^lljtOw zg|lcKa1(HZ`S%RTH2*?S6Vb+3vsN2CuRK7AK>bHee^C+DYY}@>5G9wXA$ zlGU9Ysh&X`w+@LQ>=uHsD^O#yf;^Q2wmeSM z;6i~esDOg(l)z?#7J8N42JQcoi7kyvB3-l)m$Y8V&L0+%MuB*gqo5fG2P`7s77ZAJqqLCg0@il9&KfsCs87@YRLs?m@MB|7u|x z1i___tTeDI!R%&FO7Vm;;ApKrZ80ql<7)*6za0o)-b?^|h z$ovvH;!rH*W&{hxp9VKJBZjzLTBeP1$4x!XmX6z&_g+B>HMeMhu<8=1KRU@MYbw}Q zU2>~K09@0XHDFCVYLKUE%@50bqOG_0i{hCP_2M5joNFqt%_2UmZipO@%1_GRgf;#H zZvNP8Fi>|w#T3u1;pt>8#va$b$pg!b)4sO%+?#Amu(?S+(#**55#L`P=VS2W?I8Vz z&AQnbM+t5IqCqEcC%Y|0LFlG2#7mLOhEgypV#qXD^0&~_AR!%8$l0oHOm?_1)Ud~G zXQmOK&KejBSPj=_jf*hIKO{paB_-NWxWy%R*m&0)oRRq!D1q>V(YH@g#Z1o@BFQ@7 z1t=X|2gweE;8>Pr6QLHDXNSTIVamAduFa$z=55`#4>S2k4YcIb@#F7-PQkdFroEW$ z2h(#5tUb$xY50XLeYN1*Y>{VoHNw>al?_xQilX6PIHwQ`l{KoyYyK8)m!gxWw8&0GIsVH32?+V8Ab$Eh3>6erS5)=dwi!Mq;k*&SqI z*5ghipa;OAUJ?TOxG{ZkimpISD|+6Y4V9KUH{NNmC~~O-X+zT7uye0tk~shQ(z%Cn z?b3mOgM&u;b+@hq@Vlim+?c_n8oVg+*QXU;&ZwKdD!BTwVI~atBrgY8XwV^({dWHp zaG!lSiPTh$&epy?p2;6xMhAVUPf$rV#c*Z5K$jLDy1rVYasX3KUcU`oE0RKhP&BNyVVTy2|lw3Y042ht}xT9%T?x*m z2tsC;&OBY7Kay%<+bUoV6Fc+BWtY#@!}uI}W2R|HUCD5uy~msJ^vt%JGEj(Rv&W91 zZ^uHfQ;Rvs^$E}}k$*s{u6p&UTY!r-ku{-H_bOLI<{qtQgIK{}fI1`8yi$CyFk7Ys z0@4LC*5 zdEi{mr}MGeRSAHGv{0xyi~v%Fr~P!Fp_i@nXp!a5Q=_bQy^G9Lr^-9yKDZCDmepC- z_;N?BN|2LpWZKqb4U>kVC@GjL@%%(K!hyA;%qq~z;!tebrt4N`E=xtO*reZyo6^iIB1+lQl&3Nz@8lu>3PN5Ep?eL)M)(Lbh5-;qpllV>9x}jF;hY zk5tV8!gb1l=LhmW)<_oot^Qek1~Y28rh==z&RRHP-?vK$rW>AV33KXp?_nrd|9WN* z(3?#~C}n`Rg8TD@l&(WB+v+!< zroAT7`YTJ%yG_p>ncuL{+O^P#*4h=F%CU+2oxHwB`T9y{=yRD3Bid@gy9J{jhEVys zfL_kP@DWp42bwIK{T_9)_HP^Wm|tU&;!qQ1Xv{s&B!IiTKKUTzaBJM*WOv82F{o3? zcynl&r(* zx{4A!+TP^IPsjH~%?N;Xm7_^+PmWzLj&GI+SfjbQRC>|#@$CA$mjP1Uj>dw9w-ov# zSyZ0)wZ|}c_+Nyphndw8v!g~)+iueSX)`(w2;$z+aj_6KA0P4y`+wBx&bzQS36Qs8oiN<6@O z97OeP2o|2#u4Ww9%j5pFoThccPs8s(SYKAfg-KJlGt5D=Um)|*vnBl~lrC#I27Bbv zyQh|RLN>(om~HO~fY~Y^_cSlE$zD8? zf5&~qJILFLgk<##jP&_wkZBJP*49{)J(N=%cg7AdqD4)4w6|c`CMeT^VfQt`_sT8ZJpSH8U_v>|sjf?T zX@*VpLd-Sx&!d>d*>h94`lR8&)%tXsJFdRAvzcEXcmT*4C=jTVc_8n<2r9|uY#{PY zNK7uhaVc0KA_NS381#QVpG!-(fdfeU5g6z3>1k`&X-ki06x8>TSZ!IP*JRPpP31K* zb)N=6Cm)hE;X=a^ibV-yH?})Md!tN;Kab^?Y`y_?aKd2vC-sI0Er~Yw2p1-eF<{R; z!H`%@{mL}(u9AI~6u{^s%(nNJjm1DF6vM@}AXuA72dK$*6Sks2)e|dml`BJ*;&^QW z`C{HmeG=;>#~d9;05y8>4A58V!iXxoh&R+TFd@y)T<`9UnyaDz4E1@ge^5gx$thKk zpe**8JIDrPFFvH3b{-Hnxl^){q1&zD=U!kzPQbaQoPt{DjbrsZKXF_?XSL0Z5$+Hs z>?9Lt@3Ii-3VNnTkf#E?ToZ@aZ%2r%G&IIBl2rEVIJ|SU{{H1sjnmj%qeIKIZ1V%I z(*?&`!L-H3#7_9qzN>TIcb93-N(3!_YLtsju0JDL* zHbKW}Mf%;9rNNGAKsU3aso=p0JMtY+`gy}?d!e=|@TKK&gXRHBFz8DR$mReYJk(zj z4P)N7W6yAv60x9T|5)qh6x5@#wmpV$sS6~Em(d!G;WKz~@@muTlNEm&L*){|2}@>P zAs0&(Iz6~&s9d#hZqh@ANdPA?faHBlG~%bObG9R(he1tP6U5(8cduKEXktt(=<=7s zKinu7gPRy*)1Z!%XW`B0^I1sG7TgTvU7Prv?glztDOTy7eQ5~cMxbg5?XOOz1hhD@ zI?Gad^m58q2P~Vl(L6NnxnB|?tZZ3~MIvnsM`F6(j5k>33BEjc@^6z!;V(nJI2YlV zJ0GEOhfFR$`}3t|3P9V#Y$*seRszRwoXV{j1NOaKcdNw@1?xgi=hjTPZ}zk;%5mQ@efqf5 zj0}{i!R*e0@flDF7f@#2HX5`QW*H5THMP(sk{$6Gd6*8j@_IsTGtDZGZ=u=jAGXJ z4JD|j4Ie0086F%`yGG#^QYZhr+r~FTi9k{J;E;qySkCoT;aX=8F{lWwxayQ|b^7#}2elqjd&NNKp$wMZQVqegd6M^f$Yw(2KU1OlPl{!CP$zQ3j7NL3G7R7j0pD|>E{)F(`o;6^B!>Xf%yXNqRn{Vdb+%G+yVQ6FMwp`!@~)3^c{_Esb0G(gKN zu0g|-FzyI&Q6Q2~2;ioj(Vj1M;HaPM=31JI{xjF|ECR4<2KYYRdx_nuROu5hRwCs1+P0;aj)%3J;O1u-1#1(x#(OF!@*a!nXLg_*2{ZZ z@@3Hy5CjW)+_GX`M2LN9FJh%0AV9Dc*ry>xhqoY~T#u2y!^I+q+E6y&lVT8>Ak+m= z-`I`IJi(77pdSP{)yXX&ngJZFKFwjeuNDINz`r)2RxP_hgz2XXI^_w!LA_~vpvoaS z)AEiw6#QLQIkYoy)(EB+Ba4t($z&}V| zUVHWs1=9{$NbZd6HU7P`f~l*)HT(>8Aj7u5IxFy5PKJ!%k+lFgM_UMw1K){nn;T5= zJg{n@U>D5Tb#H8o1Ph?-Q^y5NQe7U46m8GNY3>8pbh;OH1>5l0o!Dlp;!*GwNQL67 zx37R3-5AYp!6H~wHB+q&#D5JB?|DQ&_vT*jBq|Xi8FqT#jtc$9?62`sr{KwjWEBJXNn9*XMq_k=Jxq53Z_zG+u5Z?Zj!{U1}^@0BE@ z3I~T5*K4wh9=lY@HJz;-ls=0Vg!BKuX1Ig@NJ32Ve^W3uJ&58BeZMmz{+~17;7LC{ zyV|&RM!Wz6cBBE#qI&~irxCtHKscr$rkk3L$4_*q*Z_<*Q9(8 z>o4VPk*eg5MJPCReU-59;K<$J#p*Y~sqnHs_)>nyGOa1eNp|}Cb$hIg_vehNdTHOW!i> z!!}Z6n{K97ZMiljx#}%?@a^aJtgnw>OMcS^_l`#7Irg|IICfFNjMx&p;Iv1}tyc|8 z@-r0%{Psl%$D30~5=jGiC5Ld_Eh;@tY%-fjK#&c-17OO?90bx#q2l z-fCQ=tr%t|`9Vd_M7*j;ZGGEwgd;vDTF3Hkahhanyw6iVOnfRde@3zoTLwEB&|y8w zi4~SgU_O258Tn{Ec}2zcpk9h=ahj-ean5lO@uk&f9{kiDp{2(i$IE^nPznmp!*nNN z5wRlRilxl=;1QtMr4T7nTsG+bYOqOoL$gmrq{O)ce5@Q)O`h8s!*7&a0&nMnD z!jB~QROe)9z* zqq^+*m`z2k%@yBzFr`;Ict+MdQP-s{C^$T4{Dy1x;BXXvP<&(cbchCbEL&}T4l~s& z1ApExOUTC7;or=ZR2_T;SM}HnrDw7+d=(X!`sDEZMdT4Vg}LVg&(~fo-3S;rWhGvf z@UF8#HD}tg9(U`k1Ua}bbO3_wt85b^b68Y zIb)i_f=$csxn5zKlKr)8uys%`P6g9WZAIxb#sRN?*Fn0}-jUA=&oM6|kMsHT*Wve> zIoFu=onRMiI8)X)(uzm2AE>+>0|w;E2QOH-=8JFe2Mk!Q+vr5=1o%$52WF5SF8!MC z*r1)~I^sIVtz(5nXgP7Gou)Y~vLwt?MbhqJIUIs7aKDhTvHiZ_X={8?2iXlbPNn5! zBq}|-WrwZ0L=!*DYvsDWZ`6mBBY7?0>Ksj1rmyVU?*kUHgTs`;VUEFJ=i~;J6p^&& zYa2}rB9}imhIh8pG8ORNgtk>sk~30FmWh9JHm7E(u^ApN-ZW+vDVG#|ctOGBf + Wave equation <_rst/tutorial3/tutorial.rst> .. ........................................................................................ diff --git a/tutorials/README.md b/tutorials/README.md index 3abeb6c..d1b09e3 100644 --- a/tutorials/README.md +++ b/tutorials/README.md @@ -7,5 +7,6 @@ In this folder we collect useful tutorials in order to understand the principles |-------|---------------|-------------------| | Tutorial1 [[.ipynb](), [.py](), [.html]()]| Coming soon | | | Tutorial2 [[.ipynb](tutorial2/tutorial.ipynb), [.py](tutorial2/tutorial.py), [.html](http://mathlab.github.io/PINA/_rst/tutorial1/tutorial-1.html)]| Poisson problem on regular domain using extra features | `SpatialProblem` | +| Tutorial3 [[.ipynb](tutorial3/tutorial.ipynb), [.py](tutorial3/tutorial.py), [.html]()]| Wave problem on regular domain using custom pytorch networks. | `SpatialProblem`, `TimeDependentProblem` | diff --git a/tutorials/tutorial3/tutorial.ipynb b/tutorials/tutorial3/tutorial.ipynb new file mode 100644 index 0000000..06a7d5e --- /dev/null +++ b/tutorials/tutorial3/tutorial.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "# Tutorial 3: resolution of wave equation with custom Network" + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "### The problem solution " + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "In this tutorial we present how to solve the wave equation using the `SpatialProblem` and `TimeDependentProblem` class, and the `Network` class for building custom **torch** networks.\n", + "\n", + "The problem is written in the following form:\n", + "\n", + "\\begin{equation}\n", + "\\begin{cases}\n", + "\\Delta u(x,y,t) = \\frac{\\partial^2}{\\partial t^2} u(x,y,t) \\quad \\text{in } D, \\\\\\\\\n", + "u(x, y, t=0) = \\sin(\\pi x)\\sin(\\pi y)\\cos(\\sqrt{2}\\pi), \\\\\\\\\n", + "u(x, y, t) = 0 \\quad \\text{on } \\Gamma_1 \\cup \\Gamma_2 \\cup \\Gamma_3 \\cup \\Gamma_4,\n", + "\\end{cases}\n", + "\\end{equation}\n", + "\n", + "where $D$ is a square domain $[0,1]^2$, and $\\Gamma_i$, with $i=1,...,4$, are the boundaries of the square, and the velocity in the standard wave equation is fixed to one." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "First of all, some useful imports." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 2, + "source": [ + "import torch\n", + "\n", + "from pina.problem import SpatialProblem, TimeDependentProblem\n", + "from pina.operators import nabla, grad\n", + "from pina.model import Network\n", + "from pina import Condition, Span, PINN, Plotter" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "Now, the wave problem is written in PINA code as a class, inheriting from `SpatialProblem` and `TimeDependentProblem` since we deal with spatial, and time dependent variables. The equations are written as `conditions` that should be satisfied in the corresponding domains. `truth_solution` is the exact solution which will be compared with the predicted one." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 3, + "source": [ + "class Wave(TimeDependentProblem, SpatialProblem):\n", + " output_variables = ['u']\n", + " spatial_domain = Span({'x': [0, 1], 'y': [0, 1]})\n", + " temporal_domain = Span({'t': [0, 1]})\n", + "\n", + " def wave_equation(input_, output_):\n", + " u_t = grad(output_, input_, components=['u'], d=['t'])\n", + " u_tt = grad(u_t, input_, components=['dudt'], d=['t'])\n", + " nabla_u = nabla(output_, input_, components=['u'], d=['x', 'y'])\n", + " return nabla_u - u_tt\n", + "\n", + " def nil_dirichlet(input_, output_):\n", + " value = 0.0\n", + " return output_.extract(['u']) - value\n", + "\n", + " def initial_condition(input_, output_):\n", + " u_expected = (torch.sin(torch.pi*input_.extract(['x'])) *\n", + " torch.sin(torch.pi*input_.extract(['y'])))\n", + " return output_.extract(['u']) - u_expected\n", + "\n", + " conditions = {\n", + " 'gamma1': Condition(Span({'x': [0, 1], 'y': 1, 't': [0, 1]}), nil_dirichlet),\n", + " 'gamma2': Condition(Span({'x': [0, 1], 'y': 0, 't': [0, 1]}), nil_dirichlet),\n", + " 'gamma3': Condition(Span({'x': 1, 'y': [0, 1], 't': [0, 1]}), nil_dirichlet),\n", + " 'gamma4': Condition(Span({'x': 0, 'y': [0, 1], 't': [0, 1]}), nil_dirichlet),\n", + " 't0': Condition(Span({'x': [0, 1], 'y': [0, 1], 't': 0}), initial_condition),\n", + " 'D': Condition(Span({'x': [0, 1], 'y': [0, 1], 't': [0, 1]}), wave_equation),\n", + " }\n", + "\n", + " def wave_sol(self, pts):\n", + " return (torch.sin(torch.pi*pts.extract(['x'])) *\n", + " torch.sin(torch.pi*pts.extract(['y'])) *\n", + " torch.cos(torch.sqrt(torch.tensor(2.))*torch.pi*pts.extract(['t'])))\n", + "\n", + " truth_solution = wave_sol\n", + "\n", + "problem = Wave()" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "After the problem, a **torch** model is needed to solve the PINN. With the `Network` class the users can convert any **torch** model in a **PINA** model which uses label tensors with a single line of code. We will write a simple residual network using linear layers. Here we implement a simple residual network composed by linear torch layers.\n", + "\n", + "This neural network takes as input the coordinates (in this case $x$, $y$ and $t$) and provides the unkwown field of the Wave problem. The residual of the equations are evaluated at several sampling points (which the user can manipulate using the method `span_pts`) and the loss minimized by the neural network is the sum of the residuals." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 4, + "source": [ + "class TorchNet(torch.nn.Module):\n", + " \n", + " def __init__(self):\n", + " super().__init__()\n", + " \n", + " self.residual = torch.nn.Sequential(torch.nn.Linear(3, 24),\n", + " torch.nn.Tanh(),\n", + " torch.nn.Linear(24, 3))\n", + " \n", + " self.mlp = torch.nn.Sequential(torch.nn.Linear(3, 64),\n", + " torch.nn.Tanh(),\n", + " torch.nn.Linear(64, 1))\n", + " def forward(self, x):\n", + " residual_x = self.residual(x)\n", + " return self.mlp(x + residual_x)\n", + "\n", + "# model definition\n", + "model = Network(model = TorchNet(),\n", + " input_variables=problem.input_variables,\n", + " output_variables=problem.output_variables,\n", + " extra_features=None)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "In this tutorial, the neural network is trained for 2000 epochs with a learning rate of 0.001. These parameters can be modified as desired.\n", + "We highlight that the generation of the sampling points and the train is here encapsulated within the function `generate_samples_and_train`, but only for saving some lines of code in the next cells; that function is not mandatory in the **PINA** framework. The training takes approximately one minute." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 5, + "source": [ + "def generate_samples_and_train(model, problem):\n", + " # generate pinn object\n", + " pinn = PINN(problem, model, lr=0.001)\n", + "\n", + " pinn.span_pts(1000, 'random', locations=['D','t0', 'gamma1', 'gamma2', 'gamma3', 'gamma4'])\n", + " pinn.train(1500, 150)\n", + " return pinn\n", + "\n", + "\n", + "pinn = generate_samples_and_train(model, problem)" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00000] 4.567502e-01 2.847714e-02 1.962997e-02 9.094939e-03 1.247287e-02 3.838658e-01 3.209481e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00001] 4.184132e-01 1.914901e-02 2.436301e-02 8.384322e-03 1.077990e-02 3.530422e-01 2.694697e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00150] 1.694410e-01 9.840883e-03 1.117415e-02 1.140828e-02 1.003646e-02 1.260622e-01 9.190784e-04 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00300] 1.666860e-01 9.847926e-03 1.122043e-02 1.142906e-02 9.706282e-03 1.237589e-01 7.233715e-04 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00450] 1.564735e-01 8.579318e-03 1.203290e-02 1.264551e-02 8.249855e-03 1.136869e-01 1.279038e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00600] 1.281068e-01 5.976059e-03 1.463099e-02 1.191054e-02 7.087692e-03 8.658079e-02 1.920737e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00750] 7.482838e-02 5.880896e-03 1.912235e-02 5.754319e-03 4.252454e-03 3.697925e-02 2.839110e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 00900] 3.109156e-02 2.877797e-03 5.560369e-03 3.611543e-03 3.818088e-03 1.117986e-02 4.043903e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 01050] 1.969596e-02 2.598281e-03 3.658714e-03 3.426491e-03 3.696677e-03 4.037755e-03 2.278043e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 01200] 1.625224e-02 2.496960e-03 3.069649e-03 3.198287e-03 3.420298e-03 2.728654e-03 1.338392e-03 \n", + " sum gamma1nil_di gamma2nil_di gamma3nil_di gamma4nil_di t0initial_co Dwave_equati \n", + "[epoch 01350] 1.430180e-02 2.350929e-03 2.700139e-03 2.961276e-03 3.141905e-03 2.189825e-03 9.577314e-04 \n", + "[epoch 01500] 1.293717e-02 2.182199e-03 2.440975e-03 2.706538e-03 2.904802e-03 1.891113e-03 8.115429e-04 \n" + ] + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "After the training is completed one can now plot some results using the `Plotter` class of **PINA**." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 11, + "source": [ + "plotter = Plotter()\n", + "\n", + "# plotting at fixed time t = 0.6\n", + "plotter.plot(pinn, fixed_variables={'t': 0.6})\n" + ], + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7YAAAFpCAYAAACoDbRLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABP7UlEQVR4nO3df/QkdX3n+9ebGQdxQRzAIDuATBBng8SLcRYNuWaNIKJJwHNisrCR4D0Y4yZustc1J3i8R3NRz9F4jbvummQnSkROFAj54WyCAUSN3hCUSUQZ8E4YwMCMKAG+jhCBEXjfP7oaanq6v13d9an6/Kjn45zv+X67u7r70/XtftfnVZ9PVZu7CwAAAACAXB0QuwEAAAAAALRBsAUAAAAAZI1gCwAAAADIGsEWAAAAAJA1gi0AAAAAIGsEWwAAAABA1uYGWzO72MzuNbPtM243M/uwme00s6+b2Y+FbyaAnJjZmWa2o6oLF065/SfN7B/M7DEze93EbY+b2U3Vz9Ye20ytA7CQHGtd9dzUOwCNNah1B5rZ5dXtXzaz46rrjzOzh2u17g9q93mxmd1c3efDZmZt29lkxPbjks5c5fZXSzqh+nmTpN9v2ygA+TKzNZI+olFtOFHSuWZ24sRid0l6g6RPTnmIh9395OrnrE4bu6+Pi1oHoKGMa51EvQPQUMNad4GkFXd/nqQPSXp/7bbba7XuzbXrf1/SL+upWrNaTWpkbrB19y9KemCVRc6W9AkfuUHSs8zsqLYNA5CtUyTtdPc73H2vpMs0qhNPcvdvuvvXJT0Ro4HTUOsALCjLWidR7wAsZG6tqy5fUv19paTTVhuBrerJM939Bnd3SZ+Q9Nq2DQ1xjO0GSXfXLu+qrgMwTG1rwtPNbJuZ3WBmrw3asnaodQDqSq11EvUOwFOa1IMnl3H3xyTtkXR4ddtGM/uqmf2Nmb2stvyuOY+5sLVtH2ARZvYmjaa06KBn2Is3Ht/r06Mw9z9+cCePe/iahzp53C7cevMP7nP3Zy9yn594+dP9uw+0Gzy49eYf3CLpkdpVW9x9S6sHfcpz3X23mf2wpM+Z2c3ufnugx+5FvdY97aA1Lz5s4zMjtwjI23duXaHWJahe6+zpT3vx048+InKL+vX4o2tiN6GVNQc+HrsJRVvm/bH37l2l1bp7JB3r7veb2Ysl/YWZvSDA404VIlnulnRM7fLR1XX7qVbQFkl6wQvX+Sf/8sgAT4+SXLpyauwmzHTe+utjN2E/Jz931z8tep/vPvCE2n72Tn7urkfcffOMmxvXhGncfXf1+w4z+4KkF0lKobO3VK17zgsO89d/8pXdtw4o2AdPvoJa169Gr61e655xwr/25/3uG/tpXUQP3Xlo7CZ05uCNe2I3oSjLvFfu/I235VrrxsvsMrO1kg6VdH81zfhRSXL3vzez2yU9v1r+6DmPubAQwXarpLeY2WWSXiJpj7vfE+BxUZiUQ2sTk+1PMegm4kZJJ5jZRo2K1DmS/kOTO5rZeknfd/dHzewIST8h6Xc6a+liqHUA6kqtdRL1bh8lh9m6yddJ0EWlSa3bKul8SX8n6XWSPufubmbPlvSAuz9ezU45QdId7v6AmX3PzF4q6cuSfknSf2/b0LnB1sw+Jenlko4ws12S3iXpaZLk7n8g6SpJr5G0U9L3Jf0fbRuF/OQeWpdRf82E3Ke4+2Nm9hZJV0taI+lid7/FzC6StM3dt5rZv5X055LWS/pZM/u/3f0Fkn5E0v80syc0OgfA+9z91j7aTa0DsIhca51EvWtqKIF2FoIupGa1TtLHJF1qZjs1OjHdOdXdf1LSRWb2A41Oovdmdx+fuO5XNTpD+0GSPlP9tDI32Lr7uXNud0m/1rYhIcQIV6UHmiEG1mWM11Pp74em3P0qjTpG9eveWfv7Ru07BWV8/fWSfrTzBk6RU60DkIYca131/NS7GYYeZldTXzeE3GFpUOsekfTzU+73p5L+dMZjbpN0Ush2Zn32phRCV6g2dBWIUlhHQ8IoLgAA+SHQLobRXKQoy2BbYlgr8TUN3aUrpxJuAQBIHKG2PYIuUpBVsCX8ITdMUQYAIE0E2u4wbRkxZBFsCbTIHaO3AACkIcVAe8jtBwR5nAePb/d9pl1gNBd9ST7YEmpRCkZvAQCIK4VQGyrELvLYKQVegi66knSwJdSiRIzeAgDQr5iBtssg26YNqYTdaf8bwi6WkWywJdSiZIRbAAD6ESvUphBoVzPZvlSCrsSoLpaTZLAl1GIICLcAAHSr71CbephdTb3tKYVciVFdNJNksAWGgnALAEB4BNp2Uh7NHZv1PybwDldywZbRWgwN4RYAgHD6DLWlBdpZUh7NnUTgHa6kgi2hFkNFuAUAoL2+Qu1QAu00OYzmTjPvvUHwzV9SwRYYMsItAADLI9TGkdNo7mqavH8Iv2lLJtgyWgsAAIBl9BFqCbTzlRJyZ1nkfUYI7l8ywRYAo7ZAU9fctSl2E6Y649gdsZsAoAOE2sWVHnKRHoItkBjCLYYu1dDaRNO2E4CBcLoerSXUtkfIRR+SCLZMQwaAYck5vIYw6/UTeIHFdBlqCbTdIOSiK0kEWwD7YtQWJRl6iF3EtHVF2AWmKz3Urt+xt/VjrGxaF6Al3SHkIiSCLQAgKIJsWJPrk6ALdCtWqA0RZJs+ZoqBl5CLtqIHW6YhA9MxaoucEGb7Q9AFuhut7TvUdhFml3ne1IIuIRfLiB5sAQB5Isymof5/IORiCHIPtbHC7GpSDrrj/wsBF/MQbIGEMWqL1BBm0zb+/xBwgcX0EWpTDLSz1NuaSshlFBfzEGwBAHMRaPPCKC5K1cVoLaF2dSmHXAIu6qIGW46vBeZj1BYxEWjzxyguMFvXoTbnQDtNaiGXgIs6RmwHKJWOKp0sIF2p1AmEQ8BF7rr8ep8ulBZqJ6UUcgm4kAi2RcmtI7pae+l4AXHkVkewuGvu2kSNBdTdaG3pgXaaVEIuAXfYCLYZGVKHc9prHXJHjOnI6NqQ6gsYvUV+Qo/WEmq7M14HBFz0jWCbIDqY0/HdjUA3qDnDRcDFEBFq+5HCKC4Bd1gItgmgU7kczvoJtEf9gcT0ZKQth2NrCbWriz2Ke8jtBxBuB4BgGwEdyfCGEHKZjoyQqEOYRLjFEHQxWkuobS7mKG7fo7cP3XmoDt64p5fnwki0YHv/4wfHeure0YHsF1PrgNVRkzAL4RbLeua6R2I3YS5CbVpijeIyelsuRmw7QscxPgIusD9qE+Yh3CIlKU9DJtSGESPgEm7LRLANjE5jekoKuExHRhtDqU9dd4SHMLWMcIvShB6tJdSG13fA5cRS5SHYBjCUzmLuSgq4wKJKq1MxR3FWe+6SQi/hFos649gdQWtNqM85oTYvMQIu4bYMBNsWSusoDgWdNQxN7rUq5amIk6a1NeewS70E9kWo7U+fAZdwWwaC7ZJy7ygOHaO3GIoca1VOQbaJydeTW9Al3CJnXX1nLfqzfsfeaF8ThLwQbBeUYycRs9FhQ8lyqlelhdnV1F9rLiGXWom+pVgTUhutPWj77taP8fBJGwK0pHt9jN4yaps/gm1DOXUQsZjcOmycQAqlSLHj2recQm5utRIIOVqbSqgNEWZXe7zUg27XAZdwmzeCbQOE2vIxNRmlSbluEWinyynkAuhX6EDb9HlSDbpdTk8m3OaLAw/mSLlziPD4f4dhZmea2Q4z22lmF065/UAzu7y6/ctmdlzttrdX1+8ws1f12vBCpPg+fujOQ5/8wXyprqsU31sxUevSVcpo7UHbd/cWald7/phtmCWVUfQhWLbWmdkrzezvzezm6vcravf5QvWYN1U/P9S2nQTbGa65axMb8IHi/96Oma2R9BFJr5Z0oqRzzezEicUukLTi7s+T9CFJ76/ue6KkcyS9QNKZkn6vejxkKtWAlosU1x81coRa152U3vOxwlOKYTLFkLt+x95O/kecdOwpbWqdpPsk/ay7/6ik8yVdOnG/X3T3k6ufe9u2lf/aFGy0wXuglVMk7XT3O9x9r6TLJJ09sczZki6p/r5S0mlmZtX1l7n7o+5+p6Sd1eOhoVTeuykGspyxPpNErUtU7qEkpeA4S2ohl3DbqaVrnbt/1d2/VV1/i6SDzOzArhrKMbYTUukUIr5ST5Ry/+MH69KVU1s+yhVHmNm22hVb3H1L9fcGSXfXbtsl6SUTD/DkMu7+mJntkXR4df0NE/dN8wCfBKVSvwhg3Rmv29jH4OZQH6l1aCvGaG0qQXER4zbHPh53qF8LlHitu6+2zM9J+gd3f7R23R+Z2eOS/lTSe9zd27wKgm1NKp1CpCOHzlsk97n75tiNQFoItP156M5DCbf9oNb1bMh1JMdQW5dCwA0dbtueSCqFWt1Qp7XOzF6g0fTkM2pX/6K77zazQzQKtudJ+kSb52GMvUKoRU7a75nr1G5Jx9QuH11dN3UZM1sr6VBJ9ze8L6aIWcOYJhsH6z06al2CQk0f7Xu0NvdQWxd7ijInlQquTa2TmR0t6c8l/ZK73z6+g7vvrn4/KOmTCnA4BsFWhFqsjvfHwm6UdIKZbTSzdRqdIGXrxDJbNTqJgCS9TtLnquknWyWdU51db6OkEyR9pad2YwkEq/hi/g8GXh+pdQiipFBbFzPghgy3HGu7fK0zs2dJ+itJF7r7344XNrO1ZnZE9ffTJP2MpO1tG8pU5IHoo+OTyVSLpQxkyl0Q1bEVb5F0taQ1ki5291vM7CJJ29x9q6SPSbrUzHZKekCjIqlquSsk3SrpMUm/5u6PR3khGYkVLgi16Ujl2NshodalJ9fR2tLFmqI81GNuQ2tT6yS9RdLzJL3TzN5ZXXeGpH+RdHUVatdI+qykP2zb1sEH29L2NsfsaM567lI6WoTb5tz9KklXTVz3ztrfj0j6+Rn3fa+k93baQLRCoE1XjOO5hlwbqXVhDbG2lDpaO02MgBsq3LY91jZ3y9Y6d3+PpPfMeNgXh2yjNPBgW0KozWEjMNnGUoIukIK+61gONWfoMjpZCZCcPkdrhxRq6/oOuIzcDsdgJ43nGmrHJwvJ+aQhObc/1/cNEEKOn9mh4n+FoeE4yPz0GexD7LDgPZY+/kMZyDkIzpPjayPcYohy+oxipM//GXURwDJin0EZZRnkVOQcNsBD7ERy8hNgMX3VslLq0aJ720s4noppychFCnWGacjxHLR9d+dTk5mSXL5BBtuUpVDYY8sh4A75ZCkYltxqUsipYqs9Vk6ht69wS11ELEwRLUMfx962DbdDP4lU6gYXbFMdrc2t89iH1AMunTiULoe6FKtDO+15U+7sMHILIBddj94ycluuQQXbFENtDh3H2FIPuEAMXdezVGtTyiMz9balGHIJt8BsfHdtWvqYmtwXam9/0u0hFC63EyalIMV1luLOEqA0h9x+wJM/uUi1zV3XUGoi+pbaZwzhdHliqTY7MnjPpWsw/5mUNraphbPcsP5GLl05NXYTUKhUPmMpBsNlpPY6Uvn/AmO8J7GaFMMt0tRoS2tmZ5rZDjPbaWYXTrn9WDP7vJl91cy+bmavCd/U/KU44pirlNZlSjtN0E4uta7L91wKn6vUgmAoqY7iYnhyqXV9IeCkj7NIo4m5W1czWyPpI5JeLelESeea2YkTi/1fkq5w9xdJOkfS74VuaO5S6CyWiPWKUKh18T9PQwp9sV9rl/9rdvalraRaN5R6gZEuwu2yOzV476WpyX/lFEk73f0Od98r6TJJZ08s45KeWf19qKRvhWtiezE3simNLJaK9YtAsq91uYod8mIqNdwiadQ6ZIuRW6ymyRZ1g6S7a5d3VdfV/bak15vZLklXSfpPQVqXOToN/Ym9rhmhKEIWta6r91qMz9CQA20d6wE9S6rWxd5+Iz+hwy1T0csRakt6rqSPu/vRkl4j6VIz2++xzexNZrbNzLZ9f+XRQE+9uliBg0LdP0bH0YNka10bsUIt9hUj4FIzMUORtW6aWKGmlK+yiYWRW0zTZAu6W9IxtctHV9fVXSDpCkly97+T9HRJR0w+kLtvcffN7r75GesPXK7FGaCjEFes9c+obfaodT1gdHK+EsIt9TBpRdQ66ghChttldnDwHkxPk//IjZJOMLONZrZOo5MIbJ1Y5i5Jp0mSmf2IRgXwn0M2NBeE2jTwf8ASBlnr+vys0AlojnWFDnVW69ihsRhGbYGw5m453f0xSW+RdLWkb2h0lrxbzOwiMzurWuy/SPplM/uapE9JeoO7e1eNbqrvAkuYSgv/Dywih1qXc6eRoLa4Pke3qZfDkUOtA5qKPWrbFDW2H2ubLOTuV2l08oD6de+s/X2rpJ8I27S88IZN00N3HqqDN+7p7fmuuWuTzjh2R2/Ph7CGVuv6qluE2nYOuf0APXj8E7GbgYKkUuuWrUEl1ZSHT9rA8aItHbR9d7TRb+pzWsqpDBERagFgfxxPG04f6zH0tiznGQZAn5iS3F6onQOcITlvjUZsc9TXBjW3UBuic5Tbnqm+R22BHHRduwi04TEyAISVUohh5BZoj55HC6mH2vFoSf0n5cftUp//K0Yp0IWc3lc51IRcdb1uU9+uASVj5LYddgwgWu/j8DUPxXrqIFLc+McMmzmE3BT/Z0AMXX4WUq4BpWAdI0Xf2/v0Xp6n9Pc/4badEOF20ZH80t+TOSnyP5HTyEYIqQXKHEIugPD4zPeny3XNTkD0gffZbA+ftIGA2wIjt8NFL2QJKRTjXMJjam3s6383tJ0ryEcK9QthpFRbZ6EWAssj3MaT0vHXaC79rWJiYncKUwuKTaUUxGP/D4ESpfDZHqKu1jt1EkgDo7fLYdR2mIo9K3JpSuo0cmZPYDGhRr26Ciu51KdF98CvbFrXUUsAYDGcNTl/fEtH94oLtl1Oe4qxBzuXDuOixq8rVsCluABhpFqjQkwjm/YYKYZddhZiKFKtN30Zj9wScJs5aPvuVqPd63fsbVzzqcNpGHaFWAChthsxpyd3/T/l2DKULrUatX7H3id/+niOlI7B6uJ/wXRkdIX3VjtMT26OnQDDUtyIbQlS6yz2gT1dQLdK7UjGDpfj509hJJc6CgwLI7jAvoaXoJbQZ4dwiKF2LMbobamdfaBrsWtVaiOmKY7ipoLZK0C3GL1dXZvgT03Py3BTVIJidxRTwXoA0hbzM5pDeIzZxtD/G3b+IRX0DVbH9OT4eI/GV9R/oIu9wn1s1FP5GpyU9Lk+6LghZSHqWgnv8RwC7aRSwi0QUgn1KGUE3OmYrj0MbP0iowMyWwnrhil4KE3fn8scA21d7u0HSpbCsfFdIeCGQw3PR/7JoUNd71UsIbh1ra91xB5kID0ldSb6Drghayf1EcgX4fYpKYzaUk+7RbKKhFDbHOsKWF7IjWifn8WSQm1dqa8LQLoYvcVQkBhm6HKPCkFtcX2ss9z2op23/vrYTQCCG8LU3b5eX0rbGg7LQBtdvpdLno48iXC7/Kht07qdUt0domLWfi4bTd7wy2PdAfH08fkrPdDW5fZac9vxhzTxPoqP0VuUjKQwRVeFl2DWXo7rMJedLn0xs8PM7Fozu636vX7Gcn9tZt81s7+cuP7jZnanmd1U/ZzcS8MzFKqWEWq70cdrzrFmloJah5QNOdymcKxtbszsTDPbYWY7zezCKbcfaGaXV7d/2cyOq9329ur6HWb2qqaPuQy2eD2hcxFOl+uSvcm9uFDSde5+gqTrqsvTfEDSeTNu+013P7n6uamDNqJHQwy1Y0N+7QNArcvEkKYj1zF6u5ih1mszWyPpI5JeLelESeea2YkTi10gacXdnyfpQ5LeX933REnnSHqBpDMl/Z6ZrWn4mAsjbQHo29mSLqn+vkTSa6ct5O7XSXqwpzYlaQij/UPtKNR1vQ7YsRoNtS4jQw230rBHb9HIKZJ2uvsd7r5X0mUa1be6er27UtJpZmbV9Ze5+6PufqekndXjNXnMhUXd2qV48psuRuzoVITHqG3WjnT3e6q/vy3pyCUe471m9nUz+5CZHRiwbZjQ5WeNUPuUHNYFtXFh1LqaZd8/ffahCLfD0eV05AL7/Rsk3V27vKu6buoy7v6YpD2SDl/lvk0ec2Fr2z4AVpfqm7tpJyrlIn/I7QfoweOfiN2M7Hxv79NDjAQeYWbbape3uPuW8QUz+6yk50y53zvqF9zdzcwXfO63a9RJXCdpi6TfknTRgo9RvNRDSA5Brm/rd+ztrOYOsV5S67CMlU3rBlufHj5pA8ef9uChOw/VwRv3BHu8PmpdLgi2NaE7gimF2mWL9LT7pRR2h9hZS8R97r551o3ufvqs28zsO2Z2lLvfY2ZHSbp3kSeujYA8amZ/JOlti9wfzXVVw4baaWyiy3CLpVDrBmj8GRxirRqP3A4h4B60fffCI9UF1+jVat1uScfULh9dXTdtmV1mtlbSoZLun3PfeY+5sHSSVwtDOA5tUePvguziOyG7etyUpD7albmtks6v/j5f0qcXuXPVQVR17MZrJW0P2Th0q+S6EUpX6yilna0DQa3L3MqmdaWGmLmGNjUZq7pR0glmttHM1ml0MqitE8vU693rJH3O3b26/pzqrMkbJZ0g6SsNH3NhbOU6EqsD0XfgTCHg5tBZY+fLPt4n6ZVmdpuk06vLMrPNZvbR8UJm9iVJf6LRCQh21U4R/8dmdrOkmyUdIek9vbZ+IHL4XJUsdl2dpc1OvwHWQWpdIcYBd2ghl3AL6cljZt8i6WpJ35B0hbvfYmYXmdlZ1WIfk3S4me2U9FZVZ4F391skXSHpVkl/LenX3P3xWY/Ztq1MRa6EHKGL0SGM3QmqP3+Mws+U5Hy4+/2STpty/TZJb6xdftmM+7+iu9ahS7HrFNAnat1TcjhxVFOr9XFKrHGlH3e7zHTkJkrrl7r7VZKumrjunbW/H5H08zPu+15J723ymG0RbDOXYhEt5fiD0Af3A31JcSp9irUqdV3U0tI6W0BKmnxec6yFpYfbRZXSzy0RwTawvvYupl4Yx+3r84NPhw0II8VRkqGiAwWUJdcR35LDbVejtugfvRelObqxmpQL36Sc2gqkpKTjEakD7aS2/nLbZgK5qB/Lm+IOLcIfUkewDaiPUY7UOjhN9HmCKUaagLTkWLNKR51En9gRsrwUgy7hNgw+F91g65aJFM4+3FaO7afwYGgIPWnKsX4CbVCL9pdKyC0x3C46zZqanCaqRiBdFuCSPjx9vBY2hhiylHbGlFS7ACAlsUNuieE2NPqj/ct+jbc9Di2lTuA0JXYMS3xNAPJz0Pbd+/ykLmTtpMMFlCOFUVwgBZwVOYCuOgglB8Cuz/SZ4hmSr7lrk844dkfsZgAzhaxlqdWvJsF11jKMTADLSX3woDTjflVf9be0MyVzduT8scs2Ual1CruQy2tkwwzkKdRobEojuqnUTeoikK4+R3CHHARTqcd4CsG2pS5Ga4f0QenytTLVDogjdg3rKojmNG25iRg1sqSvsUJ32H6HQbjF0Ay6cqS4xzl2hzCGIb5moAttaloJHck+Q2fMcEvNBNBUX6O3pYTbUnZcDlX+PZmIQncEh9xZGfJrB0oS47McaxS1pNFbAGXjBFNxrJYVUhxgyx3BFsnookMcaucDxQdIUwrBMkYb2BmI1LHdTFOX4baUUdtFUIvTQrBNBB8MLOK89dfHbgKwn77rWAqhdizX0dsSpqADWAzhdnU51nKMRN+i5dpBL/lrMWJKedQWKFWOn5FUOx59tiv2toMROXQhx3qUI6Ymo0SDrR6pbJBjd0xSxDrB0OV45tg+P7ephtqx1NsHAGNdhNsSRm2Rp8EG2zbYmwgAceQSGnNpJ9CFVAYP0AzhFqUgoUXEyORsodcNOyNQuiF0JHMLi320N0StpD4CYFryvhap3/PqMDW2P1mv6Ryn6yFfQwgOwDL62EmXW6gdy7XdAIYndLhl1BZ9Wxu7AbkJtdeF0dr51u/Yyx5EoGM57EkmHALDkEM9WsaDxz/ReNnY62Bl0zr6qD166M5DdfDGPbGbUYwyK8gcsUfeKBjNsa6AYSsh1Hb9GqiTSEnsPlZsDx7/xH4/be+/6GO0FXJQIedR2xK2P0PDiO3ATfvQ5lyEVnPI7Qf0vnEAStdlqCqpU3HQ9t3F1lZgyPrqV9Sfp49R3ZAjtw+ftKGoeo50EWwXEKKQxNyz3rSoTC4XuzNW0pTka+7apDOO3RG7GQAiSDncLrvjj2l0GKLYO8n7CrlMS26upL5qzgY5FXloDtq+u9WesvH92dsGoC/Um8XQ+QS6F2Na8DxdtylUWEt1p15fYh87PRSs5R713fHoIozGCripdNqGfuwQytJ2Q5vK5zInBHaUrM02MtWOf6zjXBeVQxtzRM3OS5pVBK30ET5z/aCnuuEE8JRc6wua4av6kINcg2IX7WbUFrkYXC9/2b2JuYxs9NkhZHoykIaSZhIMoaZ09RoZQQfCyDHQTko13AJdGlywLVXMkNnX89JpA1CKFAM8M1rQxuOProndhNZyHaWdJfTrCRFuSx61pZ8aH1uxAqTQQUqhDQD608UGnDoCIIbSAu2kkl9bH7reNpU06yq2RsHWzM40sx1mttPMLpyxzC+Y2a1mdouZfTJsM/M2lO957KMtIdYloxKYhVoHYAiodU8ZSugL9ToZtUXK5vbwzWyNpI9IerWkEyWda2YnTixzgqS3S/oJd3+BpP8cvqnxpBqEUgq1Yym2CWhiaLUutbo2xNrRxWtmKhzmSbHWxahHpY/SThPqNXO87XJS2+6WqMkaPkXSTne/w933SrpM0tkTy/yypI+4+4okufu9YZuZr646GSl3AlNuG7AKah2AIRh8rRtaoJ009NePcjUJthsk3V27vKu6ru75kp5vZn9rZjeY2ZnTHsjM3mRm28xs28oDfKiWlUNw7LKNjEigI53Uuu+vPNpRc+PhMxhODvV8Ho4Py04nte7xh/6lo+aGRagbabse2o7aMh0ZXQg1Jr5W0gmSXi7pXEl/aGbPmlzI3be4+2Z337z+sP6H49n49quEDts0vI8GbeFa94z1B/bbwgyVWiuAjC1c69Yc/K/6beESCLX7ih1uc9J0O8WO37iapMvdko6pXT66uq5ul6St7v4Dd79T0j9qVBCz12Y+PGcNTRfHOWCKQdc6xENdR88GWesItdPFXC+M2iK0Jr37GyWdYGYbzWydpHMkbZ1Y5i802qsnMztCoyksd4RrJqQ8Oz9dtXnIe8TOW3997CaUKolad81dm0I+HAaoTX1kp98gJFHrxvp4zxFqV9dm/Qxp1BbpWztvAXd/zMzeIulqSWskXezut5jZRZK2ufvW6rYzzOxWSY9L+k13v7/Lhg+t85djqB07aPtu9soheanWutLlXNuAHA2t1qUQag/euKfV/TkECmhmbrCVJHe/StJVE9e9s/a3S3pr9QMNe0QRyBW1bj5qWzfYAYg+DaXWxQy1bcPsao/VRdB98Pgnlh49X9m0bultw8MnbchmByd1On3MOcpALh/41ZTwGhCGmR1mZtea2W3V7/VTlnmumf2Dmd1kZreY2Ztrt73YzG42s51m9mEzs35fQVqW6eAw3RToHrUurhih9uCNe5786eN5QkthdLt0s7a/OY/KN6l11XLnV8vcZmbnV9c9w8z+ysz+v6oGvq+2/BvM7J+r+niTmb1xXlvo3SBbjBxl60JJ17n7CZKuqy5PukfSj7v7yZJeIulCM/vX1W2/r9F3LJ5Q/Uz9Ggqkjx1eKBy1LpK+A1ofYbav51123bU51pZR0OzNrXVmdpikd2lU506R9K5aAP5/3P3fSHqRpJ8ws1fX7nq5u59c/Xx0XkMItqtYdlQjZOAqqeOX2mth1CqasyVdUv19iaTXTi7g7nvdffwFsAeqqlVmdpSkZ7r7DdVUuU9Muz8AJIBaF0GfoTZWoE21HRgZ4MDL3Fon6VWSrnX3B9x9RdK1ks509++7++elUT2U9A8anal9KY2OsQVQjscfXRNiyssRZratdnmLu29peN8j3f2e6u9vSzpy2kJmdoykv5L0PI1OXPItM9us0ddQjO2SxK5eFCHk8Vvrd+wd/NlKqXXpC72Due9Qm5qDN+4JMqV12eNt2xxri+VlUus2SLq7dnm/mlZ9V/bPSvpvtat/zsx+UqOvHPs/3b3+GPsh2CYstRHOEDjwvhj3ufvmWTea2WclPWfKTe+oX3B3NzOf9hhV8XphNS3vL8zsyjYNRlpKrG8oErUuE32F2hQDbd24fW2DTpuTSS0jp5NIFarzWrcaM1sr6VOSPuzu468W+1+SPuXuj5rZr2g0GvyK1R4niWB73vrrdenKqZ0+R84HZWM2RiXS5O6nz7rNzL5jZke5+z3VdLt75zzWt8xsu6SXSfpb7TtF5WhJbAkBREGtSwOhdn8hRm+XCbelj9oOdYAmQK3breq7sStHS/pC7fIWSbe5+3+tPWf9K8Y+Kul35rWTgwwDC/VhLnmvVcmvDY1slXR+9ff5kj49uYCZHW1mB1V/r5f0v0vaUU11+Z6ZvbQ6Q+gvTbs/ulFyZwXoALWuB4Ta2XJsM7I0t9bpqe/GXl/VujOq62Rm75F0qKT/XL9DFZLHzpL0jXkNIdgC6Nv7JL3SzG6TdHp1WWa22czGZ7z7EUlfNrOvSfobjc6Yd3N1269qtOdup6TbJX2mz8YDXUphxx8n1guGWleInANi27Yvs+OAmXSDM7fWufsDkt4t6cbq5yJ3f8DMjtZoOvOJksZffTb+Wp9fr74C6GuSfl3SG+Y1JImpyCmKuWFPoWOD1T1056FZb+hiqqaWnDbl+m2S3lj9fa2kF864/zZJJ3XZRgBoi1rXva5Ha0vZzoc6qVTXOM42T01qXXX5YkkXTyyzS9LU7+h297dLevsibWG3LKKgcAHDxecfQFuE2sW0eT2M2i6OQ3fiINgGFOJNTIcPAABgNkLtckp9XbnhcI/usGaRvTY7FCguAErFiAGwuNLD37Kvr6+TdOVwxuGuBqFymC6eOnr1iIbRaQDYX461kQ4ZFrXsjuW+AlbJ+grvQ5+OjP4RbBOSY2cGAACgD0xBjoudCkgdZ0UGgAFh+j2AHJUYas84dsd+111z16ZenjuXMyUDiyDYAgB6w8wUAKnpI9ROC7FNl+sq7C4Tbh88/omFdpCubFq38PH+fO0PlkWwnWKZEQ1O0rGcg7bvzuJEAQAAII4uR2u7DrVNA23Tx+hrRBfIEcE2EeyZApAydt4BiCHX4zpDBNrVHjdUwO1j1BboyyDelRxDAAAAkI4UAmsXo7VnHLujs1A7+TyhdD1qzdmR0ZdBBFuUj9EkAADKktMU5L4CbeznHOt6x0QJh6nRN+0fwRYAAABJSWFEt6lY4TLk8/NVR4vhEMI0EWwTwIcDAABgJKev9okdasdSaccsTEdGHwi2iI5gDwD7oi5iqAi1y2vbnkXXTU6j6qnh5FvdYK0GwBx6AACAxUwGo5yCUmqhdizVdi2jhONs0S+CLQAAAKJ48PgnnvzpWqjR2tTDY5v2cawtckawBQAAABpIPdSO9dXORXZIcJwtukawBQAAQNFCjETmEmrbYtQWuSLYRsYJQrCI89ZfH7sJAAAMTo6hNsc2A20QbCdwljIAAIBytB2BzDkgLtv2RdZZl8dHcwIpLIIUBwAAABQqpWDOcbboEsEWAAAARRryaC0wNARbAAASw/Q7IL6SQu0yr4WTSCE3BFsAAAAUh2DWnz6+hxiYh2ALAAAA1JQ0WjuW62tiBguaItgCAACgKG1Ga3MNgF3oYtSbE0ihK2tjNwAAAAAYmmnfTX/pyqmdPucZx+7QNXdt6vQ5gFgItgCAuVY2rdP6HXtjNwMA5kp1tHZakF1tma5DbmgPHv+EDrmdyaCIh2DbEh09AAAAzNIk0K52v9wCLhALu1UAAL3hJCAAupTaaO2yoTb0Y9Qt+jo5uzRywYgtAAAAEFDoMFra6O2ih7c8fNIGHbR9d4ctQgkYsQWAAeG7BgGUKpXR2tChtovH5szPKBHBFgCAhOQ4XZupisBIl6G2z+dYFjtPEVOWwZbTlAMAsDq+KxJoJtToZcqBsy12XiEHWQZbAAAAYCx28Oo71JYcooFlEWwBAAAwSDkfa9o23Ob82oFpCLYoQt9T7mLvGQZyluMxpADSFXubzOhpP9h2YB6CbUscw9QehQoAAPQtxIhl7FAb+/nboA+N0Ai2AHplZoeZ2bVmdlv1e/0qyz7TzHaZ2f+oXfcFM9thZjdVPz/UT8tBJ6R7Kezo46ymYVDrUJqmI+PUkGFpWuvM7PxqmdvM7Pza9VNrnZkdaGaXm9lOM/uymR03ry0EW6BHHM8iSbpQ0nXufoKk66rLs7xb0henXP+L7n5y9XNvF40EgJaodT2IOQ05ldHSNu2gX4IA5tY6MztM0rskvUTSKZLeNRGAp9W6CyStuPvzJH1I0vvnNYRgi0Fjr2IUZ0u6pPr7EkmvnbaQmb1Y0pGSrumnWehTCiOTQMeodQlrG+hSCbVAAprUuldJutbdH3D3FUnXSjpzgce9UtJpZmar3WFt0xYDKMOaR6VDbm+9T+sIM9tWu7zF3bc0vO+R7n5P9fe3NerQ7cPMDpD0QUmvl3T6lMf4IzN7XNKfSnqPu3vzpgMYAmpd+WKfNCol562/XpeunBq7GYOR0s7ZHGqdpA2S7q5d3lVdNzat1j15H3d/zMz2SDpc0n2zGkKwRVQpFQYs5D533zzrRjP7rKTnTLnpHfUL7u5mNq2j9quSrnL3XVN2zv2iu+82s0M0KoDnSfrEQq0HEhSyHnI8dDDUOuyH0VoUqOtat5pgtY5gCyA4d5828iBJMrPvmNlR7n6PmR0ladpxYz8u6WVm9quSDpa0zswecvcL3X139RwPmtknNTpWg84egN5R6/LEcaX5evikDTpo++7YzRicALVut6SX1y4fLekL1WPPqnW7JR0jaZeZrZV0qKT7V2snx9hGxohle4xMZGerpPHZ8M6X9OnJBdz9F939WHc/TtLbJH3C3S80s7VmdoQkmdnTJP2MpO39NDtNy0zFS+XYcuofCket6xDTkMMh6KOlubVO0tWSzjCz9dVJo86QdPWcWld/3NdJ+ty8wzEItgD69j5JrzSz2zQ6pux9kmRmm83so3Pue6BGhfDrkm7SaG/eH3bYVkxgR1I3CPlFotYlqOSTRqXcNhRtbq1z9wc0Ovv7jdXPRdV1q9W6j0k63Mx2SnqrVj+zvKREpiJzsPsw0ZEbJne/X9JpU67fJumNU67/uKSPV3//i6QXd9tCAGiPWtcdRmvLsbJpndbv2Bu7GWihaa1z94slXTyxzMxa5+6PSPr5RdrSaMTWzM6svjh3p5nNTMtm9nNm5mY28+BjAEhVCrVuaFPC2MEVHqPqmCeFWpcjRkTZqYC0zQ22ZrZG0kckvVrSiZLONbMTpyx3iKTfkPTl0I3sUyrHngHo19BqHdKRUrhnG1i+3Gtdm2A1tB2H6EZKNRv7ajJie4qkne5+h7vvlXSZRl+YO+ndkt4v6ZGA7RsEPiBx0IHDBGodgCGg1gEoUpNgO+8LdWVmPybpGHf/q9UeyMzeZGbbzGzbygPlhAqmfS0uVJiPse6ZhlOsTmrd91ceDd/SyEJ/7ti5lzdqYnY6qXWPP/Qv4VuakFymIafSzi4HD3LZZpAP+tf6rMhmdoCk35X0X+Yt6+5b3H2zu29efxgnZAaQj2Vr3TPWH9h945CtXDpoGI5la92ag/9V521jJwqA1TRJl+Mvxx07urpu7BBJJ0n6gpl9U9JLJW3lRAOYhY4cEjWoWpfaVHzqQhiMEKCBQdW6MY6vBcrXJNjeKOkEM9toZusknaPRF+ZKktx9j7sf4e7HVV8wfoOks6pTPKMhOnVAdNQ69Iq6j0iyrHUxR2tTmd6L+KjbaZsbbN39MUlvkXS1pG9IusLdbzGzi8zsrK4bCMzCyARCotbFR4cB6B61DkCp1jZZyN2vknTVxHXvnLHsy9s3C6VKpeOa2jRMpIFa19zKpnVav2Nv7GZkK5VaWEddHI7cal3b0VqmISM11NtucAanhKTY0QGQvpJOqEIdXB6zWAAMBdsKTEOwDYQOxXwUIQCgFgKLiL3jjuNrMUbtTh/BNjF8aJphRwLQXtupUF19DqmDZWNaKACgCwTbKZj3Hh4dVQCLKLVmlPq6gC6EGK1lRwpiYAAmDoItsIDYU6IA5KvLUBurE0VNBDAE7JTMQ5bBNtW9b6E6FqV9eEK/nrbrmRF5IA+l1cKUUReRGnaaYB62EZiUZbBdFMUxHorOU1LdIQO00fVIYSk1pJTXAeSk7XY3xxNHXbpyauwmANEMItjmiE4QAJSh63rOsVwoDQMSSAl98nwQbBOW+wepi/bTgQPCymEKas61MOe2AzEQatP20J2Hxm4CMBPBNjCC10iqnbkcOvHAMmJ2Bvuoe6nWlNXk2GYAQDtkgXgItomjYwQAIznVw77aGqIDxQ4/pCLkDjrOa9HcNXdtCv6Yh9zeT8ToutbmtN0BwXamlDb0uX2oumpv7D1gTI8CZuvr85lDPcyhjUBoaw58vNX92cZiKFLKGKUh2CIoOnRAe32PNOS2kU25zvTZNnb2ITXLvid4L4WR6xmRY9eyWVLe1mA6gm0HuviA5vDh6rKNTLcDUJdiTUyxTUDfFg2phFq0Re3FGME2Iyl/cFNuG4B+9L3X/eGTNiRTe1Jpx6LY4YcuNA2rXYVajq9FW8vW9FRHn4eCYJuZFDtPfEcjgJhi1sVY4Zq6iNQdvHHPzOC62m1AbCn2tdEMwXYVbfZkd9npSOUDl9JoyTyMSmAI2nQUQ31GYgWuvutRTvUPiGkcYus/CG/Z42u7OCNyDNRjSAMKtqUV0tgf4Jy+yiKE0t4/QKm6DpwpBNpU6iKAYXnozkNjN6Fzses72hlMsI2h685HrA4WH3oAq0kheIWujykE2tBizGTh2EcAXWlbo9veP4Vt39Ctjd2AXE9NnpKHT9qgg7bv7uV5+kSBANDWZN1qWitTDbGp1EVmsSBVQ9x5Ql86Hxwa163owTZ1Dx7/hA65Pf2B7XEnrIuAm2oHr6kUisgQN7TIT8h6t7Jpndbv2BvksULKvZ4BAGZbdrCHbUMZCLYd67tzFyrgxv6ApzIqAQCpCFkXU9jhB6C9NqO1pZw4atKi4TZ2nxfhEGwLtcz0u1Q+2Kl13lKacnfpyqk6b/31sZuBhB28cU8yJ/hIddQ2R+zsA5CLHGY6SmH7vSFqdEr9zVwRbHuQQuculdAKIG2hD79Iof4BGIahHfbT57G1Xe4w7WIbMW8GI/3iMhFsG8jlONsSMCoBjJxx7I5ip4lhMaHrYtuZLIwqAPkbyvaFADssg0prMTfGBLb5Uuu8AQiD+rc81h2AaTgTcjqa1mn6pd0bVLCNjQ5KfhiZwBB1sfGl/gEAgC4RbJEERmuBcNghU4YudgZQG1GqLo6vTXVUtG27Fp2GnMoJCYF5CLYNheoMMGqxP9bJsJjZYWZ2rZndVv1eP2O5Y83sGjP7hpndambHVddvNLMvm9lOM7vczHgDdYBR27hYV/mj1gHlo1YvVOvOr5a5zczOr647xMxuqv3cZ2b/tbrtDWb2z7Xb3jivLQTbCPgQPGUIIxJDO0tjAxdKus7dT5B0XXV5mk9I+oC7/4ikUyTdW13/fkkfcvfnSVqRdEHH7UVA1L94Yn/92QBrIbUOwaU6ijzGyVYHaW6tM7PDJL1L0ks0qnPvMrP17v6gu588/pH0T5L+rHbXy2u3f3ReQ3j3IZrUO7hM5+zM2ZIuqf6+RNJrJxcwsxMlrXX3ayXJ3R9y9++bmUl6haQrV7s/0pb6Zz821k8xqHUIKkSoHcrZkFOT2qBLYHNrnaRXSbrW3R9w9xVJ10o6s76AmT1f0g9J+tKyDRlcsG0TVkK+Kem4dKPwwlGKI939nurvb0s6csoyz5f0XTP7MzP7qpl9wMzWSDpc0nfd/bFquV2SOJf/FCF2zHT5eaIGTsd6KQq1rgCpjJDGagfH185GvX5Sk1q3QdLdtcvTato5Go3Qeu26nzOzr5vZlWZ2zLyGZPs9tiV8x+PKpnWdfCl1DigG8ax5xEO8744ws221y1vcfcv4gpl9VtJzptzvHfUL7u5m5lOWWyvpZZJeJOkuSZdLeoOkT7dsNxIy5Bo4TZd1cYg7/ah1w1H6NPdQoTb3fjOmy6TWNXGOpPNql/+XpE+5+6Nm9isajQa/YrUHyDbYxvLg8U8EPX5giB27rjpvITtuTEOe6z533zzrRnc/fdZtZvYdMzvK3e8xs6P01PFkdbsk3eTud1T3+QtJL5V0saRnmdnaaiTjaEm7W7wOzBG65mG6HHb2DbQuUuugS1dO1Xnrr4/dDCQkdM1OoL52Xet2S3p57fLRkr5Qe4z/TaPDMv6+9pz315b/qKTfmfci6K0kIIcOTQgrm9YN5rViVVslnV/9fb6mj0zcqFGn7tnV5VdIurWanvJ5Sa+bc/8ilD4KIVEXpO63AUMcrU0EtQ6txZwKzTRkNNSk1l0t6QwzW1+dNfmM6rqxcyV9qn6HKiSPnSXpG/MaMshgm8Bekf2U3rGj44aa90l6pZndJun06rLMbLOZfVSS3P1xSW+TdJ2Z3SzJJP1hdf/fkvRWM9up0XFoH+u5/dkIVev6+HyVXgNnGerrHghqXUH6DpiXrpwa9Dn7mIY8lNk9i9btAfRRm9S6ByS9W6OdeTdKuqi6buwXNBFsJf26md1iZl+T9OsaHaaxKqYiL6GrqXmlTkvOreOW4o6PklRTS06bcv02SW+sXb5W0gunLHeHRqeKR2FKrYGz9FEbB9ChSha1rh8lzmxJ5YRVQBML1LqLNTrMYtpj/PCU694u6e2LtGUYu1YyklsInGfoHbcSN7gYpr4+Z6XVwGlynH7ddocftRAl6CNwdvEcy4zWMg15utxq99AQbJfE12Csrq+OW8qhFkhBjjMQSqiBs/T52qiPQHhdhttUQi3CoQ73a7DBNvXOXs4du5zbnvr7Aoipzw10jqOa8xBqgTKEDqChj6eNYQjH15a2TSpR1GNsc/8Qd/01GOMPUC7HnPX9gafjBvSv76//KeG429w7Q+zwQ6piTnEf92HbfA1Q1/3gZUdrmYaMXHHyqJb66OSl3rGL0Wkj1GIozjh2x+CnkuW2k28sVqClPgL9WTTg5j6oM1Rd1nN2HoZDsM1Eih273EchJlFYUKqDN+4Juge+71HbsRTr4Cyl1UcAq0stsA59h2gK2MnYv6wnxLedgpLTdzyOpXDcWew25FIoOAso0I3YNWiWcbtKqo/s8AOGZZmdoKUfX5vi9gbTMWKbqfqHrI/Ri1Q+1F2FWjpvKF0po7Z1qYzgll4f22InH9AvRmvDSaW+oxmCbSAxO3ldhNwUP8ipdtqAoUoh3Er97+ibfE4A/WOHyXRtQi0njQqHPmscgw+2IUcxUujkzets1Tt9OXXMuiwQjNYCy0uh7tVNq2shwm7q9bKLGkltBDBkqdd97G/wwXZocvyQstcLQxfqzMihpyPnIse6twhqJAApzmhtiB2bKX77R+nbjVKls5u9EHQwwup6fXY1IsEUKQwJdS8ezjsAQOK42pDahtpF6jK1NiyCrcK/qejkhcF6BMLraiPK57V/OaxzdvKhC7yv9tU21A5xJg/KlH2wTbW45dDhSFkf64+9ZEBY1L3+cN4BAAiPKch5yz7YhtLFhpxO3nJYb0C++PwCQH9ijtamdOLAEEKEWraBcZX1jkwQb/DF9LW+GJFAbkLOTun6/U/d6xbrF0OW6ky9GDiuNhxGastAsO0BnZD5Hjz+iWLWExtdgLrXhT7qJDv9gDyECLUcWzsSKtQuWp+pt+ERbGu6fIPRyZut73VDIQH6+RxQ98LJcV2ykw/oRgojtaVMQ2aktixlvCszkWPHpEsxRmkJtUC/qHvtcYgGgLFQoTa10doYATPkc7KtS0OjYGtmZ5rZDjPbaWYXTrn9rWZ2q5l93cyuM7Pnhm/qbBx7lp9S1wMjFHlLvdaF1leQKfXz3gdCLbqQa60b+jY2hZHaUjBSW6a5wdbM1kj6iKRXSzpR0rlmduLEYl+VtNndXyjpSkm/E7qhfeoj3A61oxfztdNxw2pyqHVddOr6DLdDrXvLYn2hCznUOuwvZKhtO1qb8zTklU3rgodaanU6mrwzT5G0093vcPe9ki6TdHZ9AXf/vLt/v7p4g6SjwzazTEP6IMTu1BJq0QC1rgdDqnvL6rtedlEfhz6yljhqXWaGMlLb9ShqSqO09Eu70STYbpB0d+3yruq6WS6Q9Jk2jUoBoxhhpPD6KB5oaJC1Tur/M5JCXUgV6wU9GGyty1HoUJvasbV96SrUUrPTsjbkg5nZ6yVtlvTvZtz+JklvkqSjNqwJ+dQ649gdWe/RGn8wcp7eUZfKB73PDjsjFMOxSK075Khn9NiyvDx4/BPF1Ly2OEQDKVqk1j3t2d0GpiFuY1MMtV3X7JVN67R+x96gj4fhaPLu3C3pmNrlo6vr9mFmp0t6h6Sz3P3RaQ/k7lvcfbO7b15/WPqdmRgb/NxHMlJqPx02LKiTWveM9QcGbWRXnbtYn5eUakYMnHcAEXRS69Yeyk68kHIerGkrRBjt4ljaSUPedqWqSbq8UdIJZrbRzNZJOkfS1voCZvYiSf9To+J3b5Mnvv/xgxdtaxR09uYbtzWX9gIzdFLrchIz6AyxhpT6eoc4spaZTmrdM9c9EryhQ9VFqM1tCvKyobSPQNsWOxW7M3cqsrs/ZmZvkXS1pDWSLnb3W8zsIknb3H2rpA9IOljSn5iZJN3l7md12O7BqHd8Upqyl3qHrO+iQUcuf9S6kYM37onaASrtsIxJqdROOlbD1WWt6+qwsKFsY1Mfpe27Li8yLbnvMJtKLce+Gh1j6+5XSbpq4rp31v4+PXC7ltJVQY3d0Rub/BD1WWBy+gDTYcOyhl7rUlJawE2phlIj0WWtG0J96kKX6yyFPuyy6oF1HHJjj8imVM+xr6AnjypZKuG2btoHq20nMOcPa6zO2lD2JGM4Uqp3qc5aaSLFetp1naQeQiLcLopQ20zsQBsCOxa7RbBdQEqdvVlS7Ej1gUIBhJVivcsh5A61BgOTQoXbkneW5BT+U625faPGp624d2nXBZAAlZ6h/U8uXTk1dhOQgD46eyl/tuonrYvZ0UilHU2k/P8EhqaPUJvazsncpV7jwYjtUlIcyRiq2B21kvckA1I+9W5WhyPkKEPOnZo+aiX1EJPajtqW+J7qa5Q2h7qdkxD1P3afdQgItkvKpbNXKooD0J+c613OYTQU6iViGofTnKbddiXXdcA0ZOSiyHdqX3v46CzEkcp6L3FPMvLS53swlc8dFsP/DakY8jbzmrs29Rpqc90RmSp2kOajyGDbJzoN/Tl44x7WNxARn7+89Pn/GnJoQXOLvE9KeE/1HWglQm1ooUIt289+FBtsGckoS2rruIQNLsrQ93sxtc8i9sdOQKRsCNvPGIG2K0OehsxIbX44xjaQcSeCPWVhpdg5G8JGGVhNzsfcli5GzaQmYlHzjrvN9T0VO8xSl8MJGWpT7MuWimAbGB2+MCgCQHOhvi9yEdS69FA3kZsSTiyVStupx+EwUpuvoucXxNrjxzSw5aW+7nLdi5wSMzvMzK41s9uq3+unLPNTZnZT7ecRM3ttddvHzezO2m0n9/0aMJL653VIYv0fqImzUeuaO+PYHVm9l8ZTjUsPtUOchkyoXVyTWlct99dm9l0z+8uJ6zea2ZfNbKeZXW5m66rrD6wu76xuP25eW4p/x8YslHT4mqODPCgXSrrO3U+QdF11eR/u/nl3P9ndT5b0Cknfl3RNbZHfHN/u7jf10ObkUeuGidqZNGpdQVILswivi1A7kPo8t9ZVPiDpvCnXv1/Sh9z9eZJWJF1QXX+BpJXq+g9Vy62q+GAbG52O2cbrJpf1k9Pe5MSdLemS6u9LJL12zvKvk/QZd/9+l41COzl9lksRe31TE+ei1i2oPjU5doCsB9nYbVkNU5DDYKS2lUa1zt2vk/Rg/TozM4126l055f71x71S0mnV8jMN4hjbGMefTap3QIZehGJ3xpZRUgfugId/oIO27277MEeY2bba5S3uvqXhfY9093uqv78t6cg5y58j6Xcnrnuvmb1T1Z5Bd3+04XMXLaVaN/Q616Uca2gM1Lr8TNavyctdbYtj181ldVlnhzQNuatQ21etzrDW1R0u6bvu/lh1eZekDdXfGyTdLUnu/piZ7amWv2/Wgw0i2KZmaB2/3DthJYXagO5z982zbjSzz0p6zpSb3lG/4O5uZr7K4xwl6UclXV27+u0aFc51krZI+i1JFzVvOvrAyaXCS6mWDqguUusSkmsA7QL1NQxGap/US63r2mCCbQojGZNKHsVNqQPWxoA6b0G5++mzbjOz75jZUe5+T9WZu3eVh/oFSX/u7j+oPfZ4r+CjZvZHkt4WpNGFSKnWDW0nXldSq6fUxadQ68JJpW5hOLoMtanV7bYC1rpJ90t6lpmtrUZtj5Y0Hn7eLekYSbvMbK2kQ6vlZxrOPAOlvTHO7XjTSfX25/oaJqX8fsncVknnV3+fL+nTqyx7rqRP1a+oiub4uIzXStoevol5S+29W1pt6AvrLHvUuoYItYthZ2F7hNqgFql1+3B3l/R5jc4xMHn/+uO+TtLnquVnGsyIbU4mPxCpFbAhfGBTCwaFeZ+kK8zsAkn/pNFIhcxss6Q3u/sbq8vHabSn7m8m7v/HZvZsSSbpJklv7qfZCIFR3NWlXl+pjQuh1jVAqF1MH7Wz5ONrmXrciaa17kuS/o2kg81sl6QL3P1qjQ6zuMzM3iPpq5I+Vj3uxyRdamY7JT2g0XkIVjW4YJvSNL2mZnV0uixuqXeuukTHrVvufr+k06Zcv03SG2uXv6mnTiBQX+4VXbavFKnXupIPxVhGDjWX2rgYat18KdcolKePUJtDLQ9tgVr3shn3v0PSKVOuf0TSzy/SlsEFWyn9Dl9TQ/zwdI2OG0qSS60basilhmPIcqhNqRlSfQyNUDsMgwy2Uj4dPvSHUIsS5VbrSg65OXd6qI8IKaealIrS6mFfmHo8LIMNtlJ+HT50h04bSpZrrZsWBHPq3OUcZOuojwglxzqEfPUZakup97kbdLCV8u3wIQw6bBiKUmpdjHMOzFNyh4YaibZKqDux5bRDLxWE2mGKFmwPX/NQrKfeTykdPiyGDhuGpuRat0jHYl4nkU7KCDUSi/je3qcXW19iItQupu+px2wv0jL4Eduxkjt82FfunbXz1l8fuwnIGLWOjkgTuddJoASE2uZiHEvLtiQ95X5R1RLOOHYHG/OC8f8FRvgcYDW8PwDkhFCLMUZsp2BEoyx00oD9UecwiVoJpIPR2vlinfGYUJsuRmxnYHQvf/wPgdXx+cAY7wUAY6l/Rc6Dxz9BqMVUjNjOMd7YM7KRBzpnwGKocaBuAmlhtHa62IGbUJs+gm1DdP7SRscMaIepycND3QTSQ6jdX+xAKxFqc0GwXRABNw10yIDwqG/DQP0EMM+Dxz+hQ26Pe8RiCoFWItTmhGC7JDqA/aMzBvSD+lYu6iiQrtRGa2OFWwItlhU12J63/npdunJqzCa0Vu8k0AkMh84XEB8BtxzUVADL6CvcphJmxwi1eWLENiBC7nLocAFpI+DmidoK5CO10dq6rsJtamF2jFCbL4JtRyY7FHQI6WQBuWPnXR6otQBCq4fQZUNuqkF2jECbP4JtT2Z1NErrHNKhAoaBUdy0UHsB9GVWQD3k9gOSD6+zEGrLQLCNbLXOSEodRjpNAKZhdko81GWgDClPQ15EjqGWQFuW6MG2hBNIdYVOC4DcEHS7wzYBAMIh1JYnerAF0Nx566+P3QRgIQTd5RBiAaAbBNpyEWwBAL0ZyvkGFkGIBYarlGnIOSDQlo9gCwCIbl64yzn4ElwBIC5C7TAkEWw5zhYAsJplwmHoMExABYC8EGiHJYlgCwBAaARRABgmAu0wLfcNyx3gpDgAAAAAlnXwxj2E2gFjxBbIBDt/AAAA9keYhZTQiK1Exx0AAABAM4zQoi6pYAsAAAAMBaFsOQRaTJNcsGXUFtgfnwsAADB0BFqsJrlgK9GJBwAAADBCoEUTSQZbiXALAACA8hHYphuHWdYPmko22EqEW0DicwAAQOkIb08hzGJZSQdbiU49AAAAyjfkMMfoLELI4ntsx+H20pVTI7cE6Bc7dgAAGI6DN+7RQ3ceGrsZvSDEIrQsgu0YARdDQqgFAGB4xoGvxIBLmEWXsgq2Y5MdfoIuAAAASlJKwCXMoi9ZBttJBF2UhtFaAAAg5RdwCbKIpVGwNbMzJf03SWskfdTd3zdx+4GSPiHpxZLul/Tv3f2bYZva3KKhIOcgXGIAyvn/0VaJ/89JZnaYpMslHSfpm5J+wd1Xpiz3fkk/XV18t7tfXl2/UdJlkg6X9PeSznP3vYHallWtA5Auah1CqwfGVEJu6SH2jGN37HfdNXdtitCSdC1Q6/5a0ksl/b/u/jO16/9Y0mZJP5D0FUm/4u4/MLOXS/q0pDurRf/M3S9arS1zg62ZrZH0EUmvlLRL0o1mttXdb60tdoGkFXd/npmdI+n9kv79vMdOxRDCRE7a/D9yDsUDeh9eKOk6d3+fmV1YXf6t+gJm9tOSfkzSyZIOlPQFM/uMu39Po/ryIXe/zMz+QKP68/ttGzWEWgegV9Q6dGZaoOwq7A4xvC6yPEF3fq2rfEDSMyT9ysT1fyzp9dXfn5T0Rj1V675UD8HzNBmxPUXSTne/Q5LM7DJJZ0uqF8CzJf129feVkv6HmZm7e9OGACE0DYcpBeABBdqxsyW9vPr7Eklf0P4F8ERJX3T3xyQ9ZmZfl3Smmf2JpFdI+g+1+/+2AnT2RK0DEBa1Dr0qPYCGsmiQbfp4Aw64TWqd3P26ahR28vqrxn+b2VckHb1sQ5oE2w2S7q5d3iXpJbOWcffHzGyPRlNn7lu2YUCX5oXJPoLvAAPt2JHufk/197clHTllma9JepeZfVCjvXs/pVGn63BJ3606gdKoHm0I1C5qHYCQqHVAIkKHWeyjSa2by8yeJuk8Sb9Ru/rHzexrkr4l6W3ufstqj9HryaPM7E2S3lRdfPTk5+7a3ufzL+kI5VHIaWdQV3Tezg+GeZiFdw9+7wf3Xv3Xu//7ES2f9+lmtq12eYu7bxlfMLPPSnrOlPu9o37B3d3M9hsBcPdrzOzfSrpe0j9L+jtJj7dsc28ma90HT76CWhcO7Qwrl3ZS6xI0WetuPuvd1LpwaGdAN2fSThVa6xr6PY1msHypuvwPkp7r7g+Z2Wsk/YWkE1Z7gCbBdrekY2qXj66um7bMLjNbK+lQjU42sI9qBW2RJDPb5u6bGzx/VLQzLNoZ1kQRasTdz+yiLRPPcfqs28zsO2Z2lLvfY2ZHSbp3xmO8V9J7q/t8UtI/alRXnmVma6uRjGn1aFnUOtoZDO0Mi1pHrQuFdoZFO8MqudatxszeJenZqh1/W51rYPz3VWb2e2Z2hLvP3EFxQIPnulHSCWa20czWSTpH0taJZbZKOr/6+3WSPsdxGABmqNeL8zU6490+zGyNmR1e/f1CSS+UdE1VVz6vUZ2Zef8lUesAhEStAzAEc2vdaszsjZJeJelcd3+idv1zzMyqv0/RKLfut4Otbm6wrfYUvkXS1ZK+IekKd7/FzC4ys7OqxT4m6XAz2ynprRqdDQsApnmfpFea2W2STq8uy8w2m9lHq2WeJulLZnarRqMBr68da/Zbkt5a1ZvDNao/rVHrAARGrQMwBE1qnczsS5L+RNJpZrbLzF5V3fQHGh2X+3dmdpOZvbO6/nWStlfH2H5Y0jnzdrBZrB1wZvam+tztVNHOsGhnWLm0c8hy+R/RzrBoZ1i5tHPIcvkf0c6waGdYubQzVdGCLQAAAAAAITQ5xhYAAAAAgGR1HmzN7Ewz22FmO81sv2M0zOxAM7u8uv3LZnZc122apkE732pmt5rZ183sOjN7bortrC33c2bmZhblDHBN2mlmv1Ct01uqM0H2rsH//Vgz+7yZfbX6378mQhsvNrN7zWzq1yjYyIer1/B1M/uxvtsIal3f7awtR61rgFqHUKh1/bazthy1roEcal3VDupdF9y9sx9JayTdLumHJa3T6IvIT5xY5lcl/UH19zmSLu+yTS3a+VOSnlH9/R9TbWe13CGSvijpBkmbU2ynRt9D9VVJ66vLP5RoO7dI+o/V3ydK+maEdv6kpB+TtH3G7a+R9BlJJumlkr7cdxuH/kOt67+d1XLUunDtpNbxE+q9RK0L2M5qOWpduHZGr3XVc1PvOvjpesT2FEk73f0Od98r6TJJZ08sc7akS6q/r9ToTFnWcbsmzW2nu3/e3b9fXbxBo+9961uT9SlJ75b0fkmP9Nm4mibt/GVJH3H3FUly94W/8yqAJu10Sc+s/j5U0rd6bN+oAe5flPTAKoucLekTPnKDRt99eFQ/rUOFWhcWtS4sah1CodaFRa0LK4taJ1HvutJ1sN0g6e7a5V3VdVOX8dEp6PdodFr7PjVpZ90FGu1F6dvcdlZTFY5x97/qs2ETmqzP50t6vpn9rZndYGadf7n0FE3a+duSXm9muyRdJek/9dO0hSz6/kV41LqwqHVhUesQCrUuLGpdWKXUOol6t5S1sRuQGzN7vaTNkv5d7LZMMrMDJP2upDdEbkoTazWatvJyjfaSftHMftTdvxuzUVOcK+nj7v5BM/txSZea2Ule+wJpoETUumCodUDCqHXBUOsQXdcjtrslHVO7fHR13dRlzGytRtMC7u+4XZOatFNmdrqkd0g6y90f7altdfPaeYikkyR9wcy+qdGc/K0RTjTQZH3ukrTV3X/g7ndK+keNCmKfmrTzAklXSJK7/52kp0s6opfWNdfo/YtOUevCotaFRa1DKNS6sKh1YZVS6yTq3XK6PIBXo703d0jaqKcO4n7BxDK/pn1PMnBFl21q0c4XaXRA+gl9t2+Rdk4s/wXFOclAk/V5pqRLqr+P0Gi6xeEJtvMzkt5Q/f0jGh2LYRHW6XGafYKBn9a+Jxj4St/tG/oPta7/dk4sT61r305qHT+h3kvUuoDtnFieWte+nUnUuur5qXeh12kP/7TXaLTX5nZJ76iuu0ijvWPSaE/Jn0jaKekrkn44yoqY387PSvqOpJuqn60ptnNi2SgFsOH6NI2m19wq6WZJ5yTazhMl/W1VHG+SdEaENn5K0j2SfqDRHtELJL1Z0ptr6/Ij1Wu4Odb/fOg/1Lp+2zmxLLWufTupdfyEei9R6wK2c2JZal37dkavdVU7qHcd/Fi18gAAAAAAyFLXx9gCAAAAANApgi0AAAAAIGsEWwAAAABA1gi2AAAAAICsEWwBAAAAAFkj2AIAAAAAskawBQAAAABkjWALAAAAAMja/w8Byp/MmJZZ1AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "We can also plot the pinn loss during the training to see the decrease." + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 12, + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "plt.figure(figsize=(16, 6))\n", + "plotter.plot_loss(pinn, label='Loss')\n", + "\n", + "plt.grid()\n", + "plt.legend()\n", + "plt.show()" + ], + "outputs": [ + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6oAAAFlCAYAAADiXRVWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABFFklEQVR4nO3deXxcdb3/8ff3nNmzTJa2CW26Uyi0lBbaAiI7AoqAuLAp4NUr8rt6vVe9V0G9btcFl+v1ongVETe8iiIoCIIIlEUppaXQspTShbbpmqXZM/v398dMpqFrOktmkryej0cekzkzk/k0p5PkPZ/vYqy1AgAAAACgXDilLgAAAAAAgMEIqgAAAACAskJQBQAAAACUFYIqAAAAAKCsEFQBAAAAAGWFoAoAAAAAKCueUhdwMOPGjbPTpk0rdRkH1Nvbq4qKilKXgSHifI08nLORhfM1snC+RhbO18jC+Rp5OGelsWLFilZr7fj93VbWQXXatGlavnx5qcs4oCVLlujMM88sdRkYIs7XyMM5G1k4XyML52tk4XyNLJyvkYdzVhrGmE0Huo2hvwAAAACAskJQBQAAAACUFYIqAAAAAKCslPUcVQAAAAAY7eLxuJqbmxWJREpdSlEEAgE1NTXJ6/UO+TEEVQAAAAAooebmZlVVVWnatGkyxpS6nIKy1qqtrU3Nzc2aPn36kB/H0F8AAAAAKKFIJKL6+vpRF1IlyRij+vr6w+4WE1QBAAAAoMRGY0gdkMu/rSyDqjHmImPMrZ2dnaUuBQAAAABGvcrKylKX8AZlGVSttfdZa68Lh8OlLgUAAAAAMMzKMqgCAAAAAErr+eef18knn6x58+bp0ksv1e7duyVJN998s4499ljNmzdPV1xxhSTp8ccf1/z58zV//nwtWLBA3d3deT03q/4CAAAAQJn40n0v6eVtXQX9msdOrNYXLppz2I+75ppr9L3vfU9nnHGGPv/5z+tLX/qSvvvd7+qmm27Sxo0b5ff71dHRIUn69re/rVtuuUWnnnqqenp6FAgE8qqZjmqO1u7s1sttyVKXAQAAAAAF19nZqY6ODp1xxhmSpGuvvVZPPPGEJGnevHl673vfqzvuuEMeT7r3eeqpp+oTn/iEbr75ZnV0dGSP54qOao5+8fTr+uMLEf3Tu0pdCQAAAIDRIpfO53C7//779cQTT+i+++7TV7/6Va1evVo33HCDLrzwQj3wwAM69dRT9dBDD2n27Nk5Pwcd1Ry5xihlS10FAAAAABReOBxWbW2tnnzySUnSL3/5S51xxhlKpVLasmWLzjrrLH3jG99QZ2enenp6tH79eh133HH69Kc/rUWLFmnNmjV5PT8d1Rw5DkEVAAAAwOjQ19enpqam7PVPfOIT+vnPf67rr79efX19mjFjhn76058qmUzqfe97nzo7O2Wt1cc+9jHV1NToP/7jP/TYY4/JcRzNmTNHb33rW/Oqh6CaIw9BFQAAAMAokUql9nt86dKl+xx76qmn9jn2ve99r6D1MPQ3R45jlCSoAgAAAEDBEVRz5BojS1AFAAAAgIIjqObIZegvAAAAABQFQTVHrmNkJaVIqwAAAADyZEfxcM1c/m0E1Ry5xkiSkqP4PxQAAACA4gsEAmpraxuVYdVaq7a2NgUCgcN6HKv+5shxMkE1ZeV1S1wMAAAAgBGrqalJzc3NamlpKXUpRREIBN6w9c1QEFRz5MkE1dQofNcDAAAAwPDxer2aPn16qcsoKwz9zZGbCaoJ5qgCAAAAQEERVHPkZOaospgSAAAAABQWQTVH7qA5qgAAAACAwiGo5igbVJmjCgAAAAAFRVDNER1VAAAAACgOgmqOsvuoElQBAAAAoKAIqjka6KimUiUuBAAAAABGGYJqjvZsT0NSBQAAAIBCIqjmyBnoqLKYEgAAAAAUFEE1R3vmqJa4EAAAAAAYZQiqOWLVXwAAAAAoDoJqjgiqAAAAAFAcBNUcuZnvXJI5qgAAAABQUATVHLlO+ltHRxUAAAAACougmqM9iykRVAEAAACgkAiqOXIGhv4SVAEAAACgoAiqORroqLKPKgAAAAAUFkE1Rx6Xob8AAAAAUAwE1Rw5zFEFAAAAgKIgqOaIfVQBAAAAoDgIqjnKBlXmqAIAAABAQRFUc0RHFQAAAACKg6CaI/ZRBQAAAIDiIKjmyHHYngYAAAAAioGgmiMPQ38BAAAAoCjKMqgaYy4yxtza2dlZ6lIOaGB7mgRBFQAAAAAKqiyDqrX2PmvtdeFwuNSlHNDAYkopgioAAAAAFFRZBtWRwMP2NAAAAABQFATVHDnMUQUAAACAoiCo5ojtaQAAAACgOAiqOaKjCgAAAADFQVDNkYd9VAEAAACgKAiqORpY9ZftaQAAAACgsAiqORrYR5XtaQAAAACgsAiqOcpuT5MqcSEAAAAAMMoQVHO0ZzElkioAAAAAFBJBNQ+OkZIspgQAAAAABUVQzYMjhv4CAAAAQKERVPPgOAz9BQAAAIBCI6jmgY4qAAAAABQeQTUPjpFSzFEFAAAAgIIiqObBNVKSfVQBAAAAoKAIqnkwxihBUAUAAACAgiKo5sExUoqgCgAAAAAFRVDNA/uoAgAAAEDhEVTzwBxVAAAAACg8gmoeHIIqAAAAABQcQTUPhqG/AAAAAFBwBNU8uCymBAAAAAAFR1DNg8P2NAAAAABQcATVPLA9DQAAAAAUHkE1D46YowoAAAAAhUZQzQOr/gIAAABA4RFU80BQBQAAAIDCI6jmgaAKAAAAAIVHUM2DY6QUc1QBAAAAoKAIqnlwjNieBgAAAAAKjKCaB8cYtqcBAAAAgAIjqObBMWxPAwAAAACFRlDNg2OkRJKgCgAAAACFRFDNA4spAQAAAEDhEVTzwPY0AAAAAFB4BNU8uEYipwIAAABAYRFU82CMlEilSl0GAAAAAIwqBNU8ODIipwIAAABAYRFU88AcVQAAAAAoPIJqHlwjJQiqAAAAAFBQBNU8sD0NAAAAABQeQTUPhqG/AAAAAFBwBNU8uEZKEVQBAAAAoKAIqnlwmKMKAAAAAAVHUM2DY4ySzFEFAAAAgIIiqObBYegvAAAAABQcQTUPDP0FAAAAgMIjqObBMelLy/BfAAAAACgYgmoeMjmVLWoAAAAAoIAIqnkY6KiyoBIAAAAAFA5BNQ8DQTWVKm0dAAAAADCaEFTz4Jh0UqWjCgAAAACFQ1DNQ7ajSlAFAAAAgIIhqOZhYDEl9lIFAAAAgMIhqOYhu5gSQRUAAAAACoagmgdW/QUAAACAwiOo5oFVfwEAAACg8AiqecjOUaWjCgAAAAAFQ1DNA3NUAQAAAKDwCKp5YHsaAAAAACg8gmoejEknVTqqAAAAAFA4BNU87OmolrYOAAAAABhNCKp5GPjmMfQXAAAAAArHM1xPZIx5h6QLJVVL+om19i/D9dzFYspsMSVrrW57cqN+/vTrmlIX0pcunqNZDVWlLgsAAAAADsuQOqrGmNuNMbuMMS/udfwCY8yrxph1xpgbDvY1rLV/sNZ+SNL1ki7PveTyUW6r/t793FZ99YFXNKkmqLU7u/XO//27/r6+tdRlAQAAAMBhGerQ359JumDwAWOMK+kWSW+VdKykK40xxxpjjjPG/GmvjwmDHvq5zONGvHJa9Xd3b0xffeAVnTi1Vr/+0Mn640ffrMbqgN7/02d1z8pm2TKoEQAAAACGYkhDf621Txhjpu11eLGkddbaDZJkjPmNpEustV+X9Pa9v4ZJL5F7k6Q/W2ufy6vqMlFOiyn99G8btbsvpl+94yQ5jtGkmqB+++FT9KFfLNfH73xBP/3b6zp91njNmVitM44er5Bv2EZ9F0wimVIkkVIknlR04DKeUiSRzB6LDrotEk8pmthzuW5DVA/vXr3nsYmU4snUAZ/PHOC46zjyOEYe18jjGLmOI69r5Domc3zP7QP39Xuc9IfXVcDryO9xM8cGXfem7xPw7rnN73HkOAeqBAAAABid8kkrkyRtGXS9WdJJB7n/P0s6V1LYGHOktfaH+7uTMeY6SddJUkNDg5YsWZJHicUVjUQkGS1fvkId693S1ZG0uv3JPs0f72rnq89p56t7bvun2VaPhnz627Zu/WBJp1JWCnmky4/26fQmT3aLnWJKWav+hNQbt9mPnrjUF7fqyVzviyv7eX9CiietYqn0ZTwlxVL5vSHgGMnrWPm2b5bXMfI5ktc18hjtP5Ee4Lms0nUkrVXSZj5PScnMscHX0/fLveYBXkcKuJLfY9KXrlHAIwU8Jv35fm7zu0ZBj1ThNQp5jSo8Ushr5HNHVujt6ekp658BeCPO18jC+RpZOF8jC+dr5OGclZ9ha6tZa2+WdPMQ7nerpFslaeHChfbMM88scmW5e/GuRyRFdPyCBVo0ra5kdfzqmU3qib+oG96xSCfNqN/n9rMzl5F4Us9v6dDNj7ymn77UpkTVEfrixXPkFqBjF00k9drOHq3e2qlVzZ3a1Narlu6oWnui6uiP62Ajj4NeVzUhr8JBr2prvJoS8Crg3dNZDAzqQg6+zHYeva4CmY5k9tLrKpC9zZHHdbRkyRIN9/8na60SKatYIpXu+Cb2dIGj8T3HBrq+A8cib+gMJ9UbS6gvmrmMJdUbTagnmtTO/j3Xo4kDd4cH+DyOwkFv9qM64Ml+XlfhV32lT+MqfRpX6Vd9Zfp6lX943tDYn1KcM+SO8zWycL5GFs7XyML5Gnk4Z+Unn6C6VdLkQdebMsfGjHJYTCmVsvrJUxs1rymsxdMPHpYDXlcnz6jXoml1+saDa3TrExvU3hvTdy4/Xn7P0DvCvdGElm5o05od3Vq7s1trtndrfUuPEpnvQ3XAoyMnVGrm+EqdNKNOdSGfqoNe1YR8qgl6FQ55s5fhoPewnnukMcbI6xp5XUcV/uI+VyKZUl88qb5oUj3RhHqiCXX2x9XZH1fXoMuuSDx7vLUnpvUtvdnr++PzOBpX4VN9pV/jKn1qDAd0RDioI8IBTawJqjEc0MRwUEHf6D2PAAAAGF75BNVnJc0yxkxXOqBeIemqglQ1QmTnqJYwqD740g5taOnV969aMOSul+sYfeZtx2h8pV9ffeAV7e6L6RvvmqfJdaEDPqY7Eteja3bpgdXbteTVlmz3bmI4oKMbq3TOMRN0zBHVmtcU1pS6UMk6cGOZx3VU7TqqDnhzenw8mdLu3phae2Jq7YmqrTeq1u6YWnujassc29Ud1eqtnWrtie3z+JqQV43VATXVBjW1vkLT6kOZywpNrAnI47JtMwAAAIZmSEHVGPNrSWdKGmeMaZb0BWvtT4wxH5X0kCRX0u3W2peKVmkZKvViSqmU1S2PrdOM8RV669wjDvvxHzp9huoqfLrx7tU6/VuPaVJNUAGvK2vtnimaVoomUtre2a+UlRqq/bpy8RSdd2yD5jaFcw5FKD9e19GE6oAmVAcOed9IPKmdXRFt64hoR1e/tnVEtL2zXzs6I9rS3q+n1rUqEt8zFNnjGE2uC2lqfUhHjq/UUY1VOrqhSrMaKkfkwl4AAAAorqGu+nvlAY4/IOmBglY0ggz0DJMl2vrlt8u36KVtXfrOZcfnPM/0XSc26U1H1uuu5c1a39Kj+MDqP2bPv8/nOppcF9Jps8bphCm1rEILBbyuptZXaGp9xX5vt9ZqV3dUr7f2alNbn15v69Xrbb3a2Nqnp9e3ZTvyxkiTa0M6qqFKsxurdOzEdFd+Uk2QrjwAAMAYRisjD6Uc+tu8u09f//MaLZ5ep0sXTMrrax0RDuqfz5lVoMqA9NzchuqAGqoD+yzwlUxZbW7v06uZOc6v7uzW2h3deuzVXdn53vUVPh3XFFY4EVNiwk4dP7lG46uKPMkXAAAAZYOgmodSLabUHYnrgz9brpS1uumdx9F5wojiOkbTx1Vo+rgKXTC3MXs8mkjq1R3deqG5U6u2dGhVc6ce3xnXH9cvlyRNH1ehxdPqtHh6+qOplq4rAADAaEVQzUM2qA7j0N9kyupff/O81rX06Of/sFgzxlcO23MDxeT3uJrXVKN5TTXSyVMlSQ/99THVHXm8ntu0W8++3q4HX9qhO5ent2+eGA5o0fQ6nTpznE4/arwaw4eeWwsAAICRgaCaByfTzbHDGFTvWLpJj6zZpS9fMkdvnjVu2J4XKAW/x2jRtDotmlanD58xU6mU1dpd3Vq2sV3LNrbr7+vb9Mfnt0mSjm6o0ulHjdMZR03Qwmm1CnjZLgcAAGCkKsugaoy5SNJFRx55ZKlLOajsYkqpg96tYGKJlL736DqdPKNOV2c6TsBY4jhGsxurNbuxWtecMk3WWr26s1tPrG3R42tb9PO/b9KPn9yogNfRabPG6/w5jTr3mAmqCflKXToAAAAOQ1kGVWvtfZLuW7hw4YdKXcvBDPfQ37+8vEOtPVF9+z3zmJsHKL1o00Bwve70meqLJbR0Q5uWvNqih1/eqYdf3inXMTppep3On9Oo8+Y06IhwsNRlAwAA4BDKMqiOFMO96u+DL+7QuEqfTps1flieDxhpQj6Pzp7doLNnN+hLF8/RquZOPfTSDj300g594d6X9IV7X9LCqbW6ZMEkXXjcEaqroNMKAABQjgiqeRjOVX+TKavH17borXMbc94zFRhLjDE6fnKNjp9co09dMFvrdvXowRe3694Xtuk//vCivnTvSzr9qPG6ZP5EnXdso4I+5rQCAACUC4JqHgbiYmoYhv6+uqNb3ZGETplZf+g7A9jHkRMq9dGzZ+kjZx2pV7Z364/Pb9Ufn9+mR9fsUsjn6u3zjtDli6bohCk1DK0HAAAoMYJqHrJDf4chqK7Y1C5JWjStrujPBYxmxhgdO7Fax06s1qcvmK1nNrbrnpXN+tOq7frt8mYd1VCpKxZN0aULJqmWocEAAAAlQVDNw56hv8V/rpe3d6k25NWkGhaCAQrFcYxOmVmvU2bW6/MXzdF9L2zTb57doi//6WXd9Oc1umBuo65cPEUnz6ijywoAADCMCKp5MMO46u+aHd06qqGKP5aBIqn0e3Tl4im6cvEUvbK9S3c+u0V3P9ese1/YptmNVXr/m6bpkvmTmMsKAAAwDJxSFzCSOZnQaIscVFMpq7U7ujW7saqozwMg7ZgjqvXFi+do2WfP1Tffnd4O6oa7V+uUmx7RTX9eo20d/aUuEQAAYFSjo5qHgZRf7FV/t3b0qzeW1NGN1UV9HgBvFPC6umzhZL3nxCYt29iun/39dd36xHr9+MkNOn9Og/7xtBk6YUptqcsEAAAYdQiqeTDDtD3Nmh3dkqSj6agCJWGM0Ukz6nXSjHo17+7TL5du0m+WbdEDq3fo5Bl1+n9nHqnTZ41jaD4AAECBMPQ3D8O16u/anemgelRDZVGfB8ChNdWGdONbj9Hfbzhbn7vwGL3e2qdrb1+mC29+Sve+sE2J4VhdDQAAYJQry6BqjLnIGHNrZ2dnqUs5qOFa9XdDS68aqv2qCniL+0QAhqzC79E/njZDT3zqLH3z3fMUSST1sV+v1Nn/9bh+9cwmxRIEVgAAgFyVZVC11t5nrb0uHA6XupSDGq6O6pb2Pk2tqyjqcwDIjc/j6LKFk/XXj5+hH77vRNWGvPrsPS/qrG8v0a+XbVacDisAAMBhK8ugOlIMzEZLFXmO6qb2Xk2pDxX1OQDkx3GMLpjbqD985FT97B8WaVyVXzfevVpn/9cS/fbZLQRWAACAw0BQzYMzDPuoRuJJ7eyKakodQRUYCYwxOvPoCfrDP71Jt79/oWqCPn3q96t07nce110rmpnDCgAAMAQE1TwM7KNazI7qlvY+SdJUOqrAiGKM0dmzG3TvR0/VbdcsVKXfo3/73Qt6y38/oXtWNhd9tXAAAICRjKCaJ9cxRe2obmpLB9XJdFSBEckYo3OPbdCf/vnN+tHVJyrgdfXxO1/Q+d99Qn9ata3oUwcAAABGIoJqnlxjVMy/MzcPdFQJqsCIZozR+XMadf8/v1k/eO8JMpI++n8r9babn9SDL+6QLfKibAAAACMJQTVPxhR36O/m9j5V+j2qq/AV7TkADB/HMXrbcUfowX89Xf9zxXzFEildf8cKXfT9p/Tomp0EVgAAABFU8+Y6pqhzzTa392lyXUjGmEPfGcCI4TpGl8yfpL98/HR9+z3Hq7M/rg/8bLku/cHf9cTaFgIrAAAY0wiqeXJNseeo9jLsFxjFPK6jd5/YpEc/eaa+/s7jtKsromtuX6bLfvS0nl7fVuryAAAASoKgmifHMUUb+ptKWW3Z3c8eqsAY4HUdXbl4ih779zP15UvmaFNbn6788VJd9eOlWv56e6nLAwAAGFYE1Ty5TvEWU9rZHVEskWIPVWAM8XtcXXPKND3xqbP0H28/Vmt3duvdP3xa19y+TM9v6Sh1eQAAAMOCoJonx6hoQ38HtqZhD1Vg7Al4XX3wzdP1xKfO0g1vna3VzR16xy1/01U/XqrHmcMKAABGubIMqsaYi4wxt3Z2dpa6lENyTPGG/g5sTUNHFRi7Qj6Prj9jpp789Nn6zNtma31Lj669fZne+j9P6p6VzYonU6UuEQAAoODKMqhaa++z1l4XDodLXcohFXPV381tfXIdo4k1waJ8fQAjR6Xfo+tOn6knP3W2vvXueUqmrD5+5ws645uP6bYnN6izP17qEgEAAAqmLIPqSOLkuervru6I3nHL3/TTv23c57bN7X2aWBOQ1+U0AUjzeRy9Z+FkPfSvp+sn1y5UU21IX7n/FZ38tUd0492r9PK2rlKXCAAAkDdPqQsY6RxHymeq2O9XbNXzWzr0/JYOvXNBk8Ihb/a2ja29mlZfUYAqAYw2jmN0zjENOueYBr24tVO/ePp13f3cVv162RYtmlarq0+ZpgvmNMrn4Y0uAAAw8vAXTJ5ck9/Q38fW7Mp+/vArO7OfW2u1vqVHM8dX5lUfgNFv7qSwvvnu4/XMZ87RZ942Wzu7ovrYr1fqTTc9oq//+RVtaOkpdYkAAACHhaCaJ8fJb+jv+pYeXbawSbUhr57Z0JY9vqMror5YUjMnEFQBDE1NyKfrTp+pJf92pn76/kVaMKVWtz25UWf/1+O67IdP6/crmtUfS5a6TAAAgENi6G+e3DxW/e3si6utN6aZ4yu1e1pcy15vz962oaVXkjRzHEN/ARwexzE6a/YEnTV7gnZ1RfT757bqzmc365O/e0FfvPclXTR/oq5YNFnHTQrLGFPqcgEAAPZBUM2T6xilcuyobmhND8ebMb5SjjF6+OWd2tUV0YTqgNZnhurRUQWQjwnVAf2/M2fq+jNmaNnGdt357Bbd/Vyz/u+ZzTrmiGpdvrBJ71gwSTUhX6lLBQAAyGLob56MMcp1G8PtnRFJUlNtUIum10lStqu6flePKnyuJlT5C1IngLHNGKOTZtTrO5fP1zOfOVf/+Y658jhGX7zvZS3+2iP6l9+s1LKN7bL5rA4HAABQIHRU8+Q6yrmj2tIdlSSNr/IrHPQq6HX17MZ2vX3eRK1v6dXMCZUMywNQcOGgV1efPFVXnzxVL23r1J3PbtE9K7fqj89v01ENlXrvSVN16QmTVB3wHvqLAQAAFAEd1Tzls+pva09UjpFqQz55XUcnTK3Rstd3y1qrF7d1anZjVYGrBYA3mjMxrC9fMlfPfOYcffNd8xT0uvrCvS/ppK8+ok/ftUqrmztLXSIAABiD6KjmycljjmpLd1T1lX65Trprumhanf7nkde0qrlTHX1xLZhSW8hSAeCAQj6PLls0WZctmqzVzZ36v2Wb9IeV23Tn8i2a1xTWe0+aoouOn6iQb+z+2ognU2rvjam+wiePy/u8AAAU09j9i6NAXJNfUB1fuWcO6uJpdbJW+t6j6yRJJ04lqAIYfsc1hfX1pnm68W3H6A8rt+qOpZv06d+v1lfuf0XvOqFJ7z1pimY1jN4RH6lUeh/rFZt2a82Obq1v6dH6XT3a3hWRtVLI5+rEqbWaPq5C85pqNKHKr8XT6xTwuqUuHQCAUYOgmicnz6G/4wYtlrRgSq1cx+ivr+zUzPEVmsWKvwBKqDrg1TWnTNPVJ0/V8k27dcfSTfq/ZzbrZ39/XYun1+l9J0/VBXMa5fOM7O5iJJ7UquZOLd/UrhWv79aKzbvV0ReXlA6lM8dX6qQZ9ZpSF9K4Sp/Wt/Tq4Zd36pkN7frF05skSZV+j3weR28+cpzqK31638lTNakmSHgFACBHBNU8OY6UynHV3919cU0btE9q0Ofqn86cqR8sWa9/PfcoFlICUBaMMVo0rU6LptXp82+P6ncrmvWrZzbpY79eqXGVPr1n4WRdtXhKqcscsq5IXCte361lr7fr2Y3tWtXcqVhm+fYZ4yt0/rGNOnFarRZmuqb7+1n8xYvnKJZI6al1LdrZFdVTr7WqL5bQo2t2qSea0E//9rqCXlfvWDBR8yfX6JL5kwitAAAchrIMqsaYiyRddOSRR5a6lENyHaN4jvvTdPbHFQ6+cVXNT553tD58xkxV+svy1AAY4+or/br+jJm67rQZeuK1Ft2xdLN+9Ph6/fDx9Zpb76qtqlnnHtuwz8+2UonEk3p1R7dWb+3Ui1s79UJzp9bs6JK1kscxmjsprPefOk0Lp9Zq4bQ61VUMfT9Zn8fR2bMbJElXDgrqW9r79PDLO7V8U7vufX6bfr1siz79+9WaUhfS3EnVmjMxrKMaqjRrQqUm14Wy6xQAAIA9yjINWWvvk3TfwoULP1TqWg4l16G/qZRVV2TfoCqJkAqg7DmO0ZlHT9CZR0/Q1o5+3blss+74+3p98ncvyOc6Om3WOJ0/p1GnzKzX5LpQ0euJJ1Pa1Nan9S09WrerR+tberRme7fW7uxWIvMzuibk1XGTwvqXc2Zp8bQ6zZ9SU5TFoSbXhfSBN0/XB948XdZaLdvYridfa9XG1l6t3tqpB1bvyN7X53E0c3yljm6o1IIptTpxaq1mN1axWBMAYMwjEeXJdYxsDospdUcTslZl03UAgFxNqgnqE+cdrfnebaqZOV/3r9quP6/erkfW7JIkNdUGddL0eh1zRJVmN1ZrxvgKja/yyzuEMJZMWXX2x9XeG1NHX0ytPVFt74xoR2cke7mts187OiPZQCpJjdUBzWqo1HVHz9Bxk8KaOymsptrgsE+pMMbopBn1OmlGffZYVySudbt6tG5nj17b1a3XdvXo6Q1t+sPz2yRJQa+r4yeHdeLUdHBdOK2OPW0BAGMOQTVPjjFK5hBUu/rTC3VUE1QBjBKOMTphSq1OmFKrz77tmHQAW9+qpRva9fjaXfr9c83Z+xoj1Vf4VR1IL0IU8LoyRorGU4olU4omkuqOJNTZH9f+fsT6XEeN4YAawwEtnFqriTVBHTmhUjPHV2rG+ApVlXGwqw54s9+nAdZabeuMaMWm3Xpu026t2LRbP3x8g5IpK49jdMrMel18/ESdP7eR0AoAGBMIqnlKD/09/Md1ZoIqHVUAo5HjGB3dWKWjG6v0/lOnS0qvdP7qjm5tauvTjq6IdnVF1BNNKJpIKZpIyVqr+gpHPo8jn+uoKuBVbYVPtSGv6ip8qgn5VF/h0xHhgOoqfKNqwTljjCbVBDWpJqiLj58oSeqLJfT8lg49vrZF96/arn+/a5U++4cX9ZZjGnTtm6Zp0bTaUfU9AABgMIJqnlwnPd/0cHURVAGMMeMq/Rp3pF+nlv86eWUh5PPoTTPH6U0zx+mGC2Zr5ZYO3fv8Nt2zcqvuX71dcyZW6wOnTtcl8ycypxUAMOrwmy1PrpPb0F86qgCAoTKZYdVfvHiOlt54jr526XGKJVL65O9e0NtuflKPr20pdYkAABQUQTVPjjFK5bKYUiQhSaoK0NQGAAxd0OfqqpOm6C8fP10/fN8JiiZSuvb2Zbr29mVau7O71OUBAFAQBNU8OcbkNPS3J5oOqmxFAwDIhTFGF8w9Qn/5+On63IXH6LnNu/W2/3lS//3wWsUSue3vDQBAuSCo5inXob+9maBaQVAFAOTB73H1j6fN0OP/fpYuOn6i/ueR13Tx95/S+paeUpcGAEDOCKp5SndUD/9xPbGEfB5nSPsIAgBwKHUVPv335fN12zULtas7qku+/zc99NKOUpcFAEBOSEl5cp30hvSHqzeaYNgvAKDgzj22Qff985s1Y3yFPvzLFfrWQ2ty+j0FAEApEVTzlOtiSj2RhCr8bhEqAgCMdZNqgvrth0/R5Qsn65bH1uuj//ecIvFkqcsCAGDICKp5cpwcg2o0qQofHVUAQHEEvK6+8e55+tyFx+jPL+7QtbcvU1ckXuqyAAAYEoJqnlxjGPoLAChb/3jaDH338vlasWm33nfbM+rsI6wCAMofQTVPrpNjUI0lWPEXADAs3rFgkn509Ylas71b7/3JUnX0xUpdEgAAB0VQzZNjjHIY+aseOqoAgGF0zjEN+tHVJ2rtzh5dc/syRRIssAQAKF8E1Tw5Rjnvo8piSgCA4XTW7An6wVUn6KVtXfr+yqhiiRz2VwMAYBgQVPOU69Df/lhSIRZTAgAMs3OPbdDXLp2rF9uS+vTvVynF1jUAgDJUlkHVGHORMebWzs7OUpdySLmu+huJpxTw0lEFAAy/yxdN0TtneXXPyq266cE1pS4HAIB9lGVQtdbeZ629LhwOl7qUQ8pl1d9EMqVYMqUgQRUAUCIXzfDq6pOn6tYnNugPK7eWuhwAAN6gLIPqSJLuqB7eYyKZOUFBH99+AEBpGGP0+YuO1eJpdbrx7tV6ZXtXqUsCACCLpJQnx6QvD2eOT38sKUl0VAEAJeV1HX3/qgWqCnj0vtueUUt3tNQlAQAgiaCaN9ekk+rhrPwbiaeDKnNUAQClNqE6oF9+8CR19Mf1kf97Tn2xRKlLAgCAoJovJ9NSPZx5qv2ZoBr0EVQBAKV3dGOVPnfhMVq2sV23PLau1OUAAEBQzZebCaqHs/IvQ38BAOXm/W+apvmTa3TLY+v18Ms7S10OAGCMI6jmaWDo7+EsqJTtqBJUAQBlwhijX3xwsWpCXn33r2vVFYmXuiQAwBhGUM1TJqfmNPQ3wNBfAEAZqQ549Y13zdPand268talau+NlbokAMAYRVDNU3bo72EE1QhDfwEAZer8OY265aoT9NK2Lp37nce1qytS6pIAAGMQQTVPHufwV/1l6C8AoJydN6dRt79/oXqiCb33tme0ua2v1CUBAMYYgmqeWPUXADAanT27QT/7h0Xa1R3Vv9y58rBGDgEAkC+Cap48uQTVGPuoAgDK35tmjtOnLjhaKzd36OO/fb7U5QAAxhCCap4cc/hBNcLQXwDACHHV4in68Bkz9Mfnt+nxtS2lLgcAMEYQVPPkcdNBNXGYQ39dx8ibeSwAAOXKGKNPvOUoTasP6UO/WK4Vm9pLXRIAYAwgqObJddLfwsMb+ptS0OvKGIIqAKD8+T2ufvnBkxQOevW5P7ykbvZYBQAUGUE1T24OQ3/740nmpwIARpTJdSHd9M7jtHZnt975g79np7EAAFAMBNU8uTksphSJJxX08a0HAIws5xzToFuvPlGv7erRh3+5QrFEqtQlAQBGKdJSnnIJqv2xJAspAQBGpHOOadCXL5mjx9e26Nt/ebXU5QAARilPqQsY6Qa2p0mkhv6ucn+coAoAGLmuOWWa1u7s1q1PbNDRDVV614lNpS4JADDK0FHN00BHNWWZowoAGDs+//Y5WjytTv95/8va0RkpdTkAgFGGoJqngaCaSB7uHFWCKgBg5PJ5HH3pkjmKJ1L6wM+eVTzJfFUAQOGUZVA1xlxkjLm1s7Oz1KUcUnaO6uF0VJmjCgAYBY45olrfuXy+Xt7epe/+da3sYfwuBADgYMoyqFpr77PWXhcOh0tdyiF5cllMiTmqAIBR4vw5jXrngkm65bH1+u3yLaUuBwAwSpRlUB1JnOxiSoc39DfA0F8AwCjxrfccrxOn1uo///SKntu8u9TlAABGAYJqngY6qim2pwEAjFGuY3TzlQtUV+HTdb9Yoe2d/aUuCQAwwhFU8+SYHDqqiZQCXr71AIDRY1JNULddu1D9sYQ+/MsVisSTpS4JADCCkZby5HEPr6OaSKaUTFkFPHRUAQCjy1ENVfrvy+drVXOnbrx7NYsrAQByRlDNk2c/c1QP9os5mkgv3++nowoAGIXOm9Ooj597lO5ZuVUPrN5R6nIAACMUaSlPA0N/kymrSDyps/9ria788dID3j8bVOmoAgBGqY+cNVNzJlbrhrtX6bWd3aUuBwAwAhFU8+Rx0t/CZMrqxa2d2tDSq6Ub2g84FDiaSM/Z8Xv41gMARieP6+hHV5+ogNfVh36xXF2ReKlLAgCMMKSlPGVyqpIpq13d0ezx1p7ofu8fjTP0FwAw+jXVhvT9Kxdoy+5+XX3bM0okU6UuCQAwgpCW8pTtqFqrXV2R7PHmjv0vzc/QXwDAWHHSjHp957Lj9UJzp/79rlUsrgQAGDKCap7cQYsp7RzUUd3S3rff+w8s18/QXwDAWHDx8RP1kbNm6p6VW/W/j68vdTkAgBGCtJSngaCaTKa0qyuqkC/dKW3tie33/nRUAQBjiTFGn3jL0XrLsQ365oOv6q8v7yx1SQCAEYCgmqdsULXpeakzxldIkrr6979wRHYxJeaoAgDGCNcxuvmKBZoxrkLX37FCL27tLHVJAIAyR1rKUzaoplLqiSYUDnpVFfAccIXD7GJKDP0FAIwhQZ+r311/iqoCHl3146Xa2Npb6pIAAGWMtJQnTzaoSr3RhEI+j6oDXnUesKPK0F8AwNhUX+nXr687WSkrXfajp7WpjbAKANg/gmqeBndU+2JJVfo9qg561dWf2O/9B4b+Bhj6CwAYg2Y3VutHV5+ozr64zvr2Er26o7vUJQEAyhBpKU+u2bPqb7qj6ioc9BxkjiodVQDA2HbqkeP0pUvmKGWl87/7hFZu3l3qkgAAZYagmifHMTJGSqWsemMJVfjTQ38PPEeV7WkAALhy8RR95KyZkqRLf/B3LX+9vcQVAQDKCWmpAFxjFE2mFImnVOHzKBz0HrqjytBfAMAY9+/nz9Y/vnm6JOndP3xam9v2vwc5AGDsIS0VgOsYdUfSc1Ir/K6qAl51RQ40RzUdVH0u33oAAD791tl6/5umSZIu/cHf9OCL22WtLW1RAICSIy0VgMcx2Q5qyOdRpd9Vbyyx31+00URSHsfIQ1AFAEBe19FnLzxGN1+5QG29MV1/x3P61F2rFE+mSl0aAKCESEsF4Dgm20Gt8LsK+T2yVorE9/0lG42nmJ8KAMAgXtfRxcdP1L0fPVWS9LsVzfrxkxtKXBUAoJRITAUwuKNa4fMo5Euv6Nsb23f4bzSRkt/Lir8AAOxtXlONHv746aoKePTNB1/VNx5cU+qSAAAlQlAtANcx2VV+K/wehXweSVJ/LLnPfaOJJB1VAAAOYFZDlZ761Nl694lN+t8l67XwK3/Vrq5IqcsCAAwzElMBuI5RV/+eob8Vh+qoElQBADigcMirr116nE6bNU6tPVEt/tojWrOjq9RlAQCGEYmpADyOo87+mKT0YkrBgaAa3U9HNZ6S38PQXwAADsbncfTLD56kT11wtCTpku//TT9+YgOLLAHAGEFQLQC/11E8mV7ht9LvUYU/PfS3b78d1SR7qAIAMET/dOaReuLfz9IpM+v11Qde0WnfeIz9VgFgDCjLxGSMucgYc2tnZ2epSxmSwKAOacjvZhdT6tvvHFWG/gIAcDim1Id0y1Un6LRZ47SjK6LTv/WYHnxxR6nLAgAUUVkmJmvtfdba68LhcKlLGZKBob6SFPK62cWU9t9RZegvAACHq8Lv0c//YbFueOtsGSNdf8cK3Xj3au3qZqElABiNyjKojjSBzFDegNeRx3X2LKa0vzmqrPoLAEBOHMfo+jNmatUXztPVJ0/VXSu26MxvLdHtT21UT3TfN4cBACMXiakABob+VmQ6qaGDzVGNp5ijCgBAHqoCXv3nO+bqN9edrJDPoy//6WXN/cJD+tUzm9Sd2S4OADCykZgKIODNBNVMQA16D9ZRTb1hTisAAMjNiVPr9NdPnK55TempQp+950Ud98W/6IePry9xZQCAfBFUC2AgqA4souQ6Rj6Po0hi36AaibPqLwAAhVIT8unej75ZP//A4uyxm/68RtNuuF9/W9dawsoAAPkgMRXAwBzVgY6qJAU8jqLxffd6YzElAAAK74yjxmvj19+mm955nLyukSS997ZndOWtS3XH0k3a2tFf4goBAIfDc+i74FD27qhK6ZWA+/e7PQ2LKQEAUAzGGF2xeIquWDxFL23r1F9e2qkfPr5eT29ok5R+Y/nfzjtal8yfpPFV/hJXCwA4GBJTAQx0VOsqfIOOufsM/bXWso8qAADDYM7EsD7+lqP0/OfP0zsXTJIkReIpfeX+V7Toq3/Vlbcu1Ta6rABQtuioFsDA4ki1oT1BNejdt6MaT1pZK/m9DP0FAGA4BH2uvnP5fH3l0rn63qPr9L9L0gstPb2hTW+66VH5PI5mjq/U5QubdNbsCZpaX1HiigEAEkG1IOLJ9FzUcNCbPeb3uook3jhHNZrpsNJRBQBgeIV8Hn36gtn65FuO0sbWXq1v6dGLW7v0/cfW6ZXtXfrifS/ri/e9LEn66FlH6pL5EzWlPsS6EgBQIgTVAujsT+/ZVj0oqAa9jiJ7dVSjmeBKUAUAoDQ8rqNZDVWa1VClC+YeoU+ed5Tufm6rfrBknda39EqSvv/YOn3/sXXZxwS9ro5qrFJTbVDzJoX17hObVF/JHFcAKCaCagF0RRKSpOrAoFV/va7ae2NvuN+eoMq7swAAlANjjN51YpPedWKTIvGkOvvj2tLepw0tvbr3hW16al2r+uNJvbClQy9s6dD9q7br639eoyPCAU2tD6mjL66Z4yu1eHqdTps1TpPrQnKMkbVWjjFyHJNTXVva+yRJk+tChfznAsCIQVAtgIuOP0L3rNyqk2fUZ4/tb45qNJ4Z+ss+qgAAlJ2A11XA66qhOqCF0+p02aLJ6o8l9dzm3frLSzv0m2e3ZN90jsSTWrqhXZK0Zke37l+9/aBfu6Har9mN1Zp9RJUaqwOqDfkU9LmqDnhlrVXKSm29UW1p79Mvnt6kXd1RSdLnLjxGdz67RWET0a6KLVq1tUMnTa/Xto5+NYYDisSTmlZfoWTKqqM/rqn1ISVTVrUhn6KJlII+V6mUVSJlFfA6coxRMmUV8LryuiZ93Vr5PY68TvrvE2PSAR4ASomgWgBnz27Q6zdd+IZjQa+r/jhDfwEAGMmCPlenHjlOpx45Tl+6ZK5SKZsNch19Mb20rUvPbdqtja29qgx4FA56tb6lRw+/vFPxpM1+nZ1dUe3satHja1sO6/m/cv8r2c+X/36VJOmOpZsL84/bizGSY4xS1mYXivS4RiGfK4/jyOdx5HVN9vOQz1UwE+79XkcB78B1RwGPmwn+jvyZ+wQ8TvbNgJAv/VHp96jC71HQ6x60+xxLpNTaE1UyZekyA2MEQbVI/F5Xkfjeiykx9BcAgJFscJiqCfmyIXZ/EsmUEpluZiyRUiKZkowUjacUTSQViafUFYkrnrTq6Itpa0e/Vjd36q+v7Am58yfXSJK88W4dM71J61t6tGByrTr6Y2qsDmhnV1ThoFeVAY/6Y0kFfa5298bUGA6osz+eDZmxZEohn6u+WFIhn6veaFLxZEpe11HKWiVTVvFkSvGklccxisSTMia9Y0F/LKl4Kn1bIplSPJlSNJFSfyypjr64IomkovGUIvFk+iORUjJl9/s9ORiPY1QT8mbDrMcx6o4ktHWvbYQWTq2VxzVat6tXVy2erJDfo75oQoun18t1jMJBr5rqgjJSpnNMgwAYiQiqRRL0uors3VGNs+ovAABjhcd1lH1v+jDXXrI2HXAHQtaSJUt05plzC1tgEcWTA8E1fTkQzAeO9ceT6o0m1BNNqDeaUH88qWgipY6+uPpiCXX1x9XeG1PrXnvSS9LyTbuzn9/86LpBt6zb575+j6O6Cp+CXlf1lT4FfR5V+l01VgfTHWGfqwlVfoV8HlUHPaoOeFWV6YxXB72EXKCECKpFEvA6+wbVgY4qc1QBAMBBGGPkdUfuPFGv68jrOqoK5P+1Uimrzv642nqjiietGqoD2t0XU3tvTB7HqCuS0M7OiPxeR5vb+hTPdHN7IukA3B2Jq6Mvrv54Ui3dUb22M6bO/hZF4kkdqvHr9ziq9HtUGfCoJuRTbcirqoBXlX5X1UGvaoI+hYNe1YS8Cgf3fPTFray1zPUF8kBQLZKg11UiM4xm4N04hv4CAAAcHscxqq3wqbbClz1WV+HTzPH5fd3s8Ov+uLoicXVmwmx3JKHO/rg6++PZrm93JKGO/rjaemLa1Nan7ki66xtLpg749Z1HH3hDeA2H0qG2NuRVbSb01lb4VFfhS1+v8Kkus8gWAIJq0QS86R8ykXhyUFBNZm6jowoAAFBKfo8rv8dVOOjN6fHWWkXiKXX2x9XRH1NHXzwbcFesXqPxE6dkbttzfHNbrzr60x3eAwl4HdWFfKoJZUJshU91Ie8brg+E3boKn5Ipq9aeqDa396k/ltRxTWHNmRiWtel2MV1djFQE1SIJZN4N648nVRVI/wCMxumoAgAAjAbGGAUz81wbw28c4zyhZ73OPPPoAz42kUwH3N198eww5o6+mNp7976eXmRrd1/soOF2b43VAfXHk/I4Rm86cpx8rqMKv6sZ4yrk87hqqParMRxQ0OvqiHCQLi7KEkG1SIKZjmp00Mq/bE8DAAAAj+uovtKv+sqhr7K1J9zGtLsvvdjU7t6YOvrjaumOqq0nqt5YUs9v6dDRjVXa3N6nrv64nl7fqtae2EG/dsiXXmyqNuRTfYVPdRV+hYNe1VV4s0OSayvSt9VW+FQT9MrDQlMoMoJqkQwM7x28l+rA0F86qgAAADgcuYRbKT1E2VopkkiqtTum9r6Y+qIJ7eqOqrM/rq7+uFp7omrrjak7klBLT1Sv7uhOz9GN7bvq8oB0kE0PQx6YZ5sdqjwo3A58XhXwHHSvXGBvBNUiCQ6aozqAVX8BAAAwnIwxMkYK+TyaUu/RlPrQkB8bTaT3ym3riWWHJGcve2Nq74trd29M2zoienFrl9r7Yool9r/AlOuYPQtJZcJrXeXgQLtn3u3AZcjnMsd2DCOoFsnAYkr9g96JGhgG7GOoBAAAAMqc3+OqodpVQ/XQ9hmy1qovlnxjoM3MvW3vjabn4Pamu7rrW3q0fFN6GHPyAPsE+TzphaXqMp3Z+kqfxlX6s5fjstf9qq/wZf/+xuhAUC2SbFDda+ivz3UY9gAAAIBRxxijCr9HFX6PJtcNrXObSll1RxJq74vtE2Z39+4Ju609MW1u71NrT1R9BxiSXOX3aFxVOrQeKNCOq/SpvtKv6oCHbm2ZI6gWycAc1cheiymxkBIAAACQ5jhG4ZBX4ZBX08dVDOkxfbGE2npiau2JqrUnpraeaPbz1szn61t69MzGqHYfYLVkn8fRuApfNrzGe6Ja2r9GDdV+BbyueqMJNe/u1/gqv95zYpMmDLGrjMIhqBbJ/uaoRuJJ+QiqAAAAQM5CPo9CdUPr2iaSKbX3xrIhtq03qtbumFozl229UbX0RLWtLalndmxQPLnvMORvPfSqGqsD6osldFRDleZOCssxRg3Vfh07sVrjq/waV+lXbcgnl5GTBUNQLZKB/ajeGFRTjJ0HAAAAhonHdTShOnDIjuiSJUt0xhlnaHtnRK090eyiUK/t6tH6XT1q7YmqvS+ubR39+s2zmxVLpLT31FrHSHUV6Q7t+Cp/+qMyfVmX6d6Or/RrXJVP9RV+Qu0hEFSLJODZd45qJJFkxV8AAACgDBljNLEmqIk1weyxhdPq9nvfWCKltt6oNrf1vWHIcWtPVC3dUbX0xLShpVctg0LvYAOhdiDQZsNtpX+fkBsOesfkfFqCapHs6agOmqMaT2YDLAAAAICRyedxdEQ4qCPCwYPez1qrrkgiM/w4qtbu9FDjgcuW7vTH+l09aumOKpbcN9R6XZPpxO4VZDNDjgeH2gr/6Il3o+dfUmYGFk3q32foLx1VAAAAYCwwxigc9CocPPRiUdZadfUn1NIT0a7u9OJQA0G2pTvdrd3eGdGqrZ1q64nuM/RYSq+TM7gje8HcRr1jwaQi/euKi6BaJMYYBbyOonstpsQcVQAAAAB7M2bPCshHTqg66H2TKavdffsG2ZZBndoNrT3a1tk/TNUXHkG1iAJed585qtVBbwkrAgAAADDSuY7J7BHr1zFHlLqa4mAcahEFve5+Vv3lWw4AAAAAB0NqKqJ0R3XPhOgIiykBAAAAwCERVIso4HXVH3tjR9XPHFUAAAAAOCiCahEFvY6iiT1BNRpPMvQXAAAAAA6B1FRE+3RUE6z6CwAAAACHQlAtoqDXVSTTUU2mrOJJyxxVAAAAADgEgmoRDe6oDqz+y9BfAAAAADi4skxNxpiLjDG3dnZ2lrqUvAS8riKZVX/3BFU6qgAAAABwMGUZVK2191lrrwuHw6UuJS8Br5MNqJFEKnsMAAAAAHBgpKYiCnrdPUGVjioAAAAADAlBtYiCPlf98aSstdmg6mcxJQAAAAA4KIJqEQW8rlJWiidtdq4qQ38BAAAA4OBITUU0MMy3P55UNMHQXwAAAAAYCoJqEQ10TyPxpKLZjipBFQAAAAAOhqBaRMFMKI3Ek+yjCgAAAABDRGoqosFDfyMDQ39ZTAkAAAAADoqgWkR7OqqpQYspEVQBAAAA4GAIqkXkzwzz7Y8x9BcAAAAAhorUVETZjmoiSUcVAAAAAIaIoFpEQV8mqA7qqPo9fMsBAAAA4GBITUU0sHBSJJFeTMnvcWSMKXFVAAAAAFDeCKpFNNBR7Y+lFI2nGPYLAAAAAENAUC2ibEc1s48qCykBAAAAwKGRnIoo4Mus+psNqnRUAQAAAOBQCKpF5HMdGTPQUU1lO6wAAAAAgAMjqBaRMUZBr6tIPKneWEIVfoIqAAAAABwKQbXIgl5X/fGkuiMJVQa8pS4HAAAAAMoeQbXIAl5XkXhKPdGEqvyeUpcDAAAAAGWPoFpkQZ+rvlhCPRGG/gIAAADAUBBUi6w64FFXf0K90YQq/Qz9BQAAAIBDIagWWXXQq47+mHpiCVUGGPoLAAAAAIdCUC2ycNCrHZ0RWSvmqAIAAADAEBBUi6w64FVrT0yS6KgCAAAAwBAQVIssHNwzL7WSjioAAAAAHBJBtciqg3vCKR1VAAAAADg0gmqRDe6oNlYHSlgJAAAAAIwMBNUiGxxUm2qDJawEAAAAAEYGgmqRzW6szn5eFWAfVQAAAAA4FIJqkU2tD5W6BAAAAAAYUVjdp8iMMfrmu+bJ6zGlLgUAAAAARgSC6jC4bNHkUpcAAAAAACMGQ38BAAAAAGWFoAoAAAAAKCsEVQAAAABAWSGoAgAAAADKCkEVAAAAAFBWCKoAAAAAgLJCUAUAAAAAlBWCKgAAAACgrBBUAQAAAABlhaAKAAAAACgrBFUAAAAAQFkhqAIAAAAAygpBFQAAAABQVoy1ttQ1HJAxpkXSplLXcRDjJLWWuggMGedr5OGcjSycr5GF8zWycL5GFs7XyMM5K42p1trx+7uhrINquTPGLLfWLix1HRgaztfIwzkbWThfIwvna2ThfI0snK+Rh3NWfhj6CwAAAAAoKwRVAAAAAEBZIajm59ZSF4DDwvkaeThnIwvna2ThfI0snK+RhfM18nDOygxzVAEAAAAAZYWOKgAAAACgrBBUc2SMucAY86oxZp0x5oZS1wPJGDPZGPOYMeZlY8xLxph/yRyvM8Y8bIx5LXNZmzlujDE3Z87hKmPMCaX9F4xNxhjXGLPSGPOnzPXpxphnMuflTmOML3Pcn7m+LnP7tJIWPgYZY2qMMXcZY9YYY14xxpzC66t8GWM+nvlZ+KIx5tfGmACvr/JijLndGLPLGPPioGOH/Zoyxlybuf9rxphrS/FvGQsOcL6+lfmZuMoYc48xpmbQbTdmzterxpjzBx3nb8hhsL/zNei2TxpjrDFmXOY6r68yRFDNgTHGlXSLpLdKOlbSlcaYY0tbFSQlJH3SWnuspJMlfSRzXm6Q9Ii1dpakRzLXpfT5m5X5uE7S/w5/yZD0L5JeGXT9G5L+21p7pKTdkj6YOf5BSbszx/87cz8Mr/+R9KC1drak45U+b7y+ypAxZpKkj0laaK2dK8mVdIV4fZWbn0m6YK9jh/WaMsbUSfqCpJMkLZb0hYFwi4L7mfY9Xw9LmmutnSdpraQbJSnz98cVkuZkHvODzBuz/A05fH6mfc+XjDGTJZ0nafOgw7y+yhBBNTeLJa2z1m6w1sYk/UbSJSWuacyz1m631j6X+bxb6T+iJyl9bn6eudvPJb0j8/klkn5h05ZKqjHGHDG8VY9txpgmSRdKui1z3Ug6W9Jdmbvsfb4GzuNdks7J3B/DwBgTlnS6pJ9IkrU2Zq3tEK+vcuaRFDTGeCSFJG0Xr6+yYq19QlL7XocP9zV1vqSHrbXt1trdSgenff44R/72d76stX+x1iYyV5dKasp8fomk31hro9bajZLWKf33I39DDpMDvL6k9Jtxn5I0eKEeXl9liKCam0mStgy63pw5hjKRGba2QNIzkhqstdszN+2Q1JD5nPNYet9V+pdFKnO9XlLHoF/6g89J9nxlbu/M3B/DY7qkFkk/Nemh2rcZYyrE66ssWWu3Svq20h2D7Uq/XlaI19dIcLivKV5r5eMDkv6c+ZzzVYaMMZdI2mqtfWGvmzhfZYigilHHGFMp6feS/tVa2zX4Npte5pqlrsuAMebtknZZa1eUuhYMiUfSCZL+11q7QFKv9gxJlMTrq5xkhqZdovQbDBMlVYguwIjDa2rkMMZ8VukpSL8qdS3YP2NMSNJnJH2+1LVgaAiqudkqafKg602ZYygxY4xX6ZD6K2vt3ZnDOweGHGYud2WOcx5L61RJFxtjXld66NPZSs+BrMkMVZTeeE6y5ytze1hS23AWPMY1S2q21j6TuX6X0sGV11d5OlfSRmtti7U2LulupV9zvL7K3+G+pnitlZgx5v2S3i7pvXbPvo+cr/IzU+k3717I/O3RJOk5Y0yjOF9liaCam2clzcqsnuhTerL8vSWuaczLzKf6iaRXrLXfGXTTvZIGVmm7VtIfBx2/JrPS28mSOgcNt0KRWWtvtNY2WWunKf0aetRa+15Jj0l6d+Zue5+vgfP47sz96TQME2vtDklbjDFHZw6dI+ll8foqV5slnWyMCWV+Ng6cL15f5e9wX1MPSTrPGFOb6aSflzmGYWCMuUDpKSwXW2v7Bt10r6QrTHpF7elKL9KzTPwNWTLW2tXW2gnW2mmZvz2aJZ2Q+f3G66sMeQ59F+zNWpswxnxU6f+orqTbrbUvlbgspLsFV0tabYx5PnPsM5JukvRbY8wHJW2SdFnmtgckvU3pBQ76JP3DsFaLA/m0pN8YY74iaaUyi/dkLn9pjFmn9OIIV5SovrHsnyX9KvPH1QalXzOOeH2VHWvtM8aYuyQ9p/RwxJWSbpV0v3h9lQ1jzK8lnSlpnDGmWenVRQ/rd5a1tt0Y859KByBJ+rK1dn8LyCBPBzhfN0ryS3o4s/7YUmvt9dbal4wxv1X6DaKEpI9Ya5OZr8PfkMNgf+fLWvuTA9yd11cZMrxhCgAAAAAoJwz9BQAAAACUFYIqAAAAAKCsEFQBAAAAAGWFoAoAAAAAKCsEVQAAAABAWSGoAgAAAADKCkEVAAAAAFBWCKoAAAAAgLLy/wGcz6fnG/mPQgAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "You can now trying improving the training by changing network, optimizer and its parameters, changin the sampling points,or adding extra features!" + ], + "metadata": {} + } + ], + "metadata": { + "kernelspec": { + "name": "python3", + "display_name": "Python 3.9.0 64-bit" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.0" + }, + "interpreter": { + "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorials/tutorial3/tutorial.py b/tutorials/tutorial3/tutorial.py new file mode 100644 index 0000000..75ea325 --- /dev/null +++ b/tutorials/tutorial3/tutorial.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python +# coding: utf-8 + +# # Tutorial 3: resolution of wave equation with custom Network + +# ### The problem solution + +# In this tutorial we present how to solve the wave equation using the `SpatialProblem` and `TimeDependentProblem` class, and the `Network` class for building custom **torch** networks. +# +# The problem is written in the following form: +# +# \begin{equation} +# \begin{cases} +# \Delta u(x,y,t) = \frac{\partial^2}{\partial t^2} u(x,y,t) \quad \text{in } D, \\\\ +# u(x, y, t=0) = \sin(\pi x)\sin(\pi y)\cos(\sqrt{2}\pi), \\\\ +# u(x, y, t) = 0 \quad \text{on } \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4, +# \end{cases} +# \end{equation} +# +# where $D$ is a square domain $[0,1]^2$, and $\Gamma_i$, with $i=1,...,4$, are the boundaries of the square, and the velocity in the standard wave equation is fixed to one. + +# First of all, some useful imports. + +# In[2]: + + +import torch + +from pina.problem import SpatialProblem, TimeDependentProblem +from pina.operators import nabla, grad +from pina.model import Network +from pina import Condition, Span, PINN, Plotter + + +# Now, the wave problem is written in PINA code as a class, inheriting from `SpatialProblem` and `TimeDependentProblem` since we deal with spatial, and time dependent variables. The equations are written as `conditions` that should be satisfied in the corresponding domains. `truth_solution` is the exact solution which will be compared with the predicted one. + +# In[3]: + + +class Wave(TimeDependentProblem, SpatialProblem): + output_variables = ['u'] + spatial_domain = Span({'x': [0, 1], 'y': [0, 1]}) + temporal_domain = Span({'t': [0, 1]}) + + def wave_equation(input_, output_): + u_t = grad(output_, input_, components=['u'], d=['t']) + u_tt = grad(u_t, input_, components=['dudt'], d=['t']) + nabla_u = nabla(output_, input_, components=['u'], d=['x', 'y']) + return nabla_u - u_tt + + def nil_dirichlet(input_, output_): + value = 0.0 + return output_.extract(['u']) - value + + def initial_condition(input_, output_): + u_expected = (torch.sin(torch.pi*input_.extract(['x'])) * + torch.sin(torch.pi*input_.extract(['y']))) + return output_.extract(['u']) - u_expected + + conditions = { + 'gamma1': Condition(Span({'x': [0, 1], 'y': 1, 't': [0, 1]}), nil_dirichlet), + 'gamma2': Condition(Span({'x': [0, 1], 'y': 0, 't': [0, 1]}), nil_dirichlet), + 'gamma3': Condition(Span({'x': 1, 'y': [0, 1], 't': [0, 1]}), nil_dirichlet), + 'gamma4': Condition(Span({'x': 0, 'y': [0, 1], 't': [0, 1]}), nil_dirichlet), + 't0': Condition(Span({'x': [0, 1], 'y': [0, 1], 't': 0}), initial_condition), + 'D': Condition(Span({'x': [0, 1], 'y': [0, 1], 't': [0, 1]}), wave_equation), + } + + def wave_sol(self, pts): + return (torch.sin(torch.pi*pts.extract(['x'])) * + torch.sin(torch.pi*pts.extract(['y'])) * + torch.cos(torch.sqrt(torch.tensor(2.))*torch.pi*pts.extract(['t']))) + + truth_solution = wave_sol + +problem = Wave() + + +# After the problem, a **torch** model is needed to solve the PINN. With the `Network` class the users can convert any **torch** model in a **PINA** model which uses label tensors with a single line of code. We will write a simple residual network using linear layers. Here we implement a simple residual network composed by linear torch layers. +# +# This neural network takes as input the coordinates (in this case $x$, $y$ and $t$) and provides the unkwown field of the Wave problem. The residual of the equations are evaluated at several sampling points (which the user can manipulate using the method `span_pts`) and the loss minimized by the neural network is the sum of the residuals. + +# In[4]: + + +class TorchNet(torch.nn.Module): + + def __init__(self): + super().__init__() + + self.residual = torch.nn.Sequential(torch.nn.Linear(3, 24), + torch.nn.Tanh(), + torch.nn.Linear(24, 3)) + + self.mlp = torch.nn.Sequential(torch.nn.Linear(3, 64), + torch.nn.Tanh(), + torch.nn.Linear(64, 1)) + def forward(self, x): + residual_x = self.residual(x) + return self.mlp(x + residual_x) + +# model definition +model = Network(model = TorchNet(), + input_variables=problem.input_variables, + output_variables=problem.output_variables, + extra_features=None) + + +# In this tutorial, the neural network is trained for 2000 epochs with a learning rate of 0.001. These parameters can be modified as desired. +# We highlight that the generation of the sampling points and the train is here encapsulated within the function `generate_samples_and_train`, but only for saving some lines of code in the next cells; that function is not mandatory in the **PINA** framework. The training takes approximately one minute. + +# In[5]: + + +def generate_samples_and_train(model, problem): + # generate pinn object + pinn = PINN(problem, model, lr=0.001) + + pinn.span_pts(1000, 'random', locations=['D','t0', 'gamma1', 'gamma2', 'gamma3', 'gamma4']) + pinn.train(1500, 150) + return pinn + + +pinn = generate_samples_and_train(model, problem) + + +# After the training is completed one can now plot some results using the `Plotter` class of **PINA**. + +# In[11]: + + +plotter = Plotter() + +# plotting at fixed time t = 0.6 +plotter.plot(pinn, fixed_variables={'t': 0.6}) + + +# We can also plot the pinn loss during the training to see the decrease. + +# In[12]: + + +import matplotlib.pyplot as plt + +plt.figure(figsize=(16, 6)) +plotter.plot_loss(pinn, label='Loss') + +plt.grid() +plt.legend() +plt.show() + + +# You can now trying improving the training by changing network, optimizer and its parameters, changin the sampling points,or adding extra features!