bokeh.io.output_file

    def save(self, title='Training Results'):
        if len(self.figures) > 0:
            if os.path.isfile(self.plot_path):
                os.remove(self.plot_path)
            output_file(self.plot_path, title=title)
            plot = column(*self.figures)
            save(plot)
            self.figures = []
        self.results.to_csv(self.path, index=False, index_label=False)

    def load(self, path=None):

def output_file_url(request, file_server):
    from bokeh.io import output_file
    filename = request.function.__name__ + '.html'
    file_obj = request.fspath.dirpath().join(filename)
    file_path = file_obj.strpath
    url = file_path.replace('\\', '/')  # Windows-proof

    output_file(file_path, mode='inline')

    def tear_down():
        if file_obj.isfile():
            file_obj.remove()
    request.addfinalizer(tear_down)

    return file_server.where_is(url)

@pytest.fixture

def save_plot(plot, name: str):
    base = Path('test-outputs')
    path = base / Path(name)
    path.parent.mkdir(exist_ok=True, parents=True)
    suf = path.suffix
    if suf == '.html':
        from bokeh.io import output_file, save
        output_file(str(path), title='hello', mode='inline', root_dir=None)
        save(plot)
    elif suf == '.png':
        # todo sigh.. seems that png export is way too slow
        from bokeh.io import export_png
        export_png(plot, filename=str(path))
    else:
        raise RuntimeError(name, suf)


def make_test(plot_factory):

def render_tab(*, tab: Tab, filename: Path):
    res = tab.plotter()

    from bokeh.io import save, output_file, curdoc
    output_file(filename, title=tab.name, mode='inline', root_dir=None)
    curdoc().theme = theme
    # TODO a bit weird that it needs two function calls to save..
    save(res)


def run(to: Path, tab_name: Optional[str]=None, debug: bool=False) -> Iterable[Exception]:

def _bokeh_reset(filename=None):
    curstate().reset()
    if filename:
        if not filename.endswith('.html'):
            filename += '.html'
        output_file(filename, title=filename)
    elif IS_JUPYTER_NOTEBOOK:
        curstate().output_notebook()


def colorgen():

def build_trace_html(trace_events, workers, filename, plot_fn):
    """
    Render trace events into a HTML file according to the given plot function.
    :type plot_fn: Callable[[List[Event], List[Worker]], bokeh.plotting.figure.Figure]
    """
    import bokeh.io

    trace_events = normalize_events(trace_events)

    plot = plot_fn(trace_events, workers)

    bokeh.io.output_file(filename)
    bokeh.io.save(plot)


def simulator_trace_to_html(simulator, filename):

def activity_details(local_time, activity_group):

    f'''
    # Activity Details: {local_time} ({activity_group})
    '''

    '''
    $contents
    '''

    '''
    ## Load Data
    
    Open a connection to the database and load the data we require.
    '''

    s = session('-v2')

    activity = std_activity_statistics(s, activity_journal=local_time, activity_group=activity_group)
    health = std_health_statistics(s)
    hr_zones = hr_zones_from_database(s, local_time, activity_group)
    climbs = Statistics(s, sources=climb_sources(s, local_time, activity_group=activity_group)). \
        by_name(SectorCalculator, N.CLIMB_ANY, N.VERTICAL_POWER, like=True).with_. \
        copy_with_units().df
    active = Statistics(s, activity_journal=local_time, activity_group=activity_group). \
        by_name(ActivityCalculator, N.ACTIVE_TIME, N.ACTIVE_DISTANCE). \
        with_.copy_with_units().df.append(climbs)

    f'''
    ## Activity Plots
    
    To the right of each plot of data against distance is a related plot of cumulative data
    (except the last, cadence, which isn't useful and so replaced by HR zones).
    Green and red areas indicate differences between the two dates. 
    Additional red lines on the altitude plot are auto-detected climbs.
    
    Plot tools support zoom, dragging, etc.
    '''

    output_file(filename='/dev/null')

    sp = comparison_line_plot(700, 200, N.DISTANCE_KM, N.MED_SPEED_KMH, activity, ylo=0)
    add_climb_zones(sp, climbs, activity)
    sp_c = cumulative_plot(200, 200, N.MED_SPEED_KMH, activity, ylo=0)
    xrange = sp.x_range if sp else None

    el = comparison_line_plot(700, 200, N.DISTANCE_KM, N.ELEVATION_M, activity, x_range=xrange)
    add_climbs(el, climbs, activity)
    el_c = cumulative_plot(200, 200, N.CLIMB_MS, activity)
    xrange = xrange or (el.x_range if el else None)

    hri = comparison_line_plot(700, 200, N.DISTANCE_KM, N.HR_IMPULSE_10, activity, ylo=0, x_range=xrange)
    add_climb_zones(hri, climbs, activity)
    hri_c = cumulative_plot(200, 200, N.HR_IMPULSE_10, activity, ylo=0)
    xrange = xrange or (hri.x_range if hri else None)

    hr = comparison_line_plot(700, 200, N.DISTANCE_KM, N.HEART_RATE_BPM, activity, x_range=xrange)
    add_hr_zones(hr, activity, N.DISTANCE_KM, hr_zones)
    add_climb_zones(hr, climbs, activity)
    hr_c = cumulative_plot(200, 200, N.HEART_RATE_BPM, activity)
    xrange = xrange or (hr.x_range if hr else None)

    pw = comparison_line_plot(700, 200, N.DISTANCE_KM, N.MED_POWER_ESTIMATE_W, activity, ylo=0, x_range=xrange)
    pw.varea(source=activity, x=N.DISTANCE_KM, y1=0, y2=N.MED_VERTICAL_POWER_W,
             level='underlay', color='black', fill_alpha=0.25)
    add_climb_zones(pw, climbs, activity)
    pw_c = cumulative_plot(200, 200, N.MED_POWER_ESTIMATE_W, activity, ylo=0)
    xrange = xrange or (pw.x_range if pw else None)

    cd = comparison_line_plot(700, 200, N.DISTANCE_KM, N.MED_CADENCE_RPM, activity, ylo=0, x_range=xrange)
    add_climb_zones(cd, climbs, activity)
    hr_h = histogram_plot(200, 200, N.HR_ZONE, activity, xlo=1, xhi=5)

    show(gridplot([[el, el_c], [sp, sp_c], [hri, hri_c], [hr, hr_c], [pw, pw_c], [cd, hr_h]]))

    '''
    ## Activity Maps
    '''

    map = map_plot(400, 400, activity)
    m_el = map_intensity_signed(200, 200, activity, N.GRADE_PC, ranges=map, power=0.5)
    m_sp = map_intensity(200, 200, activity, N.MED_SPEED_KMH, ranges=map, power=2)
    m_hr = map_intensity(200, 200, activity, N.HR_IMPULSE_10, ranges=map)
    m_pw = map_intensity(200, 200, activity, N.MED_POWER_ESTIMATE_W, ranges=map)
    show(row(map, gridplot([[m_el, m_sp], [m_hr, m_pw]], toolbar_location='right')))

    '''
    ## Activity Statistics
    '''

    '''
    Active time and distance exclude pauses.
    '''

    active[[N.ACTIVE_TIME_S, N.ACTIVE_DISTANCE_KM]].dropna(). \
        transform({N.ACTIVE_TIME_S: format_seconds, N.ACTIVE_DISTANCE_KM: format_km})

    '''
    Climbs are auto-detected and shown only for the main activity. They are included in the elevation plot above.
    '''

    if present(climbs, N.CLIMB_TIME):
        display(transform(climbs,
                          {N.CLIMB_TIME: format_seconds, N.CLIMB_ELEVATION: format_metres,
                           N.CLIMB_DISTANCE: format_km, N.CLIMB_GRADIENT: format_percent,
                           N.VERTICAL_POWER: format_watts, N.CLIMB_CATEGORY: lambda x: x}))

    '''
    ## Health and Fitness
    '''

    fitness, fatigue = like(N.FITNESS_ANY, health.columns), like(N.FATIGUE_ANY, health.columns)
    colours = ['black'] * len(fitness) + ['red'] * len(fatigue)
    alphas = [1.0] * len(fitness) + [0.5] * len(fatigue)
    ff = multi_line_plot(900, 300, N.TIME, fitness + fatigue, health, colours, alphas=alphas)
    xrange = ff.x_range if ff else None
    add_multi_line_at_index(ff, N.TIME, fitness + fatigue, health, colours, alphas=alphas, index=-1)
    atd = std_distance_time_plot(900, 200, health, x_range=xrange)
    show(gridplot([[ff], [atd]]))

def calendar():

    '''
    # Calendar
    '''

    '''
    $contents
    '''

    s = session('-v2')
    output_file(filename='/dev/null')

    '''
    ## Distance

    Larger distances have larger symbols.
    
    Place the cursor over the symbol for more information.
    '''

    df = Statistics(s). \
        by_name(ActivityTopic, N.NAME). \
        by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME, N.TOTAL_CLIMB).with_. \
        copy({N.ACTIVE_DISTANCE: N.ACTIVE_DISTANCE_KM}).add_times().df
    df['Duration'] = df[N.ACTIVE_TIME].map(format_seconds)
    if present(df, N.TOTAL_CLIMB):
        df.loc[df[N.TOTAL_CLIMB].isna(), [N.TOTAL_CLIMB]] = 0

    calendar = Calendar(df, title=N.DISTANCE, not_hover=[N.ACTIVE_DISTANCE, N.ACTIVE_TIME])
    calendar.std_distance()

    '''
    ## Work Done and Fatigue

    Larger increases in Fitness have larger symbols.  Higher fatigue (relative to fitness) is redder.
    
    Place the cursor over the symbol for more information.
    '''

    df = Statistics(s). \
        by_name(ActivityTopic, N.NAME). \
        by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME, N.TOTAL_CLIMB). \
        by_name(ActivityCalculator, N._delta(N.FITNESS_ANY), like=True).df
    df = Statistics(s). \
        by_name(ResponseCalculator, N.ANY_FATIGUE_ANY, N.ANY_FITNESS_ANY, like=True). \
        with_.drop_prefix(N.DEFAULT).into(df, tolerance='30m')
    df['Duration'] = df[N.ACTIVE_TIME].map(format_seconds)
    df.loc[df[N.TOTAL_CLIMB].isna(), N.TOTAL_CLIMB] = 0
    work_done = sorted_numeric_labels(df.columns, N._delta(N.FITNESS))[0]
    fitness = sorted_numeric_labels(df.columns, N.FITNESS)[0]
    fatigue = sorted_numeric_labels(df.columns, N.FATIGUE)[0]
    print(fatigue, fitness)
    df['FF Ratio'] = df[fatigue] / df[fitness]

    calendar = Calendar(df, title='Work Done and Fatigue',
                        not_hover=[N.ACTIVE_DISTANCE, N.ACTIVE_TIME] +
                                  [column for column in df.columns if ':' in column])
    calendar.background('square', fill_alpha=0, line_alpha=1, color='lightgrey')
    calendar.set_palette('FF Ratio', K2R, lo=0.5, hi=2)
    calendar.set_size(work_done, min=0.1, gamma=0.5)
    calendar.foreground('square', fill_alpha=1, line_alpha=0)
    calendar.show()

    '''
    ## Distance, Climb and Direction

    Larger distances have larger symbols.  Higher climbs are redder.  
    The arc indicates the general direction relative to the start.
    
    Place the cursor over the symbol for more information.
    '''

    df = Statistics(s). \
        by_name(ActivityTopic, N.NAME). \
        by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME, N.TOTAL_CLIMB, N.DIRECTION,
                N.ASPECT_RATIO).df
    df[N.DISTANCE_KM] = df[N.ACTIVE_DISTANCE]
    df['Duration'] = df[N.ACTIVE_TIME].map(format_seconds)
    if present(df, N.TOTAL_CLIMB):
        df.loc[df[N.TOTAL_CLIMB].isna(), N.TOTAL_CLIMB] = 0

    calendar = Calendar(df, title='Distance, Climb and Direction',
                        not_hover=[N.ACTIVE_DISTANCE, N.ACTIVE_TIME] +
                                  [column for column in df.columns if ':' in column])
    calendar.std_distance_climb_direction()

    '''
    ## Distance, Work Done and Direction

    Larger distances have larger symbols.  Larger gains in fitness are redder.  
    The arc indicates the general direction relative to the start.
    
    Place the cursor over the symbol for more information.
    '''

    df = Statistics(s). \
        by_name(ActivityTopic, N.NAME). \
        by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME, N.TOTAL_CLIMB, N.DIRECTION,
                N.ASPECT_RATIO). \
        by_name(ActivityCalculator, N._delta(N.FITNESS_ANY), like=True).df
    df[N.DISTANCE_KM] = df[N.ACTIVE_DISTANCE]
    df['Duration'] = df[N.ACTIVE_TIME].map(format_seconds)
    df.loc[df[N.TOTAL_CLIMB].isna(), N.TOTAL_CLIMB] = 0

    calendar = Calendar(df, title='Distance, Fitness and Direction',
                        not_hover=[N.ACTIVE_DISTANCE, N.ACTIVE_TIME] +
                                  [column for column in df.columns if ':' in column])
    calendar.std_distance_fitness_direction()

    '''
    ## Fitness and Fatigue

    Better fitness has larger symbols.  When fatigue is higher symbols have "hotter" colours.
    
    Place the cursor over the symbol for more information.
    '''

    df = Statistics(s). \
        by_name(ActivityTopic, N.NAME). \
        by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME). \
        by_name(ActivityCalculator, N._delta(N.FITNESS_ANY), like=True).df
    df = Statistics(s). \
        by_name(ResponseCalculator, N.ANY_FATIGUE_ANY, N.ANY_FITNESS_ANY, like=True). \
        with_.drop_prefix(N.DEFAULT).into(df, tolerance='30m')
    fitness = sorted_numeric_labels(df.columns, N.FITNESS)[0]
    fatigue = sorted_numeric_labels(df.columns, N.FATIGUE)[0]
    df['FF Ratio'] = df[fatigue] / df[fitness]

    calendar = Calendar(df, title='Fitness and Fatigue', scale=18, border_month=0, border_day=0)
    calendar.set_size(fitness, min=0.1, gamma=0.5)
    calendar.set_palette('FF Ratio', magma(256), lo=0.5, hi=2, min=0)
    calendar.foreground('square', fill_alpha=1, line_alpha=0)
    calendar.show()

    '''
    ## Groups, Distance, Climb and Direction

    Larger distances have larger symbols.  Higher climbs are lighter.  
    The arc indicates the general direction relative to the start.
    Pastel backgrounds group similar rides.
    
    Place the cursor over the symbol for more information.
    '''

    dfa = Statistics(s). \
        by_name(ActivityTopic, N.NAME). \
        by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME, N.TOTAL_CLIMB, N.DIRECTION,
                N.ASPECT_RATIO). \
        by_name(ActivityCalculator, N._delta(N.FITNESS_ANY), like=True).df
    dfa[N.DISTANCE_KM] = df[N.ACTIVE_DISTANCE]
    dfa['Duration'] = df[N.ACTIVE_TIME].map(format_seconds)
    dfa.loc[dfa[N.TOTAL_CLIMB].isna(), N.TOTAL_CLIMB] = 0
    dfb = groups_by_time(s)
    if present(dfb, N.GROUP):
        dfb.loc[dfb[N.GROUP].isna(), N.GROUP] = -1
        df = dfa.join(dfb)

        calendar = Calendar(df, scale=15, border_day=0.1,
                            not_hover=[N.ACTIVE_DISTANCE, N.ACTIVE_TIME] +
                                      [column for column in df.columns if ':' in column])
        calendar.std_group_distance_climb_direction()

def compare_activities(local_time, compare_time, activity_group):

    f'''
    # Compare Activities: {local_time} v {compare_time} ({activity_group})
    '''

    '''
    $contents
    '''

    '''
    ## Load Data
    
    Open a connection to the database and load the data we require.
    '''

    s = session('-v2')

    activity = std_activity_statistics(s, activity_journal=local_time, activity_group=activity_group)
    compare = std_activity_statistics(s, activity_journal=compare_time, activity_group=activity_group)
    health = std_health_statistics(s)
    hr_zones = hr_zones_from_database(s, local_time, activity_group)
    climbs = Statistics(s, sources=climb_sources(s, local_time, activity_group=activity_group)). \
        by_name(SectorCalculator, N.CLIMB_ANY, N.VERTICAL_POWER, like=True).with_. \
        copy_with_units().df
    active = Statistics(s, activity_journal=local_time, activity_group=activity_group). \
        by_name(ActivityCalculator, N.ACTIVE_TIME, N.ACTIVE_DISTANCE). \
        with_.copy_with_units().df.append(climbs)

    f'''
    ## Activity Plots
    
    The black line shows data from {local_time}, 
    the grey line from {compare_time}. 
    To the right of each plot of data against distance is a related plot of cumulative data
    (except the last, cadence, which isn't useful and so replaced by HR zones).
    Green and red areas indicate differences between the two dates. 
    Additional red lines on the altitude plot are auto-detected climbs.
    
    Plot tools support zoom, dragging, etc.
    '''

    output_file(filename='/dev/null')

    sp = comparison_line_plot(700, 200, N.DISTANCE_KM, N.MED_SPEED_KMH, activity, other=compare, ylo=0)
    add_climb_zones(sp, climbs, activity)
    sp_c = cumulative_plot(200, 200, N.MED_SPEED_KMH, activity, other=compare, ylo=0)
    xrange = sp.x_range if sp else None

    el = comparison_line_plot(700, 200, N.DISTANCE_KM, N.ELEVATION_M, activity, other=compare, x_range=xrange)
    add_climbs(el, climbs, activity)
    el_c = cumulative_plot(200, 200, N.CLIMB_MS, activity, other=compare)
    xrange = xrange or (el.x_range if el else None)

    hri = comparison_line_plot(700, 200, N.DISTANCE_KM, N.HR_IMPULSE_10, activity, other=compare, ylo=0, x_range=xrange)
    add_climb_zones(hri, climbs, activity)
    hri_c = cumulative_plot(200, 200, N.HR_IMPULSE_10, activity, other=compare, ylo=0)
    xrange = xrange or (hri.x_range if hri else None)

    hr = comparison_line_plot(700, 200, N.DISTANCE_KM, N.HEART_RATE_BPM, activity, other=compare, x_range=xrange)
    add_hr_zones(hr, activity, N.DISTANCE_KM, hr_zones)
    add_climb_zones(hr, climbs, activity)
    hr_c = cumulative_plot(200, 200, N.HEART_RATE_BPM, activity, other=compare)
    xrange = xrange or (hr.x_range if hr else None)

    pw = comparison_line_plot(700, 200, N.DISTANCE_KM, N.MED_POWER_ESTIMATE_W, activity, other=compare, ylo=0, x_range=xrange)
    add_climb_zones(pw, climbs, activity)
    pw_c = cumulative_plot(200, 200, N.MED_POWER_ESTIMATE_W, activity, other=compare, ylo=0)
    xrange = xrange or (pw.x_range if pw else None)

    cd = comparison_line_plot(700, 200, N.DISTANCE_KM, N.MED_CADENCE_RPM, activity, other=compare, ylo=0, x_range=xrange)
    add_climb_zones(cd, climbs, activity)
    hr_h = histogram_plot(200, 200, N.HR_ZONE, activity, xlo=1, xhi=5)

    show(gridplot([[el, el_c], [sp, sp_c], [hri, hri_c], [hr, hr_c], [pw, pw_c], [cd, hr_h]]))

    '''
    ## Activity Maps
    '''

    map = map_plot(400, 400, activity, other=compare)
    m_el = map_intensity_signed(200, 200, activity, N.GRADE_PC, ranges=map, power=0.5)
    m_sp = map_intensity(200, 200, activity, N.MED_SPEED_KMH, ranges=map, power=2)
    m_hr = map_intensity(200, 200, activity, N.HR_IMPULSE_10, ranges=map)
    m_pw = map_intensity(200, 200, activity, N.MED_POWER_ESTIMATE_W, ranges=map)
    show(row(map, gridplot([[m_el, m_sp], [m_hr, m_pw]], toolbar_location='right')))

    '''
    ## Activity Statistics
    '''

    '''
    Active time and distance exclude pauses.
    '''

    active[[N.ACTIVE_TIME_S, N.ACTIVE_DISTANCE_KM]].dropna(). \
        transform({N.ACTIVE_TIME_S: format_seconds, N.ACTIVE_DISTANCE_KM: format_km})

    '''
    Climbs are auto-detected and shown only for the main activity. They are included in the elevation plot above.
    '''

    if present(climbs, N.CLIMB_TIME):
        display(transform(climbs,
                          {N.CLIMB_TIME: format_seconds, N.CLIMB_ELEVATION: format_metres,
                           N.CLIMB_DISTANCE: format_km, N.CLIMB_GRADIENT: format_percent,
                           N.VERTICAL_POWER: format_watts, N.CLIMB_CATEGORY: lambda x: x}))

    '''
    ## Health and Fitness
    '''

    fitness, fatigue = like(N.FITNESS_ANY, health.columns), like(N.FATIGUE_ANY, health.columns)
    colours = ['black'] * len(fitness) + ['red'] * len(fatigue)
    alphas = [1.0] * len(fitness) + [0.5] * len(fatigue)
    ff = multi_line_plot(900, 300, N.TIME, fitness + fatigue, health, colours, alphas=alphas)
    xrange = ff.x_range if ff else None
    add_multi_line_at_index(ff, N.TIME, fitness + fatigue, health, colours, alphas=alphas, index=-1)
    atd = std_distance_time_plot(900, 200, health, x_range=xrange)
    show(gridplot([[ff], [atd]]))

def fit_ff_segments(group, *segment_titles):

    f'''
    # Fit FF Parameters for {group} to {', '.join(segment_titles)}

    This notebook allows you to estimate a personal time scale (decay period) for the
    [FF model](https://andrewcooke.github.io/choochoo/impulse) using your times (more exactly, the speed)
    on the named segments.

    The resulting value can be used in the configuration so that the Fitness parameter more accurately
    reflects your personal rate of adaption.

    For more information see the [documentation](https://andrewcooke.github.io/choochoo/ff-fitting).
    '''

    '''
    $contents
    '''

    '''
    ## Load Data
    
    Open a connection to the database and load the data we require.  
    We reduce the HR data to hourly values so that we have smaller arrays (for faster processing).
    '''

    s = session('-v2')

    # hr10 = statistics(s, HR_IMPULSE_10, activity_group=ActivityGroup.from_name(s, 'all'))
    hr10 = Statistics(s).by_name(ImpulseCalculator, N.HR_IMPULSE_10).df
    print(hr10.describe())
    segments = [s.query(Segment).filter(Segment.title == segment_title).one() for segment_title in segment_titles]
    for segment in segments:
        print(segment.title, segment.distance)
    kit_statistic = StatisticName.from_name(s, ActivityReader.KIT, ActivityReader)
    journals_by_kit_by_segment = \
        {segment: group_to_dict(s.query(StatisticJournalText.value, SegmentJournal).
                                join(ActivityJournal, SegmentJournal.activity_journal_id == ActivityJournal.id).
                                join(StatisticJournalText, StatisticJournalText.source_id == ActivityJournal.id).
                                filter(StatisticJournalText.statistic_name_id == kit_statistic.id,
                                       SegmentJournal.segment == segment).all())
         for segment in segments}
    times_by_kit_by_segment = \
        {segment: {kit: drop_empty(statistics(s, SEGMENT_TIME, sources=journals)).dropna()
                   for kit, journals in journals_by_kit_by_segment[segment].items()}
         for segment in segments}
    performances = []
    for segment in segments:
        for kit, times in times_by_kit_by_segment[segment].items():
            if len(times.columns) and len(times.index) > 2:
                times.index = times.index.round('1H')
                performances.append(Series(segment.distance / times.iloc[:, 0], times.index,
                                           name=f'{segment.title}/{kit}'))
    n_performances = sum(len(performance) for performance in performances)

    hr3600 = sum_to_hour(hr10, HR_IMPULSE_10)

    def copy_of_performances():
        return [performance.copy() for performance in performances]

    for performance in performances:
        print(performance)

    '''
    ## Define Plot Routine

    Below the `params` are (currently) log10_period and log10_start - the time period for decay and the initial
    response (fitness) value.  We use log10 so that the values cannot be negative.
    '''

    output_file(filename='/dev/null')

    def plot_params(params, rejected=tuple()):
        response = calc_response(hr3600, params)
        predicted = calc_observed(restrict_response(response, performances), performances)
        f = figure(title=fmt_params(params), plot_width=500, plot_height=450, x_axis_type='datetime')
        f.line(x=response.index, y=response, color='grey')
        for color, prediction, performance in zip(evenly_spaced_hues(len(performances)), predicted, performances):
            f.circle(x=prediction.index, y=prediction, color=color, legend_label=performance.name)
        for i, t in rejected:
            f.x(x=t, y=predicted[i].loc[t], color='black', size=10)
        f.legend.location = 'bottom_left'
        show(f)

    '''
    ## Plot Initial Response
    
    Trim and resample data; evaluate and plot our initial model.
    '''

    initial_period = log10(42 * 24)
    plot_params((initial_period,))

    '''
    ## Explore Effect of Start
    '''
    plot_params((initial_period, log10(1000)))
    plot_params((initial_period, log10(3000)))

    '''
    ## Fit Model using L1
    
    Adjust tol according to the 'fun' value in the result (tol is the tolerance in that value).
    
    Note how the period (printed, in hours) varies as points are rejected.
    '''

    result, rejected = fit_ff_params(hr3600, (initial_period,), copy_of_performances(),
                                     method='L1', tol=0.01,
                                     max_reject=n_performances // 2, threshold=(2, 0.2))
    print(result)
    plot_params(result.x, rejected=rejected)

    result, rejected = fit_ff_params(hr3600, (initial_period, 0), copy_of_performances(),
                                     method='L1', tol=0.01,
                                     max_reject=n_performances // 2, threshold=(2, 0.2))
    print(result)
    plot_params(result.x, rejected=rejected)

    '''
    ## Fit Model using L2
    
    Using different methods gives us some idea of how susceptible the value is to processing assumptions.
    
    Note that the thresholds change for the L2 norm.
    '''

    initial_period = log10(42 * 24)
    result, rejected = fit_ff_params(hr3600, (initial_period,), copy_of_performances(),
                                     method='L2', tol=0.01,
                                     max_reject=n_performances // 2, threshold=(500, 50))
    print(result)
    plot_params(result.x, rejected=rejected)

    result, rejected = fit_ff_params(hr3600, (initial_period, 0), copy_of_performances(),
                                     method='L2', tol=0.01,
                                     max_reject=n_performances // 2, threshold=(500, 50))
    print(result)
    plot_params(result.x, rejected=rejected)

def health():

    '''
    # Health
    '''

    '''
    $contents
    '''

    '''
    ## Load Data
    
    Open a connection to the database and load the data we require.
    '''

    s = session('-v2')
    health = std_health_statistics(s)

    '''
    ## Health and Fitness
    '''

    output_file(filename='/dev/null')

    fitness, fatigue = like(N.FITNESS_ANY, health.columns), like(N.FATIGUE_ANY, health.columns)
    colours = ['black'] * len(fitness) + ['red'] * len(fatigue)
    alphas = [1.0] * len(fitness) + [0.5] * len(fatigue)
    ff = multi_line_plot(900, 300, N.TIME, fitness + fatigue, health, colours, alphas=alphas)
    xrange = ff.x_range if ff else None
    add_multi_line_at_index(ff, N.TIME, fitness + fatigue, health, colours, alphas=alphas, index=-1)
    atd = std_distance_time_plot(900, 200, health, x_range=xrange)
    show(gridplot([[ff], [atd]]))

def month(month):

    f'''
    # Month: {month}
    '''

    '''
    $contents
    '''

    '''
    ## Preparation
    '''
    
    s = session('-v2')
    output_file(filename='/dev/null')
    map_size = 100
    month_start = to_date(month).replace(day=1)

    '''
    ## Generate Plot
    '''

    def days():

        for i in Calendar().iterweekdays():
            yield Div(text=f'  <  h2>{day_name[i]} < /h2>')

        day = month_start - dt.timedelta(days=month_start.weekday())
        while day.replace(day=1)  < = month_start:
            for weekday in range(7):
                if day.month == month_start.month:
                    contents = [Div(text=f' < h1>{day.strftime("%d")} < /h1>')]
                    for a in s.query(ActivityJournal). \
                            filter(ActivityJournal.start >= local_date_to_time(day),
                                   ActivityJournal.start  <  local_date_to_time(day + dt.timedelta(days=1))).all():
                        df = Statistics(s, activity_journal=a). \
                            by_name(ActivityReader, N.SPHERICAL_MERCATOR_X, N.SPHERICAL_MERCATOR_Y).df
                        contents.append(map_thumbnail(map_size, map_size, df, title=False))
                        df = Statistics(s, activity_journal=a). \
                            by_name(ActivityCalculator, N.ACTIVE_DISTANCE, N.ACTIVE_TIME).df
                        contents.append(Div(
                            text=f'{format_km(df[N.ACTIVE_DISTANCE][0])} {format_seconds(df[N.ACTIVE_TIME][0])}'))
                else:
                    contents = [Spacer()]
                yield column(contents)
                day += dt.timedelta(days=1)

    show(grid(list(days()), ncols=7))

def power_v_hr(local_time, activity_group):

    f'''
    # Power v HR ({local_time}, {activity_group})
    '''

    '''
    $contents
    '''

    '''
    ## Load Data
    
    Open a connection to the database and load the data we require.
    '''

    local_times = ['2021-04-18 05:50:52', '2021-04-11 05:52:16', '2021-04-03 05:59:03']
    dt = 60

    s = session('-v2')
    # stats = Statistics(s, activity_journal=local_time, activity_group=activity_group, with_timespan=True)
    sources = [ActivityJournal.at(s, local_time, activity_group=activity_group) for local_time in local_times]
    stats = Statistics(s, sources=sources, with_timespan=True, with_sources=True)
    df = stats.by_name(ElevationCalculator, N.ELEVATION, N.GRADE). \
        by_name(ActivityReader, N.DISTANCE, N.SPEED, N.HEART_RATE). \
        with_.transform(N.DISTANCE, scale=1000). \
        concat(N.DISTANCE, 100). \
        concat_index(300). \
        without_sources().df
    # TODO - change default to use just the N.HR_IMPULSE_10
    fthr_df = Statistics(s).by_name(Constant, N.FTHR).df
    hr_zones = hr_zones_from_database(s, local_time, activity_group)

    hr_zone(df, fthr_df)
    ldf = linear_resample_time(df, dt=dt)
    ldf = add_differentials(ldf, max_gap=1.1 * dt)

    '''
    ## Define Model
    '''

    mass = 65 + 9
    params = {'cda': 0.42, 'crr': 0.0055, 'gamma': 2, 'zero': 2}
    bounds = {'cda': (0.1, 0.5), 'crr': (0.0, 0.01), 'gamma': (0.5, 2.5), 'zero': (1, 3)}

    def model(params, df):
        idf = impulse_10(df, HRImpulse('HR Impulse', gamma=params['gamma'], zero=params['zero'], one=6, max_secs=60))
        idf.rolling('60S').apply(lambda g: trapz(g, g.index.view(np.int64) / 10**9))
        mdf = ldf.copy(deep=True)
        mdf[N.HR_IMPULSE_10] = interp1d(idf.index.view(np.int64), idf[N.HR_IMPULSE_10], bounds_error=False)(mdf.index.view(np.int64))
        mdf = add_energy_budget(mdf, mass)
        mdf = add_loss_estimate(mdf, mass, cda=params['cda'], crr=params['crr'])
        mdf = add_power_estimate(mdf)
        mdf[N.POWER_ESTIMATE].clip(lower=0, inplace=True)
        return mdf

    def cost(mdf):
        # this fits a line through the origin (the param 'zero' is the offset, effectively)
        mdf.dropna(inplace=True)
        x = np.array(mdf[N.HR_IMPULSE_10])[:, np.newaxis]
        y = mdf[N.POWER_ESTIMATE]
        _, r, _, _ = np.linalg.lstsq(x, y, rcond=None)
        return r[0]

    '''
    ## Check
    '''

    output_file(filename='/dev/null')

    def plot(mdf):

        el = comparison_line_plot(700, 200, N.DISTANCE, N.ELEVATION, mdf)
        xrange = el.x_range if el else None

        hri = comparison_line_plot(700, 200, N.DISTANCE, N.HR_IMPULSE_10, mdf, ylo=0, x_range=xrange)
        xrange = xrange or (hri.x_range if hri else None)

        hr = comparison_line_plot(700, 200, N.DISTANCE, N.HEART_RATE, mdf, x_range=xrange)
        add_hr_zones(hr, mdf, N.DISTANCE, hr_zones)
        xrange = xrange or (hr.x_range if hr else None)

        pw = comparison_line_plot(700, 200, N.DISTANCE, N.POWER_ESTIMATE, mdf, ylo=0, x_range=xrange)
        pw.varea(source=mdf, x=N.DISTANCE, y1=0, y2=N.VERTICAL_POWER,
                 level='underlay', color='black', fill_alpha=0.25)

        show(gridplot([[el], [hri], [hr], [pw]]))

    plot(model(params, df))

    '''
    ## Fit Data
    '''

    params = {'cda': 0.42, 'crr': 0.0055, 'gamma': 2, 'zero': 2}
    fit(params, df, model, cost, bounds=bounds)
    plot(model(params, df))
    params

def generate_rec_prec_html(mrec, mprec, scores, class_name, ap):
    """
     generate dynamic P-R curve HTML page for each class
    """
    # bypass invalid class
    if len(mrec) == 0 or len(mprec) == 0 or len(scores) == 0:
        return

    rec_prec_plot_path = os.path.join('result' ,'classes')
    os.makedirs(rec_prec_plot_path, exist_ok=True)
    bokeh_io.output_file(os.path.join(rec_prec_plot_path, class_name + '.html'), title='P-R curve for ' + class_name)

    # prepare curve data
    area_under_curve_x = mrec[:-1] + [mrec[-2]] + [mrec[-1]]
    area_under_curve_y = mprec[:-1] + [0.0] + [mprec[-1]]
    score_on_curve = [0.0] + scores[:-1] + [0.0] + [scores[-1]] + [1.0]
    source = bokeh.models.ColumnDataSource(data={
      'rec'      : area_under_curve_x,
      'prec' : area_under_curve_y,
      'score' : score_on_curve,
    })

    # prepare plot figure
    plt_title = 'class: ' + class_name + ' AP = {}%'.format(ap*100)
    plt = bokeh_plotting.figure(plot_height=200 ,plot_width=200, tools="", toolbar_location=None,
               title=plt_title, sizing_mode="scale_width")
    plt.background_fill_color = "#f5f5f5"
    plt.grid.grid_line_color = "white"
    plt.xaxis.axis_label = 'Recall'
    plt.yaxis.axis_label = 'Precision'
    plt.axis.axis_line_color = None

    # draw curve data
    plt.line(x='rec', y='prec', line_width=2, color='#ebbd5b', source=source)
    plt.add_tools(bokeh.models.HoverTool(
      tooltips=[
        ( 'score', '@score{0.0000 a}'),
        ( 'Prec', '@prec'),
        ( 'Recall', '@rec'),
      ],
      formatters={
        'rec'      : 'printf',
        'prec' : 'printf',
      },
      mode='vline'
    ))
    bokeh_io.save(plt)
    return


def adjust_axes(r, t, fig, axes):

def draw_stacked_bars(df, perc, sample="", parts=[], outfile="", colours=COLOURS):
    """
    Takes a file of quantified read annotations by a pipeline (e.g. PZN)
    and draws a stacked bar chart using Bokeh, creating an interactive
    HTML file.
    
    Input:
     - A Pandas dataframe with all the required data
     - The column name of Sample IDs (as string)
     - The column names for the numbers ("parts", as list)
     - The title of the figure (as string)
     - The name of the output file (as string)
     - The colours to be used for the bars (as list of strings - hexcodes)
    Output:
     - Interactive (Bokeh) stacked bar chart that visualises
      the composition of each sample in your experiment
    """
    # Set a comfortable length depending on the number of samples to show:
    if len(df[sample]) > 40:
        width = len(df[sample]) * 20
    # Have the graph grow horizontally to accomodate large numbers of samples
    elif len(df[sample]) > 25:
        width = len(df[sample]) * 33
    else:
        # With a minimum of 800 to display the legend
        width = 800

    nr_fig = figure(
        x_range=df[sample],
        plot_height=800,
        plot_width=width,
        title="Composition of samples (read-based)",
        toolbar_location=None,
        tools="hover, pan",
        tooltips="@%s $name: @$name" % sample,
    )

    nr_fig.vbar_stack(
        parts,
        x=sample,
        width=0.9,
        color=colours,
        source=df,
        legend=[value(x) for x in parts],
    )

    nr_fig.y_range.start = 0
    nr_fig.x_range.range_padding = 0.1
    nr_fig.xgrid.grid_line_color = None
    nr_fig.axis.minor_tick_line_color = None
    nr_fig.outline_line_color = None
    nr_fig.legend.location = "top_left"
    nr_fig.legend.orientation = "horizontal"
    nr_fig.xaxis.major_label_orientation = 1

    nr_panel = Panel(child=nr_fig, title="Absolute number of reads")

    perc_fig = figure(
        x_range=perc[sample],
        plot_height=800,
        plot_width=width,
        title="Composition of samples (percentages)",
        toolbar_location=None,
        tools="hover, pan",
        tooltips="@%s $name: @$name" % sample,
    )

    perc_fig.vbar_stack(
        parts,
        x=sample,
        width=0.9,
        color=colours,
        source=perc,
        legend=[value(x) for x in parts],
    )

    perc_fig.y_range.start = 0
    perc_fig.x_range.range_padding = 0.1
    perc_fig.xgrid.grid_line_color = None
    perc_fig.axis.minor_tick_line_color = None
    perc_fig.outline_line_color = None
    perc_fig.legend.location = "top_left"
    perc_fig.legend.orientation = "horizontal"
    perc_fig.xaxis.major_label_orientation = 1

    perc_panel = Panel(child=perc_fig, title="Percentage of reads")

    tabs = Tabs(tabs=[nr_panel, perc_panel])

    output_file(outfile)
    save(tabs)

    return None


def main():

    def save(self, title=None):
        """save the json file.
        Parameters
        ----------
        title: string
            title of the HTML file
        """
        title = title or self.title
        if len(self.figures) > 0:
            if os.path.isfile(self.plot_path):
                os.remove(self.plot_path)
            if self.first_save:
                self.first_save = False
                logging.info('Plot file saved at: {}'.format(
                    os.path.abspath(self.plot_path)))

            output_file(self.plot_path, title=title)
            plot = column(
                Div(text='  <  h1 align="center">{} < /h1>'.format(title)), *self.figures)
            save(plot)
            self.clear()

        if self.data_format == 'json':
            self.results.to_json(self.data_path, orient='records', lines=True)
        else:
            self.results.to_csv(self.data_path, index=False, index_label=False)

    def load(self, path=None):

    def save(self, pm: PrettyMIDI, filepath: str):
        """
        Saves the pretty midi object as a plot file (html) in the provided file. If
        the live reload option is activated, the opened page will periodically
        refresh.

          :param pm: the PrettyMIDI instance to plot
          :param filepath: the file path to save the resulting plot to
          :return: the bokeh plot layout
        """
        plot = self.plot(pm)
        if self._live_reload:
            html = file_html(plot, CDN)
            html = html.replace("  <  /head>", """
               < script type="text/javascript">
                var liveReloadInterval = window.setInterval(function(){
                  location.reload();
                }, 2000);
               < /script>
               < /head>""")
            with open(filepath, 'w') as file:
                file.write(html)
        else:
            output_file(filepath)
            save(plot)
        return plot

    def show(self, pm: PrettyMIDI, filepath: str):

    def show(self, pm: PrettyMIDI, filepath: str):
        """
        Shows the pretty midi object as a plot file (html) in the browser. If
        the live reload option is activated, the opened page will periodically
        refresh.

          :param pm: the PrettyMIDI instance to plot
          :param filepath: the file path to save the resulting plot to
          :return: the bokeh plot layout
        """
        plot = self.plot(pm)
        if self._live_reload:
            html = file_html(plot, CDN)
            html = html.replace("  <  /head>", """
               < script type="text/javascript">
                var liveReloadInterval = window.setInterval(function(){
                  location.reload();
                }, 2000);
               < /script>
               < /head>""")
            with open(filepath, 'w') as file:
                file.write(html)
            if self._show_counter == 0:
                import webbrowser
                webbrowser.open("file://" + os.path.realpath(filepath), new=2)
        else:
            output_file(filepath)
            show(plot)
        self._show_counter += 1
        return plot

    def show_notebook(self, pm: PrettyMIDI):

def bokey_plot(dictionary_input, folder_bokey_default,mode="single", output=False, id_=None,info="",color_map_mode="continuous"):
    """
    as input dictionnary of the form
    - label : {x:, y: , label:[labels]}
    mode allows you to specify if you want everything on one single plot or if you want distinct plots
    """
    reset_output()
    source = {}
    i = 0
    color = ["blue","#ee6666"]
    assert color_map_mode in ["continuous","divergent"]
    if color_map_mode=="continuous": color_map= cm.OrRd
    elif color_map_mode=="divergent": color_map=cm.rainbow
    if id_ is None:
        id_ = str(uuid.uuid4())[0:8]

    if mode == "single":
        p = figure(plot_width=800, plot_height=1000, 
        		   tools=[BoxZoomTool(), ResetTool(), WheelZoomTool()],
                   toolbar_sticky=False, toolbar_location="right",
                   title='T-SNE '+info)  # x_range=Range1d(-6,6),y_range=Range1d(int(min(y))-1,int(max(y))+1))

    for key in dictionary_input.keys():
        if mode == "distinct":
            p = figure(plot_width=800, plot_height=1000,
                       title='T-SNE '+info)  # x_range=Range1d(-6,6),y_range=Range1d(int(min(y))-1,int(max(y))+1))

        source[key] = ColumnDataSource(data=dict(height=dictionary_input[key]["x"],
                                                 weight=dictionary_input[key]["y"],
                                                 names=dictionary_input[key]["label"]))

        colors = ["#%02x%02x%02x" %(int(r), int(g), int(b)) for r, g, b, _ in (255)*color_map(Normalize(vmin=0,vmax=5)(dictionary_input[key]["color"]))]
        colors_legend = ["#%02x%02x%02x" %(int(r), int(g), int(b)) for r, g, b, _ in (255)*color_map(Normalize(vmin=0,vmax=5)(np.sort(list(set(dictionary_input[key]["color"])))))]

        color_mapper = LinearColorMapper(palette=colors_legend)
        ticker = FixedTicker(ticks=[0, 1, 2, 3, 4,5])
        formatter = FuncTickFormatter(code="""
                                        function(tick) {
                                            data = {0: '0-10', 1: '10-20', 2: '20-30', 3: '30-40', 4: '40-50',50: '50plus'}
                                            return data[tick], " ,
                                        }
                                        """)

        cbar = ColorBar(color_mapper=color_mapper, ticker=ticker, formatter=formatter,
                        major_tick_out=0, major_tick_in=0, major_label_text_align='left',
                        major_label_text_font_size='100pt', label_standoff=5)

        p.scatter(x='weight', y='height', size=8, source=source[key], legend=key)# color=colors)
        p.add_layout(cbar)

        labels = LabelSet(x='weight', y='height', text='names', level='glyph',
                          x_offset=5, y_offset=5, source=source[key], render_mode='canvas')
        p.add_layout(labels)
        i += 1
        if output:
            output_file(folder_bokey_default+id_+"_tsne_"+key+"_"+info+"_bok.html")
            print(folder_bokey_default+id_+"_tsne_"+key+"_"+info+"_bok.html")
            show(p)
        #if mode == "distinct":
        #    output_notebook()
        #    show(p)
    if mode == "single":
        output_notebook()
        show(p)

    def savefig(self, file='MWPlot.html'):
        output_file(file)
        save(self.bokeh_fig)
        
class MWSkyMapBokeh(MWSkyMapMaster):

    def savefig(self, file='MWSkyMap.html'):
        output_file(file)
        save(self.bokeh_fig)

def visualize(
    k_mers,
    target_freqs,
    optimized_freqs,
    original_freqs=None,
    title="Freqgen Optimization Results",
    plot_height=400,
    plot_width=1200,
    show=True,
    filepath="freqgen.html",
    codons=False,
):
    """Creates a visualization of the results of a Freqgen optimization.

    Note:
        Currently, this function does not support the visualization of codon
        optimization *and* *k*-mer optimization simultaneously.

    Args:
        k_mers (list): A list of the *k*-mers to use as the labels for the *x*-axis.
        target_freqs (list): A list of the target frequencies in the same order as the `k_mers` argument.
        optimized_freqs (list): A list of the resultant frequencies in the same order as the `k_mers` argument.
        original_freqs (list, optional): A list of the original frequencies in the same order as the `k_mers` argument.
        title (str, optional): A title to use for the graph. Defaults to "Freqgen Optimization Results".
        plot_height (int, optional): The height for the graph. Defaults to 400.
        plot_width (int, optional): The width for the graph. Defaults to 1200.
        show (bool, optional): Whether to show the plot or simply return it. Defaults to True.
        filepath (str, optional): The output filepath. Defaults to "freqgen.html".
        codons (bool, optional): Whether codons are included in the input vectors. If they are, the *x*-axis legend will updated accordingly.

    Note:
        Codons must be denoted with a ``*`` in the ``k_mers`` argument. For example, the codon GAG should be passed as ``GAG*``

    Returns:
        bokeh.plotting.figure.Figure: A Bokeh figure containing the bar graph.
    """

    # validate that all the codons that should be there are there
    if codons and len([k_mer for k_mer in k_mers if k_mer.endswith("*")]) != 64:
        raise ValueError("You appear to be passing an incomplete list of codons.")

    codons_only = False
    if all([k_mer.endswith("*") for k_mer in k_mers]) and codons:
        k_mers = [k_mer[:-1] for k_mer in k_mers]
        codons_only = True

    output_file(filepath)

    categories = ["Original", "Target", "Optimized"]

    data = {"k_mers": k_mers, "Target": target_freqs, "Optimized": optimized_freqs}

    # adjust spacing depending on if there's three bars or two
    if not isinstance(original_freqs, type(None)):
        offset = 0.25
        data["Original"] = original_freqs
    else:
        offset = 0.15

    # identify the greatest y value for setting bounds
    y_max = max((max(target_freqs), max(optimized_freqs)))
    if not isinstance(original_freqs, type(None)):
        y_max = max((max(original_freqs), y_max))

    source = ColumnDataSource(data=data)

    p = figure(
        x_range=k_mers,
        plot_height=plot_height,
        plot_width=plot_width,
        title=title,
        y_range=(0, 1.2 * y_max),
    )

    if not isinstance(original_freqs, type(None)):
        p.vbar(
            x=dodge("k_mers", -0.25, range=p.x_range),
            top="Original",
            width=0.2,
            source=source,
            color=Set2_3[0],
            legend=value("Original"),
        )

    p.vbar(
        x=dodge(
            "k_mers",
            0.0 if not isinstance(original_freqs, type(None)) else -offset,
            range=p.x_range,
        ),
        top="Target",
        width=0.2,
        source=source,
        color=Set2_3[1],
        legend=value("Target"),
    )

    p.vbar(
        x=dodge("k_mers", offset, range=p.x_range),
        top="Optimized",
        width=0.2,
        source=source,
        color=Set2_3[2],
        legend=value("Optimized"),
    )

    p.x_range.range_padding = 0.05
    p.xgrid.grid_line_color = None
    p.legend.location = "top_right"
    p.legend.orientation = "horizontal"
    p.legend.click_policy = "hide"

    # decide the x-axis label
    if codons_only:
        p.xaxis.axis_label = "codon"
    elif codons:
        p.xaxis.axis_label = "k-mer (* denotes codon)"
    else:
        p.xaxis.axis_label = "k-mer"

    if len(k_mers) >= 32:
        p.xaxis.major_label_orientation = math.pi / 2

    p.yaxis.axis_label = "frequency"
    if show:
        _show(p)
    else:
        save(p, filename=filepath)
    return p

def _bokeh_timeseries_plots(varnames, time_units, var_units, phase_names, phases_node_path,
                            last_solution_case, last_simulation_case, plot_dir_path, num_cols=2,
                            bg_fill_color='#282828', grid_line_color='#666666', open_browser=False):
    from bokeh.io import output_notebook, output_file, save, show
    from bokeh.layouts import gridplot, column, row, grid, layout
    from bokeh.models import Legend, LegendItem
    from bokeh.plotting import figure
    import bokeh.palettes as bp

    if dymos_options['notebook_mode']:
        output_notebook()
    else:
        output_file(os.path.join(plot_dir_path, 'plots.html'))

    # Prune the edges from the color map
    cmap = bp.turbo(len(phase_names) + 2)[1:-1]
    figures = []
    colors = {}
    sol_plots = {}
    sim_plots = {}

    # Get the minimum and maximum times in any phase, so when we plot a variable that only exists
    # in a few phases, it is plotted against the entire time range.
    min_time = 1.0E21
    max_time = -1.0E21

    for iphase, phase_name in enumerate(phase_names):
        if phases_node_path:
            time_name = f'{phases_node_path}.{phase_name}.timeseries.time'
        else:
            time_name = f'{phase_name}.timeseries.time'
        min_time = min(min_time, np.min(last_solution_case.outputs[time_name]))
        max_time = max(max_time, np.max(last_solution_case.outputs[time_name]))
        colors[phase_name] = cmap[iphase]

    for ivar, var_name in enumerate(varnames):
        # Get the labels
        time_label = f'time ({time_units[var_name]})'
        var_label = f'{var_name} ({var_units[var_name]})'
        title = f'timeseries.{var_name}'

        # add labels, title, and legend
        padding = 0.05 * (max_time - min_time)
        fig = figure(title=title, background_fill_color=bg_fill_color,
                     x_range=(min_time - padding, max_time + padding), plot_width=180,
                     plot_height=180)
        fig.xaxis.axis_label = time_label
        fig.yaxis.axis_label = var_label
        fig.xgrid.grid_line_color = grid_line_color
        fig.ygrid.grid_line_color = grid_line_color

        # Plot each phase
        for iphase, phase_name in enumerate(phase_names):
            sol_color = cmap[iphase]
            sim_color = cmap[iphase]

            if phases_node_path:
                var_name_full = f'{phases_node_path}.{phase_name}.timeseries.{var_name}'
                time_name = f'{phases_node_path}.{phase_name}.timeseries.time'
            else:
                var_name_full = f'{phase_name}.timeseries.{var_name}'
                time_name = f'{phase_name}.timeseries.time'

            # Get values
            if var_name_full not in last_solution_case.outputs:
                continue

            var_val = last_solution_case.outputs[var_name_full]
            time_val = last_solution_case.outputs[time_name]

            for idxs, i in np.ndenumerate(np.zeros(var_val.shape[1:])):
                var_val_i = var_val[:, idxs].ravel()
                sol_plots[phase_name] = fig.circle(time_val.ravel(), var_val_i, size=5,
                                                   color=sol_color, name='sol:' + phase_name)

            # get simulation values, if plotting simulation
            if last_simulation_case:
                # if the phases_node_path is empty, need to pre-pend names with "sim_traj."
                #   as that is pre-pended in Trajectory.simulate code
                sim_prefix = '' if phases_node_path else 'sim_traj.'
                var_val_simulate = last_simulation_case.outputs[sim_prefix + var_name_full]
                time_val_simulate = last_simulation_case.outputs[sim_prefix + time_name]
                for idxs, i in np.ndenumerate(np.zeros(var_val_simulate.shape[1:])):
                    var_val_i = var_val_simulate[:, idxs].ravel()
                    sim_plots[phase_name] = fig.line(time_val_simulate.ravel(), var_val_i,
                                                     line_dash='solid', line_width=0.5, color=sim_color,
                                                     name='sim:' + phase_name)
        figures.append(fig)

    # Implement a single legend for all figures using the example here:
    # https://stackoverflow.com/a/56825812/754536

    # ## Use a dummy figure for the LEGEND
    dum_fig = figure(outline_line_alpha=0, toolbar_location=None,
                     background_fill_color=bg_fill_color, plot_width=250, max_width=250)

    # set the components of the figure invisible
    for fig_component in [dum_fig.grid, dum_fig.ygrid, dum_fig.xaxis, dum_fig.yaxis]:
        fig_component.visible = False

    # The glyphs referred by the legend need to be present in the figure that holds the legend,
    # so we must add them to the figure renderers.
    sol_legend_items = [(phase_name + ' solution', [dum_fig.circle([0], [0],
                                                                   size=5,
                                                                   color=colors[phase_name],
                                                                   tags=['sol:' + phase_name])]) for phase_name in phase_names]
    sim_legend_items = [(phase_name + ' simulation', [dum_fig.line([0], [0],
                                                                   line_dash='solid',
                                                                   line_width=0.5,
                                                                   color=colors[phase_name],
                                                                   tags=['sim:' + phase_name])]) for phase_name in phase_names]
    legend_items = [j for i in zip(sol_legend_items, sim_legend_items) for j in i]

    # # set the figure range outside of the range of all glyphs
    dum_fig.x_range.end = 1005
    dum_fig.x_range.start = 1000

    legend = Legend(click_policy='hide', location='top_left', border_line_alpha=0, items=legend_items,
                    background_fill_alpha=0.0, label_text_color='white', label_width=120, spacing=10)

    dum_fig.add_layout(legend, place='center')

    gd = gridplot(figures, ncols=num_cols, sizing_mode='scale_both')

    plots = gridplot([[gd, column(dum_fig, sizing_mode='stretch_height')]],
                     toolbar_location=None,
                     sizing_mode='scale_both')

    if dymos_options['notebook_mode'] or open_browser:
        show(plots)
    else:
        save(plots)


def timeseries_plots(solution_recorder_filename, simulation_record_file=None, plot_dir="plots"):

def categorical_scatter_jitter_SVF(df):
    from bokeh.io import show, output_file,export_svgs
    from bokeh.models import ColumnDataSource
    from bokeh.plotting import figure
    from bokeh.sampledata.commits import data
    from bokeh.transform import jitter
    from bokeh.io import output_notebook, show
    output_notebook()
    output_file("categorical_scatter_jitter_SVF.html")    
    
    # 'SVFep_min', 'SVFep_max', 'SVFep_mean','SVFep_count', 'SVFep_sum', 'SVFep_std', 'SVFep_median',
    # 'SVFep_majority', 'SVFep_minority', 'SVFep_unique', 'SVFep_range','SVFep_nodata',    
    
    # SVFFields=['SVF_mean','SVF_std', 'SVF_min', 'SVF_25%', 'SVF_50%', 'SVF_75%', 'SVF_max']
    SVFFields=['SVFep_min', 'SVFep_max', 'SVFep_mean','SVFep_std', 'SVFep_median','SVFep_majority', 'SVFep_minority', 'SVFep_range',]
    # SVFFields=['SVF_max']
    # DAYS = ['Sun', 'Sat', 'Fri', 'Thu', 'Wed', 'Tue', 'Mon']
    
    data_SVFFields=df[SVFFields]
    print(data_SVFFields.columns)
    
    dfStack=df[SVFFields].stack()
    dfStackDf=pd.DataFrame(dfStack).reset_index()
    dfStackDfAdj=dfStackDf[["level_1",0]].rename(columns={"level_1":"label",0:"value"})

    source=ColumnDataSource(dfStackDfAdj)
    
    p = figure(plot_width=2000, plot_height=500, y_range=SVFFields, x_axis_type='linear',title=" ") #Enum('linear', 'log', 'datetime', 'mercator') plot_width=2000, plot_height=500,1000/200

    p.circle(x='value', y=jitter('label', width=0.6, range=p.y_range),  source=source, alpha=0.3)        
    
    # p.xaxis.formatter.days = ['%Hh']
    p.x_range.range_padding = 0
    p.ygrid.grid_line_color = None
    
    # p.xaxis.axis_label="xaxis_name"
    # p.xaxis.axis_label_text_font_size = "25pt"
    # p.xaxis.axis_label_text_font = "SimHei"
    # p.xaxis.axis_label_text_color = "black"
    # p.axis.minor_tick_in = -3
    # p.axis.minor_tick_out = 6
    p.xaxis.major_tick_line_width = 3
    p.xaxis.major_label_text_font_size='26pt'
    p.yaxis.major_label_text_font_size='26pt'
        
    show(p)
    p.output_backend = "svg"
    export_svgs(p, filename=r"C:\Users\richi\omen-richiebao\omen_sf_paper_2020\01_ spatial structure of parks of the city of Chicago\fig\svfPlot_.svg")

    return dfStackDfAdj

'''

def predvsobsplot(Y,Yhat,*,variable_names=False,obs_names=False):
    """
    Plot observed vs predicted values
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    Y:    Numpy array or pandas dataframe with observed values
    
    Yhat: Numpy array with predicted data
    
    variable_names: List with names of the variables  (columns of Y) if Y
    is a numpy array [otherwise names are taken from Pandasdataframe]
    
    """


    if isinstance(Y,np.ndarray):
        if isinstance(variable_names,np.bool):
            YVar = []
            for n in list(np.arange(Y.shape[1])+1):
                YVar.append('YVar #'+str(n))
        else:
            YVar=variable_names
            
        if isinstance(obs_names,np.bool):
            ObsID_ = []
            for n in list(np.arange(Y.shape[0])+1):
                ObsID_.append('Obs #'+str(n))
        else:
            ObsID_=obs_names

    elif isinstance(Y,pd.DataFrame):
        Y_=np.array(Y.values[:,1:]).astype(float)
        YVar = Y_.columns.values
        YVar = YVar[1:]
        YVar = YVar.tolist()
        ObsID_ = Y.values[:,0].astype(str)
        ObsID_ = ObsID_.tolist()
        
    TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
    TOOLTIPS = [
                ("index", "$index"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID")
                ]
    
    rnd_num=str(int(np.round(1000*np.random.random_sample())))
    output_file("ObsvsPred_"+rnd_num+".html",title='ObsvsPred')
    plot_counter=0
    
    for i in list(range(Y_.shape[1])):
        x_ = Y[:,i]
        y_ = Yhat[:,i]           
        source = ColumnDataSource(data=dict(x=x_, y=y_,ObsID=ObsID_))
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=600, plot_height=600, title=YVar[i])
        p.circle('x', 'y', source=source,size=7,color='darkblue')
        p.line([np.nanmin(x_),np.nanmax(x_)],[np.nanmin(y_),np.nanmax(y_)],line_color='cyan',line_dash='dashed')
        p.xaxis.axis_label ='Observed'
        p.yaxis.axis_label ='Predicted'
        if plot_counter==0:
            p_list=[p]
        else:
            p_list.append(p)
        plot_counter = plot_counter+1            
    show(column(p_list))        
    return    

def conv_pls_2_eiot(plsobj,*,r_length=False):

def r2pv(mvmobj,*,plotwidth=600,plotheight=400):
    """
    R2 per variable plots
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pca or phi.pls
    """
    A= mvmobj['T'].shape[1]
    num_varX=mvmobj['P'].shape[0]
    
    if 'Q' in mvmobj:
        is_pls=True
        lv_prefix='LV #'
    else:
        is_pls=False
        lv_prefix='PC #'
        
    lv_labels = []   
    for a in list(np.arange(A)+1):
        lv_labels.append(lv_prefix+str(a))    
    

    if 'varidX' in mvmobj:
        r2pvX_dict = {'XVar': mvmobj['varidX']}
        XVar=mvmobj['varidX']
    else:
        XVar = []
        for n in list(np.arange(num_varX)+1):
            XVar.append('XVar #'+str(n))               
        r2pvX_dict = {'XVar': XVar}
        
    for i in list(np.arange(A)):
        r2pvX_dict.update({lv_labels[i] : mvmobj['r2xpv'][:,i].tolist()})
        
    if 'Q' in mvmobj:
        num_varY=mvmobj['Q'].shape[0]
        if 'varidY' in mvmobj:
            r2pvY_dict = {'YVar': mvmobj['varidY']}
            YVar=mvmobj['varidY']
        else:
            YVar = []
            for n in list(np.arange(num_varY)+1):
                YVar.append('YVar #'+str(n))               
            r2pvY_dict = {'YVar': YVar}
        for i in list(np.arange(A)):
            r2pvY_dict.update({lv_labels[i] : mvmobj['r2ypv'][:,i].tolist()})
            
    
    if is_pls:
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("r2xypv_"+rnd_num+".html",title="R2XYPV") 
        colormap =cm.get_cmap("rainbow")
        different_colors=A
        color_mapping=colormap(np.linspace(0,1,different_colors),1,True)
        bokeh_palette=["#%02x%02x%02x" % (r, g, b) for r, g, b in color_mapping[:,0:3]]  
               
        px = figure(x_range=XVar, title="R2X Per Variable",
             tools="save,box_zoom,hover,reset", tooltips="$name @XVar: @$name",plot_width=plotwidth,plot_height=plotheight)
        
        px.vbar_stack(lv_labels, x='XVar', width=0.9,color=bokeh_palette,source=r2pvX_dict)
        px.y_range.range_padding = 0.1
        px.ygrid.grid_line_color = None
        px.xgrid.grid_line_color = None
        px.axis.minor_tick_line_color = None
        px.outline_line_color = None
        px.yaxis.axis_label = 'R2X'
        px.xaxis.major_label_orientation = 45
        py = figure(x_range=YVar, plot_height=plotheight, title="R2Y Per Variable",
            tools="save,box_zoom,hover,reset", tooltips="$name @YVar: @$name",plot_width=plotwidth)
        
        py.vbar_stack(lv_labels, x='YVar', width=0.9,color=bokeh_palette,source=r2pvY_dict)
        py.y_range.range_padding = 0.1
        py.ygrid.grid_line_color = None
        py.axis.minor_tick_line_color = None
        py.xgrid.grid_line_color = None
        py.outline_line_color = None
        py.yaxis.axis_label = 'R2Y'
        py.xaxis.major_label_orientation = 45
        show(column(px,py))
        
        
    else:   
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("r2xpv_"+rnd_num+".html",title='R2XPV') 
        colormap =cm.get_cmap("rainbow")
        different_colors=A
        color_mapping=colormap(np.linspace(0,1,different_colors),1,True)
        bokeh_palette=["#%02x%02x%02x" % (r, g, b) for r, g, b in color_mapping[:,0:3]]  
               
        p = figure(x_range=XVar, title="R2X Per Variable",
             tools="save,box_zoom,hover,reset", tooltips="$name @XVar: @$name",plot_width=plotwidth,plot_height=plotheight)
        
        p.vbar_stack(lv_labels, x='XVar', width=0.9,color=bokeh_palette,source=r2pvX_dict)
        p.y_range.range_padding = 0.1
        p.ygrid.grid_line_color = None
        p.axis.minor_tick_line_color = None
        p.outline_line_color = None
        p.yaxis.axis_label = 'R2X'
        p.xaxis.major_label_orientation = 45
        show(p)
    return
    
def loadings(mvmobj,*,plotwidth=600,xgrid=False):

def loadings(mvmobj,*,plotwidth=600,xgrid=False):
    """
    Column plots of loadings
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pca or phi.pls
    
    """
    
    A= mvmobj['T'].shape[1]
    num_varX=mvmobj['P'].shape[0]    
    if 'Q' in mvmobj:
        is_pls=True
        lv_prefix='LV #'
    else:
        is_pls=False
        lv_prefix='PC #'       
    lv_labels = []   
    for a in list(np.arange(A)+1):
        lv_labels.append(lv_prefix+str(a))    
    if 'varidX' in mvmobj:
        X_loading_dict = {'XVar': mvmobj['varidX']}
        XVar=mvmobj['varidX']
    else:
        XVar = []
        for n in list(np.arange(num_varX)+1):
            XVar.append('XVar #'+str(n))               
        X_loading_dict = {'XVar': XVar}
    if 'Q' in mvmobj:
        for i in list(np.arange(A)):
            X_loading_dict.update({lv_labels[i] : mvmobj['Ws'][:,i].tolist()})
            
        num_varY=mvmobj['Q'].shape[0]
        if 'varidY' in mvmobj:
            Q_dict = {'YVar': mvmobj['varidY']}
            YVar=mvmobj['varidY']
        else:
            YVar = []
            for n in list(np.arange(num_varY)+1):
                YVar.append('YVar #'+str(n))               
            Q_dict = {'YVar': YVar}
        for i in list(np.arange(A)):
            Q_dict.update({lv_labels[i] : mvmobj['Q'][:,i].tolist()})
    else:
        for i in list(np.arange(A)):
            X_loading_dict.update({lv_labels[i] : mvmobj['P'][:,i].tolist()})
            
    TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
    TOOLTIPS = [
                ("Variable:","@names")
                ]
      
    if is_pls:
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings X Space_"+rnd_num+".html",title='X Loadings PLS')
        for i in list(np.arange(A)):
            p = figure(x_range=XVar, title="X Space Loadings "+lv_labels[i],
                    tools=TOOLS,tooltips=TOOLTIPS,plot_width=plotwidth)
            source1 = ColumnDataSource(data=dict(x_=XVar, y_=mvmobj['Ws'][:,i].tolist(),names=XVar)) 
            
            #p.vbar(x=XVar, top=mvmobj['Ws'][:,i].tolist(), width=0.5)
            p.vbar(x='x_', top='y_', source=source1,width=0.5)
            p.ygrid.grid_line_color = None    
            if xgrid:
                p.xgrid.grid_line_color = 'lightgray'
                
            else:
                p.xgrid.grid_line_color = None    
                
            p.yaxis.axis_label = 'W* ['+str(i+1)+']'
            hline = Span(location=0, dimension='width', line_color='black', line_width=2)
            p.renderers.extend([hline])
            p.xaxis.major_label_orientation = 45
            if i==0:
                p_list=[p]
            else:
                p_list.append(p)
        show(column(p_list))    
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings Y Space_"+rnd_num+".html",title='Y Loadings PLS')
        for i in list(np.arange(A)):
            p = figure(x_range=YVar, title="Y Space Loadings "+lv_labels[i],
                    tools="save,box_zoom,pan,reset",tooltips=TOOLTIPS,plot_width=plotwidth)
            
            source1 = ColumnDataSource(data=dict(x_=YVar, y_=mvmobj['Q'][:,i].tolist(),names=YVar)) 
            #p.vbar(x=YVar, top=mvmobj['Q'][:,i].tolist(), width=0.5)
            p.vbar(x='x_', top='y_', source=source1,width=0.5)
            p.ygrid.grid_line_color = None    
            if xgrid:
                p.xgrid.grid_line_color = 'lightgray'
            else:
                p.xgrid.grid_line_color = None    
            p.yaxis.axis_label = 'Q ['+str(i+1)+']'
            hline = Span(location=0, dimension='width', line_color='black', line_width=2)
            p.renderers.extend([hline])
            p.xaxis.major_label_orientation = 45
            if i==0:
                p_list=[p]
            else:
                p_list.append(p)                    
        show(column(p_list))
    else:   
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings X Space_"+rnd_num+".html",title='X Loadings PCA') 
        for i in list(np.arange(A)):
            source1 = ColumnDataSource(data=dict(x_=XVar, y_=mvmobj['P'][:,i].tolist(),names=XVar))  
            
            p = figure(x_range=XVar, title="X Space Loadings "+lv_labels[i],
                    tools=TOOLS,tooltips=TOOLTIPS,plot_width=plotwidth)
            
            #p.vbar(x=XVar, top=mvmobj['P'][:,i].tolist(), width=0.5)
            
            p.vbar(x='x_', top='y_', source=source1,width=0.5)
            if xgrid:
                p.xgrid.grid_line_color = 'lightgray'
            else:
                p.xgrid.grid_line_color = None    
            p.yaxis.axis_label = 'P ['+str(i+1)+']'
            hline = Span(location=0, dimension='width', line_color='black', line_width=2)
            p.renderers.extend([hline])
            p.xaxis.major_label_orientation = 45
            if i==0:
                p_list=[p]
            else:
                p_list.append(p)
        show(column(p_list))
    return    

def loadings_map(mvmobj,dims,*,plotwidth=600):

def loadings_map(mvmobj,dims,*,plotwidth=600):
    """
    Scatter plot overlaying X and Y loadings 
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pca or phi.pls
    dims: what latent spaces to plot in x and y axes e.g. dims=[1,2]
    """
    A= mvmobj['T'].shape[1]
    num_varX=mvmobj['P'].shape[0]    
    if 'Q' in mvmobj:
        lv_prefix='LV #'     
        lv_labels = []   
        for a in list(np.arange(A)+1):
            lv_labels.append(lv_prefix+str(a))    
        if 'varidX' in mvmobj:
            XVar=mvmobj['varidX']
        else:
            XVar = []
            for n in list(np.arange(num_varX)+1):
                XVar.append('XVar #'+str(n))               
        num_varY=mvmobj['Q'].shape[0]
        if 'varidY' in mvmobj:
            YVar=mvmobj['varidY']
        else:
            YVar = []
            for n in list(np.arange(num_varY)+1):
                YVar.append('YVar #'+str(n))               
    
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings Map"+rnd_num+".html",title='Loadings Map')
       
    
        x_ws = mvmobj['Ws'][:,dims[0]-1]
        x_ws = x_ws/np.max(np.abs(x_ws))
        y_ws = mvmobj['Ws'][:,dims[1]-1]
        y_ws = y_ws/np.max(np.abs(y_ws))
        
        x_q = mvmobj['Q'][:,dims[0]-1]
        x_q = x_q/np.max(np.abs(x_q))   
        y_q = mvmobj['Q'][:,dims[1]-1]
        y_q = y_q/np.max(np.abs(y_q))
        
        
        TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
        TOOLTIPS = [
                ("index", "$index"),
                ("(x,y)", "($x, $y)"),
                ("Variable:","@names")
                ]
    
        source1 = ColumnDataSource(data=dict(x=x_ws, y=y_ws,names=XVar))  
        source2 = ColumnDataSource(data=dict(x=x_q, y=y_q,names=YVar)) 
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=plotwidth, title="Loadings Map LV["+str(dims[0])+"] - LV["+str(dims[1])+"]",
                                                                                                          x_range=(-1.5,1.5),y_range=(-1.5,1.5))
        p.circle('x', 'y', source=source1,size=10,color='darkblue')
        p.circle('x', 'y', source=source2,size=10,color='red')
        p.xaxis.axis_label = lv_labels [dims[0]-1]
        p.yaxis.axis_label = lv_labels [dims[1]-1]
        
        labelsX = LabelSet(x='x', y='y', text='names', level='glyph',x_offset=5, y_offset=5, source=source1, render_mode='canvas',text_color='darkgray')
        labelsY = LabelSet(x='x', y='y', text='names', level='glyph',x_offset=5, y_offset=5, source=source2, render_mode='canvas',text_color='darkgray')
        p.add_layout(labelsX)
        p.add_layout(labelsY)

        vline = Span(location=0, dimension='height', line_color='black', line_width=2)
        # Horizontal line
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        p.renderers.extend([vline, hline])
        show(p)    
    else:
        lv_prefix='PC #'     
        lv_labels = []   
        for a in list(np.arange(A)+1):
            lv_labels.append(lv_prefix+str(a))    
        if 'varidX' in mvmobj:
            XVar=mvmobj['varidX']
        else:
            XVar = []
            for n in list(np.arange(num_varX)+1):
                XVar.append('XVar #'+str(n))                   
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings Map"+rnd_num+".html",title='Loadings Map')    
        x_p = mvmobj['P'][:,dims[0]-1]
        y_p = mvmobj['P'][:,dims[1]-1]                        
        TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
        TOOLTIPS = [
                ("index", "$index"),
                ("(x,y)", "($x, $y)"),
                ("Variable:","@names")
                ]
    
        source1 = ColumnDataSource(data=dict(x=x_p, y=y_p,names=XVar))  
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=plotwidth, title="Loadings Map PC["+str(dims[0])+"] - PC["+str(dims[1])+"]",                                                                                                         x_range=(-1.5,1.5),y_range=(-1.5,1.5))
        p.circle('x', 'y', source=source1,size=10,color='darkblue')
        p.xaxis.axis_label = lv_labels [dims[0]-1]
        p.yaxis.axis_label = lv_labels [dims[1]-1]        
        labelsX = LabelSet(x='x', y='y', text='names', level='glyph',x_offset=5, y_offset=5, source=source1, render_mode='canvas',text_color='darkgray')
        p.add_layout(labelsX)
        vline = Span(location=0, dimension='height', line_color='black', line_width=2)
        # Horizontal line
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        p.renderers.extend([vline, hline])
        show(p)            
    return  

def weighted_loadings(mvmobj,*,plotwidth=600,xgrid=False):

def weighted_loadings(mvmobj,*,plotwidth=600,xgrid=False):
    """
    Column plots of loadings weighted by r2x/r2y correspondingly
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pca or phi.pls
    
    """
    A= mvmobj['T'].shape[1]
    num_varX=mvmobj['P'].shape[0]    
    if 'Q' in mvmobj:
        is_pls=True
        lv_prefix='LV #'
    else:
        is_pls=False
        lv_prefix='PC #'       
    lv_labels = []   
    for a in list(np.arange(A)+1):
        lv_labels.append(lv_prefix+str(a))    
    if 'varidX' in mvmobj:
        X_loading_dict = {'XVar': mvmobj['varidX']}
        XVar=mvmobj['varidX']
    else:
        XVar = []
        for n in list(np.arange(num_varX)+1):
            XVar.append('XVar #'+str(n))               
        X_loading_dict = {'XVar': XVar}
    if 'Q' in mvmobj:
        for i in list(np.arange(A)):
            X_loading_dict.update({lv_labels[i] : mvmobj['Ws'][:,i].tolist()})
            
        num_varY=mvmobj['Q'].shape[0]
        if 'varidY' in mvmobj:
            Q_dict = {'YVar': mvmobj['varidY']}
            YVar=mvmobj['varidY']
        else:
            YVar = []
            for n in list(np.arange(num_varY)+1):
                YVar.append('YVar #'+str(n))               
            Q_dict = {'YVar': YVar}
        for i in list(np.arange(A)):
            Q_dict.update({lv_labels[i] : mvmobj['Q'][:,i].tolist()})
    else:
        for i in list(np.arange(A)):
            X_loading_dict.update({lv_labels[i] : mvmobj['P'][:,i].tolist()})
            
    TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
    TOOLTIPS = [
                ("Variable:","@names")
                ]
    
    if is_pls:
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings X Space_"+rnd_num+".html",title='X Weighted Loadings PLS')
        for i in list(np.arange(A)):
            p = figure(x_range=XVar, title="X Space Weighted Loadings "+lv_labels[i],
                     tools=TOOLS,tooltips=TOOLTIPS,plot_width=plotwidth)
            source1 = ColumnDataSource(data=dict(x_=XVar, y_=(mvmobj['r2xpv'][:,i] * mvmobj['Ws'][:,i]).tolist(),names=XVar)) 
             
            #p.vbar(x=XVar, top=(mvmobj['r2xpv'][:,i] * mvmobj['Ws'][:,i]).tolist(), width=0.5)
            p.vbar(x='x_', top='y_', source=source1,width=0.5)
            p.ygrid.grid_line_color = None    
            if xgrid:
                p.xgrid.grid_line_color = 'lightgray'
            else:
                p.xgrid.grid_line_color = None    

            p.yaxis.axis_label = 'W* ['+str(i+1)+']'
            hline = Span(location=0, dimension='width', line_color='black', line_width=2)
            p.renderers.extend([hline])
            p.xaxis.major_label_orientation = 45
            if i==0:
                p_list=[p]
            else:
                p_list.append(p)
        show(column(p_list)) 
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings Y Space_"+rnd_num+".html",title='Y Weighted Loadings PLS')
        for i in list(np.arange(A)):
            p = figure(x_range=YVar, title="Y Space Weighted Loadings "+lv_labels[i],
                     tools=TOOLS,tooltips=TOOLTIPS,plot_width=plotwidth)
            source1 = ColumnDataSource(data=dict(x_=YVar, y_=(mvmobj['r2ypv'][:,i] * mvmobj['Q'][:,i]).tolist(),names=YVar)) 
            
            #p.vbar(x=YVar, top=(mvmobj['r2ypv'][:,i] * mvmobj['Q'][:,i]).tolist(), width=0.5)
            p.vbar(x='x_', top='y_', source=source1,width=0.5)
            p.ygrid.grid_line_color = None    
            if xgrid:
                p.xgrid.grid_line_color = 'lightgray'
            else:
                p.xgrid.grid_line_color = None    
            p.yaxis.axis_label = 'Q ['+str(i+1)+']'
            hline = Span(location=0, dimension='width', line_color='black', line_width=2)
            p.renderers.extend([hline])
            p.xaxis.major_label_orientation = 45
            if i==0:
                p_list=[p]
            else:
                p_list.append(p)                    
        show(column(p_list))
    else:   
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Loadings X Space_"+rnd_num+".html",title='X Weighted Loadings PCA') 
        for i in list(np.arange(A)):
            p = figure(x_range=XVar, title="X Space Weighted Loadings "+lv_labels[i],
                     tools=TOOLS,tooltips=TOOLTIPS,plot_width=plotwidth)
            source1 = ColumnDataSource(data=dict(x_=XVar, y_=(mvmobj['r2xpv'][:,i] * mvmobj['P'][:,i]).tolist(),names=XVar)) 
            
            #p.vbar(x=XVar, top=(mvmobj['r2xpv'][:,i] * mvmobj['P'][:,i]).tolist(), width=0.5)
            p.vbar(x='x_', top='y_', source=source1,width=0.5)
            p.ygrid.grid_line_color = None    
            if xgrid:
                p.xgrid.grid_line_color = 'lightgray'
            else:
                p.xgrid.grid_line_color = None    

            p.yaxis.axis_label = 'P ['+str(i+1)+']'
            hline = Span(location=0, dimension='width', line_color='black', line_width=2)
            p.renderers.extend([hline])
            p.xaxis.major_label_orientation = 45
            if i==0:
                p_list=[p]
            else:
                p_list.append(p)
        show(column(p_list))
    return  
 
def vip(mvmobj,*,plotwidth=600):

def vip(mvmobj,*,plotwidth=600):
    """
    Very Important to the Projection (VIP) plot
        by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pls
    """
    
    if 'Q' in mvmobj:  
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("VIP_"+rnd_num+".html",title='VIP Coefficient') 
        num_varX=mvmobj['P'].shape[0]            
        if 'varidX' in mvmobj:
            XVar=mvmobj['varidX']
        else:
            XVar = []
            for n in list(np.arange(num_varX)+1):
                XVar.append('XVar #'+str(n))               
        
        vip=np.sum(np.abs(mvmobj['Ws'] * np.tile(mvmobj['r2y'],(mvmobj['Ws'].shape[0],1)) ),axis=1)
        vip=np.reshape(vip,(len(vip),-1))
        sort_indx=np.argsort(-vip,axis=0)
        vip=vip[sort_indx]
        sorted_XVar=[]
        for i in sort_indx[:,0]:
            sorted_XVar.append(XVar[i])  
        p = figure(x_range=sorted_XVar, title="VIP",
            tools="save,box_zoom,pan,reset",plot_width=plotwidth)
        p.vbar(x=sorted_XVar, top=vip.tolist(), width=0.5)
        p.xgrid.grid_line_color = None
        p.yaxis.axis_label = 'Very Important to the Projection'
        p.xaxis.major_label_orientation = 45
        show(p)
    return    

def score_scatter(mvmobj,xydim,*,CLASSID=False,colorby=False,Xnew=False,add_ci=False,add_labels=False,add_legend=True,plotwidth=600,plotheight=600):

def score_scatter(mvmobj,xydim,*,CLASSID=False,colorby=False,Xnew=False,add_ci=False,add_labels=False,add_legend=True,plotwidth=600,plotheight=600):
    '''
    Score scatter plot
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj     : PLS or PCA object from phyphi
    xydim      : LV to plot on x and y axes. eg [1,2] will plot t1 vs t2
    CLASSID    : Pandas DataFrame with CLASSIDS
    colorby    : Category (one of the CLASSIDS) to color by
    Xnew       : New data for which to make the score plot this routine evaluates and plots
    add_ci     : when = True will add confidence intervals
    add_labels : When = True labels each point with Obs ID
    plotwidth  : If omitted, width is 600
    '''
    
    if isinstance(Xnew,bool):
        if 'obsidX' in mvmobj:
            ObsID_=mvmobj['obsidX']
        else:
            ObsID_ = []
            for n in list(np.arange(mvmobj['T'].shape[0])+1):
                ObsID_.append('Obs #'+str(n))  
        T_matrix=mvmobj['T']
    else:
        if isinstance(Xnew,np.ndarray):
            X_=Xnew.copy()
            ObsID_ = []
            for n in list(np.arange(Xnew.shape[0])+1):
                ObsID_.append('Obs #'+str(n))  
        elif isinstance(Xnew,pd.DataFrame):
            X_=np.array(Xnew.values[:,1:]).astype(float)
            ObsID_ = Xnew.values[:,0].astype(str)
            ObsID_ = ObsID_.tolist()
            
        if 'Q' in mvmobj:  
            xpred=phi.pls_pred(X_,mvmobj)
        else:
            xpred=phi.pca_pred(X_,mvmobj)
        T_matrix=xpred['Tnew']
        
    ObsNum_=[]    
    for n in list(range(1,len(ObsID_)+1)):
                ObsNum_.append('Obs #'+str(n))  
    
    if isinstance(CLASSID,np.bool): # No CLASSIDS
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Score_Scatter_"+rnd_num+".html",title='Score Scatter t['+str(xydim[0])+'] - t['+str(xydim[1])+ ']')

        x_=T_matrix[:,[xydim[0]-1]]
        y_=T_matrix[:,[xydim[1]-1]]

           
        source = ColumnDataSource(data=dict(x=x_, y=y_,ObsID=ObsID_,ObsNum=ObsNum_))
        TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
        TOOLTIPS = [
                ("Obs #", "@ObsNum"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID")
                ]
        
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=plotwidth,plot_height=plotheight, title='Score Scatter t['+str(xydim[0])+'] - t['+str(xydim[1])+ ']')
        p.circle('x', 'y', source=source,size=7)
        if add_ci:
            T_aux1=mvmobj['T'][:,[xydim[0]-1]]
            T_aux2=mvmobj['T'][:,[xydim[1]-1]]
            T_aux = np.hstack((T_aux1,T_aux2))
            st=(T_aux.T @ T_aux)/T_aux.shape[0]
            [xd95,xd99,yd95p,yd95n,yd99p,yd99n]=phi.scores_conf_int_calc(st,mvmobj['T'].shape[0])
            p.line(xd95,yd95p,line_color="gold",line_dash='dashed')
            p.line(xd95,yd95n,line_color="gold",line_dash='dashed')
            p.line(xd99,yd99p,line_color="red",line_dash='dashed')
            p.line(xd99,yd99n,line_color="red",line_dash='dashed')
            
        if add_labels:
            labelsX = LabelSet(x='x', y='y', text='ObsID', level='glyph',x_offset=5, y_offset=5, source=source, render_mode='canvas')
            p.add_layout(labelsX)
        p.xaxis.axis_label = 't ['+str(xydim[0])+']'
        p.yaxis.axis_label = 't ['+str(xydim[1])+']'
        # Vertical line
        vline = Span(location=0, dimension='height', line_color='black', line_width=2)
        # Horizontal line
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        p.renderers.extend([vline, hline])
        show(p)      
    else: # YES CLASSIDS
    
        Classes_=np.unique(CLASSID[colorby]).tolist()        
        
        A=len(Classes_)
        colormap =cm.get_cmap("rainbow")
        different_colors=A
        color_mapping=colormap(np.linspace(0,1,different_colors),1,True)
        bokeh_palette=["#%02x%02x%02x" % (r, g, b) for r, g, b in color_mapping[:,0:3]]  
        rnd_num=str(int(np.round(1000*np.random.random_sample())))               
        output_file("Score_Scatter_"+rnd_num+".html",title='Score Scatter t['+str(xydim[0])+'] - t['+str(xydim[1])+ ']') 
        x_=T_matrix[:,[xydim[0]-1]]
        y_=T_matrix[:,[xydim[1]-1]]          
        
        TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
        TOOLTIPS = [
                ("Obs #", "@ObsNum"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID"),
                ("Class:","@Class")
                ]        
        classid_=list(CLASSID[colorby])
        legend_it = []
        
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,toolbar_location="above",plot_width=plotwidth,plot_height=plotheight,title='Score Scatter t['+str(xydim[0])+'] - t['+str(xydim[1])+ ']')

        for classid_in_turn in Classes_:                      
            x_aux       = []
            y_aux       = []
            obsid_aux   = []
            obsnum_aux  = []
            classid_aux = []
            
            for i in list(range(len(ObsID_))):
                
                if classid_[i]==classid_in_turn:
                    x_aux.append(x_[i][0])
                    y_aux.append(y_[i][0])
                    obsid_aux.append(ObsID_[i])
                    obsnum_aux.append(ObsNum_[i])
                    classid_aux.append(classid_in_turn)
            source = ColumnDataSource(data=dict(x=x_aux, y=y_aux,ObsID=obsid_aux,ObsNum=obsnum_aux, Class=classid_aux))        
            color_=bokeh_palette[Classes_.index(classid_in_turn)]
            if add_legend:
                c = p.circle('x','y',source=source,color=color_)
                aux_=classid_in_turn
                if isinstance(aux_,(float,int)):
                    aux_=str(aux_)
                #legend_it.append((classid_in_turn, [c]))
                legend_it.append((aux_, [c]))
            else:
                p.circle('x','y',source=source,color=color_)
            if add_labels:
                labelsX = LabelSet(x='x', y='y', text='ObsID', level='glyph',x_offset=5, y_offset=5, source=source, render_mode='canvas')
                p.add_layout(labelsX)
        if add_ci:
            T_aux1=mvmobj['T'][:,[xydim[0]-1]]
            T_aux2=mvmobj['T'][:,[xydim[1]-1]]
            T_aux = np.hstack((T_aux1,T_aux2))
            st=(T_aux.T @ T_aux)/T_aux.shape[0]
            [xd95,xd99,yd95p,yd95n,yd99p,yd99n]=phi.scores_conf_int_calc(st,mvmobj['T'].shape[0])
            p.line(xd95,yd95p,line_color="gold",line_dash='dashed')
            p.line(xd95,yd95n,line_color="gold",line_dash='dashed')
            p.line(xd99,yd99p,line_color="red",line_dash='dashed')
            p.line(xd99,yd99n,line_color="red",line_dash='dashed') 
        p.xaxis.axis_label = 't ['+str(xydim[0])+']'
        p.yaxis.axis_label = 't ['+str(xydim[1])+']'
        # Vertical line
        vline = Span(location=0, dimension='height', line_color='black', line_width=2)
        # Horizontal line
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        p.renderers.extend([vline, hline])
        if add_legend:
            legend = Legend(items=legend_it, location='top_right')
            p.add_layout(legend, 'right')
            legend.click_policy="hide"
        show(p)
    return    

def score_line(mvmobj,dim,*,CLASSID=False,colorby=False,Xnew=False,add_ci=False,add_labels=False,add_legend=True,plotline=True,plotwidth=600,plotheight=600):

def score_line(mvmobj,dim,*,CLASSID=False,colorby=False,Xnew=False,add_ci=False,add_labels=False,add_legend=True,plotline=True,plotwidth=600,plotheight=600):
    '''
    Score scatter plot
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj     : PLS or PCA object from phyphi
    dim        : LV to plot eg "1" will plot t1 vs observation #
    CLASSID    : Pandas DataFrame with CLASSIDS
    colorby    : Category (one of the CLASSIDS) to color by
    Xnew       : New data for which to make the score plot this routine evaluates and plots
    add_ci     : When = True will add confidence intervals
    add_labels : When =True will display Obs ID per point
    plotwidth  : When Omitted is = 600
    plotline   : Adds a conecting line between dots [True by default]
    '''
    if not(isinstance(dim,list)):
        if isinstance(dim, int):
            dim=[dim]
  
    if isinstance(Xnew,bool):
        if 'obsidX' in mvmobj:
            ObsID_=mvmobj['obsidX']
        else:
            ObsID_ = []
            for n in list(np.arange(mvmobj['T'].shape[0])+1):
                ObsID_.append('Obs #'+str(n))  
        T_matrix=mvmobj['T']
    else:
        if isinstance(Xnew,np.ndarray):
            X_=Xnew.copy()
            ObsID_ = []
            for n in list(np.arange(Xnew.shape[0])+1):
                ObsID_.append('Obs #'+str(n))  
        elif isinstance(Xnew,pd.DataFrame):
            X_=np.array(Xnew.values[:,1:]).astype(float)
            ObsID_ = Xnew.values[:,0].astype(str)
            ObsID_ = ObsID_.tolist()
            
        if 'Q' in mvmobj:  
            xpred=phi.pls_pred(X_,mvmobj)
        else:
            xpred=phi.pca_pred(X_,mvmobj)
        T_matrix=xpred['Tnew']

    ObsNum_=[]    
    for n in list(range(1,len(ObsID_)+1)):
        ObsNum_.append('Obs #'+str(n))  
                       
    if isinstance(CLASSID,np.bool): # No CLASSIDS
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("Score_Line_"+rnd_num+".html",title='Score Line t['+str(dim[0])+ ']')

        y_=T_matrix[:,[dim[0]-1]]
        x_=list(range(1,y_.shape[0]+1))

           
        source = ColumnDataSource(data=dict(x=x_, y=y_,ObsID=ObsID_,ObsNum=ObsNum_))
        TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
        TOOLTIPS = [
                ("Obs#", "@ObsNum"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID")
                ]
        
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=plotwidth,plot_height=plotheight, title='Score Line t['+str(dim[0])+']' )
        p.circle('x', 'y', source=source,size=7)
        if plotline:
            p.line('x', 'y', source=source)
        if add_ci:
            lim95,lim99=phi.single_score_conf_int(mvmobj['T'][:,[dim[0]-1]])
            p.line(x_, lim95,line_color="gold",line_dash='dashed')
            p.line(x_,-lim95,line_color="gold",line_dash='dashed')
            p.line(x_, lim99,line_color="red",line_dash='dashed')
            p.line(x_,-lim99,line_color="red",line_dash='dashed')
        if add_labels:
            labelsX = LabelSet(x='x', y='y', text='ObsID', level='glyph',x_offset=5, y_offset=5, source=source, render_mode='canvas')
            p.add_layout(labelsX)
        p.xaxis.axis_label = 'Observation'
        p.yaxis.axis_label = 't ['+str(dim[0])+']'
        show(p)      
    else: # YES CLASSIDS
        Classes_=np.unique(CLASSID[colorby]).tolist()
        A=len(Classes_)
        colormap =cm.get_cmap("rainbow")
        different_colors=A
        color_mapping=colormap(np.linspace(0,1,different_colors),1,True)
        bokeh_palette=["#%02x%02x%02x" % (r, g, b) for r, g, b in color_mapping[:,0:3]]  
        rnd_num=str(int(np.round(1000*np.random.random_sample())))               
        output_file("Score_Line_"+rnd_num+".html",title='Score Line t['+str(dim[0])+ ']') 

        y_=T_matrix[:,[dim[0]-1]]  
        x_=list(range(1,y_.shape[0]+1))        
        
        TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
        TOOLTIPS = [
                ("Obs#", "@ObsNum"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID"),
                ("Class:","@Class")
                ]        
        classid_=list(CLASSID[colorby])
        legend_it = []
        
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,toolbar_location="above",plot_width=plotwidth,plot_height=plotheight, title='Score Line t['+str(dim[0])+ ']')

        for classid_in_turn in Classes_:
            x_aux=[]
            y_aux=[]
            obsid_aux=[]
            classid_aux=[]
            obsnum_aux=[]
            for i in list(range(len(ObsID_))):
                if classid_[i]==classid_in_turn:
                    x_aux.append(x_[i])
                    y_aux.append(y_[i][0])
                    obsid_aux.append(ObsID_[i])
                    obsnum_aux.append(ObsNum_[i])
                    classid_aux.append(classid_in_turn)
            source = ColumnDataSource(data=dict(x=x_aux, y=y_aux,ObsID=obsid_aux,ObsNum=obsnum_aux,Class=classid_aux))        
            color_=bokeh_palette[Classes_.index(classid_in_turn)]
            c=p.circle('x','y',source=source,color=color_)
            if plotline:
                c1=p.line('x','y',source=source,color=color_)    
             #added to allow numbers in classids   
            aux_=classid_in_turn
            if isinstance(aux_,(float,int)):
                aux_=str(aux_)
             #        
            if add_legend and plotline:    
              #  legend_it.append((classid_in_turn, [c,c1]))
                legend_it.append((aux_, [c,c1]))
            if add_legend and not(plotline):
              #  legend_it.append((classid_in_turn, [c]))
                legend_it.append((aux_, [c]))

            if add_labels:
                labelsX = LabelSet(x='x', y='y', text='ObsID', level='glyph',x_offset=5, y_offset=5, source=source, render_mode='canvas')
                p.add_layout(labelsX)
        if add_ci:
            lim95,lim99=phi.single_score_conf_int(mvmobj['T'][:,[dim[0]-1]])
            p.line(x_, lim95,line_color="gold",line_dash='dashed')
            p.line(x_,-lim95,line_color="gold",line_dash='dashed')
            p.line(x_, lim99,line_color="red",line_dash='dashed')
            p.line(x_,-lim99,line_color="red",line_dash='dashed')   
        p.xaxis.axis_label = 'Observation'
        p.yaxis.axis_label = 't ['+str(dim[0])+']'
        if add_legend:
            legend = Legend(items=legend_it, location='top_right')
            p.add_layout(legend, 'right')
            legend.click_policy="hide"
        show(p)
    return  



def diagnostics(mvmobj,*,Xnew=False,Ynew=False,score_plot_xydim=False,plotwidth=600,ht2_logscale=False,spe_logscale=False):

def diagnostics(mvmobj,*,Xnew=False,Ynew=False,score_plot_xydim=False,plotwidth=600,ht2_logscale=False,spe_logscale=False):
    """
    Plot calculated Hotelling's T2 and SPE
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pca or phi.pls
    
    Xnew/Ynew:     Data used to calculate diagnostics[numpy arrays or pandas dataframes] 
    
    optional:
        
    score_plot_xydim: will add a score scatter plot at the bottom 
                      if sent with a list of [dimx, dimy] where dimx/dimy 
                      are integers and refer to the latent space to plot
                      in the x and y axes of the scatter plot. e.g. [1,2] will
                      add a t1-t2 plot 
    
    """
    
    if isinstance(score_plot_xydim,np.bool):
        add_score_plot = False
    else:
        add_score_plot = True
        
    if isinstance(Xnew,np.bool): #No Xnew was given need to plot all from model
        if 'obsidX' in mvmobj:
            ObsID_=mvmobj['obsidX']
        else:
            ObsID_ = []
            for n in list(np.arange(mvmobj['T'].shape[0])+1):
                ObsID_.append('Obs #'+str(n))  
                              
        Obs_num = np.arange(mvmobj['T'].shape[0])+1
        
        if add_score_plot and not(isinstance(score_plot_xydim,np.bool)):
            t_x  = mvmobj['T'][:,[score_plot_xydim[0]-1]]
            t_y  = mvmobj['T'][:,[score_plot_xydim[1]-1]]
        else:
            add_score_plot = False
        t2_   = mvmobj['T2']
        spex_ = mvmobj['speX']
        
        if ht2_logscale:
            t2_=np.log10(t2_)
        if spe_logscale:
            spex_= np.log10(spex_)
            
            
        if not(add_score_plot):
            if 'Q' in mvmobj:
                spey_=1
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,t2=t2_,spex=spex_,spey=mvmobj['speY']))  
            else:
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,t2=t2_,spex=spex_)) 
        else:
            if 'Q' in mvmobj:
                spey_=1
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,t2=t2_,spex=spex_,spey=mvmobj['speY'],tx=t_x,ty=t_y))  
            else:
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,t2=t2_,spex=spex_,tx=t_x,ty=t_y))
    else: #Xnew was given
        if isinstance(Xnew,np.ndarray):
            ObsID_ = []
            for n in list(np.arange(Xnew.shape[0])+1):
                ObsID_.append('Obs #'+str(n))  
        elif isinstance(Xnew,pd.DataFrame):
            X_=np.array(Xnew.values[:,1:]).astype(float)
            ObsID_ = Xnew.values[:,0].astype(str)
            ObsID_ = ObsID_.tolist()
            
        
        
        if add_score_plot and not(isinstance(score_plot_xydim,np.bool)):
            if 'Q' in mvmobj:  
                xpred=phi.pls_pred(X_,mvmobj)
            else:
                xpred=phi.pca_pred(X_,mvmobj)
            T_matrix=xpred['Tnew']
            t_x  = T_matrix[:,[score_plot_xydim[0]-1]]
            t_y  = T_matrix[:,[score_plot_xydim[1]-1]]
        else:
            add_score_plot = False
        
        t2_ = phi.hott2(mvmobj,Xnew=Xnew)
        
        Obs_num = np.arange(t2_.shape[0])+1
        
        if 'Q' in mvmobj and not(isinstance(Ynew,np.bool)):
            spex_,spey_ = phi.spe(mvmobj,Xnew,Ynew=Ynew)
        else:
            spex_ = phi.spe(mvmobj,Xnew)
            spey_ = False
            
        if ht2_logscale:
            t2_=np.log10(t2_)
        if spe_logscale:
            spex_= np.log10(spex_)
        ObsNum_=[]    
        for n in list(range(1,len(ObsID_)+1)):
            ObsNum_.append('Obs #'+str(n))  
                       
                       
        if not(add_score_plot):
            if 'Q' in mvmobj and not(isinstance(Ynew,np.bool)):
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,ObsNum=ObsNum_,t2=t2_,spex=spex_,spey=spey_))  
            else:
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,ObsNum=ObsNum_,t2=t2_,spex=spex_)) 
        else:
            if 'Q' in mvmobj and not(isinstance(Ynew,np.bool)):
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,ObsNum=ObsNum_,t2=t2_,spex=spex_,spey=spey_,tx=t_x,ty=t_y))  
            else:
                source = ColumnDataSource(data=dict(x=Obs_num, ObsID=ObsID_,ObsNum=ObsNum_,t2=t2_,spex=spex_,tx=t_x,ty=t_y))
    TOOLS = "save,wheel_zoom,box_zoom,reset,lasso_select"
    TOOLTIPS = [
            ("Obs #", "@x"),
            ("(x,y)", "($x, $y)"),
            ("Obs: ","@ObsID")
            ]
    
    rnd_num=str(int(np.round(1000*np.random.random_sample())))               
    output_file("Diagnostics"+rnd_num+".html",title='Diagnostics') 
    p = figure(tools=TOOLS, tooltips=TOOLTIPS, plot_width=plotwidth, title="Hotelling's T2")
    p.circle('x','t2',source=source)
    if ht2_logscale:
        p.line([0,Obs_num[-1]],[np.log10(mvmobj['T2_lim95']),np.log10(mvmobj['T2_lim95'])],line_color='gold')
        p.line([0,Obs_num[-1]],[np.log10(mvmobj['T2_lim99']),np.log10(mvmobj['T2_lim99'])],line_color='red')
    else:        
        p.line([0,Obs_num[-1]],[mvmobj['T2_lim95'],mvmobj['T2_lim95']],line_color='gold')
        p.line([0,Obs_num[-1]],[mvmobj['T2_lim99'],mvmobj['T2_lim99']],line_color='red')
    
    
    p.xaxis.axis_label = 'Observation sequence'
    p.yaxis.axis_label = "HT2"
    p_list=[p]
    
    p = figure(tools=TOOLS, tooltips=TOOLTIPS, plot_width=plotwidth, title='SPE X')
    p.circle('x','spex',source=source)
    
    if spe_logscale:
        p.line([0,Obs_num[-1]],[np.log10(mvmobj['speX_lim95']),np.log10(mvmobj['speX_lim95'])],line_color='gold')
        p.line([0,Obs_num[-1]],[np.log10(mvmobj['speX_lim99']),np.log10(mvmobj['speX_lim99'])],line_color='red')
    else:  
        p.line([0,Obs_num[-1]],[mvmobj['speX_lim95'],mvmobj['speX_lim95']],line_color='gold')
        p.line([0,Obs_num[-1]],[mvmobj['speX_lim99'],mvmobj['speX_lim99']],line_color='red')
    p.xaxis.axis_label = 'Observation sequence'
    p.yaxis.axis_label = 'SPE X-Space'
    p_list.append(p)
    
    p = figure(tools=TOOLS, tooltips=TOOLTIPS, plot_width=plotwidth, title='Outlier Map')
    p.circle('t2','spex',source=source)
    if ht2_logscale:
        vline = Span(location=np.log10(mvmobj['T2_lim99']), dimension='height', line_color='red', line_width=1)
    else:
        vline = Span(location=mvmobj['T2_lim99'], dimension='height', line_color='red', line_width=1)
    if spe_logscale:    
        hline = Span(location=np.log10(mvmobj['speX_lim99']), dimension='width', line_color='red', line_width=1)
    else:
        hline = Span(location=mvmobj['speX_lim99'], dimension='width', line_color='red', line_width=1)
    p.renderers.extend([vline, hline])
    
    p.xaxis.axis_label = "Hotelling's T2"
    p.yaxis.axis_label = 'SPE X-Space'
    p_list.append(p)
    
    
    if 'Q' in mvmobj and not(isinstance(spey_,np.bool)):
        p = figure(tools=TOOLS, tooltips=TOOLTIPS, plot_height=400, title='SPE Y')
        p.circle('x','spey',source=source)
        p.line([0,Obs_num[-1]],[mvmobj['speY_lim95'],mvmobj['speY_lim95']],line_color='gold')
        p.line([0,Obs_num[-1]],[mvmobj['speY_lim99'],mvmobj['speY_lim99']],line_color='red')
        p.xaxis.axis_label = 'Observation sequence'
        p.yaxis.axis_label = 'SPE Y-Space'
        p_list.append(p)
    if add_score_plot:
        p = figure(tools=TOOLS, tooltips=TOOLTIPS, plot_width=plotwidth, title='Score Scatter')
        p.circle('tx', 'ty', source=source,size=7)
        
        T_aux1=mvmobj['T'][:,[score_plot_xydim[0]-1]]
        T_aux2=mvmobj['T'][:,[score_plot_xydim[1]-1]]
        T_aux = np.hstack((T_aux1,T_aux2))
        st=(T_aux.T @ T_aux)/T_aux.shape[0]
        [xd95,xd99,yd95p,yd95n,yd99p,yd99n]=phi.scores_conf_int_calc(st,mvmobj['T'].shape[0])
        p.line(xd95,yd95p,line_color="gold",line_dash='dashed')
        p.line(xd95,yd95n,line_color="gold",line_dash='dashed')
        p.line(xd99,yd99p,line_color="red",line_dash='dashed')
        p.line(xd99,yd99n,line_color="red",line_dash='dashed') 
        p.xaxis.axis_label = 't ['+str(score_plot_xydim[0])+']'
        p.yaxis.axis_label = 't ['+str(score_plot_xydim[1])+']'
        # Vertical line
        vline = Span(location=0, dimension='height', line_color='black', line_width=2)
        # Horizontal line
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        p.renderers.extend([vline, hline])
        #Do another p.figure
        p_list.append(p)
    
    show(column(p_list)) 
    return

def predvsobs(mvmobj,X,Y,*,CLASSID=False,colorby=False,x_space=False):

def predvsobs(mvmobj,X,Y,*,CLASSID=False,colorby=False,x_space=False):
    """
    Plot observed vs predicted values
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A model created with phi.pca or phi.pls
    
    X/Y:     Data [numpy arrays or pandas dataframes] 
    
    optional:
    CLASSID: Pandas Data Frame with classifiers per observation, each column is a class
    
    colorby: one of the classes in CLASSID to color by
    
    x_space: = 'False' will skip plotting the obs. vs pred for X *default*
               'True' will also plot obs vs pred for X

    """
    num_varX=mvmobj['P'].shape[0]
    
    if isinstance(X,np.ndarray):
        X_=X.copy()
        ObsID_ = []
        for n in list(np.arange(X.shape[0])+1):
            ObsID_.append('Obs #'+str(n))  
        XVarID_ = []
        for n in list(np.arange(X.shape[1])+1):
            XVarID_.append('Var #'+str(n))  
    elif isinstance(X,pd.DataFrame):
        X_=np.array(X.values[:,1:]).astype(float)
        ObsID_ = X.values[:,0].astype(str)
        ObsID_ = ObsID_.tolist()

    if 'varidX' in mvmobj:
        XVar=mvmobj['varidX']
    else:
        XVar = []
        for n in list(np.arange(num_varX)+1):
            XVar.append('XVar #'+str(n))    
                        
    if 'Q' in mvmobj:
        num_varY=mvmobj['Q'].shape[0]
        if 'varidY' in mvmobj:
            YVar=mvmobj['varidY']
        else:
            YVar = []
            for n in list(np.arange(num_varY)+1):
                YVar.append('YVar #'+str(n))
                        

    if isinstance(Y,np.ndarray):
        Y_=Y.copy()
    elif isinstance(Y,pd.DataFrame):
        Y_=np.array(Y.values[:,1:]).astype(float)

            
    if 'Q' in mvmobj:  
        pred=phi.pls_pred(X_,mvmobj)
        yhat=pred['Yhat']
        if x_space:
            xhat=pred['Xhat']
        else:
            xhat=False
    else:
        x_space=True
        pred=phi.pca_pred(X_,mvmobj)
        xhat=pred['Xhat']
        yhat=False
        
    TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
    TOOLTIPS = [
                ("index", "$index"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID")
                ]
    
    if isinstance(CLASSID,np.bool): # No CLASSIDS
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("ObsvsPred_"+rnd_num+".html",title='ObsvsPred')
        plot_counter=0
        
        if not(isinstance(yhat,np.bool)): #skip if PCA model sent
            for i in list(range(Y_.shape[1])):
                x_ = Y_[:,i]
                y_ = yhat[:,i]          
                min_value = np.nanmin([np.nanmin(x_),np.nanmin(y_)])
                max_value = np.nanmax([np.nanmax(x_),np.nanmax(y_)])
                
                source = ColumnDataSource(data=dict(x=x_, y=y_,ObsID=ObsID_))
                #p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=600, plot_height=600, title=YVar[i])
                p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=600, plot_height=600, title=YVar[i],x_range=(min_value, max_value),y_range=(min_value, max_value))
                p.circle('x', 'y', source=source,size=7,color='darkblue')
                p.line([min_value,max_value],[min_value,max_value],line_color='cyan',line_dash='dashed')
                p.xaxis.axis_label ='Observed'
                p.yaxis.axis_label ='Predicted'
                if plot_counter==0:
                    p_list=[p]
                else:
                    p_list.append(p)
                plot_counter = plot_counter+1
                
        if x_space: #
            for i in list(range(X_.shape[1])):
                x_ = X_[:,i]
                y_ = xhat[:,i]    
                min_value = np.nanmin([np.nanmin(x_),np.nanmin(y_)])
                max_value = np.nanmax([np.nanmax(x_),np.nanmax(y_)])
 
                source = ColumnDataSource(data=dict(x=x_, y=y_,ObsID=ObsID_))
                p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=600, plot_height=600, title=XVar[i],x_range=(min_value, max_value),y_range=(min_value, max_value))
                p.circle('x', 'y', source=source,size=7,color='darkblue')
                p.line([min_value,max_value],[min_value,max_value],line_color='cyan',line_dash='dashed')
                p.xaxis.axis_label ='Observed'
                p.yaxis.axis_label ='Predicted'
                if plot_counter==0:
                    p_list=[p]
                else:
                    p_list.append(p)
                plot_counter = plot_counter+1
        show(column(p_list))
        
    else: # YES CLASSIDS
        Classes_=np.unique(CLASSID[colorby]).tolist()
        different_colors=len(Classes_)
        colormap =cm.get_cmap("rainbow")
        color_mapping=colormap(np.linspace(0,1,different_colors),1,True)
        bokeh_palette=["#%02x%02x%02x" % (r, g, b) for r, g, b in color_mapping[:,0:3]]  
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("ObsvsPred_"+rnd_num+".html",title='ObsvsPred')      
        classid_=list(CLASSID[colorby])
        
        plot_counter=0
        
        if not(isinstance(yhat,np.bool)): #skip if PCA model sent
            for i in list(range(Y_.shape[1])):
                x_ = Y_[:,i]
                y_ = yhat[:,i]
                min_value = np.nanmin([np.nanmin(x_),np.nanmin(y_)])
                max_value = np.nanmax([np.nanmax(x_),np.nanmax(y_)])
                p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=600, plot_height=600, title=YVar[i],x_range=(min_value, max_value),y_range=(min_value, max_value))
                for classid_in_turn in Classes_:
                    x_aux=[]
                    y_aux=[]
                    obsid_aux=[]
                    classid_aux=[]
                    for i in list(range(len(ObsID_))):
                        if classid_[i]==classid_in_turn and not(np.isnan(x_[i])):
                            x_aux.append(x_[i])
                            y_aux.append(y_[i])
                            obsid_aux.append(ObsID_[i])
                            classid_aux.append(classid_in_turn)
                    source = ColumnDataSource(data=dict(x=x_aux, y=y_aux,ObsID=obsid_aux,Class=classid_aux))        
                    color_=bokeh_palette[Classes_.index(classid_in_turn)]
                    p.circle('x','y',source=source,color=color_,legend_label=classid_in_turn)
                    p.line([min_value,max_value],[min_value,max_value],line_color='cyan',line_dash='dashed')
                p.xaxis.axis_label ='Observed'
                p.yaxis.axis_label ='Predicted'
                p.legend.click_policy="hide"
                p.legend.location = "top_left"
                if plot_counter==0:
                    p_list=[p]
                    plot_counter = plot_counter+1
                else:
                    p_list.append(p)
              
        if x_space: #
            for i in list(range(X_.shape[1])):
                x_ = X_[:,i]
                y_ = xhat[:,i]
                min_value = np.nanmin([np.nanmin(x_),np.nanmin(y_)])
                max_value = np.nanmax([np.nanmax(x_),np.nanmax(y_)])
                p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=600, plot_height=600, title=XVar[i],x_range=(min_value, max_value),y_range=(min_value, max_value))
                for classid_in_turn in Classes_:
                    x_aux=[]
                    y_aux=[]
                    obsid_aux=[]
                    classid_aux=[]
                    for i in list(range(len(ObsID_))):
                        if classid_[i]==classid_in_turn and not(np.isnan(x_[i])):
                            x_aux.append(x_[i])
                            y_aux.append(y_[i])
                            obsid_aux.append(ObsID_[i])
                            classid_aux.append(classid_in_turn)
                    source = ColumnDataSource(data=dict(x=x_aux, y=y_aux,ObsID=obsid_aux,Class=classid_aux))        
                    color_=bokeh_palette[Classes_.index(classid_in_turn)]
                    p.circle('x','y',source=source,color=color_,legend_label=classid_in_turn)
                    p.line([min_value,max_value],[min_value,max_value],line_color='cyan',line_dash='dashed')
                p.xaxis.axis_label ='Observed'
                p.yaxis.axis_label ='Predicted'
                p.legend.click_policy="hide"
                p.legend.location = "top_left"
                if plot_counter==0:
                    p_list=[p]
                    plot_counter = plot_counter+1
                else:
                    p_list.append(p)
        show(column(p_list))
    return    

def contributions_plot(mvmobj,X,cont_type,*,Y=False,from_obs=False,to_obs=False,lv_space=False,plotwidth=800,plotheight=600,xgrid=False):

def contributions_plot(mvmobj,X,cont_type,*,Y=False,from_obs=False,to_obs=False,lv_space=False,plotwidth=800,plotheight=600,xgrid=False):
    """
    Calculate contributions to diagnostics
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj : A dictionary created by phi.pls or phi.pca
    
    X/Y:     Data [numpy arrays or pandas dataframes] - Y space is optional
    
    cont_type: 'ht2'
               'spe'
               'scores'
               
    from_obs: Scalar or list of scalars with observation(s) number(s) | first element is #0
              - OR -
              Strings or list of strings with observation(s) name(s) [if X/Y are pandas data frames]
              Used to off set calculations for scores or ht2
              "False' will calculate with respect to origin *default if not sent*
              
    to_obs: Scalar or list of scalars with observation(s) number(s)| first element is #0
              - OR -
            Strings or list of strings with observation(s) name(s) [if X/Y are pandas data frames]
            To calculate contributions for
            
            *Note: from_obs is ignored when cont_type='spe'*
            
    lv_space: Latent spaces over which to do the calculations [applicable to 'ht2' and 'scores']
    """
    good_to_go=True
    if isinstance(X,pd.DataFrame):
        ObsID=X.values[:,0].tolist()
        if isinstance(to_obs,str):
            to_obs_=ObsID.index(to_obs)
        elif isinstance(to_obs,int):
            to_obs_=to_obs
        elif isinstance(to_obs,list):
            if isinstance(to_obs[0],str):
                to_obs_=[]
                for o in to_obs:
                    to_obs_.append(ObsID.index(o))
            elif isinstance(to_obs[0],int):
                to_obs_=to_obs.copy()
        elif isinstance(to_obs,np.bool):
            good_to_go=False
        if not(isinstance(from_obs,np.bool)):
            if isinstance(from_obs,str):
                from_obs_=ObsID.index(from_obs)
            elif isinstance(from_obs,int):
                from_obs_=from_obs
            elif isinstance(from_obs,list):
                if isinstance(from_obs[0],str):
                    from_obs_=[]
                    for o in from_obs:
                        from_obs_.append(ObsID.index(o))
                elif isinstance(from_obs[0],int):
                    from_obs_=from_obs.copy()
        else:
            from_obs_=False
    else:
        if isinstance(to_obs,int) or isinstance(to_obs,list):
            to_obs_=to_obs.copy()
        else:
            good_to_go=False    
    if cont_type=='scores' and not(isinstance(Y,np.bool)):
        Y=False
        
    if isinstance(Y,np.bool) and good_to_go:
        Xconts=phi.contributions(mvmobj,X,cont_type,Y=False,from_obs=from_obs_,to_obs=to_obs_,lv_space=lv_space)
        Yconts=False
    elif not(isinstance(Y,np.bool)) and good_to_go and ('Q' in mvmobj) and cont_type=='spe':
        Xconts,Yconts=phi.contributions(mvmobj,X,cont_type,Y=Y,from_obs=from_obs_,to_obs=to_obs_,lv_space=lv_space)
    
    if 'varidX' in mvmobj:
        XVar=mvmobj['varidX']
    else:
        XVar = []
        for n in list(np.arange(mvmobj['P'].shape[0])+1):
            XVar.append('XVar #'+str(n))               

    rnd_num=str(int(np.round(1000*np.random.random_sample())))
    output_file("Contributions"+rnd_num+".html",title='Contributions')
    if isinstance(from_obs,list):
        from_txt=", ".join(map(str, from_obs))
        from_txt=" from obs: "+from_txt
    elif isinstance(from_obs,int):
        from_txt=" from obs: "+str(from_obs)
    elif isinstance(from_obs,str):    
        from_txt=" from obs: " + from_obs
    else:
        from_txt=""
    if isinstance(to_obs,list):
        to_txt=", ".join(map(str, to_obs))
        to_txt=", to obs: "+to_txt
    elif isinstance(to_obs,str):    
        to_txt=", to obs: " + to_obs
    elif isinstance(to_obs,int):
        to_txt =", to obs: "+ str(to_obs)
    else:
        to_txt=""
    
    p = figure(x_range=XVar, plot_height=plotheight,plot_width=plotwidth, title="Contributions Plot"+from_txt+to_txt,
                    tools="save,box_zoom,pan,reset")
    p.vbar(x=XVar, top=Xconts[0].tolist(), width=0.5)
    p.ygrid.grid_line_color = None    
    if xgrid:
        p.xgrid.grid_line_color = 'lightgray'
    else:
        p.xgrid.grid_line_color = None   
    p.yaxis.axis_label = 'Contributions to '+cont_type
    hline = Span(location=0, dimension='width', line_color='black', line_width=2)
    p.renderers.extend([hline])
    p.xaxis.major_label_orientation = 45
    p_list=[p]
    
    if not(isinstance(Yconts,np.bool)):
        if 'varidY' in mvmobj:
            YVar=mvmobj['varidY']
        else:
            YVar = []
            for n in list(np.arange(mvmobj['Q'].shape[0])+1):
                YVar.append('YVar #'+str(n))               
        
        p = figure(x_range=YVar, plot_height=plotheight,plot_width=plotwidth, title="Contributions Plot",
                    tools="save,box_zoom,pan,reset")
        p.vbar(x=YVar, top=Yconts[0].tolist(), width=0.5)
        p.ygrid.grid_line_color = None    
        if xgrid:
            p.xgrid.grid_line_color = 'lightgray'
        else:
            p.xgrid.grid_line_color = None   
        p.yaxis.axis_label = 'Contributions to '+cont_type
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        p.renderers.extend([hline])
        p.xaxis.major_label_orientation = 45
        p_list.append(p)
        
    show(column(p_list))  
    return

def plot_spectra(X,*,xaxis=False,plot_title='Main Title',tab_title='Tab Title',xaxis_label='X- axis',yaxis_label='Y- axis'): 

def plot_spectra(X,*,xaxis=False,plot_title='Main Title',tab_title='Tab Title',xaxis_label='X- axis',yaxis_label='Y- axis'): 
    """
    Simple way to plot Spectra with Bokeh. 
    Programmed by Salvador Garcia-Munoz
    ([email protected] ,[email protected])
    
    X:      A numpy array or a pandas object with Spectra to be plotted
    xaxis:  wavenumbers or wavelengths to index the x axis of the plot 
            * ignored if X is a pandas dataframe *
            
    optional: 
    plot_title
    tab_title
    xaxis_label
    yaxis_label
    
    """
    
    if isinstance(X,pd.DataFrame):
        x=X.columns[1:].tolist()
        x=np.array(x)
        x=np.reshape(x,(1,-1))
        y=X.values[:,1:].astype(float)
    elif isinstance(X,np.ndarray):
        y=X.copy()
        if isinstance(xaxis,np.ndarray):
            x=xaxis
            x=np.reshape(x,(1,-1))
        elif isinstance(xaxis,list):
            x=np.array(xaxis)
            x=np.reshape(x,(1,-1))
        elif isinstance(xaxis,np.bool):
            x=np.array(list(range(X.shape[1])))
            x=np.reshape(x,(1,-1))
    rnd_num=str(int(np.round(1000*np.random.random_sample())))                
    output_file("Spectra"+rnd_num+".html",title=tab_title)

    p = figure(title=plot_title)
    p.xaxis.axis_label = xaxis_label
    p.yaxis.axis_label = yaxis_label
    p.multi_line(x.tolist()*y.shape[0],y.tolist())
    show(p)
    return

def plot_line_pd(X,col_name,*,plot_title='Main Title',tab_title='Tab Title',xaxis_label='X- axis',plotheight=400,plotwidth=600): 

def plot_line_pd(X,col_name,*,plot_title='Main Title',tab_title='Tab Title',xaxis_label='X- axis',plotheight=400,plotwidth=600): 
    """
    Simple way to plot a column of a Pandas DataFrame with Bokeh. 
    Programmed by Salvador Garcia-Munoz
    ([email protected] ,[email protected])
    
    X:      A a pandas object with Data to be plotted
    col_name: The name of the column to plot

            
    optional: 
    plot_title
    tab_title
    xaxis_label
    yaxis_label
    plotheight
    plotwidth
    
    """
    
    if isinstance(col_name,str):
        col_name=[col_name]
    first_plot=True
    
    TOOLS = "save,wheel_zoom,box_zoom,pan,reset,box_select,lasso_select"
    TOOLTIPS = [
                ("Obs #", "@ObsNum"),
                ("(x,y)", "($x, $y)"),
                ("Obs: ","@ObsID")
                ] 
    rnd_num=str(int(np.round(1000*np.random.random_sample())))          
    output_file("LinePlot"+rnd_num+".html",title=tab_title)
    
    for this_col_name in col_name:
        ObsID_=X.values[:,0]
        ObsID_=ObsID_.tolist()
        aux=X.loc[:,this_col_name]
        y_=aux.values  
        x_=list(range(1,len(ObsID_)+1))
        ObsNum_=[]    
        for n in list(range(1,len(ObsID_)+1)):
            ObsNum_.append('Obs #'+str(n))     
        if not(first_plot):                   
            plot_title=''    
        p = figure(tools=TOOLS, tooltips=TOOLTIPS,plot_width=plotwidth,plot_height=plotheight,title=plot_title)        
        source = ColumnDataSource(data=dict(x=x_, y=y_,ObsID=ObsID_,ObsNum=ObsNum_))
        p.xaxis.axis_label = xaxis_label
        p.yaxis.axis_label = this_col_name
        p.line('x', 'y', source=source)
        p.circle('x', 'y', source=source)
        if first_plot:
            p_list=[p]
            first_plot=False
        else:
            p_list.append(p)
    show(column(p_list))  
         
    return

def mb_weights(mvmobj,*,plotwidth=600,plotheight=400):

def mb_weights(mvmobj,*,plotwidth=600,plotheight=400):
    """
    Super weights for Multi-block models
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A multi-block PLS model created with phi.mbpls
    """
    A= mvmobj['T'].shape[1]
    lv_prefix='LV #'        
    lv_labels = []   
    for a in list(np.arange(A)+1):
        lv_labels.append(lv_prefix+str(a))    
    XVar=mvmobj['Xblocknames']        
    for i in list(np.arange(A)):
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("blockweights_"+rnd_num+".html",title="Block Weights")         
        px = figure(x_range=XVar, title="Block weights for MBPLS"+lv_labels[i],
             tools="save,box_zoom,hover,reset", tooltips=[("Var:","@x_")],plot_width=plotwidth,plot_height=plotheight)   
        source1 = ColumnDataSource(data=dict(x_=XVar, y_=mvmobj['Wt'][:,i].tolist(),names=XVar)) 
        px.vbar(x='x_', top='y_', source=source1,width=0.5)
        px.y_range.range_padding = 0.1
        px.ygrid.grid_line_color = None
        px.axis.minor_tick_line_color = None
        px.outline_line_color = None
        px.yaxis.axis_label = 'Wt'+str(i+1)+']'
        px.xaxis.major_label_orientation = 45  
        hline = Span(location=0, dimension='width', line_color='black', line_width=2)
        px.renderers.extend([hline])
        if i==0:
            p_list=[px]
        else:
            p_list.append(px)
    show(column(p_list))  

    return

        

def mb_r2pb(mvmobj,*,plotwidth=600,plotheight=400):

def mb_r2pb(mvmobj,*,plotwidth=600,plotheight=400):
    """
    Super weights for Multi-block models
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A multi-block PLS model created with phi.mbpls
    """
    A= mvmobj['T'].shape[1]
    lv_prefix='LV #'        
    lv_labels = []   
    for a in list(np.arange(A)+1):
        lv_labels.append(lv_prefix+str(a))    
    r2pbX_dict = {'XVar': mvmobj['Xblocknames']}
    XVar=mvmobj['Xblocknames']        
    for i in list(np.arange(A)):
        r2pbX_dict.update({lv_labels[i] : mvmobj['r2pbX'][:,i].tolist()})
        rnd_num=str(int(np.round(1000*np.random.random_sample())))
        output_file("r2perblock"+rnd_num+".html",title="R2 per Block") 
        colormap =cm.get_cmap("rainbow")
        different_colors=A
        color_mapping=colormap(np.linspace(0,1,different_colors),1,True)
        bokeh_palette=["#%02x%02x%02x" % (r, g, b) for r, g, b in color_mapping[:,0:3]]                 
        px = figure(x_range=XVar, title="r2 per Block for MBPLS",
             tools="save,box_zoom,hover,reset", tooltips="$name @XVar: @$name",plot_width=plotwidth,plot_height=plotheight)        
        px.vbar_stack(lv_labels, x='XVar', width=0.9,color=bokeh_palette,source=r2pbX_dict)
        px.y_range.range_padding = 0.1
        px.ygrid.grid_line_color = None
        px.axis.minor_tick_line_color = None
        px.outline_line_color = None
        px.yaxis.axis_label = 'R2 per Block per LV'
        px.xaxis.major_label_orientation = 45      
    show(px)
    return


def mb_vip(mvmobj,*,plotwidth=600,plotheight=400):

def mb_vip(mvmobj,*,plotwidth=600,plotheight=400):
    """
    Super weights for Multi-block models
    by Salvador Garcia-Munoz 
    ([email protected] ,[email protected])
    
    mvmobj: A multi-block PLS model created with phi.mbpls
    """
    A= mvmobj['T'].shape[1]
   
    XVar=mvmobj['Xblocknames']        
    Wt=mvmobj['Wt']
    r2y=mvmobj['r2y']
    vip=np.zeros((Wt.shape[0],1))
    if A>1:
        for a in list(range(A)):
            vip=vip+Wt[:,[a]]*r2y[a]
    else:
        vip=Wt[:,[0]]*r2y
        
    vip=np.reshape(vip,-1)
    index=np.argsort(vip)
    index=index[::-1]
    XVar_=[XVar[i] for i in index]
    XVar = XVar_
    vip=vip[index]
    rnd_num=str(int(np.round(1000*np.random.random_sample())))
    output_file("blockvip"+rnd_num+".html",title="Block VIP") 
    source1 = ColumnDataSource(data=dict(x_=XVar, y_=vip.tolist(),names=XVar))         
    px = figure(x_range=XVar, title="Block VIP for MBPLS",
         tools="save,box_zoom,hover,reset",tooltips=[("Block:","@x_")],plot_width=plotwidth,plot_height=plotheight)   
    
    px.vbar(x='x_', top='y_', source=source1,width=0.5)
    px.y_range.range_padding = 0.1
    px.ygrid.grid_line_color = None
    px.axis.minor_tick_line_color = None
    px.outline_line_color = None
    px.yaxis.axis_label = 'Block VIP'
    px.xaxis.major_label_orientation = 45  
    hline = Span(location=0, dimension='width', line_color='black', line_width=2)
    px.renderers.extend([hline])
    show(px)  
    return

def advanced_visualize(
    file_dir: str,  # Dataset directory.
    target: str,
    load_multi_ex: callable,
    CON_config: dict,
    NN_config: dict,
    show_plot: bool=False
) -> None:
    """
    # TODO: write the doc.
    Predict and Visualize the result.
    """
    print("[IPR]Visualize model result using bokeh...")
    print(f"Building up from container with data at {file_dir}...")

    container = MultivariateContainer(
        file_dir,
        target,
        load_multi_ex,
        CON_config)

    print("Building empty model placeholder...")
    model = MultivariateLSTM(
        container,
        NN_config,
        create_empty=True)

    # The folder where model is stored.
    load_target = input("Model directory >>> ")
    load_target = f"./saved_models/{load_target}/"
    print(f"Loading model from directory: {load_target}...")
    model.load_model(folder_dir=load_target)

    # timeline = pd.DatetimeIndex(container.dataset.index)
    # Time line for x-axis

    # true_y = np.diff(model.container.get_true_y())
    output_file(f"{load_target}visualized.html")
    print(f"Saving plotting html file to {load_target}visualized.html...")

    # ======== Differenced Value ========
    print("Building up forecasting sequence for differenced value...")
    train_yhat = model.predict(model.container.train_X)
    train_yhat = np.squeeze(train_yhat)

    test_yhat = model.predict(model.container.test_X)
    test_yhat = np.squeeze(test_yhat)

    timeline = pd.DatetimeIndex(model.container.dataset.index)
    true_y = np.diff(model.container.get_true_y())

    # ================ Training Set ================
    plot_diff_train = bokeh.plotting.figure(
        title="Differencing value, training set",
        x_axis_label="Date",
        y_axis_label="Differenced Value",
        x_axis_type="datetime",
        plot_width=1400,
        plot_height=400
    )

    plot_diff_train.line(
        timeline[1: len(train_yhat)+1],
        np.squeeze(model.container.train_y),
        color="navy",
        alpha=0.4,
        legend="Training Set Actual Values"
    )

    plot_diff_train.line(
        timeline[1: len(train_yhat)+1],
        train_yhat,
        color="red",
        alpha=0.7,
        legend="Training Set Predicted Values"
    )

    # ================ Testing Set ================
    plot_diff_test = bokeh.plotting.figure(
        title="Differencing value, testing set",
        x_axis_label="Date",
        y_axis_label="Differenced Value",
        x_axis_type="datetime",
        plot_width=1400,
        plot_height=400
    )

    plot_diff_test.line(
        timeline[-len(test_yhat):],
        np.squeeze(model.container.test_y),
        color="navy",
        alpha=0.4,
        legend="Test Set Actual Values"
    )

    plot_diff_test.line(
        timeline[-len(test_yhat):],
        test_yhat,
        color="red",
        alpha=0.7,
        legend="Test Set Predicted Values"
    )

    tab_diff = bokeh.models.widgets.Panel(
        child=column(children=[plot_diff_test, plot_diff_train]),
        title="Differenced values"
    )

    # ======== Raw scale values ========

    train_yhat = model.predict(model.container.train_X)
    train_yhat = model.container.invert_difference(
        train_yhat,
        range(len(train_yhat)), fillnone=True
    )
    train_yhat = np.squeeze(train_yhat).astype(np.float32)

    test_yhat = model.predict(model.container.test_X)
    test_yhat = model.container.invert_difference(
        test_yhat,
        range(
            model.container.num_obs - len(test_yhat),  # Last n observation.
            model.container.num_obs),
        fillnone=True
    )
    test_yhat = np.squeeze(test_yhat).astype(np.float32)

    timeline = pd.DatetimeIndex(model.container.dataset.index)

    raw_plot = figure(
        title="Aggregate Graph: all",
        x_axis_label="Date",
        y_axis_label="Actual Value",
        x_axis_type="datetime",
        plot_width=1400,
        plot_height=400
    )

    raw_plot.line(
        timeline,
        model.container.get_true_y(),
        color="navy",
        alpha=0.3,
        legend="Actual values"
    )

    raw_plot.line(
        timeline,
        train_yhat,
        color="red",
        alpha=0.7,
        legend="Training set predictions"
    )

    raw_plot.line(
        timeline,
        test_yhat,
        color="green",
        alpha=0.7,
        legend="Testing set predictions"
    )

    tab_raw = bokeh.models.widgets.Panel(
        child=raw_plot,
        title="Raw scale"
    )

    # ================ Training information ================
    hist = pd.read_csv(f"{load_target}hist.csv")
    loss = np.squeeze(hist["loss"])
    val_loss = np.squeeze(hist["val_loss"])

    num_epochs = len(loss)

    log_loss = np.log(loss)
    log_val_loss = np.log(val_loss)

    # ======== RAW ========

    loss_plot = figure(
        title="Loss Record(Raw)",
        x_axis_label="Epoch",
        y_axis_label="Loss",
        plot_width=1400,
        plot_height=400
    )

    loss_plot.line(
        range(num_epochs),
        loss,
        color="red",
        alpha=0.7,
        legend="loss on training set"
    )

    loss_plot.line(
        range(num_epochs),
        val_loss,
        color="green",
        alpha=0.7,
        legend="loss on validation set"
    )

    # ======== LOG ========
    log_loss_plot = figure(
        title="Loss Record(Log)",
        x_axis_label="Epoch",
        y_axis_label="Log Loss",
        plot_width=1400,
        plot_height=400
    )

    log_loss_plot.line(
        range(num_epochs),
        log_loss,
        color="red",
        alpha=0.7,
        legend="log loss on training set"
    )

    log_loss_plot.line(
        range(num_epochs),
        log_val_loss,
        color="green",
        alpha=0.7,
        legend="log loss on validation set"
    )

    tab_hist = bokeh.models.widgets.Panel(
        child=column(children=[loss_plot, log_loss_plot]),
        title="Training History"
    )

    # ================ Finalizing ================
    tabs = bokeh.models.widgets.Tabs(tabs=[tab_raw, tab_diff, tab_hist])

    if show_plot:
        bokeh.io.show(tabs)
    else:
        bokeh.io.save(tabs)