Example03

.gitignore

@@ -89,3 +89,5 @@ ENV/
 
 # Rope project settings
 .ropeproject
+
+.DS_Store

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)

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))
 

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)

funfolding/solution/likelihood.py

@@ -229,19 +229,10 @@ 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.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)
         else:
@@ -350,53 +341,36 @@ 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
 
     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
 
 

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