From 4a84fc3791994d5c9cc2829043175c22f7e7c5d8 Mon Sep 17 00:00:00 2001 From: Josh Siegle Date: Wed, 24 Apr 2024 11:35:16 -0700 Subject: [PATCH] Initial implementation --- README.md | 22 +- doc_template/source/conf.py | 7 +- pyproject.toml | 5 + src/aind_ephys_rig_qc/__init__.py | 1 + src/aind_ephys_rig_qc/generate_report.py | 257 +++++++++++++++++ src/aind_ephys_rig_qc/images/aind-logo.png | Bin 0 -> 40774 bytes src/aind_ephys_rig_qc/parameters.json | 5 + src/aind_ephys_rig_qc/pdf_utils.py | 106 +++++++ src/aind_ephys_rig_qc/qc_figures.py | 71 +++++ src/aind_ephys_rig_qc/temporal_alignment.py | 288 ++++++++++++++++++++ tests/test_example.py | 16 -- tests/test_generate_report.py | 20 ++ 12 files changed, 769 insertions(+), 29 deletions(-) create mode 100644 src/aind_ephys_rig_qc/generate_report.py create mode 100644 src/aind_ephys_rig_qc/images/aind-logo.png create mode 100644 src/aind_ephys_rig_qc/parameters.json create mode 100644 src/aind_ephys_rig_qc/pdf_utils.py create mode 100644 src/aind_ephys_rig_qc/qc_figures.py create mode 100644 src/aind_ephys_rig_qc/temporal_alignment.py delete mode 100644 tests/test_example.py create mode 100644 tests/test_generate_report.py diff --git a/README.md b/README.md index 7ff36f6..c202236 100644 --- a/README.md +++ b/README.md @@ -8,13 +8,11 @@ ![Python](https://img.shields.io/badge/python->=3.7-blue?logo=python) +## Overview -## Usage - - To use this template, click the green `Use this template` button and `Create new repository`. - - After github initially creates the new repository, please wait an extra minute for the initialization scripts to finish organizing the repo. - - To enable the automatic semantic version increments: in the repository go to `Settings` and `Collaborators and teams`. Click the green `Add people` button. Add `svc-aindscicomp` as an admin. Modify the file in `.github/workflows/tag_and_publish.yml` and remove the if statement in line 10. The semantic version will now be incremented every time a code is committed into the main branch. - - To publish to PyPI, enable semantic versioning and uncomment the publish block in `.github/workflows/tag_and_publish.yml`. The code will now be published to PyPI every time the code is committed into the main branch. - - The `.github/workflows/test_and_lint.yml` file will run automated tests and style checks every time a Pull Request is opened. If the checks are undesired, the `test_and_lint.yml` can be deleted. The strictness of the code coverage level, etc., can be modified by altering the configurations in the `pyproject.toml` file and the `.flake8` file. +This repository contains scripts that run after the completion of extracellular electrophysiology experiments. It serves two main purposes: +- Ensuring that all data streams are aligned to a common timebase prior to upload +- Generating a QC report that can be used to validate ephys data quality ## Installation To use the software, in the root directory, run @@ -42,12 +40,12 @@ coverage run -m unittest discover && coverage report - Use **interrogate** to check that modules, methods, etc. have been documented thoroughly: ```bash -interrogate . +interrogate -v . ``` -- Use **flake8** to check that code is up to standards (no unused imports, etc.): +- Use **isort** to automatically sort import statements: ```bash -flake8 . +isort . ``` - Use **black** to automatically format the code into PEP standards: @@ -55,11 +53,13 @@ flake8 . black . ``` -- Use **isort** to automatically sort import statements: +- Use **flake8** to check that code is up to standards (no unused imports, etc.): ```bash -isort . +flake8 . ``` + + ### Pull requests For internal members, please create a branch. For external members, please fork the repository and open a pull request from the fork. We'll primarily use [Angular](https://github.com/angular/angular/blob/main/CONTRIBUTING.md#commit) style for commit messages. Roughly, they should follow the pattern: diff --git a/doc_template/source/conf.py b/doc_template/source/conf.py index afc260a..6d52a0d 100644 --- a/doc_template/source/conf.py +++ b/doc_template/source/conf.py @@ -1,12 +1,15 @@ """Configuration file for the Sphinx documentation builder.""" + # # For the full list of built-in configuration values, see the documentation: # https://www.sphinx-doc.org/en/master/usage/configuration.html +from datetime import date + # -- Path Setup -------------------------------------------------------------- -from os.path import dirname, abspath +from os.path import abspath, dirname from pathlib import Path -from datetime import date + from aind_ephys_rig_qc import __version__ as package_version INSTITUTE_NAME = "Allen Institute for Neural Dynamics" diff --git a/pyproject.toml b/pyproject.toml index 73c3124..b460011 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,11 @@ readme = "README.md" dynamic = ["version"] dependencies = [ + 'open-ephys-python-tools', + 'matplotlib', + 'fpdf2', + 'scipy', + 'rich' ] [project.optional-dependencies] diff --git a/src/aind_ephys_rig_qc/__init__.py b/src/aind_ephys_rig_qc/__init__.py index d0a8547..024d51f 100644 --- a/src/aind_ephys_rig_qc/__init__.py +++ b/src/aind_ephys_rig_qc/__init__.py @@ -1,2 +1,3 @@ """Init package""" + __version__ = "0.0.0" diff --git a/src/aind_ephys_rig_qc/generate_report.py b/src/aind_ephys_rig_qc/generate_report.py new file mode 100644 index 0000000..7c6c7d2 --- /dev/null +++ b/src/aind_ephys_rig_qc/generate_report.py @@ -0,0 +1,257 @@ +""" +Generates a PDF report from an Open Ephys data directory +""" + +import json +import os +import sys + +import numpy as np +import pandas as pd +from open_ephys.analysis import Session + +from aind_ephys_rig_qc import __version__ as package_version +from aind_ephys_rig_qc.pdf_utils import PdfReport +from aind_ephys_rig_qc.qc_figures import plot_power_spectrum, plot_raw_data + + +def generate_qc_report(directory, report_name="QC.pdf"): + """ + Generates a PDF report from an Open Ephys data directory + + Saves QC.pdf + + Parameters + ---------- + directory : str + The path to the Open Ephys data directory + report_name : str + The name of the PDF report + + """ + + pdf = PdfReport("aind-ephys-rig-qc v" + package_version) + pdf.add_page() + pdf.set_font("Helvetica", "B", size=12) + pdf.set_y(30) + pdf.write(h=12, text=f"Overview of recordings in {directory}") + + pdf.set_font("Helvetica", "", size=10) + pdf.set_y(45) + pdf.embed_table(get_stream_info(directory), width=pdf.epw) + + create_qc_plots(pdf, directory) + + pdf.output(os.path.join(directory, report_name)) + + +def get_stream_info(directory): + """ + Get information about the streams in an Open Ephys data directory + + Parameters + ---------- + directory : str + The path to the Open Ephys data directory + + Returns + ------- + pd.DataFrame + A DataFrame with information about the streams + + """ + + session = Session(directory) + + stream_info = { + "Record Node": [], + "Rec Idx": [], + "Exp Idx": [], + "Stream": [], + "Duration (s)": [], + "Channels": [], + } + + for recordnode in session.recordnodes: + + current_record_node = os.path.basename(recordnode.directory).split( + "Record Node " + )[1] + + for recording in recordnode.recordings: + + current_experiment_index = recording.experiment_index + current_recording_index = recording.recording_index + + for stream in recording.continuous: + + sample_rate = stream.metadata["sample_rate"] + data_shape = stream.samples.shape + channel_count = data_shape[1] + duration = data_shape[0] / sample_rate + + stream_info["Record Node"].append(current_record_node) + stream_info["Rec Idx"].append(current_recording_index) + stream_info["Exp Idx"].append(current_experiment_index) + stream_info["Stream"].append(stream.metadata["stream_name"]) + stream_info["Duration (s)"].append(duration) + stream_info["Channels"].append(channel_count) + + return pd.DataFrame(data=stream_info) + + +def get_event_info(events, stream_name): + """ + Get information about the events in a given stream + + Parameters + ---------- + events : pd.DataFrame + A DataFrame with information about the events + stream_name : str + The name of the stream to query + + Returns + ------- + pd.DataFrame + A DataFrame with information about events for one stream + + """ + event_info = { + "Line": [], + "First Time (s)": [], + "Last Time (s)": [], + "Event Count": [], + "Event Rate (Hz)": [], + } + + events_for_stream = events[events.stream_name == stream_name] + + for line in events_for_stream.line.unique(): + events_for_line = events_for_stream[ + (events_for_stream.line == line) & (events_for_stream.state == 1) + ] + + frequency = np.mean(np.diff(events_for_line.timestamp)) + first_time = events_for_line.iloc[0].timestamp + last_time = events_for_line.iloc[-1].timestamp + + event_info["Line"].append(line) + event_info["First Time (s)"].append(round(first_time, 2)) + event_info["Last Time (s)"].append(round(last_time, 2)) + event_info["Event Count"].append(events_for_line.shape[0]) + event_info["Event Rate (Hz)"].append(round(frequency, 2)) + + return pd.DataFrame(data=event_info) + + +def create_qc_plots(pdf, directory): + """ + Create QC plots for an Open Ephys data directory + """ + + session = Session(directory) + + for recordnode in session.recordnodes: + + current_record_node = os.path.basename(recordnode.directory).split( + "Record Node " + )[1] + + for recording in recordnode.recordings: + + current_experiment_index = recording.experiment_index + current_recording_index = recording.recording_index + + events = recording.events + + for stream in recording.continuous: + + duration = ( + stream.samples.shape[0] / stream.metadata["sample_rate"] + ) + + stream_name = stream.metadata["stream_name"] + + pdf.add_page() + pdf.set_font("Helvetica", "B", size=12) + pdf.set_y(30) + pdf.write(h=12, text=f"{stream_name}") + pdf.set_font("Helvetica", "", size=10) + pdf.set_y(40) + pdf.write(h=10, text=f"Record Node: {current_record_node}") + pdf.set_y(45) + pdf.write( + h=10, + text=f"Recording Index: " f"{current_recording_index}", + ) + pdf.set_y(50) + pdf.write( + h=10, + text=f"Experiment Index: " f"{current_experiment_index}", + ) + pdf.set_y(55) + pdf.write( + h=10, + text=f"Source Node: " + f"{stream.metadata['source_node_name']}" + f" ({stream.metadata['source_node_id']})", + ) + pdf.set_y(60) + pdf.write(h=10, text=f"Duration: {duration} s") + pdf.set_y(65) + pdf.write( + h=10, + text=f"Sample Rate: " + f"{stream.metadata['sample_rate']} Hz", + ) + pdf.set_y(70) + pdf.write(h=10, text=f"Channels: {stream.samples.shape[1]}") + + df = get_event_info(events, stream_name) + + pdf.set_y(80) + pdf.set_font("Helvetica", "B", size=11) + pdf.write(h=12, text="Event info") + pdf.set_y(90) + pdf.set_font("Helvetica", "", size=10) + pdf.embed_table(df, width=pdf.epw) + + pdf.set_y(120) + pdf.embed_figure( + plot_raw_data( + stream.samples, + stream.metadata["sample_rate"], + stream_name, + ) + ) + + pdf.set_y(200) + pdf.embed_figure( + plot_power_spectrum( + stream.samples, + stream.metadata["sample_rate"], + stream_name, + ) + ) + + +if __name__ == "__main__": + + if len(sys.argv) != 3: + print("Two input arguments are required:") + print(" 1. A data directory") + print(" 2. A JSON parameters file") + else: + with open( + sys.argv[2], + "r", + ) as f: + parameters = json.load(f) + + directory = sys.argv[1] + + if not os.path.exists(directory): + raise ValueError(f"Data directory {directory} does not exist.") + + generate_qc_report(directory, **parameters) diff --git a/src/aind_ephys_rig_qc/images/aind-logo.png b/src/aind_ephys_rig_qc/images/aind-logo.png new file mode 100644 index 0000000000000000000000000000000000000000..04c9a5241102254951ac400d9a0c3d0a5cda6803 GIT binary patch literal 40774 zcmd43g;!MF`#uaRAX0)zcPY{^bcxd4-QC@dC>_$>-3%~vcXxMpH&VaDqo413|AhCf zN0i0IBPu-g+)~w23eCs~-muQ{TuPk?OHvBwiCY=}`!~@1|1PpG1+&oaQU6_0 zp#S~e(MJS_%+xo6aTxt;mgnrk&BIbF=4_iz8d_4>FBNF8gMi`mXQpi$2CoI~uRLR1Tt)7U8d=hT9~MWFjcL z>MiZNHIk_C2w`t&EpiS&S0bLl*ESgh0SI2d9^BGCSRC_62M3Epo>x>a>cq%b4- zHFr;sZ@J~NdQ*0yP>WndSU2GT-s+~F6oWye0=F(^$g1c4?=EH2)uKPS3c6(|g|N=U zHGFS}XZ>pro8RKIZQi?Vnc|6+wU1qj0x#bbEGEQaiZkSD;c}${1me@~wxrx5`_K9k z#q%M56jOS7H2et0G@t^gwV^f$M1`t8z z*%g8h)Y+@kXui9IDjxn8m+dEawlqidb7I!Y+jjSPS4Jbjm(y>{Z|yik$B@t{lz#{YI70QM zDZ;D4Bly>bJ$C0m`nEO#-4)p zw3T!CXSzlYH8|(b5}!=RgH0Jc0v(Jd)rp`-TYOP0@v(k1A2pVE^88&xW;L*eAi%;8 zB+>P{)<`zvIIkmfxAzIQe$`l+p-ANnZt6vz%f~n0a>?T?43X^v#;ykYMo*vBt2etC zB>s0RuKh)Ds7Ci;%z5SfVBF3%DA?@^Zs}1PS?fsM7lBp0SmIENT*lbk^%mxKUgweT zcN>i|sOmF{zpH3yCkusve^t&2#=CIbCh2JlhOSXi&xbRVIDSirG7*@FZ!~avmWUW| z_AUh8wquho18a%fo2K!vtM<~;i|9fBfHLP*^ZVs?9-*Q7$&t-5i#;%L%Rl(^Vj2%u z&>&wgAC3SE)yETDtnx3|t?NeSh7Ru(ogU?yCl8dMa-ug(QZBwGTtqdhCcmazusgeH z`-UqxHM`2D#RXz~6#;FT+0FO&flGr$5)}&zr*zs<-HUlt5xHA_8D5#}gw|LcK>ZrL zttWwQ^)DI+th3JIbwh~ym$uqwzIay!tj0z0k*9wFwgQ&IL{%87oIC^Fu4ZRGJG812 zYix>;F4W4X;wf0Pq~`)5GD1e(Y@s#e^d@@4D~pt&v8ZJ><8MrkK+uWIg^_(He`_hn zGQ4SRB4@TZOZPThL1m^Gv)W|(Dv(!^OoUh@?Fzc2mN zMHUMG${XpE=*yu@nqKPqZ^%d6QW{DHepiq8AA7u(6~49Av}EM+v_Kbh)GPUyxhl9Q zOhP^vjjab?sj$V%5<7mDq?%L;q&o@pFY{mAeF8#g>K|_SyDd{(q3J3ejfd0Cc;nmJ z%Gc=smW&()iknMP_o16H%B+~?37X-kd^++ty|y@#n06KEE$uD}mOaen8Or(Qjvj9# zN6jJR9EAVwUqFQi66-7i^c6(futZqzN4H>Dwu0_cse9DD8^6n17;y$mNdDbBDrRSw zF_&UCSu)4Ze@$eQj4cY*8#uQQ);dz^d zi?OBGz>+`LQpR79{pS#w({^`c5{@SL%PO&c#gQD~?=~s~qLj_eHB2NAp-z(hR7x8K zE+Y9d3nz$qwSc3sgSme+e$@5<_V>Gy`*W-w2}Mgf?YzjC_(&OJRKvrSjJ_7JrqUta z$gKmWWPztJkUXNaR>vD@)9hU`|JRf99n@JwE=JnEHC{=5DLgoIa`dAT87yi8Ohn?i zrP>JUPq*s=eP%3G>Gl+NWK3a8ncIh* zb|QJavcQ$ceThLRin&ZO|ALGmyF+BIU`fYdi(Eu>Pjzn?MNu}g7a>T4(|BJ&=hQZ8 zjy|(`E>7JCMaf3+Uw%v2t`fn4D26-5<}mPWv|=1lebua+-{sLlVW&WP2e#q2sAMg` zODi7s4h#I}B8=^nQ6SY6?;EAu4Shq%Y=$&$KE&0%hH_!^Ne%WV(WH_|N5H(Lk}f-s zN)$YI5&Vb$j%~~!nN}nV^Ua`UQlKwo3}8PP)@%?owOCjo-#1l|O2%X-R5HUQeGu(J_N6W^|K>X-6dZaS1?8;G|I^}i#qYa@^YT4(t+Kyc zE!by4ssD32)npP8 zJB|jNm4_s9M=X5-;0=``wF_|m^Yl0aV+gMd%lf{m%9Q+H6z#?15s{P5=q~=aH@Z7Z z(=$iIDIZXQQ!?i9xNx6tmYo+$T7ejo+pOXKtN(UJ@b9;jwq-Oo)&9jjV{4vyHvF#5T!T5UwBs*Ky zB8fn_1gg}Zn%YOmN;%W};B95_n1^iO8$CC;m-iIP9GpG{f$|~lp?{5}4lR&~O|g0- zYYKa|ZHtI@WJNAai*OG%M}6E_%Ow=!=Fz!(jni{!kryjGBn zUa|JR75<^Y!mGpa!%gHvV&LcaEl1phx-Q#1fA)!I*x!S>w%dDrlO70f)(I3`c)UU47ks;1+dvEJrUX zjPmU($3yNQaw-xok%9?kGjE#__GKzVOaXM&k3ENhe}to z4o7|dx5sI=SzQ~S%ktY_&eu0Py>bRVZYT&i|M&Svb`k>g<35_ zn!7#rq`R#h=tSKKb?w2R!JoEoIH_@fpjwxGKOqgftfXuR zh>TNfRU*NKrLt_M9*m{*M2l)y~v6?%b z%x6f5Br*k>#4{Xl9&*|j&W*@SP3r*ruH5262YY?wr&7l4r@6E-#A7r!e>Ce=d$&_g z=EX}`DUujDZtJ@uCWUfC5d}J6IqLE&xqsQ{1hMgCsR71Y{DYN-V1-D~$CxSDoe0q+ z$?-w03XHkBZd?LowUJW%p{Fn-{Vlb}$2zk-<(Gb0wU;;&q}t>iu9)*ZXIgFR_V^uT z=Bvf`tK>l8g2zL^cvk|Pwt^PG3&c#z+Yag}Q0V;>$mM%-E`CPDdC+)LdgId~ie0*XlL#iW(}9XRK{ zG1h6gN8P3|x@laWZ?b=x}zFC^_K7{K~>auvaHOd4zSZtvnu>FrF-GC(Z8WkDS+iMBXAA zt8*PL8AXlT++;D~RzE(vm*1sXsoEn#onC<-LUKozJun~+4m7TkijvTl z%{w`M7Ior&kGNhYXLU_YTxt(^Zs6^akqjmVs)6W|hOpHx))w>G7+g{I#!6W#q zD;0^3Hu#xf`{R2Wo^WvtGz$tt|poxL;1kNtWF)nMSX@`&E@Iw3<|$?)pdK&})~50s8bNlcn48oD9$#ZW1Bl8`r18C0zG&u4tcnr@nWfx~ye6 zP0aMnJs2*xS<-GYv3sEZaO(gZQ~tXf2{{pM*NW}Ca_R56&3lTtzQaF{2?4zd!d=L% z$xDG-Y@>q+=5(ysf=(!{#eMi^9j`U4dwxsXzh~#u6VJ0o5pUqYvm@LOO8BslWmW7^ zw&L{WIrR&3-N}CGwOz)jBCcfYs<~w7tPWdtUED_m&-1pvcs5V{j%vZIt_N~6zS2A( z4!~mWyEb|6XgnbvVCd`1K>FNgoi95S+WObdQAI)_MxYpm+)ZaQ)j0} zdRkZyy{02t2E7peSBPO3T`y0x#hpSAxW}ffP{86-Hs(}fu?mQ67*n)vsd-S1|Gjq^ zt07;W)5Dfpb+j|jYtCx|vsnU;a@Kmspsq(JY{#eg7JYn>;PIbGyNgTfo*Kd?{ZUU& zQYsX{;wuAW*u#cxj_4l~?h(7EE#`e&6fIFHGrJWnW(~?cI+(xPY_)2Z@GNuJgm=J| z@mfM+kAU@W&L3XyU)^Z9d@dc2UtHntM_~F?-IM3kqo<6?`&mC$%vyeW|F?{Ka5mS+ zupFf|W!I06rj%K^_;o~}t2CKjFj% zD#$Bt)Xnhrtg}}Vf#vTz74u6JvdU`*wTsR2uGinZ^jVXLjxI4b_oZCBbY>CQnA<)* zWr5#jdZ)>a3J>hr5*2^xVJomD$%S5K)(h?!BLUrT$hxXN`pmkuDY+tE0U^q1WR_Vr zreEisu5}}wI}rb4@}xMe+#RGonX{?QjmkMF;*kSB6R(AZ=RIvOEcK$p~O?%ecRA}C3PeBs1c6v$%s4(EC&A6 z$78oaesrsxh$@ol-MT$oGP3EwA7P<(8=#|Z3sG6wxR=sIW#_n`q2<6U<9&wG^tf5@ z1E<`F=`3e5$Fpn+GL2NxDp*)#fxl8L1*nTkn&4L_(|))0MCv!QKedn21+cl0y}9$u zVSfHTo%CqdAy zX>J#1vt@y|c2+#dx$-0oWw1TtMbZPf>HPXVhqNM0Rz?I+U%y0XyRQ%py5v(?9up#( zFy)EOm|ngBv}LqB;U>H(ZHq^BY%vNnreK8<&DMx2QTHYWgqufeSYGpCK}i#h%2eF_ zqE0v``_870b`|AXc(`Gu)=8Ofoya@mF6BT&q8TISEyppnz0Rx_{^oAVj02{R4YEm}n7p z;3Bf2wrEi@eFf4^A)gzolYpMGP}6J1j73J~ZS#SyCf?*WCEN;Y#!4K;RO*Q@UqT6r z!Z0=>w1(GG3{2RP7VO0g3iy9;`u)_s);N26t~g61+9R)Az@XBT=jFDODDoPLb0n6S zA)Ka8BEe2&!9q@qHk0HVz3c2pqy-=^dObjv3&9bOkpud-`@TC!^P+LF3abg&Y1ruL z6dUawKTfZyP;f)=IUHeEp08_!QK7jTy0FdZ9?X1&Bi=`pB=_dGSNL9a4Y+kVuf_Q{ z=ypqhopsa^c~a5*0g@zwl#GON$TAZ2Y9Ap}bFL3XYh?zzNjI^k^HWnk24QJ`Ylq#4 z?p(&VUvrc-blM z3qEoSTOILY2||z1tJ{cWs9l@SfYnQLzkrkclpBOL^7EjKWp=sNg7-D{ai`9j@-`(vSJy{QSC=Nn}IKoA1;on4rp5@D_^am=d}iZnX_KQsd(S0CeTh zc3Rx-Efi4gbUb^q`>!}|^W`w`l5>r!`Gg;bbo_S4-~n~ZYv4cuP4HcAHGyc!ra2C& z=O2%RwrdZkFs20^L%?-*vE0D#?PNr-U(y}qw)csfZE;Fxf7MJiLplVS4C65rc z+>!lw{fJvVWr~#Qc{SmeIaqDK;n49OLsyyTa&CsFLyM{=-We8^y`udclB9v+YpBxW z{c!1(SSNLI611o#GLt`Dt{dA=p$7u$rG<+AT27nmeJdZpRNzIOD5f|$GZ^Y6E$>gB zV*K9}A#*+TRN2l)+sfj`-8TbM*_%bu@&yd_R^A_IU0!*U-@ofbZYxEPi`D~rF{A~8 zW6OB=31-7I(gd3R1lSY!n3BF?YR{fV&%#}vKd z8WEM9q(Nl7mTJ(S1nMO^pbT*bw;-c`E3z9#=>yt;#(FlG3#TnUv+Zv{NLd{li>P(b z4nWGP^m^s3w+G;bS3zff-?rVm;%>;5h#=)QNh(UrUVW;*IigXP_1)>e3n-nWa1c;# zg2mP)Qzg5B1Xzql%I-J2Yc1||&EpGxErdpYQ{gJ4zeuwo2vRk@JW9@Cvt3Uf(yIeX zedi}|15Rk!ofB(m^O?@##jTn$5G>^(Lsa`x{MS_M&$=G%1Tr)Sh+_QuD4V))qjDjjE;Q95ZFsHOrA$tvEv+!Z?zqavOZ=~qk8B(6r(8b zwb2hSX76boNJ!xHS(!}4=YS&)e~#T{gijPB-=12ngXv>Svue~sg%XD1aRCRXAjsABu7kk9Fojw@PJDr5CxsB z1aK-uYgzE1OqpE=RU`a9_3$SpGmypzpN`);I(zalG6rch@B~7w+B-7)Ii8zUji?xn zL$v}0NxZFra@18y@)S^0`Kn^_`FPga3@Xs@3O&=#DrxHqOW$6@r_iqYAvdr~Pn#cjq**OMrUNz|=#E$d+aiV%&#C&9i(|}h?#3VX;hjDs&vGo`qCSA5~0UL>vYxQ-?~jH;&GcaY7XB`Qa4^s&)nK9 z+q(Z5p1_)H@ej=+uizoHk_u@w_3jG8oFSgq?J+Gy5KHl{CN}|bKnuj+M9Z5&L7!gz zCT08zpD=lDX9pCSo3S$)f+a%HEx6y@CF0nLNM^L1L=v<16kbvWP|^Z5rqUbeTSGkz zsCe;7=g$WXJlrirb6YK5@G9L=AfT9k<4PqhyXEfmy~xx7z*7ee;+h}XP6&OE3ZlBd z6`0?_r37vWuBXh|$tY^Vr#Guj_Fr)_+!JsZ(TZ7d% zS49HF(iNYW1E=-Lgo{vuq~yfFz?-HJAUONotFg(2_bvJJ5E?ekeE`laNtRC+jBlW9+(}fKavtO-c$Ip9Deo4Si?T?@K4*h@bD2f2a*(^ z2r<2NGVy3$WQ4Tw*fE8U08BMd=KwAe#x1gi1cWhOshNC-*%;+8OboJ(28p@Rs`LqS zSxqDWaTBI0{#Da}2;{G|s~k#Lc%X9lxG{C8hZb>Bk(2-mzW%8FRN9zJxVwJ&!&_>` zuw8?gS;g1TI+>0F`E47o|2qr3lt*F8{9md?{<6!rD7JBLr4E_oHhLm*FsDcuYG6}R zoK@l&E-*@b1nt1cV}MVfe)u@Ozaf!iOX`{0(wL#V+c(k2=6r|hT3X)dO@rg@9 ztU*C$c0Tk*<`E1B14ctw)@Lp}XpemJLeq`4G@-m|+8Rap>L`Nmm`c?X)FF^Bb5IWO zRiZ@@fd4hGv<%sfnHgC{zERRWYkoyNuk>FXeuV8!AfpLzb{SQwzRO3QDo-N(*Il%R z<3O#R2A_yAx(DN4!I72W+Qe|9i60hEv{20o`NYwrU6Cu}9)B^M4ec!fpQG!qEoTbJ zl@Z$+XaGUv@8jPhwp6sVy=2wuz4%Qop2WLszHfYGn5TyaX<3Ey_MdnPcYx$Kcf$3$ zk6ekli+sy@ucS;Mq`wg`I`^Btf%v<1U)T+tPz<`k2iffWUD?Qi(ntr615onlo<>c? zucms&Wy~npU!aGE-^P1b`#SVICem*?eu7G6>{i`j3`@lQfnx?l3VJjAYn2TffHdRl zoBPj1F~US~qVI7J!wQsstsUHR5PwO3y0}Wswk_tEZ^RzN+!s*Hw^y!ZB?PE4?c)xO zv)+7R^{F`&0P|r05S|~6D&c+shuRcFK*M#rlYTIKYfKX(7_LG!es(wM>|DouES~XO z}o$b|&zkNj%B})I8tBe@P zsujYsC>00*bQ^&#fus9xs28Kc58PndhX!Q}aUjSr&omu{(k5f1|G`ZD-?rBw@BzJK zGEpyHNzXR{JD1}S-^B8P$;!mazR;jXSSsY;!%yZO^-5+tzrX}RVjj_2%Iti*@#*6! z2aZFr($>?`YlwZTs;5eMzqVS*^r|nNXXUL6NR?)rM#WbtwSH#3hoXn+I+4oe0;vP6 zyALsE0Y}ADCa2@>cfhfvYSs8=)hqE)o$MI*IYsD4NCae^kB&p5Cl)01e)YArQySXA*RZC)$=JOo|JIK+}$>vihn(a z!vqRh+$n%dl7iUMb=2O!X{5gx@H>{$PxNoVZBD-8(ys_6{Ql)j(b)s87RcVKLg=N- z0E$^yJpcsng!jXF1A_%r$zLdi3}Uv7CbDylaLwW1n2$+B)^vhVMypz$A?cxVKBUv=v9*QQB(aLXT2K~$?Ie{K~_{q%4d;8Qj?fl8?i8~JV()% zqrvZ%31FLZt71*%BieyaEe)vw58UnS<2lZyZ^X) z5-hdsIj>D!mX1}{^z18)${tJg$Fbx2xM51${5~~H`E*0z)1d~MLaNV|BEn0c53VQX(k@bh1pH!o5GPYm?g-jY@=PW^ZiDT3i;EQV{33!Lq(v6Oq$r;;+!DG|%F>_@*@%K4g{ANB|F`|8V1Ym;z-Yzgw zxmtEL+q`4kalHsrWZgrfN9~Git3DHo7qVd0j3K6IZi;npbfp-A&z(IMf?#?WlXu(} zn*wLULfNEqI27*Lt>Al~{=GpMut^%SUTew%!7NltwEC|UQJHe;h>~PvtRY97$=-aA zZpxQjaACwoO z-F4k1*1&uR^VIA+YXzXjN`6LeDMDXK7IK}wv;tkv%Tr!iF?$3ub_W4vQ>|kl{ zx#ReCPeh6vxO=|n#51_09GV1=D=}v&)=$M_AE-4F~WMdoapiE6Mi}Fl?L}g6fszW z+5yd)!3J%!n$W&{q6sT5Ea_-(VZ{l7K#r;@m5>jA z9j9<&g$x;RTGsz83&-SfyHW9FX(iA3=)5;n#vjseb9sp2|A{csoZj(t*W`iAt4kg;1j%A&<=EW*2)&5{!e5IZ!`^`@;QQm zUdXx{Z3$q9y7mU@9mP|)rVDp2&o5h#ShdXPS-kc*Ca3Ea87vt`fKy#@za`0zo5ah{ zFYZ2v9KM~#TjF$TMF?N9pm!u8#s)+xst#N!lHqeJw5|9kLti!d7y#S>h4+RBP;Ja# zMujT37!Wr9aRIl@(rgG3QFlU~O6MyJqO6RoxUNdVz1n%rsv-N%swv3fO_vNn&u&8h zwD~;>l>lRMp1$uLZB<>Y7LI$P@CO27jVw!-?QaR48&c?%6T^F_wri&sH^mXFbax$Mt3~Yj-Y6NFy zUEI<~e3sg9iTU(pubH2M%tEmxONLh$k$Tv(Wt=a~1;95cOL4KmKe{=Ad%@YV@BqD> zKB1vLFAbn*lR_(|7OPwtW-IKGdOmf9)koxRyK zT5b`?J@bQC-_sAnG^6EGuH^ngSSZA!O@BZl4&+fC3JE(7NmNt~?aHFb4+#p@95Dk$ z@_Xf+1$bHYcdYZ6du9A5j>~MRYdE~&kzY0jQoR4@+R)Uwt>*(O%nXg+cH-cn-u{m5 z6gZ@A+A>T%8hkxchKK%U*pyTN9mURFbkztBhupW?EZ&_SU<4lIf|WK(7fP+h>{1hX za4q{|ohs=K$3{*dyJrB-R6aZ7S%5tvX(m6xPd$;A;y1*JLQGDhb|{S(2KY(b32!(< zRI2HPI5U9v&>1xWVuE*v^InZQ>HvB&ZPcQ>7c}GyFs?5kv0$8-DM;0PwjPCdSk3Qy z&oI_uR9|-d<2BRUcb8_DTxc_`o=Xe|5TP@GcV?N)8asdiR5UOwmJ8wjxyvDEdVVa0 zT6SwTJs;4(_FaY#b>@0mGupzxdAB|8DHfM#E}l%RDs0~spFg=Dgp%S82|}_tLJ8`g znS+dvBQulk;q;~=NC1Z_XU{-SY3HKP@G9@?DT2K6oQqhHzGoQM3N-2wOyc&}eE<7y zWk(zU$}(fiRU*3;B41tWE}iLSiQuF^9s3kLa1D-R+dQqeL0MBbq17PLKely5{HOOx zIRpqkbzawipBTbGl z8jG9_EfafZicQbF18Jw^7?oX$Ccx&#I zkB>#EZ`bg-ZvO0q#x|0rjAn>jZ zZ^iu*t$Y5o$zOFlb)kokBtpB}ZH2YU5pMs(BNOif>dU{i48DNSK<|d6S3jVl734|C zX>b*$zGU^7W8#PX;%)RrJ|eNWGMH6&2z>`&nc!QSd)2yHM^_sfTJdSCLgQ1I;SVAD&qSMmGF3s{-Ci!(H|i7Z|}8sC*Jew{#TH7-8 zG!*yfq4Wq|2jaVVbT(zdG*Qv)oe}kH)X!p5fECAi%%acFqHFdwTU;wCo!~3^m>$jA z9B^r6upl60_oWEOmCt*!twCB|p{U&Xj=~oVgd_U+X-f+|LBwz3@>q16L&O6&W($6LoPde2Q;6iBpuTDW8DznY`r$()_VH` zEX-z~Cl>mSMW*p-4s6%mIwVQ&Y_48herB(n3XIo)Hn3I@-L@jIetrl{X!RKJRh5I>$>jWHd>qt%WRE{`B*QCfa4VF85f_Tt<-eLfqSG_Wq z{YgetFC8_58rP&mhIe>Po`SBelbP#$Q(1VN!2BrA#3VHiD{IbF^R-xCQbMuZmybwg z#89W_Rm&$`Vp5V_rH?%ki8{_cyKA(PXTnxiQoh0BIhk9P@Z5__NX=fBL(Y#7Klv|h zMkhSMNqNly^#k>s5DH7EcmVlWv_6WtHCX#X{b2TboY5}HbKHxgTTmC+&{F^^<146I zrWQ|7H#}uJzp9p$1u9yczn1N$HjHkuv@-5ZDrqmWy41=BPgO4^`%jrf}DXMzs~!0 z|9j4j6W?F6?)B&~U+_&Qubyd7<@@P@839Y|_aVcmL(iUD5Mo*qTMM){)wD%A1}8Pz zCx$VbUV<*zhAgw6pVm9L-yi4k1EK*yk-qp*cGBc3e`;4oGfM!_t5=9ku|;Ww+SM_S z?jBNo{0MOIWygE^_A4ix7v;4uzsB+APQo3(|909yc`J5WS|)1M2QV5rQ93cLFVu*;oZY zCj*HWEo@NCi1d3C;SE7)Sx>#r9{NSp}v{l7yAZU zm{L9?T|Lw`2Kgz~_Kj=b81dxmAVYLg3X%%+q`_?VB(vv|8oeKj*5A{J5DyNGE`L0| zvYIsN%zCvr%Dum?5eCHC9;U+bEqmwdH!d(*qiQ~jM$}kvj{*HP44(vKaUKJh;|Z~ z>Y7ea5#*m#SCkqzzIqjT5sInvA?Z!4g2scj-U`uW$eE>_nJF-u;k)k6@OcbrZg-;x~s~xbO8Kr;tBh0$!D)4qw_AC-Qx1NnsAB`U_ zNDi}YTD<8gLYm1GpqAT_Y)mhJz+t>-tV_Au_kRrd;)csx``q z{<1-ELV~MBZk!oe!g8O6hIE4v>_YdFb0b0T5)P=H(4h3oJ8=FgrZwye$RuWjrQAKN zK82j%>Q6X!g4>R7FH4_9dM2StET+=1V{hh|Y$%*b_%_s_p z^n)u1DEc$B=z+18vCJHy9D|R-Y;i%+Rp!x3iW!-iI0OWir`P9S&Bfnv^gzD*J4Os~ z#LCjPZOjM~Q}p8VTK4%}xtUIVcte!)^%ucd=z?>-s17J4u%5HDY7zc_!;`(sD%>Jb znq}S`C_!R6 zuWde}1b}L<(s8fqYb^n_VjR0s4iNr1t@CD1OPK3&(U8m*sHYJ&=?_IGY^{U2|L0M0HpTSf7EGY z6tkE)N~6VA3!>-9U>PXIsm)}Q@)b6U0f9Hg!CzrHYQ$}EVs;;;oZ7{yATWKg$jkLg zDlG6vV_&^#r1ECR}mv6_^cu3 zISLl)e4tCs^Pkyot1tsI1Ue05p=o(|go^|-#R-%zI&+}9=2w7RaY{V`TF=~;o07gq zpXGOpb_($ddilju35z>>hF7$QnM?|7$|W{U383nSlz&zn3*P_e8kk_f5{tfG)Y{|J zogwEuL6ai|APL_`hJOf}CCx&u-gKFrZAk@7922h4y#<)T5=YO#1O=^9Y0epa2^`F1 zk&geJ=wpRxWNuaeA2BMD*NgOSqFewfm*{Q|oGF3J)9ZXw3$4-HD64I~8=nt?HzMR{ ziwvEv0zZcwU}_X_tr=WChtPep>q-O(P;hws=cB;M;ki;snHtSqCp@eG-fKR#i$ z5YsTt6fz`aPVAixhDga|F)qWp01`TWSWt?wDNz5setpw1vYIzb>{$cSKGmR6ozvG0q=!#c2>HVEw{H_`^zkx3*VP?)eq{T;FF^;Dq|ky9{A z3FgaD|B+BrB_u9=5eI8vxTE(N%dHrgm-D5_^mK<8H}n|cOeS5fHm_?0d)>^nI-(?q zGx!&i38F&6-sQh})Ab6~$lToF^Dap;sMKa2SPfFO*xp;SDU61A0lPK;WR*BRm1pvCOgR2DE=g{@H(6N?g0G-zQ4#F98YX0ZaNsdRDh5~Y%E(ILp6!(HShF9LI zp;1b~|F%sU*tSWP+i%v&M!e$KIMeGpFWS4wxeJWYF0^Yr&~jLV212Y2v(mk`Q`W|> z7}Maif#sh`@gh-LYY$otl%??iTJdt(z5zbTn1<~;X{VCUr?Ox{Gt>Q4F-iYn0F#_N0;;re}Ynf)8Ij#9Z?}k1KDp-Du!jR!g4(b2h?*atVuCs){-ZU#SPx)l)5V z)eLLoZ$3{jTMlotIYm+cYGu*7vpXg1u3_HbUUXWgzt34A<;WeCD-?}Mp^2(%`wU3t zsOOVVyrd>XLgzW|)u&}O^fCW_PkhA~C*4ctd|v-qSRC&YfXk+p4*i5`o$0uY%aB$+ z5;Onc#5w`u2`pt@J1CwJcv>->mZn_GHEA-1+ z_gj!1QapKLtxPQ7n+?}x3uvOUvC&S<&Ub)`-$jz0awB|+!E6yqX?BkT8$0yD zzkhb1uJ}cNq_4OpugENEq(qu0kJL2!aFF%Y3g+*uIwow@4J!8-@!QtE!7q5?C7pH2 zIAmKz>C-GHiSxja5-n=#iwqVRAvz?1Hh_?Y79+`-=>AxwcahHTdjz;qVull4?vK9$ z>q((TX*n{dY?S5h9dkFnluWB`;sHWOYvEd)91YDSKs1Y10& zs>zz_&}gdYKRBHpxN!zP`Akg%`kh}I7|qaR(5qV%zUqozGfa~;z`TliIVXp|Tf8~4 zHqsXe!3d2o;~h55sF=@4Zui2}$y_`4ow7}K@8;8|($MzP(kH2$B$yOwRDeg0 zrPu@@gSeU05yPKchj&Me=cm?cGs4eCx~d&%Yhan4S(dl10394X+pu%NE477}8Qf_j z0N8vPA@7gTVF}2$(0sb!8c&jmvR!_zhZJ1$hpPIU@TM}?pIjyBjWwj4DIix%D`s^K zR~pPZRwM4%n$A}hW@d}}0wdw&8|K(#Q1wD?I`%hEA-mgF6E|&ujZPC#b?@tS0o8cb zK9&Kq>LIU|&*Ki&=~ByFZDA$hi-pc=Z+s9mZ=;O@p&*9Sv%e-IF2eC5qy}LzIgQPc zH)|Aal;;lHbZ-Deos34F*~J;>FSqohz#jwt6JKPqEbrM5Rs7#OIVXH>O^&vzC-_Ze zd_c_#o-~1xpGx#xD&^J;;W=t#C=O0T%3_qtVj`%Voa7kB7)2=Wb^FC!RdpOe)VevM zXa!ffEbqXPb;f*y97+knX_}j{7faxb1xyW~{nY3BFth66T33E(yrHqU0_)1>hiJVm zV^jYL|M+A661zk8(zR{#Bptue++ujVYzN`XGp8XiCgg4E>j88%m%unWpG`B?Xj@OA z6EAnP>`QrdLM_Hr-vWsCU>YvYMuqQd+2aB5pZp#{yqWL@cteo{nK9JG$pEH%BgN~q zrH@uhX;g&=AuY&Lj`HDFFMokzWx>Ed=nxpaonk(5a-AB}1%_=+XIP;}HYbUQ8|wGq zURlt>5y#)TQz-gR0fJOT0?Kw^_D9-`hS1mhkJ%*vZ#+~pVFJbPMl5t2MOha$UlC^c zxiu;vSS#P%(E%j6*Uynsz=N-@$)FmF6)6W#jjM|?iTxhj`zZI^GmTXu*~_q?@I!1l<28;`>E58MMAb8TiKqeVqXOwh=oWCHtox9V7gVJ zw*i<~)Zd?X_0D$7uM-9JdawD6vHCOEv0_$Mxq0JNl8V7MHH9A;Yi1y-DgGJidTv#E zM_AJ8dWAiamE2J1dV`)%$uKh*u))0~dMyl;n3F=W6Ud9Lzs8xq_Ps@&GC#R8s@$l> zgRvKS7b6#_jQ7(m5_HM@(T|h{a9cgtr66MERz7WWxL0a`lGu3{LoftU#y0Euy~?!s zyK`=0{BQQN)1^|Fv_Lt@qBNrjs2*=-wZ8Apur zJONml3FOZ=KREWG`>vtdmEw?;+azS74pQmP z3?|^r$9OhBgnS!{otPv6Y2O#UY>i}2Bf)UF<5l&;bw)j3iU*WC#dU-}H z9JtY&1S#IP9uKM~W$K+cGglrfWFyH*zj1^x-+&Lq#*FeaiZ*D3gn)r3ud6a$$PtAp z&B{(tC=|%Mdj5?>xOXcOs<{GTX%^B}qFTs2%PY z0-Yrpr2F&o@(}U`4!(jREf;EbDOVaaVwEhm@SQI zaAbh@JV#*e&$3yiDg~R_o;$SPc<*Dy*~}kz+MnO-WpzGfz#U7iW8?phU<{-a-LP;@ zeY;(n#P2RQ%39smx=sSfNB(=nsls_hQJ&YxAPFOv5!Lj#A)c=u-u$HCusAiROeN!& zi3HjPHTTxPqdX69!~NT*QyTiQCJU$r~7TL#|El6q4_{TtDv9+3ik$d!P#Y$~ zDOHjz_mU+Yi*V}fP7<%8SAz+AdWpZnV3_r6eJ4jJ%CDByd5PX7so;^aw4`*dSJ_mir>}| zr=Ax-mg<~EG&jxT{`}`!@e*_V?e#ru@wcIbK>PRIeWgkE@QKN8qPlz!m~48ur8EfVAA zY+8P|?6Vl89l3V55n3g(h*6Q*lg9(^^ckmiBU6$&M6ci5$dWX@LV~KWO=rWZC!xV2 zi)6$UDoN3Ch`1i`Pm^Nf?JNB@Z{hmcR=KtS=2MC!UR`=K=WEu#k-yeLzlQLUTr%vs z4oAF)XA21zQICNhjH12B2_1NQXm>X2YDP?g82j2JBXnNRjLuO_)EW}NOY!danpxS< zj__Rif(=UQK)#BC91yF}1bQi3l{8 z;x$Ml;Y4nchn!w57}^d5e&?E~p1H?6d5g6|MmJpxd{`#Q!3wvV_;yDs1XI=3!C7TC zlr`x)2?tElW0*K76)&G(Yn_Qfn$WzvK_1^pn7dt7)iSRmQ>~C1w zhYU_H2MFLP&Fq#}XL{31ggYsz;d^Z2gD*9SG=6toVwXR(tT$W~6~))VZ_Xc%R`T>le_&M>eVh;cD^q!PZu?{OTU%)Ycrq_lbz~ifK3$Ht9Iy zELV*78hjP4h*~MxUf1I{WUZs5!LX$4+WnAyQ^J;D->zJ zOagB!RiOFOla4 zu0eQ~gixja=iJwAJ1l9QYC&NKL-gU!fkDH9?g}OTy+k8YYHZnQPsAdYsIl(I}hyI=LA?jQ z>dEqS9qY6)=`A%!U&YX#h|u))_}8obeKE?8BLv4|9!1AThz?M|TzV;@2X;T468 zW%aQEZ*PeYUGF=afAeaH`gjm%adTN=#1EVc{G#R>rZ5lIp7rKrVaMmm)FvDA`PyeQ zq=LzBs%-53SZ&Npp19%FFx2&fGLA(2Ma{!cvz^O6Z81TO0b7z4x3*ELj5;n=R- zm37O$)6Oq%B3(>lNb+Q3{OnSK%q+ zjI7)0PJxY;Bj&i!o^Iw4jap8b@D6FZoo^a(G;x1w@NCGg;|lBgclHcm&%rNwJ?7Y^r#LGOx z?@Nw}imI%UmgyuvQ=Myv=X-Gq{ontnlNiE`6nLmq-zuf*`8JYcV8*$Vem@EmeW=WV zZoL!mzO?qUv{%3u_q}H-Pg^HPi`@uc61epY!iyym?W9c=__|fFM(S;SLyZi+@DA^Z zdYQlFJ1qjtdQrZ0g}9H~3B-HKf17wBK7D{!W0aI-!reB`e>q}Yi8~aU7B*CG!8377 zcuIG&z5N6FysCB@y(7 z`{0=KF|KkbT!l+3t*F;Lc0HAXC!%dZbd&zbBZr*PVW$Z`TKjEOfkYs-{pEQ0i$?j~ zk&YjhQZ4D2tcRX4*J13b@Oj@iCmzI~y4k!>60arbQFH6JMiJ(q3+;gtO1!N(=bH(# zHED1B<@?bnOj3gDL>pY>zQpSZ$fGsDZ2j#e+Ckf^cH3(V?U80A%z4A()si5vczUUP zX>&lZDH|z?j3AYwotY_jH7kfs`*U0M=SFTkL0!~sma3Fu92yPsriYfY)OHB0^UvRc z58Gc@hkYP=m8+~bYY}a*U5k%ynnFJxVbp%u4(L95%ic)3WuOStbg7Q$FF}>p4k>{o*RhedL*AqDP#=O2b;8*ZWWGsHQ!X19`MeaaA zS?oXr){t>3TIdG-Ue7r`GG2|-mwmR)ubJ6lmSM?_rrhGXv4kJ@w`al+CV zrg;8&MYq2aQbcM<1YkB5aNX<)pRi9`K)6-?^ffckr)s`p$NA`TW~%qUgYvB83Th}n z+5`==183!TH;w4z4o9gzgH7)@D3+XDS!Lz=_$BHHDQSU_{3{&fb85*l9~1p%)d`Hz zEmI{=yeB^g3@29T<}ACtC|9kNWr$9MVHrJkP2l&gfz#1nOZ*GCQdIT3_#Gn^Pp`zIfrjvKpAn&1KbA_}MA z6HC00{3mwcg|?NTOV!|TB}PFmW&9xn;#K=8diW!)_?h~giZKr=BJ?@}9+qvkWDrVd zJNQ#tLZbYO&zGWu;-Z0Oq#sc`qV^AW+M<^?@8qLMGc6RE(@k&deU@AA9G+M3^i*uk zw72rL_m9>YW;*Y?Yr}R* z{iLl3g~HWGUsbdj`*|gKZZP3>i}9%Ru5m$3${lG)m0T76G3*%|w7{W{`y?E99}KFHpAb9y`o+w9j}QLrCMV)G5`FY zthPh?26n!(q4ytx+P2A=2nG#ht<%e6L_;QG{-RsdNc{6t%ykF;;DgKR>(aUCq zYKXXXc1&kj)?}aLA{qPwqeQ+Z&*y#je(^_=pX~Hi>1#E{pmxTQ3`R&(JM<)B5cWG_ z9}Fnv4E(RA__AxCZn{xp#H<)LH&i;8EO&R4@bHvXBX3KRWXA*S%Y!qeTXg_Z) z+q3P%YJx+l2;0o;_mj%;d|Qo&D)HVLc9=#kxZlPFhYiCD&xHXY0yjs+9czIO4PI1` zBfjM+jK+%ewu(&QH3K2+)UNW}h)e1AEIc;+r!X~pLF4xq9~wN`Itf{0VHKA5T!|g^ zgpmYT^Fpf~X|5;)gHx?p)b!saHYKnV`6{b1;e%m*cvN0O!V!@^fjr_U4lmJU*Ww*V z{!`~DOR6wpcS!PwbuOB$rFW|bkh9q~gM!qeUi@^WxUrWND(HDSx!4dtvzO4tDwf>a zpBENnHC#oVb00?dB4YtAdw=#P<{_!&nO?eK8JZprY(!F8Hx~dlpr9@*5aBNyu7+z~erHtqsT1KX2|COea7M&3 zo}SwDh8ca8+a{siT-k(ymte~)P2G)%-oq|=v&Da^+EYo|8q4MTw!QY&g24(C)d@dT zKc5dy5{ykH6(T`xtZB>3-WWsb$oUAyZ16@*2C zYJEvzpCk*J_dMdzxurg`1=IiiTQLCt1t(l+lEF0Xdh!~#k$Y3}t$K{FUsuPiM$JTe z*sMC5(h2wss8U@B`PDZe9Nyh`&d4lcfZ*%II3luNe0 zW6zm%&Z`-o`Qeu;+@|a#H5r-$GYO2K+IG$eku5(VOj7}0n35T!Ug&K0&Y+{Qb8#0g z9t;NPxVX@}uLzHXHo2_%@3h4jQF*B&deY5XHpY%s2y&3?B14mpUU_ndy$MK3S>koM z+~zM>B0-_U;aq#yucFn+?JPqAw(JKL6d_6os?O}d5HNbBHm6Z=c=WO|1Dc02vF(c|IuMf2-c?!jiEoTSPD6CIu6Mq5)% zW2f6QKf;rM(~hA|Mu5OrDYU3NYy zF{opR5+RqaYcZ}1GQmo{JUL|z$`F{JoqB0G-P_y~t&QZ^(0ph_xcbs62-W+c>NomhJo)!;pYxZTf@7zyd5fz<@3FE;ScWDvCD#;| zK+CcbUunEh3PmvhST$x&V0YVqwbRXZ^PaV zTh8zieFUDtfNop4A00{f!?BED9{gti>RUAj=?v z!a7nJu1Dv6WbT(kI&o->=7O5MjQ7si-%swVfCAV`Rvm+8XfN0ZTv5{H6S4Igd&F5e3 zZd?yX0Mgjt5KRH}o!xQWt~mxJCeDHb!NcJjmVLemSc=RH`j-XYxK`9N<G$v-5OR`fecxI>nOJB0@BI=D@by`dk@%1^TN}QskI~=8aPtS>k(EdE^ zHN3q>CxKNvs%oIgfbA%snsA(_Utl@9S63iL;EgWb??~Ekr8i&@b{FdJMhv8y8 zrHZptIb)jq=b>K817Tqn&Q2t2QYArEZ+`Z_MYw|2daUZUil`TOs+q=kZ|#jUkJ-s+ z`uluC_4>DZd=8SZFrq3G8n;;^opuLmS)tdIe+0!AVN8)#FT5CA%wc)5L|Fmf~wRIz>r#`vM<52v%_C+|*h;#-$7Lu=?{6XK+gdUTY8X9rwE4zG$D4OIPQ~fGPU7AQWfSZ`|c$x$>*UK zGPbk3HAB)Uw{bj0+P`Nf^bi3ZHbQfGlv^P&EcxU{gU>HCk!Ko=%6~$PROBz$T)*Io zS#5h{GxbW+n{j>%5{2v8cl(XLa)x3g2bD9|-}PDT>YdGXeiD251Q9~!O60u9BX0zh?|Qm(m#%oYC{C4$jXf(<*7SHv95E7jJyZ4$ zFtPg>46({#y$fMRY&*T{l>^(5=s!GlQ9a)Ij}9$;F#XnIFw11@58FA4J47~uxM{bM zMR*HSzkGQv))zb?V`bQ=?l&|O7D#4R&G&g#N0qgTUg?LQ$J5t3i_}oo@a+i1@o2S-Z)=xnDJ>Qr@gQNZ7tG2o`#HYuj=idY((cOVExAKt~U%z)$qwzKVrm@b42E zajmamO-C60-*h4!%ZglK%~!ic4t>azZ}3jbZz#5RYpLD#7O8O&42Ns+KbHRyE)_r) z-(?1Y3t&E!^%?t&@?%s7we$6Igyn9l6c(smxo3kx_bZzPXGBi@$6q- z4y|&rGoawQ((DKMeCfp8f13C`e;M+odoh$*&%q zrJ$kfF zmojmP7S6hWK)m<&#*w%h#Na|IOX!Y!-HodW>gKA3>VueK)=OP%oO8s{KQza=F|ZM^ z1b46Grj*ri9>wKrWD1RjYk8YZs0iZ`bAEKZ{ts;w$gZOXw3QLZBc7aH3h)bjS{n7I zJeT{!#hm&YmNyMFUHJb%i2CpC;V)9_iW1+9j^^Y?_a1B@>Y*0~g;r+K`EJ79fJ;us zN)%I*4_|cu-ti-B-8pDVypPmFrzNF2!`zDN-}Mlnn^qS1yv-bx`=#BGKFg<)utXZ4 zn^?Y$$X`?MT2A+`~YLD4P|2YASx`^3rTg23|N$LT=UJfK)-+ z^PlUGIU;4CLo8ugDs$Y1*f)Lrt+_;5vZK;IlkolXh5nmuIaY|b9Q2W?a7sZ9av*pD zS3a=~(4i5@SN~5%IX6O!?@Ls zvoG;yuZah7{dQBRqpnr_NbXZeX(E}A1RW1ugrUNm4Rk6;f5VDO)?w^3H7n7`_bomW zd+)n@=JkK-Oig)QVI`>6je;A$`c7LXB16fE+1>XvupRfgg8iMW&80<;9+7QJeQI5= ztMGP$r=LrC_AnoB2^vK%J3;za0>6eXe6R0b7+tqTxo-12lJt9aE_~ZqLTXL<{8;`F zg7lwDI>S#rH>oT9y;f{GV!w?GUI=e_M0OyBO+?BXQRQDJ{Up(@w|9#Hzmu942 zpIrudQR4yS?)Ag75;pjyWpOQ`qmoLHojJNBhYxRMb;MlqTnUwSD0w4473!KR09m7|SI8=bbd-((KBmX zawE%K=<*WkRqHRoYy|~nM)Xc!X#afw*a4j=4?o-{{RQQTounbAm*Vu1Ssjr`(5Z5* zzNliNRQg7L)Lx*YQ+F+|$D8Zb;C?$VGn$^3CR!<;PX#Q1cW=pMe9#K0P_!P>hj%Ra z7^mkko!~~n#V;3|rYE$2ZM`w$sWf^%C0@Anr{6PVU|@VluiT=S1M)!9x!SBrXt9rE)Tk#P!aH*eJ0_5tJk@4ZJuu6k{IyS`xi~ z%b6@5)>(URM5v=w=lW}KsGK(QwUL}s*6-qPYYCvx?Vf+17>?64)yOB1(p+FvQb#ZM z1sg^^iR-`ILy6N_otfCSm4F_kn;1(| z8U(Oi@Ql8Crm(bLsb9V%^BW|Rwz%6!w|II;iqtDcd5vvVU2B9YmcWwLs8Jn9ncgBC z*NkONg7}!eQ7MO9=1~CYpuee4BGAL(U=&1)J-_SxtN9G5#TH0oenaiDx#!MfUc8>5 zMxuvYE9Rr{H#{#T)w6ql1RTU}c(*y{smOgHd=J326D}?7e9rkBU~ibu6*o%%J= z5wPk41mv&?9_6gCi5^IxkHA%i*Y3o`?MjDiNQI_ z#%-jIH>Q!X%H6TvgjD@lWjHMAE|_C5)} z^b0FHYZ$g)y1g^L&E2qq+D>}T`PoT=?6dyftD(?u?_~xg5kRxPQluJ83$Ym|F2&77 zc2m1_0b`Vnfa_ZU*iS=8CXv@{Wx{oSuGH>i;Z)ycV*qB?JMI%@7;D>4tH^SKW(w10 zGbT0)>G`!qez9pG8YYq`p~2_5`7uE34c+CQ#O`x!fDbe@i7fkQDHnD|7Vb*bH@DE@ zV8D$6FN!ihZ9CAT=*Tl+8_pFFYyM8JXrFaf(1>-B^8FA>zVw^Xt%w?22b84OSw_sg zeW*Op<$cQxGETi9!C&q^W@z+xwpz4!N~4K`<3LskO02nRHr4uw*f-b5ki73~iCI_N zS*eCSV|!ukZ{i%`~SQiyEDRJ_Mo~ zw2usQ^`O9;@?I8i(@E8I$5dEkLASq9^pKpZGj67Q@<@#9Y@?3Ogb!!3yEbNE+Wm1*0p@0{0K{J&FwszPGYNE>G7$%n_M#^&Rm>Ve&LA zOOEiMI)~gWzG6;8;@-u(dDPbT{`t{woV1!$ZM&vpAWm0caXZ0lYC95_9QNHmF{*OK z><|llF01yXKeUP>NWDPE)f48=pVN9UR%U0=p@h2%dsOIq4{;<@6c z4M%Cd;!z2%@|IE4HlxqMHB@nN>UYCQ3t+>EmO=X!KfTR^p>5f`5`1ip79hr{PL0$# zpZKQ6wpsP;yIJLcp{bCIRZMlr&?B&x4^f~1Q3*9S*la3*lmdZ zPWdMT_Q9X9f_?|{7O9he_)%4`=3Uczlv#TIA-qfwcD9I!Xl5k+jD;)3*_(y0nVLkd zs#D(CRAD?_eDqz1p(%fdJ+tKSClLXRNnQ~izAJ7^``znwkK&54K&IBGoq8&ZQ?8UV z4oVEdw{ao`MpzREQY0rvWJ+^#%_G=JOlk5GHOV_)xLmx?N=5(U4sx?X3M=fqaTyPs z>%G0?Xz~p;r7jRr@Iq!vuW@z;9$-I`bK~#b1BhDBjAS=vC;=Y z;uV-gDi}f<5eeaLvE3o&i@}^7hJxLV8~b%hS!s)G1=Fhpf|$rQ`s%Lg!XVGv>|MO- zE=rFnD4ChaiV-q&F0hAfQ6njS#%BG?ALEct}(tHB+}zCOJ1eIa|0mp4-g)2twz zXuWGwk968gxb7OO&o>^PY22LE1xZ24#Cz=hKLX1c#L1YtEwO21G3c)+6zCjB^$EOP z{Jn`uG!T*g!vs5tTBN`Wr<-&dWrSr0!%9Zn5`QtQgg)!OUL5c`-OG3HQG(bA2ujgQ z+!sA()S6=hSh$0>-(rSYrZ19-%erHj_V&L0XTV5;D!II4eal;|ttgWz#OrKs0G1eM z2+9wQ%Zx|~r@6#(872gltQ?`)f}2<0PqSQ5@XAn}*2gGm&*N&`2I@6$qGMAaD>vR9 z=dqHa-&sjzC}f?Sgja~?Ik9{2Y(skPMREkK?_jRs|LQi3AWNogadAt9=Fq45;{^!L zNbuSfx5(X;Vqjm9HLC%i#?9q}I|Qlf>ktLEC0w|qBoOkOZ1+MZMb=0{4CYUG>h?lv zqINKeU{GJmy!wcsLs@7?GCCzR#>1aiUAAc)Yz1)1X6d?Id$!ujD7i2W^D|eU^$9PR zi}$0=@fK#92DTMe%o{(e<3d_6vfmh;G|qgJTN9&LsQ_LRE%NJzpC}V9cB zN1}|6eViwF>x7P@eWx{}JN^atK)^7c9^^#OEJ2RQt^-yPk+oq36+ilKXzi}IZ}du4zweO2aIP0Q(gneqH0h7Ral9i2ZJtF&gHbypv^&0>p*|k zVmDcA-ifDRt6K#fG_<|+p|%&+sygIf-+z^+xHD^D8Nxlrq{Cv6F;R- z$q%p<`9_a480zsh!Fn!)M!ax7;;P)2^9MMWDo-v3b{ z!Z^4vtjJHQ3?GT5b+4r*{fnD%$Hb*8#Eofw4fT=! zn~n4=J$~{0>uU#*IbyuSP>B`CKz}!}n7R2}$6dhoByx=O_Z56Z-f`yCGldbY4|*N@ z6hCoUZcF%aPur$nMuAIs>|ZooD0)^{4-N39LCY+21N5JzeorIzrIVH*3Ic=!-}DGffd_PiSg}O3?8(K zy!sJ771Z%%BE+P_ji2#)rq4~@2w3`dUExfyt`QNtdvU~`dI}?6>&=8lsIX{wcjz|d z{@J@&LvEbP!iv&}N-os>)OI|}qXs;1pgb^c6BCzsW590#zogAU6Ydf4>r-4 zGEF8Mcut;;zmy-XT#Tm@bhl3QnY=yBo!{jm@%(f1;V%7~x#0)~P`1YV6%OFuw3+Ev z-J@r8nJba@la|L8(O!ym!fd=WbSb+?j52$b5o5Y7s292VVrHPe77OL-%d7tL4`U8< zDQD)kgdfH;I-Sb?_|_e}Ggu)qPE0==EE&CcnYVPLP^hSzK?Hrd#mfQL6Ex_k85b}ocfdD~?O{GcF>b@!n$rw4>$GC*Ph23GKjTi!+ z%}vD#gL^43?0ONAd;DcsQkL)L3@IUk5#QVhJ6Cq6Sf*q|M!>9|_`qp<+=RP-uMNu^ zd6CW-PBg}lsvP!XW0Qz8>txam(fuMqj_kK0g(~vs*4BOK-c)xKn)38cGt)mAHo{#A zit{H7;Ld&_gm^9aY{5eaA0&azE@82xR^r$ZHw0i_!p{k)rA0AnpQ7fMk7KiRe0Iwm;$pGoVIOw-CmfW^h-0~_Y1q}@8!~R zN1sTPF#I0OO*_bSC^3wAf4&jYF~>84pwW5fI`Nu0G5O3h%?>eyjROPiX^*dyEg^p% zt$;p!_{InTI8gmvK%+S`GfPe2gdDget|(|RC5yg)cD0QWlf^VY(9&`F58Jx|JV|1+ z?Zm#%RqS}ugtU$-&nZJPOU0Pdmw07gT!FZ+gfU^(IhkN_lAXKOLSg@_owL!%AAq94@wFzLB}#wJ?|rTT zMcZY(d(?|{8O$=6)IfQ2THt48bhx_>>x8__A{cIqC!F;>@bU(5A-8urAO*u1{G}vSt`y zK%EI}rQ~#Fm4RZRrG@WB^T2vzPug7O^ zInW=#k3>*kA@#Lb*N_o+zcf7H+R+y7QoCvi5T&@;%vU&zGyzKcM52f`E3#Z9@hYH^YY?FJge?4 z_?e0u9cdU7NSsPOux7boI9Na;M^&u(Hu_Im-acAwd3D{c+Lp#9B5^5yyEu!t@I|8hn4a z;K}qh15ZUOTh(bbe47rMis+7q&f%t|uC)dWl+-6{IuA%$djb%W#=* z(At(L9OT|brh8MgL)?vCIVnSp_~;V^=6?aMRzf|~rw?q5^{2#!Wfy)wjdGB)vY}&a zt$!I$8Ax5O*jG?JiP6fs`pfie&r4gnLB9K4`R(q)3T%du42o_-LSHc|k2hUE5_8+x ze!4c*!N*7*0o6k0Y}4YT80a$f-~AP!qTPPugH8MuD{X~FH#HGh*GyM_*w>6?5dE%z zScZ#AH{e6x?l{c+;w*mc!U0w@ub8J&EaR}8ZSR$K4O;~R=r>U=%2QSbB}-`u{}W|H3Rvv_>v46 zDGFWe+m-_Cctp~C6T`l&Zi3Hv<_6najw*`%arZ_tQ0Y1m{bw!V%c^8f!iwDcL~;(N zs^^Tu{^I;rO=YE7<4x%oA?>BNO25tpO60v^7)PdGJX|7ZmyE<{ZgDLAf(oZ2!B=ju z#j-j5;zx#G-(;Tpks0glV6*UAQkK=#%lR}ntJn9_iht+C{6zo^426?|w1gJwl%END zPAG-lcPb8b1yOHxxVHf*I7=$8eErk-v(6SeKfc;jw@3o>CTF-|?>cBw`@Y$su9N^c z6}5%Haz>V!GijaVCQB+6e{pOuoxnaVHo(-q`qY$X$Vc!s#Z9Vkdq@Y8*hzsM$8dW8mb(nWZh$6aygu*-MIwwZrA-6 zxCJWqc_OxF8|%UpUNi7iZ9iKmvKcH-A@&#qPL35rPMwSl2|c&C;KXW2j@0lly}T#w z@gD^Y;7uxw7Biq(o}e{~^Np#m*8!gHXcSEZ|Ev(}GuCc}<;tTX&U2`%l~2oK3$@w! zqqH_J)n0@GatG->$JI9|^IyhwR#H{=0RQh@9&d_;@=fd7p}lPeMRwMb`2+_a>c{oU zX@`bu-O+>)PEsjz1XKBMwpi@JPQo>kz6D05{)nE${r0GOD@AYh6}z5YlaIR#@g5!Y zJz!D+ToNaHB)>>+Fi^p+r<7@oPx(cle&Mj&WzAelrz!`UeSj=y%DId^hM6Y~= zCZWThJnWeAzl$?ioYhLw{x5lOmmYE~_`W;#@ZtImxl-KJ!t9h$RC3h1tENy-!RP^` z&+bZxQE7Y*=h4`RH5P~&2vt1u`)>!o6~e5O_M7&9M~pdSpnoO7vfE+R)C!$C&d+ID z@-lcH>d&*I%s5S8tg;DHah^k5CmwSLkkHw@;ui_PEn#e=8;c~)eB7MLg3WU?v4b{F zA*-XS9ySCFCBLvWf?=AQ3t!_L!fct!2>cOl(bfYIxD0gc1dZCnnH6>EsSY4OtWCXi zQsr(CT-+Ue0G%kMsI|Z>6dSSIJeJ!VcUp;`A^)KorQl(iY@W6ctYZF6eFzUpC&=Rz z-{T7d#U#NbP~)m$h7ANJ7`nbvGy5Qop18f!8!Kk-^b~E~{;B#xzXW!jySK|Y13r(w zaaVdUkaRMczc@%l{9)vMJxE9O*ZyUC`2vJA&u0IMmOzOr;T$z<_m>u@;o>ntSRlZ~ z$F7qu5=lC%TcXTrf%*hOJf!MWSM$gS+{mFK)}>oOQ~#^e|Bf~G+KYb&5C*&hBXZS5 z1E+sfklAD94!9CHkyOs&K~JnXum^_QAO70GssZ*_v02a=5*na+mfn$k=d@^Yn@bn1)q&es%#OQ&7X9-`Y3( zti6V%_F-ZUgf9%pzVd%?zhPC<**zEzb@s?X7<%+8zeI|Cd-0fX$}lus9+Ddbcx98CsUE2NlgC#;G02It94>h+ z=mKln&0KT3x0NNI)4C_ff=n5$1#0}gD6H!!7WYbLx0za&YVW?v&{9Er|YC}Jz!DlDfTG3tz ziG7WVXDciO0xZ!-0XxIv{9b=1RDTv`#mBvP_U~{;AFW6ri=zbrZ}N6N<`iymROkfI zovAno5=Ucp(zmXLI5<=XZ)I`0_?N&ubh~L~4>ShgcL&j{+q0Cbe4mfr#n%u2kQ~ z_dg5y$4LuGYkw?VM;a6&*Eesbo9}$mMsYIK-EgGq{NV>^vt@GDHa3;U%jM~L%X6Rj z^cZ4R(t;pTdOE1%9n-diSq#dOO4LvWEcb7dkz&DZmJ1YEK#v{Ba4kFMilp0tw}``A zCIS9&^Jjh62Sq+s@gg+9AQ%(;R#}GS?tsU)jtXg*XK5E0F2JzZn`>4F^n~Mvfv8wZ zm)atv*SWvDpJMU@7-&<>Ep84kTc>?FnKV-hd)_r@Jp21LmrMY4Br-|l{4-OKm8jy{ zuEq*Xez#f-olg}Ai23RY(NR6-v9s8O_#McU+hP+!zeySHoON-ZzIE|gXo!9E@FIQ= z04-E$K`!(GYWe3S#SXkfAvJ+^(R!Lj5iX=$o1Xe0uKq&~WN1_a3RR*wvAaI;(%V~Bk6e_3i3N$(pIY&s};Zg&_p`KdNkX^{xDdD{V z=Xb`4^H)0%{mrSRofN?y((EZPfLZ(InEhnrz4VxtLszriervG!1`N(D!_>;)5C1sY zS9yzyNS8$dg4F517bK8XKdP$qC)?3X-nUY%78I!fbxE8cF@}b~7Q~~AXH2DB<7!cq zK^=&e`M=JtJD%$H{U6E5O!iF2rtI*bLUwkAWM@S7-YVI9uY<@I83~a+63IRxdxx@D z;(MLv`~051fB*d99545|KlgQy>%PYOeH}XR4#!MVifL?eP2Kej9MkB2+Wj=5^ackC zK!gnyMNAL22qJp>nRyJKxqLXMG$8nxVlIoa8dM9QXX)x~3azSUJZ$;YMrm*GW+t0i zZ#IPsHgMx#V<-w9!fE)B=zp%w*lDKCK=p99%YsN<=aEem6(A|iI=L>ECY%&BTfVD< z?bntXWQZ!<=xD6`;YVCei|g@N+N(cNztua9{GtF6CZ3x^aUHHY!uV@0jBBvxHHPdu z$aT!53YN+WmD-6QeH>DBJ7=qAuezPH&Vh|fF9>fr&vRJdhQW0Btv9ciN;R~OEj`}T zIeQT8aAuop`3YmVm!m@q0vE&FfslE!_nzTozl30l#-O>>ZGDZzAVmI3(8!N%tEa%h6u_bY zy%y^rbL_2b{@!L=>ujNn|0AOJg}-;IBinZszHK!wiLVa?m7}wZxEmU>qi9|@TmfL= zT4HIk$Y7o+F5Gui={=SMSZxSQG!HngzI1o>q5b1L8op8v4*UIChYE~-jt%dmOcW$z zP#azQ@DlJzXS~X!$_x7)_~Y?G;h~;-4YiUr(*tQz&OnDD_0u1HgLM}}D(3ytw%r{{ zcBO(+fT)a)JXq9@XQf7JK8nw#i`?jb|GCcyaF0BU?b2o>h66I{E%We4`LF9kZf0(g z(z8Jz-Um$9ECOKx?hBcENs+X2t$gB#&z0FbA+1k7>uWxVdR0{7(JO38-6g|D7#%5W zwMaN7Xz4O*dWZhdS^4Z+gh1%bO+u}>t)To=Z90m$&%W3Oo@1C4E zR0g1hSSdF1&hh31wHU^yzicnz+1}cGzKfN|Hj!~Pf~w#PP8!5F3~RYaVJYciqKDyn zYiK}DFYh#sE%{dQf1~J02YX8O5u%rdE*x#obMf1^#->yWmEd!k0%f8#!)LEw=}*}#t@rL%NHt@O zzx0=^QRNkgQG8>wm^oTl>v`wuTHra~Z1m=2GDvtfcW+wZ8rN*c+<+rZ`o?9~9H~$xp5B+@_=0=PiQN}ecb}?pi4nm7L;9c8)vl7 z(+gNh3U%kN0Tk=oEm7h(-%!yqw<5zWJGg+KeFO*~WI|IQQ@R}wi=STkUTHx?3|sz z&fYN&k)^4vw;Bp+pDD@GE8lwIFY$~B8$&7wLJPl348tZ3Zmp#@H`DRk1n1FufS}p$ z*UT+(g(ALfm9i=e3V)0sENvZmNEJOu9$6%swF2qL&|z#9W5o`lOK$~1J88t^%1T!f z$Z>IhWHlLc&OOA!Z4THuc{J4&H#3BT5DM;;A(v+ng{|NTJRKW#KqZ55S#FkKsG*t% zuFGA*oQjQe>dcL~3SpK|yrI@n3t7fQ9}TE4Oah@#B<0-PmL3po??0Ls_{)es=1UUr zdI1BX3jhYuUk@v`Wbq=X&en0*khntOj~q8%oX^E?t_n37+a#5F2lXmk)P=wZz}3O` zN6{pn2k69(qT!kgcS9$BS?F`34-Dd@aT+JY(+4cKl(=foSr`OU+McH*Lgap52w7IL zSiN(?5=;3;O66iA{^?8fcUdd)v3NlKbd*GhrZ5!FS#t{gX!T)U;UCu=I-PbLGO@BY z)##CefSsn=<`df3Y$lCoEf8O=4fFOFVST`bA8-I3k6YHy)-5^~HFUbqz85GzQ0M#t z>fuU91BwLOH=ChY%$wJQ7&E2e9eQnv5&J-K6qHIYcNFEYyCgD9*cU3OV8An|rPxP! z0T^P};x0ur0X<;<87+wqP3?u~h_W@Y4I*TFjolotCaEj2(y}1kc$feoH0RVwC;7~|Y zi@L(a9gP+v2$adz0RBi@M`fW+#moG*upqSqzJPA5_z|3vHk^`4oN5`ONPVJ2eWDNY z&(}Q`{RiWZ?;H|S@g~X1h2t!zt2|zRo834R##@E0T~7!pO|1*q)hZDHH;8+*2$0Jh zPPs_kY6k)#>Pwn{T)AiSMxV?g}4W;lr}Gg6h>=(q)jYp8A30~+Mr`FL#Oz!1h+WmU^c-v7s>kG z&T&fwN7iqR3-sV&B6TR=RYHcsl}-p$iRq4E7ew2QPhEsS6aZtABEoix(`vyEwYRlp z?yS)AQ7B26*?YKi8|j$;Z{HNkk3+s9RWE~ZB_uT_BeU;fJexRkk@})numE`4)LxB> zhD;Vf;pQ-SE&+M9(8PH_nWPvnOxWuU1-^Be*n~A+=cyOYOcZpWL7drJQALvE{ri9f zY5A{IRzMMIFV169aF}4Y3YpL58oY`Ae%dt0oCO-l=)}(L&8C|_i+7a=npP)-58t}{ zGAzeL9knh-*0q@mNZy_w;f(y;GWY9tDbNGM)!8#BzSR23oU^OD$~9ag{b#Ncn91MC z(N#?bLd_%0EChMPHuY6DSJ1JHUbgRQ9JCXA&adVo5Xi+161G_8$$D1NdiT7Ct=#j9 z#vZphhJw<9y3170Bmy0ELP8x5(9Er|ObfaaMq^N*9>?Q0(sD?x|63Uf0yeQ1ws2o| z1>ZQaqLsIlgL(+vh6-t}0QvG7e+v7&kh#4aCk_BTc4?>(0Px;CW5kWehNOxZ7->U$ z>?nJ|mRYkFeksi`IzU-sBsRO8;{cU};Am}vpCnAc>d74!idF`CtdF!TU;zKw7Vz-H z7p}3mh$wA(41rk&5|85-wsiG#f>q{xDRki=)%VeELPGRYS}*=mx)$i&$oNYF;u6-) zN;iSOWv&|e1$<)p}H%(ePQD7JQ4mDo=3NP7?N!{~Ra&X^#`3EsHAfH;X) zkH_wcGED4&;p+{var3>{%7fng^yVFCdl z0M6y)hM<+nsM&6tR1zqs9?U=dJ?l{URcADEuWQ2b^*GJ&7VjtjR8jq4AQc22i~Jc+ zAQFMhhe?L4946#xNw1o$A0I{w+^H6rXNtUHy!GF6#3a;&;h%9rQ=Z9uKYJct0950G z?Km}oK`Ao!v(vZM2bB?hkv3X z)Q>&-Lbnf)xx=P4XR?h)lHz~s`~UWcOb)p=h>*gR?((_+V<;qAww;+nGZfAtTw|IE zIol;$mhQ@xYOz^UEC}=T8RChnx)q0o@`Hcvu{3V8 zGJH?^_8;X)AsCPYQ;bPj92PT_PTmWZxxfU4k0h(E;%O>dUsCy>rhcdX-tOcq?hl$B zPel2C{b;PdW8@Zok5jA0;o>4&h%t|G4yOQn+GEokJQe<}AP#a~cKES26+^x2-U*w;fBibFazBnA}3pk>hY*`;|3z8Af#}@K}@Jt!{uN8TG?!JyCIw%AJA;x1%X%%PYXvKExa3g#D2EWcO#c25& zFJo^;Ry}{l$f(B^A_U^;T(v>2Qrxi<%dIJ2vFIOJyj}!P`l-tpx;G>g;F~!uZkc0m9y#nSp@tuN>FQE!rqdo42!JkFVA@FX}{f}#2 z1n{gUA&NU97S!L5R&NLI43}`gJd%%HzloBtW4aTaC5nYWIAn$j|D21Z&UZd&Q>d(c ztF=o(FpR$6U1QJG6FCw5&$HMDQGsYn>!SheH;WyifjI>lPoLBfFGQtxR-R%r%Be#~ z#>C`SeOj+#*ZussZ~5CG39A=#mdgo(eia{e5~uI1GPeEQ>?LuC@L9p?k(b|AaOjEO znajXkkF12ab%x7*uqH#(ce%VsVBoegxt`JcOesB+(G1<~Dy1u2D1`};{e^vlS7QRB zu-MRMJVMc##Ck9B+gaCq-C629FO{8Ro9~Fv?_<@vqit8N#4tl_1VT#Ul2uW$2hAwmA+t|;T5eIKizVn#A^5*Ab?lO_ z%go{Gw}agZ$sw!xgIt*e#zN?M;z;I!xcQ zcBLkIlRwp`Z1>(r%;^31=MlT^6o0pzN>Z!}%U@oDR`g>Z?76#^lfp8q@@<*E_b$I- z=)R`&`i{Cn-O69k>1LU1+BspmZs#`{pMLIth~+T)E7hD&CWCrJnG0=4|FLXqR^Qdp z|I^eUE&}l~Fhm$f;$}@kaAwvg7CseAr+3=Bs3eyGMR&V8Qv`wlsw=SfHE zy62)oZ8=PDO_*Abh%xn!6 zi5`V(`WhtKFs}B=u}s#v+*Q9RIG%iMtsdD8pX`VEyk6Ki_;DarFG?k=JGMy3{xr*P zXJ@!5@BW15wPdFjDFqP{H;bW%{?s2RE#Ph0IA#Ms#zCs4KhU!?bIZ;wI?C(qMGMvI zei<(Eg46idhqrM@}qxqQ>xtPBmpFCKlR+3Yb JEt4^R{vWz1e7FDr literal 0 HcmV?d00001 diff --git a/src/aind_ephys_rig_qc/parameters.json b/src/aind_ephys_rig_qc/parameters.json new file mode 100644 index 0000000..98376b4 --- /dev/null +++ b/src/aind_ephys_rig_qc/parameters.json @@ -0,0 +1,5 @@ +{ + "report_name" : "QC.pdf", + "align_timestamps_to" : "local", + "original_timestamp_filename" : "original_timestamps.npy" +} \ No newline at end of file diff --git a/src/aind_ephys_rig_qc/pdf_utils.py b/src/aind_ephys_rig_qc/pdf_utils.py new file mode 100644 index 0000000..ecf67ac --- /dev/null +++ b/src/aind_ephys_rig_qc/pdf_utils.py @@ -0,0 +1,106 @@ +""" +Helper class for generating PDF reports (using the fpdf2 library) +""" + +import os +import platform + +import matplotlib.pyplot as plt +import numpy as np +from fpdf import FPDF +from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas +from PIL import Image + + +class PdfReport(FPDF): + """ + Sub-class of FPDF with additional methods for embedding images and tables + """ + + def __init__(self, title="Running title"): + """ + Initialize the PDF report with a title + """ + super().__init__() + self.title = title + self.set_matplotlib_defaults() + + def header(self): + """ + Add a header with the Neural Dynamics logo and a running title + """ + self.image( + os.path.join(os.path.dirname(__file__), "images", "aind-logo.png"), + x=150, + y=10, + w=50, + ) + self.set_font("courier", "", 9) + self.cell(30, 15, self.title, align="L") + self.ln(20) + + def footer(self): + """ + Add a footer with the page number + """ + self.set_y(-15) + self.set_font("helvetica", "I", 8) + self.cell(0, 10, f"Page {self.page_no()}/{{nb}}", align="C") + + def set_matplotlib_defaults(self): + """ + Set the default matplotlib parameters for the report + """ + plt.rcParams["font.family"] = "sans-serif" + + if platform.system() == "Linux": + plt.rcParams["font.sans-serif"] = ["Nimbus Sans"] + else: + plt.rcParams["font.sans-serif"] = ["Arial"] + + def embed_figure(self, fig, width=190): + """ + Convert a matplotlib figure to an image and embed it in the PDF + + Parameters + ---------- + fig : matplotlib.figure.Figure + The figure to embed + width : int + The width of the image in the PDF + """ + + canvas = FigureCanvas(fig) + canvas.draw() + img = Image.fromarray(np.asarray(canvas.buffer_rgba())) + self.image(img, w=width) + + def embed_table(self, df, width=190): + """ + Create a table out of a pandas DataFrame and embed it in the PDF + + Parameters + ---------- + df : pandas.DataFrame + The table to embed + width : int + The width of the image in the PDF + """ + + DF = df.map(str) # convert all elements to string + DATA = [ + list(DF) + ] + DF.values.tolist() # Combine columns and rows in one list + + with self.table( + borders_layout="SINGLE_TOP_LINE", + cell_fill_color=240, + cell_fill_mode="ROWS", + line_height=self.font_size * 2, + text_align="CENTER", + width=width, + ) as table: + for data_row in DATA: + row = table.row() + for datum in data_row: + row.cell(datum) diff --git a/src/aind_ephys_rig_qc/qc_figures.py b/src/aind_ephys_rig_qc/qc_figures.py new file mode 100644 index 0000000..d3b7fb2 --- /dev/null +++ b/src/aind_ephys_rig_qc/qc_figures.py @@ -0,0 +1,71 @@ +""" +Generates figures for checking ephys data quality +""" + +from matplotlib.figure import Figure +from scipy.signal import welch + + +def plot_raw_data(data, sample_rate, stream_name): + """ + Plot a snippet of raw data as an image + + Parameters + ---------- + data : np.ndarray + The data to plot (samples x channels) + sample_rate : float + The sample rate of the data + stream_name : str + The name of the stream + + Returns + ------- + matplotlib.figure.Figure + Figure object containing the plot + """ + + fig = Figure(figsize=(10, 4)) + ax = fig.subplots() + + ax.imshow(data[:1000, :].T, aspect="auto", cmap="RdBu") + ax.set_title(f"{stream_name} Raw Data") + ax.set_xlabel("Samples") + ax.set_ylabel("Channels") + + return fig + + +def plot_power_spectrum(data, sample_rate, stream_name): + """ + Plot the power spectrum of the data + + Parameters + ---------- + data : np.ndarray + The data to plot (samples x channels) + sample_rate : float + The sample rate of the data + stream_name : str + The name of the stream + + Returns + ------- + matplotlib.figure.Figure + Figure object containing the plot + """ + + fig = Figure(figsize=(10, 4)) + ax = fig.subplots() + + subset = data[:1000, :] + + for i in range(subset.shape[1]): + f, p = welch(subset[:, i], fs=sample_rate) + ax.plot(f, p) + + ax.set_title(f"{stream_name} PSD") + ax.set_xlabel("Frequency") + ax.set_ylabel("Power") + + return fig diff --git a/src/aind_ephys_rig_qc/temporal_alignment.py b/src/aind_ephys_rig_qc/temporal_alignment.py new file mode 100644 index 0000000..fa4ef1d --- /dev/null +++ b/src/aind_ephys_rig_qc/temporal_alignment.py @@ -0,0 +1,288 @@ +""" +Aligns timestamps across multiple data streams +""" + +import json +import os +import sys + +import numpy as np + + +def align_timestamps( + directory, + align_timestamps_to="local", + original_timestamp_filename="original_timestamps.npy", +): + """ + Aligns timestamps across multiple streams + + Parameters + ---------- + directory : str + The path to the Open Ephys data directory + align_timestamps_to : str + The type of alignment to perform + Option 1: 'local' (default) + Option 2: 'harp' (extract Harp timestamps from the NIDAQ stream) + original_timestamp_filename : str + The name of the file for archiving the original timestamps + """ + + return None + + +def get_num_barcodes(harp_events, delta_time=0.5): + """ + Returns the number of barcodes + + Parameter + --------- + harp_events : pd.DataFrame + Events dataframe from open ephys tools + delta_time : float + The time difference between barcodes + + Returns + ------- + numbarcodes : int + The number of barcodes in the recordings + """ + (splits,) = np.where(np.diff(harp_events.timestamp) > delta_time) + return len(splits) - 1 + + +def get_barcode(harp_events, index, delta_time=0.5): + """ + Returns a subset of original DataFrame corresponding to a specific + barcode + + Parameter + --------- + harp_events : pd.DataFrame + Events dataframe from open ephys tools + index : int + The index of the barcode being requested + delta_time : float + The time difference between barcodes + + Returns + ------- + sample_numbers : np.array + Array of integer sample numbers for each barcode event + states : np.array + Array of states (1 or 0) for each barcode event + + """ + (splits,) = np.where(np.diff(harp_events.timestamp) > delta_time) + + barcode = harp_events.iloc[splits[index] + 1:splits[index + 1] + 1] + + return barcode.sample_number.values, barcode.state.values + + +def convert_barcode_to_time( + sample_numbers, states, baud_rate=1000.0, sample_rate=30000.0 +): + """ + Converts event sample numbers and states to + a Harp timestamp in seconds. + + Harp timestamp is encoded as 32 bits, with + the least significant bit coming first, and 2 bits + between each byte. + """ + + samples_per_bit = int(sample_rate / baud_rate) + middle_sample = int(samples_per_bit / 2) + + intervals = np.diff(sample_numbers) + + barcode = np.concatenate( + [ + np.ones((count,)) * state + for state, count in zip(states[:-1], intervals) + ] + ).astype("int") + + val = np.concatenate( + [ + np.arange( + samples_per_bit + middle_sample + samples_per_bit * 10 * i, + samples_per_bit * 10 * i + - middle_sample + + samples_per_bit * 10, + samples_per_bit, + ) + for i in range(4) + ] + ) + s = np.flip(barcode[val]) + harp_time = s.dot(2 ** np.arange(s.size)[::-1]) + + return harp_time + + +def rescale_times_linear(times_ephys, harp_events, delta_time=0.5): + """ + Applies a linear rescaling to the ephys timestamps + based on the HARP timestamps. + + Parameters + ---------- + times_ephys : np.array + Array with the ephys timestamps + harp_events : pd.DataFrame + Events dataframe from open ephys tools + delta_time : float + The time difference between barcodes + + Returns + ------- + new_times : np.array + Rescaled ephys timestamps + """ + splits = np.where(np.diff(harp_events.timestamp) > delta_time)[0] + last_index = len(splits) - 2 + + t1_ephys = ( + harp_events.iloc[splits[0] + 1:splits[1] + 1].iloc[0].timestamp + ) + t2_ephys = ( + harp_events.iloc[splits[last_index] + 1:splits[last_index + 1] + 1] + .iloc[0] + .timestamp + ) + + sample_numbers, states = get_barcode(harp_events, 0) + t1_harp = convert_barcode_to_time(sample_numbers, states) + sample_numbers, states = get_barcode(harp_events, last_index) + t2_harp = convert_barcode_to_time(sample_numbers, states) + + new_times = np.copy(times_ephys) + scaling = (t2_harp - t1_harp) / (t2_ephys - t1_ephys) + new_times -= t1_ephys + new_times *= scaling + new_times += t1_harp + + return new_times + + +def compute_ephys_harp_times( + times_ephys, + harp_events, + fs=30_000, + delta_time=0.5, + wrong_decoded_delta_time=2, +): + """ + Computes ephys timestamps assuming that they are uniformly samples + between barcodes. The times_ephys are only used to get the + number of samples. + + Parameters + ---------- + times_ephys : np.array + Array with the ephys timestamps + harp_events : pd.DataFrame + Events dataframe from open ephys tools + delta_time : float + The time difference between barcodes + wrong_decoded_delta_time : float + The time difference threshold between barcodes to detect a wrong + decoding and fit the barcode time + + Returns + ------- + new_times : np.array + Rescaled ephys timestamps + """ + sampling_period = 1 / fs + + # compute all barcode times + num_harp_events = get_num_barcodes(harp_events, delta_time=delta_time) + timestamps_harp = np.zeros(num_harp_events, dtype="float64") + for i in range(num_harp_events): + sample_numbers, states = get_barcode( + harp_events, i, delta_time=delta_time + ) + barcode_harp_time = convert_barcode_to_time(sample_numbers, states) + timestamps_harp[i] = barcode_harp_time + + # fix any wrong decoding + (wrong_decoded_idxs,) = np.where( + (np.diff(timestamps_harp)) > wrong_decoded_delta_time + ) + print(f"Found {len(wrong_decoded_idxs)} badly aligned timestamps") + timestamps_harp[wrong_decoded_idxs + 1] + for idx in wrong_decoded_idxs + 1: + new_ts = ( + timestamps_harp[idx - 1] + + (timestamps_harp[idx + 1] - timestamps_harp[idx - 1]) / 2 + ) + timestamps_harp[idx] = new_ts + + # get indices of harp clock in ephys timestamps + (splits,) = np.where(np.diff(harp_events.timestamp) > delta_time) + harp_clock_indices = np.searchsorted( + times_ephys, harp_events.timestamp.values[splits[:-1] + 1] + ) + + # compute new ephys times + times_ephys_aligned = np.zeros_like(times_ephys, dtype="float64") + + # here we use the BARCODE as a metronome and assume that ephys + # is uniformly sampled in between HARP beats + for i, (t_harp, harp_idx) in enumerate( + zip(timestamps_harp, harp_clock_indices) + ): + if i == 0: + first_sample = 0 + num_samples = harp_idx + # fill in other chunks + sample_indices = np.arange(num_samples) + t_start = t_harp - len(sample_indices) / fs + else: + first_sample = harp_clock_indices[i - 1] + num_samples = harp_idx - harp_clock_indices[i - 1] + t_start = timestamps_harp[i - 1] + + # fill in other chunks + sample_indices = np.arange(num_samples) + if i != 0: + t_start = timestamps_harp[i - 1] + else: + t_start = t_harp - len(sample_indices) / fs + t_stop = t_harp - sampling_period + times_ephys_aligned[first_sample:harp_idx] = np.linspace( + t_start, t_stop, num_samples + ) + + # take care of last chunk + num_samples = len(times_ephys_aligned) - harp_idx + t_start = times_ephys_aligned[harp_idx - 1] + sampling_period + t_stop = t_start + num_samples / fs + times_ephys_aligned[harp_idx:] = np.linspace(t_start, t_stop, num_samples) + + return times_ephys_aligned + + +if __name__ == "__main__": + + if len(sys.argv) != 3: + print("Two input arguments are required:") + print(" 1. A data directory") + print(" 2. A JSON parameters file") + else: + with open( + sys.argv[2], + "r", + ) as f: + parameters = json.load(f) + + directory = sys.argv[1] + + if not os.path.exists(directory): + raise ValueError(f"Data directory {directory} does not exist.") + + align_timestamps(directory, **parameters) diff --git a/tests/test_example.py b/tests/test_example.py deleted file mode 100644 index 06e9e0d..0000000 --- a/tests/test_example.py +++ /dev/null @@ -1,16 +0,0 @@ -"""Example test template.""" - -import unittest - - -class ExampleTest(unittest.TestCase): - """Example Test Class""" - - def test_assert_example(self): - """Example of how to test the truth of a statement.""" - - self.assertTrue(1 == 1) - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/test_generate_report.py b/tests/test_generate_report.py new file mode 100644 index 0000000..42e5c70 --- /dev/null +++ b/tests/test_generate_report.py @@ -0,0 +1,20 @@ +"""Example test template.""" + +import os +import unittest + +from aind_ephys_rig_qc.generate_report import generate_qc_report + + +class TestGenerateReport(unittest.TestCase): + """Example Test Class""" + + def test_generate_report(self): + """Example of how to test the truth of a statement.""" + + generate_qc_report("", "qc.pdf") + self.assertTrue(os.path.exists("qc.pdf")) + + +if __name__ == "__main__": + unittest.main()