1"""
 2This code is for Bayesian estimation of the noise levels.
 3It is it not implemented or used in the current code base.
 4The packages it uses are not part of the requirements.
 5If you want to use it, you will need to install them manually.
 6"""
 7
 8from corems.mass_spectrum.calc.NoiseCalc import NoiseThresholdCalc
 9
10
11class BayesNoiseCalc(NoiseThresholdCalc):
12    def from_posterior(self, param, samples):
13        """
14        # Legacy code for Bayesian efforts - not used.
15        pymc3 is not installed by default,
16            if have plans to use it manual installation of pymc3
17            package before using this method is needed
18        """
19
20        import pymc3 as pm
21        import numpy as np
22        from scipy.stats import gaussian_kde
23
24        smin, smax = np.min(samples), np.max(samples)
25        width = smax - smin
26        x = np.linspace(smin, smax, 100)
27        y = gaussian_kde(samples)(x)
28
29        # what was never sampled should have a small probability but not 0,
30        # so we'll extend the domain and use linear approximation of density on it
31        x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
32        y = np.concatenate([[0], y, [0]])
33
34        return pm.distributions.Interpolated(param, x, y)
35
36    def error_model_from_trace(self, trace, ymincentroid):
37        """
38        # Legacy code for Bayesian efforts - not used.
39        pymc3 is not installed by default,
40            if have plans to use it manual installation of pymc3
41            package before using this method is needed
42        """
43        import pymc3 as pm
44        # from pymc3 import traceplot, plot_posterior
45
46        with pm.Model() as model2:
47            sd = self.from_posterior("sd", trace["sd"])
48            y = pm.HalfNormal("y", sd=sd, observed=ymincentroid)
49            start = pm.find_MAP()
50            step = pm.NUTS()  # Hamiltonian MCMC with No U-Turn Sampler
51            trace = pm.sample(
52                1000, step, start, random_seed=123, progressbar=True, tune=1000
53            )
54            pm.summary(trace)
55            # plot_posterior(trace)
56            # traceplot(trace)
57            return pm.summary(trace)["mean"].values[0]
58
59    def simple_model_error_dist(self, ymincentroid):
60        """
61        # Legacy code for Bayesian efforts - not used.
62        pymc3 is not installed by default,
63            if have plans to use it manual installation of pymc3
64            package before using this method is needed
65        """
66        import pymc3 as pm
67        # from pymc3 import traceplot, plot_posterior
68        # import seaborn as sns
69        # f, ax = pyplot.subplots(figsize=(6, 6))
70        # sns.distplot(ymincentroid)
71        # sns.kdeplot(ymincentroid, ax=ax, shade=True, color="g")
72        # sns.rugplot(ymincentroid, color="black", ax=ax)
73        # ax.set(xlabel= "Peak Minima Magnitude", ylabel= "Density")
74        # pyplot.show()
75
76        with pm.Model() as model:
77            # mu = pm.Uniform('mu', lower=-1, upper=1)
78            lower = ymincentroid.min()
79            upper = ymincentroid.max()
80
81            sd = pm.Uniform("sd", lower=lower, upper=upper)
82
83            y = pm.HalfNormal("y", sd=sd, observed=ymincentroid)
84
85            start = pm.find_MAP()
86            step = pm.NUTS()  # Hamiltonian MCMC with No U-Turn Sampler
87            trace = pm.sample(
88                1000, step, start, random_seed=123, progressbar=True, tune=1000
89            )
90
91            return pm.summary(trace)["mean"].values[0]