diff --git a/src/pyoteapp/false_positive.py b/src/pyoteapp/false_positive.py index 24f7bb3..c1c135d 100644 --- a/src/pyoteapp/false_positive.py +++ b/src/pyoteapp/false_positive.py @@ -162,8 +162,7 @@ def compute_drops(event_duration: int, return drops - -def calc_sigma_lines(three_sigma_guess, slope, y0, bin_delta, debug=False): +def calc_sigma_lines(observed_drop, 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 @@ -172,7 +171,13 @@ def calc_sigma_lines(three_sigma_guess, slope, y0, bin_delta, debug=False): three_sigma_line = three_sigma_guess area_fraction = tail_area(three_sigma_guess, slope, y0) / bin_delta - p = 1 - area_fraction + if observed_drop > three_sigma_line: + drop_fraction = tail_area(observed_drop, slope, y0) / bin_delta + drop_nie_probability = drop_fraction + else: + drop_nie_probability = 1.0 + # print(f'drop_nie_probability: {drop_nie_probability:0.6f}') + p = 1.0 - area_fraction if p < three_sig: delta = 0.05 @@ -242,7 +247,7 @@ def calc_sigma_lines(three_sigma_guess, slope, y0, bin_delta, debug=False): 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 + return two_sigma_line, three_sigma_line, four_sigma_line, five_sigma_line, drop_nie_probability # def excercise(): # noise_gen_jit(1, 1.0) diff --git a/src/pyoteapp/pyote.py b/src/pyoteapp/pyote.py index 5cf9f03..1158a25 100644 --- a/src/pyoteapp/pyote.py +++ b/src/pyoteapp/pyote.py @@ -6416,6 +6416,16 @@ def showQuery(self, question, title=''): msgBox.setDefaultButton(QMessageBox.Yes) self.queryRetVal = msgBox.exec_() + def showInfoNonModal(self, stuffToSay, title=''): + msgBox = QMessageBox(self) + msgBox.setIcon(QMessageBox.Information) + msgBox.setText(stuffToSay) + msgBox.setWindowTitle(title) + msgBox.setStandardButtons(QMessageBox.Yes | QMessageBox.No) + msgBox.setDefaultButton(QMessageBox.Yes) + # self.queryRetVal = msgBox.show() + return msgBox + def fillPrimaryAndRef(self): # Load self.yValues and sel.yRefStar with proper light curves as # indicated by the spinner values @@ -7430,6 +7440,9 @@ def finalReportPenumbral(self): f'{self.observedDrop:0.1f} three sigma drop from noise: {self.threeSigmaLine:0.1f} margin: {margin:0.1f}' self.showMsg(stats_msg, blankLine=False, bold=True, color='red') + stats_msg = f'fit metrics === observed drop has probability {self.drop_nie_probability:0.6f} of being a noise-induced-event' + self.showMsg(stats_msg, blankLine=False, bold=True) + stats_msg = f'fit metrics === {self.magDropReportStr}' self.showMsg(stats_msg, blankLine=False, bold=True) @@ -7454,7 +7467,7 @@ def finalReportPenumbral(self): self.updateFitMetricTxtFile() - def finalReport(self, false_positive, false_probability, error_bar_plots_available=True): + def finalReport(self, error_bar_plots_available=True): self.xlsxDict = {} self.writeDefaultGraphicsPlots(error_bar_plots_available=error_bar_plots_available) @@ -7511,19 +7524,24 @@ def finalReport(self, false_positive, false_probability, error_bar_plots_availab self.reportSpecialProcedureUsed() # This includes use of asteroid distance/speed and star diameter - if false_positive: - self.showMsg(f"This 'drop' has a {false_probability:0.5f} probability of being an artifact of noise.", - bold=True, color='red', blankLine=False) - else: - self.showMsg(f"This 'drop' has a zero probability of being an artifact of noise.", - bold=True, color='green', blankLine=False) - - self.showMsg(f">>>> probability > 0.0000 indicates the 'drop' may be spurious (a noise artifact)." - f" Consult with an IOTA Regional Coordinator.", color='blue', blankLine=False) - self.showMsg(f">>>> probability = 0.0000 indicates the 'drop' is unlikely to be a noise artifact, but" - f" does not prove that the 'drop' is due to an occultation", color='blue', blankLine=False) - self.showMsg(f">>>> Consider 'drop' shape, timing, mag drop, duration and other positive observer" - f" chords before reporting the 'drop' as a positive.", color='blue') + # TODO New nie report + # if false_positive: + # self.showMsg(f"This 'drop' has a {false_probability:0.5f} probability of being an artifact of noise.", + # bold=True, color='red', blankLine=False) + # else: + # self.showMsg(f"This 'drop' has a zero probability of being an artifact of noise.", + # bold=True, color='green', blankLine=False) + + self.showMsg(f'The observed drop has a probability of {self.drop_nie_probability:0.6f} of being a noise-induced-event', + bold=True, color='blue', blankLine=True) + + # self.showMsg(f">>>> probability > 0.0000 indicates the 'drop' may be spurious (a noise artifact)." + # f" Consult with an IOTA Regional Coordinator.", color='blue', blankLine=False) + # self.showMsg(f">>>> probability = 0.0000 indicates the 'drop' is unlikely to be a noise artifact, but" + # f" does not prove that the 'drop' is due to an occultation", color='blue', blankLine=False) + # self.showMsg(f">>>> Consider 'drop' shape, timing, mag drop, duration and other positive observer" + # f" chords before reporting the 'drop' as a positive.", color='blue') + # TODO End new nie report self.showMsg("All times are calculated/reported based on the assumption that timestamps are " "start-of-exposure times.", @@ -7600,6 +7618,9 @@ def finalReport(self, false_positive, false_probability, error_bar_plots_availab f'{self.observedDrop:0.1f} three sigma drop from noise: {self.threeSigmaLine:0.1f} margin: {margin:0.1f}' self.showMsg(stats_msg, blankLine=False, bold=True, color='red') + stats_msg = f'fit metrics === observed drop has probability {self.drop_nie_probability:0.6f} of being a noise-induced-event' + self.showMsg(stats_msg, blankLine=False, bold=True) + stats_msg = f'fit metrics === {self.magDropReportStr}' self.showMsg(stats_msg, blankLine=False, bold=True) @@ -8034,8 +8055,11 @@ def computeErrorBars(self, plots_wanted=True): if self.errBarWin is not None: self.errBarWin.close() - _, false_positive, false_probability, self.observedDrop, self.threeSigmaLine = \ - self.doFalsePositiveReport(posCoefs, plots_wanted=False) # noqa + # _, false_positive, false_probability, self.observedDrop, self.threeSigmaLine, self.fourSigmaLine, \ + # self.fiveSigmaLine = self.doFalsePositiveReport(posCoefs, plots_wanted=False) # noqa + + fp_plot, false_positive, false_probability, self.observedDrop, self.threeSigmaLine, \ + self.fourSigmaLine, fiveSigmaLine = self.doFalsePositiveReport(posCoefs) # noqa if plots_wanted: self.errBarWin = pg.GraphicsWindow( @@ -8056,7 +8080,9 @@ def computeErrorBars(self, plots_wanted=True): self.durBarPlotItem = pw2.getPlotItem() pw2.hideButtons() - pw3, false_positive, false_probability, self.observedDrop, self.threeSigmaLine = self.doFalsePositiveReport(posCoefs) # noqa + pw3 = fp_plot + # pw3, false_positive, false_probability, self.observedDrop, self.threeSigmaLine, \ + # self.fourSigmaLine, fiveSigmaLine = self.doFalsePositiveReport(posCoefs) # noqa # Suppress this print statement # print(f'observed drop: {self.observedDrop:0.2f} max noise-induced drop: {self.maxNoiseInducedDrop:0.2f}') @@ -8174,7 +8200,8 @@ def computeErrorBars(self, plots_wanted=True): if self.timestampListIsEmpty(self.yTimes): self.showMsg('Cannot produce final report because timestamps are missing.', bold=True, color='red') else: - self.finalReport(false_positive, false_probability, error_bar_plots_available=plots_wanted) + # self.finalReport(false_positive, false_probability, error_bar_plots_available=plots_wanted) + self.finalReport(error_bar_plots_available=plots_wanted) self.fillExcelReportButton.setEnabled(True) self.newRedrawMainPlot() # To add envelope to solution @@ -8473,12 +8500,22 @@ def doFalsePositiveReport(self, posCoefs, plots_wanted=True): return self.falsePositiveReport(event_duration, num_trials, observation_size, observed_drop, posCoefs, sigma, plots_wanted=plots_wanted) - @staticmethod - def falsePositiveReport(event_duration, num_trials, observation_size, observed_drop, posCoefs, sigma, plots_wanted=True): + def falsePositiveReport(self, event_duration, num_trials, observation_size, observed_drop, posCoefs, sigma, plots_wanted=True): + # print('we are computing error bars') + # QtWidgets.QApplication.processEvents() + # msgBox = self.showInfoNonModal('Computing error bars') + # msgBox.show() + drops = compute_drops(event_duration=event_duration, observation_size=observation_size, noise_sigma=sigma, corr_array=np.array(posCoefs), num_trials=num_trials) + + # QtWidgets.QApplication.processEvents() + # msgBox.close() + # del msgBox + # QtWidgets.QApplication.processEvents() + sorted_drops = np.sort(drops) - max_drop = sorted_drops[-1] + # max_drop = sorted_drops[-1] # three_sigma_line = sorted_drops[int(.997 * drops.size)] # TODO Add new sigma line call @@ -8513,10 +8550,17 @@ def falsePositiveReport(event_duration, num_trials, observation_size, observed_d 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, + two_sigma_line, three_sigma_line, four_sigma_line, five_sigma_line, self.drop_nie_probability = \ + calc_sigma_lines(observed_drop=observed_drop, three_sigma_guess=three_sigma_estimate, + slope=slope, y0=y0, bin_delta=bin_delta, debug=debug) + if self.drop_nie_probability == 1.0: + indices = np.where(sorted_drops > observed_drop) + if indices: + index_at_observed_drop = indices[0][0] + self.drop_nie_probability = index_at_observed_drop / num_trials + if plots_wanted: pw = PlotWidget(viewBox=CustomViewBox(border=(0, 0, 0)), enableMenu=False, @@ -8548,11 +8592,12 @@ def falsePositiveReport(event_duration, num_trials, observation_size, observed_d pw.plot(x=[0, right_edge], y=[0, 0], pen=pg.mkPen([180, 180, 180], width=2)) pw.addLegend() - pw.plot(name='red line: the brightness drop (B - A) extracted from lightcurve') - pw.plot(name=f'black line: max brightness drop found in {num_trials} trials against pure noise') - pw.plot(name=f'green three_sigma_line: 99.7% of noise induced brightness drops are to the left.') + pw.plot(name=f'red line: the observed drop (B - A) extracted from lightcurve has ' + f'probability {self.drop_nie_probability:0.6f} of being noise induced') + pw.plot(name=f'black curve is LSQ fit to drop data using the Gumbel Extreme Value distribution') + pw.plot(name=f'green lines are at 3 sigma (99.7%), 4 sigma (99.9938%), and 5 sigma (99.99994%)') pw.plot(name=f'PyOTE reports a "pass" if red line is to the right of the three_sigma_line,') - pw.plot(name=f'BUT !!! If red line is left of black line but right of green line - refer to regional coordinator.') + pw.plot(name=f'BUT !!! If red line is left of the 5 sigma line - refer to regional coordinator.') else: pw = None y, x = np.histogram(drops, bins=50) @@ -8579,7 +8624,7 @@ def falsePositiveReport(event_duration, num_trials, observation_size, observed_d # TODO Changed criteria for noise-induced-event fail # return pw, false_positive, false_probability, observed_drop, np.max(x) - return pw, false_positive, false_probability, observed_drop, three_sigma_line + return pw, false_positive, false_probability, observed_drop, three_sigma_line, four_sigma_line, five_sigma_line def displaySolution(self, subframe=True): D, R = self.solution @@ -9100,7 +9145,15 @@ def findEvent(self, plots_wanted=False): self.newRedrawMainPlot() if not self.suppressReport: + # print('we are computing error bars') + QtWidgets.QApplication.processEvents() + # msgBox = self.showInfoNonModal('Computing error bars') + # msgBox.show() self.computeErrorBars(plots_wanted=plots_wanted) + QtWidgets.QApplication.processEvents() + # msgBox.close() + QtWidgets.QApplication.processEvents() + self.suppressReport = False if need_to_invite_user_to_verify_timestamps: diff --git a/src/pyoteapp/version.py b/src/pyoteapp/version.py index f5d9d2c..f1a11a7 100644 --- a/src/pyoteapp/version.py +++ b/src/pyoteapp/version.py @@ -1,2 +1,2 @@ def version(): - return '5.4.2.5' + return '5.4.2.6'