tension_explainer.py

#
# Copyright(c) Daniel Knüttel 2022
#

#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.

from scipy.stats import norm
import numpy as np

class TensionExplainer(object):
    """
    Class that can explain tensions in terms of the standard error $\sigma$.
    Can be formatted in Jupyter Notebooks.

    GNU GPL'ed v3.
    """
    def __init__(self, tension):
        self._tension = tension
        self._age_of_universe = 436_117_077_000_000_000
        self._planck_time = 5.391247e-44
        self._seconds_per_year = 31557600
        self._seconds_per_month = 2629800
        self._seconds_per_day = 86400
        self._seconds_per_hour = 3600

    def scientific_format(self, flt):
        exponent = np.round(np.log10(flt))
        base = flt / 10**exponent
        return r"$%.3f \times 10^{%d}$" % (base, exponent)

    def explain_tension(self):
        probability = norm.cdf(np.abs(self._tension)) - norm.cdf(-np.abs(self._tension))
        if probability > .2:
            return self.high_probability_explanation(probability)
        if probability > .05:
            return self.low_probability_explanation(probability)
        return self.very_low_probability_explanation(probability)

    def high_probability_explanation(self, probability):
        return fr"Probability is ${probability * 100:.0f} \%$. This is a likely event."

    def low_probability_explanation(self, probability):
        return fr"Probability is ${probability * 100:.1f} \%$. This is an unlikely but possible event. Try to re-run the experiment."

    def very_low_probability_explanation(self, probability):
        if(np.abs(self._tension) > 31):
            return self.impossibly_low_probability_explanation()

        wait_time_seconds = 1 / probability 
        if(wait_time_seconds < self._seconds_per_hour):
            return (f"The expected wait time for an event like this is ${wait_time_seconds:.0f}$"
                    " seconds at one possible event per second.")
        if(wait_time_seconds < self._seconds_per_day):
            return (f"The expected wait time for an event like this is ${wait_time_seconds / self._seconds_per_hour:.0f}$"
                    " hours at one possible event per second.")
        if(wait_time_seconds < self._seconds_per_month):
            return (f"The expected wait time for an event like this is ${wait_time_seconds / self._seconds_per_day:.0f}$"
                    " days at one possible event per second.")
        if(wait_time_seconds < self._seconds_per_year):
            return (f"The expected wait time for an event like this is ${wait_time_seconds / self._seconds_per_month:.0f}$"
                    " months at one possible event per second.")
        if(wait_time_seconds < self._age_of_universe):
            return (f"The expected wait time for an event like this is {self.scientific_format(wait_time_seconds / self._seconds_per_year)}"
                    " years at one possible event per second.")

        if(wait_time_seconds < self._age_of_universe / self._planck_time):
            return (f"The expected wait time for an event like this is {self.scientific_format(wait_time_seconds / self._age_of_universe)}"
                    " ages of the universe at one possible event per second.")

        return (f"The expected wait time for an event like this is {self.scientific_format(wait_time_seconds / self._age_of_universe * self._planck_time)}"
                " ages of the universe at one possible event per Planck time (!).")


    def impossibly_low_probability_explanation(self):
        aou_planck_times = self._age_of_universe / self._planck_time
        tension_31_prob = aou_planck_times * norm.cdf(-31)
        return (fr"The tension is ${self._tension:.0f}\sigma$. This is virtually impossible, for comparison at $31\sigma$"
                f" tension one would have to wait {self.scientific_format(1/tension_31_prob)} ages of the universe"
                " assuming one possible event per Planck time.")

    def _repr_latex_(self):
        return self.explain_tension()