From 9526156e8e7933accb43c3ad23e912f585f34480 Mon Sep 17 00:00:00 2001 From: sudojan Date: Tue, 30 Jun 2020 19:34:06 +0200 Subject: [PATCH 1/7] add Mac DS_Store to gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index d78a011..ec5d674 100644 --- a/.gitignore +++ b/.gitignore @@ -89,3 +89,5 @@ ENV/ # Rope project settings .ropeproject + +.DS_Store \ No newline at end of file From 96905cef8947ac435fcd976de1f92d200cb1148e Mon Sep 17 00:00:00 2001 From: sudojan Date: Tue, 30 Jun 2020 19:35:11 +0200 Subject: [PATCH 2/7] replace for loops in StandardLLH with numpy functions --- funfolding/solution/likelihood.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/funfolding/solution/likelihood.py b/funfolding/solution/likelihood.py index 1417b11..8a27220 100644 --- a/funfolding/solution/likelihood.py +++ b/funfolding/solution/likelihood.py @@ -229,19 +229,11 @@ def evaluate_gradient(self, f): h_unreg -= part_b if self._tau is not None: if self.log_f_reg: - reg_part = np.zeros(self.model.dim_f) denom_f = f_reg + self.log_f_offset nom_f = np.log(denom_f * self.reg_factor_f) ln_10_squared = np.log(10)**2 - pre = np.zeros((self.model.dim_f, - self.model.dim_f)) - for i in range(self.model.dim_f): - for j in range(self.model.dim_f): - pre[i, j] = self._C[i, j] * nom_f[i] - pre[i, j] /= ln_10_squared * denom_f[i] - for i in range(self.model.dim_f): - reg_part[i] = np.sum(pre[i, :]) - reg_part[i] += np.sum(pre[:, i]) + pre = self._C / ln_10_squared * np.broadcast_to(nom_f / denom_f, np.shape(_C)).T + reg_part = np.sum(pre, axis=0) + np.sum(pre, axis=1) else: reg_part = np.dot(self._C, f_reg * self.reg_factor_f) else: From 14a9a26d91d53b7031f3c4fa0d567a14c7fe4de3 Mon Sep 17 00:00:00 2001 From: sudojan Date: Tue, 30 Jun 2020 19:36:08 +0200 Subject: [PATCH 3/7] make example 0 and 3 run again --- examples/03_mcmc_pull.py | 5 ++--- funfolding/visualization/visualize_tree_binning.py | 4 +--- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/examples/03_mcmc_pull.py b/examples/03_mcmc_pull.py index f029c58..7250357 100644 --- a/examples/03_mcmc_pull.py +++ b/examples/03_mcmc_pull.py @@ -77,13 +77,12 @@ def smear(unbinned_f, factor=10, exponent=2, scale=5): binned_f = np.digitize(unbinned_f, binning_f) vec_g, vec_f = model.generate_vectors(binned_g, binned_f) llh = ff.solution.StandardLLH(tau=None, - C='thikonov', - neg_llh=False) + C='thikonov',) llh.initialize(vec_g=vec_g, model=model) sol_mcmc = ff.solution.LLHSolutionMCMC() sol_mcmc.initialize(llh=llh, model=model) sol_mcmc.set_x0_and_bounds() - vec_f_est_mcmc, sigma_vec_f, samples, probs = sol_mcmc.fit() + vec_f_est_mcmc, sigma_vec_f, samples, probs, autocorrelation = sol_mcmc.fit() print(vec_f_est_mcmc) diff --git a/funfolding/visualization/visualize_tree_binning.py b/funfolding/visualization/visualize_tree_binning.py index 68391f0..8cd99e2 100644 --- a/funfolding/visualization/visualize_tree_binning.py +++ b/funfolding/visualization/visualize_tree_binning.py @@ -36,9 +36,7 @@ def plot_hexbins(ax, **hex_kwargs) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="8%", pad=0.05) - cb = matplotlib.colorbar.ColorbarBase(cax, - cmap=cmap, - norm=norm) + cb = cax.figure.colorbar(norm=norm, cmap=cmap) return cb From 9f6013c60ff065440320224c906932342a72eac5 Mon Sep 17 00:00:00 2001 From: sudojan Date: Tue, 30 Jun 2020 20:18:19 +0200 Subject: [PATCH 4/7] make example 5 and 7 run again --- examples/05_fact_example.py | 4 ++-- examples/07_fact_pulls.py | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/05_fact_example.py b/examples/05_fact_example.py index c41f230..9778bc0 100644 --- a/examples/05_fact_example.py +++ b/examples/05_fact_example.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib -matplotlib.use('Agg') +# matplotlib.use('Agg') from matplotlib import pyplot as plt import matplotlib.colors as colors @@ -295,7 +295,7 @@ def generate_acceptance_correction(vec_f_truth, svd_values = tree_model_uniform.evaluate_condition() ax.hist(bin_centers, bins=bin_edges, - ooweights=svd_values, + weights=svd_values, histtype='step', label='Tree Based ({} Bins; Uniform)'.format(tree_binning.n_bins)) diff --git a/examples/07_fact_pulls.py b/examples/07_fact_pulls.py index d38c7cd..611aa8f 100644 --- a/examples/07_fact_pulls.py +++ b/examples/07_fact_pulls.py @@ -60,7 +60,6 @@ def do_single_pull(obs_array_binning, llh = solution.StandardLLH(tau=tau, log_f=True, - vec_acceptance=vec_acceptance, C='thikonov') llh.initialize(vec_g=vec_g, model=tree_model) @@ -83,7 +82,7 @@ def do_single_pull(obs_array_binning, sol_mini = solution.LLHSolutionMinimizer() sol_mini.initialize(llh=llh, model=tree_model) sol_mini.set_x0_and_bounds(x0=x[idx_best]) - best_fit, _ = sol_mini.fit(constrain_N=False)[0] + best_fit = sol_mini.fit(constrain_N=False)[0] vec_f_str = ', '.join('{0:.2f}'.format(a) for a in best_fit.x) @@ -97,7 +96,7 @@ def do_single_pull(obs_array_binning, random_state=random_state) sol_mcmc.initialize(llh=llh, model=tree_model) sol_mcmc.set_x0_and_bounds(x0=best_fit.x) - vec_f_est_mcmc, sigma_vec_f, sample, probs = sol_mcmc.fit() + vec_f_est_mcmc, sigma_vec_f, sample, probs, autocorrelation = sol_mcmc.fit() vec_f_str = ', '.join('{0:.2f}'.format(a) for a in vec_f_est_mcmc) From 1292eab160141e178f07ce7fb3be884c586abe0f Mon Sep 17 00:00:00 2001 From: sudojan Date: Wed, 1 Jul 2020 10:25:23 +0200 Subject: [PATCH 5/7] fix bug in 9526156 --- funfolding/solution/likelihood.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/funfolding/solution/likelihood.py b/funfolding/solution/likelihood.py index 8a27220..9f2a75c 100644 --- a/funfolding/solution/likelihood.py +++ b/funfolding/solution/likelihood.py @@ -232,7 +232,8 @@ def evaluate_gradient(self, f): denom_f = f_reg + self.log_f_offset nom_f = np.log(denom_f * self.reg_factor_f) ln_10_squared = np.log(10)**2 - pre = self._C / ln_10_squared * np.broadcast_to(nom_f / denom_f, np.shape(_C)).T + pre = self._C / ln_10_squared * np.broadcast_to( + nom_f / denom_f, np.shape(self._C)).T reg_part = np.sum(pre, axis=0) + np.sum(pre, axis=1) else: reg_part = np.dot(self._C, f_reg * self.reg_factor_f) From cd57f434290c1893027d1e10819294ad54d70d22 Mon Sep 17 00:00:00 2001 From: sudojan Date: Wed, 1 Jul 2020 11:23:26 +0200 Subject: [PATCH 6/7] futher for loop replacing --- funfolding/solution/likelihood.py | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/funfolding/solution/likelihood.py b/funfolding/solution/likelihood.py index 9f2a75c..7abe665 100644 --- a/funfolding/solution/likelihood.py +++ b/funfolding/solution/likelihood.py @@ -231,9 +231,7 @@ def evaluate_gradient(self, f): if self.log_f_reg: denom_f = f_reg + self.log_f_offset nom_f = np.log(denom_f * self.reg_factor_f) - ln_10_squared = np.log(10)**2 - pre = self._C / ln_10_squared * np.broadcast_to( - nom_f / denom_f, np.shape(self._C)).T + pre = (self._C.T * (nom_f / denom_f)).T / np.log(10)**2 reg_part = np.sum(pre, axis=0) + np.sum(pre, axis=1) else: reg_part = np.dot(self._C, f_reg * self.reg_factor_f) @@ -378,18 +376,15 @@ def evaluate_gradient(self, f): def evaluate_hessian(self, f): m, n = self.linear_model.A.shape hess = np.zeros((n, n)) - for k in range(n): - for l in range(n): + denominator = np.dot(self.linear_model.A, f) + for i in range(n): + for j in range(n): poisson_part = 0 - for i in range(m): - A_ik = self.linear_model.A[i, k] - A_il = self.linear_model.A[i, l] - nominator = self.g[i] * A_ik * A_il - denominator = 0 - for j in range(n): - denominator += self.linear_model.A[i, j] * f[j] - poisson_part += nominator / denominator**2 - hess[k, l] = poisson_part + self.tau * self.C[k, l] + A_i = self.linear_model.A[:, i] + A_j = self.linear_model.A[:, j] + nominator = self.g * A_i * A_j + poisson_part = np.sum(nominator / denominator**2) + hess[i, j] = poisson_part + self.tau * self.C[i, j] return hess From 18df8c1aa181efa9efb585e7e545fdf0bd935aa9 Mon Sep 17 00:00:00 2001 From: sudojan Date: Wed, 1 Jul 2020 12:17:40 +0200 Subject: [PATCH 7/7] further loop removal --- funfolding/solution/likelihood.py | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/funfolding/solution/likelihood.py b/funfolding/solution/likelihood.py index 7abe665..6da1bf2 100644 --- a/funfolding/solution/likelihood.py +++ b/funfolding/solution/likelihood.py @@ -341,34 +341,20 @@ def __init__(self, g, linear_model, tau): def evaluate_llh(self, f): m, n = self.linear_model.A.shape - poisson_part = 0 - for i in range(m): - g_est = 0 - for j in range(n): - g_est += self.linear_model.A[i, j] * f[j] - poisson_part += g_est - self.g[i] * np.log(g_est) - - reg_part = 0 - for i in range(n): - for j in range(n): - reg_part += self.C[i, j] * f[i] * f[j] - reg_part *= 0.5 * self.tau + g_est = np.dot(self.linear_model.A * f) + poisson_part = np.sum(g_est - self.g * np.log(g_est)) + reg_part = np.dot(np.dot(self.C, f), f) * 0.5 * self.tau return reg_part - poisson_part def evaluate_gradient(self, f): m, n = self.linear_model.A.shape gradient = np.zeros(n) + g_est = np.dot(self.linear_model.A, f) for k in range(n): - poisson_part = 0 - for i in range(m): - g_est = 0 - for j in range(n): - g_est += self.linear_model.A[i, j] * f[j] - A_ik = self.linear_model.A[i, k] - poisson_part += A_ik - (self.g[i] * A_ik) / g_est - c = 0 - for i in range(n): - c += self.C[i, k] * f[i] + A_k = self.linear_model.A[:, k] + poisson_part = np.sum(A_k - (self.g * A_k) / g_est) + + c = np.dot(self.C[:, k], f) reg_part = self.tau * c gradient[k] = reg_part - poisson_part return gradient