Make monte carlo runs identifiable

examples/MonteCarloExamples/scenarioAnalyzeMonteCarlo.py

@@ -20,7 +20,7 @@
 r"""
 Motivation
 ----------
-This script is a basic demonstration of a script that can be used to plot Monte Carlo data with 
+This script is a basic demonstration of a script that can be used to plot Monte Carlo data with
 bokeh and datashaders.   These tools are very efficient to plot large amounts of simulation data
 that is likely to occur with Monte Carlo sensitivity analysis studies.  For example, running this script will
 create an HTML interactive view of the simulation data.   Instead of seeing a fixed resolution, the user can
@@ -80,7 +80,7 @@
 
     Read all three steps before advancing.
 
-The next steps outline how to run this script. 
+The next steps outline how to run this script.
 
 1.  This script can only be run once there exists data produced by the ``scenario_AttFeedbackMC.py`` script.
 
@@ -102,26 +102,33 @@
 
 import inspect
 import os
-FOUND_DATESHADER = True
+import holoviews as hv
+
+
+FOUND_DATASHADER = True
 try:
-    from Basilisk.utilities.datashader_utilities import DS_Plot, curve_per_df_component, pull_and_format_df
+    from Basilisk.utilities.datashader_utilities import DS_Plot, curve_per_df_column, pull_and_format_df
     from Basilisk.utilities.MonteCarlo.AnalysisBaseClass import mcAnalysisBaseClass
     from bokeh.palettes import Blues9, Reds9, Greens9, \
         Blues3, Reds3, Greens3, Oranges3, RdYlBu9
+    from bokeh.server.server import Server
+    from bokeh.application import Application
+    from bokeh.application.handlers.function import FunctionHandler
+
 except:
     print("Wasn't able to include the datashader_utilities.")
-    FOUND_DATESHADER = False
+    FOUND_DATASHADER = False
 
 import Basilisk.utilities.macros as macros
 
-
 filename = inspect.getframeinfo(inspect.currentframe()).filename
 fileNameString = os.path.basename(os.path.splitext(__file__)[0])
 path = os.path.dirname(os.path.abspath(filename))
 from Basilisk import __path__
 
 bskPath = __path__[0]
 
+
 def plotSuite(dataDir):
     """
     This is the function to populate with all of the plots to be generated using datashaders and bokeh.
@@ -139,20 +146,23 @@ def plotSuite(dataDir):
     sigmaPlot = DS_Plot(sigma_BR, title="Attitude Error",
                         xAxisLabel='time [s]', yAxisLabel='Sigma_BR',
                         macro_x=macros.NANO2SEC,
-                        labels = ['b1', 'b2', 'b3'], cmap=RdYlBu9,
-                        plotFcn=curve_per_df_component)
+                        labels=['b1', 'b2', 'b3'], cmap=RdYlBu9,
+                        plotObjType=hv.Points,
+                        plotFcn=curve_per_df_column)
     plotList.append(sigmaPlot)
 
-    sigma_BR = pull_and_format_df(dataDir + "attGuidMsg.omega_BR_B.data", 3)
-    sigmaPlot = DS_Plot(sigma_BR, title="Attitude Rate Error",
-                        xAxisLabel='time [s]', yAxisLabel='omega_BR_B',
-                        macro_x=macros.NANO2SEC, macro_y=macros.R2D,
-                        labels = ['b1', 'b2', 'b3'], cmap=RdYlBu9,
-                        plotFcn=curve_per_df_component)
-    plotList.append(sigmaPlot)
+    # sigma_BR = pull_and_format_df(dataDir + "attGuidMsg.omega_BR_B.data", 3)
+    # sigmaPlot = DS_Plot(sigma_BR, title="Attitude Rate Error",
+    #                     xAxisLabel='time [s]', yAxisLabel='omega_BR_B',
+    #                     macro_x=macros.NANO2SEC, macro_y=macros.R2D,
+    #                     labels=['b1', 'b2', 'b3'], cmap=RdYlBu9,
+    #                     plotObjType=hv.Points,
+    #                     plotFcn=curve_per_df_column)
+    # plotList.append(sigmaPlot)
     return plotList
 
 
+
 def run(show_plots):
     """
     **This script is meant to be configured based on the user's needs. It can be configured using the following
@@ -161,7 +171,7 @@ def run(show_plots):
     First, set ``show_all_data = True`` to get a broad view of the data and find a time window to investigate closer.
 
     Once the data is characterized, the user can set ``show_extreme_data = True`` to look at specific run cases
-    within the window.
+    within the window.title=title, xAxisLabel=x_axi
 
     Finally, the user can set ``show_optional_data = True`` to look at any extra data to determine why the extrema
     cases exist.
@@ -171,16 +181,15 @@ def run(show_plots):
     :param optional_plots: plots additional user-defined plots
     """
 
-    if not FOUND_DATESHADER:
+    if not FOUND_DATASHADER:
         return
 
     show_all_data = True
     show_extreme_data = False
     optional_plots = False
 
     plotList = []
-    analysis = mcAnalysisBaseClass()
-    analysis.dataDir = path + "/scenario_AttFeedbackMC/"
+    analysis = mcAnalysisBaseClass(path + "/scenario_AttFeedbackMC/rerun/")
 
     # save_as_static: save off static .html files of the plots generated into the staticDir directory.
     # The staticDir will be created inside the dataDir folder.
@@ -192,13 +201,13 @@ def run(show_plots):
         plotList.extend(plotSuite(analysis.dataDir))
 
     if show_extreme_data:
-        analysis.variableName = "attGuidMsg.omega_BR_B"
+        analysis.variableName = "attGuidMsg.sigma_BR"
         analysis.variableDim = 1
 
-        extremaRunNumbers = analysis.getExtremaRunIndices(numExtrema=1, window=[500 * 1E9, 550 * 1E9])
+        extrema_run_numbers = analysis.getExtremaRunIndices(numExtrema=10, window=[1e9, 2e9])
 
-        analysis.extractSubsetOfRuns(runIdx=extremaRunNumbers)
-        plotList.extend(plotSuite(analysis.dataDir + "/subset"))
+        analysis.extractSubsetOfRuns(runIdx=extrema_run_numbers)
+        plotList.extend(plotSuite(analysis.dataDir + "subset/"))
 
     if optional_plots:
         # nominalRuns = analysis.getNominalRunIndices(50)
@@ -207,22 +216,25 @@ def run(show_plots):
         shadowFactor = pull_and_format_df(analysis.dataDir + "/eclipse_data_0.shadowFactor.data", 1)
         shadowFactor = shadowFactor.dropna(axis=1)
         shadowFactorPlot = DS_Plot(shadowFactor, title="Optional Plots: Eclipse",
-                                               xAxisLabel='time[s]', yAxisLabel='Eclipse Factor',
-                                               macro_x=macros.NANO2SEC, macro_y=macros.R2D,
-                                               cmap=RdYlBu9,
-                                               plotFcn=curve_per_df_component)
+                                   xAxisLabel='time[s]', yAxisLabel='Eclipse Factor',
+                                   macro_x=macros.NANO2SEC, macro_y=macros.R2D,
+                                   cmap=RdYlBu9,
+                                   plotFcn=curve_per_df_column)
 
         # plotList.extend([statPlots])
         plotList.extend([shadowFactorPlot])
 
+
+
+
     analysis.renderPlots(plotList)
 
 # The following must be commented out before this script can run.  It is provided here
 # to ensure that the sphinx documentation generation process does not run this script
 # automatically.
-if __name__ == "__main__":
-    run(False)
+#if __name__ == "__main__":
+#    run(False)
 
 
 # uncomment the following line to run this script.
-# run(False)
+run(False)

examples/MonteCarloExamples/scenarioRerunMonteCarlo.py

@@ -60,8 +60,8 @@ def run(time=None):
     scenarioName = "scenario_AttFeedback"
 
     monteCarlo = Controller()
-    monteCarlo.numProcess = 3 # Specify number of processes to spawn
-    runsList = [1]  # Specify the run numbers to be rerun
+    monteCarlo.numProcess = 10  # Specify number of processes to spawn
+    runsList = [678]  # Specify the run numbers to be rerun
 
     #
     # # Generic initialization
@@ -88,7 +88,7 @@ def run(time=None):
     # Step 4: Add any additional retention policies desired
     retentionPolicy = RetentionPolicy()
     retentionPolicy.logRate = int(2E9)
-    retentionPolicy.addMessageLog("attGuidMsg", ["sigma_BR"]) 
+    retentionPolicy.addMessageLog("attGuidMsg", ["sigma_BR"])
     monteCarlo.addRetentionPolicy(retentionPolicy)
 
 
@@ -97,6 +97,6 @@ def run(time=None):
 
 
 
+
 if __name__ == "__main__":
     run()
-

src/utilities/MonteCarlo/AnalysisBaseClass.py

@@ -5,18 +5,26 @@
 import numpy as np
 import pandas as pd
 from Basilisk.utilities import macros
+from bokeh.layouts import column
+from bokeh.plotting import figure, show
+import panel as pn
+import numpy as np
+from tornado.ioloop import IOLoop
 
 try:
     import holoviews as hv
+    import datashader as ds
     from Basilisk.utilities.datashader_utilities import DS_Plot, curve_per_df_component
+    from holoviews.operation.datashader import datashade, dynspread, spread
 except:
     pass
 
+
 class mcAnalysisBaseClass:
-    def __init__(self):
+    def __init__(self, data_dir=""):
         self.variableName = ""
         self.variableDim = 0
-        self.dataDir = ""
+        self.dataDir = data_dir
         self.numExtrema = 0
         self.extremaRuns = []
         self.timeWindow = []
@@ -66,17 +74,15 @@ def getExtremaRunIndices(self, numExtrema, window):
         times = self.data.index.tolist()
 
         # Find the closest indices to the time window requested
-        indStart = min(range(len(times)), key=lambda i: abs(times[i] - window[0]))
-        indEnd = min(range(len(times)), key=lambda i: abs(times[i] - window[1]))
-        self.timeWindow = [indStart, indEnd]
+        ind_start = min(range(len(times)), key=lambda i: abs(times[i] - window[0]))
+        ind_end = min(range(len(times)), key=lambda i: abs(times[i] - window[1]))
+        self.timeWindow = [times[ind_start], times[ind_end]]
 
         # Find outliers based on largest deviation off of the mean
-        self.mean = self.data.mean(axis=1, level=1)
-        self.diff = self.data.subtract(self.mean)
-        self.diff = self.diff.abs()
-        self.diff = self.diff.iloc[indStart:indEnd].max(axis=0)
-        self.extremaRuns = self.diff.nlargest(numExtrema).index._codes[0]
-        print("Extreme runs are ", list(dict.fromkeys(self.extremaRuns.tolist())))
+        mean = self.data.mean(axis=1)
+        diff = self.data.abs().sub(mean, axis=0)
+        self.extremaRuns = diff.transpose().nlargest(numExtrema, self.timeWindow).index
+        print("Extrema runs are ", list(dict.fromkeys(self.extremaRuns.tolist())))
         return self.extremaRuns
 
     def generateStatCurves(self):
@@ -121,27 +127,27 @@ def generateStatPlots(self):
         varIdxList = range(self.variableDim)
         varIdxListStr = str(varIdxList)
 
-        meanRun.columns = pd.MultiIndex.from_product([['mean'], [0,1,2]], names=["stats", "varIdx"])
-        medianRun.columns = pd.MultiIndex.from_product([['median'], [0,1,2]], names=["stats", "varIdx"])
-        stdRun.columns = pd.MultiIndex.from_product([['std'], [0,1,2]], names=["stats", "varIdx"])
+        meanRun.columns = pd.MultiIndex.from_product([['mean'], [0, 1, 2]], names=["stats", "varIdx"])
+        medianRun.columns = pd.MultiIndex.from_product([['median'], [0, 1, 2]], names=["stats", "varIdx"])
+        stdRun.columns = pd.MultiIndex.from_product([['std'], [0, 1, 2]], names=["stats", "varIdx"])
 
         meanRun_plot = DS_Plot(meanRun, title="Mean Plot: " + self.variableName,
-                                            xAxisLabel='time[s]', yAxisLabel= self.variableName.split('.')[-1],
-                                            macro_x=macros.NANO2SEC,
-                                            labels=['1', '2', '3'],
-                                            plotFcn=curve_per_df_component)
+                               xAxisLabel='time[s]', yAxisLabel=self.variableName.split('.')[-1],
+                               macro_x=macros.NANO2SEC,
+                               labels=['1', '2', '3'],
+                               plotFcn=curve_per_df_component)
 
         medRun_plot = DS_Plot(medianRun, title="Median Plot: " + self.variableName,
-                                            xAxisLabel='time[s]', yAxisLabel= self.variableName.split('.')[-1],
-                                            macro_x=macros.NANO2SEC,
-                                            labels=['1', '2', '3'],
-                                            plotFcn=curve_per_df_component)
+                              xAxisLabel='time[s]', yAxisLabel=self.variableName.split('.')[-1],
+                              macro_x=macros.NANO2SEC,
+                              labels=['1', '2', '3'],
+                              plotFcn=curve_per_df_component)
 
         stdRun_plot = DS_Plot(stdRun, title="Standard Dev Plot: " + self.variableName,
-                                            xAxisLabel='time[s]', yAxisLabel= self.variableName.split('.')[-1],
-                                            macro_x=macros.NANO2SEC,
-                                            labels=['1', '2', '3'],
-                                            plotFcn=curve_per_df_component)
+                              xAxisLabel='time[s]', yAxisLabel=self.variableName.split('.')[-1],
+                              macro_x=macros.NANO2SEC,
+                              labels=['1', '2', '3'],
+                              plotFcn=curve_per_df_component)
 
         statRun_plots = []
         statRun_plots.append(meanRun_plot)
@@ -160,7 +166,11 @@ def extractSubsetOfRuns(self, runIdx):
         """
         idx = pd.IndexSlice
         baseDir = self.dataDir
+        new_list = []
+        for run in runIdx:
+            new_list.append(run[0])
 
+        runIdx = new_list;
         # check if a subset directory exists, and if it already contains all runIdx requested
         if not os.path.exists(baseDir + "/subset/"):
             os.mkdir(baseDir + "/subset/")
@@ -196,36 +206,36 @@ def extractSubsetOfRuns(self, runIdx):
             pd.to_pickle(dfSubSet, baseDir + "/subset/" + varName[-1])
         print("Finished Populating Subset Directory")
 
-    def renderPlots(self, plotList):
+    def renderPlots(self, plotList, cols=2):
         """
         Render all plots in plotList and print information about time taken, percent complete, which plot, etc.
 
         :param plotList: List of plots to render
         :return: nothing.
         """
         hv.extension('bokeh')
-        renderer = hv.renderer('bokeh').instance(mode='server')
+        #renderer = hv.renderer('bokeh')
+        #renderer = hv.renderer.instance(mode='server')
 
         if self.save_as_static:
-            print("Note: You requested to save static plots. This means no interactive python session will be generated.")
+            print(
+                "Note: You requested to save static plots. This means no interactive python session will be generated.")
         print("Beginning the plotting")
 
         if not os.path.exists(self.dataDir + self.staticDir):
             os.mkdir(self.dataDir + self.staticDir)
 
+        figures = []
         for i in range(len(plotList)):
             startTime = time.time()
-            image, title = plotList[i].generateImage()
-            try:
-                if self.save_as_static:
-                    # Save .html files of each of the plots into the static directory
-                    hv.save(image, self.dataDir + self.staticDir + "/" + title + ".html")
-                else:
-                    renderer.server_doc(image)
-                # Print information about the rendering process
-                print("LOADED: " + title +"\t\t\t" +
-                      "Percent Complete: " + str(round((i + 1) / len(plotList) * 100, 2)) + "% \t\t\t"
-                      "Time Elapsed: " + str( round(time.time() - startTime)) + " [s]")
-            except Exception as e:
-                print("Couldn't Plot " + title)
-                print(e)
+            fig, title = plotList[i].generateOverlay()
+            figures.append(fig)
+            print("LOADED: " + title + "\t\t\t" +
+                  "Percent Complete: " + str(round((i + 1) / len(plotList) * 100, 2)) + "% \t\t\t"
+                                                                                        "Time Elapsed: " + str(
+                round(time.time() - startTime)) + " [s]")
+
+        layout = hv.Layout(figures).cols(cols)
+
+
+        return pn.panel(layout).servable()