From 02fe53ce82779f2dc4f8940f03220a429bb86f3c Mon Sep 17 00:00:00 2001 From: elizabethg60 Date: Thu, 14 Nov 2024 13:11:32 -0500 Subject: [PATCH] gpu geometry updates --- ...bug_gpu_eclipse.jl => debug_gpu_eclipse.jl | 74 ++++++++--------- gpu_test.png | Bin 0 -> 14049 bytes scripts/NEID_GRASS_all_lines.jl | 7 +- scripts/projected_RV.jl | 5 +- scripts/set_geometry.jl | 37 +++++---- src/GRASS.jl | 2 +- src/eclipse_comp.jl | 5 +- src/gpu/gpu_physics.jl | 6 -- src/gpu/gpu_precomps_eclipse.jl | 77 +++++++++++++----- 9 files changed, 117 insertions(+), 96 deletions(-) rename scripts/debug_gpu_eclipse.jl => debug_gpu_eclipse.jl (73%) create mode 100644 gpu_test.png diff --git a/scripts/debug_gpu_eclipse.jl b/debug_gpu_eclipse.jl similarity index 73% rename from scripts/debug_gpu_eclipse.jl rename to debug_gpu_eclipse.jl index f450d101..98cf0a38 100644 --- a/scripts/debug_gpu_eclipse.jl +++ b/debug_gpu_eclipse.jl @@ -1,6 +1,10 @@ using Revise using GRASS using SPICE +using CSV +using DataFrames +import PyPlot +plt = PyPlot #NEID location obs_lat = 31.9583 @@ -9,26 +13,30 @@ alt = 2.097938 #convert from utc to et as needed by SPICE time_stamps = utc2et.(["2023-10-14T14:29:23.500000"]) -# time_stamps = utc2et.(["2023-10-14T17:29:23.500000"]) N = 50 Nt = length(time_stamps) Nsubgrid = 4 disk = GRASS.DiskParamsEclipse(N=N, Nt=Nt, Nsubgrid=Nsubgrid) +lines = [5434.5232] # array of line centers +depths = [0.6] # array of line depths +templates = ["FeI_5434"] # template data to use +variability = trues(length(lines)) # whether or not the bisectors should "dance" +blueshifts = zeros(length(lines)) # set convective blueshift value +resolution = 7e5 # spectral resolution + +# make the disk and spec composite type instances +spec = GRASS.SpecParams(lines=lines, depths=depths, variability=variability, + blueshifts=blueshifts, templates=templates, resolution=resolution) + """ CPU GEOMETRY CODE CALLS BELOW """ wsp = GRASS.SynthWorkspaceEclipse(disk, 1, Nt, verbose=true) mem = GRASS.GeoWorkspaceEclipse(disk, 1, Nt) - -mu_grid_eclipse = Vector{Matrix{Matrix{Float64}}}(undef,size(time_stamps)...) -dA_total_proj_eclipse = Vector{Matrix{Matrix{Float64}}}(undef,size(time_stamps)...) -idx1_eclipse = Vector{Matrix{Matrix{Int}}}(undef,size(time_stamps)...) -idx3_eclipse = Vector{Matrix{Matrix{Int}}}(undef,size(time_stamps)...) -z_rot_eclipse = Vector{Matrix{Matrix{Float64}}}(undef,size(time_stamps)...) -zenith_eclipse = Vector{Matrix{Matrix{Float64}}}(undef,size(time_stamps)...) +LD_type = "KSSD" mu_grid_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ)) dA_total_proj_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.Nθ)) @@ -41,55 +49,37 @@ zenith_matrix = Matrix{Matrix{Float64}}(undef, length(disk.ϕc), maximum(disk.N t = 1 GRASS.eclipse_compute_quantities!(disk, time_stamps[t], obs_long, obs_lat, alt, wsp.ϕc, wsp.θc, - wsp.μs, wsp.dA, wsp.xyz, wsp.ax_codes, #wsp.z_rot, + wsp.μs, wsp.dA, wsp.xyz, wsp.ax_codes, mem.dA_total_proj_mean, mem.mean_weight_v_no_cb, mem.mean_weight_v_earth_orb, mem.pole_vector_grid, mem.SP_sun_pos, mem.SP_sun_vel, mem.SP_bary, mem.SP_bary_pos, mem.SP_bary_vel, mem.OP_bary, mem.mu_grid, mem.projected_velocities_no_cb, mem.distance, mem.v_scalar_grid, mem.v_earth_orb_proj, t, mu_grid_matrix, - dA_total_proj_matrix, idx1_matrix, idx3_matrix, z_rot_matrix, zenith_matrix) -# GRASS.eclipse_compute_quantities_updated!(disk, time_stamps[t], obs_id, obs_long, obs_lat, alt, wsp.ϕc, wsp.θc, -# wsp.μs, wsp.dA, wsp.xyz, wsp.ax_codes, -# mem.dA_total_proj_mean, mem.mean_weight_v_no_cb, -# mem.mean_weight_v_earth_orb, mem.pole_vector_grid, -# mem.SP_sun_pos, mem.SP_sun_vel, mem.SP_bary, mem.SP_bary_pos, -# mem.SP_bary_vel, mem.OP_bary, mem.OP_ltt, -# mem.mu_grid, mem.projected_velocities_no_cb, -# mem.distance, mem.v_scalar_grid, mem.v_earth_orb_proj, t, mu_grid_matrix, -# dA_total_proj_matrix, idx1_matrix, idx3_matrix, z_rot_matrix, zenith_matrix) -mu_grid_eclipse[t] = deepcopy(mu_grid_matrix) -dA_total_proj_eclipse[t] = deepcopy(dA_total_proj_matrix) -idx1_eclipse[t] = deepcopy(idx1_matrix) -idx3_eclipse[t] = deepcopy(idx3_matrix) -z_rot_eclipse[t] = deepcopy(z_rot_matrix) -zenith_eclipse[t] = deepcopy(zenith_matrix) - + dA_total_proj_matrix, idx1_matrix, idx3_matrix, z_rot_matrix, zenith_matrix) + +extinction_coeff = DataFrame(CSV.File("data/NEID_three_extinction.csv")) +neid_ext_coeff = extinction_coeff[extinction_coeff[!, "Wavelength"] .== lines[1], "Ext1"] +mean_intensity = GRASS.eclipse_compute_intensity(disk, lines, neid_ext_coeff, LD_type, idx1_matrix, idx3_matrix, + mu_grid_matrix, mem.mean_weight_v_no_cb[:, :, t], mem.mean_weight_v_earth_orb[:, :, t], + z_rot_matrix, dA_total_proj_matrix, wsp.ld, wsp.z_rot, zenith_matrix, "true", wsp.ext) + """ GPU GEOMETRY CODE CALLS BELOW """ - -disk = GRASS.DiskParamsEclipse(N=N, Nt=Nt, Nsubgrid=Nsubgrid) - -lines = [15652.79] # array of line centers -depths = [0.25] # array of line depths -templates = ["FeI_5434"] # template data to use -variability = trues(length(lines)) # whether or not the bisectors should "dance" -blueshifts = zeros(length(lines)) # set convective blueshift value -resolution = 7e5 # spectral resolution - -# make the disk and spec composite type instances -spec = GRASS.SpecParams(lines=lines, depths=depths, variability=variability, - blueshifts=blueshifts, templates=templates, resolution=resolution) + gpu_allocs = GRASS.GPUAllocsEclipse(spec, disk, 1) -GRASS.calc_eclipse_quantities_gpu!(time_stamps[t], obs_long, obs_lat, alt, [5500.0], disk, gpu_allocs) +GRASS.calc_eclipse_quantities_gpu!(time_stamps[t], obs_long, obs_lat, alt, lines, LD_type, 1.0, neid_ext_coeff[1], disk, gpu_allocs) # make a plot -plt.imshow(wsp.z_rot .- Array(gpu_allocs.z_rot)) +plt.imshow(wsp.ext .- Array(gpu_allocs.ext)) # plt.imshow(wsp.z_rot .* GRASS.c_ms, vmin=-2000, vmax=2000, cmap="seismic") # plt.imshow(Array(gpu_allocs.z_rot) .* GRASS.c_ms .- wsp.z_rot .* GRASS.c_ms) plt.colorbar() -plt.show() +plt.savefig("gpu_test.png") +plt.clf() +#z_rot & dA & ext + # # get size of sub-tiled grid # Nϕ = disk.N diff --git a/gpu_test.png b/gpu_test.png new file mode 100644 index 0000000000000000000000000000000000000000..17d3d26e2131e28da2e8b2c93eb32482a0abfb11 GIT binary patch literal 14049 zcmeHuc{r78+x}yyR3y7Z$dD-|RE9E_A_l`kRuBa%`>|)x5AP9}(#q;V2 zLXJZavLUMN@D-8fz7hEA4CcaBjE4P9jEl)_Gep@0<6vWtv9Y|)?re73(bC@Tl)!0$ zll<%!7>t9XxS*izuNMf|-!>OyUnjhS583H(QP&YcXid;RWQj5fmIxxMqImwSrt7Qe zelLvXQUr0n^Zt|jOgtyTZd}|+{+)^cXqYTjc8Cby-NN+LZu%a~wh-?UrNq7YqN`dP z{fRrao%OuPFz*_|qTly^S8t76^tlucTEUeR2={9pyrYkRb?*4?=|?Lqe3Y z@HJl!gd9OmQ#$F_P6<2QbsN8EsimWHOi-}VOW^3yOLy+vNls0j%2hTqGpk$*ejXV)IJ!iY zHeNmh2Z1W0qeDAF8r!ps<*y5_8{Ku9ftfj_tE+2w`sVkPjJD?H=(e^t+(J@PlKAbB ziwP?$Ze@suIcHC{P^U0XHvpu@R+o;A+PFd>eDMta+Lc%c|PQ?_D=T`3uCKh-E1+iO>%tXE*yE-_mg8^A2KJeu)^ z5iKg};b+{DZsm1)syz*_&S)5UuEm828VO`{roetA$HY`sG!2h8$NKrqXqA^E!oZehy<@VP*)z!u<)IcXZ78K_b?)VoGp=z1Pkbb z30?MeTFz``y~zpaVa#ZS&{1nt!d3QEQ?bK#&VP*CQ)H9p+l}H0V~u%))Kso-cNjSU zuLI-H_t*0oOV~EcpztCD=<^nR1OO#}(g#M41~+DB?dp=(Tcjo3MD>0lZxpGF!}eTA zs+NDDqd>Lp{o(x2rzc16ng68-?tloc=AUnAA~UAVr28smRS7OYfwVZ%MI$*UK923% z*@XGqB3xkqxd>C8$jNOs;_)~VX41*R#yUFFU2IIaHT`RarC-jH-=xo6BR!}OlPsm& z@3-XPcb(G$^O&@x;{2iF{LJe}kHurc-K{OzdZ54z9kSjtRdLbLXD?j10LT5(<;!iQ zcNQ=|Mn2m0+VmBwX~bWC+B+RmH>z%GYTAV#+rM?la$;=QzNl+zzM+-bxKiRg-2!Ew zXFn)?{Ay-YuV{iwL`$w^D~7na;+tlqLb|VUdS6;L?F$zk>|vEmOikriYzY^!*1U2h zivF1Pk>khHTAd#=2!?Q#t$ll0)~2y=g0!xOopu#wXRAMY^oTIWs|bMM=jW%crKNS% zlSjPv=XMWVG92tA_V)F$I7o|%YG-9-NtI`GPHJ7cbg8mveB1)c>*M6o@@jb@2Hjws zN?bz1F%FKmUe8~=aB?%Sv`o?$&?%^W%BkR1v@*VjbQ)_2@=fNw9UUE?rk&RtT}1E& zml8*6-ZCf3Ji`yI*GtSbzmfuyybs@T11qX#&!7byvT z^(u8_WF)F{bYvtkD@(1tz5UbY&(FE0yK@u(KD)cSrDg;L4xz`Fv4%Hvnv+F7h&9c! zHQti2u_o>?gd1?1JWkKRVBL|bOQ5f;4LqbFFYmL5<@AvK*pi#mcoPS9W4<={&B~yc zk?)>E>T+^&SMd{+=oU3u=i7A0`1@0je*2aj6(#$OSGA)ifO&JUpr9ZnK3*QS=tpF4 zrpZ0{V#3IdARqQ;PS#x*Luz zRq!@iD0cK(Sy>e)@)J8?9)tZ0f?T_Z>oEFS~m79O#xBZzB zQk2aybM}mA+h>Hzt`w3!c1fC<`E0Kx9MG*JJEHCzQ2&Z_+TZQ?Rl88?KPyE1hCNdO zC&u9H{llco%L5J&eg0K0ttL>cHAh~f=s~KGnqVbqW+>C%O2(jvwEwF2cpM#xfmE;Z zhk9uLRbjtf^Rh3Luc?2n%HIIt>%Xw8CgI4MjsJYtziuAt&F1oO&#x~elARawpC^%> zt@D)~GmWPnJV=HX}S;;&`U+8G6 ztH;nViyo1X&`Z%yQ}^^NhaGqx9Q;RaS>ssi{BW&^>zr|s{)Pz|f~n8uk2rnnoTa6u zpWp3`)!q28Sc;*XM(GsTYU`Q4qTo01N#ioa1(~II|hiik5Nl8id^;05@ zZ}RiCp?U9~Y@E|wK5%BzgPonnhkio!W5N?M+_Nhu?_ZR~ z7gE-R9ET?KjF^~~Ha#gohB%N#>gUsj)n{?ybo6gM*iOwDRa{(*TX1r6`r#@Wk^d+l zK+DF)=HmKO8U9mmRwM?QBXxY22ng;&^-I~k{-z55p<`papa5`lx2LG0xL!PeZm~So z4Nu#0f9JkjtM;RLxw)%4+bEDGwk#^jx2MICf9$wly&rHO$Re+Fw+Dj|^aw``^8C%$ zuh*fjSDeHIu4RUoFIewdsE}$Hpv5>5jqbKFlFN={*$9keQsHV${Ta|IVUa=c(9$CW zT9-B%a}cRUu2Om%Q{gs0uZeh$bq+l!<`KL~6FtpTZ7-fN@C1JjL8|V-LC2XiDF0`W z&oBZ4mtMAfjfP41Vzn=AXOP#1sfdWkgR~%(N)p0mt&0~p*imQz*bA1go&&H~cXu!O z))0}hZTEr1($dn_$sjhFbT~cM@BQeDXNo5m2Z@`5cwev89fxE*^)t%xR$X{Yag=QR zmnJanEoxfh-&{dluUMBca%$Id`}&A?v6_u=#k}E_8b(+xp4ZiCQ$I6Ywq)S+vf+mr z`dq`N9|k>3p?*wq$R4_6t;+iPda^NSRrJivy~FF;+S*4?oM>hvaPbKWs;a4}B_}7h z7!j9@jM}i_;o&3BB|E)DuTY>Z{G1j=*z4En_wQ2@)wH#5rWhG!=QQsDVqrliTVSU)#aPk0dbS(bTB4E0ADEoXsX$;=)gvs}ywl|6^yd>y3m$ZjD!N z9eR5&54-M6!*W{4qM}KJJj$dYW2dB%QAb&K8gH(_zKa1Ii0}Bc*j+AO-eW*F$DQCm zRyH=8PEL6zPc;&-dFG9*z(#?^yL6`)73r?6t@#B6v`=PKq!$(zDt%`YAdwoMtVJ)q zw~Hg0Rc#r~u?p5bxt8HFUggksq6!Awl-|F851l*2Yug+!%I8m?MpNwAb5Wd`5U&!E z91?QCwkKB+HsOa`8T|e}K2D{fT2E@#<05+bHfpzSW!-jkq^GBktExJeW?0rk#KLI; zAj5ob<5O%YnJrIAPClyG(yn84$RVRqdi^q|y#IJ!`OEvHBk$(EYgL_#o7<^qB8Mx0 z@x-~1B`?z_KvH(D`Zuv7;+G3oXs8pe%qSjt5WW2g916O_!4&L=DMWC{FNO$g^zAR| zw-=!q^ZQPqozMNE+wW;~Ycb-AGjGg!59Kdwi6;Gl-!FYX)LUxzUFrA0I8QiJIt-R$i#L~hUxbA4gx1A%8oD1Ja<9wI^;o-2 zb4XnC-gX*wL&Mk7`0?2C34R3OxPg`0*4KIuv5+>UFji7}Y*#ujpJq^!-#Q!ViMNEs zzQ!zSb45$*g+R&F$Su^_98TNOeXZnlc=dF z8!c_^KM=vUV%FqIJEruU0pMC|dM0rx3<>_{3`{R7bAV7$5*VQLwmW?iaVs-aREINA zxz^UENpy*@w6@mJ(|bvDgkP@WU>NL`ujxR1G0; z(UoP&9TCxx2EZUg7?yGF7T6+@Wo>nYu?;EmMuh<58}7ra9D38XoqG-au5>c1Yfx31 zz;lD3Y^!DSnj&s=n@3JGMau!@Y?=E~o##4#U5cLHES{THs)auYZ0ESlogQxK6 zt*~}Dvu^!upHCVOA7n-b?~EPVd~dYzfas{G_=w9AT8!I7TqtG`yOI8mgF$2%lI|C>Pm$e6SBpM%!@{{Px#dyHC{js8X@Y(&<{CRNOJ3> zJko5QxD%5uJ?JRZSIx|kej>L_@XMVG#Sh?#6J-Mw&G{Zyqij;Z#Q{YrEdsvPToga`2%Q-ZF7=2|U8>TF>Y@_W(3xNPdU z>*{4vhpTjbYqj0IOkRkw+KMXa;1srWC*@ldtba*0Z8A=hQxu=2zx;g_x+N!a_O0bZ-LmKuj;s zsZ$y+#q2x2RQoQ~Fx!8~=2wmo*E)3Q(3i=?FC>C1Ojk@h1mD(DQ;WX0gZAjL$tnKE ze-}W5)DHpdnudd9(wXf`)x7Ynk(hQZod*dyKBlOGUiAC)=@K^&5578*Yjkv!4?Gmm zLJUqpDC`cD>T)UsYC>rg7F$We&e>nt_wZu#>@t0+W~exkp$5C7+cDGpT;8t2qM|cH z9yT_q*G)~`J&eE_c>er3&at_<8P0ET=j3$$bg|RqUk4CSCk_0HxUewr1l~`Op5}Z= zz_2x#KbJ0EoMQv`CW4UI(xP6ZpKccf_t1rRLKkYleb#bzpbVcQvwo`Isnfv9);7;c zm6!A+m}KCWb_Kn}qLk)}IHgc4U{3u5)@_MuK-ZEpGH6eGefxY+L|FL4@I0oMpifnX zmeQ@S=M6m9MtjP!C8%L#)F*fkA%W?>WIp5nt>Ox8l2HmV1Fkw9)62dW(EZ%!Ri?KK z+foivFO)p0;JaEMWokvqiJ_kZn?%I6N4P#x#t6ktXj6WVAB$Lb@Dq0pnbrpKiHYe5 zR%~SRTn>v;mD$wHDE;Aw_P~*m8(1*N+z^&baGNsbCN~^m*TBkxA4>UUMb$*@oFfXmyA2 zoJ+vbouR@-8?SM#H4t(h>n;scmE*=U66vPo!zQ;K%E@3S+k2Mrr8dx7Vb2X$8dk}x zz8tpADUn_2Y4Ce01GWE^3?Ts&UU-`Z|KU}=fAhDj=i2~qgwpXP^jS88L zrY5C_C*U3En3(V1t@>}*9fu3I#4EqlG&D@tvoF&M&b_m6Gk{g<+WN|D=etMy)?<3X zOuVTZN|3Ku=&iP?*Hgs5;22ZEU&%Hc6F7@t@^<7o^fmArq;$S{>)Qb(Ko-nPsD87E~dsrAFRN8j^ zm*L?(qq~!OdwR&756KX@yAk9(<$=5GRHoKQh2@v32l%0Tw7E``+IlOXS$4bfIYR|L zbfZFKslWhZ_ou}HW}{%?x$tt(qASycumFuGK{ZkvNNVL0tfoQv6SwYcGkEGBY|xOu zOo=0i>{UjIZ8o|g2aaERGu>M-33t{8x`?8@@$=_4uzT~|SM2By9Jmw=@4_sOtAW3O zack`Ay86DR#;0s{{NVDo`x>u|zyMM<78>iMUjFLt8@Jnyb&Yj_6_OAgB02#p3OuO*-03R=J!k_@O2-!bGa8#eV zkV4N*H^r9cd>Z5p_5ePP-#$O(O4Bb=2bO)5pI;ODH(1PhcKyYMGkic!ZWaUvG0CDG zFYR&HC3aA{aAMkf@~qvX6%Yh@`nC%Qp2}I}A+am_Tn~WE_eKeyB^8dUYSq>lUa|Kf zYJ*({m(QCWA;*~{?%eR9O5Q1Hnx}-%b@3ER!O>CvZUlZ&{%NLJ0_j3~3(j%zmNNTZ zzTdq>eZ3I5& z3u>KU;I?ixPJ8mRBx-jV^$U&Z|LSd~nk?9aj1uJh{f*o#Z2yA~?R6QR?!#T~1TQ=3 zFRi0q9ZE7`bh4DFQ(r06h8$-6&GfrmoienS5 z*KHnHCy<$E%?^s4&igj&A)0`zK8P?gG38n`AD-^YlH1ywHpfP(gw@s6q5f=sfTzTm zm~*`a_MJdasyfxGPoSLEO<7Itsr1I2+;mT#l7!1lTRC<;*Z8gXW2+JnzC~_lMqd6r z0{srvG+-RU?@u2;HfMM)9V4!!rKMdfbcpt$J{YsQ>Iw%HEQ54fcWOlT;Z=^)RYSm$ zk`oh~DmLmWP&Whhrl1$1RRn^uv^O8TtQOSdW-#(NAsOMwXi71ZX;Qs7@=P_kxcJ&^ z=|Z4g*^+kut=dEQxoTQ&MWqn_W@zoX5P?khm*~_6v2{kc40P9)ZAaKC4+|q(GEhxd zSC`=UFB36zS|vg}j%{OT$2G`Pz|nUV9}@zeIpDcIO%@sH`gLcTQH3FjI&g8YGkJjD zPIlR+AG67$b4b-UPC^iFPtqFTzJ8wz_(7siA%7bO^t040!Mpff6(0^NsTC;6qd-#m35kv9AUE{r2tKng#}G z%gfHt?bk9se)!M?E2pKa+dYcy>FH@hEl;Qm%Z6zRi?Jtf78fzxTwJR`5G1I08UF5T zE|!!M@RL0`=0VzN*V?NS?lv$=kRMT;tss`eHQIW5l3|-{2g?mQ4O}Aj91?%?7y3m4 z`b~S?1*df4sRHXx#|$#$A@$E-q|qRNaNM9+6AFJ~*)Uco|4U5HAuAMS(NQ{n- zu1o5da$g>QrYc=DG&==yZ*wh@STx|VP#;yXL*NDmDP8WdMommW?+MU*?Ni-3c`ma9 ztb_J|hmcpzd)a~Y2zYv=co!o3?o7dxlVQ+Zxp{ealnL_zaf z?aDk04g08o09dm*&q1ImiaSq*fwW!vQ733!b$=)7I}!RNJysn8Sk9c3JmiEZv>~sD z7uMQ3I^sWnz6>({3h0kzFNjlO;ggR7zS6n_Fq;F4VJK9vu&tvD$D`xZi;s zz747b47;NxP99%U$PB zpJMbv)xnTqI~OBH7 z5R-J64_xT?j`WD}y_^`@kr+vd54G}%C;8(uUM=tgBZ(Vh#T!5BHphS_p$I9pK9y6v)S>qR!q57bVw2MagD6O@ zXV$_7d@;w+3;~&R>a{IE?YJNrLLt!u%&DH<-Yx3@%p}lADUkT=<5ih@JR?tOwec!e z#~Pp5P2t9{4jT=fpd(Fh335@if*Cv@2BF`#FZSKLccds7EGL73oU}A$TU%Q~BLFbR zq#%MsR~_qBXmhxr zon88X7jYBk=s4D(1<{QsNV=m^QU+j|waxUGXkEQ}L=h$ZV6(k=_N;XjyVWk3JIYr? ztqAy_%+9F?St>zHQcG@z@b|!j!zhz`5Iu>Ew-v#0|*n&}iVz8!m7mk6+8Z z2)19KA%r4AV0%&3d_Y`JXs7EZB_>u2)jO7je$FHryqw`d@w1A4=P(iFmKItZ|h2_-g55Q-rCd@1;o%2VkyF( zg=4SXRp#gC4f6jyY38>))6eWM^V)r;*@uR?sza4Tkz_MKv9>@@k2WZAPAV%ifJJHd z-nJV+YrwVyf{Dpy1-2zGLPtS4{Y*Mj0Ezz{#_U@CNf#hGfE=7#@mkK!fL*(GX7RtT+b=Vt@{<)Z!CX z&($xigN>6*MOmODM>VbwUITlFvAP!?@Qcs>P1XyJ{L4u{jVv|9Q&F!lLjE zz8u^*#TM|XE+bPMB$hiddVw7#5iY5y>~S-Ij=;pstOYy+nv!D6Ybp1Ug|&_7nd>a# zn2W)}nDHj1$s-#OGq8wwF$fZSho}07{;*z7Y^K1;{>O&6aSRu9NFp9?ybRzWT*$;* zdTnGM6nP3nk#auN;Kv|jD#u2ykp4ETL^;rx%I#^<+_ z8?&7%s%f;NVH0sfn&jw8cQHQwnWt}!I23xA(@IR6h;E6AU5sb)vvj?A)ew6|#k{rQz?*bv;WF$#3R;^_&SZP1$ytEr-BIy&V>>Ur;2iy~3( zV}$ffH&eG4@M3Soq8f+;m!i(VVnX>e6%sXg`qtO`Kpv{_r_j%!lx1#dX#qKRtS(y^ z*i1Vq6kl;Fz`ywR?IvXD3sY_n$2uJ)OI{CmFZrA4io@IC9aY3?kDm*PSjcH=J_jU- zhj%Cj_fCH!QoBDiH>YjP7;|-&s?90Ed|_Rq5I8_u{wpAxNe!-s8>T8aCU+jD@2$Kq z3w{-=+c*aZ++`;n{GrzoQj!(u6>Gg1|HK#S?<>q`X@+y81kGbVP}%_H)Z z`L(5v8Q=D{w(}7Q(dZaRrH~E9J>%a8LTFIw$cYmL%AySh+!cNOSJ^o^PD7P<%V*Tj zgMZKp<^>wA0)e-v@jnQNy|XYrFHZ{ublVdSx&3Tx89;9#-Ct4*;_L*%%LjJ-`^OCH zjq_8*&Qo0{Hs^e0&{qnBA!OjSx#Tt$;j%HP7ZNvG7ivA;6wM@T@dTV4aBta`Ou$2f zJd_p2uYftKgUjEdrDG6R%j(WQq;*rW;=@QbfA5aZ@<{HN7dnP&U>c+AQU-*rc)o`D zMPy_K9*-}dx%RI6|617 zlJr@fU`3x-eXFWWT6>C)HRQ!;yk-_^5WlhEk(ZypY|_YC4X_rsG}$@1HWt~5`lFDJ zO-v+2qJnRub8Z7ic%)S+dgDJ1maP}MJe`K?XkNO6c9j%8yP^x=>jQRagLOgpW>1h0 z#*BQRLk%M64&kS7U1wV&RA^jNdYcC8kyvzpLnzGWnvif&C%O#z__+?wDBlKWF#od&q1Y|E)$7FJfl z5MOaC{&cEX=Oj;^Jc-uT>#N`{s;GMJvSh-18fU529le z7ucm1L1v)>6CROl3=(vwkq=L$ziSbP7S?hhwCD#Oz;x=SKNk$l9OdDWhq)WF&iMHF zPcT1oR*PV`2^~eqFn85_Qwa+1>hq8ZN!v!8Yzvu@2lM&IrQ&_Og|M z6krQ{0_kHK;wsRYK#ImiZcdp8zqvJBGjWerx(RSa*nN3Qdd6U6y})Z-3}ttqm+=iQ z<*{WmV=PQ=O#A-b zyCaYypo%bZVa=nb%)J1lIX=_;=|?bWU<7D=c^Vb)@_{S|=dd7ox9?$LoRhTP^8adO z8qr=$B&>^L#-9*BeEI}7+*4xY9p{J=7^o8LvKxn7KQ_Je$0fxUFcL@SUCCb?KsXHr z%R-#ungG`p;-wAi>##IG);xk;J)chXy_B6|5n|S;6a^RXe~vy_z!(X(Kks3V!md?7iOT_?~A zz_>96x2gy0!lNAi2HoG{$@CJ;P$iiy%%Cxg**)hXditS)A-(7u6AFv8(yLKWxiGBi zUq0*F-sloB!laVo+>`0p&eZ8zIS<(=|CS!b(4rIL_#ENnXqvs5^zXD830yo-DnexJn75c9w@R^bTfbaewUFxA{di7%qyYZ}0#r zy=|SH5cFR{vj$N>c>Yw!tO?ziUn&x6UM(NuSc*JwT!rD5g{rKM_s_(ZK=Fq?; zcoED$YnYTny+mcP{x>L3h8CKSAHIPyrlzF795F|C+1DAlxtao@DZQ}?_8-Lyrkl{% zugG!Cs;}_YwW8Z`=m<>N>et6og#K57T9xQP-psm3V;)$UG>K$WR&U?8#3^k`h-1wD zXGn!9y0Z%DFLzR}a{sM~Eed+Rx3aK*=J7Sk|2_j4k$(_KDx3!K{x&RrapJ5O3^$@C zbbJ=qnhlFd58I8!@#b(*+wgCVkr^Nz$I}W1)M1d~W_|dm5t#?@3qaO!hy{Wh>`oJ& zMQ)r=(aGq5xvFHSn5$@|1hSlofl_<`K@@4I=l5BRP3E&1|2Bs~ECTWZ(*asw2dK%( z-9rZnP$O8!z)8z!W0$-!8!RvNEjzYKz1kN?-EaX|MB-ZBSQ&41Hl6Tc3QinEPc^~y zOqEZTDCzuucVF&X&zhw&n*ZKxroZ4{E~veh+_T4M5k?6yWy<2!U((r1Czzvc=~A z{oB%!T+A=etanlMVE*|9zZL)uHb~)!34a}L-S@R~d(gR*--B_cx@@rkC7P(8H=M0Y zaS;(_MZ=M(PGMo;#Kgpl;0Ho?M^he^*OJ|%OK8|Ua=S=e3i7&+$5LXuNnCu@V>w6& zRqeiPo{7+hEyo%ndJ60#5Aa`Y8eM_}7!T8fp3yu%V`wmXq2F%;0xK%lU%0S?b`_O1CvX8q&yHXXeFL`R&^` z+yaaM!(7w!Jhsd*LtR~+ya9D>zymeRSc#d*0-FjYi_dtJzu%G}Rl~P4D@RXD_%?J5 z{>wkthtbK<-69Q#6arbYHmx^jiG&Q`aR)bdlacN&!OGft)VX9cvV>__~x8Kx44}GQ(X!2^`Y( z;-7tOF3ZN(QqYe-H8kw$@dN^jPpQDCJjUyelXf$#CY!6ye3+{7=+OUiP>V8mb9Fcf zoz&hOSCHAIH27+3=$=2)Y@I6fznaeEJt#&piQ0rg(#@1#e&!Jj2@ise#Pf4Vm~a?x zzkiykYk2)tR?pJQQyHL@CSu}o%|z0aoh3h