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        import theano.tensor as tt
23        from theano import as_op
24        from scipy.stats import gaussian_kde
25
26        smin, smax = np.min(samples), np.max(samples)
27        width = smax - smin
28        x = np.linspace(smin, smax, 100)
29        y = gaussian_kde(samples)(x)
30
31        # what was never sampled should have a small probability but not 0,
32        # so we'll extend the domain and use linear approximation of density on it
33        x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
34        y = np.concatenate([[0], y, [0]])
35
36        return pm.distributions.Interpolated(param, x, y)
37
38    def error_model_from_trace(self, trace, ymincentroid):
39        """
40        # Legacy code for Bayesian efforts - not used.
41        pymc3 is not installed by default,
42            if have plans to use it manual installation of pymc3
43            package before using this method is needed
44        """
45        import pymc3 as pm
46        # from pymc3 import traceplot, plot_posterior
47
48        with pm.Model() as model2:
49            sd = self.from_posterior("sd", trace["sd"])
50            y = pm.HalfNormal("y", sd=sd, observed=ymincentroid)
51            start = pm.find_MAP()
52            step = pm.NUTS()  # Hamiltonian MCMC with No U-Turn Sampler
53            trace = pm.sample(
54                1000, step, start, random_seed=123, progressbar=True, tune=1000
55            )
56            pm.summary(trace)
57            # plot_posterior(trace)
58            # traceplot(trace)
59            return pm.summary(trace)["mean"].values[0]
60
61    def simple_model_error_dist(self, ymincentroid):
62        """
63        # Legacy code for Bayesian efforts - not used.
64        pymc3 is not installed by default,
65            if have plans to use it manual installation of pymc3
66            package before using this method is needed
67        """
68        import pymc3 as pm
69        # from pymc3 import traceplot, plot_posterior
70        # import seaborn as sns
71        # f, ax = pyplot.subplots(figsize=(6, 6))
72        # sns.distplot(ymincentroid)
73        # sns.kdeplot(ymincentroid, ax=ax, shade=True, color="g")
74        # sns.rugplot(ymincentroid, color="black", ax=ax)
75        # ax.set(xlabel= "Peak Minima Magnitude", ylabel= "Density")
76        # pyplot.show()
77
78        with pm.Model() as model:
79            # mu = pm.Uniform('mu', lower=-1, upper=1)
80            lower = ymincentroid.min()
81            upper = ymincentroid.max()
82
83            sd = pm.Uniform("sd", lower=lower, upper=upper)
84
85            y = pm.HalfNormal("y", sd=sd, observed=ymincentroid)
86
87            start = pm.find_MAP()
88            step = pm.NUTS()  # Hamiltonian MCMC with No U-Turn Sampler
89            trace = pm.sample(
90                1000, step, start, random_seed=123, progressbar=True, tune=1000
91            )
92
93            return pm.summary(trace)["mean"].values[0]