From ad3e5d1577a6bad8d4bb27dd5f642dab7634bdc9 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Mon, 21 Apr 2014 23:58:17 -0400 Subject: [PATCH] preparing public version --- CMakeLists.txt | 89 ++ c_tools/CMakeLists.txt | 50 + c_tools/analysis/CMakeLists.txt | 6 + c_tools/hod/CMakeLists.txt | 21 + c_tools/hod/users_manual.pdf | Bin 0 -> 92581 bytes c_tools/libzobov/CMakeLists.txt | 8 + c_tools/mock/CMakeLists.txt | 23 + c_tools/mock/loaders/CMakeLists.txt | 11 + c_tools/stacking/CMakeLists.txt | 39 + c_tools/zobov2/CMakeLists.txt | 1 + c_tools/zobov2/jozov2/CMakeLists.txt | 24 + c_tools/zobov2/voz1b1/CMakeLists.txt | 4 + external/cosmotool/CMakeLists.txt | 86 ++ external/cosmotool/sample/CMakeLists.txt | 65 + external/cosmotool/src/CMakeLists.txt | 92 ++ pipeline/datasets/example_observation.py | 250 ++++ pipeline/datasets/example_simulation.py | 144 ++ .../apTools/chi2/__init__.py | 2 - .../apTools/chi2/cosmologyTools.py | 109 +- .../apTools/profiles/__init__.py | 6 - .../void_python_tools/backend/__init__.py | 1 - .../void_python_tools/backend/launchers.py | 1319 ----------------- .../void_python_tools/partUtil/partUtil.py | 2 +- .../void_python_tools/xcor/xcorlib.py | 66 +- zobov/CMakeLists.txt | 16 + 25 files changed, 934 insertions(+), 1500 deletions(-) create mode 100644 CMakeLists.txt create mode 100644 c_tools/CMakeLists.txt create mode 100644 c_tools/analysis/CMakeLists.txt create mode 100644 c_tools/hod/CMakeLists.txt create mode 100644 c_tools/hod/users_manual.pdf create mode 100644 c_tools/libzobov/CMakeLists.txt create mode 100644 c_tools/mock/CMakeLists.txt create mode 100644 c_tools/mock/loaders/CMakeLists.txt create mode 100644 c_tools/stacking/CMakeLists.txt create mode 100644 c_tools/zobov2/CMakeLists.txt create mode 100644 c_tools/zobov2/jozov2/CMakeLists.txt create mode 100644 c_tools/zobov2/voz1b1/CMakeLists.txt create mode 100644 external/cosmotool/CMakeLists.txt create mode 100644 external/cosmotool/sample/CMakeLists.txt create mode 100644 external/cosmotool/src/CMakeLists.txt create mode 100644 pipeline/datasets/example_observation.py create mode 100644 pipeline/datasets/example_simulation.py create mode 100644 zobov/CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..2dfc7d4 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,89 @@ +cmake_minimum_required(VERSION 2.8) + +project(zobovPerso) + +include(CheckCCompilerFlag) +include(ExternalProject) + +check_c_compiler_flag(-march=native SUPPORT_ARCH_NATIVE ) + + + +find_library(MATH_LIB m) + +macro(add_genopt _sourcelist _ggofile _basefile) + + unset(_structname) + unset(_funcname) + + if(NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/${_ggofile}) + set(_ggofile2 ${CMAKE_CURRENT_SOURCE_DIR}/${_ggofile}) + else(NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/${_ggofile}) + set(_ggofile2 ${CMAKE_CURRENT_BINARY_DIR}/${_ggofile}) + endif(NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/${_ggofile}) + + foreach(arg ${ARGN}) + if ("x${arg}" MATCHES "^x(STRUCTNAME|FUNCNAME)$") + SET(doing "${arg}") + elseif(doing STREQUAL "STRUCTNAME") + SET(_structname ${arg}) + elseif(doing STREQUAL "FUNCNAME") + SET(_funcname ${arg}) + endif() + endforeach(arg ${ARGN}) + + if(NOT DEFINED _structname) + set(_structname ${_basefile}) + endif(NOT DEFINED _structname) + + if(NOT DEFINED _funcname) + set(_funcname ${_basefile}) + endif(NOT DEFINED _funcname) + + set(_cfile ${CMAKE_CURRENT_BINARY_DIR}/${_basefile}.c) + set(_hfile ${CMAKE_CURRENT_BINARY_DIR}/${_basefile}.h) + + add_custom_command( + OUTPUT ${_cfile} ${_hfile} + COMMAND ${GENGETOPT} -i ${_ggofile2} -f ${_funcname} -a ${_structname} zobovConf_info -F ${_basefile} -C + DEPENDS ${_ggofile2} + ) + + set(${_sourcelist} ${_cfile} ${${_sourcelist}}) + +endmacro(add_genopt) + +macro(configure_exec _source _destdir _destfile) + + SET(TMP_CONFIGURE_DIR ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}) + configure_file(${_source} ${TMP_CONFIGURE_DIR}/${_destfile} @ONLY) + file(COPY ${TMP_CONFIGURE_DIR}/${_destfile} DESTINATION ${_destdir} + FILE_PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE) + +endmacro(configure_exec) + +include(${CMAKE_SOURCE_DIR}/external/external_build.cmake) +include(${CMAKE_SOURCE_DIR}/external/external_python_build.cmake) + +SET(PYTHONPATH ${PYTHON_LOCAL_SITE_PACKAGE}) +configure_exec(${CMAKE_SOURCE_DIR}/run_python.sh.in ${CMAKE_BINARY_DIR} run_python.sh) +configure_exec(${CMAKE_SOURCE_DIR}/run_env.sh.in ${CMAKE_BINARY_DIR} run_env.sh) +configure_exec(${CMAKE_SOURCE_DIR}/python_tools/pipeline_source/prepareCatalogs.in.py +${CMAKE_BINARY_DIR}/pipeline prepareCatalogs.py) +#configure_exec(${CMAKE_SOURCE_DIR}/python_tools/pipeline_source/applyMaskToMock.in.py +${CMAKE_BINARY_DIR}/pipeline applyMaskToMock.py) +#configure_exec(${CMAKE_SOURCE_DIR}/python_tools/pipeline_source/buildSkyProjections.in.py +${CMAKE_BINARY_DIR}/pipeline/miscTools buildSkyProjections.py) + +SET(python_build_environment + ${CMAKE_COMMAND} -DPYTHON_LOCAL_SITE_PACKAGE=${PYTHON_LOCAL_SITE_PACKAGE} -DVOID_GSL=${CMAKE_BINARY_DIR}/ext_build/gsl -DPYTHON_EXECUTABLE=${PYTHON_EXECUTABLE} -DTARGET_PATH=${CMAKE_BINARY_DIR}/ext_build/python -P) + +add_custom_target(python_pipeline ALL + COMMAND ${python_build_environment} ${CMAKE_SOURCE_DIR}/external/python_build.cmake + COMMAND ${python_build_environment} ${CMAKE_SOURCE_DIR}/external/python_install.cmake + DEPENDS gsl cython netcdf4-python ${PYTHON_AUXILIARY_DEPEND} + WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/python_tools +) + + +subdirs(zobov c_tools) diff --git a/c_tools/CMakeLists.txt b/c_tools/CMakeLists.txt new file mode 100644 index 0000000..863d264 --- /dev/null +++ b/c_tools/CMakeLists.txt @@ -0,0 +1,50 @@ + +include(FindPkgConfig) + +pkg_check_modules(CAIROMM cairomm-1.0) +pkg_check_modules(SIGC sigc++-2.0) +pkg_check_modules(CAIRO cairo) +pkg_check_modules(FREETYPE freetype2) + +find_path(CAIROMMCONFIG_INCLUDE_PATH NAMES cairommconfig.h HINTS ${CAIROMM_INCLUDE_DIRS}) +find_path(CAIROMM_INCLUDE_PATH NAMES cairomm/cairomm.h HINTS ${CAIROMM_INCLUDE_DIRS}) +find_path(SIGC_INCLUDE_PATH NAMES sigc++/slot.h HINTS ${SIGC_INCLUDE_DIRS}) +find_path(SIGCCONFIG_INCLUDE_PATH NAMES sigc++config.h HINTS ${SIGC_INCLUDE_DIRS}) +find_path(CAIRO_INCLUDE_PATH NAMES cairo.h HINTS ${CAIRO_INCLUDE_DIRS}) +find_path(FREETYPE_INCLUDE_PATH NAMES freetype/config/ftheader.h HINTS ${FREETYPE_INCLUDE_DIRS}) + +find_library(CAIROMM_LIB NAMES ${CAIROMM_LIBRARIES} HINTS CAIROMM_LIBRARY_DIRS) + +IF (CAIROMMCONFIG_INCLUDE_PATH AND CAIROMM_INCLUDE_PATH AND SIGC_INCLUDE_PATH AND SIGCCONFIG_INCLUDE_PATH AND CAIRO_INCLUDE_PATH AND FREETYPE_INCLUDE_PATH AND CAIROMM_LIB) + + SET(CAIRO_FOUND 1) + SET(ALL_CAIROMM_LIBS ${CAIROMM_LIB}) + SET(CAIRO_HEADERS ${CAIROMM_INCLUDE_PATH} + ${CAIROMMCONFIG_INCLUDE_PATH} ${SIGC_INCLUDE_PATH} + ${SIGCCONFIG_INCLUDE_PATH} ${CAIRO_INCLUDE_PATH} + ${FREETYPE_INCLUDE_PATH}) + +ELSE (CAIROMMCONFIG_INCLUDE_PATH AND CAIROMM_INCLUDE_PATH AND SIGC_INCLUDE_PATH AND SIGCCONFIG_INCLUDE_PATH AND CAIRO_INCLUDE_PATH AND FREETYPE_INCLUDE_PATH AND CAIROMM_LIB) + + SET(CAIRO_FOUND 0) + +ENDIF (CAIROMMCONFIG_INCLUDE_PATH AND CAIROMM_INCLUDE_PATH AND SIGC_INCLUDE_PATH AND SIGCCONFIG_INCLUDE_PATH AND CAIRO_INCLUDE_PATH AND FREETYPE_INCLUDE_PATH AND CAIROMM_LIB) + +SET(ZOB_LIBS zobovTool + ${COSMOTOOL_LIBRARY} ${GSL_LIBRARIES} + ${NETCDF_LIBRARIES}) + + +include_directories( + ${QHULL_INCLUDES} + ${HEALPIX_INCLUDE_PATH} + ${CMAKE_CURRENT_BINARY_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/libzobov) + + +set(computeAverageDistortion_SRCS computeAverageDistortion.cpp) +add_genopt(computeAverageDistortion_SRCS computeAverageDistortion.ggo computeAverageDistortion_conf STRUCTNAME Params) +add_executable(computeAverageDistortion ${computeAverageDistortion_SRCS}) +target_link_libraries(computeAverageDistortion ${ZOB_LIBS}) + +subdirs(libzobov test mock stacking visualization analysis zobov2 hod) diff --git a/c_tools/analysis/CMakeLists.txt b/c_tools/analysis/CMakeLists.txt new file mode 100644 index 0000000..4fd2153 --- /dev/null +++ b/c_tools/analysis/CMakeLists.txt @@ -0,0 +1,6 @@ +include_directories(${CMAKE_CURRENT_BINARY_DIR}) + +SET(voidOverlap_SRCS voidOverlap.cpp) +add_genopt(voidOverlap_SRCS voidOverlap.ggo voidOverlap_conf STRUCTNAME voidOverlap_info) +add_executable(voidOverlap ${voidOverlap_SRCS}) +target_link_libraries(voidOverlap ${ZOB_LIBS}) diff --git a/c_tools/hod/CMakeLists.txt b/c_tools/hod/CMakeLists.txt new file mode 100644 index 0000000..ef85f62 --- /dev/null +++ b/c_tools/hod/CMakeLists.txt @@ -0,0 +1,21 @@ +include_directories(${CMAKE_CURRENT_BINARY_DIR}) + +SET(hod_SRCS header.c main.c utility.c sigmac.c transfnc.c transfunc_file.c + nonlinear_power_spectrum.c least_squares.c hod_functions.c + xi_matter.c one_halo_rspace.c two_halo_rspace.c + input_params.c dFdx.c mstar.c halo_concentration.c growthfactor.c + halo_mass_conversion.c nfw_transform.c + pair_density.c tasks.c wp_minimization.c + kaiser_distortions.c + tf_eisenstein_hu.c jeans.c + populate_simulation.c test.c + dark_matter_statistics.c + meshlink2.c nbrsfind2.c i3tensor_2.c + mcmc.c halo_mass_function.c halo_bias.c + nrutil.c qromo.c midpnt.c midinf.c polint.c splint.c spline.c + zbrent.c qtrap.c trapzd.c cisi.c complex.c amoeba.c amotry.c + gaussj.c powell.c linmin.c f1dim.c mnbrak.c brent.c gasdev.c + ran1.c jacobi.c splin2.c splie2.c ran2.c sort2.c +) +add_executable(hod ${hod_SRCS}) +target_link_libraries(hod -lm) diff --git a/c_tools/hod/users_manual.pdf b/c_tools/hod/users_manual.pdf new file mode 100644 index 0000000000000000000000000000000000000000..f9bd01d07bdc77831b2cd32ebf12b748a79a2795 GIT binary patch literal 92581 zcma&tLy#`O(k9@xZG3Ipwr$(CdD=W}+qP}n#%bF&=H7{!)xVfoWks&DBC4KvlT1-m zoQ{c}1%_*XTO(I9Q8N<< zQ!{>k7#CM(Gb1|~&+HZ5%UVEN!jFHhalPOHe*e+2601b8=AQG~O{`yH`Atp|0c5=74#|1IL% zG!Jqpt{*IlvY|M~Hum|3*4NOT1x`^mGCxA!dT>0T(ho7I$WM#}u+p|EIh~l^Kr1r$)Rgh(Ldlc2lAz6L^h< zO}jd}sKw7i(cTqFNIw(s3P~jg5c>Rns#szF2!pfuaw0w-cyr9dQYHDZtw}`n8cxeL z+qFFw>STmDW%lRBSTMLD0uBES+xi}mWWy@LsN zM;4OuiG(lxW?PWiG|(1=8ywx2Kfqo@A%-(4wB&8=&3M+EUub5H99h~Ow5lVj;QpfO zffPHvx{KyiBwBnaMwa*H>7JYxB^oG%Az^3fDb~8Xy4w(XU!cJ&HOs`(RWe?|3LO#0 zg&(4vj2b}hOj2v(-KmH(`7_)|JIs&`=&g3>tBq@Wz86Xx?D_C)kwzlQ zJUbfx_-9>1m;hFsL+b<8!yMT{J6%QH|Eb0wxgROI1CAbIx5jNDq_`iYmg)maK2i~h`4qQXbsaPVX83xPOKq_@(MaE%4sYTE8`GR?KY6=z|io~I6Aci_%K;x zB{;hj2c&+srb`@UGESjexQcEi&F6u1~ zXjdQ)CXbio9jP}JN&Y}YDp`}6U36q4Fpc~0P&OOP8{3_U&giu3AT%U-)8Fovur4=v z-aF@$^f}!-UX^{PKt%wNVoBzj&fC7V!beG91v;aflMfydkUERTfNhpf%%@MHeF**^ zz23$xnY@?jvExfEnO-xd8?zDOxD1sW--{Gx=%;{A6nk6@Ik2Ru1@gdn@WegAcFB>X z2~%Y03g@#|^QZYcKV5@-f$M}s@U^5U)`7XFAiI}+Eldx`rEa|%7q)M=+}!EIssd## zad)@|KjcF5kt`tJ`I*ZsFGtv$^#Ie1MKz*z@Hwk}nBIDLi!r95X(rJ-6Rd5GPi}lB z@O#7Nfb?67A(Zho6At8YR#x76>uuOESUg%TFa1;}BEUfQ!~+#Gcpdv7Bge|AmZi2N zeqvj7UVRO-z5~`O$;($I>2$!6R%{Wt!f5Za=+3;7Q$79bpRh) zB7@EiOw@v$*t0p+@+>xZ-2B;jw+X#bqtK{^J8|7CXUCbCW@##MGS$MhW(5inEdNM1 z+?jc$pggFeX8+!n2Ib)uZ$uiYU}nI=WY4sHpgRk@{-3p!@re-ret{Pp2m$md7GUp| zBSp?8G&>jQSrx&?Lf%1bwy;^qwlzy!=C}3_Ctb-hmy}Vnz0CKzC3~_9&uYD}PZenj zTD<8v6V7`lsmN9Y@AtJmEjhekynCePIore7weXZ^osLH`%J0W241zft2%N1dz>6oDeES9kQXFSl>#orl^EpTx3fT+xO{~{ z?SipmW>>*!1|KZ$LFh;=?~^pXkt)$03Z+MS4j8+aY;7VU^)<#I+8>?+<|X@nY`?iv z)8-_al%3%4dO;=+Z08Ai9moQT6?bnlf{6$8iV0ucCj=1|0TUGevd$s;8_tmcu1RfN zI~2fgGPOnVGYY`skY+-6+y(J5mEXG1{W*Qg*~#@dgwk#0QqIo39j)EfQjb5?o}0K31;WA z52ZIAA>D<0>p@whYdzhdCxly?nO6p9lgG zT6~#GfXy%W%rCVPqoaP*_2-G&!b(sgieic!^aNvPi zywoM2CppQ`vmvB(Y`J|ywPVS#`9AVIVsVQ1bohx30|gRNBg=pB8&TheP(*W~oj}pn zZran-#9|{d>^FX0lzL?*f*P_$rT);nC8YyWsO=B{d{W-(tL4I;t+{AZo*2)@aKt=F z*N|~6Mehs`ox}obVRT8i)+;Si_FzGfq__e<%+L8;nJu=TN>7e_-G71NZo{FFUDt+02BC6r~;T(To#rk&;n@bG3)U! z0SI5(IxXevTXoQdq}>Nr+OICE>0SQ#0m8XMZGsoOpP>qWl>lO@bfOH>#bidQMPvM5;L@zRkABQ&c0=0dQKNavX-K?(Qf^JZ z`4BU00A;r7l8(B_L(eR#M?_Lmmf^k*CF!r7!TybgL5)j~DGTm?w%6Sx$y$8#jlB+i+RX2R3sYQs$-gBYSEiB zTq)RH_zVxQgdp>WPpVRv?S-k6wu@giap>(`Yy$>N9}i=tL8_r4_&bMt?*x!uoPeRf z*bJ3JvJf|kn)`ski^hhGZ#^Q{)~y!}n-Ir&{6)5B=?5SG5KDc7IDtL0k@HV zT`t(Sq6K(4H?aZ#me%u-18S`fz`I_YkW~sQEVaD8heDGhF0WfTiNxxy?bZHh+=utW zFl4??PBW7ml4R_VuBgE&dzIm;I+pp&SC}C+ofExq>pi_Zy5xYw4Arw^H-ZUTWeIC# ztv&Z~ak&s!{e%Lq$E_%Uc-ke;0@v78DbEg< zP3W`OwliPb0Mx)=Yk?RvR;D^}=*?l%G4a0U-aKSWd4LBh>8_?^M_Hhci(pS%+n$(l zaL7Dh>U()Cx{L>mwTR#c(@3~@Qz1~Pf{>4l#_rWKshOojUyDXz+>M)Pld}G0;VRDt z>+(FhxYdtOL&qV}RZmsL)niNw3EdK#suS9L-Bzh}}Qv|M{crn|k87~Dgv3B&M)ha@q1pCn=de7JN zY`)gbtcv&Q%hSE8)cw~rh3^Y`0$>n_I0^=+rWCzR&7J7x;Uv9XXN`4YOGPs95SYNE zc>LjM_d37yn@UQ!(~$h^Y)}PNpT+ZnlZ%Qv&mp9v*0B(o*|B~SNWG@%e8^SRuEK9v z{X6qzdPOclB+E)x?8^^!4SKf6u81xi=TuxBC7Vm2k^Ht8G;C()9DZ1{^7_-$ukfd!#=+)T>m>5l>Ao3w=e?hf)J-)K;MI=p=8wvT z_I6@M<8@vdgfKz)u7x>KrM~4HG*MNBtROJ=^K6a1e^l|W2?5m~DC&2|RsB@lr9xM6j=AcXn0Hh84>vp8KxaPNlzDOc zf9J<_E=dNxoZG`N{0oKO@~-rdzo9k6OAH)RaiK|mkGehzoXV9*hX-wvnv{~S%4ojE zhisT_;pd|e&VfmOx<~(*`ECD&nSu%ILscj)zG&*LwCsDEA7kIUB|U*cLOX)zVleFX?TPR#^Q^9 zXwGG%;_~WDzcW_A%pz^1o+}O#C$%N_&k1i_TC;{Uo-V&=)Dn72;Szq~*(K{>-R3<5 z1|=j0Yv%c$(~A9C9j(3@c!5i;-~KfodHwM%p`ccGeHNzJMRbl1ThG#r7x9nmT1tX# z7!vzM&eZO`v~K7~O6+m?M8X46cU?x)_}9=wI8LHi?u!oVNIl>(y!MnJyG!IAOH=#N z@6S6gOyApghHClmzmZ>H$mohuGkeqj155uo{THMDOHlu`J-uEI%1c(`nkC)N|yve+Li6euJg6`ub(c{z?Aju^!2x;E*jmPUhZ{h zdA9vNd=ZkH@nS6CPIQ8iJpUMdta5W2bjmp7?DV<$F%HUW)9^pMK2?N9IQ>FxMK~+( zs4wvv-*H~qT$qHmir>^h{e@w=Vt|TB?@s=dxj+hWn^n*8@&UFtb_6OM1egG01qton4|08 z>`ft-1-1Ei>$ige#l3_HK}~n;zeAq-NIT<6!=*&Bi|HA1kNeh_RbB>-E#F|#fGg4v zhz7v?mWZdo~7?8`{h7FtX$&R3-Q6^f>r=MuOXPyZ`F3{w#X*iWkc*KaRJM#c}E0 zp=@BJK;07ZC>FK+HZw%LvdWy>l%7N7DX4>onVKs+VCWHH`(f-8(W=6+0mu8*nARi$#5&jX?B=X2>(lB8b|$f2#qPFUj@(?Gw+M6qQm$nmf> zl_JJTk(S~`aLu)~V~slFnwwh$Uu&QhNb^D3IqbOb{iZ^vvc@^${Ib z1ePsrrpFZrxxtNyQfx6m7c_}gs=Bo{Sr)u~gvXnM7|`^3Q0@Bo1uHa%sQzp+F>Xn2 z+BwgGAGyN?k=k~2Ktlf z{-`=0mU~oWPJ-qmDTVxOo{nAoStgCM5v2 z0?Jdoj<~E=WU?22AsT}qq4zJ9(%%GP_J`u?SvdIZs%T?1pX0TW*;K+kF}^hfHGXVd z+>L|+;O*Oo+Ijv&l#~LYORyx>vgblU7V~jVl^Yb%no#sh`b~!_s#s$!mSq!dK_ERw zpBAIeW+Q`8B#7MG@q|;ox7YQt+;9w3-aLtG83UFRa-vf7uf+A&{L6gIkHORWmaHkd zPIakLr);~j=PT5KJ<{SWygL?cr9tOD()b-LTa$7rf@94*HoYv%>q0`yIPlGSg7S|A zdWeO-Kpnh-sACod4RAf8lpbvpon0rj$&N-YW!gs$x*q{IoSFm-X&9<7q)VgP3)bF0j$cNtJf#p!Gx(YzS;(T(&K;YUw=g-?2S6AW(R z{3$=$^LxA_)x7+A!CZd*0ZkD2L;Qkk8uW+ix(oRElzT($96v*SLb-iS$LN^up?c)sI8ka8UZ|fnRJuta|9aKbQ(h>p%XHKExE$`@ zTJbp)Ti>DG`qYoJok1Og%%SmE)l{Mh#o&N`cfZs#Vjty}o(7nAOfmMcSWoFk#->Nj zd2}79&o+&PA-8|cv6n1$!TXnbX4E^Nugt+$bCW%e!Rm(^+{XIc22pg7#F!xHa?MWl zvW|gs+g*YmJC3(Dh#(f%AHSJ3OWC0{$9X9gWS>3V=)t(@hCi!tKQ~@5V-o9T{%O}O zXZ!Ac@k$hU>~Biy!5g%mzXKZ%bK{=de?WU1@XvO&z_E>1ZcykF)H|Aa`99y=GmL=J zL!7~Y#Eqsq_Ve;&>dQx*dPJIuCm1|mzYOWBSuewJY?7uj^|ECX)nlAN7#j7Loqpv? zXEagWYD7?334BY;Kb4~8Cj+PsxQ|?XHMUdbWY>yl8fIR}35dF(Gla!^_~lBoQz!CU zM0g3Hn`a@FLY*Qk*{p2hk8$j05kxqd$g|vwzN$t{S;EJb3WY;OX?9->>U2!a!&Bup zLVV(+;p8Hj9H`44&?4{s7E29rU=cn-m`I^f@FiND|Immrdu^F<>dM^y_!milb7sLj}NOd^cxMkntxJFmC#_S}yZW~>2 z1}YSo4~WJe%1m9aR150J5>hyrsyFF&u9zW1D#%)$n!FB67A*ri^ve21EmaaXXT1H3 zr{cP$doj`Tb;C_X+b>VBD>SXP z)Y%VFLX;aT9y|vV8sPH64C3TkXDd3^e_yk5X65=#P3?MQII>*keq@N7^Xs|Lju5$} z%0V+X$+LFw9ChM?odxWgi^`N3Fr0q(|J&+nJUV|_r5x+cMA$yy1luPadGX4$St;l9 z2#abGiHR|Fr?PkZ!*GDU(YC&-Y){>(PHl1VOTUA*KK&{8Y#i2hfBP5J3&b>a-iPj> zT%|l-mHRX_MkOf;5cL+6#z21QPEK_Y-|ZWLT9cZ}*hrIh@2V7Mx)Y)BV>m+ z=|w26c{|354V9wqxrXZ&GhGy?Ms>Syh*#?b=kx?%{~6-;?rLY?rMM-SZEl)B&T>z< zVhG#Ryqyc_7Gm257O?n|0V!kE_W~@T2>6pf>Fs&K!Tr!tPAdk)&nF;>jQ2AU;pNzA zSEF=~z2-Yq>T0NSuDO-kl<1>X4aXB59Ag5-#rs6fJ~XxQR4~@h{pnySgpgUB>)9=lOB%->SmuoI=`MihaW- zbZH%eBdZu`XJ}MxM{25x>sGZBxh_o^5}Z)};=c^XG_&qLzaR7OFVAmZMs4RIDpAl; zS}X5L$Hi3Q{jm!q=c@_?@HKRKwd_jdD7!#;6f1Wut(zqUjWaq?zH%B{hTC{^4R?%B zZgqU#jL`L8O1~eq(PY<@Z-GjJ-t0P`i^(w1{ixuChy#x30H$(=5vV5pk9A-QUvV+4bf(oE%i4dpmkN{U!$SIpQn;r z_dCWYtlTRp>+w2^gS003)J1k?l~g&_D9YDOIqAS25Tt3*4i7yZyItvWQ5T=V?#}3vAGW{sAI2Klh(R}MSJLxZYei!C6*w|a z*zRpX6zCgiEf;8$_es1du=#9K0daB_ga5hdN@T4Vy=~|*_eNB3SjZUV-g>6D%}E7| zsQ>(&MQ^9Hs%ybhZfT;D98H+yH)N<8YU#agn}3}t#JDQ$3Eq_Kp8guDAYWZRk;RbS z6Ggd|7I?q;ekZb0s$sg*Mz54VBtEL;X4Bpkh5V$VO$QAfOC`4m%lemE%+n_8^a=Y)9Dag>7~-De8I|GI`Pt#f?bOLtK^y zt1nZlE38OzZV@vu8Yyhobi<9EoMIrOLF_fs=bgsCF)@-c)1qI7R@+!nWUPzvI|K7d z?9Z4IRX!ji7{uKnzQE?^@dv&JrwG75DpXlg(*h{Me}q~)3=^|D6drvtAuon$FYjjR zn-P@Kb`U46va-Oe!ufMu*>JWc0LyPLo0PGRtiAw>L*mQv z&tS&CwTWrZbVP-yl?N4KP8bjL>oUid|W!XRO-YryPrQfOL<8G_Qs z?7qtwV=N}M9)qIN(?mB|(232+5%4bQ2e&k5(e13;1OF7#I-!8Nxjihl1fi+pbk+l4 zh3JPK{)n1wRfUilyIwh~bObqWDT{(!8T*i3r`KTA4;rrA>=dz02G+mrQ^k#J+`*JV zLZL#3UcA|(R=*p5^Rn5SjI(IUZm?UIQm%qWI+b-%k8*PkySaC+mC(oBoROXfwRCjl zC3p2@(SXHS?~SHPWG?LN)&w-~8}1qw=eZlwuynT$pjPqy|CnP zQ-nU=c^x;u7B(&7GH!=C<7F&31YAIb_K-p)SG#R_1jWi)8l%F=%;C z?*Qm+N(O{z&OEkyiIivw;IillKDsseF$chSRhX^HY&(_VC(xt-BIq-8ZE3PfCjQ() z{reb8`}dvEbv}Q|b`h*==xal2yM`K}F@3!s&m+!=w}rc0iiVVGRF53~EvJnX*h%_r ze+7s=LZQRM3ETHf+a1p?ZR<#m09G!yElYS$jy1DAOd=lXpKZ|1SpFGKp*)a92)*io z$R9k;pfiNC;8z^gGlxc1K}mbzOSU>agO8R^K7?>^(NCVxA?`l?h2};ZV zm>fS?r+>_lmjBVjLZ862ktYQBt-8^&r_Ikg(B|E=wbuKzPb zmnY;6aU+!rxX$v(Gml>3&eF(lnWHx?N)-S!qAJN6)>?Ex&93_8Oqs=ji+rRpo|1xNgZL|qVSO2-KW z%W|^4XVT|!7u;N4aE_ebc)RI_CEeswSIR)~@E=@%v8K$H?4U6r4&2#f@`)}Sa;!Fq zE(8{OGkts;q_kBPXEnK;htA`FI5GxU;5;FhNu8e&;cxm^jLDp-AOlvS7Wd0@?${pA zln||=?~?$%@QbFr4KhlYm*Qe;w(iBYGl@2#a#CK*-lded^s#{RF|W3r5ygvB*JR?O ztCHpTcGsGww36uel>jznBR(VYJ7sF}kxd+UZ}i%Zs)WpQukv)Pbp9fInC+PFq*E(0 zXnAzk6($o7jJu|8>Z;paTgj>UQRl9+Vu|vT7e6JSl5h;`-e1^@n{ZN+J>)Xt1Z{Dk z@8JiZhQ+1)n=(I1%8drj#piBs$re>1T20ya_M?7DTzfRYtT@l%vXxumHsL=!VP#8N z%04M1Gu#WJ3~T<;bZ1M^$Hv~CNzR7bqpQ9xjN8(eDB(+Oc6Lh$%UDkR7Em#_twDAO zfDGhvS~m_}-c_4(po~~9V>_=x+^6=f0$qIm4$-zE5tG1W=*^RQD!!&)%z1TwYjMr) ziCtbqw5v=GbX1_zzc*m&<9ujsePr4h4O&RDh(Y^@2(W*DM=gwk=;%l<^k<0!DfMqo zfEOC?@9l24QuE9Oq=N}Dh&(54dZ{ID^%n}X*GM>Lq#+H+&pf@5c@n+yBCUDjeC7 z_OI%3k7$uh@p7{0oIHK@;si^b7^NjyFK1rr9z6V_x^o@P)m%FlXy;5xbD4EYmNCX5 z7|Uxk9Q3rlDQi{6KXpMZ7GnC-FMx2b$_&*wz)-^M0~^P?B62`Zs`Ui4?y$~x{p{Q% zX9h)0fNuR-oGyZG3sIcyb57QSt#6rmK#DO;3f=u4d;JIy48K7ivTJGAfmIXIqJ6{o z|K9%F!j%?M3EBluIebtYI{`l@4PPqBwm+}VDO8n!%M6eSEN2VNtja|NIQX<5ek$%rg3_t_2p+#DqmJU6H)z_c=Y~L#(+!YRG znH_`AQNlKY<4FnP$(a8+?DHl&NVJ|nb>^D+b!EfzF&5^Y+<3Y0d|t1Tv$)@C`UNjX zkt+Qkmi)hw_`jIU#?Hk4e=(Vb^?zV8%YT!Q|IgfGkDitT-h`XqXK`)kta{P;lB2+N$u zoI&nL+`GJIiz9_c>buxg?r!e4_3!5+EY&e_@y680+oww_<8kQn7suQ$Ugu%-iw98% z6rP@B4`GO2GD^Yv>yNVIXWarkvWLKDL;5%J?~ zN^(PSTcvHw_0ujvQbRK)B+8qBb>}r#5e<>zVA>P>45;BDU>-4u`oD{ZYId0~^DqUD zZ?GfZgG19o_SkTG=Cgp)LdMMJ$`s-3$h)ZU8LtINd|9CKf%ulVihwB>B!n?Rz3#tx z1Dws*0b0Sap4emTwdhX||Jv4Z!xcE5#kol9U`l^0C0C*(^(S8YU{O=^?H-VXTn{)> z2q2jONB=75xX;KJg|VI~C%&=D$b>KzGk3{0twc_hF?B50Dh0Wvj}Z11-TJ5qc1V#^ zqZ-9R`p7_Vmcj1m=Sj373(_wo#Q(0Z+n8!{-$ffQ6C**C5erfXcVbWUHc4)_2YvtM z3H)A5Z^I)1#Sz|MT;)aD^guJuCf^5^G#jdq&_%~t(JoOTAfhauVS^i&A4zga?s0V4 zK#(PSqo+^X1nL;96WSI&=7*aVNtIvdHmTc+1vki*XY~PfM)@x{i}5(Vu|z|E-{<|UU=QkuTU;p!aFC@ z9!aH(417}XylpSO^8n2m-d1B(8zlab&mG17Tlyl5(F8p_TmiE|65PHWwk)L1BB!HO z!tAAB`jE*R#}y;ZRq`VX)EiRa1lorZqpK_>^6wU42<*~d%)_VyL|$wDGpuiPRAEESjRpq@uMPKA|g5!5=m~x<9pHgG*?6@hj!UY?2bh7;T`3rHe@)c z_+5nt4nd+!_iY%9gqO==WIBNmOxB+R(b3`{!A^4TufNp65k$*4B0zJTOriT)V$tn@ zcE@{piuezIzfhl075fOjPQ4q6qop8LPZESvLAZ{g~^xj$>pZ~~6 z3fw4Ltr=spUlwp){;EF&lB6rCtZva%LP0|jWk}~Q*^e_OuC)iy0&4;g999tc6)`Yc zO;UNkG=0!@v3F3!1tiET&VG{d0!*+e7tRY)PjYX=#7cznkh}!Bc<k({WNV%x+vIXF zE282C!(LvcYY9Pc6$jHLmUg>}xqX&M+Fy!hD>t5bPtY&|U@!;#7OEJh?s=!Hqw#fT z+YjSX>?CA58^yw!iQII%CwKAjF$A)}Zjyfo(7cBJ2KVczs}jM2iu$O4TOz5*{UM01 zSY@4JR8{Zv!YGSyJZ^LQOC!+zmx8fTsnCnBm5bzTl5O2?Aa2yg4Q@@k)T&*b*h;uP zRp-zXZY^x=dS4s1+#XO1&RgIyQ(kiq3$2)sztRZ1AS6U1p^@7y&!QsNt6?>by@B`E zerED;m5>+k;g8XM(0*D(i`gExrxYddnh`Q+6C{Bvouf&@Du2v*SVeS@J3gHO#3MgXU|5QK_rmjD64+UcbssOC0+o(rxZ zTc;qvbz?$-!vs$m5hv3eVW@No(q8+y{8qWUV45ucHkko5y! zUF7BEeS^G!yuF9BSTW@iduRvoz;wC!kmHprTu)_XHPw|pBjU{{1A{#QF;YAMu_!ix z8I#04XE}uGLz>d}*O~b-!OyGMSc-TS)Ktf(e*o=%f>eirK%O?Fpu3xxfnRVD_K!}& zSlSMEHVz(4BZO$*i8vA#WjIBsq5=$JoH4sZ){5Mk2W@7jS#<;poWz8~gQ$nSrhBWv zU$-2TV?#EI1ym19 zX3G_lq1v#QJ%}wB7?0&Y=NIVxH&8s1;rvo1$gXuFudbavo(F$B<7cX0G5j#(9?70S z`tkPrBgWzNiVOqAk;LlHHcxRNy%Tq8IC4h)5(Q_B37kF_O$MN)V?pgMU6XM0XS#ZH zEyxI()07rheImg@Vi;Rl7<8nFQ<|TPKL|^hNP4z_%wmD=(utD|gR=DG(1YYdjfKw9 z`o-kzaGrI`lMPU{u26Isdk_B+7Ff)7MTq7qIuI2wIsC=_P^It`%~|mDyn$2ged!-4 zghoHIauk53TLi8+k~<({(F>x5DT3vzDf%iLqGzPBvM5cTy`!cRda+kqhjq%oW1Oth zYPP0+3l_uZ(% z4~HHH)mj*4f`}i)JaVI7GqvMd?ksiIYh41@V9Zy(c2eRd@-sANLnw|jD!_#Z|6_0j zlQpDpET4_eTD{)d#RcSre3P|k-9 z+zf@L4}?`SZNe-0qckB8Hrz-HVw;k)UoQ0kr=nD(dT0?Ry_IH1k4KtFiGy+jwh8VXOIH=PwsDq6_FTgm<)>;(#AGsuGla-IlyuOF|{>l(@b-;KH(ZKq2yEL~sS2lBAM7@~_J)$oS(LMqDu!0d5- zH?nON8H&0NJZe|i;iF~#STK81+udEwu$z`vcRBU>@KD?9tD#T-d`zZiRK^BxNmImm zf5`TO@p2@f^u#4@7qJ7%h82)>Tvyb3~82V2(snBzpL=;pVtaZ|L_y50+64vO z2F#(9_pEYvBvkUSY{lLql8Z;-Q8b!E#!*bZe_RIxWX+VC!iefG=8k#mY6m}Pc3vra zxb^xU@_mt&`yqtbNV3ie}`Q0^}S?^fR%89gU7<6?#5EOEF4m+7*jtZN+Gf8Ef zvj0+fNp$#gK1DTSaA}2#Ues@W7Y#K<@GQ*?ZMpfSKpttP9XznYm8H@l(D`vv ze1PiH1vH`v0+Be6hMRgpfYUC~#x0&%2Phgh{*&zBntG(K2O&gdv}T$ziBT1da!L*N zT3UBS(gY{`xumYMJL_0;Du4z*6EYZ@1$WenCS4H$PH?BUQq;n#wB+#-syvY1gpNrF zT+5FMGdhKbuuH!HoAjRWpdizVbJV@O1<}#G`7Cz9^#xg^;__fj5YJWE_tOB1m3V`_i29^{M~)t8lL#lXdKUk6F6p2~Mb*>9FVDRV+>YR7g%_(VgB zwfruiL{#4Gv`V&6{^aLb6dxGc%s^by?i56Ny(@CdsgNJT?rETF8%7lg0vejr@nT7& zwqfK=YeX|?nZPXZYFO8s2v2yo&*Um14tW(9#>jk_PrWj0yfMCX_ytN9PBKY|QDWkqkUdmC}}} z4Bfw#WqDJaj~8N<_qSw^E!yAtOMMfdSNlRW)NC`1$%$b?N3wmnXF?Ti^C^00@rBv1 zUuFNaLpf`{lQqmo|K~N^L;cX1GZ4h$k5)2A(>HjSI>~S>4iG2*v`j`=$J|S|SK>aP zi6VJ#{Yy?DkRGz&lDN|B)s54LKiX8UI-i!SKKB5Q>jmaofGH zO{8$e;l8)d*rap~6Fg+C2OkhhV4PD{aLipLf*YO~&NmqlC`3#tn0bT{wS#{R44Id{ zecX86T&FlqS;>gvBB%=Y6<0RNpX(SEX;wIO2l(4D zG${xA!fls_eQ(vwSJbiz1amhg@AG2C4(&s}3Z#7o;x7h0i1qsJLU*jfy{Mdp$W}K{ zG|iPsG90H)^kqm2Qhp*)j%4nqM&wJII%G^c#=;bL?Ac%NAxUY`V#&LZ@R$G9KV2+> zk=E&NlI@4EfpNP4ST6MOSQQ;~Ao&E?x^f??cAyVA;~*FsXp8a=FX7-YA3n!cpiQ^s zfPhYj)SgnS95DHQzjFvsQc{?=bPh|4=AD2S9Mc9I%}*_2VeK5nvI!e2z8OAF+@iD3 z$~`GQTi0Aol!Vs2y=Gl^#f4^paFE*eqMdHQS!`rabX^X-TkBhd#21xdB~WFBRi)k+(_(DaN)6gN z-r=K9&K*mLa#Lf7l+aK=x19VynSzT`ah03G86T*F<8fv+0-ubW{9GRuxi@>CNLmq? zk8E3lE|vr#1m0(H?x&B*pdUNEVI?;N`cfY7mS=#4^H+RP$^Z z6=5bci;6WB>}P2D_njEC1@0R08ond@Y-q(7b=)_zY%6l5)qvL?PwZYG7=WH@W!uJ` zEYcXCT$NHilj3Y-({6{Tsjx_`hc4h}b*h0DY9QmsIiTheaG6*!u|k|=y1+O}@g(d# zx>e|2U*H*Rfyqs0q(GXMK@r$R>du09H(Jdh=cK9vzK+4XxY_40IfeObRzU%{57O1G zN&i}2;gTxBoEwn_rIG!U5jmMns?vA%S-%fqJZRd~B&`vdX^3}(QnbZK2vUK0HGoQxhv! zA(IB>rWM>DWw91o-%e)du1%yd_F4Yw&dDzd@(z`mh~3F!885W>m%8d2elg>5ojtvBaU3&o_uw`vwyNcpaxkOP)^C4 z`-dypmy46Clju?_1vZME+>GhY-&vjAxj-qMkWd$@Tmi(?uf14$&S(f%DrdI~%^^u| zH_88CL;>gsx$1^Dr_{-j*m=4oGCGUV&kxrG_2@q3OvF$9(D?nO+{OQ5 z?3{u$3)XF0w%t|T<@(FEZFJeTZQHhO+qP}nw(IuYH|}}aFXt_vG9n{Z#EKkqek0+! z7M^c3137`2Ib}s7;OT$`NQ~_hq{1a72PQ)1D_im5Tz74Z4uMHij#x?Pq~GZ^ZAL9~ zzRPm=!PaUv@c~azh&07h;Gt0T>WV<(*CO+UF zizr=>9kP#FFM(bk>MncTNeL3p?oT=d~BKE+_&C zEkD}EJW!q9n*d=*VZWLw0(sCgry8F%=;50DlEo8?1A;J0D3~9whF5zVilBw7@XcZ6veY4>l^9BM~sl^Y@ZGtoaP=<19i(&MO|AGRCs zAL{EJa<^2G;QRA1Pe%1;{3G{^B`%~lf1fa)>gg+I(p9OeF3pREtxqFL_?cVu^z^zV zGK|j#k*GBYl;+;bmr9G*aoU2fml_^oD;51n6cT85!dNs*QI>xN);)%>>N!nkj9o)| zU1I>(H$xvka13DCRK!P*qI$m~eiH1o&Q7hm)UENySQG|;emORh0<(bWjzz34h%NMj zW(ABA$sw0F*ATsFP6V~{PBk$78y*OQBKoaBF)sG{=TDS5V-{CcQlHExlB;!>HHRkh zQcrl6T_m1=Ak>h$ExzPrC-8AKX^reX3pxHyC6gLxO=x!l_aqM3x!|g#iCJ1g@~Fop z1qJe}3H-HA0DZ{61O~l?LP#mDHh*ETSC9|P91}&BVRvt#Y?%3B78qI$?O86PPxgRo z#LU6iYti`Z00HTHdTxlz*cjUmY4g2g* zCO4VF)CMmKv}O}0xrU(Olc^<9Y_HTI^DMe0 zB_~+hhrdEzoR*36mNG^i1ME&=_bi4iEE-dtCcXs8`j0}A-B5`vKn(vit zjSB8RRPJf0)NxjRTd0?9fLTrMxvoxmBGffKnAjqCT)ZDo48a zExrkc1>uTRnUe3UDB(~c%Z7Y&`$XiE3dW`S49Yb0RN87s@hQ5ErGcvcI%G?uqes-M7xp?c|KLdk!)DGJv zA*KOe0(}qF5pcTa!DGmm6$-47*DXhnmyPUyr`5oH$*yn6D%}jHfGRKW5M&`IEtS%t#X}S<(TqV)r+61tQ~%+~ztV zEdxnWRK+1fo|p356rodH8l_H?8GFqtLuQN>ywb_K4&OU#AZ`;`v(_$4bo>1NA|b5| zU86z*QQcptYXHKmdH$-FcNgXr)Y;73$r(_U>Q~hRO?4S(fYXYLS*L_Et=n=n(*tff z^wA6vh8YgbU_0tXhn-yg6Ua}68t2p#^=XV9fZuB)#LLcRhs--ysS1Tfpj`r$p`bKg zc(gLj1xBiM47Q5&X@Nb%{4Gx`SP198sEdYqy&b`Z&4gOujWZ^cNkms08XFIGxT#7k zQ(?pfTOxgKBxc0d_qt_1ZufVVo9yOIK+g5NO)AhgHc@+@W=Z5fU5%-4 znKM)=?Rn-Hq8Q%^VC{>#LNr(ZOZZrH9-9wn`GKD{CTV>!{^{_2-|hv<8N#)P`%upI zRzx_zr>4tu*cvB8lM&d4Pvl<@>_&MOHf&Gh!C(>4G|t#c_ZJBHCD*zmw4Ki9t!Vn1 zG5yZvBEl&<2wd1w)qsAAa*d0hDg=1E)r*HgXqQe384QlFoLnFVkFU{U1;J_XJ-6pa zLk}Rxj3}X#Gc7z~d_q*PehQg`jJUK}dILco@8#_q^0m2=$i+@<&&fHY0?AD zQ)ng0)}3yRf|uX8x0s`o`L!nn47l`WI{=EwVB##Q0wvhi5{wuCQ#EVj_&9H@N5cYw z%kKw9k`UB_nPNe;@w<=GCr}e3NRY~cDT+i4XH6hA zVGXO?47be9$1qkbRprB$%!lY?6&GN8?*t?OZ0aN{S$EGZ!N`}{^#K*OUE7>_87BF| zQkH9lZaJK3wigHDYe!Sw_w54+2f2<3s^y`CN%LUI7Wg+ME>gwzuf-0HWG+SY3KFQ3 zQvz;Z_Z(2Kz@P<3(L^3hNKDb%05n?GvgT)@OOQCct=pBlfD2>`8f4Fdhul$ zUV~POV7XfT20hmGX9@*aVI~xW3WI>i5nFE9V-DF#J4$z44BYILs|lAEao!wG9B1VU z3n;2DOukk+Dp}A`zt3Vl@w;&R2l}GcU&6bQF{_H1W5D~!ip{AfG|xvEm=GHp!WMq- zA9%%1L^zgm!*6fw3}FX-Cte+hZI*UQHn^M{8EkY~`mibJ0(@Xj>xKwo*aOOl`EqxD zV7OTTbO#+|w{4O|xsf%S2W~u+P*upmtsO*Z%&K%@|tOlK@J=|)?> z^k+?{_>5N0LejOo#ngqx9!`hPxJ1+|3dAMPj4z9`&jqlGt>bP9M_E&T$Fj5vgBUSuL>mHxy5^X=xE`?6zmD(t7hzE074*Gq7Su(%Zy54G9mi3^6;oO z+ac?F3O-UN1TEAgE(I0HnaVvzRV}~~I9W*pelYbANa^XnUps$GHcmmNqBD|z(1J2> z(1fKM4UAytg@LLLrXd&$0(#GW&G3igsbtrnq%e+w9>|zzeldk}xV957-FflspgP{S z%`RjSUSz)lHyfVBu)?#cDuA|)HZL6WYz@;(>i!n$!i|~#Y@Mjs{XmIWI~;C0T?v>9 z-%0SIWJ%Fi$DZlHPOEJeiwALBM=UtqTP4B*`HS5*>8(L(*@|JRA;ca8R?&`BD-Flu zwF2I)wM<1fLA6I!55iX6&t5(c(z&FpC?jD`8P!I(w0h%KH{G$iW+#&@O7en7xszL! z3I)xoi2yBQspRB^EM`Au00P86 z>XZtoV*}be{=^zVcv9`+h6r#JZpYQSQW%F9k3VBj-Z2sgzkJ-y(OAnI47Fpp*a zd+%h9Bh*afln>bb*LvkCneJUFpsC`1jxlo?dX&3<;`~`BsadI#***4D_|Zf#kH+kZ zHl>}Id6`7Io8%3kmRrUM-QK2V@z)u^8*T0|!7ox~UXuKxVui4WBx#2=mfa(j7f+v~ z+VIo)^TnQ@}IYFVuKd>#l_wNre<@>=o)~g4E&6mfXZ4fbU-!9+~aONHc-}AGNF9gNFHbfdLl* zI^kW4L&y%Us&-Zp^>)fXxS#uh)EUPlc4F@_wYh_8t1h~H2ym16Q|*>zDS38g)!u=7 zjk3VD=TzD!t4qRLS(tF@?Jq}uN%pIK>FgZ8?s0nb9$~Yb?5D6GzB8FdEq>uSUtGy6 zvO#U~zJjN4YTrjfzzO{tQ3Xyi`!*mOo|>!u&fkCu{zVgQq7qCD1Sk+>chuE!S)A_Ws>KKmiGabmy?2lepSl6YQ(+{+l)cBl;^ zltM5=x2`PXKw&j?A*d&j`H0f?%Co2$+6W%NTCS~yV+`kd_g>mLY}2Gnttz3&&i>_F zTtbMglJ_mnUkz8a#T?Zz$EumS#2hT0!6x!|00mc60oldG5k;t$mK~_|eEAypP`uev zjDt!S)nJAathuM`2nKIS;>R3zZ^!(r{j!(@RLaN|`!-z|pp5pZ@zgoYP5-~?m&uu1z^@XOJ9wXFzGh`*bRKkT7qu%|m zj!Ivn>7QHHq{yA6nVD%sYz~H5FmF7&(}(o;=yCcWj2AlPKuI}fO8SF~t8gp+=OSQd zH)%s1>NQ@!8}(Wgf4UwA_C{6SW#QH5YivJmq@BUH|fiR6Ln?!wwZZaW~l` zXKQzs?9@d#jRT5JQRmP1?2|xLUK!m7a=zmm{}DhOFv~dWA;vnW$_X?0WrKRU$| z!sOCq(bx1g-eeDQhaon0a*RBdEVO&z1R==e{m{Ulc^)%v7kcApI6!M$d*YZK$*>|0 z4r)Oows)2fc9&NjqGNiR%;FTXWM=dd59_<5$*QZaF~odEO2wt+CdYJ6-ooMGBrpAQ z83tABEYJ8*BG8(B8kl`9m})07rLqR3bw)$ZN%*_nl{aQu>*#__zV1S$@<%=&H3$ry zBsKcbM(HdK?<5aLtiP5lTi7(&+G^NByH;`AN@>be#rP?iG94l=oFJDxHf!)>%U$^h*zrhBPck#^e zJx{(Z_27?P?OFiTDZ$x^fCT>=Y}vk_jO3url&MLi!^?AI(pJF+<9qM&2MC(lXXC%X z?|*>qe~RhMj12!D{4)J_F`bF|e*=E6G&kb^h$ea6)Ueqsrp=Pyd(-y?&|@$LLPl<; zLIvUlF1FLftBdvH4y;UKBK1*nP5o1crYWUyR@BSnEW`(erP{J3W@RO%WO>whO?@By zN66cAwi7~2=ZnMl( zLTy?a>!6^_!Lirz|QG1n=={Plz4_V=f}hgIk1DeDUita{lYRlbMw9%ax3fuD?o%nJ*8L8344IInm4WZ zyfycsQ`Y~Y+Z3)r88_$TlJsEe6XWKZu|gz#=&JcEp_nC8?oVgNguG2LM(#Kz(j)UA zo~M<1UBDe^ZxCq*cwl8x-$y%o-kgq^gH01b%_~&t2>6}#OmG~dJtr#S>rLacSg9*= zBLs$e``4Rw*1BkU{?XZp4|wytv%(RAJRqUGg_2RED3jAeMb*2m!V6#$Ywo&O^tDlM zvH_Dd-*%)A;ZfJFvsapeuTmMQ!LRnonIayl0r^t`kGoRCH^NMVazkw5CTE zcmnhtEii%zm==zk#~|V_tVi<@W`VkzCtQO~t&=}a6faJP zP7+S~>nfy0_#f0ix>7?~3B!@GUaRy1q4bN{4_8?HK+-z~^l4M6a7O(cCOBtcGHQr_ z;TYa>KEk0gQvhTE7?1Gr-C|IczEMi7bc2`U#D|-G=%Wm8F$h=0(zMbeQk%YT-32L) zM5+aTuuzzJq#xb>t3WRrd(kSf#cGc|qOAjB_T2pSbMV4!MccSr3DE`7C-_jO&p%N1 zcC)v1{t+gEqlGghQldj}y*q_O?pC_lEHLkk@BA8haR0swe4#h)Iw;NH!B zYz($guocmb+ueQ`s)cXgr=GM*`=aFSdHiDy=>eM85Y`KpB!O&5E+^*O0T6Sq zn-2+x)oL&>^HV=>RVfIN4^2f=F3M#VZ_yBvGtRTI zIeZZYOUKvxQyZT$Irp0q1x-~|3uXZs!2ZY1rA*=d^=}=^{JwSwQ8&|KcW#?8T9xZK zw4-YySlC!}1-{*ouU{OQ_`d@|c)^HJfeBdT;#J2k_|iU>La4n8z*WAb!{w`89`*Kk zCi&A+_zFAL1-%yI)-Tl@Z&A?F*JO7EU@Kw{Ti@kKxI9Y~%U z;y*DHA&Yd)NP(O{phObOaWZ72o*^>8QC%VtTlr`JX3LFs^|1~Fuy}VPtYSkb7c^VH zir1S)*$rD=$LwIi0_-K~{&sK$R*;%}az1S(LODmA8J&KJi?(?evBIu&Ft%wkW4-Ez z31WeKZv%auz!mgNm7@FmVUuv;R7a3A;Cuc(aV=7JnnKRfI+;=@CTMJ?i}vUCpIw0q ztGKFtG_QYStLt8S~42>7p$8Psu?)Nc}Absm9J;BkxpUojMHz=vo6~wBD8y?u9_=zA_y?E7xj^n<-Nr zbTvQ+UGm5hZMEW3)>8Mo->|6e5a_mPd99le)vSZ+*BYrYJ<5qjrhEb{WtkPrrUoRif#y+fqSj+$E}yumVnw%O zF}A);qJRvpx@00F*)FzUM(%wP9me^111O3pq>ll5ADBzV=bS&bSAyiZVk@+ph1!Yg zm%Mz%Zb-VXIp6$N_(kzprdZ?!%t$d=tzI($ZQ41VYNtabm%kF5TUd1*J&O+vjS+(8RHmt`Y znh|My#yCC;w~Jl!f{J~fJ(Rtv0oIm}t$E%?YaC>?r;i^B0%*3c>d7TYOFS;ZBQ? zsIkv68cpadm8<*3u;c`##}edG*#t0mOUi04&t8csrp>GO$!PKawP89luFQ9Fk;x&c z1T=FzWYc0Lt2Um{Y<1DR@)pXXdO<=5At0z86WdkH9AM2^r#hV93nMj zhW?8KF@;WYab8Nk-bekT#L!3%&%q6yo9>l9Cq9SP+P48+CSS&cw|E?}-bpbPkD-hK z(Z1=I4x_l<#Ro|<5jgKaMY(BRsd(E@1d?P1(c4{aRK7r?&D^1iEyZ^!n>av$1*=8E z`nacTxK|EiP(i995StG-@!O%h%P7*Gp@P@^?Hiw*oPUo)7$K!TP~f&o5X(lXXZN1 za=g-^wvlo#*E8*MS=wH$%>VeNK0Qh?nkp%j3+J?T$1d``-wwzvDt~*_Dk{|LsE`Lo z^6(YO+;*}gwt5)|^8HKpH@k$#HG&e>Snv9*Td0Pl1#Dk(y+i(}kywXDJTYo*Ag%z_ zwEWG1xqa4P;NNt*SJ(k*TM0 zDyBl9RV9P>nq`PSN^g`#Rv1sd6;h;*#Xc)GN2wJ(nF7>WQ1KmY{RK*H;m|m(Uy!wp z+k7EKBQQc^9p&ZhrRc5;x}cGhC=eUe%o7i7@nF*qHt0rSi&PpYU^=Jl_Rp(K0~b@8 zmToIch=3&uThOMU49mC)f8hq1JV(m5P+FO(0*vN{fbPBU6tqAm+(c_vR zcMd%BP zVXVKkrFgcjs%kAI;lcqX{Fi8%)9!ir@Jf^TY6+FNko`%qHpD?a%Z(%brQFM_O-k28 zP@gT-yw%PRGbzFIac`{mtDwq|XKdF}gkkSwi3gJv>i9*%dd6i|*;i%9U*}@OK1FMG zCY&_*)QY57s_0QiDZ!!n-!YX8crINq%q^)Q3O2v^@UKBnuW5KYT0dJpNU)GXe#6#s zOpB$7`^9c85gve$Pus=87N}Zx>ED2=Viy20M3`B~roEy$WZY9Eu4Bq7HMt z>X0r@#+b{HlMq45`fg9t;zrB1LzATwh7n>OSeGhPW9A&siPT8A)QJs3|3r}Hw`ypu z*A$aiS{Q40(N?qoNWsy-5EFDNNdk8y(cTB`f|_j zvmZoJ)M_@PS4r&Q=8f@c_hgMuzLuiYdI>lS(H0v%UXT|eaI5R6Pm10(Lpwj*;3mee z0;W={i@d0ppeqV1@p7!>D%^k~{Bj;a5Q-`ERm1G9-UYrr-ONqAYVw)Rs+GRX&bP3? zOaKC3YVo2hZYB;;YlBuR1n{NInhLE&(rt|aQ^jaBk?Uw4Tjd{n1}9Q}p`@MSAI|n* zbp?pz5@>n-QP!EACty#Lh~NojCohlaQ(2U##f2&AHNY6O=xQy)N7ynv z>L9+&w`5-70RT~>w^8E?|G1WuQvkReG3r9id<7qDgfD4&Ie5;uKUk+clgQPm|4Jg~ zOVPgzfebMp&n#HMiu>1mSp+qEImm)eY0%pfLZB=wiheZfjNgWT-~NRyM#Rn;M& zLhzgiw{WWKAbRr6l0eIMr0|!$eUpU;Ze{nH;o|-rj=6neHKmi({6$q-cF8f8xQ)L z%r8&9=Ber1nUVUv6~)W5 z(%LB&M(5u(pH~a*zwZBzDx)_{CAr76YsvSe+L8!(sYxeHxy)|fE6L(6yBM5K#Xf@7 zzwqtnmZs{o1WhghG`M5*lO8Q%_NJJ&s(yaQL>4e!7vmSmccU6^x=TQSj_JM$$KG+e zN;E&l$xvaCate#qN{7Q(awEdQ;zi(acYc8`XIBR}*D2X|*1IJTvA-e2a@vE0`GEwI znzY7itxiEA*_&(VSWO#m(N?IvWX!*rU42$MdN()vlsOndzf=iOiW++#C#tPo^szBu zr-09K4QSnfin46$u`GO=#j%)ZK@IA4jP>ZEz_#5GLu`yg_lx6b)E6yOpC+eM_ETq! z)GIM3IuTus>+r#_1;dxS1=fZ`PEsWzTs3>ORdDrVwV~&wf^brwcTMKJJrf}1w8CE$ z&do6kt5jB6Peg6LdjLdHXy2H+&)26Qu|50S5PDk%+%hTTw4+5bUSeRv7=ye>6-NnC zHd+AW*9E=ON(EnKGPOw!>Rj1K3UmMTy#^vDkNuhFi%C)fnf)3J$x!j@k*dzPKPB{N zWTB7pMg(_X=H-g;5LHD=1}n|U2JZ&k)#S|9IU~jBc|~mgOf$r%FSsFU%eZ=G=o@^; z%D9*2>O=3+eQI1_YAG1>T@c@eJ$1;Lmit93B)c0v%h@+UxR%z?+kbA8b0tXJ%r~8p zA=JyZlBZB(V8+vvO2<5phO>Y8dcej*n`Ll41ch4AjbiR3bc01qyV#*oi=Px&Q?hrA zl&7}uT(Zml7)BW|<1leowF9vjSkSSlq1LhfPt3>K zeTcWUw<+w571m12kGDXweY(5PSMAJ4j`pSJHi-HOkbXZ$J~zA_uWOj!w@Xt&Y@bA* zKTAU`nDN-Z?|ENx{(OI9V}rZa2_qlvz0V21tKup$Wx&1f3U4m%y{`>r9puo|G%jb% zU!R`;TMhxF6xYXjHGkgsCCG*Gw9PueP5SrVfjvLqG{DJ|ICi;)1}#*K5=)jt&(_os^RSs`^3CF$A?{0kIkjIy!{() z0q$q_k^hvvxH_8JXM9smpO$gI?H%rqu$KV#YFN)LNa;i;lo zn#8hF2ElZz>p5|om>tc2DgR~XoUpJy=BgcCQhQkbKQ?y$DBk@956`a$F^v`K6u+Vo z=RSG1GL*V#0z0I424oNx{IsPRmk;ij1~#h@ zKVLj+v~>tM&5T+~iE4@K=DPiA8Dx4InJ9lxMp9?@OmmUkFDihZ5Bz(2i_eW=*`7D)K&TJJ$xci~=Y_dP3X4kdXXnqFB6ccc3NB&cP9vk6Pz)8% zM`CI6s>-j!NK989O|78%=ce06t$n{!W%;!FDALxZZm15-ijLxU@ z3vrqAjdp&Qg#r))_du(lxm>ch>FrqzkQ0b3oN++F=dmt{&Lkus-kvT8EC!@j0t zBIk5OQM+6<^(4i1i!v>3pQKi$%LP6b#LQmH{ItZWW9>3BeuKAxB&J)GKn&-%SD%y> z9YF;Wgjx_}kBL}QO@ z7Do|v-?%z>pGJ8d$3aoP0wg`HWSm+NJxV;53MEFiuHd1vWUE#h()TxqS{~!OQ;>-YK&p{4Q4#}v#)t1ay+|SW83P5*N5@? zox&|puNCU9%O`!uLppR{_@K3>zKJs9!-6>vO1T1=!xr;lVz#%YdqAXv$kQJGzB_tm z<*|!+R!8$c&rlu0syTU1;~7n*)WZbD)g9_*{@36WK(InU9gqJ&nW7FMn+aw9mLEji931baBx^Z1!(&?* z1#*uGDJx)6`5yX@nZhjv&?BCSE|%lH*Dy6|7-`X@9qCCmEe8KCupZQ2TA z66~>7jqdC4NHH5&K53-n=}ynfQycGQdBkL)*{3vSV*_aD@oEqQT$lW8se1v`m2cS| zK9AGJa2!ZVp@}I37d=&G#)&7vo#uRSng>LIry6(Et5{HE6gx;K>;;7M(%)=R9q~$g z7^Hy9$l-Vfy{oth+hkKB#qQ|M&MLlBug2FCxD%p5Ev~QM3(&1u#I<##g6z7f;}cl} zo(kphiSA9#^*Q{$Jw~l1BPFnm);c9KvfrSz!&s0cw6}XA5$C3(rsaptjygh$7^inz z%x)t0J^5Cn(ke)bZc^C}l&RIdSze1@G)j!+)%yLp|86r)7?+-K`i@dK&u;tDDmm2B zrqEo*J2AQ-#EoR!mk%^S@c1m>ICvc3D9PdQC%vmXAFpzDdGugszl8y0Yah<20G3Ef z?qQ~`=~)1hc}~%~0ps>!8|zup-?6XNLewSa#nRH25>Ec2vTjfAU6#wD82a_t&iA4A zdU}GR1R7gGrnO;iwNze%NqNs&bd)l3$ZU}i-j+R}Ik;>3Fp!|%QF0HNWxs3em!F)H zE3Z9v?ZDjV;+DoAX41lg==203X1$Y4XxD_CKxthZ4*PzE&3cO+ zmPBmhnglXT9jA+dwHFV-ANR+2O_En;v~pI;-5>C3Y*>t=bO}ab*G6$hY08MpZhFw= zIP^zuw~&G7!7Os2o9#yYud@&ad@!t;I&lK3SwMVN{5A&?1u5yz*e-~0PKymrKvcEB zvg9}6{d~2*&kQiz#VZoCgG^QS_DaGy>LPCFpBOfvVZHHlOQ`&QI8*p;0YH}81kg;r z_UViu23QzLQxOx~mIAS%qOyFS>rs)p1nFL^y(jv?b|k6_XS?WlPDfb0lM$H4=sc-J z;VEOS0zaUKB<_uVDJLHmJ9@($wL55j0C>m z9IKUr5A>#JICRvz4eLN%rk@PD^klv*K>Ms_CW487;EGprIqh9nL9kibrzLaQyJYM{ zpWh1*$SijGmVUKf8Vj9hAdMqkB_-=}g1mLiTk8t7;&Ng7@#9G#WV1hF3Z(jhVSidM zkSv3(bsj6XEiD{91uV8PQWiy7!pT&b_QO%GTy-Ss#V!XMxyBBuJ}3D_yQvWzVVQ}= z57>QTos=FBzbc&OLFy#3EKI%Q61q8!ZR&T6dx1ijj;csgM_^s)F~ei!`X)lU*>@N5 zzHH05(Y=;2Yi-86%4FN~hDCP2qzv~i_qPo8CPOHJdU|OEnXm+3Z(byV*wt;NOkPZF zz-a>YG2^F_!lf>JtytzeF~&D!uy!WxOAuf|S*{k7hGnSDMyq+ z@|H5b0PGFDO9B{dCX~XI8uy2a@Ew{Ls*d&q8pZ{f&(hJTtPSSTehm!%g}H1>&RZmC zjZ_LSJxV|va=JMPw1h{Lxw~1nMK)J3(+a1d(A_JmU>>YeaXZXpvRU*v$$tU|M(LtC zK~NdC7PM2f^T|LgZI=GM}Mv^%=3zs`tdYb?)!nOFu^_5GLp zFeM!AI$kBa^~LP!aLCA>kQ$1yRObjtV^A~imV5{fSVCR(zk2&OM9el|}FGuFeQLU;N4kha7{XtK%+D0J>HehL8Ajm|yKj|3BH zG@C_pePUWF2SP>&sx8a=Oxt)l3xkBR85Qg!Mk``xvuj~xEps?1q*~)f10FD2E>J00 zZpc~UunIWIEK({Wpe>76cOezV?d#4ihji`ou0W+`}9wGFf0f zfr=IJaN13S;LnBo8RF>l)P%)7xxMmv{(a`SwO-5;4dK0T)^77OYa0@cW=2_(wkL!| zg}xg2!+2@3&wVg~Zh9zSBbo%X%vcj*@$aKX#LdF~BwMEG`uycSSLydjTv*3y!VavO z5BLLCb|9)s7l?RES#AszVOw$QL(bU!GT=fXlu^I4$c)WtH zR2w))bdVH5ZnA&dl<|SOI`JWbKum}=KH=I^2gRvMBE#vgfEM4fs*X3nf`1P}P&M|j z2VEZ@?E`5NIYWafv4IP@v}XzQScS^(!DET*R(G2Hc!fC~+zR0JvPtQ?549j2nZ-{4m1WDhScMJ9*}`%5LPAA=M0aIh6836&Z(-LDv&B9_D=?&vm7l&+b&-e5SjsPj^wD5d$FeK zR`S}<0vu=5j|l619dLU}ISmo5n^X?uQy~p8fGKoKe0vCtk zzt^dkjs7{^kW^$3ZP=NW$U^p3)}bYNt;*-~pc}k>!bn=9Wbv<7T#=0Yk%m6RxOO1@ zggZ3VQ{mrfDl4JgBi-gY5l3oHW0VA@Fg?Ud0}oLT6gN_pS%GswL9C5dclEw3rv7cb zEj$b}aY~glJbIS^CJ_u}5DAT#LRDcdxwcaNY(>!U`2_ufb8^_E>3h)-2SriMI_gQt zmQc``=Gx{cV==DfEoXFUyc_T?-CS2b2dHrwVVpcvbvz2zjBe}J-CoK#bTCu>Mm(eE zaLdjXN2M3EELeT#ZnbOojX*>5`T0iItdKXQ^Xjnl3XhLFclS}kUn0pQGl_U;~&6fDX$lsTN|@{b@bVotV1O)Y6dzqgupx?=+z(b z;mb1V1A6-A9wJ`;#rWQ&la%Rfx52I56fSQ2M*~FJ98T3}up!3PXE&@R(7~-T8btYB zsNGd(PYnhfhmx;I+P+5wpA$X~14||*D%QFv*L6pV-wtl!H|146dhTN1qKkX5MT%&O zBA9;F1Wd|rNN*Z8IRbUvJg@=X((RmFj(8c=t>f&jOnNx-6!;9F#%O8r4E?G)LK8#) z%(M4C`kXj*<%k@}F`0|EiXc}6;KGCqjvezVQ}~41cd1#G3g^iZ7z|#UwXv}2)w#32 z$GCiB>uE)Y3$skXgDlDx{hXOrTD8%LFTBZSrRfi9NwN!-PXb-VC zr{028t!zc|PUBV_k%qjse+m+d}GY#xe!3dn0KM#!edn8%Y)>)Z7k&u2}!lA)<>zg3aut=d=HG7f#Kk|;N61OrX?NwoR?L{aWV2WSCFRwWtXLdWVBi&^C0_CPN4@N<_=0VGe)0&d z{wso@R4B)ql0h35u)VuBGAn$g}-$q#I@6t(g(=31h#E^EuMNYsy>;+s+Fh!wjH51DH??J$h6H98~3LefY1PK*f zN-KKYR?rp||4EQ8qEzyOZgbNhoDis5**=s*J&NOGFiWJ1W+tXIb35AYu*X?_<}u1 z7vJ5S)JbgTGnccWF0^UH$j`(%b4)R_i223OE{4d`J<%g zNOOFjd_o1m%UDF~tBj=F1|}6}67qHr*68Y)V6^%bgzH5dby|E-$?F}fsbV%a+%M;H zuhSgow~mB7zBBdz8i{e)cBqQxDiJf+3DyBF;70kQPU33Z6e8VXvU}P2aw5S~J@heL zVSPi1vT^BNkVJb+QDaBTwk;`Uh98NtX`D!`ccD43w8Y;#&ZS#|r-lg2G`G1-Xy0cw z4n~m*zAVh5ZV1m93YouNVb>Hb8JSEWIbeE8COn(Tok^K&c_D>^yNdWjzSOyWA-lBL zn^X#13<6Z&X|dW{Z{khmHUrWtQf1z(;B zp=RQ-i^KX%DdxVQb(5R_Cc{gWF+c&D7AV;8Y%iS8vtc}v`#b)nXaFvV3Y<{#& zEqx^MyZzl_$RApG-4XKJZ*E>qVTOl@5<^T&*TtIT-^F^3Q_D!6~&+D>Bhr zR`=~Z^%_#wD2(t!=h1AS!{5{?eC3v&Q%<^Y_GAV4u+$T`Hq{H45qco(2F_q%ONyRq(_=qV(2Pb)tbMPy z_QHWTLN&xDXgIyADFwbQ-wX&g^U=tgdQz=+O%$YOlthCAAIkr`Z+9*P(5Um_{kyLR zEe2Lo-lY-^QQu#)L^DMUrtMQp5w1Q+3bo+hVJxyq)w`g05s>tlrlTlY744pK17Iy` zCJERDz-UY>_2`POuBI^weL-2YD_GE496X9Od7W|MMV=IS6;e4RPM;RXcq{}ey<^-1O=gZq2@uDZETES^x1XgVQ}`C+Dh2M+QZoyDw@~^TRu+Vs*QU)GK`dQdQC2-XSf`XeZ=WQGh2=BUBH+lpFMKpY% z*|+U?O2n@1)8yxPTJiG#dg?Z+^o^)U9hg5!dt)&Yst`q6#QWg62-HU3Y5A=Jjz=$+ z;x+UY*?Vb8KGnmn;ci5wO=mj^DPH9#{1qqSK+!HDvEYJ$DpXV9xKoRmvxuQ)$42$W zD2;R}-AKpYY&t}`yHn|qZV(A+DQS>y5RfiuM3ImZ5s;P! zkrDw#xU)9~pMKx@&%NjT=ef^*X3wmC=Us2FH8bzl(FIjj1HXb-n zUpW;_<6S?kTN|%yZ1aB5Df2cXoh9TgFGs=+LwI*jf(1jj+UyrKU+*QQgUd3fnPeM}J@KYv5{rdrBSxhFf- z#6$AR;u*fwI0+oxeKj_qjN=1!gpVp(YlR&Zgff^l6Cd2BObxD8A6WZ{2D?ziXcN7R#J5F@NMi}eOEvmC7-S)=il+P-2G#zzK zFEj8JwcQC!(2AxMtHg0D(enE@b9$r5{WGZaxdITd%V#rd<#j`^%Vm%YKUY*}r-!dD zA>~$5Q_N(Ry|dxAy7A=YfwS8ULgwX9IO%V~w-h8v$i^&06Kv&l9YPkHwf2rS$Ud5b z(<^nJ+~D}K#kj`R-1xyxWJ%!yuO$H%}sn7w@! z3HfpooQhr)Htg)JQ3?~j3_r_Ty^G3!Dp*mGHN#nAS82Csl1E6lR+psDVC5ZZ{-I`O zSt{km^lEv3qV-p#!PLlM>XX)YJYT*FMt?es!{WbB_Gx7X&CfQGMabw_$8FBit` z|JLtqI9`bow()D7il0-L^K)JQc?bX~4a5i<6gB3~)YWzt*5 zd_hA$cE)!Czh_3dpL{#HbaO@1F(oA#L2SWfUZhQPVt@Q~*gVB^1{{R*0e`=f5AWNC zPZ#^X0joIpXZ|jGvmpi8VW%%shOuw9get?6HZbK08~VG)-AYEiHpmbp;)TcO-sI=8 z*Yax>=ji24qX?zclDedGzi(xctc6$Gyb8AoGtk1e$h?jY=-kHC(RX54+7|PjWjIBa zWMu~4-Rzi|_jU3v%y^dpe(A_j`MfA9EQ!()qy@YW?d;fP=WD8M-2LIICAd*Y!-SSkNrAdNv$fAWf(02wiKI8!$~Xe*4*CHbO| zNKSQs;|N=Wfya-4gW&x5=7YF2N)L@WJsqYj&8Ps4)EKJ(5w)hz;G`=aF3DgBuIk%3 z5CIzPF4MfcfDZ0E&;niqy%=bvR&}n9tfSD?d9I`&I(3M!dyYan{B{M*H_B&_28%NDI4yW+c6+{ltDG z{m8wWDo-lXdh0wS%QQghyp`9lhVM+n{$3NGpCOrl8=s#(+59}&d}VYgukxn{7{<+? zSLj7G_*^Z>%WpNaoeA<5OPj`Il+Re^FSh!n=u<4<&3UJYR!}r48^6ChCW{yqO1g^V z)y1K1q41(h!K}1|kLXh?eDNz|Ne@lwF&(;^GD>|)D^mxF1@;bRPMMqeA*eoGCp}Vc zHdl;4RyQWU)ko$Iw4y7-k=yCDZcQce>ueEezT4Es5OBjrcRgKu!uqW^Mj@H^N^tul z{MGc+Jc_G4y2Lx;`N8^HPs;*bv{}>kEoRyuogKF{2IZBNzCUpwI_WKo&JNKnP4=$x zd=!Ye>yoj{j_=XagN$Dk#BiNW0x_6k=St^`ySW1;{yq<)f_MFRStfB1({kLBN2XX7 zqaT?n(&Q`R8dwf*I$DQ`BYa8wi2LdR*Zv2JkGkjQwqd2yCR>tlmHfNFo4AJEXosX- z4F>RYe8slzj@iun7vZ(xEqul*R{?*ZuWCrSu0^#(XQ#nhG0UHD(}9Px{fhWx$-vh) zn7qMTUZHK2#d?~b2wH^itJ-**m{JEeRZfh%`vqKs#l9esC^P9Tq_Kee)P!1WEhzvYaZa0ThfP2Bz0 z`1B^LNXp}M^TXzoBsB5OKJQ9IXz80CzWPck{E&M;;huamyzwh7I)*W7g_<1Hp^jUR z@$H-!uOCZDAJTD1n|?{Eu}p|Hw=``C6d;oRNI-)j+eSigw+ac7c_Lr`!xrn_6?Z4^ zb&I&^8FHWIQc*9KT7D6@2(8HnPv!1S&f+{BW}WZ@Uf;vsvgly2Dt}LSR?^=b95%kl zUWi-wc(fCH*L0t_WduWp@#w_ZoZKDXozhO|AQ{Eu(FuGN(Tvp=Dy}*eIs9F6Vc;+X zj(hN{@Sh&8Ci21M(LDK>C>hOk9o@o+l9pIt)=f^Ns%mJeDAP6=TL%7V!va0g1TEcN zm#dB|*4MCIrhP4Gi6hU85F^?NIxeLh8q=ZQ|M7=ciqKjG8oiO z`ntPST0geH5+S6ipIeI!KW6=pz?f;kK_Ja=y8iC!*p^lN)0HZ3UhF>{a<>V3ci0ve&y1 z2K_fLbhbv7o=uo1K5h2O`jC#!I{wI0?3u0?&E~C}#+b4`*w~Pu$SMWeb8lqryKkhs z3hNtReo*gcB55J96(R?dBR!bRb}Sy`W4^s2L9W;s>WXz1EHY}WelPL~Nd)Dq4cw(= z0%h@Lt;SCA)+=jFaq^D^xN2NiMLY*m=x0s@cShb>inWNCFx*!e98T&sI6=a#L)>jm`b`!pCCoLaEuK?BGX4hPu>H?-?9+4t#3C9m}8? ztX`oJ9bOe?kmDW}eelB7c$n?OgivnIB+`VI4!*Y${4m{)ZBp``MfPL(L~ZJ(_>-hO z?bo_Xu~Q49dW6kX; z6=8yh+Y-+ks3vG=f{P}6UU#fwt;?!ehCY_djL)-c$cx#Q366R5;$=uvOqAFl+wAM` zy&&HYB&nZ_-`$JG^5vEa?3ko#{I-Wvyg=phn$l&H+%`xL6jYknP%2(~oYZ3uj`Uot z6j3%at1zVz*0j+q|FOHsf^hY7EhbI+SKHU) zc%56eyz{=;d|o4VTq=P=I~Xk)dAA8glQ-?QC{~t{Q!9H+2NR+<);o<~$YmslrRc7* z(ro9=`X=aobbYW9H~h^*!jls&8m&WvhBr}VLkJuiz+#$X;@~p6GOEiyBK>*?WmWO% zz)Kqnj(1)8jJAyx1Wc9n5XHG@ySeveOX5%CGmD>4=512emvb-P>K&`|RGgV%KqNa| z^*H}rVG_|Q@f4}WhHoz;$Ees2-@eCK_I|6)onjTqc8L4a-NahY>I7Q^yZ$GE4n}nq zvEevvf*3I?9dqsoq537beU?d<1(BR$%ZmGNLsluxopUj?`Xajb8*bc6&wr=I@j5?V z6c=0cY4WxNpTrs&TU>Zmg9u^RlZ9z#ra>{+y_3ioWc1@u#Oe141xc&Mz%y@Mk^=|d-*%yV}g9LWnef>$$1 zq|z#5F6YPNB`Oh^uZ{hZ4=Y}$uSPQRB_MlTYcxyOdyF>MF6&xh zcxG=^UXJ@(puyuC5z-9&Z zR`VbRpRfXywl$&%2d&y&FUDSURf1cD#Cv5!6X^QG^v6gAMP|v#8Ip_lb4}(uiasW$ zf&x*Sbwh4BWFA}{gjZ*@oD&wwk3AW2?E5H45J(>i+Hzfw$sMl*n*(iiu%(^1Hz#Nr)|2193)|5)F0_GT8=@c1J?^V^&W--1oL%q(C zc;=d2suZ4kU6hqaFNMy$ERqmImG3Y#ioBz8JbF3tv{KFzqBYPq5dnR^kSzt zV)FL0JmcvhkzUAR^3VGcJ;ch5;&G)DOWQrrpyY^(K=ue8yLf+P zqhuhq+KSVC)xO*?3;0U6ih87X<12#IyxrR~y6ozLQ(Z?-CDC;7$r5_O-!_6sI-NC4 zvXkwVQXZ5Jo%W|?CmpR?Hla2{Qk8@vMJJ?phMBhV5sW)gau^Wi9>s4^AmHK;$@Nh3 zOO4FjeWcW^WYoG`&0Rs-PYiK6_aS-)`rPy2Q6IG_RYm>Y%n_d#W668k5!2V>_Gn~f z65<=-U+Y++95XCttTFDWO?T|`P!&rt?`XGw4xp-CDQa66xu;y_v=DJ-@Lr|aq}zQ6 zH~0bYl2_i>(sYn=Y!QxpHv>OB$6g-2$5)Aaz*247HS#KYpED|_MyZ#lB=NJio{kG8 zd9udv;OW=ziK7H=MP)cMVJvrf^~{xL>L_f)-H!`ju&)w*4Fcx&Vu91iue&z2#Cl+|4V zZ2T|DUnfsT3QnCyn)0IDj1g+c>_hEHQ0KQYJe@*xaFwYyD9YRN5{`@u*XkOy?%(so z5$Q6Rxh9_{i!4_zi;U{Rk@cuO+ba)LhUT03*|bHlq7stAgIDwRFYdRuJoTN(eleLX#*nNKpw=e4`25B=t*9?;?Cp9~cjywWN;?uiy zL{*Y&LX@;zHba;Z*GKaG^IGMC#;^2#VtSaySu$7nmCU-(B)3$cBOP7LKm&|wq1f8Q z>17nq`Q?mZbm)vJC_%JXGmA*L=dnz5DRb&cgUYKr%aC^1{FbeXZpm6vfIhw1v7<{Vbo#`~2h!=d+*qseBU zr-m|n8dbY2H(ejh|5|VEsz>H#!S05V`}JG<1~2l!H|CGi1raY}6A)Me?`>yyFTA^t z!&ln&?X^b*_A(boIU_5%u6Qk}_@Z>lq}H_NMx*mra?z-^g?hYt1`Dq|BsDm44|`-aS{Fu!EN3fu}y|>nP4#RO(t>shaGvoi)&KvZ$+Uu1(j~Y6Ha^rO?}wUbARcH-VjCKg2X)**?d1_ z<8GB%ci3TD^Of1ax396;jvqe1L#gCqAKJDXN$b2BsE6pUCNqpD4LHxa^%8>-3ciZQSLjBf_4@HZ(jx++n$qhgCMlG%jSlwOWYhi>;|~ z%;(68=Q4siLLc}Rb6vP}g!)*jC`kv;qA)+Wu*FU<-_Co%xM@Nl{oT2&J=_$!!G<|`>ayuR~sts7ueauXQ>&e=ot&V^80FY_}C;v(YD zklveap~2D#jG9+HVy#~Fa9FbihM&~lm_syCH+e9#eJGh9J3{C7P#24hxi}i`|2KH!(ENI z9=*XCxYDw}lU%8hCYsra($kEk;z+MNlVTNvp0Y5(wy7)aW`9kNQ2`~xYfYboQly9T zC}Dl$&QbL)? zmtXSHSfte2YK+ws$5bkkUFQ{<8z9sWWu#baUprRXt zjh9@ddf#6E@tu^r`nvs|a*|31RnF_i_}dN7Fjj*X2ciZDdl{Fn_91_{ioBv9G@G{` z*hNhA(uG7hd(;vFqDyGHHtkRSSmA}-jZF<*_i<#kp0b_NFPrAbI(+f^@WG6eCObno zwX82I)m({hFw|NhB%xFfMo&^21tQl6h#m;5PxDT&+)0am_^n4iv__^_nOcb!98N-6 zSi80Em5}Z)aA9L4-b%VNQkAl~e4u=KASb1D^k{E6r*UpkWg(|B!~N|I zwq~#u#m33v;R{C)*wO2F@foXg=fd9V;o?+CjoNsJNSjyF=@Iapn{KNBmVb>uE?fr{ zJOXlNH{saLE7V+))ihSI8T_3Ohc;r{o7hc^&_Pbf_islQB)o8>mtIHCaF^IIppfrt zZxkwO0xXmKiY3;5dD|?nafVYf>5*Ov$%PWz|E9@R{#hPEp-FKegpGG~I%`wQ%%C zfd1GRO`N8p>fig*@TR=)i7cx`Cdmc`U@Vr^fs(Mfvce+WeBHAzRu)nr-{jz;h#(!*93# zTy*&H+z9aHPW~Uyjqv`+6cz?4S~ys_SyS`!{yT3-d0{YBLs^g4exCilBX70f9`L^A zkJKh;xM3^7^s$jaSUtvHQ*m=JcamCfH5-sVrjKTLNgSVNNQp!0UN%wLa_H|vSSE~d z3Oe`Q+1pq-`gYRg$>{8dcqYv6>l0tMka6C8bbPvbjSYR*G@7>mU_7PxQ?+oJwKmD8 zsmP|%gYl)Y;(=kJ6#XqNo5tWT#|I(0ArIsYPL9eNl~d%IzwD|{^Y9Z4=^QOi;}|Dt z(}V)=P*H%V+FUr@56_B5H655}5?`nOpzxMQ+j+ZX&Y>k8C z(muXO>G(<0X>@x6Z*vV*&17ur?bP?IXL0R`49BY{G+*N@{R4)zt!>Q@&Xi4i6zWe& zZlhKPf|HC|WE%vPvnaMXq9<04o$_?nPX%>Xl8e32*^@8_40`jfj0Pg~yWpm5@XgTF zu}^Hf3}*S3-6<67%9;06*6COxJ-eebasIqfU2w$ANK|sLRLo?jTY*0Cjg)xRzWt17 z54Bx+CHus|xYS6Q>n(OZ0tyj^MFpaAd*ZtgC%k07g-mY!$BNi*Pz`6KRh)HH4EH=A z4Q-kvFl2^5qZxHoITDDEg7SRbpRirGv1jtRpZyJ#;s>1sB(CqwjIUEHL|_|d^nHwNTBTtUhzG)- zp<*=;Nb?#8H9iKLvADJ;_J+pDrS9a-B#oC!iZa%sB_U0&N7CN6%iK}uuK<=>i0KdR z44LhG&<`gj+iEe9+xg0w;d$T9czFG4_aff@jan3ru4YG%!(Kz02SnTTgFR-NVI_N( zTTTx`zzPF;3#CPH9VK#CSBJm&$02w*tnRfjQI1`EdOg@{XL~jN?mmSN1HqfOV!Kgm z1MDqJFZLSJe4?o9^ZZVqOiUb9VxcbR@AHab5;xf#?S#5 zwjja+4#O#4Gma?2bs1H-`JGcj#(!<{IbFP#EGUq_ zr3j73+$hGclE77B8_HgdLEN?P!EXxXx%<^P7`Gjb}xCDzBAy&Mur{u#Q;|uH%B2?91WW>&f90iZ{we zizePTG2@hdDXQz~6UPylx=n-sQU%`Xo+{@j6UTz?X>Zg|trY?YiS;|_bV|n(@1{ zCR6gybg-QL%6$to3<6VTRtVxpAr;7~=ei`O-n4rkTc@^V)I`=hqt(hZ{0>l96(0m# zGg(lRA6H}~NlMZ_SR*Ojz*8pFjHkEC8NEq55DZ`R((v)q&Vl#vu}G+e4|d?>9K1Wv zt_EOsvMjjB5*xmZO*WoSJ$bEN&ndzex!_$4QDNtpk<8B`1F7DWD*zWYRjDYrpPZs&=fEbNI;N!}X@3NEAeww+%{^%9$?%a^)nZ7G?T#b6IH_ zWJ|gW)2{>DedCN-91#4hOx)@M$r2_@4ZYpgNd+EQ_VLPBNPiiW!Px8Z6x_+92z-Yd z-Y~jGiuVco5doI|bc09J?b%boqy#uPuS^|&q0b*Ar+MlsZ3RS=(}WP)2<&C(84<>v z>ZG%(Xs=GVr>L_npB2Z?g>teT;TA>TOc*eNxe zHxvf#vV3qoaUD|2V+1h~9{VTu9aAFW`KT^O4VL$3>Pf^lE>yT$VBQbU+-YEBj40g@ zt?!p#zLUcA(x9uBK>=*YAY8Q@n{>wwBlY$6_-Q=cOGZ-?^O? z+(T)q$-v9;95Pi9*}HK}kv2>Hptg6<1NrRWi!+tOaA-LyqM(Sp#fsVHR70Q72fh}R z;)&h*G3J4T;y28-LeVuQ2&d+Wjgc~$5SFLr8k$2Re2(+czA1t|x$l`9(A*xDN418& z8q;7*mvxpH<;u&tlatytm0o=uq?%5U_NGBueLM!|gD@Y#j;8osCT@8&?=5LuY+7** z<)PIDDgNd?%roKSJ@q;+jsqJxYd@wR#c$B0 zqrtpq>QFcl=Z(&a&nL7nerM2IACjE+ov0bew91ZSon*6Mf%P_Af86QMRtbjX87T&X zG1j(xdD<)7?<*u;#Xo&txTZ)rGeK2kNRC9C<(z!CqfbfhPVQPzum`{C+)_+R^b<+h zvaa2^LnRQTMRXH`zI?TdoMY)y2SnItkrYBN3 zvkWV9HHyyrR)sa2%PP0K(1f_>cs*$ms8E9SZ~`-ZCiNt zy55GXTRaDGy=fFqQ{u2TB#Oll ztqPiCA5${TIrCVTJJ02INx)f3pX*M)zTey;U~@qja)5uCbbh$8NZnIZG(Y`HBNgds zG}!4&opzQ1F+!;4^_5qnEBPM^J3NY6_RMgEw4MZoRtY+`6!a?8k|>gwezj8|&}mKh zXg*$s{nale&dPpbTJkJ+JOtP_@k569uO6RkJMwXo;;z4{{Ux^ma zRqs?sn+FCWY@X8mp{YE1o;@6i-?H47#XdMtSF#7pYkYs@9Ei#ZMNx&$x{{?mlZ)o?>Xa-;UbI|VA|Z}S$TPA95G`Ko;4jvirQxn&(nvY(BWu)rs$(sm+*=s0H9Eg?EAFS8-wwL|prhfA zJUuVYZC6I_)NbI(lfI%3@KYR-8AQd7)@Grj8upeds%>*P*5)sOwPMG(A8B8GUpCCwv zp>e)kG4C}PO?=8;@wD(!ul7d8VXN%l@DSf!uYcoOeghchzSv2>} zZidGMATg86uH-Xrx1Y9SS+z*-*sI6zl|}Bc=`hljG-ASwiHHjNzTPDGgu@dRSl-+5 zv=p`bAsYhe{Y>9ZTJ7e#^?Guuj90b=pEhUMPWY+ji}5wRS%Ys^?5I&Y?LL|37UOkI zcp@Ql`U0(rT@%xIi?vT1E%<|NoNM*Yk@7jb;Dj^RA5JaIS%vPThV0Exg=-c2Zt$G8(X^5a@V76KO}FIRv(44_!nA z9V7vyfgW9{e`W#TLt%wLJ#6R&(D(lWKrYTp0CGWKiN62{;ro&I9YQWB7XFtI!eS2J zQONZZnqMg7y1+32^$!U7|F;lwU&bN#-{6o7`UVbl=)`}Ako$j(L+;BsQSVZqKguhS-$Syz}Qp4n-D>X>L#@v;>7rv839F-XU!-7@ zp^gr_V|8>RYG6!vK;rRVqIu{9UlJQY6%#iTJ4dS@WB$bGcOHS>%%%Fjgep2nMOBHK zQ$*x@MVO5JUQ?Y99fpKotOgUhKT07GbO>x3LFy1NI^^OlIO?4I=$x?I0I75GqH|uv zIMq41(K%uF59C770i6J#6SxGUbAE5ej}B|}Ybcm3{V~)ZfWU<7k5Y&LIs``0|91en zq(uJ22uz0lZvcTw#veodff1Of{89Sfngp1r`_bv2LtdgP_5Zun0F!9{2nbBX{V4sv z3kXbb{TS*GjKCz>kJA4dQ-IWY!RWl-8HI-*o##7$@$jPa!1(JI*TAIK4*>oI58CiY z>35F%`_6y$2;;2pZNyDnEuh91NJ>vpOhblMLP<)O6ATlPk`}IJE;deXjxN+(7m&gX zE@*AQ+OTnPb(63*aiJCfgA`4^Uvu(vgLG`n-K<>=sCi&F0&NKV(7&HAQh)yc|EQdg ziH1hbdu@hP$a7)Ha5q0}!)6j22kwceKJV^`jjKisYXO z9UBKR2UnY)X-OMPOAEk!2JB@6bO4z)faPoDVq#}*(einaE|Eyqi2}G!r{tlwFjh&l?3rN}yFySRFfIGDKv5Ez( zafyhyy17`G*rR)87(Eoi?xx0-^il4@0vo)Dms;ZHz?YKa2r#K;!o0FX+lWIYBo(1& zGx=ydD0a?qq&En!OH4fIA{HCj=(tw*%r2&2IS)M`sP}C+OF(<&A(t-4wSt6z z(i_rlk6m;so|#dcwG-ZBVo#H5&s)&}Lnq-o$UoHHrF;EJ<=S#G3X;+Qm22qn|4!r( zn1%&uI%+xCK+!?Xf1%U0b!yV#t3Ys9)yCg|G9I zX8wJGg2HKUc8MlnTK;!wpsA!JE2#?LN&QbWT`j{}KahgX1E~HV2NX5X>%T$Ke?bad^vul%ow|$viyVjm)RP2LL(lyDJk&e_FeqRn zasze@H>Ut_aKZ9n5c2*8;eP=(Fm&Xm0v*muD za3xh`6=@}4!jqOc{uriGXOQL9`A2*0u~1ON5umz z{SOrnFpW+YE;f$lATvjMdlPg(Y`EH)xLN~ZgBxUG=VWaHC?Y%I?mVt!?TOtH3KJ zIO2!pc*th07T>Ng!8vIq@Bsf76l%s4=`JJklMTDQT-i^`~U6c5>kmL|a2;%VSPo)GSyx)0GTT($%M;FF>e@I(^ zDqyVlODVz3!++pKflJ~Sn71FM*bm|NOXU4Mm2-kG80jBK^=BhkIXCPH+lxFN@CDml zqy%6e6aoWM|3c<+3;ZTa|4VEIjP@@f<>iM4DSt>!YW@pU{*$!5)U@+Lbs&?A;^KQ_R6>n~^!6s4XiF-jJV z84ya#PK}7J>ktQgRW^wNBIQt4*9&L2sqvOhqcoK9hD(bIv{WJ9||Hsw8q_+PPHC~j%)ZmZ){yf7*fU3*y3iHzWR}Pq(6AT9b zg|sge^god{wDdpJ4^Tv!I9OfsH7+G#t_H{yO7NG=O@5H6yPchd8=yGAAPZ-A8xIpZ zsG9?%xj=RnuC73yi5*O>xZ2p;*g=~EGXyQMcQ>=X#4@np((h7JLsdpmLK7eoP0ha{ zlJB-W)J^-{Vdc6kB_ThJRKR8c%);*nbONB)-~Z8niD`UhHrmHvmc227B% zyCdM7LS0;R0p5!g)W5w<6HqXIrwMT#F>M`nfF{KMCWyfa)k6PZ#r}mRE_riYFw^&9 zo}t%V{9M$qJm@_00rc|QaghtH|Dz?p0L*#?>htpo01*U!YN*Z25BQdRK!Ah~7Mv9T z7@ivlIY8U-3P1q%2L7O#{Jh-gd@zfOABaQo0_!CF&|JWOfp&HQ)}?VL5H59J`>;z0g-xcU=;+K1FZ=K7Qa1@KaSO$Vy6zLv1%8^8I^2 z;X)|?tzP1~td{tHMGr1SKQP6>VggiEK-JR^Rdt~NU^zf3^bFKShqm~44rqPo@l%Ze zlfwtROoa~&_-H`TpAT4@fcjv(z)}StKw!`^F2D)q0i0tfB)ovE5df@FSQ|cmbZA|m zH8fBLqhlyNU-SrAkwBN<6TASudxuEEQl|e@ZZsdF1?xJsC9H_1YAO(K^Avs>Ohl28VL=PJp zkPZMS=umKe)%h_yuzC6mL0`m|{)wQWrT-~Gb6?0#b5pwuVfveBh5RgrN>^YJhYJKm z&j7a%WML0I0s8!V`~Wazxj-uR7FH%8OG|()K$cL=4FQ>$xw~0_%uHMy9e^|s$lk;Z zxV$8%VL|TS<+PN(oR}Ugyd>~D)4<{`|0Ja^a}DH|F$bm5U--U!1pz1-;E_w^m#=@- zMyi(VetGD zEkR5F+q5L`SL(Sm=djg`-{)LUMO{$}8hTad{e8}1VV-}ObM8yl8t*Sb36(1s=>FS( z0QLL=mqV?9;s4jL`#)d`#<{=Qrnn`Iat{NJAu$a z8g3T$+SHufoai7W8&_9g*#s7apyuM@0y+W~jDX!hOrT4C1Z+v;_ce>(S3&;e3J5do zCIMgC0e(ZDln~h+oc6cniFHjsE%Z zk}zgq9>z(R%e7$a?wtEx@EfKSv(7`SmM@nq3fa3yP8!&yS(-bGYH^I-yx}gms-rF@ zYc7T&Pwk5*mbWYLnd8GN4yx1LR^g7Llbfer;D{J;;wAnd&NA?!altJ8&T|MEJ!9FG zmM^QtCj6l#T9OpLj!zs=y&5Gn@VV_IFnbu4GnvbtjuLWDvQp%|Fl2CU!SZR+Y&!#| zsL_NcGukak?9R(U{CTbzbk44CJJlf!BgqaS0k0F<*>p8f1u>DV;o; zWD!j=CgjX&8>hDAH~5ClhgFem>b@an(wSi8!gCqr-W7C?+I&6cJvaXT4h_5FlKtz^ z4F)44ecj*IR4cdzY*>aH4G31Tcjco;~!8&)EBf~XPx#f5T9%Qc%q7|Mf)N11SF5FDuvP5QSf*F@HC zcekfFU6Y=-z1h{$2o^%7LhXWFv7)+PwDf2nGw$nk-8Fj-PooM3S`Et9*JETKHnUlX z`)_8+aJI!(+CP<^a$=4CQrIGSG`&9RlJ>q=SjII5w}Y~u5GA9a%jo0obBdK#?> z@+dsqPped(Xy0KK?|?79VXE_I&;_>A_3y}xKbFgWhzuYc09cH?fJ?yrr^o;@FUt&A z@~6-MTK<#F0JPpO3Gs8O>SFooGLHk%azNt#;&E8||22PRnvtFHNMHSN@j@z) zKzglKj+>~Tp|6>ohR`r2aI-O|?s~m!7vD|#dM1!SEt3Ev4}C}Q{rmTA7Vn~jT0|LP zzKV)c)K@q;<1#4Zya$@^J6&t<}dEvJD@KxQBdzz+7X09k=-Kz1N| zkORmOKxf2d&T^m9+i0UdBOv0BB!N*bl~;HiBy)q~eatnQbp8pjlW<=O|aS2!sz=pS6|4>{E` z<@C>Al4kMGXVSW^Li9pN9JTMd(Ce$T0yh;2LGcPVuRgdYj@d0-=nRPyrqPN^WFWMZ zCj=cDya_h(?-Y5?Yuu&F-#{AKNn5RIM?e^gY^pYS_Z5K;>OrUFg4E=+oeA;%$5oYk zxHFr}i*DKZOEYet-L-|X-A}@Wl6)d`#iL@?3Lf0sCx3X@Rdmx+zk?L9Y<4yZPbk8V9PMY{QWz>LqbXo1<5jF_Dn_cvBLlo{T#9w8(G$B zZRD5U)>xdMU!EzX+-Z1{PqN})lW|7gvGQ=GK6fY#%f89}D!s`2$;|Cq=BX@kh?ny*)4?7#(+hz2FWm0>j)6k}$v!Mh|q%nQz_=C^h%Z4M~=M`s|!6SMS1;6i3;mojn zPOCY^6XT25UqE^!HnyH@^~$bi){0|kNcBqkaM6ypOY37zjkts|t)20G{CmlF5^=E* z@(XlcW9Z~0vGJ;{g!?x>MI5=xWMN6fN-jn^LqTnFga17y%@7mn6ElfU9s4UjhBb{~ z5Y^R!CV}Y@yK?*LLS0jFvxeZ}F)1@c%b||@AHzSgad#s#-4`Q_#iUg=@IKahNo<0l zf5WQ{n!LefLcKtz}Do%)zhy1N|*4GU#dd?+Lnn^H{5pc+4<&pQWt2 zqe)dm6;7|9MzHuy)W}=INZIFAg$hOpg8nV4+l|*`D?1nF?Rn@)xx8%lPp`^`=ms@8 zPpp2J^umo5n-9x>)Qc#PNlGHSpUJO9F13HD>NRcNVBkZS z$UB1xb0P^6MA}3fTB4R>ehd$t8EdF3@7{KHmbmd`LPvY+@THlGwUQ|9a8@qpOTlf3 z2ATH>{o=I{YL4kU!TIL4W=EdAdL4$^;zVmHHvIwRAvGv66^9DsZkd)+XDH+XSowPP zYTXm2r$SGl(|=U+yyPnp=a3kdX17$8 z00U2*s(-N6%kUK1E|2{uS=o{ZvyG_gH}P0wb(_6s`wK%}D!*PIMq1j+nE_i99uQN} zULnkz7ZJRdZ>U1lF)THAjYFIEdPjM7SWo!G#@XZnWPJC5YQ&!3;AGz>A3+%+B4%eC z-}U7e;RdlYcHiEbclmsz-q=+od-?I9vYN_0l`kii@_KZfEK12JK~KFsjnc9{^7F&l z4Gd>S-X0<#pt9Ar6^ZA-+x|iv-&g2A@9n$pH8Co!4}v_n8&KY4@C9q9oIfypHZ!T6 z;x6k(fJWl*ex`G}nVs)#!}ieAPj(%dry{Ob0}Pzq&{k9l4HSIJ#_!ZV5&>BrL`_&v zEph~4QNNArlMHImL>7qE$WgbbG~jW4K{pQ1n;Gw!Yj=B4EX+59f!sw5Ej!@3i9c3E zq;jPr(L2-Wka#4yKz=ojX-cvj3xD_Ze%%}F_Z6IaM`JVMXCNhEUrG@LzjbN5mOXxv zVJFAQNOqO2WjwZWNSU3MLjE2(scwtL$Y#?u6z?yY_T{Zkudh3PfP1{OqR3riyP)!5 zjpX!3m!J*dlh*1{TyzoYD??LaXB`SA1?LA{{VKF)1M0`~?6*yEJV)kxvUK;j%w+c{ zP^U(COB;%&^?erhPcW<&$?%mj{Uf!!8?X6AHj=ARfg8r%vSJXf%|hG`jc=CG`LZT| zV*v*^TNLf3lDrZsHEht^_uLkCz(9A#bLAm9Q7e7bb^;&!f<*#+WHIvTZFMozWMD~) z7RtSuh`MwkQXkCvyPbPKN%DOdBiuDM1 zSeSB5fk{lO<+e`qZ515=jAo3rPBy0Das?(lQDHk!MB)o``S;J$sjt5EU`BD5`kF&5 zt6=795LPEB@hHC0e$@V^B=tipf2oh*hM#XNk=FCD%6%&3+ek*!!I1=B|&Rw0^SHg(CPL-?6V=-jSYQ)z`khx^( z#hmQF=zjF}dd|_cwQI4Ww{+7%hdkJ$seNEU=NddKP4E-PI2K`FWhajsGrIPGo;x?q zw@s21d7g>sFbXr(Hh7^@kTo)o6Y0{XP05|TQdHK)o#tV=Ey($3y08WH<<(oU<2Ji^J@#u#Wov| zC$P2VFJ>YvHhjs)xfBil_sv8g=)!&3Ow>kmBCtw!pvmV(A@`^uwd39J2+c(DKtxMV zlf!$KiOA?hy-57bf_$JSS!D{6993`^<3l-Rsl-yw;hTS04=RP5f*s2+x`NX4JEE5-Z0+gN6p` z@h8$QGFq!-Ct4>-sQX$Pc;z6_lajZqI+m6xW&5lH`{j^M?24rw8zXy5e@U$@m~V6@r0qc1glLw(8XX zhqZSM(xhLry}P=)tS;NOZQHhO+qT(d+h&(-+qTtJ^;YlanLRyc@0l|v;{BHYj634~ zmYLUDzqPI>)SB2O-Z4Tfq664y*!qY>--!5Yx}q^bgRNho!6C&KhY?wmPYM^3u0;EL z5|Ma^xQNev6wl>lB@lAT`DXIeWAE2z&W8r$sHv}yre(8MYsrSz_%j!iN49V=D9MBd zdZkZ^Gt5Vj9smYw5R|<5OLL>BGP4EM6RP0({a+Y`us_=cw@eE0p|K(Eb&OF5a41ub zq^)IG$Ed%kuF|xS@rApFIc@Qnyc;XTU4Vaq!8q3AAc8k%H^~XhoQOr3QNu|lf%XL~ z*8m)}*#aWIx&RH`U8(V-8R#RVMQvqXa1nhEh7c>9Cv4{JTgZ|iZ@JhdI_ zzMI+n9bWH(Nx7~m*q)I{i0cDCrnmE}RsNz$?!vhasmjftcX%yGi8G`%U;P2I5hm=l{U?3OH~K&P|7#&%8H}#620) z=%!>1hhlx{1OIa5g<|mH;}CpxEmS>aZ8Jbp?;C8-4tOMcZNgU6e4F*ndMko=2-N4h zjeDQquKl&X>2}{d=@gw=JTMb#MR5u54b{&{Fg}(n%pNigDi$v5#1h zoyK?vkd`Ea`jZ}ls_^Wo65yG%y$0;|J69%~FQvY=b2SkObk}aJ`LOC^ut~E!N5~{* z67s*CwZFk<qt7|@o+1X~c6?$>p!0zs9UOIHOpUJ=lKh_>4^(&76jxFaz9X_oVe6{IW)5(jU+L(8X=zo}7ut-^y$mjZLsFOdvb6iBs3b+BTY|=@F-E(OL_0bs z%%T#l45iIJkj;$GRr~`m%O4+zJG_z{DpA^bm8m8QvE-apo+Da|GUfVpzU8mD#5zOF zXT5^o3vg4L)qaJ~%$~#fNnvk-xV$CzJ|sFiC2I25IRyjIZ=BRWeQj@wEoz-&o9KV; zPK~udnd?IpM650sY!Y>&&2U_z2^B3{J-D{NAo6v?=3GbR4fo+FKLT!B9zjtdI+g;H z2hbl2hNFxG2~mj2U?HY zyGOcg-;sO169&)ovuBaTc$uZ$&~y$_bvU8>MV;h)OV90({=1R|dCc z!VWGIE}mHAsvX!SYIm=(FaLTYIL0MF5nkPm8P&>sq0`S|L)9mckJkaK}eTN-C ztzdyBNDE9jvaNqQEx$c|-zJt8Gc909Bi<+1W^-9Rv=h$2ZACA50?`+aT4zBN|H3+_ z{&2HEM_QmqT(Qwy#)2Ies3M!@P-rb=Xzd-+NYC@&90!b{kQBYuajJ`k7-hFz#Weq|#CCrTm@d)@_!dzzHJbsglS#NgB6iN0B zy)3&+fG(!DjUu;OHC^9H!I(}0_hDBw9 zPALnK5Ub{CDjLEdv=7wQ)>Es1T19C)(MxPf*FFj@F; z*AY8|#R7hb@%~mM9DK4ntOVNy8XD!SHsIMXmtDbJi3Pl#iYTDG_X_%%I}ZfWlzaR>Ge+vegfo}WVAYTy6dl=g zZ2|G3ShVQz4Bzo}fqq$6t4qhZ|E5gNxhG*SPZYrflWVzDS)eWjWW{C+c6wn*`ONT%)Op;oRQ~TDKr=7D`u@!X z`hSw&|D6l*FP#2A22jdRgpZQ$(+{JiVfn96T;@xeE5x?|Cd?vA8g1Uv*aH%$p3zpq-CO~`D2!}Ofa{BUql>q zrc|fsOW`0JunJK!7&1>AtRIJzNrPLIkf7<&$(E|o|K1pLu9j4e9A_w+q&#zRa?~=T z27|wi49oO&RYzZp1Ql+CoK)1tC&l2UB$aM9!IXV=!*KzVjQA}k)oCKxoqNsR`|#Jt z;g*UmY|w6Vy;zBdeecY%>b9tciaY@aA#E5)=d*>bI#_~OcD;*5q{E4=!UxJqEx4rz z0(ti}v)RMBb>+;lWe}5Ondr60JVY3_^;)x*Y>6WJN}6H@c)`eA=|ZVQpm?xeQ)P>! zWo<`!@!65g9AyU$*~y$xLwBy#^TRD!hLSZ?2^4rirJ&`Jg5#08qEP&ZkM5FPftLOp z)|!JEeNaclc~7IBdo`vCVk6w+u1^Jtf)8bHXHT7A7fCDLbFGTDM{R!V-UPo!3o^Mo ziK<0E4)rhgXzx007o1&bkln1WwCO$;BVbH&38MEueQvLIM*5>iz40{igx`266+)!N zETyX1;%H*Hp`BJPjnng=`t;`RPpNK|8hcJ}R@-62U|@oEqS`@BiJ%2opq2%!r})jj zu8^eriI4yR)C;l5W5sxNgz!D^Sc8+v{k#|6?IDMlx$)7~XmiPv9+yj$(RK7qgl>r+ zuF=5}^djki!G{D9V7R9=W7hZAVKxP;#SP&nK!FknsSOq^yVTpMx;YYN^&Q*piEX?vWK)LVgr{vH`VZ!o#@9n-)Kqwg9YXH9sC`$8a@xDKCAPNx7Fk6UDuh;jI8>_ zyHbt}%5$E)M&>7wGkZfc*;&Bis zj7Ze?0|WOY*+3TzzD{axGL<5D7rR|o@Z}LX-srBv9&kRAP8u@L@h3NE8gm+~9X9#H zmi7kGBm1v05QNZj9@b(zoG-g-dSbuX!BwO^p9G%aN^N|%-qr~k!g%golrM5`o4N7x zlpaH5N$S&lS|*NV1v3*sH^^b&MidKvIueP25xT)1Kaj>69kk<-z+&vOyi0TYAtJ5( zs8rQi)VQ{IKJ6DxV_0H9lp$hsBg5cCN5}H3=;pFEpNdz9C9BKE>w)MyynoUFCz)+( zgRg@#Emk$0RB_fnJlte{zW_ulQ_WPPNHS$;ncM)2 zCfXou+0~$O;)wIr!Sl2`SbG}#xOVbAZgC8h#Slmo@ylCD|2AJ}SUoOVq%Z?HibRmD zr$u!==x*{HJh1#$5$x!usl?10Cj|sfhaB1Tp?t7RJsrso=gSYi{ML0iZ59G+Yu!9~ zTtunWa*?{Sh=QYnf{CMAFPI>ufE=V>ArdAXYM%9xf0LE9fz@Fbjq&n)>4}O%bx3Ul zPo?9gI?x;Aq=oE1Gxw8wVv}If!^YGAZ6o+PX)nrw zv4OPJ+_abDnk@+fp)`62&`%mq=BS>cU(oM6ajsUB-aS{LFjZ$hW8V4Mm)=LdK!{%m z=JesKqy6XbQ95Az_b%(SX+@&-&y1e>%D7}4Qbje@679NSFhMwQiGzB}IVxw}T_eQm zL}|q8J)K_40SD^<-YB447~ytRxyc&B!*rutJqNtkD4Uzqh%i^9QspV$0-d=oZ__jT zEjFg4L2n?3a;zW!<{kYz_VB+^iIx5(I`pSX{2Oxn|5YUZjk*03YWu(EWByLW{F#RN z&m_!0lP`Y}FP{m^e}rQG%bWZ+=lUl+@caD#Oun$t{~5k4h|xAtSHu{c6n!B+=1cbv zCv9Cjb=dj|b!cFZC(^=CQVC88S{vFnez9UQGx5uGW|@K6`CM{E+c}=-*b*L>d3wBq zl@ZhbxtjNz?6s^mFr56?jbsbDhLeqB#-_wMPT*kZ9ZIN`uvCOLx7(S5jalc_XLYR` z+a5N0w9VwTd-Gc`4E7GNpWe_VJW(?)E*%wcV8}hF4SqaV7fR~JAk??7Yc^$gV^HC` zv4Z837ma+y%uwd{J(ws%pQ)?h3lX3(KCQ*&@87;tEzw7xYYWH{j<=;MLY4 z&Kx+r=3LD_H+}?cyH$56k$PxgvSN5#`N83tNUF@UxErH0&dwRi(x-<=Qg_#QtSfdt z1Zn=hK|dRnpyC$X{zfvRDzCTBFk2kQ?XZ?Dy>zTDVEE2vW1j2_)7c!p>T*og&FtNA zYFUP3eDZ`Wb~etrCjoM12U4D~xQTNiq+TXsMXM+nFP|DM7p}Fh`o8qT20@4#>etf% z^QB^x|88ndQ5MYM3i@-l|)+S3i22h6EHH`F_g2cZ`thPRjV@3YH%rIg>O5tKet?&nM;5W5+=)E64(_8Ysxws;b7v3Jg%EYQ4IoZdkS}PQtvD#R4Ty4TV>SjZbKoDm z1E`=ie>Oy97)(n@POLi)tY7>!*UZ6QWWZA!q$Hz+@Df^cF>IJb?`o5f)R4IGIK{FE39crlX{kn;5UD_^KU8@@1sI@N1hngA z;a9Kwk;h-jgd&NK!r$#Buiszh8K%yoKf0Qb%z0%Mi%KV!jk7!}L+F@t1@K}W0{!D> z%Trgli9z)RdtQZlyhG$|o2PCDwEQIchWeVXJqbg`2Xlx#s4YE3;eziRlClG6&8ZJ(D?qo@S~pH zuY%e2oBI$y=fly27B0jRdhS(>I{VOr%Mc)GSgootI&8ChIDfvyRyX& zBFB>ZELvdnt}$+>lv9q^kKSMF>Aiq!nnu#;49++8sPx93Fa1%P2#gsoF)Kor*-CN{ zC8%@Ab3$f_hWwrqo&ugy6@^A}yCsryXl6)^zIOTPB*%Vu~UqE(y*3fEAh8l9rgWj1aTv7%M zT8GcKSA(4+RXBEjhP)PzN>E3n)70AUc}I((--oBJt&MP8cZ6Qj*6AR|`i7hv1tB%>gd~Dp31|2Ss4sSFOgk83dsCksAo}B(o7AK4k z49mr(YD6v+frOQlij@_sp9a9OIKo{tc`*jRHwIbzt?k>v?EEf)mhN6Ln!-d^{^*I! zu5ZFy4XcNp5@?tZq!#YCy3OrN!lgJO;T@ZXsazyXPaLEq2!>UWV{LKh+w7pVW=?%s;>(--N;PSSGZpBn9(z`L5O0>)Df{D zY_0{u#@S=(9u7VwktE5NkrAMF=F98L_K*}1buFEX{K+e#mJb?$QfLsj*gu6_qAi68 z+gkP=64FyD`TKv-t~JPW`>{IEy6VKIYv>S{=Z;6;ubFXsHmWrG_UgLld!C*&oED;2 zh8JpEu(m$YD`#EVoVN^{J@HbumiTU2)og*`HR66J6z(qyCbS)<-CnhXrT6Y{$dYFA zK^Q3QY+xrfH=PSIrd0tNP%RKG$&>GMbqsoEG{GDPUgfGeQxud{f`16AjsD>pC6#h7a_C zlJ4C*0Vg=LD!&4cP(M+y%P>`pq{$A`aM*975K~>F1?D^KD*=q@tp|-za5%6y6bnZi z)rwvsL%&&;@Dli--d+@LLvphlm@Q$t9W}Z{%1)$h`xBp-l(~ka1Qe5}|M|IJwdDh> zMbJ2PYy;;d<459*k>Upgrm1s6sRXDVupiaQ98a(ig+&l^c_|$|L?SIEE=hW?w8>_l zfVNjK14dwdJ)R<=GRSNG6Zi!r)Tju5tP+=Eu^+s|ajAwvs91`10b+4jzJdwLp>J8f zfA?EBuvw@l^_8*iT*PKIEmI~fSai=pmjoa_&C5P&itQ5Aq`oDl75DS=Na86o)#_yW zcJ)}&CH;1OP?=dw@p$JP50$SG!!O&gd%lt0HObccq;@h;Qdec%WGk_a%UdQj<8HWX}zVd1H(sRX;`@BV9qHn51P@kKWp|74f=;37rC9b0={0o$g~jj z==@a;lM|cLbkGIW5)Z2T_G(?Sm1mpB#hl{9rvc#Q64whK!M+wn#DEVhFzF3oPm}R> z3ZmUC^z?dlFbPwL?B0$zzY+8i?r(y*x6uW9St597#nbOtQ;mY5st9BT=yqP^!hpr( zh{x2AuO;yXi9(vAbiP(nKu^sD4!hRS-Vc(Je^g3v-I$KLwx5|gy)9@Wo65CYcr5#J z!aO?@rx(4`a9hL#(CH_1j{CK&hd=`|hK{BpC&qR)`-Jf!P7B$C!kOvamkfT*BeU&m zK;*Qri=5SB#+<$q6LwZpPCFj*U8X@+;J^u?F&wN5(`*g=PLeOHA0u#`rOK#$Fi}4< zfWw1zpl)?8I8yv7Tqw>tyTM0{sL$s2PNe)8%SLBKV5}~A zA6q2}HU@TIYj=+kXjvA*)PLFXQRk zMs&oAVs@9#mXjsFSb{4_i5DV9?6?U$16DHibIs#ZvOmYC*stO1>Znmh)!AZJvZnAs zN*Jg*HWpT9XTMbAOe6@Oa0x+kdgz|#vzMsL;CYBXG+X!?^b=+%n72|rlDy!c1lANJ;}+yQ zY#kGY%rOF4>%V4)H)+dmMyC-iFbBh5_JUPZzz0=3Bz*2?xNY65Za!r3K7yd9r?-_D zpUk|+uuOhR?G~00z()I}20U6%MR#pH^)oF_8sCtc*xcxMnU~4rK_9%&xVw-9E2GAd z=NxrUuQFVU|Hx|PSI!oiqbLIx#2tcK&hr89s`dHYNZY4T`drzOI%U@W70BQ!^f)t? z&f(V?VM3P7ct4gxaY7*5fiJ7sQyM55)Z9Iy8BAkgYt>b3G4!nPD=21oWRs<{P{g=g zwho@e$2GRI8Z{}6;=pD2SY*X%G~tzJH>70fZRkNex^r@KtWKIoz@Q#``-QUwQ1824 zeobo}46iO5%c+gM;#;&>EK)|D=+1D%a>q}cHsrND5P7|A4NpGqKGVCoGA(B7b&>|NVnh&NvBRZ&XlSLSXw3jSOzieRVpT30(`c6gZjHFlMKlb_w#yB+H5Tel zZ&ni(T=ybEu0pCpE!x~>HHO$0Wp?2|AewstPKrmrHI79FuP1jz^&jm!?M;VCsbpMmxMBt+Njy&O*YbhKlY_KdmWXF04EYF5bU50_l3zEvEblyLVK%*z zrc6$_RFX)YD(mN%*(ipCqI8I+pMyd|F?RCio7a#(053Kfe^3ep`7uK?#voXa&j9+B znODPt26Lw^OK>BS2}k5?__<$BsdMFnq!_}VHlbje9b~^us)j=xm85O`z{Gq-Uuo_J z4iR`aU{<;7e^1omKzUuh!wJ)IxO6}h#EZtc?7(`@hk+2NdCh% zy{&B{_pLSthh+jB{iWJHdZZ+_$J*@e7Hx_den#n6b*PZ}Tsd4RinFN~+0}1Ad7yYH zMpkBUjt#R zlyijM@h&Em3{4U$PT>q_ggG7j>g>}r*SC0Ti4?dfaB5415sR}WH**#sxYSZc zRB<8`_o};{GjRbqRXMHMb&nxL*Lj$OUiT6!udoYsNiqY3cU| zmlW{OP4Yw% z$&G!KWcf`T>1;D6+wPbKK@8l(kkWYLc!3wKy2ksY>q=?&#|$^4g=Y>3p727cU?5+R zN=)n#GBk@t4KHD?^UiDNX}kuV0m1ula#L?(#9DwEkQfdrQqYcQnC0NAJcdYEPI^RG zzO}|T^WICfaK~*k9v}}Q&%7qO!L}fRctwCxfmiann@-=9vtEMvTFtI**wka!8@2@D z19KTEbMI1EKh{|&LbU70&04U~VN#Ldr^Q@3XslBQ56vzam}FZP-%gvC<@3B?>I`l# z)VCp@;eb__Z&RtxZQCQ(TCTQVFTCqfQUYx zzBm%j>LGuTk?VvZ&hwWQ6OurI{nY^}-Wrz-_e-v{>n=H9M5`MJsk<68x?P%eQ^Sc1 zh4yBzy6>nbKb-*!NkUS5VFJ_TLL{kxMdA;O(0fEEzHZUU$p zb$AYAc%;)37UB=IETmYzN-Rw{)rDU-=rnSd&t`{Kf<9S;yw5pn?PH#%NeK}#Th6uZJDBJT+@g(7Jvrn@DYxv-HyCBxwg`|gvV$qv0@Tzw zLkIAHB#%SeJ799~u(Y;KEo{+KuLWy3-u<3`tLO#U5kq`e!1^$on;lDN> z{nPKDq5VV>KC6G2KfA*ICD#9W^G_tepZ}k|l+U)AKV!rhDqt!Ii%FG@1*Hr= zV(MSh6n4YRchKrud`)&bclh#rzRP}J22xP+z2I8MVGpSV3&W`&!b@nrq@^{Bm`2uy zOKvyGr5*@n6{JcXqm-$e7-PPj2U(YCI@xf4|M+#tt=#W$-czhVHsM$~dPAF+SjE5R zMfSW;bwh53lgX{Vf~4ZxfL7kN&~^5dgSl}Oc!@$zaMarHzAeSC z8O$2#dOIE}HOurASgF2?xyP``zP4%Y@67UsKV^D7z1qa>P3+TE0bF$yCN`UqQNBAVV*!!e8lMqf8H z!eY{m#)6ZR?$3{{22EPCnt6*c8<%&Tgc9E??`@Ccr;4(r&9RrgN&F}Q#nc!OUHR|1FG=eupN4`f^_Zc2^)hXTPJM6> zm>yPUS=FBR>b}Ed5r{KJEJIs)DBjO*yyo)|#OV7eBqF(jBq?$J;Uuz=E1mY$BR{Hw z9+OGAo0kRXb9rGO>HvQyDdnKS3>GgsIfd?IEUps%ado|Q4%o>y)!dA3FY65t^!FG- zMJ0l;0&(s86kN7&4d_Uk~rwSg=?C^hN%=N;D~PGqtm$;Mqx=*p2S% z1KA{-`*D*k=*Cp(8~8KN`j=#@ULJN8{m`mzuDIq^rdCPcFq;+~CNzYdT^6brhs5cp zqWR~~7huf9@QOn6Ee5nf_H4PSIFC4bJg0BI%i05?MoCE4&e$T;bt{WgP2rSzd8Ue1 zCcDp-om4Zo>QFN_j^YO#jR%Q{>?tn_otKnSnYK$lMFbigFtmzivT&hL_S z_oVq>#WoAc{Y{gKKZg)hA3mA6DEHS1>6UmJ1%HWUM}E8@L#H7a^ZVMSK*Q2FBflWy>h`oe+Ja4q4#zd z{u8-f7xhWTmF&=|TWXhzzhJ{40JPIH1y4R+Qlg#UgQA3LheUO=)N0fM9Up<$$sxD9 zno`umx_b@(@*+)C)ru%@h-OE!J=QB1^)L_dfv3{t2PqK~tX{;5Wj=LHD)Tp@?!NWW ziBqXT&ETqS1vENrkGpMQXS-$QX4AHJfOICM?m7<5O1GA$)<-X2FU;#1QeF zilbZ!0R^npMP$z_VfZVLloxGAGE)Q5*aB{`otrR&U@VQt#}k^eF>v}K%jxBjcPa=p z;ZaA{0)u%gZ<^yEl97b3;ev zBiP(t^WM?Ad>pS!)LamCpII}~PpL;ravR9;UQygZv6@$hma)pfXFP1TZvzBfV~IF zqxfdU4JmYh^-#UY(PtG$!#}u?mNrf)BdlCxvh3pFm)XmpoZJ+GGvFU)W5@ zMGBtlQngYX<~dIY6>0v17z#gvxdtv+VqR0Z*=$icm8c9rogLxQ@P@8pVxC1TwFFig z{AfUK=vlfN-3j#t(N%nvYui<_aQQA|VVub2na4*8*F->ZO z49PRL>4-oB+yTaA{TKLRAeB{J0KAZ9TZ&ZiYs}CF;^j-b_Ciy_^68pcxQ?K3aT*Ze zS{48@4e(GX>Kfk^i~Ddf)qZ3LL`$c7fif#3GV94fUvVii3-7SThx*oxBYynd5lvFz z)l0!VS~vP{*Ie^I$jHmYPok+SG!7b#`qOP>XCss_Tu+@1V``sZt_#S-I4D_e+RBh1 zt*gttwt~#kqs*F5omRDOcI*3?cum0A%}9z;>5ult z$>omTDFatL0|^f_;&J<9K!-h!CXT|wG6dU3^~^Y%j_%&e`j0p+<4hUg{F5Z?Zb=_X zpck#~q_tpTs|h^~J-<|$>g(N8#oE76akC;i9y#5fC3=ctx2xk*ZLnht@qoJUg_$}U zd0?LTrQczqs{R71RrqOTfvDEK*3c;~imlgCTd=c0O5Gh22>3|(MH6#6*~!YyK~gbS01k46euY--cNpc4>KRVKeWfLFh&Wq zMY};$Gw)F`BCiWUWg(ZKJ&Tnbnc%(0lcZG(x!$17WJOPm6u^twf_m0#g48V4%!fpq zmSB@xAL_>{ja~#F@Rrk&z5r}aSwQ^G- zCpL`o9R^{N7a)X=QpfrC??6;nSmo`t7czZ;3S`IJ``(%=Vw~#vw73h!lG8wf7Cz50 zE{p}!F+UuZn>|kx2^F;h!>}-a6c0UB!AX%SapvS4Q{C??`|h*VjJU)UQioNVbw2ZF2>UUVOEc z-En;o1XfB_8;`7pi9(}bNV8`h-^F;2voW!ZY)bCxN>;A$md`6Vn4-t#J3aSiTGq6< z_6SMFma=MG)QArr(g&hzHj{kVL1VdGs(g)_5|G*r%6;K2hU13xV48g`>(T*d;!V~68plfv`h6YqqD5cbJ^kxdZv4@m z`Swk@bbu&EtD%ahMzHlP!)QHg%G(bzs)G%UiwincZ7;q)5bx%?MEh36>&r0dEu(Og z)dil}qXSgMynxlCgSg_9fF?-rNYG#3gJ*(GBiX$o^x4rxqO_6B47LF#zl*!Cy3cJ z=`mQ3Ng<c?k4!x5P&Ww#cTFH3-zQoSOVocE|5NM`6WZH;Gk4RUb9~xGrvp zk8X_j;a5kS6*pg9K$<@w(XfRtpiMKrqY;TCSG(wN>?J!tiFtI{D&AgJz{ii4hV5c^ zxe5ZZ9W$O+`CN{7Ygg7d{HTOxG7bOm&fLzN;HBeGc32W?a@m7Mq{hP{;|#^`phc7p zDC-#$7&&KDZSyEqi^S;nD*yyPs+%`Y68eBql2Mta0$qPF&(a#V#39z7)AE(lU0~Ms z>3%sN>-B)-62b-zm_KmD?6pyFe;iA>3$U}H2);`n%aEp>+$J4hkO67Tq-9?}aSQZ5gdU)9@i_~Uz#`I}Xx)c}*jI3LO`ffdA z@8~wNz&IuM^OZYghQbv_41T%TxPuw47VKkkT)Pc&!&XGDYi7<{NdWnEuXahg;+QU! zR!N*1YGAVoyF~wTb128JUxMPV9}v41vZ%+OTvfsu02xKe5=Z2O6-AI(Uz4zFSu}4c zdqd49yX4$@a}B5x8hTqkDM`OZ|Qy-3t;hP}dyyS{2NLc1|Z_06)f!xrD70+>d~-lN!=Al zK2A5brlJR%ebEzoOW&pDV;#|uDyE@9daf^PKOY5pp!>VO?jM+*X%m-7?c${S=BiQs zlrihN{{d%gy^fT{fygAkdpsrv63xhXDf}=R4(Bn7&wM4nQn(tC#%|QTP1#md1=tLJiE9OUxqEBhn)bY)05t1A5i+jcLu5*1oi@9))$Rft3fI4hF+Y?x) zob%KLjb*&$(Sx8zpw|pi;J~HuFWF`dyK{!M18a#E9^{vq{AFu*MW98-WMyoYItJ!> zf`tpt1Vx}<$fk9RgJ_{|KMg9PyLBXJiBO=x=KRf9cJu6Yl3hXyE}WL^Fd?x}YAHOk zf}OmbTv-{e#sL(|I^O(~qg`3Bfzn;iU%fLRTO}g{Yq`Xj6A9W7SGO6@1RKKa&THDg z3q!yZY~TDi!|HvAj1}5VVu}4V^`4msjb7~Jq$qaV#RiRW$-!|fp9KkJw+M11;IQ>jxQAv&^1$E8TXrMC zZ^S!03T&gH>Q4Kjb!QZ(?XSG5=|;e0I95RSeZ(@G5<9m= zicMrYyKcw^j%&SMUl!}lG^|vsF{ZsMwuLN%nbxdZRWGkrR_0x#vwZW6J|STNK_TG? z98FD8ZXABsCixgy63Y(}i|s9&m2C;qSGnQkl_29CiM6a*lQ^_RWC&3K#!*f>Hdeu^ zqVHrBE7!Z0b6)KicMx!j$d?lG0kW*#Ldn6U9vpcCW%U@`1ja$!`os?HcJ*lRR^XEF z4O){KhVqt{@lSf7Z0#sS&d;YTGgJ3flVoc#@df-00=7?aXx?V>5Ct*taMnJSp*wqm z{wCss)-#B8fi^P^jUm?5IBJ;AB98C^wcULd6seuD$T@kzF8Fqa9FOY6KE{S0x7=w;YhuXA}~>bGJR9(hCUC+-y@zWDkiNZ znenfi_R7|i_2Qa3lum@o+G<<7fV^6ai~dl1+=InP&6HS=frX4xpkz3YZl>n0d@Zac zAZ(O~5QUENWo)-L4B6)~bnUH~M;1k%iHE!)To*VjX+7NoqfEbUv7~Atz(-e3tnQR+YRziV3lE2gSOt^0R-wXWr7 zGmE{hmBIf7T>V{g{ja+9?_S)$>$3k@p8X$S>t9XT|B$mkLD%1y>wi(S|ATk^t2Fzw zA;!w|vp?{!7Qz3Xef?d&`yULfjfJB<_1_@uXUX^P7Vm#ndjEyP{);?r`MKQgxAgv( zP;PJP@^=Hxe?g*zT-dAgnz^*|M%o93j+)NALML>8i)hD@&XoadFpquYPuRI z0|5w3J`nk@h0eixrS2e)L2Z;4qh!i)?$otwv5}V5lC=E}aG>ZoJU-t*0~*0h3HBmg zrEa=vAcrrE>OPh~xY>jU<2J_HE;lZ{9^0HLa@`)XKe<_h?Gwj%B5?1>Es(G$Ue2ec zek_#_zUf#4fR6sDHsjuG+ANiu@YlREFOvoIgg}vcgMB6B$w1!j>6}$hk6w^yH5Ruk z0`*l?|MRcyd`<4nK+sy}NEfGwTS!t|PcT~8do}Ci9Bw^2?A{;Urvc4x8Uo? z`^yc_J=s26AvRT0-DUn99ejZ2fbAY@NUqPpZL$aru8Eos*@Fl35PE(bP35X6XASdS zvLTPJ;ryOS0ne4{Xq}#KwIP`3Ntbu_DzL|wjkbmgKCgwt_(c(EY;G%eHS`1H7@uzQ zLR%eiC%tK(#vng`)48dIT^J>IHFu?N?tnahNksx@Lr84_kxXx~<=pF>o*B|T6=S-Z zIAxBJJba1sCrT@42FK{9CQX?qB!^zSKZjt3u+LMGT_~4V5;-SkhQ=7QB2bZAD6d~K zK4)eI!03-AY?-ql15uKikf2x0T(U@#bRyUxs4Ck*YF-Ls$bVk|Eb~_EJBL1pQ_4-4 zLqVFhQ+%F~9+rO)!l3X4QKUl{DAY$z_BugrE_{V=M-Yis{!vhtlfum-Z$(}E(01bedP;)RIY$ngyCr!-l zZ$QkxqX(748V0fcsoi-XZ7@`55kUb^nj8w;WBMxRSJSQbyAxg*eF zUNQ6`!|gW-Gi1Nodhu9!(pk$MpW-Q&x^3)DF||LBTer^XFSaMFOQek75;K}Y_r+T$ z_J&WMq4~BF=uJo17#+V!OAo3sE<||LYD#YRM>Rs;CXj=OKRbDw5Vg%yj0yGJOl>VZ zHjA?ewKU83B=0iiB&d2=N2<~Ad?{x{m5kgt(<+dMkWTZxLLI%Gintr8#`)7jCRdiq zDmL{8HH)swtUBhWVa(fh3qJ}n6T-lyuT-wLuyMtgmC0HD#MP=hQ5NC8hv=2a z(_ch@ocrYUz3<`*Y2OJOCW&XHaza{`*WltMaP0#n>qehv0VT zGNcWJ)CU%RFwRIt?GTAUtsyfXl9mj#3!!H%Fu}_SDG8&8kNPSJYiD!C<3{R9xW;c{?xc!8hwP7?=_Ta;Q!Y=E0p z!2dUHTSWR-=g<4Z=GLSt0Yccr397^ zn#eCCf1>rX^3pPWj~~`rF=3ou{$W2MyVkEa1Xmp*2HvRb-m~bq(?=&oVswA5oB_{KqasR~3u4ISK5s7ycjakQ^VL9c! zzka7~*zHI~_9ivnoN(HGEY)|f9;yaW7ntisEkoH7FvhEUnvoDMk>iPiZL#a2Yf2Ej z9L;N;+-|=$gQvmxIlh;lsm%$ZlJ$C}En;Pz`*SOx!$0xq5)cd{75$?J#;3#fGg0`3 zaoHIY_B96vM5*w})JmhgMx%M$(Wwor5Y=;`l1fhp^)}4uN%eeFb^wrC056xzQhq#k86nfVAk5p^`zOwZ?o0LrsPppx#EgHR4cbiMIK@ev{Qa*E4tok{| zko8ISl51}($wF5yhXntCm)$d~hUj4EqdlVs($(z~0p*q*tc$U3R{`UV;x5F^8~jid z1r@?#{%N%tL>L1f(FDr`hH|R8Gl*c=Bv`YgPNvS}IwhO0=3CZZ2U~u{ZVK*+F<1~j z_e;5OMTl=}CGYkPrAQb*L;;DP=fKvLU>1kmzjYynlo9fPi9%VYxC+U4icbT7arI9qWb z8-KY`{OdY5b)A$uzBMM?$mgeqSYZ}(c+eJo#k<_(vB{2JPu7gg(qH*>H})9{6?GU5 zsj>#}q;)F&Bspa{M|d?e>UFAZ;FM#vj4HeDW~;c22k8=V<2or)V#(b7X`)4$XI(9Q zL>oNYwSCca>hxq_s6(TjrSG@98r9%DZ2U}r-5QwGNSYjwc|V6_P%`15Y;=H8l6d0R zV5G)Aa)X)GAD_orrz6>qZr_9w%l1+r0(5zc;eF{Ky;4RIVOl4?fWh_CX2zvQc9Ik) zh#2;E%*i|kXNY-}KO-iWk;ryofB(6T=R&L(lxU#oiD%x`X|o6%cF*XPaS~uzE*xE>35sNVr7~CCr;BpDJfe!<5g~4*@#TL(!TlO zly!=XmeV5h-;fpHZc+#{xFBy;fT}R-(0Fhtv2-Ce~1c2|{>v`kbj;ajbBx>RW`30~h7W1niFer@rT0CYKJ$zeQ(N@&!X^A7Rz1Mc+^g+`NZu^c z!NQg?T;fRu$H&4EdL^l1XP5Z zt7(-}^@b{TaI#kFVu?jBY-`OSOAJx?e@y#{50n67LQ%V+6Z-IE=r8v7A1rloZWtXS zcJetCQSe8+cOGjyTxhDgkGIU1Qd+iPhnXxFD5~<_&JrwC>=o;zjb!NU7v9^h>x{h0J z4y$lc0$6WD#kdMmw)#SyGh^PvyK;J>YAqc6Sy@pP+K zD0z=7pnhtA?cD-#CD!K6>rzMdG`Ed!^YMq4h6rXJL%#Xrq)Tcx@?>`FTqh~rQv^}X z=G>Xc=QQM96Emu7MTO+;YV#3uvabt{!nljFyU|=<3Fs$SavXP)dt8X`W(85Re^VwD zz9n}4l}dEUIcOVA_Ji;4AUz2pv&M2j!v=qN7un7h4j(GvAe==4QZ~1k%+q|cC`oH3 zDVRjOX7h*i-5dAZh;k)&9&Bttnit>p*eTbe7-rG%Z2l1>HH31m$y{orU;J;UipKi* zsXak;WslCI?82|rcwW13eM|*KR7d!xDZFAkBB!WL zFb3FeqxIi*ZeR!`# zCd>yN!MvTJ$bB!N#1Ex&_|Xd&=}zw)#(F~Y;rSO7DqB0Qs2HwPaU_8f4_XLi2WZe} zxvr|LJfD@8kk3qs7->sX_rTGq!Qjm{*?vjy2_sxupu|wSj6s&#?}s-w0R4|cy1cPH zKKB^tHxfh>a*(tEiNs>{`5}XQK}VmxLfkapF?8njmG$l98*#93gxk&<3%)NIm!iW@ zF^aKuqO}U#{CK8-&R-YVUftj=QbO9g!YN597z+`YdgCgan4c=^v#L(BTNk`D{}5st zrDIS_R;Cz)J$^qqDM0Sftz9QXB$yhuwLI}gXpH~J+wM>PjM!n$*HLYjIVK}Tl8TLGP*YFekMGXr@^aPh61 zJ9B_UUp@_H&CV!g;2rO`)+bNam~4e{BX6W|7=-UPwAw=72W|ih&I=gL@Z`hxtSB<~(z`6MvImd)c^&p0B*MhwlDls!5sIB{Q)GXWWwM$Kj&2~> zyc=;y7(izpV40-(Mc4A-Q|pMZLrDBbY->$U3Ol{Kn26sC7^MiQ;&SqUVm;IRkLu#4 zpPZXg2+gnV?luaKa6E~lzj_|k*(yezB(7u{vbaX>ieM~!%#^7YK)oP(K(TTqpgdEZ zN25D=d-2Ta-wF4fh-%!YXc)tFK(@+(}Yq**pj+PouEoZ|i zDio_#65kO}S8lJMQZz_|SyflE07LbK9};-*YpnyteM}^ke1!yQV33FbCTY}Tj^`65?a|rfvGR?WM%ru5Vudcsyg)?YO%C0>7+~LOdkS!{3nOa z=r=KE;1q5**kj8d!Balc&sU;G)7>`aGr%(a$!E5^?jZsZtyz(JIhC(~u=3K0%wGvY zOyE<6K-Fbb2)%(o(^$&31TJ@O*=|MV#6FRMQR1DD&#Jl3ge%Vin^GkBC_;$T9Dv2FiK;xN`h- zf*GkSXlB|`!p_6I5##jgCs0uHO(FiRZKC$X9!Qk16bO0PJLC4XR; z;xW&rt>t`)&>}H>R2uxSbxS7mQ0MH?%CG4g{GOt1_4Bm(Bb1GF9D@G*X23Ec9@~g( zNC-J9MknL3M%PnXC&{d;Fx6gun%*xjLwuvJtYjUv`@z8xi4`1;eBx~Hbt~#)`-dmm z6y||ZpT4mh4D-|8=(x~3@wKt5kgeqKCXbOxq)xAT@yli2H)2{t^}B>_oWBu~Rh0|P zp=WMbmSQNg;f&RW8I76Bj5Arzw;JIn!!A>MJGw&L45yQhCDa}y#F$zx|L4A+pGMJu5X&PM|a$0cI0XV%x#xw;>)A=jeI zm5Z_5Om^fgTQg_r;G1Ks{s2{6klrVv652qivt{$6Q2kR+W~CsEC3Ts7AM1c)-vVkr_bEq5 zu;|SuHs^XjA~`JYE1LF3VUi7jLe!9D>CNXB@sN6`Ls`kYy>JMHmIy>qY#|O{A6LuV%CO?v8(D7$)w1nlFSYn}zTYkpo2Nd3i{MeaK zb?~K&xE4BYh6jn3Q2NL&;y&S2j=Ds7A5A=TddXc=S5?xq>rbYK%HcU-BgnAyRFChR zSCN-L-3#tyJn8Wee!&9=!w;T#bSPB48GL;^w2tfXt*{s!xWec~o0m`NUftiB1&^x8 z;-MLSgkEq&7{btL!K(tK(pfW(XWZzsvZX$lq2e!ky(?`-H?zPA(RDuV#wyA|Joo$> zUQts%g{-i#h2KVvarta|6u={x2)UY{Up^5frdH)HRlsJ(*j zjSyN~X{rR^@79A_6H9mLSBydCDmu3kG^t~6F3g`97w5INtu&Z=2u_ikSO}5p)+|C( zK<=7&sn0rJ)bRQSGdcz`5qNg@35W>Ml?i{F(+&=imZ}#0`s6nu)H=y*vIZ`w?G<%E zz>15ll>OOuLAT$B)+7{09t{k{~=o<|$k&Sj3BP*OVJ$F7%xCDB?%E z&%b^cnxr=+bl%nC93bV|{tTJWP47m}-!%%e_c7iJo`J+iG4BFHxn7ApUIeD*n*kC$ z`YNSlYA(HO*p6dD`Y7LSZ?E&y)bw%vK3NVUGyqIqrU)OA8hQAJy0U5HHWGwaYDjfd z2{kg(XVCMyR&Th#b5~}>hpKTd7i=`|;9Ym*6)4D3+*;Q|0F2^Pgg<=w>1^{&Fu~R^ z0lMU7mhq{}2~)<(R?Gw2W>_f;X|S94D(w%NML(8wf$-b?N%MJ*Joybe-?leDFZ#nz z3FZJ@gmhx50E(Tp0ZETUo*<7+GS=TF%|<~Uo+6M@<;NPtk1lC&Vtm!HPqKnLy}7g2 z9X-?Rvl{v&$kw??VRx9y1Bx9ZQRi$bdvy2?vnra#lQV-R`6N!mPK+`$WkYJC? z)+n5SlsLxwmv|OBHcC^WKXZuaGir0gTb(nj17q4rhx_?QnMS90dP7f#nm1U@t|+4< zEmib&G?makG*K*b^jRzvu!^p;2zW;1*S@#B&DnA#24OhyiDpKr^VY|@kXv4IX$$2c zUW#wjgsJtpP~hx9j6Gr5_vbFn-;XsBW`?fxLa~O67*FEN`6i1gMBkKk#2w_vb@Zir^7#okzj;m#zX^9EZ|msG zJ9t^3fr(Mp)!$Xv+sD&k;v>lheynC_6fF(&DUN!cDjg)0;jlC5>@@vYPGak?$pR8t1c#4q3>DEtrkik(Y3=;p zAKap54^)dS{L{72X>lNdgnQmon-$%s?&@t5PjdNbD|gt)QH`CAVreclxzV@h38Z8V zQ5{s1)-mL=(2a7EUeaUmx ziSm{+d;)7GsvoW=eWsst_YbrqOzo#|Im1A-8VIS(5&7>vLi_UiQz^JM}PC14v)Rc^pTv+>GevarZIqm2X$p z%Sb(C917ebvNr5nA~Wa2R6XJ5d++VG@9e zT4U7l4cDZRo;aATW5^+XFP;=Rtch}AfTffmv&FE{Ux8p3WP{YO98!(;RnoTAwx0Kj zOiMk~6z)WA(#1WP$X;2zQ62?PSAHQOH8yv(kx>`B^$m8<{xH&!QXdb(qi@PL^LO4z z*SE$qS7k*c*_m(08|!rc#9U8Ii{D^IrT`}E@ekyK^{J_-1WQ*J)AQ2_b+v4;)rS-e zD3W&aG|m4=+?otT%_H?q7Y3#oPqsKLw>xfTJ&9+3ix<1Rzs7l-YSd9TTuYOMfn|c$ z$H;4pdn_SN z?VSyWg(p_w`X^HYieaJ&-+q98I#~6GW>XIyp}`buM_X_thmYs(OnbNq5(+flVTH>$ ztspYnDBwoXM1}S&Rh(*(rQf~ep#O{TWLNu<|1QQ3|5H5WUCta&xZ0vM?s{TPwIKip zZ|xX=b>|Qh9T9nSdaVbC4s8 zY@6rIbT#UjGrOe|BjQ{~GqBkHL)`VSbJ+)f-xR#3%if*O_zLx|B=GaVJh&gZABZ;k zh!7^stkz`y96e}F{%8yYJTbJ!9BprLYVNF1*)zK*e^1`Yq~KX!vljX$0v!p)42V~Rm1 zUED7Eo!BVzszMF*LwJYaM`o&&g!$65=K^{PwMTWMu~D?UZJdo%v10W)ToXgSMc(bc z??pduj`TO0Nl1Dm9;(DAO3Y?IO(ZT@{_G`cG;hvU6^UPyGAyYn-E8O_P29`^5b?X! z97&H*#r~EIg2~Nvd!HF)T@BsAS7ai!NPgUT-fVq)FtpH+JuJ;6`(B|!x`Og;6|urM z!4Fj3ALmD&^n>)Hk9aYh5GvWs>G#jh@~X0!dogxkej?Owt^a)y5hIAS=j|d3r1Qn6tiJEr@&RG2_5mAZ!@zh$U{t9l43w{z$j5<< zlq<`YPj;>wm3B}&m1@vJjCFbRIrXir=GXH`6)qU32Wd4SJ+VG+f;VP7*#ws)jAGIO ztt0XzMq3~Bq{|-;MsD>jNCE>d9jVOgr8D&Xc%)QJX|XFvsP9*-BRypqmxabIHm5{p zYs<(=f@GBYkNr_Y#=;oXSe{`6X5kPEohWN=tfPuLA_+G9ma(5gp-u3_IvM&}(M@Xl zF_|5M>v~D+6fwVe2JeF10s^Z^>kuO- z8GcXTvp#!B*s}AmG2omh$L_}NayIGYY}(oKwsS7zo`ipYZ>irKrr!%c^Muc!vG}n> z#P!W&1o^u6Nw~ss_$^ZRhj$C2-9IFW32k_8iWC9}yspcAR{CU8hTpPRD4@JEB{XF> z-1ph>-h3n>Cc%KXO&7zyoLa0PAu0)le6NN@H_DB*FuqTz3P~wM9$?YTOxJYr^?Nh! z+imIU&&1aI&U?<4>Ob*oI+tWbAX{IGM!PZsM~3cK(}!`!mFL3*idiLr6L zTNTr_0vjjMS8!aA^ouU3A)M4QG%Tbc8agQgG}L=vU+g50vII5T1`X{DwZ;90aBJJP zD};;Nw5N8c{Y%%P$BpG%q)}8;hq3GSpZB~!?>%`xeTC`UIO7eNwSnK{hyI?m(Uf&< zr?3qbNlc6%c@E9Y?-sG=_9S&3SL{u`uo+;gkbT3u95!ypJCggHa$QyHYvnv*aVkh!MIGmZsxll61jS+IJe|M~=S2{^6;QuJ$xA zBj|Bcm^BX14Yj&8=agZ&g_WAA378m-vW*N}>^4oj{HnQF@^TZYzE$KAQH}uTmksf^6{9lDE!RQzTP)(~zPi)h1Vw+o+f`T{f+@5BfxXh1?JX zpEsjCQtYh40t!9Pp;m^A7f2A2qo7p;HUu4O@NDsXu6Rigd8pN&scy@8a&EINldXpSj>P7H6`v9CkKrlX)Ea5yZ4SfZWS|AwQ$l>@F?j! zp7Iyv=5LqKmqN_(-8lIDDRNbBO$vs)48UmU;r;8DnHo}Ir>N9Q8dlEw&gY{XH^+Uu zGK~u};~;a=qQePEG8ALdUn7Y5LHta;B5<@aDpY|G-#c-lgq-^t6p`Nvq#-a0)3oWP z?>lYKwNFeZc=U8VDu%)DxN=^gDQq9o*;0$rTLgW$=RGB(Gx3cm28!58?nIjle2gBE8I?8mY^*}h^H zIRe0D9m~jzN1k?X`PJpPBpTy-V>8T`BF8*57)VLomu^_6`G-e9!+RI&;gZpZ7p?Bv zG`mJ}Io*KJ|mvtCC*n&OI(K1Iks{hWlwIph ziA04c%?%+SyagUz=aOvYg$ckZnJ|XVZoA|ekPrEnf_sMM{vuz*NgR%X2potMVSa64 zZDEEyD({y)I0PMq7|MpxrnM4zqQUsb9C!(&b|sSLcig4(~sHvZl>j4!J)Jx*Up1M zRlS>nb8JZRpyDgwjjzzkL3Ep46sE?Bd2U2CY58Ru@i z-Gz_x=GX^4f{6CjO2l*sC_ho5*rY`Yem&-5^*K9IR_m?jh3ojB>?I18MWE5(Km>R_ zm>FOzJkY&~xX=GKFiEN?oIrEr%h9jQPvo*gEIVt6i(WoIIg^Y95mKsj{%NtSTIQ-q1F%q3G`hjdE;J3@n)wii_8v}kG zo-9HtSa%ac(eerk)p#@N3!BX^r(KNlu(8;fK>x;V@)bTTlNz}#gIkjeSLjbwT+FNB zukb!?A+0gjXZJs_K7I6){-BbJK_&ZHvjxXRuwHj21DrL~o*R1!8f1Jb1|Mr#)c z@U4T5P+Hq)?Rvt6Wjm-I34QL0hvI_pAKwuAVK2FNpCfvt6sBTi4Ajs@R+(k&`i9rp zWA3e@{3Nk=;PR%U!eyq;yCi9Fn05@sGg;*?06lV9n*~*wVgVk` z&5G*Ao4Yph3N>mnZN?n_dUGT+x`Ft-)j+i;po_?&->m zu=pWuGM(B%3Iz)Z6S4S@Wom61d_@ZCNGtlArYU8c^kJM&rB3~x1`ataGGDR^>vB4c zrw<4SzL@myr=Fi5@A5|^m=6~1oChu=e%7j^T@&=|n_P%LHv!oNS=b)8M-+!~Y`)9k zCkeCkEeQ-0pEzm6G(+)eELQ7w33X!26InqJzv2Y^g1LH*N-I4i-3e6SnjkP;t`gbR zgdwTeoM!tGT9g9<_}Taw@5>@lJi~vRbQ^W<8$uuoNcv)Cti4+>d2^F+>I|jX)H>8e zfqy8j7qTuoC~>faGS@hrC|V+%89&YF+SpKE4|f*s#>Cp$HEh%^H-b3W90#^g|ERf{K=FT|-+$-1?7 z5ov*{yH|zqhH?S)C!L`vGd`8v^3?lYk45)plTn3@p}F!BD#0^tMrU4A=>S-X6ojI8 z9ioS!nh*`$ld^3KK!8g4Zq1H@2d|RTJki74S<)1l2Muv@!qcgWR;E+A#q!*hEbIkjVn+a%NS1iJEY5atlu8&vMm-zl zbEQiiILd+rQ;D?tom=f1zE=2oGq70Rjq>J{nyk`g<6{2Mv80(JVUPJo&77FNF;%Q) zS%v8IhFLjA0&E`SezN01%A7P?c{au$PnLIiJZs|9879n!7UaKRK4KsL=FTVyf!-ZW zZ}k^`z~LR6%@BUhIgP|Gzq8qeEUUb(IrT1tQ6LKW6lE24S*EaxspRhFq6;@Oq}?cO%sD}mlf_6PZ;apCt! zt1W=Q{Tdo{51c>iW;&(hY9{$iHi(K8`*~kGQ@^y7qN}pJxqLrs*Y||^MlhvRRujhC zSy4+>t&bP|fu-uR%AsYD$h-ikn146;v3VoP-UGL@oOWs9&N?{ZdWh(9d)x$Og%3ef zBeqxvP1hdAW!f6bw9Db_Qwi0gsV#8PgM7r9yyZy4gnA_&&H$oNzn<=tys_O;XIpK?qhhjm4}!&e9|5mMbBT` zgk7X#KZq3^Xyg)kIk%OaxIL2P&QX-G*DwuN`;~8I^l5)En*1Xw^Gz7Kw=ti1w87V# za@HP0){%5s7xR^L z6g@<+RY|av5;r*)N?_`Kbi$^!yiQ0|#SgSz8qkt|<4AJ)+l{U~K`bP;Z+%p^8bXa& zmS{3SNgLv8Q$yiY%1m_n`;n-W>phu}sQzaM967-a{d=?v7ptG4Ng>35c*eXcv_1*L zqPZCrOCzIFUb-5p25Vy+KC5o-hx`QBP3lsrYhTezI0PIZ2bDT`hwkrMw3(ju=|RoU zYQ;F+MzEmKL-(@htDg>ct0K4#PS1Hq{%7zyAFDBaMyZAEtf|jc+E6&)0+_pZ3a7Sh z-x26+);jU$+1p$#R>RIjA3Q+*3UXr9D2?qAL0$Su0nY@w828+K#5>>4bGfFT$+Sq? z+)BWtl`@*xgH_aGQS+X5Io7BeG@xY)bSP~7o-28Ehkn?wKz)2{UcVSZPqR9 zPM>;`yzTWN2k3uVnWU+pEGs6?AR?#C{x59(|9!~8LdppQBVe)q8&WQ?@Si6f;LHC+ z!1`}Vp#d*Rp@IK8Cowm4_+a8h%PhofYh&^ckg31)|NEK4-=61xUMvFy0RC1h^U>JA z(8$2p*x30G5&()-#aoi?DcbOkY~=USxv=&G?i|rQGbb@)84c3yaFt`oL z@^i7BY3Xybf-B?OPqZek2R11))n7NRZ%N|BI%J;OyY4>o6Rh|Y8slYSc#p2sd+~2v z7oencJUPkXjPgap)W%;`=GqJLQy%|6X1Oi^4*MWt9_vU|H{--yC zZ~SkM>t6aZ+g|{0|A+Fy8q6l={C9umWcjB*|Lgky`wWc(0Q#Fhn?w&OU_f97n`zZU z0aaOX;kYxPV8%5lV^M+c=M4OFXw$wcH%ifs2ihDqAiQ0Wa9)^Hj!u*x1_SA@cBA_l zWCsPn?Iws~pkr(Zu!wjFqT2pmJ~`dGKUbuT49KP(#*i6?udO7`w6LNd)?2N{X);*2 z@50V_!;f=gAJY~@xqXMXy+9}Kguj-%+tzJnl+Y|kk@V{z#L50$OhQ_k$R?d+L~Em# zlw_&%Rer9PD@$sKvT8;8bVF!JYaw;+X7o5ZRes_$?qD-TokoG;fpdOPbn(7q1iBDv zrCOC>G++}=`H~rZjSL$Y7uOgStB=~nzrP|Td)bdnv7nWsBwjy8X>B*}$jIDOGh6zO z&D5bGSfU?hC+X8NR`0syd-dOU3m}UgoY8Po+5!Au9Fv!Q)1P44rMo2TjD+%5vOhGa zV~d%rq{ujaHLw65f|$NhZ!#uM9_Y@g$OzSR%s|08bR3gFo<2-~vfkc4mFO$7m`zu`qeSt8%MuCZ;hvVJsT;e!JZV+im1yHn9^s2@w zZ5|=MV}P*bQ7vN&gJ9zf3hr;+$gH+=I)cK0Q93j`+r4<8Q7(b} zW}^VSHg^`^+-(TE^T$ zjfP+5a*NymwddCvcYxZ_h8TJy`q`K7=}O+be4RTKgZgYx%Kfi5GQfYGH7h#U8asb9 zaiEn@lqHogw{>*-4^|jb&i5b=I`HaqFfsh=q+8U`$%Iyv8^FTC$-)6(XJrQffUH{J z26DE>|Lx!Yb+fGlt+9)Roudw|r4Ah_=sg!Zsfewe`@af8Nx}VrNRVq^T*DtB`*F@*fQGBz#%=tVp5O{}kE zKn@PDL0|O8#__Tph?R}wr7j#GR^Uq+hy`4<{GvZFzYdtw^PgjYSOBcxj|Tr3mmR!0 zgTXrgAp>%;v%h>cFncE8=qkAs!#rG6kFE8ELw1A>LG?gasXuVg?L&X?n|fR%o!Gb;=5uMOj$y08QP+8X~M z1F^BNzq}U+0=wCZc3>IjD;W^*YJR~ms9-SEe~!V*&hg6jKoC3EEB!zmujdU21Q)u$ z?2i+?uKpjNx%Wdk7CEH8Bi_s8~9KLB_$ec28O;sU=X{8JaO?Dd)f za)7`8{j(jozgO+p*f?L$IXer-OB(>e&ihIRAReq?5rs2fPvcHa3Kz86O9|Ezxj-4PN8}OyCfj}&;^aF8#7uO42*jd?Ljls?acoh%Wf$Xo= z8#@=rs~F1xcGH)04)zw7mo@;0x>xH6>>0q9{tRU21n(R#?8nZ<{?bpuuFLgW_HXRW z`8vJ;IYF=d0LTU4dKGuM*k0KM#KQ7wj{<>1_Deg1SOFX_Z2$s?-IulkbFBT{|C}5Q zEv!r&;CXpTnLk)~n1J_MQf6gaTPITRA}3WQWtOrrwIu~RHmNfBonDkwn^zpb4rJpH zX9I#~2fXe%S%f%QxY#+xS;e_n*jdCx`TyT7bm75p4u7#mJN;{e=VAc?IpC?N#1zEg F{}(THA!Yyo literal 0 HcmV?d00001 diff --git a/c_tools/libzobov/CMakeLists.txt b/c_tools/libzobov/CMakeLists.txt new file mode 100644 index 0000000..ead0a51 --- /dev/null +++ b/c_tools/libzobov/CMakeLists.txt @@ -0,0 +1,8 @@ + +add_library(zobovTool loadZobov.cpp particleInfo.cpp contour_pixels.cpp) + +set_target_properties(zobovTool PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) + +SET(ZOB_LIBS zobovTool ${COSMOTOOL_LIBRARY} ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY} ${NETCDFCPP_LIBRARY} ${NETCDF_LIBRARY}) +SET(ZOB_LIBS PARENT ${ZOB_LIBS}) + diff --git a/c_tools/mock/CMakeLists.txt b/c_tools/mock/CMakeLists.txt new file mode 100644 index 0000000..8761175 --- /dev/null +++ b/c_tools/mock/CMakeLists.txt @@ -0,0 +1,23 @@ +include_directories(${CMAKE_CURRENT_BINARY_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/loaders) + +IF(SDF_SUPPORT) + add_definitions(-DSDF_SUPPORT) +ENDIF(SDF_SUPPORT) + +SET(generateMock_SRCS generateMock.cpp ) +add_genopt(generateMock_SRCS generateMock.ggo generateMock_conf STRUCTNAME generateMock_info) +add_executable(generateMock ${generateMock_SRCS}) +target_link_libraries(generateMock simu_loaders ${ZOB_LIBS} ${HDF5HL_CPP_LIBRARY} ${HDF5_CPP_LIBRARY} ${HDF5_LIBRARY} ${LIBSDF_LIBRARY}) + +add_executable(generateTestMock generateTestMock.cpp) +target_link_libraries(generateTestMock ${ZOB_LIBS}) + + +SET(generateFromCatalog_SRCS generateFromCatalog.cpp) +add_genopt(generateFromCatalog_SRCS generateFromCatalog.ggo generateFromCatalog_conf STRUCTNAME generateFromCatalog_info) +add_executable(generateFromCatalog ${generateFromCatalog_SRCS}) +target_link_libraries(generateFromCatalog ${ZOB_LIBS} ${HDF5HL_CPP_LIBRARY} ${HDF5_CPP_LIBRARY} ${HEALPIX_LIBRARIES}) + +set_target_properties(generateFromCatalog PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) + +subdirs(loaders) diff --git a/c_tools/mock/loaders/CMakeLists.txt b/c_tools/mock/loaders/CMakeLists.txt new file mode 100644 index 0000000..55818a3 --- /dev/null +++ b/c_tools/mock/loaders/CMakeLists.txt @@ -0,0 +1,11 @@ +SET(SIMU_LOADERS_SRCS + simulation_loader.cpp ramses_loader.cpp + gadget_loader.cpp flash_loader.cpp) + +IF (SDF_SUPPORT) + add_definitions(-DSDF_SUPPORT) + SET(SIMU_LOADERS_SRCS ${SIMU_LOADERS_SRCS} sdf_loader.cpp multidark_loader.cpp) +ENDIF (SDF_SUPPORT) + + +add_library(simu_loaders ${SIMU_LOADERS_SRCS}) diff --git a/c_tools/stacking/CMakeLists.txt b/c_tools/stacking/CMakeLists.txt new file mode 100644 index 0000000..fce1f3d --- /dev/null +++ b/c_tools/stacking/CMakeLists.txt @@ -0,0 +1,39 @@ +include_directories(${CMAKE_CURRENT_BINARY_DIR}) + +SET(stackVoids_SRCS stackVoids.cpp) +add_genopt(stackVoids_SRCS stackVoids.ggo stackVoids_conf STRUCTNAME stackVoids_info) +add_executable(stackVoids ${stackVoids_SRCS}) +target_link_libraries(stackVoids ${ZOB_LIBS}) + +SET(stackVoidsZero_SRCS stackVoidsZero.cpp) +add_genopt(stackVoidsZero_SRCS stackVoidsZero.ggo stackVoidsZero_conf STRUCTNAME stackVoidsZero_info) +add_executable(stackVoidsZero ${stackVoidsZero_SRCS}) +target_link_libraries(stackVoidsZero ${ZOB_LIBS}) + +set(computeVelocityProfile_SRCS computeVelocityProfile.cpp) +add_genopt(computeVelocityProfile_SRCS computeVelocityProfile.ggo computeVelocityProfile_conf STRUCT Params) +add_executable(computeVelocityProfile ${computeVelocityProfile_SRCS}) +target_link_libraries(computeVelocityProfile ${ZOB_LIBS}) + + +set(stackDensityField_SRCS stackDensityField.cpp) +add_genopt(stackDensityField_SRCS stackDensityField.ggo stackDensityField_conf STRUCT PARAMS) +add_executable(stackDensityField ${stackDensityField_SRCS}) +target_link_libraries(stackDensityField ${ZOB_LIBS}) + + + +set(stackVelocityField_SRCS stackVelocityField.cpp) +add_genopt(stackVelocityField_SRCS stackVelocityField.ggo stackVelocityField_conf STRUCT PARAMS) +add_executable(stackVelocityField ${stackVelocityField_SRCS}) +target_link_libraries(stackVelocityField ${ZOB_LIBS}) + +SET(pruneVoids_SRCS pruneVoids.cpp) +add_genopt(pruneVoids_SRCS pruneVoids.ggo pruneVoids_conf STRUCTNAME pruneVoids_info) +add_executable(pruneVoids ${pruneVoids_SRCS}) +target_link_libraries(pruneVoids ${ZOB_LIBS}) + +SET(makeAHFOutput_SRCS makeAHFOutput.cpp) +add_genopt(makeAHFOutput_SRCS makeAHFOutput.ggo makeAHFOutput_conf STRUCTNAME makeAHFOutput_info) +add_executable(makeAHFOutput ${makeAHFOutput_SRCS}) +target_link_libraries(makeAHFOutput ${ZOB_LIBS}) diff --git a/c_tools/zobov2/CMakeLists.txt b/c_tools/zobov2/CMakeLists.txt new file mode 100644 index 0000000..38ad864 --- /dev/null +++ b/c_tools/zobov2/CMakeLists.txt @@ -0,0 +1 @@ +subdirs(jozov2 voz1b1) diff --git a/c_tools/zobov2/jozov2/CMakeLists.txt b/c_tools/zobov2/jozov2/CMakeLists.txt new file mode 100644 index 0000000..81c1f1d --- /dev/null +++ b/c_tools/zobov2/jozov2/CMakeLists.txt @@ -0,0 +1,24 @@ +function(omp_add_flags TARGET LANG) + if(OPENMP_FOUND) + if(NOT LANG) + get_target_property(LANG ${TARGET} LINKER_LANGUAGE) + endif() + if(LANG MATCHES "C(XX)?") + set_property(TARGET ${TARGET} APPEND + PROPERTY COMPILE_FLAGS ${OpenMP_${LANG}_FLAGS}) + set_property(TARGET ${TARGET} APPEND + PROPERTY LINK_FLAGS ${OpenMP_${LANG}_FLAGS}) + else() + message(WARNING "omp_add_flags: target '${TARGET}'" + " link language '${LANG}' is not supported") + endif() + endif() +endfunction() + +add_executable(jozov2 jozov2.cpp jozov2_io.cpp jozov2_zones.cpp jozov2_watershed.cpp findrtop.c) + +if (ENABLE_OPENMP) + Find_Package(OpenMP) + omp_add_flags(jozov2 CXX) + add_definitions(-DOPENMP) +ENDIF(ENABLE_OPENMP) diff --git a/c_tools/zobov2/voz1b1/CMakeLists.txt b/c_tools/zobov2/voz1b1/CMakeLists.txt new file mode 100644 index 0000000..98fa941 --- /dev/null +++ b/c_tools/zobov2/voz1b1/CMakeLists.txt @@ -0,0 +1,4 @@ +add_executable(voz1b1_2 voz1b1.cpp vozutil.c voz_io.cpp) +target_link_libraries(voz1b1_2 ${QHULL_LIBRARY} ${COSMOTOOL_LIBRARY} ${MATH_LIB}) + +include_directories(${COSMOTOOL_INCLUDE_PATH} ${QHULL_INCLUDE_PATH}) diff --git a/external/cosmotool/CMakeLists.txt b/external/cosmotool/CMakeLists.txt new file mode 100644 index 0000000..7faca8a --- /dev/null +++ b/external/cosmotool/CMakeLists.txt @@ -0,0 +1,86 @@ +cmake_minimum_required(VERSION 2.6) +list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}") + +project(CosmoToolbox) + +include(GetGitRevisionDescription) +include(ExternalProject) +include(FindOpenMP) + +get_git_head_revision(HEAD GIT_VER) + +option(BUILD_SHARED_LIBS "Build shared libraries." OFF) +option(BUILD_STATIC_LIBS "Build static libraries." ON) +option(ENABLE_OPENMP "Enable OpenMP support." OFF) +option(ENABLE_SHARP "Enable SPHT support." ON) + +find_path(NETCDF_INCLUDE_PATH NAMES netcdf.h) +find_path(NETCDFCPP_INCLUDE_PATH NAMES netcdfcpp.h netcdf) +find_path(GSL_INCLUDE_PATH NAMES gsl/gsl_blas.h) + +IF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf) + SET(FOUND_NETCDF4 1) + FILE(WRITE ${CMAKE_BINARY_DIR}/src/ctool_netcdf_ver.hpp "#define NETCDFCPP4 1") +ELSE(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf) + SET(FOUND_NETCDF3 1) + FILE(WRITE ${CMAKE_BINARY_DIR}/src/ctool_netcdf_ver.hpp "#undef NETCDFCPP4") +ENDIF(EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf) + +find_library(NETCDF_LIBRARY netcdf) +find_library(NETCDFCPP_LIBRARY NAMES netcdf_c++ netcdf_c++4) +find_library(GSL_LIBRARY gsl) +find_library(GSLCBLAS_LIBRARY gslcblas) + +if (ENABLE_SHARP) + SET(SHARP_SOURCE ${CMAKE_SOURCE_DIR}/external/sharp) + SET(DEP_BUILD ${CMAKE_SOURCE_DIR}/external/sharp/auto) + ExternalProject_Add(sharp + SOURCE_DIR ${SHARP_SOURCE} + BUILD_IN_SOURCE 1 + CONFIGURE_COMMAND ${SHARP_SOURCE}/configure --prefix=${DEP_BUILD} + BUILD_COMMAND ${CMAKE_MAKE_PROGRAM} + INSTALL_COMMAND echo "No install" + ) + SET(SHARP_LIBRARY ${DEP_BUILD}/lib/libsharp.a) + SET(SHARP_LIBRARIES ${SHARP_LIBRARY}) + SET(SHARP_INCLUDE_PATH ${DEP_BUILD}/include) +endif (ENABLE_SHARP) + + +set(HDF5_FIND_COMPONENTS HL CXX) +if(HDF5_ROOTDIR) + SET(ENV{HDF5_ROOT} ${HDF5_ROOTDIR}) +endif(HDF5_ROOTDIR) +include(FindHDF5) + +include(FindPkgConfig) + +pkg_check_modules(FFTW3 fftw3>=3.3) +pkg_check_modules(FFTW3F fftw3f>=3.3) +pkg_check_modules(EIGEN3 eigen3) + +include(FindPackageHandleStandardArgs) +set(NETCDF_FIND_REQUIRED TRUE) +set(GSL_FIND_REQUIRED TRUE) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(NetCDF DEFAULT_MSG NETCDF_LIBRARY NETCDFCPP_LIBRARY NETCDF_INCLUDE_PATH) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(GSL DEFAULT_MSG GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH) + + +set(GSL_LIBRARIES ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY}) + +# CPACK Configuration +SET(CPACK_PACKAGE_DESCRIPTION_SUMMARY "A toolbox for impatient cosmologists") +SET(CPACK_PACKAGE_VENDOR "Guilhem Lavaux") +SET(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/LICENCE_CeCILL_V2") +SET(CPACK_PACKAGE_VERSION_MAJOR "1") +SET(CPACK_PACKAGE_VERSION_MINOR "0") +SET(CPACK_PACKAGE_VERSION_PATCH "0-${GIT_VER}") +SET(CPACK_PACKAGE_INSTALL_DIRECTORY "CosmoToolbox-${GalaxExplorer_VERSION_MAJOR}.${GalaxExplorer_VERSION_MINOR}") +SET(CPACK_STRIP_FILES "lib/libCosmoTool.so") +SET(CPACK_SOURCE_IGNORE_FILES +"/CVS/;/\\\\.git/;/\\\\.svn/;\\\\.swp$;\\\\.#;/#;.*~;cscope.*;/CMakeFiles/;.*\\\\.cmake;Makefile") + +add_subdirectory(src) +add_subdirectory(sample) + +include(CPack) diff --git a/external/cosmotool/sample/CMakeLists.txt b/external/cosmotool/sample/CMakeLists.txt new file mode 100644 index 0000000..969480c --- /dev/null +++ b/external/cosmotool/sample/CMakeLists.txt @@ -0,0 +1,65 @@ +SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS}) +include_directories(${CMAKE_SOURCE_DIR}/src ${FFTW3_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH}) + +IF(SHARP_INCLUDE_PATH) + include_directories(BEFORE ${SHARP_INCLUDE_PATH}) +ENDIF(SHARP_INCLUDE_PATH) + +add_executable(testBQueue testBQueue.cpp) +target_link_libraries(testBQueue ${tolink}) + +add_executable(testInterpolate testInterpolate.cpp) +target_link_libraries(testInterpolate ${tolink}) + +add_executable(testSmooth testSmooth.cpp) +target_link_libraries(testSmooth ${tolink}) + +add_executable(testkd testkd.cpp) +target_link_libraries(testkd ${tolink}) + +add_executable(testkd2 testkd2.cpp) +target_link_libraries(testkd2 ${tolink}) + +add_executable(testkd3 testkd3.cpp) +target_link_libraries(testkd3 ${tolink}) + +add_executable(testDelaunay testDelaunay.cpp) +target_link_libraries(testDelaunay ${tolink}) + +add_executable(testNewton testNewton.cpp) +target_link_libraries(testNewton ${tolink}) + +add_executable(testPool testPool.cpp) +target_link_libraries(testPool ${tolink}) + +if (HDF5_FOUND) + add_executable(testReadFlash testReadFlash.cpp) + target_link_libraries(testReadFlash ${tolink}) +endif (HDF5_FOUND) + + +add_executable(testEskow testEskow.cpp) +target_link_libraries(testEskow ${tolink}) + +add_executable(testAlgo testAlgo.cpp) +target_link_libraries(testAlgo ${tolink}) + +add_executable(testBSP testBSP.cpp) +target_link_libraries(testBSP ${tolink}) + +if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) + add_executable(test_fft_calls test_fft_calls.cpp) + target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES}) +endif (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) + +if (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) + include_directories(${SHARP_INCLUDE_PATH}) + add_executable(test_healpix_calls test_healpix_calls.cpp) + target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES}) + set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) + add_dependencies(test_healpix_calls sharp) +endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) + +add_executable(test_cosmopower test_cosmopower.cpp) +target_link_libraries(test_cosmopower ${tolink}) + diff --git a/external/cosmotool/src/CMakeLists.txt b/external/cosmotool/src/CMakeLists.txt new file mode 100644 index 0000000..8901027 --- /dev/null +++ b/external/cosmotool/src/CMakeLists.txt @@ -0,0 +1,92 @@ +SET(CosmoTool_SRCS + fortran.cpp + interpolate.cpp + load_data.cpp + loadGadget.cpp + loadRamses.cpp + octTree.cpp + powerSpectrum.cpp + miniargs.cpp + growthFactor.cpp + cosmopower.cpp +) + +IF(FOUND_NETCDF3) + SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc3.cpp) +ELSE(FOUND_NETCDF3) +IF(FOUND_NETCDF4) + SET(CosmoTool_SRCS ${CosmoTool_SRCS} yorick_nc4.cpp) +ENDIF(FOUND_NETCDF4) +ENDIF(FOUND_NETCDF3) + + +if (HDF5_FOUND) + set(CosmoTool_SRCS ${CosmoTool_SRCS} + h5_readFlash.cpp + loadFlash.cpp + ) +else(HDF5_FOUND) + set(CosmoTool_SRCS ${CosmoTool_SRCS} + loadFlash_dummy.cpp + ) +endif (HDF5_FOUND) + +SET(CosmoTool_SRCS ${CosmoTool_SRCS} + bqueue.hpp + config.hpp + dinterpolate.hpp + field.hpp + fixArray.hpp + fortran.hpp + interpolate3d.hpp + interpolate.hpp + kdtree_leaf.hpp + load_data.hpp + loadGadget.hpp + loadRamses.hpp + loadSimu.hpp + miniargs.hpp + mykdtree.hpp + octTree.hpp + powerSpectrum.hpp + sparseGrid.hpp + sphSmooth.hpp + yorick.hpp + growthFactor.hpp +) + +include_directories(${GSL_INCLUDE_PATH} ${NETCDF_INCLUDE_PATH} ${NETCDFCPP_INCLUDE_PATH} ${CMAKE_BINARY_DIR}/src) + +set(CosmoTool_LIBS ${NETCDFCPP_LIBRARY} ${NETCDF_LIBRARY} ${GSL_LIBRARIES}) +if (HDF5_FOUND) + set(CosmoTool_LIBS ${CosmoTool_LIBS} ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES}) + include_directories(${HDF5_INCLUDE_DIRS}) +endif (HDF5_FOUND) + +set(CosmoTool_LIBS ${CosmoTool_LIBS} PARENT_SCOPE) + +if (BUILD_SHARED_LIBS) + add_library(CosmoTool SHARED ${CosmoTool_SRCS}) + target_link_libraries(CosmoTool ${CosmoTool_LIBS}) + + if (BUILD_STATIC_LIBS) + add_library(CosmoTool_static STATIC ${CosmoTool_SRCS}) + endif(BUILD_STATIC_LIBS) +else (BUILD_SHARED_LIBS) + add_library(CosmoTool STATIC ${CosmoTool_SRCS}) +endif (BUILD_SHARED_LIBS) + +install(TARGETS CosmoTool + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib) + +if (BUILD_SHARED_LIBS) + install(TARGETS CosmoTool_static + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib) +endif (BUILD_SHARED_LIBS) + +install(DIRECTORY . DESTINATION include/CosmoTool + FILES_MATCHING PATTERN "*.hpp") +install(DIRECTORY . DESTINATION include/CosmoTool + FILES_MATCHING PATTERN "*.tcc") diff --git a/pipeline/datasets/example_observation.py b/pipeline/datasets/example_observation.py new file mode 100644 index 0000000..e7518f9 --- /dev/null +++ b/pipeline/datasets/example_observation.py @@ -0,0 +1,250 @@ +#+ +# VIDE -- Void IDentification and Examination -- ./pipeline/sdss/sdss_dr7LCDM.py +# Copyright (C) 2010-2013 Guilhem Lavaux +# Copyright (C) 2011-2013 P. M. Sutter +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; version 2 of the License. +# +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +#+ +#!/usr/bin/env python + +# does analysis in real space assuming LCDM cosmology + +import os +import numpy as np +from void_python_tools.backend.classes import * + +# does full void analysis. Also generates 2d/1d stacked plots and hubble diagram + +# to combine multiple samples: +# For example, combining dim from 0.0-0.1 and bright from 0.1-0.2: +# 1) add dim+bright to sampleList - "+" is required to detect combo! + +# if True, will scan log files for last known completed state and run from there +continueRun = False + +# stages: +# 1 : extract redshift slices from data +# 2 : void extraction using zobov +# 3 : removal of small voids and voids near the edge +# 4 : void stacking and combining +# 5 : 1-d and 2-d profile generation +# 6 : fitting and plotting stacked voids +# 7 : determination of error bars (if requested) +# 8 : hubble diagram generation +# 9 : liklihood determination +startCatalogStage = 1 +endCatalogStage = 1 + +catalogName = "lcdm" + +# where main data files are present/will go +workDir = os.getenv("HOME")+"/workspace/Voids/sdss_dr7LCDM/" +inputDataDir = os.getenv("HOME")+"/workspace/Voids/catalogs/nyuvagc/" + +logDir = os.getenv("PWD")+"/../logs/sdss_dr7LCDM" +figDir = os.getenv("PWD")+"/../figs/sdss_dr7LCDM" + +ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" +CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" + +freshStack = True # "no" to continue building onto existing stacked voids + +errorBars = "CALCULATED" # "ESTIMATED" for error bars from multiple null tests +#errorBars = "ESTIMATED" # "ESTIMATED" for error bars from multiple null tests + # "CALCULATED" to use fitting procedure +numIncoherentRuns = 100 + +ranSeed = 101010 + +useComoving = True # if true, convert to real space using LCDM cosmology + +# try different portions of data. available options: "all, central, edge" +dataPortions = ["central", "all"] + +numZobovDivisions = 2 +numZobovThreads = 2 + +dataSampleList = [] + + +newSample = Sample(fullName = "lss.dr72dim1.dat", + nickName = "dim1", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72dim1.dat", + zBoundary = (0.0, 0.05), + zRange = (0.0, 0.05), + skyFraction = 0.19, + minVoidRadius = 3.479, + fakeDensity = 0.01, + volumeLimited = True, + includeInHubble = False, + partOfCombo = True, + isCombo = False, + comboList= None) +newSample.addStack(0.0, 0.1, 8, 12, False, True) +newSample.addStack(0.0, 0.1, 12, 16, False, True) +newSample.addStack(0.0, 0.1, 16, 20, False, True) +newSample.addStack(0.0, 0.1, 20, 28, False, True) +dataSampleList.append(newSample) + +newSample = Sample(fullName = "lss.dr72dim2.dat", + nickName = "dim2", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72dim2.dat", + zBoundary = (0.0, 0.1), + zRange = (0.05, 0.1), + skyFraction = 0.19, + minVoidRadius = 4.723, + fakeDensity = 0.01, + volumeLimited = True, + includeInHubble = False, + partOfCombo = True, + isCombo = False, + comboList= None) +newSample.addStack(0.0, 0.1, 8, 12, False, True) +newSample.addStack(0.0, 0.1, 12, 16, False, True) +newSample.addStack(0.0, 0.1, 16, 20, False, True) +newSample.addStack(0.0, 0.1, 20, 28, False, True) +dataSampleList.append(newSample) + +newSample = Sample(fullName = "lss.dr72dim1+dim2.dat", + nickName = "dim1+dim2", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72dim.dat", + zBoundary = (0.0, 0.1), + zRange = (0.0, 0.1), + skyFraction = 0.19, + minVoidRadius = 4.723, + fakeDensity = 0.01, + profileBinSize = "auto", + volumeLimited = True, + includeInHubble = True, + partOfCombo = False, + isCombo = True, + comboList= ("lss.dr72dim1.dat", "lss.dr72dim2.dat")) +newSample.addStack(0.0, 0.1, 8, 12, True, False) +newSample.addStack(0.0, 0.1, 12, 16, True, False) +newSample.addStack(0.0, 0.1, 16, 20, True, False) +newSample.addStack(0.0, 0.1, 20, 28, True, False) +dataSampleList.append(newSample) + + +newSample = Sample(fullName = "lss.dr72bright1.dat", + nickName = "bright1", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72bright1.dat", + zBoundary = (0.0, 0.15), + zRange = (0.1, 0.15), + skyFraction = 0.19, + minVoidRadius = 6.763, + fakeDensity = 0.01, + volumeLimited = True, + includeInHubble = False, + partOfCombo = True, + isCombo = False, + comboList= None) +newSample.addStack(0.1, 0.2, 13, 20, False, True) +newSample.addStack(0.1, 0.2, 20, 28, False, True) +newSample.addStack(0.1, 0.2, 28, 36, False, True) +newSample.addStack(0.1, 0.2, 36, 44, False, True) +dataSampleList.append(newSample) + +newSample = Sample(fullName = "lss.dr72bright2.dat", + nickName = "bright2", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72bright2.dat", + zBoundary = (0.0, 0.2), + zRange = (0.15, 0.2), + skyFraction = 0.19, + minVoidRadius = 10.844, + fakeDensity = 0.001, + volumeLimited = True, + includeInHubble = False, + partOfCombo = True, + isCombo = False, + comboList= None) +newSample.addStack(0.1, 0.2, 13, 20, False, True) +newSample.addStack(0.1, 0.2, 20, 28, False, True) +newSample.addStack(0.1, 0.2, 28, 36, False, True) +newSample.addStack(0.1, 0.2, 36, 44, False, True) +dataSampleList.append(newSample) + +newSample = Sample(fullName = "lss.dr72bright1+bright2.dat", + nickName = "bright1+bright2", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + zBoundary = (0.0, 0.2), + zRange = (0.1, 0.2), + skyFraction = 0.19, + minVoidRadius = -1, + fakeDensity = -1, + profileBinSize = "auto", + volumeLimited = True, + includeInHubble = True, + partOfCombo = False, + isCombo = True, + comboList= ("lss.dr72bright1.dat", "lss.dr72bright2.dat")) +newSample.addStack(0.1, 0.2, 13, 20, True, False) +newSample.addStack(0.1, 0.2, 20, 28, True, False) +newSample.addStack(0.1, 0.2, 28, 36, True, False) +newSample.addStack(0.1, 0.2, 36, 44, True, False) +dataSampleList.append(newSample) + +newSample = Sample(fullName = "lss.dr72lrgdim.dat", + nickName = "lrgdim", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72lrgdim.dat", + zBoundary = (0.16, 0.36), + zRange = (0.16, 0.36), + skyFraction = 0.19, + minVoidRadius = 24, + fakeDensity = 0.001, + profileBinSize = 2, # KEEP + volumeLimited = True, + includeInHubble = True, + partOfCombo = False, + isCombo = False, + comboList= None) +newSample.addStack(0.16, 0.36, 54, 70, True, False, zMinPlot=0.2) +newSample.addStack(0.16, 0.36, 70, 106, True, False, zMinPlot=0.2) +dataSampleList.append(newSample) + +newSample = Sample(fullName = "lss.dr72lrgbright.dat", + nickName = "lrgbright", + dataType = "observation", + maskFile = inputDataDir+"/healpix/rast_window_512.fits", + selFunFile = inputDataDir+"/czselfunc.all.dr72lrgbright.dat", + zBoundary = (0.36, 0.44), + zRange = (0.36, 0.44), + skyFraction = 0.19, + minVoidRadius = 38, + fakeDensity = 0.001, + profileBinSize = "auto", + volumeLimited = True, + includeInHubble = False, + partOfCombo = False, + isCombo = False, + comboList= None) +newSample.addStack(0.36, 0.44, 76, 92, False, False) +newSample.addStack(0.36, 0.44, 92, 108, False, False) +dataSampleList.append(newSample) + + diff --git a/pipeline/datasets/example_simulation.py b/pipeline/datasets/example_simulation.py new file mode 100644 index 0000000..242cc18 --- /dev/null +++ b/pipeline/datasets/example_simulation.py @@ -0,0 +1,144 @@ +#+ +# VIDE -- Void IDEntification pipeline -- ./pipeline/datasets/mergertree.py +# Copyright (C) 2010-2013 Guilhem Lavaux +# Copyright (C) 2011-2013 P. M. Sutter +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; version 2 of the License. +# +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +#+ +import os + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- +# CONFIGURATION + +startCatalogStage = 0 +endCatalogStage = 0 +continueRun = False + +startAPStage = 1 +endAPStage = 2 + +# directory for the input simulation/observational particle files +catalogDir = os.getenv("HOME")+"/workspace/Voids/catalogs/mergertree1024/" + +# path to HOD code +hodPath = os.getenv("HOME")+"/projects/Voids/hod/HOD.x" + +# where to put the final void catalog, figures, and output logs +voidOutputDir = os.getenv("HOME")+"/workspace/Voids/mergertree1024/" +figDir = os.getenv("PWD")+"/../figs/mergertree1024/" +logDir = os.getenv("PWD")+"/../logs/mergertree1024/" + +# where to place the pipeline scripts +scriptDir = os.getenv("PWD")+"/mergertree1024/" + +# simulation or observation? +dataType = "simulation" + +# available formats for simulation: gadget, mergertree +dataFormat = "sdf" +dataUnit = 1 # as multiple of Mpc/h + +# place particles on the lightcone? +useLightCone = False + +# also do peculiar velocities? +doPecVel = False + +# common filename of particle files +particleFileBase = "mf_4s_1G_1k_NNNNN" +particleFileDummy = 'NNNNN' + +# list of file numbers for the particle files +# to get particle file name, we take particleFileBase+fileNum +fileNums = ["1.000"] + +# redshift of each file in the above list +redshifts = ["0.0"] + +# how many independent slices along the z-axis? +numSlices = 1 +#numSlices = 4 +numAPSlices = 1 + +# how many subdivisions along the x- and y- axis? +# ( = 2 will make 4 subvolumes for each slice, = 3 will make 9, etc.) +numSubvolumes = 1 + +# prefix to give all outputs +prefix = "mt_" + +# list of desired subsamples - these are in unts of h Mpc^-3! +subSamples = [0.1, 0.05, 0.02, 0.01, 0.004, 0.002, 0.001, 0.0003, 0.0001] +#doSubSampling = False # do the subsampling in preparation script? +doSubSampling = True # do the subsampling in preparation script? + +# common filename of halo files, leave blank to ignore halos +haloFileBase = "mf_4s_1G_1k_bgc2_NNNNN.sdf" +haloFileDummy = 'NNNNN' + +# minimum halo mass cuts to apply for the halo catalog +# use "none" to get all halos +minHaloMasses = ["none", 1.2e13] +#minHaloMasses = [7.0e11, 1.2e13] + +# locations of data in the halo catalog +haloFileMCol = 6 +haloFileXCol = 0 +haloFileYCol = 1 +haloFileZCol = 2 +haloFileVXCol = 3 +haloFileVYCol = 4 +haloFileVZCol = 5 +haloFileColSep = ',' +haloFileNumComLines = 0 + +# adjust these two parameters given the memory contraints on your system: +# numZobovDivisions: how many sub-volumes per dimension will zobov process +# numZobovThreads: how many sub-volumes to process at once? +numZobovDivisions = 4 +numZobovThreads = 4 + +# simulation information +numPart = 1024*1024*1024 +lbox = 999.983 # Mpc/h +omegaM = 0.2847979853038958 +hubble = 0.6962 + +#galDens = 0.000225 +hodParmList = [ + {'name' : "LowRes", #BOSS: Manera et al. 2012, eq. 26 + 'Mmin' : 0.0, + 'M1' : 1.e14, + 'sigma_logM' : 0.596, + 'alpha' : 1.0127, + 'Mcut' : 1.19399e13, + 'galDens' : 0.0002, + }, + + {'name' : "HighRes", + 'Mmin' : 0.0, + #'Mmin' : 1.99525e12, + 'M1' : 3.80189e13, + 'sigma_logM' : 0.21, + 'alpha' : 1.12, + 'Mcut' : 6.91831e11, + 'galDens' : 0.02, + } +] + +# END CONFIGURATION +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- diff --git a/python_tools/void_python_tools/apTools/chi2/__init__.py b/python_tools/void_python_tools/apTools/chi2/__init__.py index 5d9f6d2..c85f49f 100644 --- a/python_tools/void_python_tools/apTools/chi2/__init__.py +++ b/python_tools/void_python_tools/apTools/chi2/__init__.py @@ -17,6 +17,4 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. #+ -from velocityProfileFitNative import * -from likelihood import * from cosmologyTools import * diff --git a/python_tools/void_python_tools/apTools/chi2/cosmologyTools.py b/python_tools/void_python_tools/apTools/chi2/cosmologyTools.py index 8337a9b..9b01534 100644 --- a/python_tools/void_python_tools/apTools/chi2/cosmologyTools.py +++ b/python_tools/void_python_tools/apTools/chi2/cosmologyTools.py @@ -24,7 +24,7 @@ import numpy as np import scipy.integrate as integrate from void_python_tools.backend import * -__all__=['expansion', 'angularDiameter', 'expectedStretch', 'aveStretch', 'aveExpansion', 'aveStretchCone', 'aveWeightedStretch'] +__all__=['expansion', 'angularDiameter', 'aveExpansion'] # returns 1/E(z) for the given cosmology def expansion(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): @@ -41,114 +41,7 @@ def eosDE(z, w0 = -1.0, wa = 0.0): def angularDiameter(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): da = integrate.quad(expansion, 0.0, z, args=(Om, Ot, w0, wa))[0] return da - -# returns expected void stretch for the given cosmology -def expectedStretch(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): - ez = 1./expansion(z, Om=Om, Ot=Ot, w0=w0, wa=wa) - da = angularDiameter(z, Om=Om, Ot=Ot, w0=w0, wa=wa) - return ez*da/z - -# ----------------------------------------------------------------------------- -# returns average expected void stretch for a given redshift range -# assuming a cone -def aveStretchCone(zStart, zEnd, skyFrac = 0.19, Om = 0.27, Ot = 1.0, - w0 = -1.0, wa = 0.0): - #print "assuming observation!", skyFrac - if zStart == 0.0: zStart = 1.e-6 - - h1 = zStart - h2 = zEnd - - r1 = skyFrac * 4* np.pi * zStart**2 - r2 = skyFrac * 4 * np.pi * zEnd**2 - - # surface area of a slice within a cone - def coneSlice(x, h, r): - return np.pi * (r/h*x)**2 - - def coneFunc(z, h, r, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): - return coneSlice(z, h, r) * expectedStretch(z, Om, Ot, w0, wa) - - aveHigh = integrate.quad(coneFunc, 0.0, zEnd, args=(h2, r2, Om, Ot, w0, wa), full_output=1)[0] - aveLow = integrate.quad(coneFunc, 0.0, zStart, args=(h1, r1, Om, Ot, w0, wa), full_output=1)[0] - volumeHigh = integrate.quad(coneSlice, 0.0, zEnd, args=(h2, r2))[0] - volumeLow = integrate.quad(coneSlice, 0.0, zStart, args=(h1, r1))[0] - - return (aveHigh-aveLow)/(volumeHigh-volumeLow) - -# returns average expected void stretch for a given redshift range -def aveStretch(sample, zStart, zEnd, - Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): - if zStart == 0.0: zStart = 1.e-6 - - if sample.dataType == "observation": - stretch = aveStretchCone(zStart, zEnd, - skyFrac=sample.skyFraction, Om=Om, Ot=Ot, - w0=w0, wa=wa) - else: - ave = integrate.quad(expectedStretch, zStart, zEnd, - args=(Om, Ot, w0, wa))[0] - ave /= (zEnd-zStart) - stretch = ave - - # if in comoving space, calculate stretch for fiducial cosmology - # and take relative amount - if not sample.useLightCone or sample.useComoving: - if sample.dataType == "observation": - stretchFid = aveStretchCone(zStart, zEnd, - skyFrac=sample.skyFraction, Om=sample.omegaM, Ot=Ot, - w0=w0, wa=wa) - else: - ave = integrate.quad(expectedStretch, zStart, zEnd, - args=(sample.omegaM, Ot, w0, wa))[0] - ave /= (zEnd-zStart) - stretchFid = ave - - stretch = stretchFid/stretch - - return stretch - - -# ----------------------------------------------------------------------------- -# returns average expected void stretch for a given redshift range -# weighted by an actual void distribution -def aveWeightedStretch(zStart, zEnd, skyFrac = 0.19, Om = 0.27, Ot = 1.0, - w0 = -1.0, wa = 0.0, dist=None, bins=None, - useComoving=True, OmFid=None): - if zStart == 0.0: zStart = 1.e-6 - - def weightedSlice(x): - return np.interp(x, bins[:-1], dist) - - def weightedFunc(z, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): - return expectedStretch(z, Om, Ot, w0, wa) *\ - weightedSlice(z) - - ave = integrate.quad(weightedFunc, zStart, zEnd, args=(Om, Ot, w0, wa), - full_output=1)[0] - volume = integrate.quad(weightedSlice, zStart, zEnd, full_output=1)[0] - - if volume == 0.0: volume = 1.0 - stretch = ave/volume - - # if in comoving space, calculate stretch for fiducial cosmology - # and take relative amount - if useComoving: - ave = integrate.quad(weightedFunc, zStart, zEnd, args=(OmFid, - Ot, w0, wa), - full_output=1)[0] - volume = integrate.quad(weightedSlice, zStart, zEnd, full_output=1)[0] - - if volume == 0.0: volume = 1.0 - stretchFid = ave/volume - - if stretchFid != 0.0: - stretch = stretchFid/stretch - - return stretch - - # ----------------------------------------------------------------------------- # returns average expected expansion for a given redshift range def aveExpansion(zStart, zEnd, Om = 0.27, Ot = 1.0, w0 = -1.0, wa = 0.0): diff --git a/python_tools/void_python_tools/apTools/profiles/__init__.py b/python_tools/void_python_tools/apTools/profiles/__init__.py index 14c0fa1..6015350 100644 --- a/python_tools/void_python_tools/apTools/profiles/__init__.py +++ b/python_tools/void_python_tools/apTools/profiles/__init__.py @@ -17,10 +17,4 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. #+ -from build import * -from draw import * -from fit import * -from inertia import * -from mcmc import * -from generateExpFigure import * from getSurveyProps import * diff --git a/python_tools/void_python_tools/backend/__init__.py b/python_tools/void_python_tools/backend/__init__.py index 998fc76..db12e12 100644 --- a/python_tools/void_python_tools/backend/__init__.py +++ b/python_tools/void_python_tools/backend/__init__.py @@ -19,4 +19,3 @@ #+ from classes import * from launchers import * -from catalogPrep import * diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index f961d65..1a56372 100644 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -95,8 +95,6 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, log = open(logFile, 'w') subprocess.call([binPath, arg1], stdout=log, stderr=log) log.close() - #cmd = "%s --configFile=%s >& %s" % (binPath,parmFile,logFile) - #os.system(cmd) if jobSuccessful(logFile, "Done!\n"): print "done" else: @@ -118,7 +116,6 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, if os.access("mask_index.txt", os.F_OK): os.system("mv %s %s" % ("mask_index.txt", zobovDir)) os.system("mv %s %s" % ("total_particles.txt", zobovDir)) - #os.system("mv %s %s" % ("sample_info.txt", zobovDir)) if os.access("galaxies.txt", os.F_OK): os.system("mv %s %s" % ("galaxies.txt", zobovDir)) @@ -214,18 +211,14 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, file(parmFile, mode="w").write(conf) if (prevSubSample == -1): - #cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile) cmd = "%s --configFile=%s" % (binPath,parmFile) log = open(logFile, 'w') else: - #cmd = "%s --configFile=%s &>> %s" % (binPath,parmFile,logFile) cmd = "%s --configFile=%s" % (binPath,parmFile) log = open(logFile, 'a') arg1 = "--configFile=%s" % parmFile subprocess.call(cmd, stdout=log, stderr=log, shell=True) - #subprocess.call([binPath, arg1], stdout=log, stderr=log) log.close() - #os.system(cmd) # remove intermediate files if (prevSubSample != -1): @@ -329,9 +322,6 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, else: maskIndex = -1 maxDen = 0.2 - # TEST - #maxDen = 1.e99 - # END TEST if not (continueRun and jobSuccessful(logFile, "Done!\n")): for fileName in glob.glob(zobovDir+"/part._"+sampleName+".*"): @@ -343,11 +333,6 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, if os.access(zobovDir+"/voidDesc_"+sampleName+".out", os.F_OK): os.unlink(zobovDir+"/voidDesc_"+sampleName+".out") - #cmd = "%s/vozinit %s 0.1 1.0 %g %s %g %s %s %s &> %s" % \ - # (binPath, datafile, numZobovDivisions, \ - # "_"+sampleName, numZobovThreads, \ - # binPath, zobovDir, maskIndex, logFile) - #os.system(cmd) cmd = [binPath+"/vozinit", datafile, "0.1", "1.0", str(numZobovDivisions), \ "_"+sampleName, str(numZobovThreads), \ binPath, zobovDir, str(maskIndex)] @@ -355,23 +340,11 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, subprocess.call(cmd, stdout=log, stderr=log) log.close() - #cmd = "./%s >> %s 2>&1" % (vozScript, logFile) - #os.system(cmd) cmd = ["./%s" % vozScript] log = open(logFile, 'a') subprocess.call(cmd, stdout=log, stderr=log) log.close() -# cmd = "%s/jozov %s %s %s %s %s %g %s >> %s 2>&1" % \ - #cmd = "%s/../c_tools/zobov2/jozov2/jozov2 %s %s %s %s %s %g %s >> %s 2>&1" % \ - # (binPath, \ - # zobovDir+"/adj_"+sampleName+".dat", \ - # zobovDir+"/vol_"+sampleName+".dat", \ - # zobovDir+"/voidPart_"+sampleName+".dat", \ - # zobovDir+"/voidZone_"+sampleName+".dat", \ - # zobovDir+"/voidDesc_"+sampleName+".out", \ - # maxDen, maskIndex, logFile) - #os.system(cmd) cmd = [binPath+"../c_tools/zobov2/jozov2/jozov2", \ zobovDir+"/adj_"+sampleName+".dat", \ zobovDir+"/vol_"+sampleName+".dat", \ @@ -420,16 +393,6 @@ def launchPrune(sample, binPath, mockIndex = -1 maxDen = 0.2 observationLine = "" - # TEST - #maxDen = 1.e99 - # END TEST - - #periodicLine = " --periodic='" - #if sample.numSubvolumes == 1: periodicLine += "xy" - #if np.abs(sample.zBoundaryMpc[1] - sample.zBoundaryMpc[0] - \ - # sample.boxLen) <= 1.e-1: - # periodicLine += "z" - #periodicLine += "' " periodicLine = " --periodic='" + getPeriodic(sample) + "'" @@ -474,11 +437,6 @@ def launchPrune(sample, binPath, cmd += " --omegaM=" + str(sample.omegaM) cmd += " --outputDir=" + zobovDir cmd += " --sampleName=" + str(sampleName) - #cmd += " &> " + logFile - #f=file("run_prune.sh",mode="w") - #f.write(cmd) - #f.close() - #os.system(cmd) log = open(logFile, 'w') subprocess.call(cmd, stdout=log, stderr=log, shell=True) log.close() @@ -505,11 +463,6 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, sampleName2 = sample2.fullName periodicLine = " --periodic='" + getPeriodic(sample1) + "' " - #if sample1.dataType != "observation": - # if sample1.numSubvolumes == 1: periodicLine += "xy" - # if np.abs(sample1.zBoundaryMpc[1] - sample1.zBoundaryMpc[0] - sample1.boxLen) <= 1.e-1: - # periodicLine += "z" - #periodicLine += "' " if strictMatch: matchPrefix = "" @@ -572,1278 +525,6 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, print "already done!" -# ----------------------------------------------------------------------------- -def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, - voidDir=None, freshStack=True, runSuffix=None, - zobovDir=None, useLightCone=None, useComoving=None, - omegaM=None, - INCOHERENT=False, ranSeed=None, summaryFile=None, - continueRun=None, dataType=None, prefixRun="", - idList=None, rescaleOverride=None): - - sampleName = sample.fullName - - if idList != None: - customLine = "selected" - else: - customLine = "" - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion, - customLine=customLine) - - logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out" - - treeFile = voidDir+"/tree.data" - - if (freshStack) and os.access(treeFile, os.F_OK): - os.unlink(treeFile) - - centralRadius = stack.rMin * 0.25 - - # restrict to relevant ranges of sample - zMin = max(sample.zRange[0], stack.zMin) - zMax = min(sample.zRange[1], stack.zMax) - if useComoving or not useLightCone: - zMin = LIGHT_SPEED/100.*vp.angularDiameter(zMin, Om=omegaM) - zMax = LIGHT_SPEED/100.*vp.angularDiameter(zMax, Om=omegaM) - else: - zMin *= LIGHT_SPEED/100. - zMax *= LIGHT_SPEED/100. - - if dataType == "observation": - obsFlag = "observation" - else: - obsFlag = "" - - #Rextracut = stack.rMin*3 + 1 - #Rcircular = stack.rMin*3 + 2 - Rextracut = stack.rMax*3 + 1 - Rcircular = stack.rMax*3 + 2 - - if dataType == "observation": - maskIndex = open(zobovDir+"/mask_index.txt", "r").read() - totalPart = open(zobovDir+"/total_particles.txt", "r").read() - maxDen = 0.2*float(maskIndex)/float(totalPart) - else: - maskIndex = 999999999 - maxDen = 0.2 - # TEST - #maxDen = 1.e99 - # END TEST - - if INCOHERENT: - nullTestFlag = "INCOHERENT" - else: - nullTestFlag = "" - - if rescaleOverride == None: rescaleOverride = stack.rescaleMode - if rescaleOverride == "rmax": - rescaleFlag = "rescale" - else: - rescaleFlag = "" - - if stack.maxVoids != -1: - maxVoidsFlag = "maxVoids="+str(stack.maxVoids) - else: - maxVoidsFlag = "" - - idListFile = "idlist.temp" - if idList != None: - idListFlag = "idList " + idListFile - idFile = open(idListFile, 'w') - for id in idList: idFile.write(str(id)+"\n") - idFile.close() - rMinToUse = 0. - rMaxToUse = 100000. - centralRadius = rMaxToUse - else: - idListFlag = "" - rMinToUse = stack.rMin - rMaxToUse = stack.rMax - - periodicLine = "periodic='" + getPeriodic(sample) + "'" - #if sample.dataType == "observation": - # periodicLine = "periodic=''" - #else: - # periodicLine = "periodic='" - # if sample.numSubvolumes == 1: periodicLine += "xy" - # if np.abs(sample.zBoundaryMpc[1] - sample.zBoundaryMpc[0] - \ - # sample.boxLen) <= 1.e-1: - # periodicLine += "z" - # periodicLine += "' " - - launchDir = os.getcwd() - os.chdir(voidDir) - - conf=""" - desc %s - partzone %s - zonevoid %s - volumefile %s - Rmin %g - Rmax %g - particles %s - extraInfo %s - densityThreshold %g - centralRadius %g - edgeAvoidance %g - Circular %g - Zmin %g - Zmax %g - %s - %s - ranSeed %d - dataPortion %s - barycenters %s - boundaryDistances %s - %s - %s - %s - %s - """ % \ - (zobovDir+"/voidDesc_"+thisDataPortion+"_"+sampleName+".out", - zobovDir+"/voidPart_"+sampleName+".dat", - zobovDir+"/voidZone_"+sampleName+".dat", - zobovDir+"/vol_"+sampleName+".dat", - rMinToUse, - rMaxToUse, - zobovDir+("/%szobov_slice_"%prefixRun)+sampleName, - zobovDir+"/zobov_slice_"+sampleName+".par", - maxDen, - centralRadius, - Rextracut, - Rcircular, - zMin, - zMax, - obsFlag, - nullTestFlag, - ranSeed, - thisDataPortion, - zobovDir+"/barycenters_"+thisDataPortion+"_"+sampleName+".out", - zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out", - rescaleFlag, - maxVoidsFlag, - idListFlag, - periodicLine - ) - - parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par" - fp = file(parmFile, mode="w") - fp.write(conf) - - # continue stacking if requested; pull files to local directory - if os.access(treeFile, os.F_OK): - os.system("mv %s %s" % (treeFile, "./"+tree.data)) - fp.write("loadSearchTree\n") - fp.write("getTree\n") - fp.write("doExtraction\n") - else: - fp.write("dumpSearchTree\n") - fp.write("dumpTree\n") - fp.write("doExtraction\n") - fp.close() - - jobString = " "+runSuffix+":" - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - #cmd = "%s --configFile=%s &> %s" % \ - # (binPath, parmFile, logFile) - #os.system(cmd) - cmd = "%s --configFile=%s" % \ - (binPath, parmFile) - log = open(logFile, 'w') - subprocess.call(cmd, stdout=log, stderr=log, shell=True) - log.close() - - if jobSuccessful(logFile, "Done!\n"): - print jobString, "Stacking voids done" - else: - print jobString, "Stacking voids FAILED!" - exit(-1) - - else: - print jobString, "Stacking voids already done!" - if os.access(parmFile, os.F_OK): os.unlink(parmFile) - return - - minDist = None - maxDist = None - numVoids = None - numParticles = None - for line in open(logFile): - if "Minimum void distance is" in line: - line = line.split() - minDist = float(line[4])/3000 - - if "Maximum void distance is" in line: - line = line.split() - maxDist = float(line[4])/3000 - - if "Selected" in line: - line = line.split() - numVoids = line[1] - - if "Final zobov stack size" in line: - line = line.split() - numParticles = line[4] - - open(voidDir+"/num_voids.txt", "w").write(numVoids) - open(voidDir+"/num_particles.txt", "w").write(numParticles) - - if os.access(voidDir+"/NOVOID", os.F_OK): - os.unlink(voidDir+"/NOVOID") - - if (numVoids == "0"): - print jobString, "No voids found; skipping!" - fp = open(voidDir+"/NOVOID", "w") - fp.write("no voids found\n") - fp.close() - return - - emptyStack = False - for line in open(logFile): - if "EMPTY STACK" in line: - emptyStack = True - if emptyStack: - print jobString, "Stack is empty; skipping!" - fp = open(voidDir+"/NOVOID", "w") - fp.write("empty stack\n") - fp.close() - return - - summaryLine = runSuffix + " " + \ - thisDataPortion + " " + \ - str(stack.zMin) + " " + \ - str(stack.zMax) + \ - " " + str(numVoids) + " " + str(minDist) + \ - " " + str(maxDist) - if summaryFile != None: - open(summaryFile, "a").write(summaryLine+"\n") - - # write out individual normalizations for each void - if os.access("tree.data", os.F_OK): - for line in open(zobovDir+"/sample_info.txt"): - if "Number of real" in line: - numParticles = line.split()[-1] - if "Estimated volume" in line: - boxVol = line.split()[-1] - nbar = 1.0 - normalization = float(numParticles)/float(boxVol)/nbar - dataTemp = file("centers.txt", "r").readlines() - fp = file("normalizations.txt", "w") - for iVoid in xrange(len(dataTemp)): - fp.write(str(normalization)+"\n") - fp.close() - - #os.system("mv %s %s" % ("tree.data", treeFile)) - #os.system("mv %s %s" % ("void_indexes.txt", voidDir+"/")) - #os.system("mv %s %s" % ("posx.nc", voidDir+"/")) - #os.system("mv %s %s" % ("posy.nc", voidDir+"/")) - #os.system("mv %s %s" % ("posz.nc", voidDir+"/")) - #os.system("mv %s %s" % ("z_void_indexes.txt", voidDir+"/")) - #os.system("mv %s %s" % ("z_posx.nc", voidDir+"/")) - #os.system("mv %s %s" % ("z_posy.nc", voidDir+"/")) - #os.system("mv %s %s" % ("z_posz.nc", voidDir+"/")) - #os.system("mv %s %s" % ("redshifts.nc", voidDir+"/")) - #os.system("mv %s %s" % ("indexes.nc", voidDir+"/")) - #os.system("mv %s %s" % ("kdtree_stackvoids.dat", voidDir+"/")) - #os.system("mv %s %s" % ("centers.txt", voidDir+"/")) - #os.system("mv %s %s" % ("z_centers.txt", voidDir+"/")) - #os.system("mv %s %s" % ("sky_positions.txt", voidDir+"/")) - #os.system("mv %s %s" % ("check.txt", voidDir+"/")) - #os.system("mv %s %s" % ("tracer.txt", voidDir+"/")) - #os.system("mv %s %s" % ("normalizations.txt", voidDir+"/")) - #os.system("mv %s %s" % ("boundaryDistances.txt", voidDir+"/")) - - if os.access(idListFile, os.F_OK): os.unlink(idListFile) - - if os.access(parmFile, os.F_OK): os.unlink(parmFile) - - os.chdir(launchDir) - - return - -# ----------------------------------------------------------------------------- -def launchCombine(sample, stack, voidDir=None, logFile=None, - zobovDir=None, workDir=None, thisDataPortion=None): - - sampleName = sample.fullName - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion) - jobString = " "+runSuffix+":" - - sys.stdout = open(logFile, 'w') - #sys.stderr = open(logFile, 'a') - - if os.access(voidDir+"/num_voids.txt", os.F_OK): - os.unlink(voidDir+"/num_voids.txt") - - doneGalUpdate = os.access(zobovDir+"donegalupdate", os.F_OK) - - for comboName in sample.comboList: - if not os.access(zobovDir+"/galaxies.txt", os.F_OK): - shutil.copy(workDir+"/sample_"+comboName+"/galaxies.txt", zobovDir) - elif not doneGalUpdate: - dataTemp = file(workDir+"/sample_"+comboName+"/galaxies.txt", - "r").read() - file(zobovDir+"/galaxies.txt", "a").write(dataTemp) - - sourceStackDir = workDir+"/sample_"+comboName+"/stacks_"+\ - runSuffix - - if os.access(sourceStackDir+"/NOVOID", os.F_OK): - continue - - if not os.access(voidDir+"/num_voids.txt", os.F_OK): - # just copy - shutil.copy(sourceStackDir+"/num_voids.txt", voidDir) - shutil.copy(sourceStackDir+"/num_particles.txt", voidDir) - shutil.copy(sourceStackDir+"/posx.nc", voidDir) - shutil.copy(sourceStackDir+"/posy.nc", voidDir) - shutil.copy(sourceStackDir+"/posz.nc", voidDir) - shutil.copy(sourceStackDir+"/z_posx.nc", voidDir) - shutil.copy(sourceStackDir+"/z_posy.nc", voidDir) - shutil.copy(sourceStackDir+"/z_posz.nc", voidDir) - shutil.copy(sourceStackDir+"/indexes.nc", voidDir) - shutil.copy(sourceStackDir+"/redshifts.nc", voidDir) - shutil.copy(sourceStackDir+"/centers.txt", voidDir) - shutil.copy(sourceStackDir+"/void_indexes.txt", voidDir) - shutil.copy(sourceStackDir+"/z_centers.txt", voidDir) - shutil.copy(sourceStackDir+"/z_void_indexes.txt", voidDir) - shutil.copy(sourceStackDir+"/sky_positions.txt", voidDir) - shutil.copy(sourceStackDir+"/normalizations.txt", voidDir) - shutil.copy(sourceStackDir+"/boundaryDistances.txt", voidDir) - - else: - # append data - dataTemp = file(voidDir+"/num_voids.txt", "r").readlines() - dataTemp2 = file(sourceStackDir+"/num_voids.txt", "r").\ - readlines() - dataTemp = np.array(dataTemp, dtype='i') - dataTemp2 = np.array(dataTemp2, dtype='i') - dataTemp += dataTemp2 - dataTemp = str(dataTemp[0]) - file(voidDir+"/num_voids.txt", "w").write(str(dataTemp)) - - dataTemp = file(voidDir+"/num_particles.txt", "r").readlines() - dataTemp2 = file(sourceStackDir+"/num_particles.txt", "r").\ - readlines() - dataTemp = np.array(dataTemp, dtype='i') - dataTemp2 = np.array(dataTemp2, dtype='i') - dataTemp += dataTemp2 - dataTemp = str(dataTemp[0]) - file(voidDir+"/num_particles.txt", "w").write(str(dataTemp)) - - dataTemp = file(sourceStackDir+"/z_centers.txt", "r").read() - file(voidDir+"/z_centers.txt", "a").write(dataTemp) - - dataTemp = file(sourceStackDir+"/centers.txt", "r").read() - file(voidDir+"/centers.txt", "a").write(dataTemp) - - dataTemp = file(sourceStackDir+"/normalizations.txt", "r").read() - file(voidDir+"/normalizations.txt", "a").write(dataTemp) - - dataTemp = file(sourceStackDir+"/boundaryDistances.txt","r").read() - file(voidDir+"/boundaryDistances.txt", "a").write(dataTemp) - - dataTemp = file(sourceStackDir+"/sky_positions.txt", "r").\ - read() - file(voidDir+"/sky_positions.txt", "a").write(dataTemp) - - idxTemp = file(sourceStackDir+"/z_void_indexes.txt", "r").\ - readlines() - idxTemp = np.array(idxTemp, dtype='i') - fp = NetCDFFile(voidDir+"/z_posx.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - idxTemp[:] += len(dataTemp) - fp = open(voidDir+"/z_void_indexes.txt", "a") - for idx in idxTemp: - fp.write(str(idx)+"\n") - fp.close() - dataTemp = () - - fp = NetCDFFile(voidDir+"/z_posx.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/z_posx.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/z_posx.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/z_posy.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/z_posy.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/z_posy.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/z_posz.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/z_posz.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/z_posz.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - idxTemp = file(sourceStackDir+"/void_indexes.txt", "r").\ - readlines() - idxTemp = np.array(idxTemp, dtype='i') - fp = NetCDFFile(voidDir+"/posx.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - idxTemp[:] += len(dataTemp) - fp = open(voidDir+"/void_indexes.txt", "a") - for idx in idxTemp: - fp.write(str(idx)+"\n") - fp.close() - dataTemp = () - - fp = NetCDFFile(voidDir+"/posx.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/posx.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/posx.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/posy.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/posy.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/posy.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/posz.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/posz.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/posz.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - - fp = NetCDFFile(voidDir+"/redshifts.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/redshifts.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/redshifts.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - fp = NetCDFFile(voidDir+"/indexes.nc") - dataTemp = fp.variables['array'][0:] - fp.close() - fp = NetCDFFile(sourceStackDir+"/indexes.nc") - dataTemp2 = fp.variables['array'][0:] - fp.close() - dataTemp = np.append(dataTemp, dataTemp2) - outFile = NetCDFFile(voidDir+"/indexes.nc", mode='w') - outFile.createDimension("dim", len(dataTemp)) - v = outFile.createVariable("array", ncFloat, ("dim",)) - v[:] = dataTemp - outFile.close() - - # add NOVOID if necessary - if not os.access(voidDir+"/centers.txt", os.F_OK): - fp = open(voidDir+"/NOVOID", "w") - fp.write("no voids in combo\n") - fp.close() - - fp = open(zobovDir+"/donegalupdate", "w") - fp.close() - - print "Done!" - - sys.stdout = sys.__stdout__ - #sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print jobString, "Combining stacks done" - else: - print jobString, "Combining stacks FAILED!" - exit(-1) - - #else: - # print "already done!" - -# ----------------------------------------------------------------------------- -def launchProfile(sample, stack, voidDir=None, logFile=None, continueRun=None, - thisDataPortion=None): - - sampleName = sample.fullName - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion) - jobString = " "+runSuffix+":" - - if os.access(voidDir+"/NOVOID", os.F_OK): - print jobString, "Profile no stack here; skipping!" - return - - numVoids = open(voidDir+"/num_voids.txt", "r").readline() - numParticles = open(voidDir+"/num_particles.txt", "r").readline() - - #Rextracut = stack.rMin*3 + 1 - #Rcircular = stack.rMin*2 + 2 - Rextracut = stack.rMax*3 + 1 - Rcircular = stack.rMax*3 + 2 - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - profileParmFile = voidDir+"/params.txt" - fp = open(profileParmFile, "w") - fp.write(str(stack.rMin)+"\n") - fp.write(str(stack.rMax)+"\n") - fp.write(str(Rcircular)+"\n") - fp.write(numVoids+"\n") - fp.close() - - #if sample.zRange[0] > stack.zMax or sample.zRange[1] < stack.zMin: - # print "outside sample redshift range; skipping!" - # fp = open(voidDir+"/NOVOID", "w") - # fp.write("outside z range: profile\n") - # fp.close() - # return - - if sample.profileBinSize == "auto": - #density = 0.5 * 50 / Rcircular / 2 - density = stack.rMin / 5. - #density = stack.rMax / 5. - else: - density = sample.profileBinSize - - sys.stdout = open(logFile, 'w') - #sys.stderr = open(logFile, 'a') - vp.build_1d_profile(base_dir=voidDir, density=density, - rescaleMode=stack.rescaleMode) - vp.build_2d_profile(base_dir=voidDir, density=density) - sys.stdout = sys.__stdout__ - #sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print jobString, "Profiling stacks done (N_v=%s, N_p=%s)" % \ - (numVoids,numParticles) - else: - print jobString, "Profiling stacks FAILED!" - exit(-1) - - else: - print jobString, "Profiling stacks already done!" - - -# ----------------------------------------------------------------------------- -def launchFit(sample, stack, useComoving=None, - logFile=None, voidDir=None, figDir=None, - continueRun=None, thisDataPortion=None, - fittingMode=None, minVoidsToFit=None, - minPartToFit=None): - - sampleName = sample.fullName - - runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - stack.rMax, thisDataPortion) - - jobString = " "+runSuffix+":" - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - if os.access(voidDir+"/NOVOID", os.F_OK): - print jobString, "Fitting no voids here; skipping!" - return - - numVoids = int(open(voidDir+"/num_voids.txt", "r").readline()) - if numVoids < minVoidsToFit: - print jobString, "Fitting not enough voids to fit; skipping!" - fp = open(voidDir+"/NOFIT", "w") - fp.write("not enough voids: %d \n" % numVoids) - fp.close() - return - - numPart = int(open(voidDir+"/num_particles.txt", "r").readline()) - if numPart < minPartToFit: - print jobString, "Fitting not enough particles to fit; skipping!" - fp = open(voidDir+"/NOFIT", "w") - fp.write("not enough particles: %d \n" % numPart) - fp.close() - return - - #if stack.zMin < sample.zRange[0] or stack.zMax > sample.zRange[1]: - # print "outside sample redshift range; skipping!" - # return - - if sample.partOfCombo or not sample.includeInHubble: - print jobString, "Fitting sample not needed for further analysis; skipping!" - fp = open(voidDir+"/NOFIT", "w") - fp.write("sample not needed for hubble\n") - fp.close() - return - - if not stack.includeInHubble: - print jobString, "Fitting radius not needed for further analysis; skipping!" - return - - if not stack.includeInHubble: - print jobString, "Fitting redshift not needed for further analysis; skipping!" - return - - if os.access(figDir+"/stackedVoid_"+sampleName+"_"+\ - runSuffix+".eps", os.F_OK): - os.unlink(figDir+"/stackedVoid_"+sampleName+"_"+\ - runSuffix+".eps") - - sys.stdout = open(logFile, 'w') - #sys.stderr = open(logFile, 'a') - - badChain = True - ntries = 0 - maxtries = 5 - while badChain: - ntries += 1 - - if fittingMode == "mcmc": - Rexpect = (stack.rMin+stack.rMax)/2 - #Rtruncate = stack.rMax*2. # TEST - Rtruncate = stack.rMax*3. # TEST - #Rtruncate = stack.rMin*3. + 1 # TEST - if sample.dataType == "observation": - ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, - nbar_init=1.0, - Niter=300000, - Nburn=100000, - Rextracut=Rtruncate) - else: - ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, - nbar_init=1.0, - Niter=300000, - Nburn=100000, - Rextracut=Rtruncate) - #print "TEST", stack.rMax, args[0][0], args[0][1], args[0][2] - badChain = (args[0][0] > 0.5 or args[0][1] > 2.*stack.rMax or \ - args[1] > 2.0 or \ - args[0][2] > 2.*stack.rMax) and \ - (ntries < maxtries) - elif fittingMode == "inertia": - ret,fits,args = vp.compute_inertia(voidDir, stack.rMax, mode="2d", nBootstraps=100, rMaxInertia=0.7) - badChain = False - - ## rescale to expected stretch if using comoving coords - #if useComoving or not sample.useLightCone: - # exp = args[1] - # expError = args[2] - - #np.save(voidDir+"/chain.npy", ret) - np.savetxt(voidDir+"/fits.out", fits) - - plotTitle = "Sample: "+sample.nickName+\ - ", z = "+str(stack.zMin)+"-"+\ - str(stack.zMax)+", R = "+str(stack.rMin)+"-"+\ - str(stack.rMax)+r" $h^{-1}$ Mpc" - - showRatio = (ntries <= maxtries) - showContours = (ntries <= maxtries) - - if stack.rescaleMode == "rmax": - rescaleFactor = stack.rMax - else: - rescaleFactor = 1.0 - - vp.draw_fit(*args, delta_tick=0.2, vmax=1.4, delta_v=0.01, - nocontour=True, model_max=1.0, delta_model=0.1, - plotTitle=plotTitle, show_ratio=showRatio, - showContours=showContours, - rescaleFactor=rescaleFactor) - figure(1).savefig(figDir+"/stackedVoid_"+sampleName+\ - "_"+runSuffix+".eps") - figure(1).savefig(figDir+"/stackedVoid_"+sampleName+\ - "_"+runSuffix+".pdf") - figure(1).savefig(figDir+"/stackedVoid_"+sampleName+\ - "_"+runSuffix+".png") - # record the measured stretch - exp = args[1] - expError = args[2] - outline = str(exp) + " " + str(expError) - open(voidDir+"/expansion.out", "w").write(outline) - - print "Done!" - sys.stdout = sys.__stdout__ - #sys.stderr = sys.__stderr__ - - if os.access(voidDir+"/NOFIT", os.F_OK): - os.unlink(voidDir+"/NOFIT") - if ntries > maxtries: - print jobString, "No reliable fit found; skipping!" - fp = open(voidDir+"/NOFIT", "w") - fp.write("bad ellipticity fit\n") - fp.close() - #return - - else: - if jobSuccessful(logFile, "Done!\n"): - print jobString, "Fitting done (%d tries)" % ntries - else: - print jobString, "Fitting FAILED!" - exit(-1) - - else: - print jobString, "already done!" - - -# ----------------------------------------------------------------------------- -def launchHubble(dataPortions=None, dataSampleList=None, logDir=None, - INCOHERENT=None, workDir=None, figDir=None, errorBars=None, - ZOBOV_PATH=None, continueRun=None, voidDir=None, - doPlot = True, setName=None): - - for thisDataPortion in dataPortions: - print " For data portion", thisDataPortion - sys.stdout.flush() - - # merge the bins from all samples - numSamplesForHubble = 0 - maxRBins = 0 - maxZBins = 0 - allRBins = [] - allZBins = [] - - for sample in dataSampleList: - if sample.includeInHubble: - numSamplesForHubble += 1 - - for stack in sample.getHubbleStacks(): - alreadyHere = False - for stackCheck in allRBins: - if stack.rMin == stackCheck.rMin and stack.rMax==stackCheck.rMax: - alreadyHere = True - break - if not alreadyHere: - allRBins.append(stack) - - alreadyHere = False - for stackCheck in allZBins: - if stack.zMin == stackCheck.zMin and stack.zMax==stackCheck.zMax: - alreadyHere = True - break - if not alreadyHere: - allZBins.append(stack) - - allExpList = np.zeros((numSamplesForHubble,len(allRBins),len(allZBins),5)) - allExpList[:] = -1 - aveDistList = np.zeros((len(allZBins),3)) - - iSample = -1 - for sample in dataSampleList: - if not sample.includeInHubble: - continue - - iSample += 1 - - zobovDir = sample.zobovDir - sampleName = sample.fullName - - print " For data set", sampleName, "...", - sys.stdout.flush() - - logFile = logDir+"/hubble_"+sampleName+"_"+thisDataPortion+".out" - - #AVEDIST_PATH = ZOBOV_PATH+"/mytools/computeAverageDistortionPMS" - - # compute the average stretch in each redshift bin for cleaner-looking - # (not not fully accurate) plots - for (iZBin, stack) in enumerate(sample.getUniqueZBins()): - - aveDist = vp.aveStretch(sample, stack.zMin, stack.zMax, - Om=sample.omegaM) - - aveDistList[iZBin, 0] = stack.zMin - aveDistList[iZBin, 1] = aveDist - aveDistList[iZBin, 2] = 0.00125 - - numFits = 0 - fp = file(zobovDir+'/calculatedExpansions_'+thisDataPortion+'.txt', - mode="w") - - expList = np.zeros((1, len(sample.getUniqueRBins()), - len(sample.getUniqueZBins()), 3)) - - for (iZBin, zBin) in enumerate(sample.getUniqueZBins()): - for (iR, rBin) in enumerate(sample.getUniqueRBins()): - - runSuffix = getStackSuffix(zBin.zMin, zBin.zMax, rBin.rMin, - rBin.rMax, thisDataPortion) - - expFile = zobovDir+"/stacks_"+runSuffix+"/expansion.out" - if not (os.access(zobovDir+"/stacks_"+runSuffix+"/NOVOID", \ - os.F_OK) or - os.access(zobovDir+"/stacks_"+runSuffix+"/NOFIT", \ - os.F_OK) or - not os.access(expFile, os.F_OK)): - exp = float(open(expFile, "r").readline().split()[0]) - expError = float(open(expFile, "r").readline().split()[1]) - else: - exp = -1 - expError = 0 - - expList[0, iR, iZBin, 0] = exp - expList[0, iR, iZBin, 1] = expError - - # TODO - # compute the per-stack stretch for liklihood calculation - #runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, - # stack.rMax, thisDataPortion) - #voidDir = zobovDir+"/stacks_" + runSuffix - #centersFile = voidDir+"/centers.txt" - #voidRedshifts = np.loadtxt(centersFile)[:,5] - #aveDist = vp.aveWeightedStretchCone(stack.zMin, stack.zMax, - # skyFrac = sample.skyFraction, - # voidRedshifts=voidRedshifts) - - aveDist = vp.aveStretch(sample, zBin.zMin, zBin.zMax, - Om=sample.omegaM) - - expList[0, iR, iZBin, 2] = aveDist - - voidDir = zobovDir+"/stacks_" + runSuffix - numVoids = open(voidDir+"/num_voids.txt", "r").readline() - numParticles = open(voidDir+"/num_particles.txt", "r").readline() - - for (iZCheck,zBinCheck) in enumerate(allZBins): - for (iRCheck,rBinCheck) in enumerate(allRBins): - if zBin.zMin == zBinCheck.zMin and zBin.zMax == zBinCheck.zMax: - if rBin.rMin == rBinCheck.rMin and rBin.rMax == rBinCheck.rMax: - aveDist = vp.aveStretch(sample, zBin.zMin, zBin.zMax, - Om=sample.omegaM) - - allExpList[iSample, iRCheck, iZCheck, 0] = exp - allExpList[iSample, iRCheck, iZCheck, 1] = expError - allExpList[iSample, iRCheck, iZCheck, 3] = numVoids - allExpList[iSample, iRCheck, iZCheck, 4] = numParticles - - fp.write(str(exp) + " " + str(expError) + " "+ " "+ \ - str(aveDist) + " -1 " + runSuffix+"\n") - numFits += 1 - - fp.close() - - rlist = np.zeros((len(sample.getUniqueRBins()),2)) - for (iR,rBin) in enumerate(sample.getUniqueRBins()): - rlist[iR, 0] = rBin.rMin - rlist[iR, 1] = rBin.rMax - - zbase = np.zeros((len(sample.getUniqueZBins()),2)) - for (iZ,zBin) in enumerate(sample.getUniqueZBins()): - zbase[iZBin,0] = zBin.zMinPlot - zbase[iZBin,1] = zBin.zMaxPlot - #np.savetxt(zobovDir+'/zbase_'+thisDataPortion+'.txt', zbase) - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - sys.stdout = open(logFile, 'w') - #sys.stderr = open(logFile, 'a') - plotTitle = "Sample: "+sample.nickName+", "+thisDataPortion+" voids" - #if doPlot: - #vp.do_all_obs(sample, zbase, expList, workDir+"/avedistortion_", - #vp.do_all_obs(sample, zbase, expList, aveDistList, - # rlist, plotTitle=plotTitle, plotAve=True) - #figure(1).savefig(figDir+"/hubble_"+setName+"_"+sampleName+"_"+\ - # thisDataPortion+\ - # ".eps",bbox_inches='tight') - #figure(1).savefig(figDir+"/hubble_"+setName+"_"+sampleName+"_"+\ - # thisDataPortion+\ - # ".pdf",bbox_inches='tight') - #figure(1).savefig(figDir+"/hubble_"+setName+"_"+sampleName+"_"+\ - # thisDataPortion+\ - # ".png",bbox_inches='tight') - #else: - print "Skipping plot" - print "Done!" - sys.stdout = sys.__stdout__ - #sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - # now do all samples in one plot - print " For data set combined...", - sys.stdout.flush() - - logFile = logDir+"/hubble_combined_"+sampleName+"_"+thisDataPortion+".out" - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - - if errorBars == "ESTIMATED": - # replace calculated errors with estimated ones - errorBarFile = workDir+"/calculatedErrorBars_"+thisDataPortion+".txt" - if os.access(errorBarFile, os.F_OK): - print "REPLACING ERROR BARS!", - sys.stdout.flush() - fp = file(errorBarFile, mode="r") - ignoreLine = fp.readline() - for iZ in xrange(len(allExpList[0,0,:,0])): - ignoreLine = fp.readline() - - for iZ in xrange(len(allExpList[0,0,:,0])): - for iSample in xrange(len(allExpList[:,0,0,0])): - for iR in xrange(len(allExpList[0,:,0,0])): - line = fp.readline().split() - allExpList[iSample,iR,iZ,2] = allExpList[iSample,iR,iZ,1] - allExpList[iSample,iR,iZ,1] = float(line[1]) - allExpList[iSample,iR,iZ,3] = float(line[2]) - allExpList[iSample,iR,iZ,4] = float(line[3]) - - fp.close() - - rlist = np.zeros((len(allRBins),2)) - for (iR,rBin) in enumerate(allRBins): - rlist[iR, 0] = rBin.rMin - rlist[iR, 1] = rBin.rMax - - zbase = np.zeros((len(allZBins),2)) - plotZmax = 0.0 - plotZmin = 1.e99 - for (iZ,zBin) in enumerate(allZBins): - zbase[iZ,0] = zBin.zMinPlot - zbase[iZ,1] = zBin.zMaxPlot - - if zBin.zMaxPlot > plotZmax: plotZmax = zBin.zMaxPlot - if zBin.zMinPlot < plotZmin: plotZmin = zBin.zMinPlot - - aveDist = vp.aveStretch(sample, zBin.zMin, zBin.zMax, - Om=sample.omegaM) - - aveDistList[iZ, 0] = zBin.zMin - aveDistList[iZ, 1] = aveDist - aveDistList[iZ, 2] = 0.00125 - if plotZmax > 1.5: plotZmax = 1.5 - - - shortSampleNames = list() - for sample in dataSampleList: - if sample.includeInHubble: - shortSampleNames.append(sample.nickName) - sys.stdout = open(logFile, 'w') - #sys.stderr = open(logFile, 'a') - if doPlot: - if INCOHERENT: - #plotTitle = "all samples, incoherent "+\ - # thisDataPortion+" voids" - plotTitle = '' - else: - #plotTitle = "all samples, "+thisDataPortion+" voids" - plotTitle = setName - #plotTitle = '' - vp.do_all_obs(sample, zbase, allExpList, aveDistList, - rlist, plotTitle=plotTitle, sampleNames=shortSampleNames, - plotAve=True, mulfac = 1.0, biasLine = 1.16, - plotZmin=plotZmin, plotZmax=plotZmax, - plotOm1=True) - figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+ \ - thisDataPortion+\ - ".eps",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+ \ - thisDataPortion+\ - ".pdf",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+ \ - thisDataPortion+\ - ".png",bbox_inches='tight') - - if INCOHERENT: - plotTitle = "all samples, incoherent "+\ - thisDataPortion+" voids (systematics corrected)" - else: - #plotTitle = "all samples, "+thisDataPortion+\ - # " voids (systematics corrected)" - plotTitle = setName + "(systematics corrected)" - vp.do_all_obs(sample, zbase, allExpList, aveDistList, - rlist, plotTitle=plotTitle, sampleNames=shortSampleNames, - plotAve=True, mulfac = 1.16, - plotZmin=plotZmin, plotZmax=plotZmax, - plotOm1=True) - figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+\ - thisDataPortion+\ - "_debiased.eps",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+\ - thisDataPortion+\ - "_debiased.pdf",bbox_inches='tight') - figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+\ - thisDataPortion+\ - "_debiased.png",bbox_inches='tight') - else: - print "Skipping plot" - print "Done!" - sys.stdout = sys.__stdout__ - #sys.stderr = sys.__stderr__ - - # save all expansion data to a single file - fp = file(workDir+'/calculatedExpansions_'+thisDataPortion+'.txt', - mode="w") - fp.write(str(numSamplesForHubble) + " " + - str(len(allRBins)) + " " + str(len(allZBins)) + "\n") - for zBin in allZBins: - fp.write(str(zBin.zMin) + " " + str(zBin.zMax) + "\n") - - for iZ in xrange(len(allExpList[0,0,:,0])): - for iSample in xrange(len(allExpList[:,0,0,0])): - for iR in xrange(len(allExpList[0,:,0,0])): - fp.write(str(allExpList[iSample,iR,iZ,0]) + " " +\ - str(allExpList[iSample,iR,iZ,1]) + " " +\ - str(allExpList[iSample,iR,iZ,2]) + " " +\ - str(allExpList[iSample,iR,iZ,3]) + " " +\ - str(allExpList[iSample,iR,iZ,4]) + "\n") - fp.close() - - fp = file(workDir+'/calculatedExpansions_reduced_'+thisDataPortion+'.txt', - mode="w") - for iZ in xrange(len(allExpList[0,0,:,0])): - for iSample in xrange(len(allExpList[:,0,0,0])): - for iR in xrange(len(allExpList[0,:,0,0])): - if (allExpList[iSample,iR,iZ,0] == -1): continue - fp.write(str(allExpList[iSample,iR,iZ,0]) + " " +\ - str(allExpList[iSample,iR,iZ,1]) + " " +\ - str(allExpList[iSample,iR,iZ,2]) + " " +\ - str(allExpList[iSample,iR,iZ,3]) + " " +\ - str(allExpList[iSample,iR,iZ,4]) + "\n") - fp.close() - - - # save all void distribution data to a single file - fp = file(workDir+'/voidDistributions_'+thisDataPortion+'.txt', - mode="w") - fp.write(str(numSamplesForHubble) + " " + - str(len(allRBins)) + " " + str(len(allZBins)) + "\n") - for zBin in allZBins: - fp.write(str(zBin.zMin) + " " + str(zBin.zMax) + "\n") - - for zBin in allZBins: - for sample in dataSampleList: - if not sample.includeInHubble: continue - for rBin in allRBins: - runSuffix = getStackSuffix(zBin.zMin, zBin.zMax, rBin.rMin, - rBin.rMax, thisDataPortion) - voidDir = sample.zobovDir+"/stacks_" + runSuffix - centersFile = voidDir+"/centers.txt" - if os.access(centersFile, os.F_OK): - voidRedshifts = np.loadtxt(centersFile) - if voidRedshifts.ndim > 1: - voidRedshifts = voidRedshifts[:,5] - np.savetxt(fp, voidRedshifts[None]) - else: - if (len(voidRedshifts) > 0): - voidRedshifts = voidRedshifts[5] - np.savetxt(fp, voidRedshifts[None]) - else: - fp.write("-1\n") - #fp.write(str(len(voidRedshifts))+" ") - #np.savetxt(fp, voidRedshifts[None]) - else: - fp.write("-1\n") - fp.close() - - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - -# ----------------------------------------------------------------------------- -def launchLikelihood(dataSampleList=None, - dataPortions=None, logDir=None, workDir=None, - continueRun=False): - - # check for comoving from first sample - useComoving = not dataSampleList[0].useLightCone or \ - dataSampleList[0].useComoving - - for thisDataPortion in dataPortions: - print " For data portion", thisDataPortion, "...", - sys.stdout.flush() - - logFile = logDir+"/likelihood_"+thisDataPortion+".out" - - if not (continueRun and jobSuccessful(logFile, "Done!\n")): - - #sys.stdout = open(logFile, 'w') - #sys.stderr = open(logFile, 'a') - - if dataSampleList[0].usePecVel: - biasLevel = 1.16 - else: - biasLevel = 1.0 - - vp.build1dLikelihood(dataSampleList[0], - workDir+"/calculatedExpansions_"+\ - thisDataPortion+".txt", - workDir+"/voidDistributions_" + \ - thisDataPortion+".txt", - nbins = 20, - OmStart = 0.0, - OmEnd = 1.0, - #biasStart = 1.0, - #biasEnd = 1.32, - biasStart = biasLevel, - biasEnd = biasLevel, - outputBase = workDir+"/1dlikelihoods_"+thisDataPortion+"_", - useBinAve = False, - useComoving = useComoving, - OmFid=dataSampleList[0].omegaM) - - #sys.stdout = sys.__stdout__ - #sys.stderr = sys.__stderr__ - - if jobSuccessful(logFile, "Done!\n"): - print "done" - else: - print "FAILED!" - exit(-1) - - else: - print "already done!" - - -# ----------------------------------------------------------------------------- -def launchEstimateErrorBars(workDir=None, nullDir=None, numNulls=None, - dataPortion=None, nullDirList=None): - - IVALUE = 0 - IERROR = 1 - - # open first null file to get sizes - if nullDirList == None: - nullFile = nullDir + "/" + str(1) + "/calculatedExpansions_" + \ - str(dataPortion) + ".txt" - else: - nullFile = nullDirList[0] + "/calculatedExpansions_" + str(dataPortion) + \ - ".txt" - numNulls = len(nullDirList) - - fp = open(nullFile, "r") - binLine = fp.readline() - numSamples = int(binLine.split(" ")[0]) - numRBins = int(binLine.split(" ")[1]) - numZBins = int(binLine.split(" ")[2]) - numBins = numRBins * numSamples * numZBins - nRTotal = numRBins * numSamples - fp.close() - - # load null data - nullData = np.zeros([numNulls, numBins, 2]) - - for iNull in xrange(numNulls): - if nullDirList == None: - nullFile = nullDir + "/" + str(iNull) + "/calculatedExpansions_" + \ - str(dataPortion) + ".txt" - else: - nullFile = nullDirList[iNull] + "/calculatedExpansions_" + \ - str(dataPortion) + ".txt" - - fp = open(nullFile, "r") - - binLine = fp.readline() - numSamples = int(binLine.split(" ")[0]) - numRBins = int(binLine.split(" ")[1]) - numZBins = int(binLine.split(" ")[2]) - numBins = numRBins * numSamples * numZBins - - zList = np.zeros( [numZBins,2] ) - for iZ in xrange(int(numZBins)): - line = fp.readline().split() - zList[iZ,0] = float(line[0]) - zList[iZ,1] = float(line[1]) - - for iZ in xrange(numZBins): - iRTotal = 0 - for iSample in xrange(numSamples): - for iR in xrange(numRBins): - line = fp.readline().split() - nullData[iNull,iRTotal*numZBins+iZ,IVALUE] = float(line[0]) - nullData[iNull,iRTotal*numZBins+iZ,IERROR] = float(line[1]) - iRTotal += 1 - fp.close() - - # use 68% limits to get error bars - errorBars = np.zeros([numBins, 2]) - for iBin in xrange(numBins): - binData = nullData[:,iBin,IVALUE] - validData = binData[:] > -1 - binData = binData[validData] - binData = np.sort(binData) - - if np.size(binData) > 0: - errorBars[iBin, 0] = binData[int(0.16*np.size(binData))] - errorBars[iBin, 1] = binData[int(0.84*np.size(binData))] - - mean = (errorBars[iBin, 1] + errorBars[iBin, 0])/2 - sig = (errorBars[iBin, 1] - errorBars[iBin, 0])/2 - errorBars[iBin, 0] = mean - errorBars[iBin, 1] = sig - else: - errorBars[iBin, :] = -1 - - # writes out error bars - fp = file(workDir+"/calculatedErrorBars_" + \ - str(dataPortion) + ".txt", mode="w") - fp.write(str(numSamples) + " " + str(numRBins) + " " + str(numZBins) + "\n") - for iZ in xrange(numZBins): - fp.write(str(zList[iZ,0]) + " " + str(zList[iZ,1]) + "\n") - - for iZ in xrange(numZBins): - iRTotal = 0 - for iSample in xrange(numSamples): - for iR in xrange(numRBins): - fp.write(str(errorBars[iRTotal*numZBins+iZ,0]) + " " + \ - str(errorBars[iRTotal*numZBins+iZ,1]) + "\n") - iRTotal += 1 - fp.close() - - - - # ----------------------------------------------------------------------------- def launchVelocityStack(sample, stack, binPath, velField_file, diff --git a/python_tools/void_python_tools/partUtil/partUtil.py b/python_tools/void_python_tools/partUtil/partUtil.py index 37309ab..4f695e1 100644 --- a/python_tools/void_python_tools/partUtil/partUtil.py +++ b/python_tools/void_python_tools/partUtil/partUtil.py @@ -34,7 +34,7 @@ ncFloat = 'f8' # ----------------------------------------------------------------------------- def loadPart(sampleDir): - #print " Loading particle data..." + print " Loading particle data..." sys.stdout.flush() with open(sampleDir+"/sample_info.dat", 'rb') as input: diff --git a/python_tools/void_python_tools/xcor/xcorlib.py b/python_tools/void_python_tools/xcor/xcorlib.py index d70ebad..535c90d 100644 --- a/python_tools/void_python_tools/xcor/xcorlib.py +++ b/python_tools/void_python_tools/xcor/xcorlib.py @@ -2,72 +2,12 @@ import numpy as np def cic( x, Lbox, Lboxcut = 0, Nmesh = 128, weights = None ): - if weights == None: weights = 1 - wm = np.mean(weights) - ws = np.mean(weights**2) - - d = np.mod(x/(Lbox+2*Lboxcut)*Nmesh,1) - - box = ([Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut],[Lboxcut,Lbox+Lboxcut]) - - rho = np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*(1-d[:,2]))[0] \ - + np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*(1-d[:,2]))[0],1,0) \ - + np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*(1-d[:,2]))[0],1,1) \ - + np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*(1-d[:,1])*d[:,2])[0],1,2) \ - + np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*(1-d[:,2]))[0],1,0),1,1) \ - + np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*(1-d[:,1])*d[:,2])[0],1,0),1,2) \ - + np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*(1-d[:,0])*d[:,1]*d[:,2])[0],1,1),1,2) \ - + np.roll(np.roll(np.roll(np.histogramdd(x, range = box, bins = Nmesh, weights = weights*d[:,0]*d[:,1]*d[:,2])[0],1,0),1,1),1,2) - - rho /= wm - - rho = rho/rho.mean() - 1. - - return (rho, wm, ws) - + return def powcor( d1, d2, Lbox, Nmesh = 128, Nbin = 100, scale = 'lin', cor = False ): - # Fourier transform - d1 = np.fft.fftn(d1) - d2 = np.fft.fftn(d2) - - # CIC correction - wid = np.indices(np.shape(d1)) - Nmesh/2 - #wid[np.where(wid >= Nmesh/2)] -= Nmesh - wid = wid*np.pi/Nmesh + 1e-100 - wcic = np.prod(np.sin(wid)/wid,0)**2 - - # Shell average power spectrum - dk = 2*np.pi/Lbox - Pk = np.conj(d1)*d2*(Lbox/Nmesh**2)**3 - Pk = np.fft.fftshift(Pk)/wcic**2 - - (Nm, km, Pkm, SPkm) = shellavg(np.real(Pk), dk, Nmesh, Nbin = Nbin, xmin = 0., xmax = Nmesh*dk/2, scale = scale) - - # Inverse Fourier transform and shell average correlation function - if cor == True: - dx = Lbox/Nmesh - Xr = np.fft.ifftshift(np.fft.ifftn(np.fft.ifftshift(Pk)))*(Nmesh/Lbox)**3 - - (Nmx, rm, Xrm, SXrm) = shellavg(np.real(Xr), dx, Nmesh, Nbin = Nbin, xmin = dx, xmax = 140., scale = scale) - - return ((Nm, km, Pkm, SPkm),(Nmx, rm, Xrm, SXrm)) - - else: return (Nm, km, Pkm, SPkm) - + return def shellavg( f, dx, Nmesh, Nbin = 100, xmin = 0., xmax = 1., scale = 'lin' ): - - x = np.indices(np.shape(f)) - Nmesh/2 - #x[np.where(x >= Nmesh/2)] -= Nmesh - x = dx*np.sqrt(np.sum(x**2,0)) - if scale == 'lin': bins = xmin+(xmax-xmin)* (np.arange(Nbin+1)/float(Nbin)) - if scale == 'log': bins = xmin*(xmax/xmin)**(np.arange(Nbin+1)/float(Nbin)) - Nm = np.histogram(x, bins = bins)[0] - xm = np.histogram(x, bins = bins, weights = x)[0]/Nm - fm = np.histogram(x, bins = bins, weights = f)[0]/Nm - fs = np.sqrt((np.histogram(x, bins = bins, weights = f**2)[0]/Nm - fm**2)/(Nm-1)) - - return (Nm, xm, fm, fs) + return diff --git a/zobov/CMakeLists.txt b/zobov/CMakeLists.txt new file mode 100644 index 0000000..9c749b5 --- /dev/null +++ b/zobov/CMakeLists.txt @@ -0,0 +1,16 @@ + +add_executable(voz1b1 voz1b1.c readfiles.c vozutil.c voz.h) +target_link_libraries(voz1b1 ${QHULL_LIBRARY} ${MATH_LIB}) + + +add_executable(jozov jozov.c findrtop.c) +target_link_libraries(jozov ${MATH_LIB}) + +add_executable(vozinit vozinit.c readfiles.c) +target_link_libraries(vozinit ${MATH_LIB}) + + +add_executable(voztie voztie.c readfiles.c) +target_link_libraries(voztie ${MATH_LIB}) + +