From 47999aeaa8399dd844f99362d887aef545b94890 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Tue, 15 Apr 2025 22:15:30 -0400 Subject: [PATCH 01/17] update initial values and bounds for GS functions --- ogcore/txfunc.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/ogcore/txfunc.py b/ogcore/txfunc.py index 88a83afa2..078f626c3 100644 --- a/ogcore/txfunc.py +++ b/ogcore/txfunc.py @@ -596,6 +596,8 @@ def txfunc_est( output_dir (str): output directory for saving plot files graph (bool): whether to plot the estimated functions compared to the data + params_init (Numpy array): initial values for the parameters + global_opt (bool): whether to use global optimization method Returns: (tuple): tax function estimation output: @@ -802,9 +804,9 @@ def txfunc_est( # Need to use a different functional form than for DEP function. # ''' if params_init is None: - phi0_init = 1.0 - phi1_init = 1.0 - phi2_init = 1.0 + phi0_init = 0.3 + phi1_init = 0.3 + phi2_init = 0.01 params_init = np.array([phi0_init, phi1_init, phi2_init]) tx_objs = ( np.array([None]), @@ -815,8 +817,7 @@ def txfunc_est( tax_func_type, rate_type, ) - # bnds = ((1e-12, None), (1e-12, None), (1e-12, None)) - bnds = ((1e-12, 9999), (1e-12, 9999), (1e-12, 9999)) + bnds = ((1e-12, 1.0), (1e-12, 1.0), (1e-12, 1.0)) if global_opt: params_til = opt.differential_evolution( wsumsq, bounds=bnds, args=(tx_objs), seed=1 From 6fc372c130b983c134c3fdf87e993a84b160b1c0 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Tue, 15 Apr 2025 22:48:01 -0400 Subject: [PATCH 02/17] add sampling weights to HSV estimation --- ogcore/txfunc.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/ogcore/txfunc.py b/ogcore/txfunc.py index 078f626c3..cf0ef1a91 100644 --- a/ogcore/txfunc.py +++ b/ogcore/txfunc.py @@ -840,15 +840,16 @@ def txfunc_est( # set initial values to parameter estimates params_init = np.array([phi0til, phi1til, phi2til]) elif tax_func_type == "HSV": - # ''' - # Estimate Heathcote, Storesletten, Violante (2017) parameters via - # OLS - # ''' + """ + Estimate Heathcote, Storesletten, Violante (2017) parameters via + OLS + """ constant = np.ones_like(income) ln_income = np.log(income) X_mat = np.column_stack((constant, ln_income)) Y_vec = np.log(1 - txrates) - param_est = np.linalg.inv(X_mat.T @ X_mat) @ X_mat.T @ Y_vec + W = np.diag(wgts) + param_est = np.linalg.inv(X_mat.T @ W @ X_mat) @ X_mat.T @ W @ Y_vec params = np.zeros(numparams) if rate_type == "etr": ln_lambda_s_hat, minus_tau_s_hat = param_est @@ -859,16 +860,15 @@ def txfunc_est( params[:2] = np.array([lambda_s_hat, -minus_tau_s_hat]) # Calculate the WSSE Y_hat = X_mat @ params - # wsse = ((Y_vec - Y_hat) ** 2 * wgts).sum() - wsse = ((Y_vec - Y_hat) ** 2).sum() + wsse = ((Y_vec - Y_hat) ** 2 * wgts).sum() obs = df.shape[0] params_to_plot = params elif tax_func_type == "linear": - # ''' - # For linear rates, just take the mean ETR or MTR by age-year. - # Can use DEP form and set all parameters except for the shift - # parameter to zero. - # ''' + """ + For linear rates, just take the mean ETR or MTR by age-year. + Can use DEP form and set all parameters except for the shift + parameter to zero. + """ params = np.zeros(numparams) wsse = 0.0 obs = df.shape[0] From 9b4fdbc81f0573bacad5e1cdb3bc19de1926f560 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Sat, 19 Apr 2025 23:19:10 -0400 Subject: [PATCH 03/17] test that recover params with estimation --- tests/test_txfunc.py | 102 ++++++++++++++++++++++++++++++------------- 1 file changed, 72 insertions(+), 30 deletions(-) diff --git a/tests/test_txfunc.py b/tests/test_txfunc.py index 4aabc045c..1dfff96c0 100644 --- a/tests/test_txfunc.py +++ b/tests/test_txfunc.py @@ -228,47 +228,89 @@ def test_replace_outliers(): # affects results from scipy.opt that is called in this test - so it'll # pass if run on Mac with MKL, but not necessarily on other platforms @pytest.mark.parametrize( - "rate_type,tax_func_type,numparams,expected_tuple", + "rate_type,tax_func_type,true_params", [ - ("etr", "DEP", 12, expected_tuple_DEP), - ("etr", "DEP_totalinc", 6, expected_tuple_DEP_totalinc), - ("etr", "GS", 3, expected_tuple_GS), + # ("etr", "DEP", DEP_params), + ("etr", "DEP_totalinc", [6.28E-12, 4.36E-05, 0.35, -0.14, 0.15, -0.15]), + ("etr", "GS", [0.35, 0.25, 0.03]), + ("etr", "linear", [0.25]), + ("mtrx", "linear", [0.4]), + ("mtry", "linear", [0.1]), + ("etr", "HSV", [0.5, 0.1]), + ("mtrx", "HSV", [0.5, 0.1]), + ("mtry", "HSV", [0.4, 0.15]), ], - ids=["DEP", "DEP_totalinc", "GS"], + # ids=["DEP", "DEP_totalinc", "GS"], + ids=["DEP_totalinc", "GS", "linear, etr", + "linear, mtrx", + "linear, mtry", + "HSV, etr", + "HSV, mtrx", + "HSV, mtry" + ], ) def test_txfunc_est( - rate_type, tax_func_type, numparams, expected_tuple, tmpdir + rate_type, tax_func_type, true_params, tmpdir ): """ - Test txfunc.txfunc_est() function. The test is that given - inputs from previous run, the outputs are unchanged. + Test txfunc.txfunc_est() function. The test the estimator can + recover (close to) the true parameters. """ - micro_data = utils.safe_read_pickle( - os.path.join(CUR_PATH, "test_io_data", "micro_data_dict_for_tests.pkl") + # Generate data based on true parameters + N = 20_000 + weights = np.ones(N) + x = np.random.uniform(0, 500_000, size=N) + y = np.random.uniform(0, 500_000, size=N) + eps1 = np.random.normal(scale=0.0001, size=N) + eps2 = np.random.normal(scale=0.0001, size=N) + eps3 = np.random.normal(scale=0.0001, size=N) + micro_data = pd.DataFrame( + { + "total_capinc": y, + "total_labinc": x, + "weight": weights, + "total_tax": (( + txfunc.get_tax_rates( + true_params, + x, y, weights, tax_func_type, + rate_type="etr", + for_estimation=False) + eps1) * (x + y)), + "etr": ( + txfunc.get_tax_rates( + true_params, + x, y, weights, tax_func_type, + rate_type="etr", + for_estimation=False) + eps1), + "mtr_labinc": ( + txfunc.get_tax_rates( + true_params, + x, y, weights, tax_func_type, + rate_type="mtr", + for_estimation=False) + eps2), + "mtr_capinc": ( + txfunc.get_tax_rates( + true_params, + x, y, weights, tax_func_type, + rate_type="mtr", + for_estimation=False) + eps3), + } ) - s = 80 - t = 2030 - df = txfunc.tax_data_sample(micro_data[str(t)]) + micro_data["age"] = 44 + micro_data["year"] = 2025 output_dir = tmpdir - # Put old df variables into new df var names - df.rename( - columns={ - "MTR labor income": "mtr_labinc", - "MTR capital income": "mtr_capinc", - "Total labor income": "total_labinc", - "Total capital income": "total_capinc", - "ETR": "etr", - "expanded_income": "market_income", - "Weights": "weight", - }, - inplace=True, - ) - test_tuple = txfunc.txfunc_est( - df, s, t, rate_type, tax_func_type, numparams, output_dir, True + param_est, _, obs, _ = txfunc.txfunc_est( + micro_data, 44, 2025, rate_type, tax_func_type, + len(true_params), output_dir, True ) - for i, v in enumerate(expected_tuple): - assert np.allclose(test_tuple[i], v) + assert obs == micro_data.shape[0] + print("Estimated parameters:", param_est) + if "DEP" in tax_func_type: + assert np.allclose(param_est, true_params, atol=0.1) + else: + assert np.allclose(param_est, true_params, rtol=0.01) + # TODO: maybe the test is that the true parameters are in the + # 95% confidence interval of the estimated parameters @pytest.mark.parametrize( From d3b40e97a4689a3cb87911312f8a0dda8a37ef1a Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 10:08:03 -0400 Subject: [PATCH 04/17] remove comments --- ogcore/txfunc.py | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/ogcore/txfunc.py b/ogcore/txfunc.py index cf0ef1a91..1d98edbb5 100644 --- a/ogcore/txfunc.py +++ b/ogcore/txfunc.py @@ -288,9 +288,6 @@ def get_tax_rates( elif ( income.ndim == 2 ): # I think only calls here are for loops over S and J - # for s in range(income.shape[0]): - # for j in range(income.shape[1]): - # txrates[s, j] = params[s][j][0](income[s, j]) txrates = [ [ params[s][j][0](income[s, j]) @@ -299,12 +296,6 @@ def get_tax_rates( for s in range(income.shape[0]) ] else: # to catch 3D arrays, looping over T, S, J - # for t in range(income.shape[0]): - # for s in range(income.shape[1]): - # for j in range(income.shape[2]): - # txrates[t, s, j] = params[t][s][j][0]( - # income[t, s, j] - # ) txrates = [ [ [ @@ -642,10 +633,14 @@ def txfunc_est( # ''' # if Atil_init not exist, set to 1.0 if params_init is None: - Atil_init = 1.0 - Btil_init = 1.0 - Ctil_init = 1.0 - Dtil_init = 1.0 + # Atil_init = 1.0 + # Btil_init = 1.0 + # Ctil_init = 1.0 + # Dtil_init = 1.0 + Atil_init = 0.001 + Btil_init = 0.001 + Ctil_init = 0.001 + Dtil_init = 0.001 max_x_init = np.minimum( txrates[(df["total_capinc"] < y_20pctl)].max(), MAX_ETR + 0.05 ) @@ -680,15 +675,6 @@ def txfunc_est( ) lb_max_x = np.maximum(min_x, 0.0) + 1e-4 lb_max_y = np.maximum(min_y, 0.0) + 1e-4 - # bnds = ( - # (1e-12, None), - # (1e-12, None), - # (1e-12, None), - # (1e-12, None), - # (lb_max_x, MAX_ETR + 0.15), - # (lb_max_y, MAX_ETR + 0.15), - # (0, 1), - # ) bnds = ( (1e-12, 9999), (1e-12, 9999), From 33f0ae3a40be9584babb94658ec1aaed5fbda5f4 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 11:07:52 -0400 Subject: [PATCH 05/17] test with HSV funcs --- .../tax_func_estimate_outputs.pkl | Bin 5596 -> 5424 bytes tests/test_txfunc.py | 17 +++++++++-------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/tests/test_io_data/tax_func_estimate_outputs.pkl b/tests/test_io_data/tax_func_estimate_outputs.pkl index eb9a071a61135f1f4a44ac3f255cc9d6dcc04456..9bf2c6d2089dc98d48aef2758d0d843395a3bf5f 100644 GIT binary patch delta 467 zcmcbky+Mnmfn}=dM3zEErio=Hq5_9a6xQmj`fZ;v$wl@`-3I$9B~y|n3o_1`yqeLF zk!i9si-xE`@{X%dT=ssk|K}&YV#@uM5H$fz(G@;f0pFn1=CSx$LsTLS|_4{u3gMSNOmUUGa%WkKqc9%hf=uqi!U#U+VFCGnN1 UiA7U<-?8^FIfYO0W-QeM0J(Xd+5i9m delta 708 zcmdm>bw`_}fo1B&i7bVTJQK@IQWIWxs`uqrT3cVu2uzpUa@|pe{bL#v`{WbROAA=YF zYEFa8Pbry_H1Xh^$<@Sm3M+ASS2iavE zpZ;a5#v=QR^Q(jkCCW-MPq~iSJnC^~|02vu1+8 zOMAC-(!Q16*X$owsYd?zRRyQ-{N?j~9p`6%H2>GLBlZn&agfCsY8bE{Q;4o2JjgDZ zT+Ce0XgK*Zvl)=oVgZqfEOtP00gEN0;pRl4dyE`GJNsJ Date: Mon, 21 Apr 2025 11:17:13 -0400 Subject: [PATCH 06/17] test with HSV funcs --- tests/test_io_data/tax_func_loop_outputs.pkl | Bin 29763 -> 14163 bytes tests/test_txfunc.py | 7 ++----- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/test_io_data/tax_func_loop_outputs.pkl b/tests/test_io_data/tax_func_loop_outputs.pkl index be5e711de0ac53f930d11fd86ffeccf45ff71fa0..c21b7d14bb6153dc7d237a701262a7bf650d95ca 100644 GIT binary patch literal 14163 zcmeHNc{J5s+csrZqNGeoWh_HTDR(2Ms5>GO2XyVkST_dV}^*ZaL|ecvC?`Qvx3z0dW#&UN;=_r3SNeXgor9ASI1dt>y`24%&L<@-_)VEdVnDF3>mNOA zZdyt4aj{7jNpi;Q39MU^%*@P;=H_3&FUBkZNiMqFN0zKN=7C!%DQ2EWVoBh4@^>Zp zcpnV#KST(aS0S+{aCs9MA9p^{k99lt%2jfF+OM{A!w#Shjv za_!M~34!{;ev0~;NCnHOU53q_KX9?8t>&7jpjcudexQ=?KpU~M)Zp6%2@aqDLW*!FAlWU zo=RQb90%LeY%;i(pMbSTemNd0jfd4!mPW_?;$hFDnbH2Bcqp@HWi@>n2US;wM#-z< z;E|iivjq#{prNv7L^doIwAA7z6J=vSFWzCJYu#}mUl-ya8Xp6%wwCaV{85ngQ$weS z^C%4ZR#eaZI0EAt$k}mi7#M7=vujrgh2A$6@uj(k!8)-hqM`jD^qU?$-?7&hZXIjM zj=JOt1|rUD4+*+Mq}m2?2{j^^9ZEd?TF3=F?j8&px$gizTLqd971}^Xoy#?)7E5qC zyvX2KgE=&OWR-W&*ac!42i}MK8iQW1aRa%02UtJ2AY;Yp0*p}n>+9Z@^zPB3-B{jV zEBx=?NtzF55*vwIm%F%EJC4jcxs4Ihzq`y=Q+%$ktA;;y1-NxA;7w@wbJx{ul)UBM z%RhJh!_)5@_a(`Z7LgYJF$qY*|Cj=7zo&q%;1hmwn}8qFo!$j56r#-Qle5~AmkOblCXf!9y+ z<_RcyyS;$xP(Fq|vJTyzT}gwT6E*q)>bV$pn&PGbM@Je+Jil(`9D`$XFF)yB*-e8& z|G`oVavnyzTvK=fRb?D9+1&c0nzJx$&m(dHiOa@e=djLEkq8`HrTyHo)&(?(PYS-( zB!K%ZolZ#SoYI5`=||wS;{+VezObAG3!&FRcKtZ*eXln^PH>n6 zQFYR(B2L`5bdb>Q_GP1k=-xYDTNQCmsJWtNyXZ97+Y}!6s{}XRLfVp*i`UFM5k|$* z-*CRef^(vuKAM0`GnqU4b#ZL(T}6FW#S^e^$D~rd7LJ{3CbCntl?GXgq%)WDalS(W za}d#%2BmrVPS4wMF~}z5sm9Yd4w{BGIR}5?a#h9~J^efMao`QR^uF8~r)_-k!lCIS z!v($+T~XzJ9fuY;o-5#F>`0DhDN)&4i(HUkrb-m-NLXJ$2wkL{B{gpozUW` z9>0lU7ai_wo_RopiDGKk!!DeYJp11J;eu3<+v(t37KC#mAwrofn;ZrELgh)0z8e@P zlJ=vK&Tqzmz5ASLax~8O_pf!MuS>_^XM>cfA{&mqRGaIu-CO1`g8uTRh zKd4xNb8;&BvX__W6fBa`aJ{<*mwh~Tq!-vXPXSc@l6-}58F{26FL8h14EXP-)HXTd zeA`)k=QTH+gG%$8L#r<0#-y~VV#+CW3gkw_ww!*1^PRtnbWPE70<0U0J5>+hGQ@{R zYow-5PC?c+kqv{@IPG%1s!ksEIXGMNAvxU%=Ub$?cS+-RIwZW*H#X|QyG~869&>^zs;7F7@&bLfOaxZm;4zxmp#cRTF>^;6tvJZKtfP87W z>dWUi_9w$2p1G_nEdQB{>wn+*|LO2M z7EG6vvNuDT`42ZnNZO+s-?tyn)!CzCoZ^c{cG@EwKH3IB0}qrzZ*P{_+CR=cGulTc~45^wVa0d4=XWTsMsh>}E~_~@Nx81*9!3D)NZn=@e?Tl9DMwjJ0az;-QWIV}Zj>stQAmWkUgParh zhZpAWL9TPjSxITOsJL!6=%cC?YPZy|I`+&0C2(aDmvcBHyGU^&7pEg?sW3OZX>Er@ z0*?zSFS0|G=bJlko9sqi=4tkkiYCZfzGuhrqXsB}b%FiN8hx~Ff2J~vnlAFYm0RE2 zyA5qzU?X?3S`$TAUVhv`Q$)m|`(4y8nOM|34o#{}iCVKOpbiw%<5>nSuK=80dJDfgw4}b{+%mnYg-u z+0JKRL?HvqnJ8b(Xum3Apw?Xm8kaCI`yK=DGuwI0_S^f6{tRaSS!TQZA*0`&iI2({ zZ5Jk1GqHz>1SanT6BWxDeCGIqqB18Evm*2ButOpxkK&k`fu{8_z%%CaN4^ zwnG^BBbC`^!MqawMIJ`x-+23v0R5-s6V9W(btWH1A*L`#?#AUoz z;WpH_@_$@^_l5$CD@u>w5X1S-{JhOOriB7i>UNeK+i`5ajlhL5?-idku4A^ErGr<;o}C?nD<0koNjz6D?83Jl zhJ>5LaMH|IOGNP+hMkivW0fa73_{L{m!90Yiea-GguV^W7zT-3Mb*3IaN2tJc|`jJ zhQV6;`@7HXIBgl>aM|HiL!irj>DycObc~ZRq|Kq?I|yF|!5uFi8)#2atBNrDZFBVy+o+ph!d)vFz>z{FKWx*8zt>AvRo}t+ws_mt^toaL2mE&3DLUGyy0unB6qN5OfnEmF=A>5n~D15(ff$137 z7V@(^|Ax!4uh@2X$BT`Drrk!#?QS@BZNPeK1(`AE58jkm5`trk)Hkwk6rsXgT?Xfw z_1Ln{j;Z(k0aW1RJLU6c634!7*l6G~Un=0rE8Ds(6IYrsim0Y-(4fJ~B{jBv>%8sA;g}h!m_(MZECtW?KA1M%3*m_7=01Po@R?<+b@%`HyIrV zu^;BIgzIrm&^8}?@yk>YsPf<@+u*dKc-xVe=!1epOIFe{1e zI)2aCj$trv)@+s9iPKJN6)CkWq<}47V{zCnobQ^a7uZIgP(Z)v!cgioF0*Y_krVtZ zKLVnAn)6kTp zy6*OF+_&V4<&9WzoC4w7a(0?E;GE2zH~#z+D8SFv>Pj`ku_rPTv>#86z~y?bsI+=q zR1N3XUpn<*eoaf#yrOG@WA~LUoNzig0;Tq@W&O2`70^?MeZGA1AG!SV^2uEh+v3T4 z#;9qNkM`zE7ASDK%Qr#73GFI~Tci*~LRn<{7p$uZD82BQ^U+FYwEA*&N2!Gi8cdmr z3<-2W84-GfYrM{gKFwz+G`$Zs=q4@D`)P;zguZJ(_H#gylNP!A@9#nCvU*;h<@O>g zS#SS#svUA(L}}8Ewnk5byi0E{vqmA?WFFK9TcN{j5>$n18hUo-f7lfnougz7XuF~)h2jrw#R3e{m$VxFP+);|F97-B5p^5yevSX>R0mk&uF5i)=w(C*Qg-ngQT19kP?#2 zjK9%(d_7vtF2o%sEr+;loy!jWT#bkqw7%s&kwG;^wzj71@<^R_!rc0nHnuF%duC9_ zuS6Eh`!^1r|I#ggmrwo}roYQ4m|*#TTt1<#pAbAuqQl-RPiE&rahotozIDCo_X!Ae zIcc|}6Su80Re19*;C7|0 zcBVegiNqDtoqnz~IJK$3Xlw{~lrpRH|0q^9U&yiGuye;gTb3l8Yach9ck(XZE<6^u z!Mpl(SDVQEGx5$^^j)+W=j7`H`okS2S(>$Q5<(LUEhk_L}oG>H*_m+57eRXGfg&z-*)Kfm1Y4Nxh#Xr;pS2me!Zw z9!!JYE}46?VmS7Y#ljC~jA-DzG^jTFJ#GZaqdfMu<>QceDdK2*KhAdocba8N+&C!A z#NFC(7{?a;sCY~~g9h%DvPE|e<9yeKuBoghPJoJO#H{#hT%5$ayGDjIPJqt~^G`nJ zxcL6GC7o63#{>{Vd|bnWaO^s}4O0h=Cn0`}ecvk<+)+xA(mj3FdlIHnU$GMQ;nY#)&j$L+P`H67f zDe(3b5tx<19nb?NEv)Tr>ChWBk$A=rr(GFoe!xz03XX>M)1U~KHIME0O_1cBhUxdu z^SXcGoRDh9&#T;Qjk~Vl%E-~XeeGBTr=eGR_e5_yPMZ?^%xr*f7IeJ zCUa2KMy8EN;GCRrW;;_wo(7WFJ*vA2PWys=4`IitX;`_$hRwwf%PxxFXkR`H*MF>N zxTuVCa&W_tR95>mY#_cG%df_1=W5S6a;DH>?SZL{Vc&7HAvI?>Jl1&yQk>S?O}KBFz-hsnvX~A_>@Ey+MdP%k%kQ=JDbV5B*2^KR9ys>m zXPIqJcGE#zY4Z<}N*tTDI_dNIAUXuDXth=n!^Plm?eQ<9t8{qB9rPRp;MkpS1G=_W z(?O9YqH2pRYc{+meiIs?LrhP}6|c37{O$RUW4?UyAG!SV^2sl!3U}9P3&eX}rtdh> z4tcQ`oyZg?pyZDu>!)sdp$7((iJ^RVq%$05QV>l1G`e|-5fP=_ zUnjKQ%>`*5j;u{#cR>?^epcsgozdZ?EXocG2`Ex@q2uXLR}?|4ypXBrhBWQItq>D) zN8jG8HTS2GkgP%?=Oro$Ijz}@7G-;)H^y`=nXjH`xqwEla2*N#TC&cqqs_QQ{7^rKiaKwZ4=TR4 z?E+2H3z1$gR3jboMANG!jq`%}%U*mB>K{n){p?rJM$F>43oFDhRjpup%9YHp=6#VqL3*hQ)Y<_ z$xL&U=QNz-bKiIO@p(PZ^Bey7o$>*Yc5dW11*}PazPJK5U$VJq z;lOKPZgr0LqK%pLc|9vDJ%_M>)t3$Q%=E0nZiEF;{QL+`8d^J8ocp0!Wh4;h@8=t) z6vib(7CjT@R7VA;0FTBxw zeru}H57;-oFXk}y1Ma7nn8KxTJa2a4ejz`o<7A9Gu)`NR+;-2fYWsn}uAnHUPCqyv zqZUuebQN~FOzPMa`2)v-vb3kAKWLVGdpBgkC^2HFprG)B z<|iX_yPbVO;Go}^z}-F|?tkX!xo2K*F_DhkP{tF?db?c1X*?ih;t0xJ?GDoxWi?A< zF7Pp4R8MD@Gf0X&)9B-Kg3;$?{zd6_pc`tYk1|NZ@RSM2b+T*VbX_wV2Ee`X28 zm9xnz6Dl#Pyip;)i=Ze!}(EOjO+o__rgbQl_Q-wdL zpZ`!hjLT%b$;Kapz=R$@1ZaK^0kN~)bXO%FFe!=9Pj+4yzTc0k#uGCX`(B~E);7;g zJ>w7+`Q!B#?#NxeZ^BRo+gf`#!z*7w;UeY5w$ExXt9vownTZ0lugm%HPRk8MMwq2;ZDlYQwkV0gV z!Z`XGVV%6^VI^8+{q0uu;4pzM;Qe%ZT&fT|_Xf;&g4D9v1QHb3y{uJgtgfdk|-_itqfsLKN~$&Y9!N0H|&_Ib(39oNyMI zsq%H4O?@zG6RMw6uK@C$vpd7x4#D1*7M9&LD!}!k+Vt^RH-erxeM!TRlay9f$Sauk zWV}N+63)1D;0w-*iP>Yr=%2{| zdhXeG^`sV@$u`G&aT&-FZxgV!7g!e@bFbTvPiq=X2bp*Nw51v z>Z9`-Hm+J4N$hw;plhm*Y^}z63Hiq7zx^52b;W5Vz~>*YC+>UK?&C}yyH$m__aqKy z*|nl1m(IUQ*Pm2vyR`8}=EsA-&ZOr``6&Ba4?=ss@*Pq2D^K zy=tP2K+m>(R_OSq7rMyJM9-Z)4IcfjH?(IEBviApcJEe!?e1!lFLd1qdgAnxg}Lwj zt<;3%F9zMv&}l_tiaXftOCF;cbF0rw0e!IDtf=^_R0&)&UmDry@)j1!JePdSdw@mh zuy22;4(J&l)^_`N6c`QoM7DD&LP!2k?hZLOpnJ*8xv0swLf)C1#}n%$>xFxVLuFLX z)}un6?42a%D^O4nSLvJf9w3SMa?@Ss5pc1zc=cm%Au^wzD6=V^GX_O-l-5E4z^PkPOhM@;OG_RR*z~n?9Ymns+%& zyInCKt>F8{wJZro#7+qHu-OQ}^-APEE_O3wG8s902=jJc>;z?6=16%|LZItT727m1 zbU`eQ;hKna27>Iwn0%)%D1rg)weoIU)#2^Lec#i1CwU)bv4|&A&xcVVBVnq!$A?_( zQqWzycXso0?I6Qc>Lhg_fIxQ>p7szcY=YS>!Pl8;4nSa7!n9GVCeYn>igMpQpQ1!cLGHwAvy!A1lKg>FuRJobTf1Gj%e-9}#iiGg#*K+#ydF zd3o4+Tb%7(v_7zGvT4svbc;Uq&5fgNK<;qI%Tqs$;5)Nc#j6CLW;k`VdZXl-O$3Rx zSiMU8wg)Q2nZrhKbiAE7UFe}}(D8~Ck>L{6{hWslg(u7$zw6XnA?rrceS?zi1gie$ zcojM8=TNlW`r`NG6CUsP!F)Zq=2O^vAD_}0OR zQBU+P~)@}~b zSQm!2C)q46nzg}4L&m)~f-0cVrAk!Z=PmFJ&0+69y@Xk9b_efMN>CY3vrQ=IFqDW~ z8b2bf2DbJ}atbPL1U+&3rS&N)3EvBaZ&(dhs&?0-{f1+iIR~v!XWz~y0h<;m&i!gK zcw-RGODnmCoo)uc59b8u;(Ng(ter1&fC=`VcC07ue~0uKqkEfNZK2Y8#-m`ykf0|{ z5BMM-l(TzGxSlbMZsW&i=$b>4mE?W)q;Nb z=!Rzg9yqyS)!dVKO}HbNy@okk8hrZB=rQAZoY`vf_$o#nf}S{?7jthpQ1L|gQB>aB z$dVS6cVI)#^`Hw8lJ(4l4~J7YJG6lm&e&f@iHWT z)N-_vQRPtX&mGhC}Pm}@;IR3IXlwB}$Zl4kW@Eguq!(iZQK z&Uv)}&3;u|A@}9gMMPFAnVv;IxcUh{YcUfeNDmo9tJfP;KsH`<-$C5!0&gcyzgsSJ zu1c<5c<-?@uMh6%LOV?aYWPhe(eZ~Zv~<^6fScw+T$2=TRdMzxW5eAikVvBCVs@h+ zXoT2hMd$UQm8UL;j#3gfbLJm#V0VMJr?tkuma7x=#Ob$EgXr$@Ul!)LP`k}$z6A-! z?EZ3I-2=@r+I$=DYlfPvttF{ZtzgKI`H?KB6At+GW+znkK=Z4suNcP_zJC#pJ-T;oHQb^<-aQrD3pdV=ubKL`gz&QTVaF>DeK)AuxG73{T8m&OPG=K3`z1q$ zN`%ICOBz~F;IT@eXL_&np{5VSN91?<-;;v8+$zU{GDdhkY_#(XZk*q0Wu{NKr9sdW zryo1KR#94Qi%8g=i@A@fT2PrV)3@uS*HD?*(RsiAMndFy_ti4*Q>cUTi|NMM+q=NT zTe;y*HapM?J~3kF`-~Wd?Hz+~bK8NDkeh7>&l2>+>6}L|x!+J=6?yqVi{rpjAG&(R zQJTZk6A8%JPu8?8(@Sbz*{P^KffrI+++07r1eaq51|KZlA-v77b&}}>lq6_=zI((C zB6PMM-FI4@peIg0ZL4~DXnKMYdnKcDOc&o7JW?dw@I?XhJ67E=6Q%l;dYO+ZuC3Ynl%Y};&htG2;*aZbRzw> zH}fRQ+tBGdC9;YM_9%&Vvj^S!X9Riy^I#2!RV@^}-gaeIYd3h<*$8~jF@;O!rzq_d zlyO(j4iyXFcF^H8$x)|Tv9J*76Q{H85mGXy+AQ*2<-jq!kZwdFcGOm$ z#UHIszA&Ivzf2#froH)9tO``sk_G6u_rYMxjBl2a5qJe1nT|ec2n!t9@t>64faQrn z9h1K%K~J16Hg)W=a@1B4!-=H;rxU&Cz2yB5hP7d+OFU<9XY^Bo?;G7~9XSmu;ro%( z2HZ+7VdL>@0|zVJVC6FdYllAn1wH)p9tNG4!;) zs76mk6flMF5MKnZ$}25n(>p%R{$yS(Vc>36qwja_u?TLFiQS zCSSi{Q06ilcrDNkRm`td4qn_!kPCH_W&B1ekaH`I{rCw-f}S{i*M6ofedS|9eg}$j z)53-j&-0Bcv|d+G?z^DpstlE&V>6;kZdXUh06R;oij?^XDB36$%9ULOf75i8O;zrY z!`B!SWUT{tDdQfs9(E(>iPOU`MpkU11R>t+=FqBNuW@Twfkfrw9Uvmwdf%*0-ny)f$kj%)>%OPsA!YbQ=@2lo_~(t{VQ<6ZhSaHb?5!GeIGR@RNC?p}pu$$DUig+{TCt zec2TFsu(=41WMT-THed-N&RHMOK$+m)W+V_IZ=Z0`II9QYG>h{RK-pttqxIleJi@u z-3WH#^h)>I>#aUQLfo(3Nt$@}p(jCaN-k8GAdpDisM4JS3?HqxEZ%4(MBb#0e2#%` z92lv$Dvrlr2ctuk7Z*(pVE?zBX-_?kAZ{kya6r|KpeIiMreX7-Si4H3C!RTu&!-QW z_puGfy%0lOC-+My)n>xsfN#4}LZq#JCkis7^`n}FbT4R<8eQeY!HN!*;F z8+?^He9dgy3M{AY$A7BHL3llJdO_o-uySQlp=YP)=c(&@(Ato-a-&?TXenAVDY+~Q z9$esXaX;2V@SXqMwlDWry@hI~m~U^p{b5K6tL%Jg3#23T7cY~Y1J?!8x}6ek1U+$j zaF%zCg*czk=}jbdJ0EqR+Gg*V2@!o%eR^M>Ic)~;>l}$cz_T3rD>ieihL&THF&o_U z#n=ygdh_*TuUr6Ut=A{HrHx>GosE5gf*V0koZhXr@6f9)l_I=z${KmIP3YaJD$^0N z%V_9ZRsUVHOlVlBxkx1!o1J~W*$X0e^$S=$J`W$)1`1Y&YCul7 z_gouprp4=t)3+Q5WW2`+Lg!g`e%)&F6760Av&(Co`AH~vvvmZvSm^q$ZnTE0 z{@=DT=X-*<{XscWJ3Xkby6yC8zZ=0$+;^rs$1C3}){5j*l6)vuy^Mt84m|zbb{5&& zNbgC1m_hJ)F1&I73}+=62W_hzvg?GsOy_hA9W0^ZsXBAsReQL-^kQ%Z*W+6D%usRr zsuS$Q>9-EEjH;gH6H*d9HW*a)0tKX8wVq8kKtb)_X_`Dzp`Ockti-O3a26xe!81I@ zZ(;gof}H#3FbGUzB+4o-bqE*W0>0LAQHY7HQ zttJiZSYID|I}kUqJN?7ckW?kwxqtdD?B5^wn&X2pTz2Tz+LC1gmssLrJ}bGw zx-1?(5!`Q6#OWa%r^@z5V?x|p-?23my+rbA7Y_>s>!D}6JkEJ3rouYwCk7{uEYmOd zt{%N9It@0Va`l`}vEY;Yur{RA6UuuQJqzkgK-WudhX=0ePN*}l618w$f!@5D^^(n% z148m$MPIDWwxgCs&K$91OQd5=al-M$eOM@C8uS!h9&snA%$puEOhD_q7mJ@SM1Xy! zKhvH#Z;&YNJ1qOf7*rWtbU*HOBj}0y&S2o8SLdW5lJ+iUW6ACwl&^;sAI-H!rMFWK z3$D%tRUePgf(HYz+ijPH9WJsX92Fl=eVu~d)GKlIuTp_Fu2G}B<Ea=tgZWK; zo9!&Dxv@XjL+37Bf8$82UK|N|R8?{`JM;*%4DDpC+q(##6(Y~fc-5Dhl26E#HNkCu zcRz{=SYIh|=RCTrwb5KTBpu#uoTJE#U58>h(7YMnh2JqK}VD9 zMCgfSx^-^II_d5SK!e*2#dn0k0N0vxJY)LcJEZ1LrD*~vQkvzknj1k+oGx&WgE~Wr zQ)tm#)O>^VAo3GR(_XV)7fDg_=&?Lk06-TFy5QB*>6 zV09=wIIVmuNbE8^ykB%)_n0y0Ole$7I_XBx6Q{pB^MOQilw8QDJnhJ#=wU=Vpv}6; zff9*aRXf#Tn+#`b%iS1emRA=_&PvigY?JU|X-2$VDFmKBJswo(XA3(NyY-yAjo`5k z)hip^+6J%xK{rh&B_a90^LFl^cM|-^%{1d5YEB)MxsB0m$t8WQo`PK(zSGFEKOM8i zt(EgHvl$dAvWNh`=$qYE=wl9#(inPQt!4536FyJZ^84i>cfk6~|8IVjV1d)EqzO zVs^^R@4O~*vHnj|-L1_z7|XTPb~4Xw%xz(z!Ynuo^J<|^P*}*oE*>*^G^3M+ojNJG zb<4UO?4I1%6OvCESZ}cx{X|tdR(Se~;X9XfteDQmi_|(5ldAb3?G}@QwQgpm6_iQG zl(go?UHy}>_Hv$?rrq~2Hfk|;(`R=v-Y}oNg%9px65^tYx~zAxF5XoX&tmUjb5uK( zz0wmfTE|tkArn`cp(Ja%rMinSMeG~ z;iz_hV%iT|G}UtY^27(D@!ehI$>NE*9+osZ>*$Uh6nwd>mDv{4*Hz6v)Zm1fN^D|( zFLxPh=WTJFD!Pn4seEM4)n$e~Ia@UGWziBNQ7Y!We%c$u-&XQ((vWw2DE6b{il6(p zfB%1~|MO7$j{^M($O!#)g0s=TEnOvLw+BZ7XQ*4|CS{We;w=pmNmcS z$Zu)z`*ZETC3oqs&kL3OlI_1gr~F%f`7Li1{jyX4{&%{+kMI7LDZkHa%kT4b{e3L_ z`&jO`zpH=C?B7!AcU(@T|N4IZ)L%03&Mz7DJMO}yxTn*R;+j{+rnJJzV%2QB{#g=mBG;_=ynGDz;} zPyS*|28TFax;@gKhqDqE3IXCxP%V68v|Xx+umtKFq_dU~Z32<*BWtcpG=cZTuMIX{ z;_ZJPA^4lYoT%z*Ew<+%jUxg!x{n?~Ms_DJIYt8+mfmZrb>D&1r^bpU<|f#kUGW!s z`3;MaqAg8uTgLE@biAE79Z^N+ul<~YNag3JQ?HkRQsJjgTBdhUvR))g%Xye7WE5ycQeD?F z{*&(0Et6P59d9R2r=p~*{=(KrlhuQ$HD#AvaLyRWM^AZN zrurwn_Gs|EeH8yr$J>e1(ZLNTw;W7FM!S6D4R;j6%n+wjF@?$a-ikmfSugL8L>IupM| z^_riq&@{nSw%NZ_ov`sD^D*isSoJvck5z}a6Q_%6R%(CS6@iLAT$i@Y%Y!|AQimt< zE79@XX8*{V*ZWbb8>{=4sZJ8`-+`|-GMH-ph;y8VwH$Yw(O0@qR2 z=Z#3U^1Sv=ihf8D+cY44qKP2!bIp!{LbC5-O@Np9h~p(=^i9v5v@4Cc#7*s7b`$P-i}=SF0YO>0o4HuqWEb-h5<;OaU3hi>C$%$A%?q}$-_#CxsUc zo7mZNDjmpqqz+Bh7NgY?A_)yAd*FBoPkj=uqyBxoKjdki(!q6;e;e<3J8?SKncRh@ zn?6WXbC+{od>SmiO`n+`3P)+|~k9^+HtY0be zx{YqbS8sTsExA28bc%1?j;}4AtFBYK?@}i;M&>k>s zZ*lO#)h@4Lcs4cL^4g})(mGaTZ4=xhIo9@M+dtPftx*Zb)5%t%--)tm8k+J(rAFn& z&*V~JeSpwAZtieo=w5rSXs{CmZ_K!9;PQv>xcFX|w?FjmM%?W4Z@uoysU-I>%9ZG3 zot>Z6uX!S2=^SC3gcMlZb0f_TQqZR<<8W@PPGC4Te8hsO2|gI;zdVnYeGh8)(cFWp zhl?Z>q&sn$!RH$<@pj^ouS&ABEev%>i?i!l4fdwMA)y7MwIy_39A;cG@3n0YliZ8lkj0d< zUrBW`*v(cy6EwJwK5Z>)3DWBUpWV_&3~(douQTi{e6O6;V8Qiy{7j1X887j6;*tLz zv){Rb+!Yc^{gqbr~RFNe7xdkAOAPHc3a~Y z^%X~4?QBjmtGWx?B1ZeQO*|RaFi<_xdG{E_GRv)xUFsxI*Ts)byvL0&e77M`g|Ds8 z!p#GCp9%7}p1AL=*LRJv$hsp{_PN-BOUZEFoA=6?`aL9g_7=b1#!eWDZ#r*>tL?x0 zepd}H@3g|&pqrhr|RtyxU)LP zSMkU7Jl;;6e!0D%_ez5^V$zF3x0sXRrUDjCRh^7(A9tEAHhKYh(TeZKaMkm7-}ef< z-;%iM-+jm1iPN12#TS<}T~WYI4F&zz$xs-hwcW8H3w>yc?!H6b2H4R9$@RFN_IG+w zTU)*7x__q!v?{(@TCvwz(WP+pXIP+BpMxjljwORm)BOkZq;W{yM!ItIfi~EXQ^sqD z`xXUe$H2vd%T?!hN?(Z&m-nZIEV=jh|Fi1wcH+L@5NjK*Vs=K)90y4k7Ls9)lw7>O zej@7WYvxbd)(U%}`cpk_miha&i^N{ztzNjX@Nd^H@OI+#+0~mk*+}frQq8Lho#+%e zXA~Ez8E_d@EvnCd7JUxvL+#w=xE?n_O1GgEcdY_H;-aqD%X8s+oIzZ{p=-EV=$WO&PP(Sll`7_H~ijQ{+i8TL+CYkm5|9GB?3#Wjb>9EZ0P zr<-`NMKRB~p~oidtJ8;4VC0$}Iko?kNMC_Z-(n~3CU2f+Oo})k|K4fo>Mx0`S@Z9m z)?s(*>7y(5dHylG8P&T^$l$=%2QwHSOjutNH>W)zgJ+)M*apPN2<$+Q}-kp{mNCyd96p|Eafp}wO1J{~Y zLxu>Cmen8f&e7B4!c8Racsp_5D^(Js*RJzN1|?br_pfKb$4QRXhM915S~J8sK(ZQC z63aYja5e8|AjMNswmjmVA^m&RxIREaGNE8y`p*$Zb@4^*J6vXn(~I54r%&v@f*!~9 za2kZC1A{?iC7WOXvYIyktd~~-&$xTnefz^4C;rssKp~Mi4sRz;kEOk8Xufm>1$sr- zsjta`dkJmo`?-S92EI33_B0RR@rewZ*SK|^W53uzGu%(s_-?~LEEgG#+vChTj;eBY z{L^<0#WxW@?8NE8**Ap~>by|qm2{K&3t6Dmn$9oHWrEZSyFOghDS*8)-a}UVnjp{K z=`S;E;+x%+`}Q=!-p$+B{^=@lUL%*bG_D5{r=vpWIpFX@ytef1J3X_3b+u|#bBzNU zr*`_dp6Wi34ZCaHz+EN&Xp{Drns;*Ja6X0WyqD-?|5)>56|oH^xVj)tPg>1v(?4c~ zOfP#LiTsuc+b(V1(j(x5!boBa@A2h>U!TXu;fqb6-5vKIbX@!$IFH1dHbJ|M=pX5L zJ#l(r@Cj!MT@RGh%alMxk_%07qg<=^-a>iSA(s1h6oB)R!G(Reo<_UtqIRrsGa>Kq zqgP)a!L50utib`E;i?vOq2*vGy6a0lE%q>< zu&S7h&UvwcWSNd%&&j_aYwTHEuDZLkzkVfn>MEvlb;Uib_NV6@uj6dc_aG^&T}!!O z&Sj!+;d338eHgj*=4Kw8y?ATWkI3Whc=?O(UDSdCKO*mxsrpCX@pj_A*S_1Pij9XO zwlkQKljeO$U@%{ry5fbzYwmBezLyV4I_{*aaXaJ#A+rCG_gUDsu#>awXcMIJ%u@fU z=J9so^vtiN9x~LI(J3nJLyR9Y;q^x>sN|zBD!NYbb>1o$ciZg??l|23{qK42*k~uO zj@#q>Ti)?@;`I0mkD4k4H*{!XyT)hTJW#jboy0!Hqv>Fh`|HCB2-KHXyxvgm#^;@c z1W#ppuIgaO`Om70erHtBzZ&;2fgcYc_}4$upO<{1QJ1P7Qi#T2Hc1w1atS71;UmebnmIqPg zIp4(2d84fx9&Hd-%q5I>@dBg2bXxp)FIAvc`%|aI+lkYU2r?g_RkudL)Cw+IWLdDk zvDQA??Ixo7o@iJXok!@fd#X|G`So1wbls~%=ak{i%YMJ~-SLC=kB08%r z4|w&Xzmgw`M6Fqc2JAigU=b1hlgjV%`7iX0QL~>^tA@Hi(#@@n`Jb2Hsj{Zlt+1_~MKjt_R61;EW=ez#c=ka}j zIDPE;-rc)bo1;0e_s7oXWe>A0i^`Wjz9@OpQe}4K7ZeH)-{TI4uUDw>E z_$IJ-tp6i@I%mUwk$3X79cD^(CW!w`@$D}t6YM9LAJ2|^AamMt_30(K5O-MDzpbkY zhNhYRLT6|p%?xO7g45SZ|42V{_`rWb-{Iq)=MtcgoDaEDTz!@f4{7Lp-Mt;q>|^Np z#GMPrwu>inQZ_^N(%Qe!ZOk?w_Mm76k(jJM((NB8{ugxnRa&!(w>4138nFw^T4_Mz z-M3TB)D&5nGmneKn)QxCIGwjHW|0Dg8pw)jtKR2iM zNkK*mt$t0}lJ+JAcCC?*JQT4Hy=VRKQi0^r3+3%QskzG371Vqf=PSg#fu5C+0DZ?2&Xj|0ca6|f4c&CkLs zSACks!cUI9Wfh;r`oFo|D{-8~q$k9Z3sYw?L$2#Drnb#vG2T+<)Z%m4hJjl;41sf4 z5$(HfPMbN5$B6S;Zs#mky`kWITEiU1u%Ev}%xoSz$#ddN8QVMT`j}Dux2k!p$i6KMe5^zGLwwD(~WMNHBY&Ry7Ps+b6UsTWxMJ)>iS_nNRL=Z-x4sM>%fepHjS3VDi+ z6rCr}rm4mb@njf@ohieHFY2)zn0$;;bH1%VpZWkZRMKrqSS-MvItuShcgV-2_Ku$) z>&n9zFHcu5FlS&UJ`D6q0oB;AO~t>zb8NWVH`f09Zk>Ps*#BdeKeyQbxg87p8S($Q z9qV`D{r~uOESaI0u6G~H(HeVlz1WC6pf!)7{#y7Atzros$Q>F*8r$j@Zs(N2oK@k2 zV|2I$nY-)Q!w22CZC?)VsaeG^kiPIeeQiE?+d1ty!n+?fPygC*-8^t5h-@l%g$3AY zeIB0@!y+_({jHaTWiG@RMmu-LQ^2vCrDx&$I6C{zXq~#sGw?8-n$WJv0%|7GJf`v= zbd4t;xJuz{pPl3J;w-TA|5AHgM2SGh>;F7H@i&9fCCo|X4CkTgt*`3ER_6eTt^K0) zpb-+e*3tILELPZis4L<5m6zb)#q({#rySg4=Z30qe;0FtqS}5p4aiY7M$Ah zB=2^lIq;KGKmML1NzfCgs~_ikyJJ^A3Z^ALo=EEsPHC4IcSfB>p@FTz^XQf^{<&mB zliy@2EGlu|`s>wed5kpkR37g5;7>YNO11L(P8SFqT1Td`RgR!1PG8WDcAoJ{L?X#T zYb#nafOl1}+*etB^hE8TVb!Hb0^Na4l{1L80{o|SAM9S02cM4FTzuD=3%I8hj{opo ztz@5G>01}bRmJR%l1c$hh&BH;$BY&7{*K|~Ls|cHr2U0)N%ws!FmRBm_A^f++5Pe^ z*KY`;q^zxFUnuH<&+{VRjD8*{FeoKHiRr@mcI#)=v6kl9>7B^~1J6N~)l&+DypxdZ zoBVNB#C>n@Nn}ciPDCxNUoF%3r9kYF3=RkJL8LNH%X_!b0+_z5-*aro{p}Y$(q)r= zM2P$+^Vf%UpJx0x{uz0Et)D98ih}wBbcn_m-g~pgnqP4OgEr zGG4#Q6sg}JL_Si!@?L0sIe34+X*}wC7Yf;K$&)Dl$h+ZBy1$(#mrjQZU<15Hsc(K% z-7h_HdXlLRtN-IWs6k(NxLPC$P9J^{aU-=4oqDo$sG*`6i9a1b-?*U?5<>5poTSKs zZT`CHFAk4_mV%W^bA17@-Rjh}`<4wr32PC0!^4EALArbgl55(!^O(!UK2_~Ww zs`+d$*4+i|o%%~n!UITd-94_1@u$clv(Fp0rRA~6TTf!|E%1mYb);xB?iim(}E=d_qPW>GeVr6uUOGK zZxe<_?gl%b4~c-sCif1iY|}vZ=l0y|a}O5w-I%`_o9qJS2puJDnlg|Wr3_!ATn%ek zEKLqw%O>R4P=58{oF3e>`IN-oJlZ6V)BiX0HIZZ5y_9hQy&gj74oj*%yT#I;>^(Yi8d@VlLw!lq?OC=BXGw!;S{2&w5cR0Q7Yj=U1eJ3vs(;$MLIGvK-b-bG) z1l>KMCzdL78wLz0TOINx5Y?{sW{rhFg72qTThx>fmP4QNQR^)dxF40O=o|cJiy+jx z)&IL(9#HR~{QWr080Phbk3@DTGVB{eSo3&iUkBI+#S@SZdJH+LI(u;~!{6(3B4r}stu(3eu zB=}d@9iF@6yk4=A@279GWe(g5Mq8e%?@!N;g%(k&zDw5PDB9`} zS$?FS@S!``rxz4)s|$;b{ia!uA)`N)@4$uUkY>WqrLUC@Qqmavwuf0TPIi;INYnx- zzZRRRM~i_lb?D9m;i)_+qZWrmXmFoC_m#A>uV}aqz^4tr1s|UxcH| zsH(p3g7SgMuO1Sp%EMzv&49nK7)zV{Ud~2njXmadouvZO_*|Z|M6|;IZNqcb)!De! z#nGR$Nol9i)yKFeLO#xuB&Xw)g0ZNlbXqI=zT2I7T9PdYk?Yn>&5TDt2`{_ZQwuTl zHFlA2yRfhD66Ku_%+|xu_d2sIxT+LVSEaUnGpi#|@%_Bx4)ueL#dh#0*C57mk0|)K zJna64TPqRwJ=dZv!AKtWe46pOloI`DP_Ecc!`>~4Zu%R9*Y8{IHXoxckG?l4h1QEv z<->|~P%jdw+Uk@6+HGNyzPJ@s|NI^fF&9_((CwgSx$`(dPn`aWJaDe7&<7a{=yOVr z(Za={+>~H5c{IkQK-HV=C%h&0Tlsd%9w@FjmA;dr9Ox=fe?5J$2A9{EpAjb|3IEyW z72AAM+>FeCEbgsSxuyicPMjX_hGf3W(g%eqJ#u{77X?|54J4vTB+)?6mya$*e!?|e zmgPRUXYD>Ip0Tr0Ed`m%^(B|5YN6JWrA|;i3*y=AUP$_60^6!*DeUU5xWC2~&}7pU z1@jEOw=K>qMgBb{`jky_6`8g-P_UOo!tI6ibqtwGsGhEt~upRpy#DENavg=hh5!A`AzAHaaC9Q^V{(fzq)wpPAjsVT#{fXPXAV(l-{`4 zA5mQE(srAS2K{^UjqNOFka&ji*&SZ=khOEmmQL>sXtQ12wvTmrHxY0%A*sYZgD{Wd zdtG2zWu=p}Ea2_L=}w;vCLgHyq55@4P6{f=K+Ggkz96lHM!G~U4uA0>_#d813l7E66Z=!q8opUu`?cP};wnmRa(fEEVD?I5kk23}R^H@qGabXxf9>Bk7%V`sBF-)TSjYdOH;A@VOwVNE);7OZmB=B@Ep1T>Msvn_FHf#i^OI(4 zY&1XlAU{Kv(_}|uAv4dJv&~u!VbtpByKZ?BB0q6VApE3VIf$CB3)OJB59P>we|E+)vJPdP)lL zdg76%+*R^|LEIZz_3XG${Voa?E*w885T}ETcyinir+NtQ46-?Mu{sN^2hZ(u5`BWB z{(9lZIbPyxp4BwYbe>KI_qqROP3z@i;cH;*jzSbVf(j#JAU5)|sP_YPWVMT$nvVO5 zaI;3{YzVH`8LSyg*zH>eeou0DufF|&Fw#0sKQ(m2jl4?M&mm^G-NdB8=(st%B;fVL zeeeD5KQGPfiI|nwp6A#Q1L_PVEa$j%P_@d8Ml!9N@HQ!VVQ-WT8=qfcIkK?q`$BPd z{rLmgaKThRTE8HJu&(2+sIQ3P*p@X;4hg_gL;$K(Vzr$#N)@jY%sPTMJ+piJCN0<$H4&vMY3sc~%X{wdt;d%yI= z>HQcxpJ%*3;*Q^!M0PM1F5byNk77nep6lIrI^WL;>_H;hKl1+O`lau~vCF*59GXTt`&k zJxL#=$ySmdQWOhc=@oB=8OEaGg%A(>5!~hUpv==hd=GF-n<FDBT1+T81}i#o4LZ;Cad(VM~L(_bgc2+*eA7k@Yh7W~*IA|y@-G!QP)ZZ@#3S!fes^Qo z`%vVxULU)(7!N(4U+W&1*Fm2)H|t&=w-H`@UoEuu?J&Gm_P85I`3SU&EEs5#O5kgc zZrA=pv|{9M-{WD-3iU;~{2Cj$OK!t@ z7VhwIyR&E`-^q<}C#-~nrKeAsGUfpDSJvil?<&CS%I(~vdz~SpcYjDR?&n?n^%h^Js2H+sM{?g*iYd{gZHqOQlLbu2X-3HvwIr5uS;E$P9)y1^Xd{`55i(;Ot z9AfyIK_VOc?=T+a^&s^5&o70S;)I`pDd#?_ z6QeC~Ig9_9CGKfD87X!Y2Ya&N236qA^JiS4pFF5v`kf{ykyJJat{C|f+hvORoPy9E zi7UD$PKn^Y+oofE>1A}8|D65iv~$8s`Zt?Waci5*8nVXL=QVKV?B!%zs(N_&f>kF7 z*I|$I3DU9NDc7Kp3*opg|=r8+WEY=wVNgHX}q3^>hGM$$i?L!dM43%g#dC=EI6&#bL*nITSZ zDgEYms_`09+|SbdY4#4(2=R7Mm4={`S`P7oC5FQ6LEZFFlL405a*4D?%XD0|CGITF z0OixyyQFaA8b8Mk7K`7lG&==19$tU8a<|dh5WdjQ5QuyOM<43oM$}1fvYM~6ej-U% z#pXqn3<;}?u%vSi$NL+A&wEC2+rby$6>a%c09Vsn%|=x_qBEhHYoA=boh!^sU!9db zs0C>sF(V0FW{CTKH)8V(#xriHx%5$)080|kGv*0AV)sJ}#{*S2a~cSHYOQK&GtPo% zSM|1jo_GQ(j|*)YH}ZqtvjFWPtt@am`-&-P?{eh5SC85TtQkl#@*XR425ku4Ol;&%mHdE7T4ERm7u$}gN1!O8P4gX z7@Zf(fmMn;HGC1vRd;Bll+qTYU=8=Atl!*9zSB5dQ(buOj4weJ_R+GOJp6R>%jIMCip0_nALvXm_Ge;Ar37XATU3weEeZskdblu9&vSR@Y60j=4ZR*NE(AMqI`?))cAg!9=)-#_&Z>dC(C|2~H)ypYvfL(N zvG1k6u&s!}PBPrp{|{QPSu8r1YyO;Vim8oi9(=p;HF4nRvhSSI)rHi!|FvYiA-Zo6 z_e?k9bQT@{&2pT9$lNNO?}3t;@LRG23W|61MM7!1r>D2xLwRZ1=S~VPKgqwR;X&T^ z<^jNXochln&Vv?T(@gz?g>W_b&Ae7Ef~dwvJ1-u^U8nzTCr*!9w;_}MsVnMNQm9sl zN`QiC)iPDvU?fNyZNo@lo6Y@T^>1cBZ=LFPxHI&~~$Oj7BOxZn6MZkM|oDr+xgS|6qo6Md$Lxl+K;^LQ_ z6(YZBDT+Q}RVu2eHs30D5ZAMA?vYNAyN1jXnhP&C>p<9PKCcSzJUAGBHM8_&3wU1n_e^C^n0hQ5l?SY~z{TAu$v@^E!fHQhLV>14| zD-q;Lf*w(xltg4eC8K6$ARHnoyhoKM3xo?ux8Lr34Y?l~-nfm&pkk zfePHdu7kRecWMzmZfI0rv)$&u>^78#uw5<=c4&H0^GbrIV&=G9JSB|$;x-> z3-gwWN1nKq3-_IT$l?wRLyla^tP@m{?51;>;QU#-l74p{gr1>1ncB(>g$rUM zH|A{NxywvK3H?gzy3p}6*&-EIX!~@*x<|Vbps<)!Vr)C0+oSXfQsVl;^r{oF!&?e~ zS8VOum}?`Dr7X5-wsQn##Sab9vu1;}pPbN4KrXbhNm`K|*$r9T53ROU+k?O0j%(|( zR-y}tzq44}Xolpbv{++g<3Tk1eX7bi0c7?tQb{yNPx$jxNTPd0A*8GcejYtC3Yi+! z#oqENSu`nr?iMe439WIl>XdKPg!2O09I2ST28`r%Qr9 z`qo9N|8Xo1o?E@4tw=wB1cO*r1zGfjU-}3nZ)Yk7w$@DrQGk2urS=Kkb#Y_B_R1&k z(4!nUm!isjPa_AcyheQZjE=zVYxe^A%y7?I&Ue?&eZ3OBR{YuCZ+c2dyHnP?fIk-U z;|ktr4RazY{Ur56erJWBu)az!iCq3aQ_}S7$XQ6o!N%oVD3wYsJhaG@By-6I{>L>_ zl&1vYRtd}2hu})k{}}J4E#$b{vHq{T5x*U4*Q02f`J)8 za>5wJ9j+3L;gN=*it!W7;!Ml=Ox8+_A6sOUsI0=S7Tuv8m9NCSMiM?&1y*4-&E6z* zLRHxQ3PaB`<<(ePLo%PObrm+_m}8VDSB>d-WX#R(tj2Oq+EP`->rYSV5uv&%`Gfo_bw)$bb zT3Z$7ejAJQps&TuzGU>ZQB`4H-D;78hLzYT)2`DU5s$I=c6sWhY86;`VgG?qvL_g? zKOfWjg>q~{R+#nPp$aTufRwGmycFA1A-y?~x&#~E{Jq<;u>>oZihN?J`2h2=o~vI) zaUauY$+)6Vejm%K^L5%gmxmoycWAn{BOB9+R;m#7&BAOv_{MV>vM}|cmQxGQ?_f&K ztObs*RdRm6jvKX zPi%(y+b!8zUo7X%lQ#cr?pW|xeQe#*6^w5~<1vro-w!RqUq$(Y^w?=+=FJd-{lWDA eit5jk+yC5-^(!I&b34}WBK-gMb}Z{K-v0%@uvfnT diff --git a/tests/test_txfunc.py b/tests/test_txfunc.py index b0ffa541f..b9e1e25be 100644 --- a/tests/test_txfunc.py +++ b/tests/test_txfunc.py @@ -406,9 +406,6 @@ def test_tax_data_sample(): assert isinstance(df, pd.DataFrame) -@pytest.mark.local -# mark as local run since results work on Mac, but differ on other -# platforms def test_tax_func_loop(): """ Test txfunc.tax_func_loop() function. The test is that given inputs from @@ -434,7 +431,8 @@ def test_tax_func_loop(): numparams, tpers, ) = input_tuple - tax_func_type = "DEP" + tax_func_type = "HSV" + numparams = 2 # Rename and create vars to suit new micro_data var names micro_data["total_labinc"] = ( micro_data["Wage income"] + micro_data["SE income"] @@ -482,7 +480,6 @@ def test_tax_func_loop(): output_dir, numparams, ) - expected_tuple = utils.safe_read_pickle( os.path.join(CUR_PATH, "test_io_data", "tax_func_loop_outputs.pkl") ) From eb0708186d5f7dd36e4f38a1f9d3e2e5cd559bb1 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 12:58:37 -0400 Subject: [PATCH 07/17] format --- ogcore/txfunc.py | 12 ++---- tests/test_txfunc.py | 96 ++++++++++++++++++++++++++++++-------------- 2 files changed, 69 insertions(+), 39 deletions(-) diff --git a/ogcore/txfunc.py b/ogcore/txfunc.py index 1d98edbb5..b1f73d0a2 100644 --- a/ogcore/txfunc.py +++ b/ogcore/txfunc.py @@ -633,14 +633,10 @@ def txfunc_est( # ''' # if Atil_init not exist, set to 1.0 if params_init is None: - # Atil_init = 1.0 - # Btil_init = 1.0 - # Ctil_init = 1.0 - # Dtil_init = 1.0 - Atil_init = 0.001 - Btil_init = 0.001 - Ctil_init = 0.001 - Dtil_init = 0.001 + Atil_init = 1.0 + Btil_init = 1.0 + Ctil_init = 1.0 + Dtil_init = 1.0 max_x_init = np.minimum( txrates[(df["total_capinc"] < y_20pctl)].max(), MAX_ETR + 0.05 ) diff --git a/tests/test_txfunc.py b/tests/test_txfunc.py index b9e1e25be..eb7174f30 100644 --- a/tests/test_txfunc.py +++ b/tests/test_txfunc.py @@ -224,15 +224,15 @@ def test_replace_outliers(): ) -@pytest.mark.local # only marking as local because platform -# affects results from scipy.opt that is called in this test - so it'll -# pass if run on Mac with MKL, but not necessarily on other platforms @pytest.mark.parametrize( "rate_type,tax_func_type,true_params", [ - # ("etr", "DEP", DEP_params), - ("etr", "DEP", [6.28E-12, 4.36E-05, 1.04E-23, 7.77E-09, 0.80, 0.80, 0.84, -0.14, -0.15, 0.15, 0.16, -0.15]), - ("etr", "DEP_totalinc", [6.28E-12, 4.36E-05, 0.35, -0.14, 0.15, -0.15]), + # ("etr", "DEP", [6.28E-12, 4.36E-05, 1.04E-23, 7.77E-09, 0.80, 0.80, 0.84, -0.14, -0.15, 0.15, 0.16, -0.15]), + ( + "etr", + "DEP_totalinc", + [6.28e-12, 4.36e-05, 0.35, -0.14, 0.15, -0.15], + ), ("etr", "GS", [0.35, 0.25, 0.03]), ("etr", "linear", [0.25]), ("mtrx", "linear", [0.4]), @@ -241,18 +241,19 @@ def test_replace_outliers(): ("mtrx", "HSV", [0.5, 0.1]), ("mtry", "HSV", [0.4, 0.15]), ], - # ids=["DEP", "DEP_totalinc", "GS"], - ids=["DEP", "DEP_totalinc", "GS", "linear, etr", - "linear, mtrx", - "linear, mtry", - "HSV, etr", - "HSV, mtrx", - "HSV, mtry" - ], + # ids=["DEP", "DEP_totalinc", "GS", "linear, etr", + ids=[ + "DEP_totalinc", + "GS", + "linear, etr", + "linear, mtrx", + "linear, mtry", + "HSV, etr", + "HSV, mtrx", + "HSV, mtry", + ], ) -def test_txfunc_est( - rate_type, tax_func_type, true_params, tmpdir -): +def test_txfunc_est(rate_type, tax_func_type, true_params, tmpdir): """ Test txfunc.txfunc_est() function. The test the estimator can recover (close to) the true parameters. @@ -270,38 +271,71 @@ def test_txfunc_est( "total_capinc": y, "total_labinc": x, "weight": weights, - "total_tax": (( - txfunc.get_tax_rates( - true_params, - x, y, weights, tax_func_type, - rate_type="etr", - for_estimation=False) + eps1) * (x + y)), + "total_tax": ( + ( + txfunc.get_tax_rates( + true_params, + x, + y, + weights, + tax_func_type, + rate_type="etr", + for_estimation=False, + ) + + eps1 + ) + * (x + y) + ), "etr": ( txfunc.get_tax_rates( true_params, - x, y, weights, tax_func_type, + x, + y, + weights, + tax_func_type, rate_type="etr", - for_estimation=False) + eps1), + for_estimation=False, + ) + + eps1 + ), "mtr_labinc": ( txfunc.get_tax_rates( true_params, - x, y, weights, tax_func_type, + x, + y, + weights, + tax_func_type, rate_type="mtr", - for_estimation=False) + eps2), + for_estimation=False, + ) + + eps2 + ), "mtr_capinc": ( txfunc.get_tax_rates( true_params, - x, y, weights, tax_func_type, + x, + y, + weights, + tax_func_type, rate_type="mtr", - for_estimation=False) + eps3), + for_estimation=False, + ) + + eps3 + ), } ) micro_data["age"] = 44 micro_data["year"] = 2025 output_dir = tmpdir param_est, _, obs, _ = txfunc.txfunc_est( - micro_data, 44, 2025, rate_type, tax_func_type, - len(true_params), output_dir, True + micro_data, + 44, + 2025, + rate_type, + tax_func_type, + len(true_params), + output_dir, + True, ) assert obs == micro_data.shape[0] From a93e6eb8b51751e249d2f3c6f4d909efab41bab8 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 13:06:39 -0400 Subject: [PATCH 08/17] pin marshmallow --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index 68fadd400..b89fb92ce 100644 --- a/environment.yml +++ b/environment.yml @@ -15,6 +15,7 @@ dependencies: - distributed>=2.30.1 - paramtools>=0.15.0 - sphinx>=3.5.4 +- marshmallow<4.0.0 - sphinx-argparse - sphinxcontrib-bibtex>=2.0.0 - sphinx-math-dollar From e6c91b8d55c24d9425f0e448f6d0f3b2fb9db0a2 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 17:28:49 -0400 Subject: [PATCH 09/17] remove unused tests --- tests/test_txfunc.py | 61 ++------------------------------------------ 1 file changed, 2 insertions(+), 59 deletions(-) diff --git a/tests/test_txfunc.py b/tests/test_txfunc.py index eb7174f30..179443939 100644 --- a/tests/test_txfunc.py +++ b/tests/test_txfunc.py @@ -228,10 +228,8 @@ def test_replace_outliers(): "rate_type,tax_func_type,true_params", [ # ("etr", "DEP", [6.28E-12, 4.36E-05, 1.04E-23, 7.77E-09, 0.80, 0.80, 0.84, -0.14, -0.15, 0.15, 0.16, -0.15]), - ( - "etr", - "DEP_totalinc", - [6.28e-12, 4.36e-05, 0.35, -0.14, 0.15, -0.15], + ("etr", "DEP_totalinc", + [6.28e-12, 4.36e-05, 0.35, -0.14, 0.15, -0.15], ), ("etr", "GS", [0.35, 0.25, 0.03]), ("etr", "linear", [0.25]), @@ -349,61 +347,6 @@ def test_txfunc_est(rate_type, tax_func_type, true_params, tmpdir): # 95% confidence interval of the estimated parameters -@pytest.mark.parametrize( - "rate_type,tax_func_type,numparams,expected_tuple", - [ - ("etr", "linear", 1, expected_tuple_linear), - ("mtrx", "linear", 1, expected_tuple_linear_mtrx), - ("mtry", "linear", 1, expected_tuple_linear_mtry), - ("etr", "HSV", 2, expected_tuple_HSV), - ("mtrx", "HSV", 2, expected_tuple_HSV_mtrx), - ("mtry", "HSV", 2, expected_tuple_HSV_mtry), - ], - ids=[ - "linear, etr", - "linear, mtrx", - "linear, mtry", - "HSV, etr", - "HSV, mtrx", - "HSV, mtry", - ], -) -def test_txfunc_est_on_GH( - rate_type, tax_func_type, numparams, expected_tuple, tmpdir -): - """ - Test txfunc.txfunc_est() function. The test is that given - inputs from previous run, the outputs are unchanged. - """ - micro_data = utils.safe_read_pickle( - os.path.join(CUR_PATH, "test_io_data", "micro_data_dict_for_tests.pkl") - ) - s = 80 - t = 2030 - df = txfunc.tax_data_sample(micro_data[str(t)]) - output_dir = tmpdir - # Put old df variables into new df var names - df.rename( - columns={ - "MTR labor income": "mtr_labinc", - "MTR capital income": "mtr_capinc", - "Total labor income": "total_labinc", - "Total capital income": "total_capinc", - "ETR": "etr", - "expanded_income": "market_income", - "Weights": "weight", - }, - inplace=True, - ) - test_tuple = txfunc.txfunc_est( - df, s, t, rate_type, tax_func_type, numparams, output_dir, True - ) - - for i, v in enumerate(expected_tuple): - print("For element", i, ", test tuple =", test_tuple[i]) - assert np.allclose(test_tuple[i], v, rtol=0.0, atol=1e-04) - - def test_txfunc_est_exception(tmpdir): micro_data = utils.safe_read_pickle( os.path.join(CUR_PATH, "test_io_data", "micro_data_dict_for_tests.pkl") From 3388a048c29862ac37e0de7c004472aa7609dc85 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 18:06:25 -0400 Subject: [PATCH 10/17] format --- tests/test_txfunc.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/test_txfunc.py b/tests/test_txfunc.py index 179443939..0d9735dd3 100644 --- a/tests/test_txfunc.py +++ b/tests/test_txfunc.py @@ -228,8 +228,10 @@ def test_replace_outliers(): "rate_type,tax_func_type,true_params", [ # ("etr", "DEP", [6.28E-12, 4.36E-05, 1.04E-23, 7.77E-09, 0.80, 0.80, 0.84, -0.14, -0.15, 0.15, 0.16, -0.15]), - ("etr", "DEP_totalinc", - [6.28e-12, 4.36e-05, 0.35, -0.14, 0.15, -0.15], + ( + "etr", + "DEP_totalinc", + [6.28e-12, 4.36e-05, 0.35, -0.14, 0.15, -0.15], ), ("etr", "GS", [0.35, 0.25, 0.03]), ("etr", "linear", [0.25]), From 38afaba2696ab90fa1b1e0775f7ac9857bcfebec Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Mon, 21 Apr 2025 18:54:52 -0400 Subject: [PATCH 11/17] make test local only --- tests/test_txfunc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_txfunc.py b/tests/test_txfunc.py index 0d9735dd3..6f0a34869 100644 --- a/tests/test_txfunc.py +++ b/tests/test_txfunc.py @@ -726,6 +726,7 @@ def test_get_tax_rates( assert np.allclose(test_txrates, expected) +@pytest.mark.local def test_tax_func_estimate(tmpdir, dask_client): """ Test txfunc.tax_func_loop() function. The test is that given From 76564351d91037b6fc15b2bbaccb84e42d7ea1c1 Mon Sep 17 00:00:00 2001 From: Richard Evans Date: Tue, 22 Apr 2025 13:31:11 -0600 Subject: [PATCH 12/17] Updated version in setup.py, __init__.py, and CHANGELOG.md --- CHANGELOG.md | 10 ++++++++++ ogcore/__init__.py | 2 +- setup.py | 2 +- 3 files changed, 12 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a93fede26..69632ffe1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [0.14.3] - 2025-04-22 16:00:00 + +### Added + +- Puts a ceiling on the version of the `marshmallow<4.0.0` package in `environment.yml` +- Update `txfunc.py` for the estimation of the `HSV` and `GS` tax functions. +- Update `test_txfunc.py`, `tax_func_estimate_outputs.pkl`, and `tax_func_loop_outputs.pkl` files for testing + ## [0.14.2] - 2025-04-04 12:00:00 ### Added @@ -373,6 +381,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Version [0.7.0] on August 30, 2021 was the first time that the OG-USA repository was detached from all of the core model logic, which was named OG-Core. Before this version, OG-USA was part of what is now the [`OG-Core`](https://github.com/PSLmodels/OG-Core) repository. In the next version of OG-USA, we adjusted the version numbering to begin with 0.1.0. This initial version of 0.7.0, was sequential from what OG-USA used to be when the OG-Core project was called OG-USA. - Any earlier versions of OG-USA can be found in the [`OG-Core`](https://github.com/PSLmodels/OG-Core) repository [release history](https://github.com/PSLmodels/OG-Core/releases) from [v.0.6.4](https://github.com/PSLmodels/OG-Core/releases/tag/v0.6.4) (Jul. 20, 2021) or earlier. + +[0.14.3]: https://github.com/PSLmodels/OG-Core/compare/v0.14.2...v0.14.3 [0.14.2]: https://github.com/PSLmodels/OG-Core/compare/v0.14.1...v0.14.2 [0.14.1]: https://github.com/PSLmodels/OG-Core/compare/v0.14.0...v0.14.1 [0.14.0]: https://github.com/PSLmodels/OG-Core/compare/v0.13.2...v0.14.0 diff --git a/ogcore/__init__.py b/ogcore/__init__.py index e733effbf..bbcb58e1a 100644 --- a/ogcore/__init__.py +++ b/ogcore/__init__.py @@ -20,4 +20,4 @@ from ogcore.txfunc import * from ogcore.utils import * -__version__ = "0.14.2" +__version__ = "0.14.3" diff --git a/setup.py b/setup.py index 4ad5c2b44..f04f572d2 100755 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="ogcore", - version="0.14.2", + version="0.14.3", author="Jason DeBacker and Richard W. Evans", license="CC0 1.0 Universal (CC0 1.0) Public Domain Dedication", description="A general equilibrium overlapping generations model for fiscal policy analysis", From c995cf2ec19f3fc63bb9c63efbf8e153bd2b7e7b Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Thu, 24 Apr 2025 08:43:36 -0400 Subject: [PATCH 13/17] remove num_workers from compute --- ogcore/SS.py | 2 +- ogcore/TPI.py | 2 +- ogcore/txfunc.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ogcore/SS.py b/ogcore/SS.py index a8c58865e..13773438d 100644 --- a/ogcore/SS.py +++ b/ogcore/SS.py @@ -272,7 +272,7 @@ def inner_loop(outer_loop_vars, p, client): ) ) if client: - futures = client.compute(lazy_values, num_workers=p.num_workers) + futures = client.compute(lazy_values) results = client.gather(futures) else: results = results = compute( diff --git a/ogcore/TPI.py b/ogcore/TPI.py index b9c123ba8..318c39d68 100644 --- a/ogcore/TPI.py +++ b/ogcore/TPI.py @@ -795,7 +795,7 @@ def run_TPI(p, client=None): ) ) if client: - futures = client.compute(lazy_values, num_workers=p.num_workers) + futures = client.compute(lazy_values) results = client.gather(futures) else: results = compute( diff --git a/ogcore/txfunc.py b/ogcore/txfunc.py index b1f73d0a2..a6df207fb 100644 --- a/ogcore/txfunc.py +++ b/ogcore/txfunc.py @@ -1687,7 +1687,7 @@ def tax_func_estimate( ) ) if client: - futures = client.compute(lazy_values, num_workers=num_workers) + futures = client.compute(lazy_values) results = client.gather(futures) else: results = results = compute( From d2e23f1c15b1a783ebfe9cc681ef11fd7d0ba7b8 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Thu, 24 Apr 2025 15:21:54 -0400 Subject: [PATCH 14/17] update scattering --- ogcore/SS.py | 75 ++++++++++++++++++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/ogcore/SS.py b/ogcore/SS.py index 13773438d..6705b2e4e 100644 --- a/ogcore/SS.py +++ b/ogcore/SS.py @@ -34,12 +34,6 @@ level=log_level, format="%(message)s" # Only show the message itself ) -""" -A global future for the Parameters object for client workers. -This is scattered once and place at module scope, then used -by the client in the inner loop. -""" -scattered_p = None """ ------------------------------------------------------------------------ @@ -223,12 +217,6 @@ def inner_loop(outer_loop_vars, p, client): units """ - # Retrieve the "scattered" Parameters object. - global scattered_p - - if scattered_p is None: - scattered_p = client.scatter(p, broadcast=True) if client else p - # unpack variables to pass to function bssmat, nssmat, r_p, r, w, p_m, Y, BQ, TR, Ig_baseline, factor = ( outer_loop_vars @@ -247,35 +235,71 @@ def inner_loop(outer_loop_vars, p, client): tr = household.get_tr(TR, None, p, "SS") ubi = p.ubi_nom_array[-1, :, :] / factor + scattered_p = client.scatter(p, broadcast=True) if client else p + lazy_values = [] for j in range(p.J): guesses = np.append(bssmat[:, j], nssmat[:, j]) - euler_params = ( + + # Create a delayed function that will access the scattered_p + @delayed + def solve_for_j( + guesses, r_p, w, p_tilde, - bq[:, j], - rm[:, j], - tr[:, j], - ubi[:, j], + bq_j, + rm_j, + tr_j, + ubi_j, factor, j, - scattered_p, - ) - lazy_values.append( - delayed(opt.root)( + scattered_p_future, + ): + # This function will be executed on workers with access to scattered_p + return opt.root( euler_equation_solver, guesses * 0.9, - args=euler_params, + args=( + r_p, + w, + p_tilde, + bq_j, + rm_j, + tr_j, + ubi_j, + factor, + j, + scattered_p_future, + ), method=p.FOC_root_method, tol=MINIMIZER_TOL, ) + + # Add the delayed computation to our list + lazy_values.append( + solve_for_j( + guesses, + r_p, + w, + p_tilde, + bq[:, j], + rm[:, j], + tr[:, j], + ubi[:, j], + factor, + j, + scattered_p, + ) ) + if client: + # Compute all the values futures = client.compute(lazy_values) + # Later, gather the results when needed results = client.gather(futures) else: - results = results = compute( + results = compute( *lazy_values, scheduler=dask.multiprocessing.get, num_workers=p.num_workers, @@ -1197,11 +1221,6 @@ def run_SS(p, client=None): results """ - global scattered_p - if client: - scattered_p = client.scatter(p, broadcast=True) - else: - scattered_p = p # Create list of deviation factors for initial guesses of r and TR dev_factor_list = [ From b05deb8237c1a7beab01acadfb0a904b471e3c21 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Thu, 24 Apr 2025 21:44:24 -0400 Subject: [PATCH 15/17] scatter in tpi --- ogcore/TPI.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/ogcore/TPI.py b/ogcore/TPI.py index 318c39d68..c7241e6bb 100644 --- a/ogcore/TPI.py +++ b/ogcore/TPI.py @@ -46,13 +46,6 @@ level=log_level, format="%(message)s" # Only show the message itself ) -""" -A global future for the Parameters object for client workers. -This is scattered once and place at module scope, then used -by the client in the inner loop. -""" -scattered_p = None - def get_initial_SS_values(p): """ @@ -578,12 +571,6 @@ def run_TPI(p, client=None): results """ - global scattered_p - if client: - scattered_p = client.scatter(p, broadcast=True) - else: - scattered_p = p - # unpack tuples of parameters initial_values, ss_vars, theta, baseline_values = get_initial_SS_values(p) (B0, b_sinit, b_splus1init, factor, initial_b, initial_n) = initial_values @@ -779,10 +766,15 @@ def run_TPI(p, client=None): ).sum(axis=2) p_tilde = aggr.get_ptilde(p_i[:, :], p.tau_c[:, :], p.alpha_c, "TPI") + # scatter parameters to workers + scattered_p = client.scatter(p, broadcast=True) if client else p + euler_errors = np.zeros((p.T, 2 * p.S, p.J)) lazy_values = [] for j in range(p.J): guesses = (guesses_b[:, :, j], guesses_n[:, :, j]) + + # Add the delayed computation to our list lazy_values.append( delayed(inner_loop)( guesses, @@ -791,11 +783,13 @@ def run_TPI(p, client=None): ubi, j, ind, - scattered_p, + scattered_p ) ) if client: + # Compute all the values futures = client.compute(lazy_values) + # Later, gather the results when needed results = client.gather(futures) else: results = compute( From 3a8dccda1e53e01300c70c1a0e9dbfd173280398 Mon Sep 17 00:00:00 2001 From: Jason DeBacker Date: Thu, 24 Apr 2025 21:50:24 -0400 Subject: [PATCH 16/17] format --- ogcore/TPI.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ogcore/TPI.py b/ogcore/TPI.py index c7241e6bb..d7c066c2a 100644 --- a/ogcore/TPI.py +++ b/ogcore/TPI.py @@ -783,7 +783,7 @@ def run_TPI(p, client=None): ubi, j, ind, - scattered_p + scattered_p, ) ) if client: From c58f3a0167922c486c6c7d2e86a89afc29e6d269 Mon Sep 17 00:00:00 2001 From: Richard Evans Date: Fri, 25 Apr 2025 08:38:31 -0600 Subject: [PATCH 17/17] Update CHANGELOG.md --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 69632ffe1..6a37f1c4c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,13 +6,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.14.3] - 2025-04-22 16:00:00 +## [0.14.3] - 2025-04-25 10:00:00 ### Added - Puts a ceiling on the version of the `marshmallow<4.0.0` package in `environment.yml` - Update `txfunc.py` for the estimation of the `HSV` and `GS` tax functions. - Update `test_txfunc.py`, `tax_func_estimate_outputs.pkl`, and `tax_func_loop_outputs.pkl` files for testing +- Update `dask` and `distributed` client calls in `SS.py` and `TPI.py` to allow for updated versions. ## [0.14.2] - 2025-04-04 12:00:00