From d724716413b11d2d748019b63ec14232a362383e Mon Sep 17 00:00:00 2001 From: xingwei bao <2122533937@qq.com> Date: Tue, 29 Oct 2024 19:44:44 +0800 Subject: [PATCH] add new model scGen to mindspore --- scGen/README.md | 146 ++++++++++++++++++++++++++++++++++ scGen/batch_correction.py | 65 +++++++++++++++ scGen/batch_effect_removal.py | 43 ++++++++++ scGen/demo/AnnData.png | Bin 0 -> 44934 bytes scGen/models/vae.py | 63 +++++++++++++++ scGen/train.py | 37 +++++++++ scGen/train_model.py | 36 +++++++++ scGen/utils/data_utils.py | 12 +++ scGen/utils/model_utils.py | 31 ++++++++ 9 files changed, 433 insertions(+) create mode 100644 scGen/README.md create mode 100644 scGen/batch_correction.py create mode 100644 scGen/batch_effect_removal.py create mode 100644 scGen/demo/AnnData.png create mode 100644 scGen/models/vae.py create mode 100644 scGen/train.py create mode 100644 scGen/train_model.py create mode 100644 scGen/utils/data_utils.py create mode 100644 scGen/utils/model_utils.py diff --git a/scGen/README.md b/scGen/README.md new file mode 100644 index 000000000..fb4845b80 --- /dev/null +++ b/scGen/README.md @@ -0,0 +1,146 @@ +# 目录 + +- [目录](#目录) +- [scGen描述](#scGen描述) +- [模型架构](#模型架构) +- [数据集](#数据集) +- [环境要求](#环境要求) +- [快速入门](#快速入门) +- [脚本说明](#脚本说明) + - [脚本和样例代码](#脚本和样例代码) + - [脚本参数](#脚本参数) + - [训练过程](#训练过程) + - [推理过程](#推理过程) + - [用法](#用法) + - [结果](#结果) +- [随机情况说明](#随机情况说明) +- [ModelZoo主页](#modelzoo主页) + +# scGen描述 + +scGen是一种生成模型,用于预测不同细胞类型、研究和物种中的单细胞扰动响应(发表于《Nature Methods》,2019年)。 + +# 模型架构 + +scGEN 的模型架构基于变分自编码器(VAE),包括一个编码器和一个解码器。编码器将单细胞基因表达数据映射到一个潜在空间,生成细胞的隐变量表示;解码器则从潜在空间生成新的基因表达数据。SCGEN 通过在潜在空间中施加条件编码,实现对不同条件下细胞状态的生成与转换。 + +# 数据集 + +使用的数据集:[pancreas](https://drive.google.com/drive/folders/1v3qySFECxtqWLRhRTSbfQDFqdUCAXql3) + +- 名称: Panglao scRNA-seq数据集 +- 格式: H5AD文件 +- 路径: ./data/pancreas.h5ad +- 数据大小: 176MB + +支持的数据集:[pancreas] 或者与 AnnData 格式相同的数据集 + +- 目录结构如下,由用户定义目录和文件的名称 + +![image](demo/predict-demo.jpg) + +- 如果用户需要自定义数据集,则需要将数据集格式转化为AnnData数据格式。 + +# 环境要求 + +- 硬件(Ascend) + - 使用Ascend处理器来搭建硬件环境。 +- 框架 + - [MindSpore](https://www.mindspore.cn/install) +- 如需查看详情,请参见如下资源 + - [MindSpore教程](https://www.mindspore.cn/tutorials/zh-CN/master/index.html) + - [MindSpore Python API](https://www.mindspore.cn/docs/zh-CN/master/api_python/mindspore.html) + +# 快速入门 + +- 通过官方网站安装Mindspore后,您可以按照如下步骤进行训练 + +```shell +# 单卡训练 +python train_model.py +``` + +# 脚本说明 + +## 脚本和样例代码 + +```text + |----data + |----demo + |----AnnData.png + |----models + |----vae.py + |----utils + |----data_utils.py + |----model_utils.py + |----batch_correction.py + |----batch_effect_removal.py + |----README_CN.md + |----train.py + |----train_model.py +``` + +## 脚本参数 + +预训练参数(pretrain.py): +```text +--train_path 数据路径 +--model_to_use 模型路径 +--batch_size 批次大小, 默认: 32 +--X_dim 基因表达矩阵的特征维度 +--z_dim 潜在空间维度,默认值:100 +--lr 学习率,默认值:0.001 +--dr_rate dropout,默认值:0.2 +``` + +## 训练过程 + +在Ascend设备上,使用python脚本直接开始训练 + + python命令启动 + + ```shell + # 单卡训练 + python train_model.py + ``` + +```text + + Epoch [1/100], Loss: 981.0034 + Epoch [2/100], Loss: 939.3733 + Epoch [3/100], Loss: 922.7879 + Epoch [4/100], Loss: 913.1795 + Epoch [5/100], Loss: 905.8361 + Epoch [6/100], Loss: 900.6747 + +``` + +## 推理过程 + +**推理前需使用train_model.py文件生成的模型检查点文件。** + +### 用法 + +执行完整的推理脚本如下: + +```shell + +python batch_effect_removal.py + +``` + +### 结果 + +去批次结果保存在/batch_removal_data/1.h5ad中。 + +# 随机情况说明 + +在训练中存在以下随机性来源: + +1. 数据和索引的随机打乱 +2. 潜在空间中的随机噪声生成 +3. 样本生成时的噪声引入 + +# ModelZoo主页 + +请浏览官网[主页](https://gitee.com/mindspore/models)。 \ No newline at end of file diff --git a/scGen/batch_correction.py b/scGen/batch_correction.py new file mode 100644 index 000000000..4e07d56e7 --- /dev/null +++ b/scGen/batch_correction.py @@ -0,0 +1,65 @@ +import numpy as np +import scanpy as sc +import anndata +from utils.model_utils import give_me_latent, reconstruct + +def vector_batch_removal(model, data): + latent_all = give_me_latent(model, data.X) + latent_ann = sc.AnnData(latent_all) + latent_ann.obs["cell_type"] = data.obs["cell_type"].tolist() + latent_ann.obs["batch"] = data.obs["batch"].tolist() + latent_ann.obs["sample"] = data.obs["sample"].tolist() + + unique_cell_types = np.unique(latent_ann.obs["cell_type"]) + shared_anns = [] + not_shared_ann = [] + + for cell_type in unique_cell_types: + temp_cell = latent_ann[latent_ann.obs["cell_type"] == cell_type] + if len(np.unique(temp_cell.obs["batch"])) < 2: + cell_type_ann = latent_ann[latent_ann.obs["cell_type"] == cell_type] + not_shared_ann.append(cell_type_ann) + continue + + print(cell_type) + temp_cell = latent_ann[latent_ann.obs["cell_type"] == cell_type] + batch_list = {} + max_batch = 0 + max_batch_ind = "" + batchs = np.unique(temp_cell.obs["batch"]) + + for i in batchs: + temp = temp_cell[temp_cell.obs["batch"] == i] + if max_batch < len(temp): + max_batch = len(temp) + max_batch_ind = i + batch_list[i] = temp + + max_batch_ann = batch_list[max_batch_ind] + for study in batch_list: + delta = np.average(max_batch_ann.X, axis=0) - np.average(batch_list[study].X, axis=0) + batch_list[study] = batch_list[study].copy() + batch_list[study].X = delta + batch_list[study].X + corrected = anndata.concat(list(batch_list.values())) + shared_anns.append(corrected) + + all_shared_ann = anndata.concat(shared_anns) if shared_anns else sc.AnnData() + all_not_shared_ann = anndata.concat(not_shared_ann) if not_shared_ann else sc.AnnData() + all_corrected_data = anndata.concat([all_shared_ann, all_not_shared_ann]) + + corrected_data = reconstruct(model, all_corrected_data.X, use_data=True) + corrected = sc.AnnData(corrected_data) + corrected.obs["cell_type"] = all_corrected_data.obs["cell_type"].tolist() + corrected.obs["study"] = all_corrected_data.obs["sample"].tolist() + corrected.var_names = data.var_names.tolist() + + if all_shared_ann.n_obs > 0: + corrected_shared_data = reconstruct(model, all_shared_ann.X, use_data=True) + corrected_shared = sc.AnnData(corrected_shared_data) + corrected_shared.obs["cell_type"] = all_shared_ann.obs["cell_type"].tolist() + corrected_shared.obs["study"] = all_shared_ann.obs["sample"].tolist() + corrected_shared.var_names = data.var_names.tolist() + else: + corrected_shared = sc.AnnData() + + return corrected, corrected_shared \ No newline at end of file diff --git a/scGen/batch_effect_removal.py b/scGen/batch_effect_removal.py new file mode 100644 index 000000000..ab4a6b6eb --- /dev/null +++ b/scGen/batch_effect_removal.py @@ -0,0 +1,43 @@ +import mindspore as ms +from models.vae import VAE +from utils.data_utils import load_data +from batch_correction import vector_batch_removal +import anndata as ad + +# 设置MindSpore上下文 +ms.set_context(mode=ms.GRAPH_MODE, device_target="Ascend") + +# 配置参数 +train_path = "./data/pancreas.h5ad" +model_path = "./models/scGen/scgen.pt" +output_path = "./batch_removal_data/1.h5ad" + +def main(): + # 加载数据 + data, _, _ = load_data(train_path) + gex_size = data.X.shape[1] + + # 初始化并加载已训练模型 + model = VAE(input_dim=gex_size, z_dim=100, dr_rate=0.2) + param_dict = ms.load_checkpoint(model_path) + ms.load_param_into_net(model, param_dict) + print("模型加载完毕。") + + # 批次效应去除 + all_data, shared = vector_batch_removal(model, data) + + # 后处理 + top_cell_types = all_data.obs["cell_type"].value_counts().index.tolist()[:7] + if "not applicable" in top_cell_types: + top_cell_types.remove("not applicable") + + all_data.obs["celltype"] = "others" + for cell_type in top_cell_types: + all_data.obs.loc[all_data.obs["cell_type"] == cell_type, "celltype"] = cell_type + + # 保存结果 + all_data.write(output_path) + print(f"scGen batch corrected pancreas has been saved in {output_path}") + +if __name__ == "__main__": + main() diff --git a/scGen/demo/AnnData.png b/scGen/demo/AnnData.png new file mode 100644 index 0000000000000000000000000000000000000000..e6aa27577a0d28323dcadf6deffab18faf0dbf1c GIT binary patch literal 44934 zcmeFZcUV(-*ESl*aTrBlupxZ}6%hfc5kgS`QE4Jgz<`Q?5Q?}!Zq)IO#M1;^GQUigI^V>Mj^F81B{yFDd@B7zzxvm%_W$(TETKB#7iyQjd90&Lg zz+f0N@s-bb---`xlP|F_N6n|d&q-)R^u zlpRsTs#IsS&XaW6)S^9)ApVvqbjEuF*fGP56B z^yfSnb{Q&2B=lVwybIUUD@5vj3`%wgi$@>2RP1@gSmoqFpNkidU;XPyXyV;aCFahZ z(8)yK$%hT~4P(q8qR$r(;gT(X^{|i6=NMMzdtmkGgt)lMlOmiQo`NvJTf@D*@+W(= z>Pt={dZ0^%%{#iwAz`rkw5*4NU5QeA{?_WtFq1|8HXrl*Sr`m)!_>t}U^PzI7Jo1V zgSmw!!R5efW#=~cfX`mN`t~>Q^y<~cU8$XT1WEbg<=>8ir;w-6=Ycf#198J(hYsaI zPbW{__@58{cMIo+OAb3Y29nIpeEDNeOt^b-qbkiu($pny=4Wkrb>w)6Xjs1GbF0KL*AB1RB4?PsHcQa?XPmV-(KUZ0WqahqdX~nIB2muiheccmM-Xb~AJR5Px=y(*(RPku zOY-M`4%ZUh$M+Ve@pH>~QyCxDC+@n9b+P3H?#yT~nxF?b<`*;jSTEjT*!&`7lp#Vt zKgU!qL{@T04l8E_8@5b=OL6sXudgNA42w|)9Nu(WfFr+I|4OgqCrI!JdEFU4e~OF{ z?2;ZDgb_p`6J8CmuDv1+SFdQ#kl3oF*7qS zSr^}J->*pOjSIgev)Pim??2_!~Q}cG>_#Jwb^TTe#W8zKv!rj3uOQwIUFQm#QpA=!br}rz9m|`N`e+aVS z9AZduuBuBG&)e;SH&!z~2b#ms_P-e6>?f*QQoaFvtCK1=G-CLNM8e7-4XscbyZMDT=_TL$Vj#)1O-{1fW~YwFUmF z3Ku@1yJzt0238iBM(J~_Kt!$Q_9gKqHI=O3?d7;-bZcPH=5^oa;{siEXKDv>zB+%N zOp^d&BbDTSt+e5E7??j6vY-KlN>wLnXwn;kJ5qpBFz{83i=r-MU3U7^&*=i*g=g9c z6^Qj>^}7DCL+>Kd(|g`dnJ@2!!AO@#SIV_lrNe?Z>pzM&1#Py3={Ifsd=(^~&Z@bu z_x$d4%J%JYBN+N`=%St2AS*F;xt=Q1?MF|C^!wrzN?9;7uN~JUq1cP|S*R4_>apK< zhYgCz_2^|EF-vask%;1$o_(@2mtWITgER@nNqYusH3wm^($ofC{h;=#-YEg)dwH>af0^IouxMJkTg@EpC7~ zJlUsU;?OiwnqQ0DD9hlNvx99sFVf`!S0&cm(qeIln~g`%?e?$iCV+MmKK{Ww?tSF~sq)@aaaf|*|70+b$p3rO`>zZ(*dyOJV2l49 zBokX(JUFK6B~#P>7|ewWJzg-+I}S!hJ+E|m_7D!~1VOdMp|S}!*d4YD^B7D~O6^f1SIz2s-2vc-uu`K!OL|*9{kl$? zAGRXPp`pgZwb@e@dZoAw)KLE6GKo^qk^I8`_j;&0Q7Qy!8A#pEztVcC=fNWE8@U21uhKKoN%|~&am~4y zXgh+zs6D|l^ycto?k^=QBE42e*XE|e9rMoer;+_{_@aZZ`R?o-PS>HQadZ1_r74x$ zZML)eZyBlO(cenAF7^d)7sj!R1s3Ml-ot|9QJ+002cWbmMltjfq!(*{L?!jqHI^^9 zj#uHxpg~9eNxcv6f|&D$(Zj;o@>t3T?U#`$3}XCVuemhdTS~z{!Y~-{HP&zi<0J4s zvmgfZcOA z$JIhL@bXH&47XuYoaJpp`Q@FJVYSTJH3{s4TE*-f+yK0e9>^_oTB7Pg(yWAth=@eh zhQZM0lNSDhme4Ad!p>J;jVx$E#zdcG9!E3%)ZT(L45mBf6E}NURN!hEj&{>psV>iP z?3A~Be|$nkv4l6a%Yz}Lz)V{Bi3Og0$gXyM^ll-%mji>jFq?kB%9rv;>u7N&%UR2I zsiw|MSdiOoY-#Ux(kQDfChXD6f&PPmRGNu4*TR7R(g~)YG@t*5+RJX$xy9zB%mwyxqz| z=B&t%tS2wiXYxug7;`KUKHpxc&~OY}%&S~sP5Y9`HM^xtqGeDCe8X7!53?OBuhpV~ z{s>jSJ0B^hqY%@qASLgpW|9dt=S;ui^qi)&3EY_{yj& zHqO49$v3jS-X7wPj_ZND)A49wN>9iQQnz=?*NMKsglio@Ois2X*nkTj)*HoPK<{m9gR%( zFLCk>vZ`~7QZ$Z>c(5r|wxSb6zG*Z^?;z3JzCz9bryx#CR8)Ht!IaMt~$wzh_ec`p$mm+>UzK652=(YZ)cT^&(3G2qK=R85LA zx@|~b%pB({ICju=b4NL8-(DDjkccIE59$S*yV+{q-WeKSM^HVQXz9z(w}qM2+ar*C zRYhEq-G-W_D+PifLC)DpCD%B59QW`ov~=MI%86kqc@M4?H*J!_T7IBdrK``m-8Odu zS4v)D9(d54x9IJbOqQC*#zb*o!ctWyi-P$O@StvV;ELusD@(bowGA6Cl7D|+Zd_Jd zS)kpuFZKMavLCYrdK~+lDVYmv1KTcN>)N*b`#OhJ)^W_O$zm6#)H9;Hi{H--QLGCS z8OWsTxIS13Cy9uZFupy4=dxS7_%mqZ!bx@OML30{>n-`U^I3ac47|DOJF=SMw~dKG zUrg5K(zh<=SYz({344Mnlr4ZkJxT8^!Ly#u!b|v$a@+WEy!mV2^U(_5qDy94)b0I1 z_jh!Gb7#7c<$5gGPD+MFCJ6!cC`MH-cMN{Ge_iwQo5|9mFoCGMC)kD_t$;N_%ZuYi8kYtEp*Kw-hx&Q(vkMOX#@$?MCI_ zwrpAP=T8Bv%KX`!ZyD)&dpR`^spO^0WDFTrEhavaq&Ry#G&;oHW#8`CI)b^`DGYRR8d6dcq`R9)UzuHmfn?N2A+6*H3~3buF*xz})u$6LvAZP&q52**s# zoE`VB7Q^54Aoj^WIVpE?y&8VgQP5?LtN;;FDYWE?Qx)+%DRJ~;47{b6@+{KNYmshRd5&8RNabtc=# zsJwD78VK3X1wO-!#~kmuWn^`T8;f4AfimzU?VN3r^UA?Dk`%pUd(DNGd~R-=q`E@1 zs+D6rC4ib!88tfCjmlQKJv6I+JJw77_%;loy2)g@YV7fzKCo2Up&p-uLO;}5*&xDXH zkk)=3NOsqv;ogURHdw5spdibej+5ggZ7pF_P0Ofoj@>X2z{m^YB7c*k1(9qsxjry^ zv|X))TCCB9 zdza=qB_C;%k}Qm+cRF|!+matPKyM?OSu7iSU}Gm~xlSAX{f5>u{+{w|#}!&#BU?5P zS8YE?``QDEpc%?(O3xa2G&LPGOon(0kLf4jT{?cMn-YVlM`s+_LdeZc1pl$YhV?Yw z~R9D1qoE7!=|s!g|1Xl2OxFJm;+1aKA4Uzt_Tgt9I|pv%^vbVS@rmS%lh3NOY+TtD%LRQnD;! zGB5b-$e|qhU^cha`p1}!rgLv`lT6}l9w}{eYXiiEk>#lvU?Zo2tmSjsfO(o zXfWd~Rg`B3Ivu=g?o7v8I|a-&h7HGgn<^?eoe67~lWXk7rmGhG$aZ5tz)4)16)b(C zH|w6#-kR!Cg!Hv1c&HYGN-$Eta*N~2(t9~+lok32H@Du}QDJh0A43~noTm29eJm^X z^heNR`9URAI;#&=2x8fn`Kq&5v_FSKF$9Z{6Oj@Qd`Jn_hl^%p&&+|y^St8sGwpKn z$_k9{wdKO|%LT<%!Ia1Z&TMwS-O?KeJAJkw|WMN@^6sb=o*1zU)kQQl~~ifEAAH z&B-TsqwTa$Gh6q*tPSe1Hmz{6A{q6&E<}iNh+Az!5qNaHjba9uuEqxdyYYE2I!x-Q zMvdLlb9jV&bkVjGXJRB?zM|6tr+qUz1{5XU8f4K~(A*(hQL{jdP@r~oy$5ZMc)1=p z$(_qOe|{TLr&eP=lw{v^(&o5Rl>jQgG~V{|o4_2DmKhOSO^Duj*M(AXP+N7DSgYsC zi!{oEqQA=VJSo3m09%uhz2vmA-WpXoiu3K6ojJyIY^X1rcpnwZ3lS)~1>|*3qftSZxYp z8$_W{&CgVjV7$XsooIt~uGljzosn2pkOQcj$|W0UecN^FPCIEkn?|Brs=ZxboU0zV z>~x5!gBbQ5ge0LL5DcgL;=J~P6nEPT$n%IMJh}U~``4rPohANp`=J4L56{T!86J_u z8ltw?4(YYk-qX5R^B|9Fn&CqqCJ>hBO@{(3CO)))+i(~ZxojdX`qfqMaVm<%X+p$5 zbTOEA#6Sg@IH6cC=vv!0@!BAFygqBi2gh#=YD=S3ub^q|kMuP^M;CQi zQe#6hljBoE6s0yhCLYM`&4t>YW-@9urJutOWjf;4*m%n5G$_$*D@NOMea+s1u+(zR zIf#VR5{0{3wwKtg;<)xWaRoPf?r~b=m!;86WI&lQ(|-4p^RB(IDyvJf9USjsQ#RHU zpBQx%58p_b8JNa_1|U4Lopi3LPxyN3wuL#Vw0*peR-X|DYSJhS<_ikN>A;`ld4Cc` z6T}4@j*@7I(!30%Tjspna!GsNwu&m!tF8?UaU~I)iKGecNeJ<$w!sZ4#c}JJSF%F1%V!1!p*{Jq%zeZze{MqGEK^{GhqlQ)M zT|$BSJN4I$lD*Z=n#1ikgjd!D=?{!FwA{csWPT$*G=EPoX|7#afi9Nx!NrW+4Sa;c~JXd`r=d%A53){9F3})f9pw-3l^r} z*;**3#>ZH*U%aJib`ERoX)9jX-wY3G^>s>{zk> zFyHGT%M4tUiex3|n!yOibbNTiGnxMJo8OHeiN_Dvr9#u zHR`TcUrzho7T)IzOdmXFo;KIH6#C%YoJ%37eKzaVHtgmpRc6NhicWKf5(Ql@8jj8a z7k=;KtQRmgIsoDX@m^o3o~m-|Bq}#f%0m{v7pHL41S$U z=eefSqIOT&Tvyw`D5A*sS|=60W2^O1ow71K8$TDk_6eb3F}s*A7!AYmKsn;UM#bE0 zzp95lWIueURvplcwyZEBl84Z|`lF*_#YLbJ@dbtI(P}MGi8+a>PKQP9LqkNiI9ID$ z*)sA*adMF&2$0_E;j>>To0q-kI_BtI>eL0VPqVv@55_!g4*Bob6a-TxROTlw!B-fd ztXZG(!1wMg#=82_`kTV8jHWMFcfO+Og6Q)>Nuuh5$PZzJgkop^rdE3xq5~Ab(l0R0 zl0vHX9tLqQ#A!X-$lFHH#KUDF(_Rq_X0)9ay~>rYRx5k3Vqo?94^}@0+N*;ik0gR~ z${pz#(3kN2ZY^|6-g9X`LF$P_vLz6=1{PeU&sG2*ksTW zj+r<+o_W9&f~x&-XdHs9!A(9VAD_fYsWq*oU)V?jJ?lNFr_cajw}mdOpZrz4_FG;} z_uWv#vfLL7tX1V{Qj#3tJg0_jxX$6{>(qjNcFG~;VEN}Z?^b;Sql&{tDBtrv z)ZbvAm!2@Jl4kgz>EZw9pF!k(=hpVuVjz#_Si?dh^%7;NJc%V^lgX|-y%S$w=ibM2 ze)#^FeUZpqC$SJlo%F5KckHW*>s<59W!6?fAtaWKF!P*9rOdtHOAlE7`q8zHxxR;Y z!|!eCh;wyyK|#TEBfsnRwH)i;Le?2;b8xswQ1pzlh=`)o;6u+GR0!q+{54nIB7W|) z_UDg@JBMJ(rSg8?=^S1t`+%y2_bC~ftXwg6Ra(tT8u>2efaiB6Etuts6Th)19RHjBGb!c#!zxHBF+;lXw1- zl|Mo13hSX(PV3e8lr6TVcT>XPJunfGdrPiXkyTUYIy`wTsV$vr5zO_k;P?5vwOcD; zgSUE)aH}5AJ2ElZ2r`2$gA5Mn$E29W%p8J=-1B%_CP4_eY?z~ShiiqLZsxREeF4n> zVSoP%;ZExQt+ZX4#EEVU=gpKu3eW#(pt@%Ke9mWfcGfR&qd~aC^96w9X60}AQsd0R z)R$fuY}aP4c#pU>jV1;G>*z$6{yH%iA zy}0T2+L~^{KG^D87q+wGoC3$y8TCuzO*<=J62`^YdF&X1u)!ArKTFlb(|NX$KReQ< zLSBbn*|B16RDQcf0k#`L26JD?3jc4I{TcdFld$k`*Vm~uf0qX_A7O*4fwM;G_)XW~ zEt1Skhk~M_owOzHT+gXbXt+&dxTsLT&yVE@W!98BYxAk@z@M=B3)R*v#(=BO(ifi2 z(e-z*!Bj~kepMr}WRL7Z2ZFX8l!F3_;3SVy*?@QMcJC%qlvdBZQ(!E3A`2V5q z{~%Pr63|wn)qloJE8MvZM%nmExkh`u-SoBHaBv6(05>6Sn=*f{$fA)SoLJO;izL3Y z47>mJj@CuhJQb1a0HVssugz$fo0YKB*3>lBUkcS7IH-CJ62irM#ar);JJOGP(K>XM zC^t-wnP9CB0X$DrGj?)~UtP3ELBpUs_873f;uvU+z}`S_ht5H~ke zw{NlofX2fj(te7|A`J}Ml&0!7eCCNW^pN5czm5#9!N(dj$-`S7biTDr%!kuGv^|y@Mi$_+xa+3+*5x{bHi?jw z>rKT~+-I}En9WR66R1TK4D>+hDhL|5E^z&ztf1$$1Ki1XD-w~rz=YsPkoN{HcBPKH z-dC6~!VL}WCz^$Qt`Szi>xN3Mm*tHmF=`sMa$*o90O>#+;@`kt4dy7wMLsCNp6O_r zNmkT-F{pyd26JD~=Mk1WbYjA}5QtsY+6x8HaYkWe^f>r`DGTV*VyDVA2NjyO%US#T zeNE7~9gwT$XVB4H<1YZ6tt#kyz*KyG`%){Y>sC*yLXx_AuIj;kYN+p2N4~YX80c@+ zau&(jqPY(>Cyqp7LJP_s6qyTwvzAOqSm!%G_6U6b5}t3;oJ~A!`EF#kb$@mA;qT9fLpQFYkUaI#id;d9PavKKlc&YkX|0x zopI}S&(v={y9_+1*SmcNWs-W`3)54_ch z>mMa8v>>D-Ff!&!ffBo5F$SZDvJ)_5Hk!aeGZX-ZM*o^K{_@#h9!09KYvcwB2m{yt z{aZ5DNvYQDt2d}6-(b$P3tyPYfcj5hu?;9Qa`5H<&eMczkYt~`x~m8b@M8dG=Y)pULjDp1 zq4*)6frFoo!gPlg>}#-$>?;oRysfDTLD?={%`F&AJg-D4VOLPEZb7I!K=^>Mjh=V? z#IXyUE2)UL&%qbF;^!bsI3aR%$HfJ9!ndD!v>%$lf zfU#9Bj}_q4Ub77Nq>ig~09>OtfW(^Ns+Irf!>O;Uc_LdyktCM2L&c5M->Zrm|N7fIIKWiDuIsAh^7a~^{idY0qu{o!a2*+ zx_e{`27Y>E4n*7M3KUvy-|DTEm+Z;jQ?KZd7eKWxqpPj4C6?@!z4rFgVj=!)&rz8F zhxs3SMXT-`x4bVxDno;W)h83s5Jx(PVTE@1QMw>cJd(R*O% zqRen5)4XF$t)Yw0xZx*(4$&*=xJ%UwrQW|PjeEtFGl4f((P zk?awiG&8VgaA7ylBi>Q!)&}xQu-BB(JKvUM-gedsj(ULK!~v@GnV96Dw3UeRU5KCv0EzGbzEZ8THK0ynIih!RjI#XGBwZo!`oKn@&gA%{C@%l8eanO zs9s-RwEskOw0BE>bYE~?tjwIX&xxdhLhnI!z=oM}pF-6|uzV<4u{Kj&g2Gkt)A6HtMNwdouR2LRm!2Ky*ra12Un-VHODTGisxRTv8(f#@7tGeBjkJk z^Rx7q_}sMcBtq~|v5W_&@|k7zt?G6jqEiDFWUZr!94J>zXZKE;wTk4n$TBBi|<^MNQ)SWF^mpw;}Nm1 z8|%7dD5TIi%C|=K5D~fAq0nclNlq@$Sf>)^t$rtUl^0HF0YkX=S*-ylfC&7)f9Z3} z6z=!acDb;5pWlwko-KU)_kr4_;%0CiQrLF&t(0cf#f`k3k97Reko#O}XBYuXW%}5! z2x@P$P{54@o5@sn#{`?oT3mKVXw~Q}gQ&Jutxr)Wu!0PKx|%gMrc;Ex3{7Ey~Eg(7me*L(vhQLR@|FSRwIkIjhwZ5+p14FTGtkm z#@G0QT%b=iOG(X?b_^#h&OU5Y6W;D@R{i-|GQH3~Q1w8gRr<2(k+|@%7B=<6Ltjpb zPXUxpR%uI+*ChBu&2!d%STT8PGJIon$yy&+`ib#zOOIej^nv+M_XFCR*gaEM^v5ug zO<=Z~0a8>ZcvShZPmB<+Xai1;2!bUHMg#uzt~|@F2VlM^#wP*jH_PjSJb)6S&x_6O zV62rNj;b!?&U+0!sFaO|AJEmBf19?Q96YBU+eMCL|RU+w|Zz{n(dsTHDLRbcAE;`{=EMFTn%FlcT9Q z((S19*J?%0o+-6g@b4pRiowe?3xk#p=@+2D4g;deRZSmoy$vuA01G2B0fp&Bt@AVE|)?N9g2uqKdX3C1?Vg!xd8(wmSGUTuMX))VN0~gVXVSd%kt( zzrEbJT0OsvzI(Z$kUev9(b{D{hkF!44YIb_dXcHOPefSr9vSO!@|NWo@V~Gm0!tcr zw+dJA_>aWuf_IyQEGXG;rm?{g#vAFsO^ z)N7S>Og*QhrRGU*6?dPB^g#WIgF-3pdNdvBIny`4I}sP3iWgG_-R+J{@@P|^q|fxP z2|g2dZaRmV_UzaGmcfYVdZgIZ>@Is}+s8MA#l%ES^-CuD$n zj_j8^`UD#WO0*0>@*+^3jy^iTyBim`cE)c;_SJ1%Hl_p^&A#WEyt2kM`r@%`fW)cZ zX$6-6Ad9OBkFM`w1y=C=4bW0_%F~Wd9Z$c=+Vk(C>ssg1#0^w|M zV}*3DCo{@vCc12Cos>av3|ea)lphzKOPRBN!_UFaqxliek~?Mm6~vNe)#YXowI3pO zdLF5oxqQv5xO;t})msb1CeUcmUcK*{YxWS8sB>)zBx(mLI8Kgq7f(R~98Sp93fwG* z!X*d}v1RjiD;^z?ZAa*_3q9}-9z9`?V>fb4Bc7UVh12hT9dh-#&y%1DoTq{pUf7x_ z55>n8HYp4Us%L7?cqQI7S_sYifbw`)GeDKjO+hgOVWcLYgT})W={kwg-Jm$eY^>cQ z8V$PbR(#QMeZ?q`m)h-LAq+gR#Md6+=_1f5di}OW-{WR+Cl-$^Km)YWLus08{AO#n z1fl-{XjH4cEV(79}GSK?<5y*b!Zb z?Au;WaKJkI%oXP`7oK!&m9uW4%YC)qjuYm^Zz4pQu4f#BSH4f~h%>ibmExgvXL944 zRF(rFMU3`gJOgA*hlUgohKCOgdU~U;8+mzMf78wIS#T8!e;iDZrTPHxU^l*l5Wlbz z8r!K?$)B7XYooXyv*Ea3hzGq(Z7h1y3QTs ztg((s0-qaHA=%DzM8{F!Ko8lirz<7til(vJ!JG(Wq_B3C*f zhU0;#nsV||5)}^6Lqt{E#Tn3luJS6MDs$gi+KK{@q|?l!?zo&#j^u>zMI+lqfU+_z z=Xrf%AsHX=0dPERdJ49BI?gx{?mqw!?80xU&>K`Y3Sz$lD2R5Sff|+C<*VkTh~D&W zakzym>RAY#OLG!VnGhOUtJv&7H(v9WLaM(GH8W~>35t@ajAYjrqk{6Ui~VQ6Y9f~9 zrRF|}U1zd;WV-3ihPU7-g5Tu)IMl1ivhl!dvyqM6KRX`h`6_xTUN;SLd(PoxS!k&u zxA3M<4m(uVo1+Et;qa-;&i;uLtF(lXFfi?io^Qc z%%ir#0?P3^%k>1#&g9fj`bF|R`{s)jr7N$c_Puxe8k=X1HfnKC(<{u?rq-;*NApcye-n?2P9D$Uu^)S=>f zd$W@_?mkB)?(eoR+oK6=NdBB*vZ`Cy#!*tR7!2EGC+R-}c7_d0JKxVNXyYeC{cF|Xo zq)^qbINdDCdsJ&`r01?d93F0erH@9aGU6>=*-ldFc~2&$e&RHPxUj(Mg74sk z0r}4J9%%E@wKW|;=H3U8*5UbB`QcC=a6#f9wNaosf|v)fHqPEdGtD`DGYLnWZ>v<3 z-rx8&V)9CzQAt5fTK2{*vYxVYFbMiMVh%o43$8PsJ6{jKCXR(F_nK$wmg* z@+r=epfFn2S5g2EFOR7BbV}NTYx-s%4Va&v*=5&rzVg2_>ngz4e{SA=Q@B2`yO69oy)uW7PZr|YsPW-ppjhiNw&Vr{TE=3#xnY%0f4ukfqLgNRhcuEB zi)kn99GP5yQNV~ul!{ip%A;y%dgZgnT)MEAE5IgyLJeeBA#1}-4tU8)&`q}-O2bUM zD{7MoQMfsqm;u*-=O=VEHJJ=0K}9U!c}9V0(ap=TK-GgVn#?<%du{rtf}md31&Np!30aLZzy^^t&a zOL3+Zj;-AI?NlX)6#nqVidFt=8RZvGk$q<3w|_;Fo=(rNpfnkQQ``76HZsIA@64}7 zBS+RK7JCSh?3}#m>ppi~(lr0A9wR z(i%PZn4BnylpLCIY%7wV;wMboA9`Os;{=S)Wys+SZDX9|yho5bnl)3Nx%gEK z;f&+&T622sIYl5m`RkkNj9rA2^ZD+c2ulP%cDJbh_@v9f*G$xP8X(2U6p%av)i%AVj%qk0IrkX`F zm-3|^rz_Gc)YJL$JQd{Fu1X=_c6YVC${MJME2sSA)QIua@xWX7^@2tR1kx866eu6$ zL_rB57%D+D*mmqCDz&<(mcCR5-OYPzTJMmUK~aEIF1J9V#m< zn|L1nfN<0eoVEzPIg%f{&PNF~<#_eH+u_0B*5UNQ;F$;s8pDNSxN>BXOeO4)O%%X44k zevHCNw5u2L`Q^;-AS9eOp%Jn0QNYfs8Oex_OO$f7f96%MOPAoYfrbW$0?AYIz(*k-yKp#iM$5%1om?;QI; zf5wv&NE_wBxw=@=d#wp_CtuwlrHM;2wUuKmWebqg-Kn(_sr&edp;a_utV8ZJ<%TROY<#M)TYqeD$>hTcvEWGc2cC6apYm_T{fm59u$#z!hSt^SaOtM!r z&SrV&0^)s_|4);J$f~CjChTHg0O?KJNMqJ|_CEEPHwra{j5wu!UlW9XL)H&=aiTV4 z@mg@oAq8)iwfwn}nq-TPXl+d}JD^AvU4`?XOa+#IW_yPVG%;LP8Z2G4)`OU$ps5*t zc8;SpUJH0Sj!@TWn~)vQ0d?Jowcfn+t;@ZY2RxS}5o*i&dK`hE;UN$M-Jv$;BfR!# zvYpDIf5sl52htuUdN-8WzQEo81CL$LCJZaMECl{IYDk}Kt7vsuj&i0q9Q>ll9f8sWzpX0s^f2s zmrN@M${;pI84wPNQ<*_Q{uwk;db213x!!A-da?6C(U~YccHwX8f&*f4%k#=kIc8>KU@@;59*>Bdydz2eQphP*zg=;x>Lwe9E zk$4lRiRI2t88}s)JQMu=)Ph|y7{3!+Qy+xh;*F&1gBU8-+*&;9P#(O4G*ch|Y1 zg{JJS5@-cY7$j2pn)1Js93{G?${h;GfAO;zWEb^Z+qRD0+5ofRkkjP0H^8$~rC1JK zAV4bt`^okTlG4^t)9h*LQn};v0qoDy;^zDdyXK^hf2@}`K+1beS2b! zKf&GAs38=pnXA%zWG8b>=+J@;=&O~%ZpU%+>U_O*l$pLW{PjuDhS)YXvCi{_HH8I!7h2E6;5%YH*|jdWxL5T+v#CJ>F}ERPfh1pbnSJMT{2m|kp)1){`2 zFYv`nMr<3TzpkzOTF4G;6F77xV=@_}5m6`|c4>x8W33@xFK$tlM}MUjMDV zJQ(6p_ws0tiMhBA2xn)Gih@G=-uP?ELjTPqXiAuf>`n|-U`LFMewd{vWEW!*Hm1y- zGEF)xaUgX_Tv9nHd-=r^wKrudfBJQUv<=Y1_C`o7fYS8~Dl&!84lt_GE%0#|dHTaf zcr92{Fwf$GI70{_XUh*QU)^Q1F#slUJ4+5rAAKdzI9wgD81sTgiX8V)( z4vnJXHgFR&t%gs}9v>O}$S2Wa(|}bw2F(!s(l`3_vCg;3&llKR^%S4ac54XB3;6zc z*MA2Jvnz)aXkkf8Rh}P^?qCFKtoHEIgs{Kigs_L(@Io|cxZ>x>!HoeLx_X~nsR!>j zdx}e^zJr_2km}Nk3`JjQ{l!P&a4t5Y`)3F6o~F|%9~}4}v^(=4&*FRz2xSRi z%?E&Z$>szA_dyU%W#OXw2|z0Bf1BNP%Tg$A<+SFa}k0-lHDcT4UTs2(n0ZF=%`N;`hmA0FnLd;m!WFu9$ z9cbG!(j*X$MS@C-uWq-eSVFxx=x+<5IafyQy^a^$$;thuQXhOw z$ppDvh#wpE=cDUvP}%^si=>>IMW?9)VCSbm@h&0afhFXy0N;vg9A4&ZDn=?P07SfZ zrYr?ODtRe&zDfpQEa}Bp1bQNjUg09my$zWP059clpmld1KAQVJqoyGh4J@dn@u$=!Rn?Lj;sHNdTtKF7k(>ws zT}6=|1=>6`VuoUQ`Fu|E99+u%0Bb2uq5mA{(aH%5BpY(ilVQO^5SDsNYUH2^^AHi9 z?wWTNHIRps0Kk|wm9ZamaRSibIsZL1z`G%B_iO3;)}yrVFAR-Puk z=jxPZ2Us@nYP?rDSDqiUMFOI>9^uXgYi;PdK#!qJfMFwE{ZS6EA^jxp&7m^m>^K=< zr`trT2B0tV?$gNK7J-mD*04i>)RFRTmD*1=uv=xwq1z$A0U1(6;a>>_8FIBE;7F*4VNkLFnLeWD05NzD{W_^)C@rqw@enTE2?io-B|rlD zQ9gh7oYIK;7wAj--cudX>QaEg31t6CspkK|UHIqXuzGc7fzqx1h-i=J;z!6A+#qS<*`K z_DCQ81%dPa?+`fH$10hg#{)mpO*lY$_6r?%2MiOiX$@vASF#Zhe+62hf<~8szXA#R z-kb@HmHdbNglYjq6mvjCx(4Z{jCyOd*1p5hG4IW5OBFvfG$bQ3 zqJkCf-Otf^CbWC&hC79ooJJxdVW9>Dmz7STW>1Gwe-5Ll6ViopO3soh~FwQQ-FAKIXS)~tc9p&hDXzqJ^IEK*`s{ud>s z527f=r_&+Z1S@4*z1P47{I75_;(H&O148@u1CR~=hnS&Z#yyXRH%Hx_plw!YI|qbk zXeq*fO_=34T01TTs2N@~HD#%P8$xMcnLEAU!KN=5ra^mb3{8k};Iy*&Bdl0}5AFDi zpV2BV@c!>zK2vCw04niW3}iN~?;&P{s7RT682tYK;y%qmw;6!6Cn_>3A?Z6y?2gH` zh~EN4i@ScGCyQSv^AE)To8a^htATMaRknIF*Cr`g6XLQ+wgwg)a~UOzfz@#pXcxg2 zZ}kuP3U?jiT1cq;=+FVy!BT-p5YX)uOdwK>MF(x|krK}YL;?BPjF4`JJj)%r$rTJn zD6wKT)$wt+pJ@4no??ggkw3@kChS-3fNpjJ3;8d;X;th$1dO%5qBnkw##66zD^$~r zJXa;YG;PmV1LE*Ms#dktO|rq>$DbJ#`mdzu{V&4ut^)t7ukL^D<@MiTGygx5>Mnz| zbn2Gn&~a~7;)xzxzz$TLdw-ESUM|$RN&uW10^?C7xGA0NBDpMXsmWl+Ag2{JNLAgu z>>N&bb1b4Wq+8?#NXRf!eB4G1| zrp`q#;`;iyxn-ONOPkh5{b^o;NP=B^l5PS&^FaAHu5-N;YzQ>&_&4Hj;c;s2?(WCP zUIR!TFZ9Af6(#Mx^2_V%UkA$pkjrcYwJh$eh{zute5(^XTeLcSCMrF6XVW?O4jLYHB;YGY=Sl;xH- zGISH_T^1vA6)S!}8GId=JX`nn*J+xFSq2KDsTQBsWz2@lBuyi`$97*g#?AQc8GNC-+1Pldx}RlMAM<=C9+ajb3q_-xf?ueNZ^Tm5k z)$8AN9%|7|06bEAQjN1LtOlE91_TYq<&?>o78QcE!4n#xD(^ehwuuY;^ai=GsKE7{-_C5Z6D?|D0eSb&S)z->D5_lIho&!nX{*CKzWyd;$ z+rd46<|m|ff@$~hzc+udmE-wbmI%npf52W{vINGtwkCt3w)}Gsp&rzJ<9Iv#Uof;) zZYp{>-2kyrp>B`oY_2c}_;I+(4*y3;N6l)t4lsI1UH`|=9s27yh}4HZip5%V&xr-s z&b1OsoksTtDZ)P&UUW5)PRKDragjJTwn2pT>dZ9+>Ze9OJ zhU^Ay&K9 zH<0y%X)26kz0y-v+R5;a=k)hy&hCvAlwbZ*%;dDPFvZdNY3G#AiQKbiO6?h3yT5z$ zC#1QD(QsvP6xnYyHmIG$>2TuuYG5uae_mSH7y7)2*VJ2rs0gRv!=L;-sq#xaxk7$F z$~@^5IGrxGJ<(aWaedC_+2+2!D66H`R%QQ1o*=iiSS{K)!Qo}!|BJo1ev7h;{zgYp zQS=c}krGf4F=!AZ#~@T%N+gvMk&@0qP(eWuK|*2>loSx@9F^{pW&r6Px`vrK>*jgB z=e^E(uj^dr4|wMXW9E(>D?aP9*53O!5iLr}V6WJ6iR*Hs9L2<{>FR0yDQKbEU|!m& z=2kA7OX5`3YwIF!ifI7qw9=E!$>uf=H*vOUc_ca}q4yo26YSBN|2&^OClUH{)E` zC7)C&qQxmCF%lj&4hhe1&k-vG$wrhpYGLmUc7u1AhVQ51cd}KDh^vm=DKcKeGcM8& z6gm7hA8oJI&jIxHxw+x_`siz@c_yR+xGsTv*wxw2;f}oS=R%KO(Tr2i0N!p#h7@Yp;=+4%hE^r(xHv9fFwuWj%ezj;iZsy+`b)ef-AKg=lqYGvxi*X%6w8 zh<6mrKccb7COl82ig6sb#9V3M{}vSBib&$#wW9pjN(49(e-xyrx^*j$nmt6#7b7N^ zX~R#|fY=;JAy`q%%6gjY;>bM;ArfA{`0w)tOu@Sx(7kcn%68$@h(BMN2GP<`nmkc?AptfD!nV) z;g7hha{T10biAd-xiq6iJ>2_8QiMSV(3EUw1pA(*YY$fv$fQDwm)lI@)EFnH%-Wrj zwv^%)UabFu%hD)?yoMX@p7a`e7jSUT9$=)_DreLCXMK+lD7|Z$Gj3Err+U#5|X%X+3J*qoecltLFQbg z?1nzM7B~U#4nLWkcd3CtC{Q&%O2XS=DQsk}h{!z!2F`I`Ee5kQY%6KrwuEyMk^DO0 zVYz*GB6-#~n9rRu5=tPCm(5i?zp?CX>(cGRfRaF>{@x>=YpI4}a#UZ&%g!ijETG3E?1>opE0{e~{XgTpCbIIM};us%4hL)1}xqHW_*LxjVEU&L`g ztRtU<<)}fN7cE+vQi`6P)^{-a$V(t?>~vAo2%DkCH=RoG zURSI#2AX|~DX-8M(!Un}wJLlzw^5B4+^k8NE15Xl(W|hPDWY8(f$^mtd+(oU0}C{% zP~q9CEprHkL^ySM1j-CIws^FnzECQEx_H3E;ly`ksJ`wuyu=_=JYB4zCo{ORm8Y7l zKb%as7l;VWIp_0XDquXbV>e!5t0upqT~dn2e>e1}H?l`lrb2>n>#G#Wg&Cqst79{8 zEl62=T>^%(fH5N9=`(t0`zAI~0C#R@lA&gRo^XMUh=_U}rx~VC$Y!;>u~6<=AKs$z zIW{Cn%DP6URBRrPFXeG)ktCmLm8+Vnprl|uC+p=}wF-x9ia#QqeB3X*$>2o^#j+U} z(4g)j`y306@7WMPCbZ%a4Vu-XqoZLLYtkxJ0>F>Fe@;TixU8&(90duwC}hRcf-Lw%aCnKWI)cL&wR1T@8s*JXe&db_b9tRjh(=O#cdPc^GpE{BGu zNT1C7?I)(q{|@8}`&@bTWOQD~`eQd!vQdr2PU&17E!%b{& z$K5#DirTmbQQo@ulM~!}0AaqCXZ>jNu#Q7{PjY1^_x$Z9Zpq(S_-4H4Sl_QPThcBCRB`4r;tVOhp`ROI%&(vSqwKL$(K(#unzFtq=YkcQf648l8Hbu@1L%pGPEJ!DcG&^( z*^q0mE~qJ#vG2*X4q+93vUQDaIjOrA5hom}gc(BxA^?J7v8QT}GcR`xUme;`LDbb4 z)Fb7O>z+X_m_LH*Ko%~Yi#}TD70`S!&U5G#^%u*X9WbHG;2**TzJhLo#+;OSlpDh4 zwB!4G*LPObpCAMNTGh!V$ku!9>ZS-k?Y-!<`?`8g+u*NlJWbX+E;kP=&k?2;ZG$L1 zl-V>wM#j!3i_tuM_xA7Jvsr^%@O&dHJ(EY?zgjU%X;!#m&Zh^AHxco=Gr4HZKFZEk zG+$u#DU$3cvWTxe#_Zj!+EU6#Y1>C_6lOZ#t-|X3+h~A|IP1?(HuxY|vE80M=+J)m z6yZt}|CSrWa3y|EAM{#dR=Px#tqWtY@f19}7vj}*J{fp+bYFjl^UEgIyY+Z|P2x@3 zmmsW*>=A0LvoiW)n`qPQCAa3?+kon`qz(`Vv8Zlv(rs+acWNM5!K5lS$0 ziHpH~{k$Qn#Bq?LExX1klSBXDCR3&`n+o`1u{h!+nfww^SUIMB_N#y18mD-IE3VN0 zpX*!A8vb!hrY^qW1%y*MZf@<}-GACQ`^p9{KYY699At!6@|dL8c~Bclz!dpN|J>3; zQ*a;iwIanZWki-S$~w+@e4Y-~i{JBUEZZyDlvOJYYz1MhyEYlYl$jR3P%*A~R0Cx^ zfaEEhPpMWc&yOH7H@wC`zW&*)bOZr%M*w_c=HlKg#=8#I5+GxZ`ks-C?At144-jy8 z$+=q6$ch~GgGIGb6~w`tJWY=)aL(6vl8Q!H)$jkw$*NzRP>jRIR}!JuoJ_+y|Sr_7_{iDQz&O@q?`CA>lDyJ2JEP|+&#M>PHh4E;u_bU z0FS=*eTe#)igeo_Q@#7B-ezPH+#d|nam%X^e)D2um{}4hLygQ!cwkTM8Dt&i|F&(v zJ3Pd(GKfhy_|ny@EJN+vMx&hGbn=SC zd~e-CsS6Y^a2Ge6f!~ zZNO*0q@o1Xa9&+DdaW+YRxM52pdZzZkQV&u%mv2F=XenNI+;<>RIMRrrORrxq87Wj zR-xF+%u;>upq@-?JmsO3Frn*sHoKssostx3oKetGZC>la3C`EYgW*Asy|kQ~;x``- zQnosCJl8+*%}Em~xfSYt$d_v$E=td&az6g;uNoYB`%pDBWh~k)y9QokvQefpMiwOq zKb{_mdS1V+C&T~Xht#6S?o8B1GkB~bDV zruTm@MI`kG^~RO?Q#zQ*Z{a4JXtHDEi-vGTzB}8V6Co+dM+%|b~r zcol5P3k@2bZH=WD8(Os+bR|^A$+(O%nLJ*rtJDZ?hR*7{cSD*o%tRv?>-dL`E=T-j z@s6-X^|+)rtMFu5@9xk+y%q)8ptBN1cvAJNgZ8U9<+tlvEGGS8xWTN5uwV%eK546G zgRZJm*Qn;xkvQ1r_~7PNr@DH+;KYLgrmNQZ$2pK8&3CSqKgN1C*M_75=LW=OpdGXa zcWuK?;wtlOy^9SPJGvA0k;ZD{f(-N@wEZo1d~`4H4ZDoEeL1R>2FdC;CYp822h!c9 zPT{e(wAxKI(byZqh39eEUX$PQZpz7Z)kH~C@-G>c7vPImWjvqv|4bH+dv!h&zg!v} z7jN|4)J=bvgst)AYV7PrMyzIjJMC9HytS0%=>z3BuxM4nzLsE~FX=p^%i6m2gvHhI zkK{y#n5yy-xM^rpqa=ZNQOSpO(hiFS@#Fm+G>Q zsVGtATX;a2T8bV&Up<*To0^&)Ok~4Iv~1VVLhO!<=`&mHjp&a^i5Ax=a9#cOaU!F{ zCidum0>#kJ=g%aPK#54s$$ZBi$Lt)QZ8$Kq@yjeQvYKuRQ`~v6aK=3Ng6@U}%x4QZ zggyBmrG_*7fesbRl|_`Clu~I!3ErFcT%0-f_ZlE(E8*zo*CS{7lsl*Mt=sH}$G#c_ z(<+Bf>!*SXQYWLB|b zw2T&eHVLFl{y0;(^JHT|RGhj0c~k*XX@M&FTTvUYazyFv(W8$+)@}O^PIFrCk08bn zKRQk#H)`;e8<#G=>{Qp#M-@QTC5-JqAk9T)Incp>+1^WhJOSl2NAEjtUIxMHj6NS3e@*nF_tv;`Ed!4( zTWgcYweG^`zUS))EzEYh_C47Wn)G$wiT|dc;A_rm2qHngYFh2Y$Bfb@oVe)aw2I8*uhRf@4IsVVLN9aj=NB zZk9C%jHgVI?f)-`n6w0{48DJ~cs1DV~wehqW0Q{$3iF zT^n=6Nb_PXn%r^*nxBAe9GjAa>}LUdy#>qE+b%-b{3K+a&558SeP#nVv`7hz0M>_-@=d z6W%eLfPivw}0UwxR0O5HA~ISknOSDeGn63-1;5XL;YNoP1s-h z=AfC`K^D6X7A=up}x2e|RsS2YMFB?5!cCWSX1fxxSJyx3H`$0RuDz};9F&Bc_9j8`h(;zq7 zQl`R%Yu|ZDcMc2w>~3gisdqf`lP!H6gK=M<&A;5ajttgszDJ9L#UxISr%#`{<91gX zYs*$!BuRQL<_{)c%_9;iysO(%uKyL4&k;6Vn(ckIlt6NhI6GzSHJUmN0LM>KF?_doO-}RO!CK@DK2i zFAiKW9n7=rRKdWeV8fJ1-0p-QY!KJtL+dqt`(S>4KBpcv7yNa3GmkVFrk#~Xyrr|V zu|V_>yl5b`UZtKiDJv^=o%G@j&TZ*e_rCfxT$HnpoUBy%A{o286>(M1fH>Ya6}vj} zj~L$7YR-x+N`|y$q!cG=`SIHhAX@% zcCB7GnQZ0UUOoOO*P$1!(!}|QzN_@~^wmcj8ff}OKgJ1rZ0tj1K3u*Kc7h2AaNt&o zVxe=uNt!iGVFVYluA9vF>h7dIi4_OPbG? z*(`3wTW`>;1-Y{Yu1R6MxP!c0j2SvfKE#{&-tV}L*YqDNoOAQLNTZ!I0>V}3glNM-VUlbCpe+%TVaP6*UFp{!viV}onT}Pd{%Bhd^I)`atN7+%akO`~v^p~SrHdt2O~~eWsz#p6)wnoO>wc%` ztRv0`Jy<&z{^>-y|4CQC(6{$8QnNa?{U9#R6}_N%Z7!6QPa9hR;uvT~2FC91D+Y49 zn@*UYumP=Phz9{j=;sWPdVVu+7(ev~I4)1hoDOE8=Js)kt?c2L=lCgM148^$e2G)- zAX2-)irHQ6RD_R^0brPUr2|)L)Xq2emadJ$FkrvSQBC_1G-r_L@_pH+oGhj)70?PLvC*F)v}P%WIH> zGOGiqbAAV|tOe-in*S4X-F_5vSi||Bqi1*+)Ah@I1Ju}TOmVjqSyk<`mw&`O+bL(u;*&?b@J9&N^6{?%dO()s0RMb_liH&gue}F|!MmW)FDHdX=CP~ps&*A@Vr6J(C}goQPu1VDKl*xTAh{BR*@?Zi$Wf4jAT(j1NLvFgc&A7KGj6-_rN6ubCklB&#q-_H z)I~nbax>mp8^f(tH{OqEpMQehUgp&)kiz`+J(%J>n#9gRUeCaZdd(LJTd=u~yqKs| z2uk&KZi%I0izWXorj=Q`F(<>S8Z5oqq|4umjZEvXOsc5x?q!7yGiS0fq&bFGG-X#D zU}-_&@>nen*(ud0&i}@%OlL$!b~g%fUN?K2EmG_jUmkIw6`vjnH#y3{!{}XilbW)Q zV4QE4`1y-Bta?r>9I6bT>=UfEK91cPp`wjN{s0H!23OXCV@LOQ_aWek#inGrTd7OM zeb0dD4exZjOtoSuGohRcu)7!=;CuBowMQA@5fKrCu4C_Xdb+#WfQ8)hze--n5L1}# zu1t}yN`p2O0n8d5YsO-g3gwo;-9FniyRm1+bZ`}$Wi17FY}S*i!@P(aq)O>c(oPxc zu2YAyZN>H^Q>2miHstt^!@ELZ{6AtSn;1&=zX~gHXT9rcwr6T_=v)}9{>nI$No|X} z8@Wt1Pg)&vgmwy>c{}Z_afaxpYtw~B@o2&+o!tQ)+rj` z@!mT#g8wRX5AcN!-5j^*`0nkId`#e^M>@KRGVP@hU%M2>(c*XS4GVYqvWeLH&4?5Mn~Y!I{b(~TFPZ&NB!0%D!RW~w_YAeQ-vD}aL1;DxVWzI<`q zUYS6)uZ8-PxDKfE!C*gr?8+;J*z!r6#@DZ3xdzC}H`5VGp?pA2JME+d+t#uxdhc{% zy`?&bg+Y-Ho2S0NV`*twe>o*~MRR^jK2MSQ7~}Vyw%fE`Q(>T0O-ivY$pz8!P$pP_Gy zYighM&hS1A_W&P`E8ukAD4tFXyr>!@x=^Bo6P~w0-U`6a*P{{q`7?rr*Gd?+Z2hl4 zJNaqyFq88|ZrJWOTgZUwh6N%Gtu~D$7mD?E6s~N2(Z3_EbTy1XG)x|*b^mnZl}$9b zv6--o%I=;XXp5=MDwQ_QWLrAG;C;#ES67TO;mw0*c ziqdg+4LCFzGDSggpmJ;{DHU|LRN2lYg%?TySBErs1%1`{hDND%SFknlAT5FkcsybQvAc`#DAQ+CxXXJmlzl(_ zld10{|5_at%BE7}z7uqGg{1Z5sKU`|_+RM-jNGlr$(FbQ9AtA%CDGbB|_tEV$M?eLRLLyL`?ilPt8GC@2 z2oK(2V8Q~3k@w#vI;K)wNM)Vt(U`n6QbgH7gmNq=)9d6IbMSJV?(iwsr8?wd? z3f&ugU{;fYHDg40XUzzZ(O*VlhhJieslxBQe2Uv0_9YtOprs_Fm`iu zTL%XQ>LUjMcR|cGtL`=AM8Lg~MfKWDwqnuG$-CQfupJ$Q(x2F_IHu*~Kye?cH1(@u zq<0yRATF|+<|Su(wfBgYw3K+}&h~8tP)`55k@0?diG1gUey>8W)5aEBIyzGn>T@sH zcSk6NE|f>=-P`2a;T^7<-s9gB&OAQWz*nS9bB|3ict4*7K6v3qX@r5*!GeLwko#7- zyaVt13&Qg6##jDsZ=cM`r$5&ex51?)wgGBrvi+eI7(Uz^xiZ3z0AP;~EmByUD_ zeE2dI%JJ06lh5x%8-VVoJq=b8zeOa9Wwan|nSOEWV1pNE;2~%W@@eyhFXFD-`zGl@ z1E8@quV})ptgY8I;`T8HJ$(3Z7lK<I8o4E*koxiFacv{>NF}Ro1V+`1J&!<*>NN z#_VfP2V~70aKxp#QTlU>xpK58hK=NH)H{lluhWc;ax)qxaCv?Nd?i$kD zYuB#H3_aMOGB$qu=1or8|J}%PK?NshO&AV`$(=n0m;10?0&=*-(#B>(2Qkm6HmX-N z$a7FW@<=C<53@Q{OP-@YaiVR0ZZ1Ot8G$zbf4^g9ZU~)*|KDI&PrXrWl~p_$G<)o) zv55U>mEZqN(<>VHB(>HU1xN&>vYitM=vyljue}IjYMhV=O$7vJ?T$8?{gX!23^~OZ zG&(l+@8Xp;x>q!V>8YuQ^9B~=%NNm((`s+tylG$O9!Hfhu(2(+t&+{m%)+!n=#U0M zvj?AAS}sd^?d09fCepG>yITis7T}=~@tc(P5brewPzDRZmav4~uP?J#N+>8J~FB&%S&Bl@vot z^s^q6FJcz*dmcve>P2L_$D&$1MqW^R6Odtg1Lx$fsozpjQTg%TjkiAgHl99xT8W2; zCpKwHxTB+^KvqsJKaK6G&nkp9%Z9SDFAcJZg#`sQQo#1KX0UO?`$jIGzA|kJVy`|k zGh-Vs6h<8@$_1ngj}=6h;j}N=kmqjwl;w)=g&V%`tOA|KKmYpvi1)M z!=miYw+^guy{h*L@#V_{{++*Dr62Ry}7>W~cUz^6+4p&IuTXzMC z_qhsMRNvSbbB3Ni->W@QE_=oxIfxzYT}wqPwItbPa;PmIK@)E6_7QIb7+IK^las%G zwXSb!QYu7@Q7xiIYp+6PQAs>nJ^5UTx$JY@bKPt<(E{!Nd&yu^l zyBYbb^3Y2^sc;sA(<8I_#=m?M-Bu?j%884m=$=*a6jiJ!;s9BoPxnA9pV`~nPY4PM ziefO7$std#J9h$vxwzz7`$XWaK6X1ep zm5&Fo2Cd%ZKw2H}uMADt%sqQt!MF179YDPP$=OE_A5N{Suh%VsR&bcHpN2;sx)!l} z?1VxH8>`2)$jchZO52?pa1uc5>{*oi8FqF;6LjyiS0q_Kv@=CDCOF_&KNf}J8@eld zQ}e9&b3eEYg<>qbqxL zU(u|YV}P*0{4%7zK?(c#@gvZsE)z%gezs$xp`WfNuXU`tMCmgHRtrT&Bi8GN#px?c*ZqDOdyTylTTD^fcPJPD z4)+;;!`)A6>2Lvx18%2{5i)DQK6h{+j9t<&9?>mzUQh&&Gkd$DqC&Jw7rCtKlIDmY z6BAQ3cw(2B=1ITwTvS`_2QEdLVXsYmDJ*@E-;jsq{7jJlygwvL#;ZhXdpryWQNdG< zxPv|wP3Eh=kaY~=UE;{&KEmVZfS@@AdYLUjjUoUjR94$Fz-n)99b~hgZac4AWcPQ6 z{;D?5eRzQPl`B{FgFJXszzgtd#{pH2+PiY>dGTX}jH!h;T!lt0l`RYuD!}4`VWKg- z0DXooX;um@_PGZmAu`sTsfSZd+kcD)VvN75Q>~Pw%D={G$kKpwr?0e8f29%dfe`R@ zYrHhk@Ad1~Y4Ak&=p2P+`D-{F$q~yZcTI_>pU$ zKY!-N8@K@K%KGtj0arO1}a_MNqI2|)Yv<#VQxi8`3m)z#JY3SOc- z1eTa5ql-CXRpv-dY-}tMA`_Pj7cN}a#2BJVBH^XvCgIlU3fs58jK8-;37Y8KyVtNW zJ^e(_&hC4ou|N?xaQCi_EW`4L>k2STJ~xAInVk@dR#N8ZxWJ)8Cw_^e7=*49geT1c z+}dd_36f_2U%Lz63`(7I0;E3g2&3Ba!O$)dJf|c>0)ulkJUu<1dN290?^=TK%jSZ% z(9qCm(y|~U93-uIj)P-$vYA+voP5!2cWe0?$K#$k@fcC-lH7y@<3z-0K2%X@@kNXk zc%-p0jn#3$y^)}H2XecRXJFM+6T&`y`h<)fzM@_nyh#z#bFGk5KPMnE+lwoU|8+I~ zk^Dv#ugk;XK7kWzsUt$*g$$kvFn)8#Q-C;+H_v%+3+}sN;juAatkqf`jwA|nvJJ}t z3L$tzH9_r{8KOOTd3n6zd;)0JI z25-{gd+_M<{rmS{Su4TZy!k?RopGP1HFVqMtfVu~&1_LtaI{*%cN^#`%Y)bOfj10tP z(&_F&aEDI{Op4ycu(Z3amOGi5nSy$SHlzv|B|rR7xrc4{DjBXNDZyLeV$~$TKke=9 zUG75|Q1!x8Hco+F5a<~R4XB6&>|FU|}1)UFAO|-kgr}y*YE5A#{vQ8}oe*b2>&| zUkAvS_wC#F9O>E~K1sk0R5;<$vg3lF=e7`dg`1c>yuljNe8ypMP}~69-qA6d0+aA> z4}WG(PMTxl*y$$ht^;tDlYwNBZTcR(T~j!HMQt&+Ri^Z&(;?&UJ6Os#*7jfIyZ2Li zW~Oqiibd+|*`LNvlTFc&{==tqa9p0w=MNuFCUrf1=Txl!I?yQO-Mf4sa(c%`NA-tV zGz|@Hu0J0SXn&etW4K$={v?m(oxI26$B71A7Fzf3e{fy4I^_3@o`uCq>mxFEKYuEG z-ak?HT~!hHJ|-sSfts3eOKYpD^1XW}=>(k}(xd^6#((-$GC>@skh{Npxln2t4OxkQ zY|GcX%eS|;E#bf;(UL|H)4CIPm6e;tVww*EcHnaN)^3dC#}vU0_0W=%?Qiq9Umpvt zs;WAf+o%8iS>F#&X<=k?8TFy0aUwkysa35l0^<%@s^KORdgLhL^jh8Bxl9VLUm*V4lH7?Q#PEJknH`dklo?v3q zjfLK%DkQ^lxIb_!f8oQ-G-*$5c{)MJsD1%0_f)pCv)fMB%IH|%`UVX_9g{((WVIEI zQ*k<^H99(4U}|dmY+CM7JFCgvyLU(5pTGT7TI-Cdg+=wqL8wioa=FJ%h%&0lLBx@}G94uSqTz znFKaYtkhYU6%L%@R*l!McvShi>S`Hs4v70BybRHrQm{NWTv1TCfnWmFM74co-^$7c z-@^#bfA$@DT(i5o`==t}3Mh$Hma+D|!f1Uyoq1qPEnq_H)K@UTJgOiXJTJh^q!7B*GN&pr%3pT8E7U?|BAw5 zu_cga)#NqGgXuN{SxFn+9=!!in&3rX??e{Xe?X&U;~K(`@4cevLmv_=Z$0dKBaizi zSHM2JB0nRkO0R=1%^v#5#o^FD%%>SgyCkBr{o*+|PY9*UxGrC|ynXlXo10U$jOTCQ*WGnb++{eE{FRP+J{&Em;i04yMRc(Dbr1+V{zt3<#)td8i?|SYwjBowcPvh!b&+oK^D~0+4rR_^?}WWZ+Rc zI+L3>Z+=&;HmBPG5%IM|vcv!vIhK){AV7NnEVT;MV3&Vt)`Ad6+ zt-g>oP%G}AHkMa3BL~wfbZxOS?U}y5`?-n-z1R?g;Wat|&k6z>=dThy7VP#7%*Noi zBbd|&;jYBrTLJ5jZhGTR(Oph=BJZSe4CEpMgWkPa&B;Q0{iNp9fwdc-_q)D$8-0$o zM3VWz-x<-GjgBvhH0CR?g`GI@Y1E_H`G34ly=BY&_Rh&qZ$}fwRKOe z*Z;6M%OhlpfL#-_*ovLa(`V+@iB$trs%KYggLSX2tk0l-`uxhN_hIN=a88;%MhZ~x zau^7VKiS3ZBX;oML9f2HwgiNWVD#~hUuu@3tIv>VpVrZwQ(`+cu9GcsYUC~$pN866 zh1ay~E$vMw=+>*lz;Wnm}AcE&Y+cmNb zA65oq2>b!HP)3);U8mPW%t5bejE!^qu6whheD>LZOy6q%duj}BDFe)I!z%6W_$?j* zA|Kzse^1epZh6sUEZ~_Y0aOR-Q1Whis|ljxm`5(nW3;vzp4E9-<0Rx;{t9zje42W? z!+zi;ZWba)KRTVjBy7>bqrUcp*c2@sSnAHn!SS1q!O;eB)>_V8Pw&D3q$!qC3 z`MhefrsOKYmBuemj+~w33dCr zJ@u_AAq;e(9K?EE4vw@_S%o$JObu{_fr(<&Y?BW+*n7kG5fSFBf{$M|!r8o=nkspb zhn?xM3ND1SE7%{oxw+*eB~VT>&-_E2JKg28LgRBiA2@Hz$$l}Jd>c^L{D_mx$MGYZ z9EeemV_LdEL+{ww+KOOhQ~%ksRMkJs8V%|A84g0_NA%)QYfDSZ#`AU)qY}rNtBZVh z6cuybJtmc0JfKz`mEfUKdqC{K@m0$o;xSc|Qs`XiI@n)tj=43@AVcD60yXk%Q=w~k zNUKrQ;!F}A#Q1l%S?N7+x;&fO7wK{e(XerKaK$ztD*@4p;a|P`VuE+(s%`N(L>RZ zhX(rbh?kE~rDNrqUcs}wEHE`%wVwwFrG$#Y2!5mZ;Q{_jmoL}1xm^zS0{q35p`F># zwaj!Pw5?5z22x3fu)!wbaGA9-lpI7zqv_IzFBc z$tIR_JYHK`Sg`8;rkl0CHC$A5y(fZKPY6Du5Dhcr#J7l8c3iMhB^EoHK_ZU{@J0>9 zkH1I$wl>HhzrYu-{cC}P5J=?-3VFcKF<0t=QF?FWHT!ZdGcDb3zXU=OA~hoxz!G99 zyD8Vt2HVlur?O>;O#lkz(<_+&1N#;H(3H|R-KzkTcm?VF!(>N9JR%|^xgo^17Ds@i zM)4vt2s`(a17i8zWiRlFiaIyGprWqQh;A}>@b>nWAtndZ0Xm*Y*%SsM66N#4YWYtc zF0}=f`uhWW(bvrnScj>WySlo4rd?y2NZ!%|LK#asIy$cBW{N)Zh^eqI*0=ZK)sU#F zt}X@$^`Y-s3$8G7 zd3l*ZLCN}m0SFL{;Sw&3v8Us)90FC+0uGXzv+y7E5Z6OR#ga!lK^_i$h= zrIb(L93UB9@CrhV$Fw6G=`+Smf{;*~heI`$jnFlHZ&Qv*Y3p0wWgedSE+g-(Dxd&E z(h>Hx&ERv)=H%KabAl1V$3qBtnTxCX2ZW0cq_JUuIUj0htOlm0rv7pX7E1WppKmoK zwbpj~g*%%2;>DU|cYyjg(DZ>PZ3&Czh&!X0Y<^Bo7PejmJ5F!MapDlO=vTQA8vrYK zcXSdGC?jEbD%(|o2W8})IW&Fi+1m@j(xRf>>3h$5hhRdnC~mWWJC9y0M2X{OS05mF zI8p6p2!!bOCK`g<7CwTute@{PIe5qq2t7-8_DCKMAn{Hg_2Arx5-3t;Wcu^rG}A*K z4hQFc)r%P|R7gba9{%CO{V~X45HGVH_viEkGEy~zjbETm8p>+;-;#4RE{Ka4KUjYN zizVE9juSCPZv#XM3f`&39Yn0&+C9h!R4}lactZ8tdyQ_c+vf(8|rj)eK-uJtlrcemzVeF zDq4>EkRKFM3`1#C2)P8txQ3;*bh(D*>K^fG4+tdRkQ6^Nad!0!k)cN68}^zFWHQ|L z1~nTY=-)A4OaG-9A`B_#mhE4(-b+)U=#E_`)H}I|rsnqTHZ**^igkeaYsh`kx)8td zVm5~g<&#nQNATA_XQP8&1l&4UzCl0=XeXMdu8D{!k-D};Tp{}U`o4kEy%*`|qa#nvfKCGT$SVbg1~bv{ zu~iE?wCwy^JxTsrRc)=L76ypQ8)#9B&Q?tel;)~BD~hCW@7}$4^`<)tW@Z!DZP5c0 z&#M=Eb2yz5`t&$0Eq|Ptt)bRCZTPf}QCHI%HT@Hb24PmeJAe)!9CIe=i z7e3w)I_}CNmKr$0$cPBhKA`GS&kp<1$`r|5y{c%2FjO}HBejkHqGsP+xJ@fuKhxPK7SfB`=Knvz)-V$+L3q*2TyWx zu<81Cc6TScKZrZVDzUaCy*aolx|hH^Btb*}1SIprEl6H7Kz>I{VWvuBSQ8RDy&-6F zk^H)=fPgr7GZkOzb8-9dH!wX4vGD+lYYq3 zS_U?)^p!t!NA1at$W0n%rgNcTm1Ylqem^R{Kx8!-qDCE}qM_FhMas)Ux`F;(KvI7Q z=Pp51#q4wr|M)}w^TsJ@lUSu)8^@HMr<;c-zSOU)BYD+3#Xc z@^wzBh9-pW!w+Rw%ofdimwv0|u-MtXd{Xq0W=c`~nl!VhCw@2qE4 zKJL?RPdo8uyT|?{rdBx;Z?z+>!ADe*I6G)wIxm9~V1%_IIv1{7N#_N4u>DSK=s~K) zQB?L1e#5e@Ybxve1Wq_<`NqzC<~bZAu=As5HY3<*<5M<8AiCW1b(5h?rFC?NeJde5 zwrir?;_W>xzcoWo^hMppesFeb>@X7KH@UR^S?6iy~va-#ZjTjDQ0;4 z*LZ?|vOHPd>C6pCXY1!c_N{g<8&_cm{IaB`#l8{x@XM$}c{Qg}CdMt|EMCEsl=2>{ z(HFY?rQNNhF`6qkvh{tPV#HRFLxs<_PP3U@2UVxL6aMSi{Jo|iduhW{%yfep&ks>m zDko)I&u)E{IGS)$(R+mM$dMz#pb|ldN_c_(o&}>e{aWbY&jdR=J*SlWzk4uIvNisJ z);)juvT=@Q?k;#MMD{EU?7oTOzrw~Gm-5vacYBs%uRoi?eBlDegOiKVJmL7nY8I;R-f;j!Y<_XW!T^AdG$7%^1X3GAdOPH%tz2Oa?$*XlfO8=|3yFV^nL&ccn!H`;oSLzQlBIYl?4a98Pr{uVVi$u1(oKaL+vz_BrIP>;Sws%5_fXwKy+&F)_FjYKYPgf{T z(HsB+j#=q72`RtsINesM{#Uw} z0Z~Aj4OP}qkfEJO$zLU(Zx^mHD6m}kTFW*dGrMhZ4IV|U z)S2ecX-J0EKtT7VP8CUZ73qI} zNSpcuBE;3V!n51iyR&1}94kL&rF>EVj9T?$8se^U6rTy?%1ceXRQ&B5Ps`(+Ds~CS zY031}eSYt~pzJEbvNM_FwK86B24zXT;KsYO)zlIs7*P7xp0PXeY0AoSqy&`xBn|gt)w| z^Wa!9M6yIVCbd}uP)Cs|8pAsYMPwv%GqbPt;OWy_jvTuH{JjL6?{0vkL?GE-@a7aZ zTM*p;kmpQSJ8C4o%XM@&r-p%lE_w_9M7*fkx$ikvxY?pqT3dH~_kNV3YIX6J6@T23 z#KfqrACW0!6zd>s%rghk>fTu8_6T>KC?59799 zv)F8X8{0&3yBrK;*8M0~AM2h+(K48~$4Z}c_D(a6Nmy%N9VL&)w*`zeW)(go%Wha# z?7%*Nc>={zh~2E|!mGD2I`ZwZ%3K!8b5-}ODmlXGx3jj%KJrxVv#<9=V&Yf%YzwvI z8~h8z!_uE28Goa+pQIqkN_$pZyiB8A&8W++k0&h9Kmr(-FD=bV!EbVtv}lVvACyX( z7iYvasWIupY!BM+S9=vumM(EA0YCykmjr6+LM@jontXBje{G!ax$wzc+}Ev z?!B9wR2R%^iW)cmulBAqs;MLk2N!4-mq^5Ri8!b%BB+Rf5dkHL<3dn^2r2 z^={o;b?d(GyHzh(66PDTOqqJ&bJ5W{n_p`Sp0Wd>ga8`gv3>fYHAxg&Nd@@jsZOzg z=u%yk)o2~%PdaMV=I3b}5E<4Rea=24DyuPPLg{(_+qQ;(3c)tVE~TyHYnn+cZ0#;H zyT6rBdwh%WI5;%lw9qC$5q5$Dq4Bn-=5==gSFy(cGm;3E?v}2z?{|o{PzdzO@rptB2 z;iqRCUOiFKzo@m|v$Lpl_?0p~B4$sPf5qj`N%j6y3GaS*cV&I zRwO~{M9?iYTNB?!3rGb*%)k7|A??_V^}~zHx#}00k6!v%Y30oJ6z)B-ebg&j5x#?- z8h&Tp#HJfdB?X*YanwiR7IYh7BZWGc_L}*HkBHK{k}!^C&*PYO*gP*CA)f-?v((+-{pJ z#bf>M&92*2Z$sF6th0Hxrb;EWNim~X=W?gM_(F@hJ8HHH8cL*1ro(aG{u|#LX&X2E z>Kt=6<5CyA&f40UJ~x!75>~a50-X1>^NjHYnN3R&;s9oCIUK5F;%jGl_|um#g?CZv zTvte6ZHwL{=!APa_}qqBVxN>4HPB+fWDu6F9_L50`Wj5}c&8torH8!L6gRO?oc>5) zPVg3S;>0NM@haphgCdw3Ry9{<2c4S~fin2X_Y55YuGyk0L6D$p1iZToXRvSWI-vq7 z!{g&4Wuw*z1V16y?tj{p-nmsFXh)KmdFI4|1gNLP8f9g7(+3K|xOp88)l)gU0s5}U zCn&q!rF1vTR>h@}wjrb*N#ZIA*B4K7gc`eNR+*VAiRp%5L^Q5NNxFCA>st*ox!t!U zmpljp84A4Lk6Uo|i4jfV_yrtpdBjA$b{1SZ>hU$9P$@jrsx`O#ty;jhU6xG z{~NZg-cbbgr;DH)bZ;U>j1hccBevQ25D!H$~vYDx{6E0V;?eyzGGhmvRC<0EG3Vz+;z0q;gi_b$uJrNM`6{q81Z@AtA{ zUquz7@JC-EZ?7yHQJ9c>I83J9XC%3xd>(v8S6h@GWfa6!mYH~&8c59dRkg&0HCICS z%AE?+uZR38Qo-<0YM7*pDKJ{|Wu}SE-6S{9y1}mC^`}ak)jeMNZxIZhDI80@KMMTU z&47u$t^NE8v zA_TBzmGWLxxjTUgQU`7a4!WqTge;Zxn~rJMRyU7dIjSko!A2c!uJu~PJ~iKuCTj5{ z#fM3~#lX&<#jAu#O6H@q*FJQu&s4Jr5yx)ewLduz97+~2GzrRVk{Ag3T`|>5C*~dkpX3I_dQGrix_5FfKg922EA*x! z;+UIFt(t3eRTQ{78=_BILGZj+r7+&(U9HS{qHkISQ1K;TDD7Jf>qOLrZS&vo5#{EH z5?L(@DUIq0UXjc-%=R%sseVSeF}eTR;^{f{dR-?^&x1YBLvx&eb|vbN4883Cm`nFK zyGpJ-oqm^G)mw~`>>nUx5+Z7*W$S0%;cZU*;){U=_&77FuSe`Rewi7vqOrE^@uvX; z89SMx_t|7k`C%_D1Im1CYXq{!S>d=NB!vtHHQte@TF7Drtst(La;o>@G^M)UU@V%R z@k@eG@q0LIM26RaKpDkgAchq}2KUv!-e^(CQX_FX_V*j6BU9+TZ)$)?S3psGK8Wz1 zUVQpewIg6p+MIoAVjj&MJPNEL-~_f;c8 z-M361eGB4>7BSJ0PcTmBwfSRE;0!j3;v^0iIv-K~yO19h%O7yO0EuwAyYq*aW|DWR z^fYWv)^@rETT*{ph)S(dAKacST`%-_2g6NR#^bLwUYsL`ms%;K^Pj`=f7aknsri@C zpL6G5cKZMABnsH#%!|U$6b@{AV6Vhx*{DKCft1fcya3x#Vu2l&i6+JxbU-1FEnT(i z@|D5waNreul7n7uoG-hHUzG MBNKA&ew!=*297huFaQ7m literal 0 HcmV?d00001 diff --git a/scGen/models/vae.py b/scGen/models/vae.py new file mode 100644 index 000000000..5c41fe832 --- /dev/null +++ b/scGen/models/vae.py @@ -0,0 +1,63 @@ +import numpy as np +import mindspore as ms +import mindspore.ops as ops +import mindspore.nn as nn +from mindspore import Tensor, dtype as mstype + +class VAE(nn.Cell): + def __init__(self, input_dim, hidden_dim=800, z_dim=100, dr_rate=0.2): + super(VAE, self).__init__() + + # =============================== Q(z|X) ====================================== + self.encoder = nn.SequentialCell([ + nn.Dense(input_dim, hidden_dim, has_bias=False), + nn.BatchNorm1d(hidden_dim), + nn.LeakyReLU(alpha=0.01), + nn.Dropout(p=dr_rate), + nn.Dense(hidden_dim, hidden_dim, has_bias=False), + nn.BatchNorm1d(hidden_dim), + nn.LeakyReLU(), + nn.Dropout(p=dr_rate), + ]) + self.fc_mean = nn.Dense(hidden_dim, z_dim) + self.fc_var = nn.Dense(hidden_dim, z_dim) + + # =============================== P(X|z) ====================================== + self.decoder = nn.SequentialCell([ + nn.Dense(z_dim, hidden_dim, has_bias=False), + nn.BatchNorm1d(hidden_dim), + nn.LeakyReLU(alpha=0.01), + nn.Dropout(p=dr_rate), + nn.Dense(hidden_dim, hidden_dim, has_bias=False), + nn.BatchNorm1d(hidden_dim), + nn.LeakyReLU(), + nn.Dropout(p=dr_rate), + nn.Dense(hidden_dim, input_dim), + nn.ReLU(), + ]) + + self.exp = ops.Exp() + self.randn_like = ops.StandardNormal() + + def encode(self, x): + h = self.encoder(x) + mean = self.fc_mean(h) + log_var = self.fc_var(h) + return mean, log_var + + def reparameterize(self, mu, log_var): + std = self.exp(0.5 * log_var) + shape_tuple = ops.Shape()(std) + shape = Tensor(list(shape_tuple), mstype.int32) + eps = self.randn_like(shape) + return mu + eps * std + + def decode(self, z): + x_hat = self.decoder(z) + return x_hat + + def construct(self, x): + mu, log_var = self.encode(x) + z = self.reparameterize(mu, log_var) + x_hat = self.decode(z) + return x_hat, mu, log_var \ No newline at end of file diff --git a/scGen/train.py b/scGen/train.py new file mode 100644 index 000000000..21031c382 --- /dev/null +++ b/scGen/train.py @@ -0,0 +1,37 @@ +import numpy as np +import mindspore as ms +from mindspore import ops, Tensor +from mindspore.train.serialization import save_checkpoint + +def train(model, optimizer, data, n_epochs, batch_size=32, model_path=None): + model.set_train() + data_size = data.shape[0] + + for epoch in range(n_epochs): + permutation = np.random.permutation(data_size) + data = data[permutation, :] + train_loss = 0 + + for i in range(0, data_size, batch_size): + batch_data_np = data[i:i + batch_size] + batch_data = Tensor(batch_data_np, ms.float32) + + def forward_fn(batch_data): + x_hat, mu, log_var = model(batch_data) + recon_loss = 0.5 * ops.reduce_sum(ops.mse_loss(x_hat, batch_data)) + kl_loss = 0.5 * ops.reduce_sum(ops.exp(log_var) + ops.square(mu) - 1 - log_var) + vae_loss = recon_loss + 0.00005 * kl_loss + return vae_loss + + grads = ops.GradOperation(get_by_list=True)(forward_fn, model.trainable_params())(batch_data) + optimizer(grads) + + vae_loss = forward_fn(batch_data) + train_loss += vae_loss.asnumpy() + + avg_loss = train_loss / data_size + print(f"Epoch [{epoch + 1}/{n_epochs}], Loss: {avg_loss:.4f}") + + if model_path: + save_checkpoint(model.trainable_params(), model_path) + print(f"模型已保存到 {model_path}") \ No newline at end of file diff --git a/scGen/train_model.py b/scGen/train_model.py new file mode 100644 index 000000000..d3ae43924 --- /dev/null +++ b/scGen/train_model.py @@ -0,0 +1,36 @@ +import mindspore as ms +from models.vae import VAE +from utils.data_utils import load_data +from train import train +from batch_correction import vector_batch_removal + +# 设置MindSpore上下文 +ms.set_context(mode=ms.GRAPH_MODE, device_target="Ascend") + +# 配置参数 +train_path = "./data/pancreas.h5ad" +model_path = "./models/scGen/scgen.pt" +batch_size = 32 +z_dim = 100 +lr = 0.001 +dr_rate = 0.2 + +def main(): + # 加载数据 + data, train_data, input_matrix = load_data(train_path) + gex_size = input_matrix.shape[1] + + # 初始化模型 + model = VAE(input_dim=gex_size, z_dim=z_dim, dr_rate=dr_rate) + optimizer = ms.experimental.optim.Adam(model.trainable_params(), lr=lr) + + # 数据预处理 + data.obs["study"] = data.obs["sample"] + data.obs["cell_type"] = data.obs["celltype"] + + # 训练模型 + train(model, optimizer, train_data, n_epochs=100, model_path=model_path) + print(f"模型已保存到 {model_path}") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/scGen/utils/data_utils.py b/scGen/utils/data_utils.py new file mode 100644 index 000000000..284a0f81b --- /dev/null +++ b/scGen/utils/data_utils.py @@ -0,0 +1,12 @@ +import anndata +import numpy as np +import scanpy as sc +from random import shuffle + +def load_data(train_path): + data = sc.read(train_path) + input_matrix = data.X + ind_list = [i for i in range(input_matrix.shape[0])] + shuffle(ind_list) + train_data = input_matrix[ind_list, :] + return data, train_data, input_matrix \ No newline at end of file diff --git a/scGen/utils/model_utils.py b/scGen/utils/model_utils.py new file mode 100644 index 000000000..cf7ce24d3 --- /dev/null +++ b/scGen/utils/model_utils.py @@ -0,0 +1,31 @@ +import numpy as np +from mindspore import Tensor +import mindspore as ms + +def give_me_latent(model, data): + model.set_train(False) + data_tensor = Tensor(data, dtype=ms.float32) + mu = model.encode(data_tensor) + return mu.asnumpy() + +def avg_vector(model, data): + latent = give_me_latent(model, data) + arithmatic = np.average(latent, axis=0) + return arithmatic + +def reconstruct(model, data, use_data=False): + model.set_train(False) + if use_data: + latent_tensor = Tensor(data, dtype=ms.float32) + else: + latent_np = give_me_latent(model, data) + latent_tensor = Tensor(latent_np, dtype=ms.float32) + reconstructed_tensor = model.decode(latent_tensor) + return reconstructed_tensor.asnumpy() + +def sample(model, n_sample, z_dim): + model.set_train(False) + noise_np = np.random.randn(n_sample, z_dim).astype(np.float32) + noise_tensor = Tensor(noise_np) + gen_cells_tensor = model.decode(noise_tensor) + return gen_cells_tensor.asnumpy() \ No newline at end of file -- Gitee