From f4dc9f7783d2b07d51d48f5669359e48049d7791 Mon Sep 17 00:00:00 2001 From: He Zhou Date: Thu, 19 Jan 2023 11:37:21 +0800 Subject: [PATCH] First DOE. Only C is significant. --- anova.py | 81 +++++++++++++++++++++++++++++++ anova.xlsx | Bin 0 -> 12796 bytes env.py | 31 ++++++++---- experiment.py | 67 +++++++++++++++++++++++++ result/experiment_result_25.csv | 26 ++++++++++ result/experiment_result_25.xlsx | Bin 0 -> 11647 bytes xv.csv | 4 ++ 7 files changed, 199 insertions(+), 10 deletions(-) create mode 100644 anova.py create mode 100644 anova.xlsx create mode 100644 experiment.py create mode 100644 result/experiment_result_25.csv create mode 100644 result/experiment_result_25.xlsx create mode 100644 xv.csv diff --git a/anova.py b/anova.py new file mode 100644 index 0000000..0a0f382 --- /dev/null +++ b/anova.py @@ -0,0 +1,81 @@ +import numpy as np +import pandas as pd +from scipy.stats import f + + +def do_print(lst_value, str_col): + str_data = '\t'.join([str(round(e, 4 if 'P value' in str_col else 3)) for e in lst_value]) + print(f'{str_col}\t{str_data}') + + +def anova(lst_col_seg, n_level, oa_file, result_file): + df_oa = pd.read_fwf("oa25.txt", header=None, widths=[1]*6) + df_res = pd.read_excel(result_file, sheet_name=0, engine='openpyxl', index_col=0) + assert df_res.shape[1] == sum(lst_col_seg), 'the column number is wrong' + assert df_oa.shape[1] == lst_col_seg[0] + lst_col_seg[1], 'the column number is wrong' + lst_head = [f"{idx+1}_{ind_name}" for idx, ind_name in enumerate(df_res.columns)] + + n_col_input = lst_col_seg[0] + lst_col_seg[1] + n_exp_row = df_res.shape[0] + n_degree_error = n_exp_row - 1 - (n_level - 1) * lst_col_seg[0] + + df_output = df_res.iloc[:, n_col_input:] + + print("Source\tSource\t" + '\t'.join(lst_head[:lst_col_seg[0]]) + "\te") + print("DOF\tDOF\t" + '\t'.join([str(n_level-1)] * lst_col_seg[0]) + f"\t{n_degree_error}") + + lst_report = [] + + for idx_col in range(lst_col_seg[2]): + str_ind_name = lst_head[idx_col+n_col_input] + + df_y_col = df_output.iloc[:, idx_col] # the y column + # sum_y_2 = np.sum(np.power(df_y_col, 2)) # sum(y^2) + df_y_col_repeated = np.tile(df_y_col, (n_col_input, 1)).T # repeat the y column + big_t = df_y_col.sum() # the big T + + # generate T1, ..., T(n_levels) + lst_2d_big_t = [] + for level in range(n_level): + arr_big_t = np.sum(df_y_col_repeated * np.where(df_oa == level, 1, 0), axis=0) + lst_2d_big_t.append(arr_big_t.tolist()) + arr_big_t_2 = np.power(np.array(lst_2d_big_t), 2) + arr_s = np.sum(arr_big_t_2, axis=0) / (n_exp_row / n_level) - big_t * big_t / n_exp_row + assert arr_s.size == n_col_input, 'wrong arr_s size' + # sum_s = np.sum(arr_s) + + # so far, the first table is computed. Now, compute the second table + df_s = pd.DataFrame(arr_s.reshape((1, n_col_input)), columns=lst_head[:n_col_input]) + do_print(arr_s.tolist(), f'{str_ind_name}\tS') + + df_s_non_error = df_s.iloc[:, :lst_col_seg[0]] / (n_level - 1) + ms_of_error = df_s.iloc[:, lst_col_seg[0]:].sum().sum() / n_degree_error + + do_print(df_s_non_error.values.tolist()[0] + [ms_of_error], f'{str_ind_name}\tMS') + + arr_f = df_s_non_error / ms_of_error + do_print(arr_f.values.tolist()[0], f'{str_ind_name}\tF ratio') + + alpha = 0.1 + arr_p_value = f.sf(arr_f, n_level - 1, n_degree_error) + do_print(arr_p_value.tolist()[0], f'{str_ind_name}\tP value') + + lst_sig = [c for c, p in zip(lst_head[:lst_col_seg[0]], arr_p_value[0].tolist()) if p < alpha] + # print(df_output.columns[idx_col], arr_p_value, lst_sig) + if len(lst_sig) > 0: + lst_report.append(f"For indicator {str_ind_name}, the sig factors are {lst_sig}") + + for s in lst_report: + print(s) + + +if __name__ == '__main__': + # first test + the_lst_col_seg = [5, 1, 5] # 11 factors, 2 for error, and 6 indicators + the_n_level = 5 + anova(the_lst_col_seg, the_n_level, "oa25.txt", "result/experiment_result_25.xlsx") + + # second test + # the_lst_col_seg = [3, 1, 1] # 11 factors, 2 for error, and 6 indicators + # the_n_level = 3 + # anova(the_lst_col_seg, the_n_level, "oaL9-3-4.csv", "result.xlsx") diff --git a/anova.xlsx b/anova.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..a9a53ca6c25fba6aa8037b5f1ce00a027fecc2de GIT binary patch literal 12796 zcmaJ{Wk8%uvW4L8?(QzZ-5r8E!QCB#ySp>E1q%>75Zv9}-3c1(LGJFAyxr{|%m7v0 z^_}jit~qrSq(MNTf!;=l+@`?W>%Rxo>l-6mLj?z0J4XhES2NVtKS2I6ORna<=LQA> zvUmpsg!p$e13NoTNZ240wX0*Bv<&lLoSomI=?g{iEQDc zhQ{%rb%Mr;cEt;pVVN8`(f6mPL*C1oIGv_PWu$N%RTYv4W~mDA#1^zOjHMXRivvic zc%%}BCb1xa5GZ6Y?W2JNjX>jr;bwOLRWg&2^~NUQ|xlTfcoefO>QuJPvEobL$oIu^#b$F=UYpz ztL^pCimO{I^DO>;8xZR4gM0z@PoOKJ!vT`7Kx@AOjr4b*4Q(Ba-;j=q7nJK_L=?T1 zSSGsU`;t%sSx%Y+=qFROB?MtJKM)^me=OWvy(bfZFS!g?$I&iz_7|S7u7fHYf zVy{1!a`~=fQpe_G_SmCt5mbSYoB-gmm>mWutu5Q9_)b&l+!)|bqv#+Er&tveT$sF` zG5JSq&stQvl)g4>`i9Hc4Z<-N%i4F9t774+EFP&0H#hNL)qvVJ(YXNhT5e-tUsQR@ zq*s3kfbzkdBoHL_T^)Kb35yEFvy$3l>d{I_B|wpFH{p3rubCXdIPgEObvp@74TQP@ zq3WD&y(oBr{1bO4^OswSSKQ%WaYy+F?#_-*w$^XR7slyGcQK*{pMkvy_j{+uEZX2% zi}MS%+O02pOHNgclNL%A3qPKj3{&a4<)?)FxIAy&>poKB^vzF|sZx1YDo%(PepIja z0Lo`;Sr_W|v21#9{IDcLRgJ1D2g8|!y# zgE^E#zM!Ehihh*2oZ*jVe*N)?$+daWHHUcq7q}*fC1TBow>_g@qFvY`8G7Z+v4RXV zAT$&i@>>}D9mKj?$w`OG>;&_vd)F2)_3)m}xAd&>v+;EXTz1(354kI02O+p5z`}_h zr}_BT7~L-;-~f4NccY)cZ@B*F`~vcdtDDuY%<`WN;QugiwRNy`G&44K`aKK1puEpz zh5!OWzy<=s{C~E;T$p}+Ml_b9x7d-r5n_E9udqdSWl7o0ZV6{cV61!zQ&vo)Do)XR z>v1(!%JLkwfXeY{9D`UTGm;w@f+YxlWjB2EE+c`ixqaV@*GcC^-Qwp+`b6*hBBzTmoX z%f4X|hZ*2XI+?p2idnyVx7{|~*Zj~@%nh#WoB&ulL0_exv-4>a-yjVwbDkkPo2%q3 zG*$y2lYfGGoU8-@;Mt5OaG%O5blE%Puco3!fq|alU#4eH4j)^YO`i+R+pyR5LSzLZ zYY){=WN*~S9&Z)NW-iQl1Q#ATIgNV{7CoGO@By*6Nte6R^R_+58Pbb$3rmJA*URu3=h)J;9KU%$@qi$GDgiZ5d4~xK=v% zdI|jUc2@EkQ4sNnrAmzyzq8iK1gr}^#;|807 zfjN$n1jGx%q&}gjuH9qZHs;^@W+mG0t3i(p>+wTR#Xd5#<5VByv1!%qM$CSkXzo&t z!IprM9@#c%^w9T~83CDaC$pE23RaR&VyXKKOk@!0YXzjzRtr z0Nw2Em&(L%vYUT*3CyKTRSe7KKm8G@5ORpNmr|I|F+z|*RH&$09rXv)4-NnMFSH|^ z2z}x@aDiJ39+Yii;7-A@zG*@Ri644ZjKIM)0Uk%rl)exF0n?et#v<=5=&&>yg8Vt< z*eIws8T+uBshb+`gY!ZQv4*M9&u#JL{bpf$1RE%+Sy|7_+CJ@7%#HgNL?#i+NV|n% zv;*#pp(J&(IOJPHh>OA&ELp~cshmk0_0zkxJ^f(k2DIQ4@GOv%m10Gt==Y0Zt5R~T z8j*CAEWaeWy~jnDA`ivd`Vwu%us7a;p6vs<7Ajf-B3f7Ah5s_DH72tcqO`@>S`rz# zo&NX`_Fl^RxaUJxo__WzTeF`QK@S~_f@76=31_-SX^;f#%6Gn(1xTb1tdPUvjDd%> zPo-;W>JsvbKVhScwDk8S`^;_6oTYhDozp?Su;62}A0ujxLt z@{KwzWYT5)Fn0W6wK_74FS-leAC||ogdWJRXXA~!rX-NUL7!1;?AT^SpFGl%wy?K9 z*jK8aMnQ1F^pmnJ$`?DDwfm6b8wqKFTe6WWCA1*z9a!ud5jt}azKVJXq31C@05HM$ zLvfG?FF3pLTG0y&7_TT2sL&R`h5EVYBVsLIvIcB3>SBskBvW*@nmX)i4Nqg_aFPi5 zlQ%0r)(_>f zA(%fIQQc1$H^9ovGF!?jG+#R_Uh5YEAI|^=Ichg!xHBRD8R;pfErxT=mAj{k!&QdQ z((a{;Ylkl#x2)q{>Gi=?5#gM)X>O{sp$h_8j@pgFx^s}aU0&~5EIi0u5|7LV_y&{ilO4V0?>XCmbz`odA0WQv zMb1UIrBx&dp35tN>dE`OoET;7@Q2@$Qh6d8=)n5l)+~D++(z7c)`UkeE|FRa=+=6n zrJSd=Uq{GhcSl%^TP&Cc_CHh0%?aXdNI%ms)y&N~Wy99q_fc!4ILGp6^pb7swH9Y| zias>m6}m<6tO`}o160GiuowFBHwoKSxW;OUJ?XmK>ksEj@l?%!_I`F`NYC z@Ix_01ggmdy*GorE=0%spT)tPAKBYu}ik`6QDET4m4)lX6VS&UuGEsrS(^ zBH^?+R3mM_v4_dDMxUUlW*EbFA~>Zc`-7=xc3|(A|nLcF(4Jk;2${ zAY}BbP?KA|0#eRRC*5&zQGJ{$v0f zpRINDo>92uT^qmrfaD9xW%^llEqjsIach)00SSg-T%{d-fEYg%ghal*E6h6w*vuz{c?iK%m4RjzsNEM8 zW1OxsE5lv}qxoLS_@6-#N{MT7L+L zxS#ojUgp9P5B9MKQvrAo75O{~^D^rMML4ZrJ$ha}hCy(oW|W1U zTaTh{bi^``j`s6RL7GC!bF)Jv8cA1!ST_vFfpT{c5y@p5vx!u|MUr)i<&T(>JBlF5 zF?%+7xI?mzza*t@CWGGnCZafLn)e}(JvM<7#vqB6MbligfGH#$XI8KPomfc-yHwwr z&8Fxp3S6mKO%=JPG6^XvlAr@N_r5|ZcQ;^?X3th3#}eZEa<2O~bBxczM3DyC)FV`i zF8suWEA=Z_R(jxt*cQZv1z9-TpVUygdZWghyw8-72Xh1w;JHvu0&+CG9735&zMF=K zfSwHuz+al$FPeMX8&C>rS?wa`=Yd})VmZ^G40{^0lhCMSie?Sg`DHRYV8RIH7OD2D-8mKb4* z!R!aJEzo#UxklY3yP!o{WIl|EQ`Sp(;mOkqHWX48bnq)@79XWpQT^BsE5bH(|_aH@Cvnh+%K_dB-wY+{GM@BDLsax1ykrh*% zHLdxO&mKEPiGYPttvyZvS%Q`f0~qS!3MG%18f`Q{RxTS)*PWJQ-F~;jScb$g3n%D+ z|B7tpE3#8G;tuORC&Yq|BFIhVYB1Jvm{1s4d6B+8CFpoT4c4x2H#mDBHr~XcnIPG6 zQ=;3^@aLah-=xEEmX8T=o#y#)3XlZRP6v>X!lAuf3hc#0kJDmc>w@+1g$XL;b?dEDy!i@Eqk8}Kd9(55T_5H4#sR?m+=)| ziZOq2byJwb=?|JV*v#C!!_n$BN+uT$>Kf?`ss^pcAITce|=ZGa$=d#uE zia2mIaVcv+gbB1DWjvDT}ci7o>!cL~t%+z9kobX*jjiok0pHc{+F+s zEWfV6gdv-yH#!FFCj#g#(hj@zp@5RYJmaFMUykGy8(kD>Cus;-pW~x_~8e+?$!?Cw^0u{(#qmdzlKJhDO2SWCVgp*C^4DIQ62i8 z%Ov>^0KaDzT9PN5GkW5AQAe#pG+v69TEZL>Q_Ns1JEG4g5#b5Ibd{jv76gzZn(^Bz zaW?C-5oKn|nVW9c!@^pZRWSGwrB6BfkQ7WfXPyy_I5=w=@p~kYtkYZDX^{YNb-&rL z&Q1VI4hQ1rP0XJFP-6`m6+ZhGDOEQ2Zzx^F}*%4LczOjokCw6#v*717C zRaT9kx!IPUyGu_RK&~Klyzr1>v|hp<-Ep=rDX@{+{U)v!J#Ke9$P=Q@byyPOS?(!{ zfv`-dYV)!%i>66fUY8nM{Cx*|81e^uH_{pY|4vi9(V(w1g_D`FwefE>=-o3eV$^G4 zl#25QOZ2DNFXLa_!->Y4Ee;#5kNBof^{KbzB)QH9x~z0o4TB?9wuN#0tO&=w9#_5v zxt3Phr7cp@ZUIk0q1JO>U&#zmCEm zdz*#aNWu>y5`$YsYP@MFY`V!}ZN6&wc8KH z@OBV)4>c&&&pS{V?$cUqYz&;65|yD-*HWKj132k4>vQO`^N2F`r>VAl?G38PYy)Ow zi#c-QI+0^5HU*l$75qRW9cNO}qH1l-3Pbq7A>X1?s%AmoCg5W%w=4E}rtH0v2CX7U zP0Ux$M@{3e^-ag?aQ&M{j(3`twJW;L=io*cpB?W`gL*ASdTswjKb1uFb;c>mj{o0B6``SByu7=$xk-qzc$ji-l@+qQ z_;46K@23F*y=BpM{`No=t^HI7Q*+@Z$T~^fFd%8xhJir*aHW|p7|bn&y~3!rOn4)2bCJ_^?i}VGsy1O3)Q90g6)I3Y`>Q%c+8%-GDgsN?u8jd8y=2r*+Uts|og~ z?+jDSSsgRc)-E=A=%&V+wsGDVP*8pvVdvhRxGcK8C_(wm@!q&!@;a(Mq9sEuWI(rv z9!Y>@!$WSBu#hSxS2WQJtoct>>6fS#+N>E&^@~D__8GHQeH@MGr{l8;;gA}v=EV*X zJ%(_AB?_OQS2I4xJ2eUY1i`=VmfK>hwAuJ@CuA1@9xYJ|G{J(RxH&Ww$bGStQG-NY zG$DYR0LyP9-T)Bu&Cu5oEaKs6g-_tJ^{Er{ZR-mg4$EpJt2-xY1~ zVN7jbw8`tn*IP7D#@U_yew7|wd@(40SLB>fkG9FseC*^cRV7# zGdTzd1O}wlEMdotxN^`pM2hhSHC|cLSgHeoyg+FVhV%`kMQns`?KQJ<&{A2fd(A{qIH+oqV5k(wPh?Y7zA8}csR8K zS7O8N&IXQx3-W`i>IO}CPS#flXgoK{`oCCI3DiPtH z!`y>}gSyKw(RZ%UhMWmzeQ`joGx+JtfK`px?ub}XY6hrE<5RT@;KO_=NLZ`|%+Ppn zaGI2`Z%=wcl1~%72`>$pRN2Hpfx1S3{#{=I5O;GWSw7N*%EUqtk~Q~3da@FqqR3U} z1~7JY*nHReC%m$pu~@0b3q~}-_V3r3_4R2E5sv9jJMq7drD!09t+ zGhCSyUo$@%-PGV5|IZdb==$59>p4x$Q|p2g3KsK$kx_3*(4d7}a%t1{xr5c&R$6E! zlaZN`R1nq3RgZxeD9WoW*UMIycu|E)GA!bgfdja)eNH{o1p@@$bM|;SQfsu`F5=#` za1&$<=Ok){$ODVTldfF-?@iJM~T4<;O)pn|mOH4h(iE`|yTGGsB z_6}XkaGYlK5kCGCRj1DKq5cigRUbhx06%F>+w@)b|Z%;Nlb z0yQXC(Iu8Cnu$LZ&1ZNsU}=4y?ruT9GqJ7?&TLEmpEb>{Bn3PCEAuWU3j_rB$0g_J zbTi=ecn(w0c_&MqwdsaDSHBQ-YvMTktc%5?Y zagOocYOftx>2Tw_4O&y+Hv6q9eH|$)}`srELGryF&t7PkzOLv`{?z$VR zY_7ij`D=xN(eW; ztW`uPxYlspidmX~=%dd|TF++s2zj}aaAKdPl}5JXh1%}r5a=z+gfuN7co}z+FM}O+ zwA3n3uFlm0mGGy|wKVn%^;%8|G7|H~Rxq5ykxTEd&K+B7) zw6cX;Q?zmASqB6YPs6nmDih}*iD2nnaD#*_Gx0slM3^F+>x5+4K^~L%2ZvP}oKA2< zR1yU4e+=4hB$b_sID$gB(JHgXIX9nHWPw@n@HSR$fnF&$GiFttUmkL>2;8LF024)8 zE?c}Y=BgJkKtW|^nCW_gMW9L-qu;BWGvMo)tVc?3n%*%gSl2}PO6zne?J_a~6yH2o zJyFbfYa!Y~;Tx7qgDfM6lM>T(Y>W!ATQ2nAVb~zxAceZPzh*pXEQhU??OMf}y?eat zh#R0N-m+XhqGQ0Alvh7%bFxV5BvOn#oV1(sq2r1_<`*E$Pqe6@^GdwdAY^)($v1-LmH#Kf2SSN4m z@|5j(EP$lbmqL??;tnTwB;OnDvQ!hh26mXYyyI2jSCu^R+GuXn?qNhP0o*T|ozbnR zfoNM#Ta2We92BXt&W2Uhuh=m3&0%kyJ9rni`Zp;qhp6VYf1!S|FL%b_1M5paY}OS9 z_&ZNDYQoO#?247Oo4uzhIQf~a0Yet&oEeSF-!c162D6K6qGFc6OY zrmTSPXKr@tVPu=_63EG6`l2O;Vj*s#Av}3`Dj=Tk)C3!p@Hok;dFL}T-|4$qxynIp zZXew@=I0|w)wp|alV0BULv3HhCG;PWa1a^_L~egSXvzn!4hH2YHGrf;R73A_6GW5j zWts?puQgD`F_Eo88UE3X@j#SMrI|)F#{MHXHlZIwNRH4Hk~2*0`U^(n4KGC)+DW8# zK>TQR93Bf_KKRu^&^ikJ-s-FKIb zX`;_XeNW1@0=F! z5uHp@EmDq2T(H)zh)0u^2)807j)8?XF^+7jxXFr8bTeuY$}EEvD22FmsG}8$_K0OZ zoDChDai^8p#dpA5QI{L?1pe* zT$453fSj1Z^;LCoQdWyOlz|7?YB9+`2dCLO)P>!u(gw7l=yG7YE;3iinVaUla_rHZ zBxSdi(5u;PxhTa4-z`Zao571tCk#U6_o|Fn=cm;oAC_|ppf8KN9fN*O<{J^=jOEK2 zZ=6KA{)nzvlA#)cCOt`+I!syZ6|FdQQsKq3(U>V_p|Qh|DDD`meYhT{QYy%+jE>SY z345Iy-cE_Z*>rswy-X|L$$rhY&ryhH?(`_h4zy1f zS+qTlpU-~=p>HH@J&USsO#}qFA;ZxQc@$yN#Eux%dC?rr+uJ2AeIm5slkim>Z)vA< ze!M!kx@fW3%&g&e;GG|k)eeHgG50%Y$BzL#KZeW@^f8CV$$+0}%MR_&*izyQ-VHGn z`LuK2uN>SBbp4R$&wjpi1QYY#Wc>(h*-dDh+lxb4QaZJV(0hyDgIuDokE}p%vWP7W zB+LobYQST?fA_qCF?JfW=5zL#kHCpB5j*I#!Ta4@(lIx}tluyS)9r&%bJB&+=ts+m z^)8#a!p6(pUI#NCeTQwhY1G*3p!z(o5p1TN_2-Ax`+&U!Jqn}_g!|M))?kL0bKLQ) z$CIvk?Z@FLsEr0j6UBQkkDHnMuG$oWZ`@fE)Jo#eN811&PW~rCMJgqT;ba906cPn- zN*IQd3qmX_GEmbS#yG8(hhB(++fVFxx|^ zV?qt%S%^3t-y9K_%b2Nwte4)K&G=|xjYW{&m~67!2-3+yP`>sb5g76SQexsDC#oRk zxS&{MBVZ`Lv{a=LWFwU4E(Ymyxcy0J2OuYfX@Ntq%xC>_{1 z+wD1)$t{*#`X9w)zQaSzV-43u?(o>1N1S3dH!FzCh{>%bD86I^OVt#{!;)1ka zZ%qw8Cl_FaaYJcM?9#XK#lA~_TV%VODu|8csm4TPG<}ku(w;)na(oHYwZNNV3U4;D ze>R=9Kkc;IM5KP!E_NXhTE#2PGdMz3Ob4iTKIVFlSpKOL;fbuF&cUPqZY7dIJT%UO4oogCXt z58sTQR$8lETRfP?j*TNk>hR7Erx9S`G_EjWf(#-?#qB=bwdREG2Os?<*Nb`AUYpb%|bS z4@WZt2V8{5B}J3W!5nNQn2;++ z=&@t9QwgAPsqL1GJ8X)sG@!d|b5eqw6ojsb59% zLZ9ETHWF<8*mYbS{bEM(C_mN)^iQ}C?2gwZUf~vh)gmGM0r#6E=A*5xbIY?LAq}B3Epq7Fgbyr3jyRb4L)X z&Cc7sXQs1ASMFtB{$S}EXG!OCW5C z-fYpVhYxorL$evhVOp)XZOvV#0F$xb0|Z@BI$(aVuj(?G6Q0L^JT#}0{E=3pLex#w`>S`Z>DKn^l z!J8kn`@my%baBS6wkff)K~l^%D>HkxGURrg0>u9etZZF*yTxFq$iEj@lqOR4&kOsFrPS3 zRuBMivxydS@*9!_uz#~ruo=7*>jmH|$nx~8QipoBfZKBfwcgv>mK}|*^QPqskFn5f zE5~IyTg9S926*Mqw)uXw=VQ~8aJZ!V?E2ra{1==0pWr0CJH>i< zg>m{dt-LC_-qOk+W`95JztH`bFl6yd^ceIUkm^mhWG)O9S4$OSwR8m%$omOVi#h$E z%V(`4qA_^gMuy!#P4G5*LiC|aey)rGcG}I9FzlnkzJnfXCG)wxhqoKJsFjX^*4g-X zU(k@_XPcMc3gadb=cxo_G;wf?-A~1AJ8Uk=#B2Epk_W=38iNW4tGlLJxpG@>pS$$n z`c2Xg)iELaKEfQbUTr?M*E0-}_hA5hT=wt6gNhcxX;(S>JxDv*mHWX>P9#DF^+Mb= z(Sk!%bsh6CkH{YS*iJNTlBI6< zGNz@(3liT%c$Q)@FmkJ9Y8eR1d4PcxWfjupLGk6A%abg%aZnGd!B0ZgXb4k#)G-h zRt%^!oYa6u)|IhNMS9EVa^2E)5Jn;=vBiEtBgBELssc|)0+A_M)b9;PTn5>rgc*PL z;pqFgjLkQ6zu-@US+OikV5?^?XN`5_zED{Tmem{UP!3F0LYJ)tYnBxF)J6ka9tI?P4Cx zgwj*Rh1n4N1|d&@tTi>^jtC9ooAcR8+|yGTzDM@~@-sBO4PkgTAw;wC_pOFGE?-K= zcX!55$5V?Mmh+z$6FtgA%}3Y!=!6!%%gooVE908mDVN8^UpoCBXZlWj9@a1H2|m3_ zNrB&?{<<0eP)q;5ORpOR1Pny@`m>7WKknHd6Z~gUPx0>mdgxb=e~q_?)(8EUO7-n=|MQmo6D1Yuzoo9< z1N}Y6y)E= uniform(0, 1), self.p.is_FH_ratio >= uniform(0, 1)) + f = FirmAgent(self, self.is_RH_ratio >= uniform(0, 1), self.is_FH_ratio >= uniform(0, 1)) self.a_lst_firm.append(f) + self.running = True + self.t = 0 + def update_e_revenue(self): self.e_revenue += 0.01 * self.e_revenue @@ -79,9 +90,9 @@ class Env(ap.Model): self.update() if self.t == 200: - self.stop() - # self.picture_out() - # pass + self.running = False + else: + self.t += 1 def update(self): lst_salary = [] @@ -123,7 +134,7 @@ class Env(ap.Model): worker.working_firm = None self.a_lst_firm.remove(f) del f - new_f = FirmAgent(self, self.p.is_RH_ratio >= uniform(0, 1), self.p.is_FH_ratio >= uniform(0, 1)) + new_f = FirmAgent(self, self.is_RH_ratio >= uniform(0, 1), self.is_FH_ratio >= uniform(0, 1)) self.a_lst_firm.append(new_f) else: if f.s_profit < 0: diff --git a/experiment.py b/experiment.py new file mode 100644 index 0000000..0daa25a --- /dev/null +++ b/experiment.py @@ -0,0 +1,67 @@ +import numpy as np +import pandas as pd +from env import Env +import datetime + + +idx_start = 10 # first index is 1 ! + +n_sample = 30 + +df_xv = pd.read_csv("xv.csv", header=None, index_col=None).transpose() +lst_xv_keys = [ + 'alpha', + 'percent_search', + 'is_RH_ratio', + 'is_FH_ratio' +] +df_xv.columns = lst_xv_keys + +df_oa = pd.read_fwf("oa25.txt", header=None, widths=[1]*6) +n_row, n_col = df_oa.shape +model_para = { + "n_worker": 1000, + "n_firm": 100 +} + +lst_op_key = ['out_w_avg_salary', 'out_w_gini_salary', 'out_f_avg_profit', 'out_f_gini_profit', 'out_w_percent_hired'] +lst_2d_op_avg = [] +lst_2d_xv = [] +for idx_row in range(idx_start-1, n_row, 1): + print(f"Running the {idx_row + 1}-th experiment at {datetime.datetime.now()}.") + lst_value = [] + for idx_col, k in enumerate(lst_xv_keys): + oa_idx_level = int(df_oa.iat[idx_row, idx_col]) + lst_value.append(df_xv.iloc[oa_idx_level][k]) + lst_2d_xv.append(lst_value) + dct_merge = {**model_para, **dict(zip(lst_xv_keys, lst_value))} + # pprint.pprint(dct_merge) + lst_2d_op = [] + + for i in range(n_sample): + print(f"-- {i + 1}-th sample at {datetime.datetime.now()}.") + the_model = Env(dct_merge) + + while 1: + the_model.step() + if not the_model.running: + break + + lst_2d_op.append([the_model.out_w_avg_salary, the_model.out_w_gini_salary, + the_model.out_f_avg_profit, the_model.out_f_gini_profit, + the_model.out_w_percent_hired]) + + arr_op = np.array(lst_2d_op) + lst_2d_op_avg.append(arr_op.mean(axis=0).tolist()) + + # these codes below should be outside of loop. but temply inside + arr_op_avg = np.array(lst_2d_op_avg) + + df_final = pd.concat([pd.DataFrame(lst_2d_xv, columns=lst_xv_keys).reset_index(drop=True), + df_oa.iloc[:, len(lst_xv_keys):].reset_index(drop=True)], axis=1) + + df_final = pd.concat([df_final.reset_index(drop=True), + pd.DataFrame(lst_2d_op_avg, columns=lst_op_key).reset_index(drop=True)], axis=1) + + df_final.to_excel(f'result/experiment_result_{idx_row + 1}.xlsx') + df_final.to_csv(f'result/experiment_result_{idx_row + 1}.csv') diff --git a/result/experiment_result_25.csv b/result/experiment_result_25.csv new file mode 100644 index 0000000..be7782d --- /dev/null +++ b/result/experiment_result_25.csv @@ -0,0 +1,26 @@ +,alpha,percent_search,is_RH_ratio,is_FH_ratio,4,5,out_w_avg_salary,out_w_gini_salary,out_f_avg_profit,out_f_gini_profit,out_w_percent_hired +0,0.3,0.8,0.75,0.5,0,0,3264853080.7861195,0.07448611065757892,27068180827.948566,0.028471074467998493,0.9926999999999998 +1,0.5,0.1,0.5,0.5,3,4,3637864305.7383804,0.07297697986009385,30978395677.329544,0.025035570583021458,0.9932 +2,0.5,0.2,1.0,0.75,4,1,2179839879.06716,0.09840417843256091,14870250970.957703,0.09615406368119288,0.9886333333333329 +3,0.5,0.3,0.75,0.25,1,2,3294909225.056388,0.07458035794885493,27529372104.086937,0.027549550048933245,0.9904999999999996 +4,0.5,0.6,0.25,0.0,2,3,3716906250.022679,0.08529021609092702,30937889107.04533,0.025723945674350865,0.9942999999999997 +5,0.5,0.8,0.0,1.0,1,1,3869912345.136378,0.09546657392456424,32080217202.03705,0.02533111889073941,1.0 +6,0.7,0.1,0.75,0.75,0,3,3323780234.219529,0.07324360974831064,27481709495.610508,0.02827211821311022,0.9934666666666666 +7,0.7,0.2,0.0,0.25,3,2,3890360391.8347297,0.09312614979722612,32025419370.64081,0.027331743681722087,0.9998666666666668 +8,0.7,0.3,0.25,1.0,2,4,3724730633.529055,0.08585608644546336,30966818797.24322,0.02556013680036385,0.9997333333333333 +9,0.7,0.6,1.0,0.5,4,0,2059080141.3863091,0.13011386291850205,11919309572.784859,0.17383571880300247,0.9850999999999996 +10,0.7,0.8,0.5,0.0,2,2,3558046321.01991,0.07852769877145456,29326434059.91672,0.02645991757420794,0.9920333333333334 +11,0.9,0.1,1.0,1.0,1,0,2217056469.2871513,0.09085566892663295,14019663732.042246,0.10640474650259805,0.9943333333333331 +12,0.9,0.2,0.75,0.0,0,4,3288880578.2219296,0.07786092688764147,27583988860.902683,0.028509652392809604,0.9885666666666665 +13,0.9,0.3,0.0,0.5,4,3,3879623567.573062,0.09407352725749332,32151544223.643482,0.026003986467554037,0.9999666666666667 +14,0.9,0.6,0.5,0.25,3,1,3577219455.771972,0.0785081831642713,29903884548.088444,0.02491509062583037,0.9924 +15,0.9,0.8,0.25,0.75,3,3,3710277766.924487,0.08632310052044662,30554824167.139378,0.026183286845803396,0.9967333333333335 +16,,,,,4,2,,,,, +17,,,,,2,0,,,,, +18,,,,,0,1,,,,, +19,,,,,1,4,,,,, +20,,,,,4,4,,,,, +21,,,,,2,1,,,,, +22,,,,,1,3,,,,, +23,,,,,3,0,,,,, +24,,,,,0,2,,,,, diff --git a/result/experiment_result_25.xlsx b/result/experiment_result_25.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..9b07d16925ebddd0ec46da04de5be264a9c0fd26 GIT binary patch literal 11647 zcmaia1yo$ivNrDS?(Xgy+?}8Sg1fs9?ry;~5Fofa1a}R=U4y$nIC=L--uds}Yt8Ig zd#bv2e_d6(yZTd<0R@8rdTSx_8$xf_|1{{=A4ax@iVn7Rj*N=0VrZ{7p#BgWd#3}; z3=9Os4FUv&^mj1>J39t9Ypcw-E@@C^q~H_436evFc|>TURd~*Z7@bC2|4%F2dAWBc zI$GJi5s+P+r(26uw^$~eIofAi))X(}f~yK!UehyWL7{Uk00o4;41d=4S{t>~^zFO7 zR93PhX{?OILj=bqN>^w{vYTOjN*Ufwjozt1(2r7NEsJ>5B;$}BXggh` zz1_a~l2-fivFRS8krhPt(EHP`Je!yGO)<)I_>0JZ_otYzM?4D%58xX? zpu5L1={7+W)KWb78@`-|@@L(+ZG0I$Jr`M=njPWDcugDu@zezb{(i&i&H+oGPu;ci z$-(+Ksx8II-C-1#I-OlS6Rc`KHP(Z2b`lpzt#>V$2UVFP>ag`WTJCeSP32YD%y^dF z`2D-D?`!*4#%ouF2YGj6c?<&W#eBG$&h${2Nt@Wf<>A0)D;H`ag(gqPa03Jyet7zQKG29E5Iszv)Jk?1%-Vrl;G#S%}U zw;ey$E=n21ObFj#{r9v2`U+8$11;d=OL?;CtbKfC6ZN zfUy3L>>n!3e*}q^we2Ps!n053QwY)O_t{c%_@bJ|8T*VBP>h>zEsXxY>C?cvQ$cFrw+`a`qla;EP|Juc=<79>9*Az^R-Ly4OgY-uA< zXnM=_c6H^c?!SD(t2_l%$I<2qF;5wLDxq^eL-$wN3T+(6X() zP(I+dVt;K}eV%)|$tN6jSJEQqn$AnqS!6243^|a`Ff3R4fY!%jY=kvi-kM#g-2r(y zyt}z>2x1CK^!cV-5M~3_NuM~2=VB{5EGSFp=`6^*{lQ;X5F#wZ5;{vp^*ahlN~)gB zfKFa}vSQOW66ss^gN*N@D`ZCiNxBO1A)WpMt?n(9Lrds%8LO}4I`NGG#PRH{Jy&yU zm~*%3l%AGLm7sVv`50hNA;=!+d*;t$D;McA_Y}_d2@991`_f%9mW9-Tbzbu(oTDso z^r>|VNzzpLg;UYDy)=6$-D%32*1G;LvG89L7?f~IMYtmZt~JH7->hKt2&!0tuU^#s)EwYj(1X_2{D()f!Dz z#`7Nh6h`NSV$EQsgdOEhv(0pZ){w{UO_~|9Ke$vH;pPWvZH8~9lvo=SFtGr;m81++ zi<3XilxQZ7Y&MOhXY`;s0~CFhww9q}pP)G()))*31=seAcKG@Am4P1i+O~BkKR1YctF3x85W&6?rJ4fa&WVUuBF4)yVV?bW zoha{8Dzh%n@2D3=nZa<=blUL>SWr_JDNeV{Kmd~?V5mrZ1SFY$x?()wo&tTvo)D8Y zjZ)SgG>;lghh^2^K<)EZgl-$MBZ=}FUC70SG5p*Pr8%%x5fNyR1+hCcTa>yW$^sjk zW#&ugO<(2E9Oc6VYXKFLHd{Pdc;!MQgCo%z#8Y7*uX9nxeXzB{M;LY8QEinCN^L%$R4xP&|?=wnVUK{s(_2U&+*#_DbmfDoML(-)MVgsBA^+UZJug-+m ze~BwJ0zi_O-oSXq=g09A?ZHmqXYbS4-aU~!Pg=K?zyBSKxZ2GImD zxGjt=6FQ9rOEvxkZyq7kY($SD@zuOiPs|#C-Pb5WqDqua{#-*)3tiqQ+}NEto`~&t z-0GaL^=Kd4!)VRb5W+ZR^No}sf)kPN$2{B!vZp9Sm41@A;@WPQqjB5^On)Eb?P-^0 zT~2?xpo)yO`%0*4w2x!Pcs~sjKnM{^?UU(#tkDb*4-mUoH)cl(5`G(;3*azFTK`~b zo{_(6mavt98c^rJ8kyMH-VwKfLOB|$`MpUfCXZ;IpA5o2kA`_5^bRW@0S>UW!O;Y5 zq)C#*m}JzI1lGp{Af`XK{XD1250J1p< z_FxRC&0LP?1-CpyU1R`kuu(X10Z||6Y=Qd|%~&~|gHob?Fgn7G!dGkW-105SM48KV zE&z7bKE)unW!zkmo^lGVQGv!z$}9{24y&C6918XpNy!vc_#>aanY^}yyq3G^q8L#^ z<184<Q17Q@92ui5m2E8bG5w0Ankp0g6akCi8m z&UF@z`9U$58;Rw#C5}~+cV~<(SDK35NGcu*SN%d0J>bfMXoOOCAVe^SHzr(>DaomX zvvSZVgj>7TNIw9RY!4nPrV^5LU|QfPC3uCGDagK+!+so3=uX4TF;Gqm0Am9%rMS_*GIyNSvh`puc+-h^P8bUgn$r3 z->}ahTc?wSSagd!CSqErED@1zvv-tL-06>wRhh0uW0A2yVW{B<$7Mqz{1C?<0J{^L zNf22EW3f?B1ulUj0JB3x8rSRVYV@cR+X#$&Ppl>N@k(tbm}MdbwNlF;R4lpr)w#QL zSwzK#*AUp)Y~u0YK!ywT7T%d4F)Q1;UTWA*OGIKfHo(VYfkn_TLH&Gn?hRplZ1dTt zLVImF)B%uDNH7nl1ema`INAD09{YTzF=FCVzI^jg&0g-);5D*i%uT{*^WZvy#5)-A zO>EG8to$JkDLCWR%B&?YQu6vnx~l#3b+{#DR8fc~)>1ZrK7ULk-XTwHh#Dm7Qdvtn zQPe$*T!=D^mru8Sz%EaUJ~Pj-5({ zuWNpn61dg$iX8|p{mrFJ!p=2TNnNjTIlWmqyp1g)%bmzaB)Tu1L|b>(4)_CPE_5~t z=|bX-;&s|zrq-9S+lp(~#}bE$f@VLBEKy4fm=GN78DbS*>b*lR=JBkiu35 zN{kRs;?i0`W?Uj;&hGp0IF->vFAig4uW2{UYYuT7N>36%*Dg(JcEA>sQQMBt%o7`y z^_o%BIsEC)Zey7IgO!1Yxd0a|@NK! zv<4Za!$RnvGLu8NqzJpawqofM+Z!HbR(6RPu1*&K8R^7lmM}u1d`&^srS2amL`xCy zK;~X^9kW#)8;%4QwsbzT*{SF!8T!q&EC0Wtu)f8;J`AevW$62?5+s86!e=M$P-*)2H4Rk$PT4bIZysIEzR0#3WB(;YUy z`Sr?IC)er{;mCT|e2fuyvC0b3dTPb)<9fMDdcAtlX4_HcyGqb;v?8DI{Tf@ZHTNQ> z$Npyr$Ilmgfz}6A*Yhj8zLccpY{uK>x?W1(y-Wql@sJ+wXV7el>(*w{`*FhQ-5R}T z#Q*m`>a9HUx{q=)GqyJVtvqvi3A1Vg0R+U4`)e2br`R9De{7D9G!<+YInX)~=Us_j zH@1}ZxRB1_^Hu5Pa=71YLI}*$k~LKFgr=mtaThlIfP0z^Qh-~In5<|gKw6x#DYQe% zLGRJ8C8evehGnK8aOY-2L`i(TURF0rov#)GLUTn($hPelOLf0`gnknsoaQp~bPaxz z!6k@?={Fp4lU-QJkUUwKdpObUO{XD47OoUXPH#L4k|s}23cobK8Xo2&W*(7i46Ovb zM^^0vAp~9yRh#=VYeL@5cC8pmO62Luk}DuAA~&@-6ncZz+9!PgoE_#*A(w#;!J9S5 z_HN8L1==Z8a}iX)(1(o-F);)y8HQrv!|Y@0(A+G{iYz8`-A8K!-POk_=Y7DL)95K$ z7nP~j;dkx#E)Mq)oii+7zaKh#)@5V?>#st}mD7oRx_TBbPne|h`?<+-WK*Y9-f|1_ zbVZ(~ohko$Qrfsm6n;QpF69(=*2LR$vqWX$)AbAj zjRi*`Wnv-s$w!($M3!rjKFAixEB{kBCRJ!KUo%?-bj<7*jfhikhI^n1pe*2d1uAWLS@iEZ>)dJC)0WyM)b<@o2(l zr^mg?vmQcvdee1==zxvxlvSVW-S33?f>}Nn=VzGPZkpd-?t3UZ9v|Jgkyy-#CnFfM zAFqE7MpOtqJ+<_cpIZA%h~37Di{@!>LE)c5Ik5)YlT(t~x1;wUshENz_`5*wXTjjF zm2;aIHzL?axxN=VqXPYqz3gg)ANaAz40;#WU<>Gm%M_nM8s1Lif=0uQDJ ze0?vCzLp5X{<=0z5H(g<>Hh2GSgw!9ta!9|v8hxAySIt4C`ZRi2lLJ`_W&6Ir9A>M zd|4SemI6}BSd0i4b1<C* zM2}C+Iln0>2jfq0(R~z*Re@59yl$ciII-sxx#YbC?x$V$PC@Bk{)`FWWwR620t?Mb zgFcd>nvs8sGJACGO;sxHh$d*H+P#{+UJ9RciOr)!h!MXj^ zHHBSSP(KxLh0OA;0Y>L?Y(*gagG#~=qr2|(a#?0hXbU?PgQNC>N)<>H;#q{`TzGia zq+p;3c1qe8Ezl7UC75|hUwG;{P9m&lcIwXPua<*G@4VEUTnt1P$J&bu6UWhz1MIv4 z#jPk1m>Ts$X!cl(P;4k=#U>c-l8VZ$i_nMj3BJCv%h4I9C7x@QdQ zufS~h+Z0ygZjV)jOWnaL7h?6ju&Q19fP4o(G&}+(v~d6z=h0gvz@*ty42qAv_Xds% zD8Wafx5I~>MpCuQwt9wjD1KSYs9}@U>M2pcp8$e}d5D7nCw@>&mc7h;SC)LNn@lDH zofb+HN{dQ+9g2s-F~eS2a^r#+gmh5!sld2%O};PtARqGv?-UYJ?bvJi|mUdM#+XqXd zZhy=a_Ta8oV(GOh$bG@e7o~O92m5!n#UegHb?99U5OQ$c5PA3U)_k!Hu`YICQ3?39 zT?JB_(?#V8^;-3sdUl`FFZd0t=znNFddnN`j;0Nd4faN`xmXzV&$(C?mUPuns-Oi^ z{1EQr>hqeV@p0S8e^}eYZ`{acCa1RsUATw+XOZ=jPNI?WYwfRF9ta5j*IUfd$=%A> z@vUk%t7C1uIEVf$p#S2pU+IlyB9m8)S(IEO$>9%jq)3F#q>l@tnP&x(xF@Bm?}`=6 zA-C*2;77b%jSfCnH&Bfj(cs*J;Usr&O49p0(^2WDH{i+OOwAldt%QuEGQ)o|-nP4h z>2Nrtv^s6W(-U`yh~r^iQ*1ewv+26rFk;R)iOVaf3Ke3ca81R(6m7B$MbbV zy?9sEt`4tVh7c%5fRczH6rg7~sJqy4?Q$_+?qveIlab$d-WKxIxls=bA0gAlRW$Ve zLl+fBX^_t%{@$W?fc@$*Bx(6(45Y9KXTCPh!}b-!65T$|=G=Kp&Jg^O=R)C$UbAXz zZkza5-)8>~@0(piK{xAwr(Aj;JCSMFT&j`pp%=!Yt6oojy86tlf(9k5g4MVsiB$bb zYs#Dv%n#q-!yb}oFQV8sA59YrS74jQABT=)AHAblzzLR*%CiKKggM|sS(<9v*gdpc zeq00jlB48CpWKPi17Qg#?q>?OoK#{YOzR1HI>glTm?28|35aWO4DyHbuT-BPJQ6o) z>a)zhcorm)9c{P`d)3r#B$1;}yUUQ9Hg-0OUgX;K%?q+vuKNgDmd zQ%O9^^t8{es}0Ho-THw9UYn>A0pt+VV*LXkXeoVLwJ&MvdX+o0)_P-9bwim%nWbEw zOD6LkAm@mRH<+tjGf2sLL8F;Y$(_*5AY&>$rKm2P|DK!IbPT4pV$OvH?DM-Exo2H6 zogL@4Ld~0ZYQjS{0{u(_JcSYPi+Bd=rCLo?2DN+G!rt2+A5_&-h6^?Tyo=f)M?`kP zpt$X3kN5Xy#53BUR2R?|FlEbtxrd%>`@tMT2_8V)Nt13~uMg~dcQsvD9pyVrPRnHD zS(#anHk-6<6#Bz!H$&G6=V-iC7}fnF5Ehd1Zu>hgzRDx@;p)(3?4gZ$v3m7R^-`UN zFA%9Z1vVp!?Ui+N3%bqa#!PkUB{Q_-!%Q-o4?NXLqAuL^>R`bsmbCroBd6U0wJ3LW zYJ-S|*7K)>Tyr|FE0#MWx`_T+6d&E=`GIXyYYJEM(q>Cj0lFZ?D%bS6_+W4iA=s8} z>FyL~h)jHMbt99loo$oZ;uO*I_a6m^1Zb-|nlrj8>*<(FqbP<7X~QUVe-2clZD&G6 zMbWab!dys87nPOjADgHOmwg+Re5P$qlV>r|YS~3>sWdMKoHkqo-dCQ zB%J^jnTP(pgaa1f5b6L1U)bP4ur3<=cC?_MWscSNr!!ju%)0e=pk7PVk?1orLfN`beLcH@)6%D|Mgg<0A6DCeg5F%xaKcmUEq07a= zL5L%z7k=W79%+oxUzO6L_s)}If!Fh+%H5~|3l9zLMx?z!$2pU>qfZ1G-KCnBC^m~F=z+y&Y2HGOg`nw|_8arE2z#h!>l(CN9XrB< z^kivZhgXij%v&;ycglODr8YjuaKSUgXK9FA9EswUw^2b*XKBEOM7I#PVMVw$`p=hg zbt?k)=ZwY30RE=)NR`g?L#tB!qw0kI{=z4RP%o~2@)YROL<#;Y8{Rx|vARHtMndHp za~N$v$fDEK)RyDqDX*3em6YcJ4fTinmg;j5DzwS({V zvYs}tFK;LlByPq9AGUtJ*E3_vV}?ZDgX%KA{%Pf=3@G)uzFgTF3y1}gL&hto_^gwZr+TU6bwx$H9Ox|bYx|FVTC2;f?Ck( zr0Y zKrmd=2nt&oirfR^-4GvaR;O;isKOG6RN z@3vE`c_z_o4 zg`PHwNI*q%Fasyc5IWlF7!jZq!hOw0_UkD)ErxhO-&Ef?s8|Cy4b2}mI~^NHoE;n#GQvWeCIsOgH3vCx`~uT}Z-IVq3Pc z%G!3tfH%=ZDL>WE4fYk0Ici$n8^9NH)J!{ey&Z4H7Zt(KqV_HpUaQMFLXyvi*Deww zgNT-5^7Dgn>r?8`x<{agql#SFj71R!ecQ0x|D`1NT58HP8>^7wJ-Z%jyiZ%nD3b~#scm z((ZQBSG|l{aF>ZbFc(F)0-TiqLbd}g?MkO!ZW_w~$+EUtP&O#^ljn@y2~AJ$c6Y^1YS#`+@&uCS9lnsU^ncq*b{AmRlXN(z zFLnL@mN)+m&h67YtfN;Ln_uC?dA+`d@Gr5yJMe$deGO!|zLS8qfJtxrsVHXg%nyAz z-^)FL3~7DMK0?Z9r#EB|rL4lgb@oy;IoSQs)X8gWy>G2C;^K(r&A zATs8JTA$WO0K-R0R??oF)DzBYCdPT<8|4@fe5j-Nu5U{n0Sd9N%6ye9R-}bE zhi@I_k+3_?^K|9suwHV#BDC;8sTg7(yYU=gB^<|&JkRp>w z;#At54*ZmDpwZQK*ys2BIE@0%@&jFVTs!NJqw>PEgqbmz+hAlk8qXc@yoKb0zjUZY z&tFcLBn#od${?#w)Rze*`I8Y-zzb%*7|9T>JKQql*WZ8K07+;_VMWAA*4&Z39kR0$ z1-r^3$&~yl+m9qnWd#(^KAUxYT#Ogdb`>`V@w_aA!7nWDRN}X5+u;bLxdV1O+<*7^ z=3=a}r?bpj>lJq888UiAMcu9=GB9N zfPsi!zlSJ+fc|Bkes%EgWU055^Q*n_Fj>g|N|2F71zS^($Ht6?Y|I1!`8}7eP$v;u@q4-;Z`u#zF&tz|QexX$OceVeM z*M9Hk_eITL{cMo@`x58(9)7RM{ndl zS?v3}+FwUjzjyHa3g)j44oLs7`*+y>+e5z*4u7?`8Fu?M!*9xe(hvXX+rJYK-*kV? dTEYL({fmyMC