5.4.2.5 adds 3, 4, and 5 sigma calc to NIE in branch

Signed-off-by: boban <bob.anderson.ok@gmail.com>

src/pyoteapp/false_positive.py

@@ -3,6 +3,7 @@
 from numba import prange, njit, int64, float64
 # import time
 import matplotlib
+from scipy.integrate import quad
 # import matplotlib.pyplot as plt
 # from pyoteapp.autocorrtools import *
 matplotlib.use('Qt5Agg')
@@ -43,11 +44,9 @@ def simple_convolve2(x, c):
         i += 1
     return out
 
-
 def noise_gen_numpy(num_points: int, noise_sigma: float) -> np.ndarray:
     return random.normal(size=num_points, scale=noise_sigma)
 
-
 # @jit(float64[:](int64, float64, float64[:]))
 def simulated_observation(observation_size: int, noise_sigma: float, corr_array: np.ndarray) -> np.ndarray:
     """Returns a numpy array of floats representing an observation that contains only correlated noise.
@@ -124,6 +123,20 @@ def max_drop_in_simulated_observation(
 
     return best_drop_so_far
 
+# Locate the index of the righthand bin edge of the bin that includes the value drop
+def index_of_drop(drop, bin_edges):
+    for i in range(1, bin_edges.size):
+        if bin_edges[i-1] < drop <= bin_edges[i]:
+            return i
+    return None
+
+def tail_count_fraction(x, slope, y0):
+    arg = np.polyval((slope, y0), x)
+    exp1 = np.exp(arg)
+    return np.exp(-exp1)
+
+def tail_area(x0, slope, y0):
+    return quad(lambda x: tail_count_fraction(x,slope,y0), x0, np.inf)[0]
 
 @njit(parallel=True)
 def compute_drops(event_duration: int,
@@ -150,6 +163,87 @@ def compute_drops(event_duration: int,
     return drops
 
 
+def calc_sigma_lines(three_sigma_guess, slope, y0, bin_delta, debug=False):
+    five_sig = 0.999_999_427
+    four_sig = 0.999_936_658
+    three_sig = 0.997_300_204
+    two_sig = 0.954_499_736
+    # one_sig = 0.682_689_492
+
+    three_sigma_line = three_sigma_guess
+    area_fraction = tail_area(three_sigma_guess, slope, y0) / bin_delta
+    p = 1 - area_fraction
+
+    if p < three_sig:
+        delta = 0.05
+        while True:
+            area_fraction = tail_area(three_sigma_line, slope, y0) / bin_delta
+            p = 1 - area_fraction
+            if p < three_sig:
+                three_sigma_line += delta
+            else:
+                three_sig_area = area_fraction
+                break
+    else:
+        delta = -0.05
+        while True:
+            area_fraction = tail_area(three_sigma_line, slope, y0) / bin_delta
+            p = 1 - area_fraction
+            if p > three_sig:
+                three_sigma_line += delta
+            else:
+                three_sig_area = area_fraction
+                break
+
+    five_sigma_line = three_sigma_line
+    delta = 0.05
+    while True:
+        area_fraction = tail_area(five_sigma_line, slope, y0) / bin_delta
+        p = 1 - area_fraction
+        if p < five_sig:
+            five_sigma_line += delta
+        else:
+            five_sig_area = area_fraction
+            break
+
+    four_sigma_line = three_sigma_line
+    delta = 0.05
+    while True:
+        area_fraction = tail_area(four_sigma_line, slope, y0) / bin_delta
+        p = 1 - area_fraction
+        if p < four_sig:
+            four_sigma_line += delta
+        else:
+            four_sig_area = area_fraction
+            break
+
+    two_sigma_line = three_sigma_line
+    delta = -0.05
+    while True:
+        area_fraction = tail_area(two_sigma_line, slope, y0) / bin_delta
+        p = 1 - area_fraction
+        if p > two_sig:
+            two_sigma_line += delta
+        else:
+            two_sig_area = area_fraction
+            break
+
+    # one_sigma_line = sorted_drops[int(.682689 * drops.size)]
+
+    if debug:
+        # print(f'  one_sigma_line probability: 0.682689492')
+        print(f'  two_sigma_line probability: {1 - two_sig_area:0.9f}')
+        print(f'three_sigma_line probability: {1 - three_sig_area:0.9f}')
+        print(f' four_sigma_line probability: {1 - four_sig_area:0.9f}')
+        print(f' five_sigma_line probability: {1 - five_sig_area:0.9f}')
+        # print(f'\n  one_sigma_line: {one_sigma_line:0.1f}')
+        print(f'  two_sigma_line: {two_sigma_line:0.1f}')
+        print(f'three_sigma_line: {three_sigma_line:0.1f}  three_sigma_guess: {three_sigma_guess:0.1f}')
+        print(f' four_sigma_line: {four_sigma_line:0.1f}')
+        print(f' five_sigma_line: {five_sigma_line:0.1f}')
+
+    return two_sigma_line, three_sigma_line, four_sigma_line, five_sigma_line
+
 # def excercise():
 #     noise_gen_jit(1, 1.0)
 #     noise_gen_jit(1, 1.0)

src/pyoteapp/pyote.py

@@ -64,7 +64,7 @@
 from pyoteapp.showVideoFrames import readFitsFile
 from pyoteapp.showVideoFrames import readAavFile
 
-from pyoteapp.false_positive import compute_drops, noise_gen_jit, simple_convolve
+from pyoteapp.false_positive import compute_drops, noise_gen_jit, simple_convolve, calc_sigma_lines, tail_count_fraction
 import numpy as np
 import pyqtgraph as pg
 import pyqtgraph.exporters as pex
@@ -2374,7 +2374,7 @@ def demoGeometry(self):
                 fig = plt.figure(constrained_layout=False, figsize=(6, 6))
                 fig.canvas.manager.set_window_title('Disk on disk geometry')
                 axes = fig.add_subplot(1, 1, 1)
-                illustrateDiskOnDiskEvent(LCP=self.Lcp, showLegend=True, showNotes=True, axes=axes)
+                illustrateEdgeOnDiskEvent(LCP=self.Lcp, showLegend=True, showNotes=True, axes=axes)
                 plt.show()
             elif self.diffractionRadioButton.isChecked():
                 majorAxis = self.Lcp.asteroid_major_axis
@@ -8479,33 +8479,72 @@ def falsePositiveReport(event_duration, num_trials, observation_size, observed_d
                               noise_sigma=sigma, corr_array=np.array(posCoefs), num_trials=num_trials)
         sorted_drops = np.sort(drops)
         max_drop = sorted_drops[-1]
-        three_sigma_line = sorted_drops[int(.997 * drops.size)]
+        # three_sigma_line = sorted_drops[int(.997 * drops.size)]
+
+        # TODO Add new sigma line call
+
+        counts, bins = np.histogram(drops, bins='auto')
+
+        counts[0] = counts[1] / 2
+
+        # Calculate the tail start drop as the drop that encloses tail_start_fraction of the counts
+        tail_start_fraction = 0.85
+        tail_start_index = np.where(np.cumsum(counts) > tail_start_fraction * num_trials)[0][0]
+
+        # Find all of the non-zero counts and their locations in the tail
+        debug = False
+        x_tail = []
+        y = []
+        zero_count = 0
+        bin_delta = bins[1] - bins[0]
+        if debug: print(f'bin_delta: {bin_delta:0.3f}')
+        for i in range(tail_start_index, counts.size):
+            if counts[i] > 0:
+                new_x = bins[i] + bin_delta / 2
+                x_tail.append(new_x)
+                y.append(counts[i] / num_trials)
+                if debug: print(f'x_tail: {new_x:6.3f}  count: {counts[i]:0.6f}')
+            else:
+                zero_count += 1
+                if debug: print(f'zero')
+
+        loglogy = np.log(-np.log(y))
+        slope, y0, *_ = np.polyfit(x_tail, loglogy, 1)
+
+        three_sigma_estimate = sorted_drops[int(.997 * drops.size)]
+
+        two_sigma_line, three_sigma_line, four_sigma_line, five_sigma_line = \
+            calc_sigma_lines(three_sigma_guess=three_sigma_estimate, slope=slope, y0=y0, bin_delta=bin_delta,
+                             debug=debug)
 
         if plots_wanted:
             pw = PlotWidget(viewBox=CustomViewBox(border=(0, 0, 0)),
                             enableMenu=False,
                             title=f'Noise Induced Events (brightness drops) found in correlated noise for event duration: {event_duration}',
                             labels={'bottom': 'brightness drop', 'left': 'number of times noise produced drop'})
             pw.hideButtons()
-            y, x = np.histogram(drops, bins='auto')
+            # y, x = np.histogram(drops, bins='auto')
 
-
-
-            y[0] = y[1] / 2.0
+            # y[0] = y[1] / 2.0
             # Plot drops histogram
-            pw.plot(x, y, stepMode=True, fillLevel=0, brush=(0, 0, 255, 150))
+            pw.plot(bins, counts, stepMode=True, fillLevel=0, brush=(0, 0, 255, 150))
 
             # Plot red observed drop line
-            pw.plot(x=[observed_drop, observed_drop], y=[0, 1.5 * np.max(y)], pen=pg.mkPen([255, 0, 0], width=2))
+            pw.plot(x=[observed_drop, observed_drop], y=[0, 1.5 * np.max(counts)], pen=pg.mkPen([255, 0, 0], width=2))
 
             # Plot black max drop line
-            pw.plot(x=[max_drop, max_drop], y=[0, 0.1 * np.max(y)], pen=pg.mkPen([0, 0, 0], width=2))
+            # pw.plot(x=[max_drop, max_drop], y=[0, 0.1 * np.max(counts)], pen=pg.mkPen([0, 0, 0], width=2))
 
             # Plot green three_sigma_line
-            pw.plot(x=[three_sigma_line, three_sigma_line], y=[0, 0.5 * np.max(y)], pen=pg.mkPen([0, 255, 0], width=3))
+            pw.plot(x=[three_sigma_line, three_sigma_line], y=[0, 0.5 * np.max(counts)], pen=pg.mkPen([0, 255, 0], width=3))
+            pw.plot(x=[four_sigma_line, four_sigma_line], y=[0, 0.3 * np.max(counts)], pen=pg.mkPen([0, 255, 0], width=3))
+            pw.plot(x=[five_sigma_line, five_sigma_line], y=[0, 0.2 * np.max(counts)], pen=pg.mkPen([0, 255, 0], width=3))
+
+            pw.plot(x_tail, tail_count_fraction(x_tail, slope, y0) * num_trials, pen=pg.mkPen([0, 0, 0], width=3),
+                     label='fit to tail')
 
             # Plot a nice 0 line
-            right_edge = max(observed_drop * 1.1, np.max(x) * 1.1)
+            right_edge = max(observed_drop * 1.1, np.max(bins) * 1.1)
             pw.plot(x=[0, right_edge], y=[0, 0], pen=pg.mkPen([180, 180, 180], width=2))
 
             pw.addLegend()

src/pyoteapp/version.py

@@ -1,2 +1,2 @@
 def version():
-    return '5.4.2.4'
+    return '5.4.2.5'