From 3d6ec7b3979b6f6dff20652b1b4c4473b9a2678c Mon Sep 17 00:00:00 2001 From: ppetrak Date: Mon, 19 Dec 2022 14:03:45 +0100 Subject: [PATCH] fix: flak8 & pytest --- examples/my_db.sqlite | Bin 0 -> 8192 bytes pyerrors/__init__.py | 1 - pyerrors/fits.py | 112 ++++++++++++++++++++---------------------- 3 files changed, 54 insertions(+), 59 deletions(-) create mode 100644 examples/my_db.sqlite diff --git a/examples/my_db.sqlite b/examples/my_db.sqlite new file mode 100644 index 0000000000000000000000000000000000000000..880c4275a7a77d3ddbbb866b9ecba01ecb3bccd4 GIT binary patch literal 8192 zcmeI%dpOhm{|E35bC}bdA~Z9{l|rp-l_o4~l4H(?InH4iBh#I;+&Pm|A(W8J`B0A4 za+cWK6e&f~iE^lruivQq`*Z*P`F;QXUf$QX&u8yz*FM+vxL()u_4;gPCWf9AH<-I$ z0Es|>9RWZ9Kp;R11_J;9Y`>qZ!}bEQUO9fBZymPY|9^`Oa7gVqKkHLKK}!JZvAcE! zb_I3?b_I3?b_I3?b_I3?b_I3?b_I3?{+k3M6@gshva-M^aSFlN$Bje?40QZIXKXsA z*wf}%nE7dKLo5vbzaN3al)-R{V<3ePK!L*yjLflmSW^`+47R=L=l{<=R|19L@8{`D z3H;~23(+%>97G6kb9AOUy85}pVXVs;Zhy_zrQk4gtd%(!v2||_X2SynZt-Jz--2U! z;9QS5T|wOJ<*d!!wJY$S3q;w0xqv_*S1b%71HshFk5t&w7iAK-0ZvYZE^M3jDg9n_ zH3&AmlPEPYpf9P}I#41Cwx@U+8hZ^(NKZWLOBf1YK)n>>k$-ak$oCJ366()ubZc|x z(8pBhtin}=%xdpxiA0Xi1LykR4Gc-YyO17TJ2~ax+`mkn0I8 zKvrItLe0$QOC6hibt%60UhwW0!YLQ}OsDw~|A*d27j3FbNog-{(=O`Kj?v5$ji6VSeYO*ESlk{O`h zeP^Lv0~HhaGj8HpX|)mO2cI(tCOQLQQNNN?!HxE;KMcyGEIA7HA(^bdPs`H4ljZigF z%<5uj7fNYJn)pZ>4L{u{O5pv)UD52;RR#yRT(#2D zgUHrr~FbrNogGaJK9tbw$?a*8vzbgj96=}HR92w`gQdkg^`13Pq6&TZ*-cHnv% z3{qWf8!vux$qGjnuj3q7M7_uKwXe_29P_E(2s_7Ex+`tEdG*nKyKMTAl9AGy z!Y*XB*EFr7&(_do%#p(-da+$*-Dgsej4Q+v>7rO!Sh-1~s&N%?w5Jmc(NgIS%X_@| zt^N5Q7M*$*?w|5v!jQ%~%&zDcsUQ7)t;YnrFQVrsDHWig*QOp7#b?7M+^m$xCM?mK zaVJ-n)WE1%kMJiV0S|0Z_8?M`thSy&q;!p}%Z0#ax7)otPj!;VpPE{fpz@j*Kn{=4 z|Duylo+XGjcj5w~O{nNmAFC5C!GC3y9~aSe%RB#4$j19xY*{R^6JQv8rZ@xA208pG z?YT`=n{0l|HENH`0e%gMNdI#W_H`e3S#^8J#*z(G41GwKU>;YrDa=Ct2t{@bG zsc6b76GNTYVVk^FxsaGm+)zW-^S^9Qy)y;JLPL>?$_+%mWRmM!>|FT5`1~Dxm}mR& z;K}D1T??Xv70yh#5@ZQg31|NZFI83HG?w$uSDM@9%5r;tS7&IBvr(ZwQ75^W8xkc) z>d>(OniGRnBLJ@Gq(H^~j}u{=)}LOsybhRc)>V$wdEO;AIrN1U3vVuX2vA-$|M4cK zS?vo5FL?y&0p_hJFHpghD6CAQMCRx&ttveDhu(9lF<#2CT6Z8ZS4SmgwT^mEkhI3q zG2LjH^LLyOMwh-fJe19w#=Ph{6?cu81QO#q*m~rV1nq&rY>?T=%_-g@E7Iv(D}8Cb z2lajRy}WL=LL&Nb#o7&P>#oBH(!MK$$_KB^9o1vADN5ofqBSr?)D8fak-fZkGe(`- z)Xpx|(H!%&S68-k!JhM%E?huH7*maFdj7*XVu$XI*zF%kbjZqno}pcoa43lfcTh6g zMaq^cw6?IXy9FLujJ_~$r`d^2!33fbaE2+Ur_-&j<~eM^q}F&lyXjWS^aBarneU&G zl?=KllkYk&Nk8oXhpsSnoi7MvT=XP30-%CU!aMg)tcPu8ddN)r-^ZK5L}5k5t<-fx zHJ$Ij=qdQei)KoCVzXM0*n!5qDZGNbH_Hp)nk7LWX32!1xb6^Zr3zqGQHNT25lq|T zsC68%7fpF`;IsQsL&Z6=L7XK*JTPx<{$jKasyqB%(j_o2Q!-FMU_9bNAs+OT47%;P zcs<}gT)Azma>U3J?vm$Rd`wrxyl)XpK4)OKd`8rZwogd8ir;3?CU*>0D_K8#gMs!C zt@I84c>3GYpMxPB#f*vV+7=9}l00Ab;4Hr|y6`4@ zvbXA2!X-A(F=ZM)c}`skcBo!fuTnqN6z14*x*X_*%oS&~HR*12-OLK9WkW?{VbZ^y z$6nSPR_HYyDTjpT+337^H*llCJ3A%HM4Yq6sgu-2A8;=!%hgL;V7^%+=j>g6ed>mI zoZAEZ+|xwKvSM9qA6jkPQ(<`i2D6gCIco(yNZnM)#KxPwxcy-rgWq_V(BQ&Rz-at?Xqe%5P_HMsO)1lC<%X$<&nU`q#z8IG=p}CS>vk zidHrw>~NhjxcDir>(1wt_!`{}$8Q_(NdZsjL$6y`0vQ58lSmsvY%qWInWj)6hjDLv zfqQ0aJLy=#>W4PX8kbd6Q@u=<8CqdJFn z_u{cvmlalAG!S!dXqbaM{ZX^e@mS*PbJ1;dUmjt%Ea-xiQD=lOcFLpiLaGC&B$@!~Hb^kv9 zM9=6_d-~6T(yhr1UA6dLXqK9nZJ1ROV=q-jO>_K}5%{f)P<&5_|HyY;gk&Y>A+L`} zZF*{K2@o@FWan_|^|U*m)!Mf3-4Xk3m_NHkO+lN7t?w9g^|Jh<`w9o=lLe%GCqso< z8S!U~ST-voG z(E%C=PJiqJHST$v_U91wZh9QpF&*FW298T(+Xu?ojtGIh4^2n58ywp&Unrl9D5ly` zeCFl?+YaAi2ICx@mZ&-ZI{3J9YB4ulM}xVC8EFhKo&07*w3| zQb=NU1O52a(Pi;7M@4Z6j(S&|u6BJR^=ymSoXu$Iy<+--%mg&6#kKw2;$mYxlhA?! zeW=g*!A4i%0nNSMnf;ltj2EMsgR8(mlhlte{x0e*jut9QmAFEx2FTXmGd;!a55ICQIKeQHo7K6yMe)3Mw&=ZH3-EaS} zGQskEt+yc_lN3{<9s0_MxnFZYUCG#dPCMCL2*#l|-mYcAK0=IV*%}xz(+#$5GduM# z4)b=F++zM^%^}57lF-3?lb`bgk;%3N=kknY&a^ulkV>Yt>1=m%mHJbt`D;`=Zob#* z_M<%B54{i*e3>X3AZ1MNX;Xyx6m-cL(g7p;r= 3.7: Dictionary order is guaranteed to be insertion order. (https://docs.python.org/3/library/stdtypes.html#dict-views) Ensures that x, y and func values are mapped correctly. - + x : dict dict of lists. y : dict @@ -117,16 +116,16 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): import autograd.numpy as anp funcs = {"a": func_a, "b": func_b} - + def func_a(a, x): return a[1] * anp.exp(-a[0] * x) def func_b(a, x): return a[2] * anp.exp(-a[0] * x) - + It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work. - + priors : list, optional priors has to be a list with an entry for every parameter in the fit. The entries can either be Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like @@ -160,10 +159,10 @@ def least_squares(x, y, func, priors=None, silent=False, **kwargs): ''' if priors is not None: return _prior_fit(x, y, func, priors, silent=silent, **kwargs) - - elif (type(x)==dict and type(y)==dict and type(func)==dict): + + elif (type(x) == dict and type(y) == dict and type(func) == dict): return _combined_fit(x, y, func, silent=silent, **kwargs) - + else: return _standard_fit(x, y, func, silent=silent, **kwargs) @@ -688,39 +687,40 @@ def _standard_fit(x, y, func, silent=False, **kwargs): return output -def _combined_fit(x,y,func,silent=False,**kwargs): - + +def _combined_fit(x, y, func, silent=False, **kwargs): + if kwargs.get('correlated_fit') is True: raise Exception("Correlated fit has not been implemented yet") output = Fit_result() output.fit_function = func - + if kwargs.get('num_grad') is True: jacobian = num_jacobian hessian = num_hessian else: jacobian = auto_jacobian hessian = auto_hessian - + x_all = [] y_all = [] for key in x.keys(): - x_all+=x[key] - y_all+=y[key] - + x_all += x[key] + y_all += y[key] + x_all = np.asarray(x_all) - + if len(x_all.shape) > 2: raise Exception('Unknown format for x values') - + # number of fit parameters n_parms_ls = [] for key in func.keys(): if not callable(func[key]): - raise TypeError('func (key='+ key + ') is not a function.') + raise TypeError('func (key=' + key + ') is not a function.') if len(x[key]) != len(y[key]): - raise Exception('x and y input (key='+ key + ') do not have the same length') + raise Exception('x and y input (key=' + key + ') do not have the same length') for i in range(42): try: func[key](np.arange(i), x_all.T[0]) @@ -731,41 +731,41 @@ def _combined_fit(x,y,func,silent=False,**kwargs): else: break else: - raise RuntimeError("Fit function (key="+ key + ") is not valid.") + raise RuntimeError("Fit function (key=" + key + ") is not valid.") n_parms = i n_parms_ls.append(n_parms) n_parms = max(n_parms_ls) if not silent: print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) - + if 'initial_guess' in kwargs: x0 = kwargs.get('initial_guess') if len(x0) != n_parms: raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) else: x0 = [0.1] * n_parms - + output.method = kwargs.get('method', 'migrad') if not silent: print('Method:', output.method) - + def chisqfunc(p): chisq = 0.0 for key in func.keys(): x_array = np.asarray(x[key]) - model = anp.array(func[key](p,x_array)) + model = anp.array(func[key](p, x_array)) y_obs = y[key] y_f = [o.value for o in y_obs] dy_f = [o.dvalue for o in y_obs] - C_inv = np.diag(np.diag(np.ones((len(x_array),len(x_array)))))/dy_f/dy_f - chisq += anp.sum((y_f - model)@ C_inv @(y_f - model)) + C_inv = np.diag(np.diag(np.ones((len(x_array), len(x_array))))) / dy_f / dy_f + chisq += anp.sum((y_f - model) @ C_inv @ (y_f - model)) return chisq - + if output.method == 'migrad': tolerance = 1e-4 if 'tol' in kwargs: tolerance = kwargs.get('tol') - fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef + fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef output.iterations = fit_result.nfev else: tolerance = 1e-12 @@ -776,23 +776,23 @@ def _combined_fit(x,y,func,silent=False,**kwargs): chisquare = fit_result.fun output.message = fit_result.message - + if not fit_result.success: raise Exception('The minimization procedure did not converge.') if x_all.shape[-1] - n_parms > 0: output.chisquare = chisqfunc(fit_result.x) output.dof = x_all.shape[-1] - n_parms - output.chisquare_by_dof = output.chisquare/output.dof + output.chisquare_by_dof = output.chisquare / output.dof output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof) else: output.chisquare_by_dof = float('nan') - + if not silent: print(fit_result.message) - print('chisquare/d.o.f.:', output.chisquare_by_dof ) - print('fit parameters',fit_result.x) - + print('chisquare/d.o.f.:', output.chisquare_by_dof) + print('fit parameters', fit_result.x) + def chisqfunc_compact(d): chisq = 0.0 list_tmp = [] @@ -800,38 +800,37 @@ def _combined_fit(x,y,func,silent=False,**kwargs): c2 = 0 for key in func.keys(): x_array = np.asarray(x[key]) - c2+=len(x_array) - model = anp.array(func[key](d[:n_parms],x_array)) + c2 += len(x_array) + model = anp.array(func[key](d[:n_parms], x_array)) y_obs = y[key] - y_f = [o.value for o in y_obs] dy_f = [o.dvalue for o in y_obs] - C_inv = np.diag(np.diag(np.ones((len(x_array),len(x_array)))))/dy_f/dy_f - list_tmp.append(anp.sum((d[n_parms+c1:n_parms+c2]- model)@ C_inv @(d[n_parms+c1:n_parms+c2]- model))) - c1+=len(x_array) + C_inv = np.diag(np.diag(np.ones((len(x_array), len(x_array))))) / dy_f / dy_f + list_tmp.append(anp.sum((d[n_parms + c1:n_parms + c2] - model) @ C_inv @ (d[n_parms + c1:n_parms + c2] - model))) + c1 += len(x_array) chisq = anp.sum(list_tmp) return chisq - + def prepare_hat_matrix(): hat_vector = [] for key in func.keys(): x_array = np.asarray(x[key]) - if (len(x_array)!= 0): + if (len(x_array) != 0): hat_vector.append(jacobian(func[key])(fit_result.x, x_array)) hat_vector = [item for sublist in hat_vector for item in sublist] return hat_vector - + fitp = fit_result.x - y_f = [o.value for o in y_all] + y_f = [o.value for o in y_all] dy_f = [o.dvalue for o in y_all] - + if np.any(np.asarray(dy_f) <= 0.0): raise Exception('No y errors available, run the gamma method first.') - + try: hess = hessian(chisqfunc)(fitp) except TypeError: raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None - + jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f))) # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv @@ -839,29 +838,26 @@ def _combined_fit(x,y,func,silent=False,**kwargs): deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:]) except np.linalg.LinAlgError: raise Exception("Cannot invert hessian matrix.") - - + if kwargs.get('expected_chisquare') is True: if kwargs.get('correlated_fit') is not True: W = np.diag(1 / np.asarray(dy_f)) cov = covariance(y_all) hat_vector = prepare_hat_matrix() - A = W @ hat_vector #hat_vector = 'jacobian(func)(fit_result.x, x)' + A = W @ hat_vector # hat_vector = 'jacobian(func)(fit_result.x, x)' P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W) output.chisquare_by_expected_chisquare = chisquare / expected_chisquare if not silent: print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare) - result = [] for i in range(n_parms): result.append(derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all), man_grad=list(deriv_y[i]))) - - output.fit_parameters = result - - return output + output.fit_parameters = result + + return output def fit_lin(x, y, **kwargs):